Development of Land-River Two-Way Coupling in the Energy Exascale Earth System Model

,


Introduction
Floodplain inundation is a critical process controlling water and biogeochemical cycles at the land-river interface (Scott et al., 2019;Talbot et al., 2018;Tockner and Stanford, 2002).
Specifically, periodic flooding replenishes the soil on the floodplain, provides nutrients for vegetation, and creates habitats for animals.In the coastal areas, flooding in tidal rivers can add saline water to the floodplain and regional groundwater due to the seawater and freshwater interactions (Yabusaki et al., 2020).Two-way coupling of land and river components in Earth System Models (ESMs) is necessary to accurately simulate the impacts of floodplain inundation on water and biogeochemical cycles.Simple macroscale inundation schemes have been coupled to global river routing models to simulate floodplain inundation dynamics (de Paiva et al., 2013;Decharme et al., 2008;Getirana et al., 2012;Luo et al., 2017;Yamazaki et al., 2011), and including such schemes has improved skill in predicting streamflow (Decharme et al., 2012;Yamazaki et al., 2011).Macroscale inundation schemes model floodplain as an impervious reservoir that stores the inundated water when water in the river exceeds the channel capacity.The inundated water stored in the floodplain is released back to the river when the volume of water in the river is below channel capacity.However, the land, river, and ocean components of most ESMs are oneway coupled such that the runoff generated by the land model is transported to oceans through river network (Golaz et al., 2019;Jolley and Wheater, 1997) and the floodplain inundation does not influence the land surface processes.Recently, two-way land-river coupling at global (Decharme et al., 2012) and regional scales (Dadson et al., 2010;Getirana et al., 2021) has shown significant impacts of floodplain inundation on land processes.Decharme et al. (2012) coupled a floodplain inundation scheme (Decharme et al., 2008) with a Land Surface Model (LSM), and their simulations at global scale showed that two-way coupling increased evapotranspiration (ET) and improved the streamflow simulation.Decharme et al. (2019) further reported higher ET with a drop of global inundation extents and an increase of soil moisture due to two-way coupling.Dadson et al. (2010) evaluated the effects of two-way coupling on ET at the Niger inland delta.Although infiltration on the floodplain is not implemented in their study, ET was found to increase significantly due to the open water ET from the inundation water, suggesting the importance of including two-way coupling for understanding land-atmosphere coupling.Simulation by Getirana et al. (2021) that included infiltration of floodplain water showed a significant increase of soil moisture and ET, a decrease in surface temperature, and improved simulated streamflow, which are consistent with the global study of Decharme et al. (2019).Over the tropical regions, including the water exchange on floodplain through infiltration was found to be critical to improve the water cycle and land surface flux simulation (Miguez-Macho and Fan, 2012;Schrapffer et al., 2020).Additionally, Chaney et al. (2020) showed an increase in the spatial heterogeneity of ET and soil moisture in a two-way coupled simulation at high spatial resolution (~1km).
The simulation of floodplain inundation in the coarse-scale (~100km) river component of current generation ESMs requires parameterization of fine-scale topography.To simulate the dynamics of the inundated area, a relationship between flood water volume and inundated area is required.In current large scale land models (Luo et al., 2017;Yamazaki et al., 2011), this relationship is usually represented by a simple 2-D elevation profile (Figure 1a) based on subkilometer topography data (Figure 1b) by assuming that inundation always occurs from lower to higher elevation.Although some studies have demonstrated satisfactory performance of the subgrid topography schemes in capturing basin-averaged inundation dynamics (Decharme et al., 2012;Getirana et al., 2012;Mao et al., 2019;Yamazaki et al., 2011), the simulated spatial variation of inundated area at global scale remains highly uncertain (Mao et al., 2019).This is because by ignoring the hydrologic connectivity, the use of elevation profile in coarse resolution ESMs with structured mesh can result in unrealistic inundated areas that comprise of spatially disconnected flooded regions within a grid cell (see an example in Figure 1c).Another scheme to predict inundation area is the Height Above Nearest Drainage (HAND) methodology (Afshari et al., 2018;Maidment, 2017;Zheng et al., 2018).The HAND data, which is typically derived from high-resolution (e.g., 10m) DEM, aims to capture the role of river network (Liu et al., 2018).
While the HAND approaches accurately capture the hydraulic connectivity for regional scale simulations in which multiple rivers are resolved in a grid cell, it may not be suitable for coarsescale ESM simulations as only one single major river is resolved within a grid cell (Wu et al., 2011).In this work, we implemented two-way hydrological coupling for the land and river components of the Energy Exascale Earth System Model (E3SM; Golaz et al., 2019) to address the following questions: 1. How does two-way coupling impact water and energy cycle compared to one-way coupling?
2. What are the spatial and temporal patterns of the impacts of two-way coupling?Specifically, which regions and periods are more sensitive to two-way coupling?
The default inundation scheme of the river component of E3SM, Model for Scale Adaptive River Transport in (MOSART), uses the elevation profile approach.We developed a novel inundation scheme for MOSART based on a log-linear relationship between the river flow volume and floodplain inundation area.Satellite-based inundation products can be used for model calibration and validation for inundation dynamics due to their growing availability and global coverage (Di Baldassarre et al., 2011;Huang et al., 2014;Papa et al., 2010;Pekel et al., 2016;Prigent et al., 2001;Prigent et al., 2007;Schroeder et al., 2015;Wu et al., 2019).
In section 2, we present a brief description of the land and river components of E3SM, the two-way hydrological coupling scheme, simulation setup, evaluation datasets, model calibration procedure, and the novel inundation scheme.The calibration of river geometry and parameter estimation for the novel inundation scheme are presented in Section 3. Evaluation of the simulated river discharge and inundation fraction is presented in Section 4. In Section 5, the impacts of two-way coupling on land water and energy cycle are presented, followed by the discussion and conclusion in section 6.

Hydrologic processes in E3SM land and river components
E3SM is a fully coupled ESM that includes models for the atmosphere, land, ocean, land ice, sea ice, and river components.A detailed description of the hydrologic and biogeochemical models in E3SM version 1 are provided in Golaz et al. (2019) and Burrows et al. (2020), respectively.In this study, only the E3SM Land Model (ELM) and the river model (MOSART) are used.ELM-v1 (Bisht et al., 2018;Drewniak, 2019;Liang et al., 2019) was developed based on the Community Land Model 4.5 (CLM4.5;Oleson et al., 2013).The hydrology part of ELM is equivalent to CLM4.5 that parameterizes canopy water, snow, surface and sub-surface runoff, soil water dynamics.Readers are referred to Oleson et al. (2013) for a detailed technical description of the hydrological processes in ELM.
MOSART is a physically-based river routing model that simulates transport of water from hillslopes to the river outlet through a subnetwork and a main channel (Li et al., 2013).The current one-way coupled macroscale inundation scheme in MOSART predicts floodplain inundation when the total water volume exceeds the channel storage capacity, and the excess water is transferred from the channel to inundate the floodplain.When the total water volume is less than the main channel storage capacity, there is an exchange of water from the floodplain back to the channel.Given this scheme is one-way coupled, the floodplain water is not lost to the atmosphere through evaporation or to the land through infiltration.The inundation fraction,  !" , is given by: Eq.( 1) where  #$ denotes the volume of water that is in excess of the main channel capacity and  represents the relationship between  #$ and  !" derived from the sub-grid topography (SGT), for example, with the elevation profile (Figure 1a).A detailed description of the default inundation scheme in E3SM, hereafter referred as MOSART-SGT, is provided in Luo et al. (2017).

Two-way hydrological coupling scheme
Similar to the one-way coupling of ELM-MOSART, in the new two-way coupling MOSART receives total runoff from ELM, routes the runoff in the river channel, and simulates floodplain inundation.Unlike the one-way coupling scheme, in the two-way coupling scheme MOSART sends the fraction and volume of inundated water to ELM, which then simulates the infiltration of the inundated water into the soil column.The maximum soil infiltration capacity ( %&!',)*$ ) is used to estimated floodplain infiltration in the two-way coupling scheme, which is formulated as: where  +*, is the fractional saturated area, Θ %-# is an ice impedance factor, and  +*, represents the saturated hydraulic conductivity.ELM calculates the infiltration of inundated water using an approach similar to the infiltration of surface water storage in the soil column (Oleson et al. (2013)) that is given as where  !" represents floodplain inundation volume,  %&!',)*$ is the maximum soil infiltration capacity described in Eq (2),  * is the available capacity in the soil for infiltration in the floodplain, ∆ is the time step,  +*, is the saturated soil moisture in the topsoil layer, and  -./ is the soil moisture.The coupling scheme does not allow floodwater to infiltrate when the soil is saturated.The use of  !" in Eq (3b) ensures floodplain infiltration only occurs in the fraction of soil that is inundated.The total volume of infiltration on the floodplain over the model coupling timestep is sent to MOSART to update the inundation volume at the beginning of the next MOSART time step.

Model setup
ELM and MOSART global simulations were performed at a spatial resolution of 0.5°× 0.5° for 1981-2010.The simulations used the 3-hourly 0.5°× 0.5° Global Soil Wetness Project version 1 (GSWP3v1) atmospheric forcing dataset, which has been dynamically downscaled and bias-corrected based on the reanalysis data (Compo et al., 2011).algorithm (Wu et al., 2012).Land cover and water depth were used to estimate Manning's roughness coefficients for the hillslope, subnetwork, and main channel (Getirana et al., 2012).
The elevation profile for the default inundation scheme in MOSART was developed from the 90m-resolution DEM from Hydrological Data and Maps Based on Shuttle Elevation Derivatives at Multiple Scales (HydroSHEDS; Lehner et al., 2008).
Four sets of simulations were performed in this study, as listed in Table 1, to evaluate the impact of coupling scheme as well as the new inundation scheme.First, a MOSART-only simulation forced by a pre-built 3-hourly 0.05° runoff dataset (DLND-MOSART-1way) was performed to calibrate the river geometry (i.e., channel depth and channel width, see section 2.5).
The pre-built runoff dataset was developed for a long-term global flood analysis called Global Reach-level Flood Reanalysis (Yang et al., 2021) using the Variable Infiltration Capacity Model (VIC) land surface model that is calibrated and bias-corrected (Lin et al., 2019) against machinelearning derived global runoff characteristics (Beck et al., 2015).The GRFR runoff dataset was very well validated against >14,000 river gauges globally (Yang et al., 2021).The calibrated river channel geometry was then used in the next two ELM-MOSART simulations with two-way and one-way coupling, respectively.These simulations used the newly developed inundation scheme (described in section 2.6) to investigate the impact of model coupling on the simulated water and energy cycles.A fourth configuration, ELM-MOSART-SGT-1way, was performed with one-way coupling using the elevation profile based inundation scheme of Luo et al. (2017) and the calibrated parameters of Mao et al. (2019).The fourth simulation was used as a benchmark for evaluating the new inundation scheme.

Evaluation data
In this study, we calibrated and validated the simulated streamflow and inundation using the Global Stream Indices and Metadata (GSIM; Do et al., 2018;Gudmundsson et al., 2018) and the Global Inundation Extent from Multi-Satellites (GIEMS; Papa et al., 2010;Prigent et al., 2001;Prigent et al., 2007;Prigent et al., 2012), respectively.The GSIM dataset includes monthly, seasonal, and yearly streamflow estimated from daily streamflow measurements of ~35,000 gauges worldwide.In this study, we only used the monthly streamflow GSIM data.The GSIM gauges had different temporal coverages within the simulation period.We used the first two-third of the available data during the simulation period at each gauge for model calibration and the remaining one-third of the data for model validation.
GIEMS is a 0.25°× 0.25° monthly inundation dataset based on multiple satellite observations that does not separately identify lakes, reservoirs, and irrigated agriculture from the river inundated areas.The modified GIEMS data (Mao et al., 2019) for only river inundation areas was developed by excluding the water bodies that were identified by the Global Lakes and Wetland Database (GLWD; Lehner and Döll, 2004) and the Monthly Irrigated and Rainfed Crop Areas (MIRCA2000) products (Portmann et al., 2010).The modified GIEMS dataset was upscaled to 0.5°× 0.5° for model calibration and evaluation.Monthly GIEMS from 1993-2002 was selected for model training, and 2003-2007 was used for model evaluation.
We evaluated our model at both global and basin scale.Three basins with different climate characteristics were used in this study to perform evaluation at basin scale: Mackenzie (cold region), Mississippi (subtropical region), and Amazon (tropical region).

Channel geometry calibration
River geometry is a critical factor in river routing models (Yamazaki et al., 2014) and inundation schemes (Decharme et al., 2012).In this study, river channel geometry was calibrated by minimizing errors in the simulated streamflow compared against observed streamflow at the basin scale for 268 major river basins using level 3 basins dataset of Linke et al. (2019).Due to the coarse resolution of MOSART, basins that only have GSIM gauges with contributing area less than 100,000 km 2 were not chosen for use in calibration.Based on this criterion, 59 of 268 basins have at least one qualified GSIM gauge, which in total covers 45% of the land surface, excluding antarctica.When multiple GSIM gauges are present within a basin, the gauge with the largest contributing area was selected.The shape of the river channel cross-section is assumed to be rectangular, and the channel width and depth were calibrated with the following equation as proposed in Andreadis et al. (2013): where  represents the 2-year return period daily streamflow, and  0 and  4 are curve fitting parameters.The 95% confidence intervals for the parameters  0 and  4 are [2.6,20.2] and [0.12, 0.63], respectively.The  for each grid cell is estimated by aggregating daily runoff from the corresponding upstream cells.A set of 5 values for both  0 and  4 were sampled uniformly based on the suggested 95% confidence interval for each parameter, which resulted in a total of 25 parameter sets of river width and depth for the calibration simulations.
The metric of Normalized Root Mean Square Error (NRMSE) was used to evaluate the performance of the model during calibration and is given as , Eq.( 6) where  % and  L % denote the observed and simulated streamflow in the i-th month,  I is the averaged observed streamflow, and  is the total number of months used for model calibration.
The values of  0 and  4 for each basin were selected from the set of 25 parameters that produced the smallest streamflow NRMSE during the calibration period.Basins without a gauge for model calibration were assigned median parameters values (  0 = 7.2 and  4 = 0.27) as proposed by Andreadis et al. (2013).

A novel inundation scheme
Uncertainties and biases in the inundation scheme may lead to larger uncertainties in simulating both land and river processes in the two-way coupled land-river model.A novel inundation scheme is implemented in E3SM to simulate floodplain inundation dynamics more accurately as compared to the default inundation scheme of E3SM.Specifically, there exists a log-linear relationship between the satellite-based inundation fraction and simulated total volume according to our preliminary results (see Text S1).Therefore, we developed a new inundation scheme that estimates the floodplain inundation fraction ( !" ) using the log-linear regression (LLR) when the total volume exceeds the channel capacity ( -9 ) and the relationship is given by Eq.( 7) where  #$ represents the excess volume in the river channel, and  and  are parameters.
MOSART with the new inundation scheme is hereafter referred as MOSART-LLR.Once the inundation fraction is calculated from Eq (7), the excess height (ℎ # ) can be estimated by assuming water depths in the floodplain and over the channel bank top are the same (see Figure 2): Eq. ( 8) Eq. ( 9) where  -9 is the channel fractional area,  is the river width,  is the river length, and  represents the grid cell area.Next, the floodplain inundation volume ( !" ) can be separated from the excess volume with following equation: Eq. ( 10) And the channel volume ( -9 ) is updated by: Eq. ( 11) The procedures to calibrate the slope and intercept parameters  and  in the log-linear relationship of Eq. ( 7) are described below: 1. Estimate the LLR parameters a and b at each grid cell using the DLND-MOSART-1way simulated total volume and the GIEMS inundation fraction.The estimated LLR parameters are used as initial guess for subsequent ELM-MOSART-LLR-2way simulations.
2. Perform ELM-MOSART-LLR-2way simulation using the LLR parameter values obtained from the i-th iteration.

River channel geometry
The monthly scale correlation coefficients between the simulated and observed streamflow at the GSIM stations are greater than 0.6 for all the calibrated basins (Figure 3a), suggesting a satisfactory performance of MOSART in simulating streamflow after river channel geometry calibration.Over some major basins (e.g., Amazon, Mississippi, Yangtze, Mackenzie, etc.), the skill of MOSART is excellent with correlation coefficient greater than 0.9.The calibrated river width (Figure 3a) and depth (Figure 3b) shows similar spatial variability as compared to other global studies (Decharme et al., 2012;Yamazaki et al., 2011).Note that Mao et al. ( 2019) calibrated the river channel geometry for MOSART-SGT by minimizing errors in the simulated basin averaged inundation fraction, resulting in unrealistic shallow river depths for the Amazon and Yangtze River basins, and low spatial variability at global scale (Figure S3).By separating the calibration of river channel geometry from the calibration of the inundation process, we estimate a more realistic river width and depth (Figure 3).

Log-linear inundation parameters
Estimation of the LLR parameters for ELM-MOSART-LLR-2way required four calibration iterations for the NRMSE of average inundation fraction to change by less than 0.01 between the last two iterations.The LLR yielded correlation coefficients between the simulated total volume and observed inundation fraction greater than 0.8 during the calibration period over the wet regions, such as the Amazon basin, lower Mississippi basin, Southern Asia (Figure S2).
The performance of LLR is satisfactory, with correlation coefficients around 0.5 over the high latitude of the Northern Hemisphere.
The default inundation scheme for some grid cells requires unrealistically large floodplain volume that is 1-2 orders of magnitude larger than the channel capacity to yield even the minimum observed inundation fraction of GIEMS (Figure 4).This unrealistic volume-area relationship results from the assumption of elevation profile approach that assumes all the lower elevation locations need to be filled before a higher elevation location is inundated.However, a higher elevation location can be inundated before all lower elevation locations are flooded due to the flow path connectivity in a grid cell at coarse resolution (Figure 1c).The new inundation scheme requires that the floodplain volume needed to produce the minimum observed inundation fraction is of the same magnitude as the channel capacity (Figure 4), implying a more realistic relationship.The excess height in Figure 4 is related to the inundation fraction through Eq (8) in LLR scheme.For the default inundation scheme, the excess height can be estimated by inserting the inundation fraction in the elevation profile.
The floodplain volume and excess height relationship is not always monotonic in LLR scheme (Figure 4b).For example, an increase in the inundation volume can lead to an initial decrease of excess height when the inundation volume flood water spills over the bank to fill the depression or natural levee (from blue line to red line in Figure 2a).This is because a small volume increase leads to a significant increase of inundation fraction.After the inundation volume exceeds the threshold (e.g., the depression is filled up), the excess height increases as the inundation volume increases (red line to green line in Figure 2a).However, the use of elevation profile in the default inundation scheme ignores flow path connectivity, leading to a monotonically increasing relationship between excess height and inundation volume.3.75°, 51.75°), respectively.Excess height is estimated from the inundation volume and inundation fraction using Eq (8) for the new inundation scheme (solid line), and the elevation profile for the default inundation scheme (dashed line).The red squares, blue circles, and green triangles represent the minimum, mean and maximum inundation fraction from the GIEMS dataset, respectively.

Discharge
Using the calibrated river channel geometry and calibrated LLR inundation scheme, ELM-MOSART-LLR-2way shows a reasonably good skill of simulating streamflow seasonality (Figure 5) and interannual variability for 15 selected major basins, with correlation coefficients larger than 0.8 (Table 2).The Congo River basin is an exception showing a lower, but still a satisfactory model performance with correlation coefficient equal to 0.6.The simulated streamflow captures the observed seasonal cycle, but model biases still exist in a few basins (Figure 5).Table 2 summarizes the location of the GSIM stations that were used for evaluation along with other model evaluation metrics.The lower Nash-Sutcliffe efficiency (NSE) that less than 0.5 for Mackenzie, Volga, Lena, Kolyma, Murray-Darling, Irrawaddy, and Congo basin indicate larger biases in the simulated streamflow.Those biases can result from uncertainty in the ELM-generated runoff because of 1) a lack of accounting for water management (Voisin et al., 2013); 2) atmosphere forcing uncertainty (Li et al., 2015); 3) surface and subsurface runoff parameters uncertainty (Huang et al., 2013); and 4) poor representation of snowmelt dynamics (Toure et al., 2018).ELM-MOSART-LLR-1way simulated streamflow seasonality at the 15 major basins is essentially identical to that of ELM-MOSART-LLR-2way (Figure 5), suggesting that two-way coupling does not impact the long-term streamflow seasonality of large basins.While the mean relative change for all months between two-way coupling and one-way coupling is small, the variability of the relative change for certain months in several basins can be large (Figure S4), for example, winter periods for Mackenzie and Godavari.

Inundation
Compared to the default SGT inundation scheme, the new LLR inundation scheme in MOSART significantly improves the simulated inundation dynamics at global scales (Figure 6).
The ELM-MOSART-LLR-2way simulated inundation explains about 96% of the spatial variation of GIEMS during the validation period ( 6 = 0.96 in Figure 6c inset), which is substantially superior to that of the default inundation scheme used in ELM-MOSART-SGT-1way ( 6 = 0.07 in Figure 6b inset).Furthermore, the new inundation scheme is able to capture the temporal variation of the GIEMS inundation with global spatial correlation coefficient larger than 0.7 for each month during the validation period (Figure S5).We acknowledge that Mao et al. (2019) performed model calibration using atmospheric forcing of Qian et al. (2006) instead of GSWP3 that is used in ELM-MOSART-SGT-1way.However, the biases of simulated inundation fraction remains significant when the atmospheric forcing of Qian et al. (2006) is used in ELM-MOSART-SGT-1way (Figure S6).The default inundation scheme underestimates the observed inundation significantly even with an unrealistically shallow river depths, especially over the Amazon rainforest, South and Southeast Asia, and partial high latitude of the Northern Hemisphere (Figure 6b).The model biases are reduced with the new inundation scheme (Figure 5c), though some regions show a slightly overestimation of the inundation.Overall, the averaged global inundated area during the validation period for GIEMS, ELM-MOSART-LLR-2way, and ELM-MOSART-SGT-1way are 2.31 × 10 ?[ 6 ], 2.46 × 10 ?[ 6 ], and 1.54 × 10 ?[ 6 ], respectively.
The new inundation scheme captures the seasonal variation of the basin-averaged inundation fraction better than the default inundation scheme over the Mississippi and Amazon basin (higher correlation in Figure 7a).Although the default inundation scheme shows a better correlation with the GIEMS data at Mackenzie, it underestimates the temporal variance.
Additionally, the default inundation scheme fails to capture the spatial distribution of inundation fraction for Mackenzie, with a low spatial correlation coefficient of 0.05 (Figure 7c).The performance of the default inundation scheme is relatively better for Mississippi and Amazon with spatial correlation coefficients of 0.48 and 0.70, respectively (Figure 7c).The spatial distribution of the inundation fraction is improved with the new inundation scheme significantly, with spatial correlation coefficients higher than 0.97 for all three presented basins (Figure 7d).

Impacts on the water cycle
The two-way hydrological coupling of ELM and MOSART affects the global water cycle.In the two-way coupled simulation, the total infiltration increases as the floodplain inundated water infiltrates in the land during the flooding period (Figure 8a), representing the driver for the changes of other processes.Since the high latitude of the Northern Hemisphere has broader inundation extents (Figure 6a) due to wetter soil (Figure S7), more areas can be affected by two-way coupling through the infiltration from the inundated water.For example, 50% of the inundated cells (i.e.,  !" > 0.01) distributed above 40 o N based on GIEMS inundation dataset.
The increased infiltration leads to an increase in soil moisture (Figure 8b) and a shallower water table (Figure 8c) through soil water movements over the affected areas.The higher soil moisture in the two-way coupled simulations causes larger surface runoff (Figure 8d) as precipitation and snowmelt have less available pore space to infiltrate in the soil.Higher water table also leads to an increased subsurface runoff (Figure 8e).Additionally, the surface water fraction increases with the two-way coupling (Figure 8f) and the increase is mainly distributed in the high latitude of the Northern Hemisphere, where the surface water area is more sensitive to the change of surface hydrological processes due to frozen soil (Avis et al., 2011;Woo and Winter, 1993).In the two-way coupled simulation, the river loses water to land through the floodplain infiltration, but it also receives more runoff (Figure 8d and e), which results in a change of streamflow dynamics.Specifically, two-way coupling reduces the variability of daily discharge by decreasing the maximum daily discharge (Figure 9a), while simultaneously increasing the minimum daily discharge (Figure 9b).The reduction in daily streamflow variability implies that floodplain plays a critical role in mitigating extreme events like flooding.For example, the 30year averaged annual maximum floodplain inundation is lowered in the two-way coupling simulation (Figure 9c).The averaged floodplain inundation also decreases, except in some areas over the Northern high latitude (Figure 9d).On average, the global total floodplain inundation area decreased by 3.3% during our simulation period due to two-way coupling.Table 3 provides a summary of the number of affected cells at the global scale for various variables of interest.At global scale, changes in annual runoff averaged over the grid cells affected by twoway coupling are greater during years with higher annual total runoff (Figure 10a), suggesting an increase of interannual variability.The larger changes of total runoff in two-way coupling simulations are caused by higher infiltration on the floodplain (Figure 10b) due to larger inundation during wetter years (Figure 10c).Indeed, the more inundation water infiltrated into the soil, the more runoff will be generated attributed to a shallower water table depth and higher soil wetness.Additionally, the changes in total runoff are dominated by the changes in subsurface runoff, which account for ~92% of the changes in total runoff.Although a similar pattern for the impacts of two-way coupling on total runoff are found at basin scale (Figure S8), regional differences in the changes of surface and subsurface runoff due to different climate characteristics are noticeable (Figure S9).For example, both the Mississippi River basin and the Amazon River basin have negligible changes of the surface runoff, whereas the Mackenzie River Basin has comparable changes of surface and subsurface runoff (Figure S8).The difference in Mackenzie River Basin is because of high latitude basins are characterized with lower baseflow index (Beck et al., 2013), hence, the surface runoff can be more sensitive to two-way coupling than other areas.Here we find that areas with lower annual runoff (e.g., less than 500 [ ] ⁄ ) have larger changes in infiltration due to two-way coupling (Figure 11a), resulting in larger changes in total runoff (Figure 11b).Although the wetter regions tend to have larger inundation extents than the drier regions, the infiltration capacity may constrain the floodplain infiltration such that the inundated water cannot infiltrate when the top layer soil is saturated (Eq (3b)).Both the number of affected cells and the change of infiltration and runoff decrease along with the increase of annual total runoff.This suggests that relatively drier regions are more sensitive to two-way coupling than wetter regions in terms of the water cycle, therefore, land-river interactions through the floodplain reduces spatial variability of runoff over the inundated areas.

Impacts on energy cycle
The subsurface thermal dynamics in ELM are impacted due to the changes of soil moisture in the two-way coupled simulations.Increased soil moisture with two-way coupling leads to a higher soil heat capacity that consequently modifies the soil temperature (Figure 12a).
The Northern high latitudes show an increase in soil temperature, while other regions with changes in soil temperature show a decrease trend as compared to the one-way coupled simulation.The spatial pattern of surface soil temperature change is attributed to the changes in the annual average soil temperature (Figure 13).Specifically, floodplain inundation results in warmer soil for grid cells with annual soil temperature lower than about 8 °; otherwise, the soil becomes cooler on average.Changes of surface temperature affect the growing season length (GSL), which is defined as the number of days between the first 5-day period with average temperatures above 5°C to the first 5-day period with temperatures below 5°C here.While the average annual soil temperature is warmer in the two-way coupling simulation in cold regions (Figure S10), the GSL can be shorter (Figure 12b).The reason for the shorter GSL with two-way coupling is that more energy is needed to heat up the wetter soil with higher heat capacity during the transition from the freezing season to the thawing season (see May to June of 2001 in Figure 14a).However, the GSL of warmer areas is not affected by two-way coupling as only the hot months are cooled (Figure 14b).
The partitioning of net radiation is modified by two-way coupling as well, with a higher latent heat flux (Figure 12c) and lower sensible heat flux (Figure 12d) compared to the one-way simulation.However, many modifications on energy cycle due to two-way coupling are currently omitted in our implementation, such as energy exchange between land and river, estimating surface albedo with floodplain inundation, and including ET from inundated water.Thus, all the changes in the surface heat fluxes are driven by the changes in soil moisture.Notably, the latent heat fluxes over the tropical regions are not sensitive to the change of soil moisture (Figure 12c) because ET is not water limited for these regions (Brum et al., 2018;Xu et al., 2019).Note: A cell with relative change larger than 0.1% or less than -0.1% are counted as affected.The relative change larger than 5% are counted as significant increase, and less than -5% are counted as significant decrease.For the soil temperature, absolute change is used instead.The cells with absolute change larger than 0.1℃ or less than -0.1℃ are counted as affected.1℃ and -1℃ are used as criteria for significant increase and decrease for the soil temperature, respectively.

Conclusion and discussion
In this study, we developed two-way hydrological coupling between the land and river components of E3SM.The default inundation scheme in the river component of E3SM was inadequate in capturing the observed spatial variability of floodplain inundation, thus, we developed a novel data-driven inundation scheme that is able to capture 96% of the spatial variance of a satellite-based observational dataset.We calibrated river geometry parameters for MOSART and parameters for the new inundation scheme and performed ELM-MOSART coupled simulations with one-way and two-way model coupling to investigate the sensitivity of land/river processes to land-river coupling.Our comparisons reveal significant changes in the land processes at the global scale, with two-way coupling producing wetter soil, more runoff, and higher surface water fraction.Two-way coupling also has impacts on river processes, as evidenced by the lower peak annual streamflow and the higher minimum annual streamflow because of the infiltration of flood water into the floodplain soil and hence, the influence on runoff.While riverine inundation is mitigated by two-way coupling during flooding periods, inland inundation is more extensive.Overall, the water cycle at about 20% of the global areas are influenced by the two-way hydrological coupling.
The water cycle shows different sensitivities to two-way coupling in regions with different climate characteristics.At global scale, our results suggest that wetter periods (e.g., years with runoff higher than the long-term annual runoff) and relatively drier regions (e.g., annual runoff less than 500 [  ⁄ ]) exhibit higher sensitivity to land river two-way coupling.
Larger changes in runoff are observed during wetter periods because the water table and soil moisture are more affected by the infiltration of inundation water, which is larger in the wetter periods than in the drier periods.However, the relatively drier regions are more affected by twoway coupling than wetter regions, where the inundation infiltration may be constrained by the lower infiltration capacity of the soil (Eq (3b)).In contrast, the floodplain soil of relatively drier regions is capable of absorbing most of the inundation water, therefore, leading to larger changes in the water cycle by modifying the daily and interannual variability of streamflow and runoff.
The energy cycle is modified by the land river two-way hydrological coupling as well.
Since only hydrological exchange is implemented between land and river in our implementation, two-way coupling has negligible effect on the surface net radiation (not shown here).
Nonetheless, partitioning of the surface net radiation is different in the two-way coupling simulation with higher latent heat flux and lower sensible heat flux.The higher latent heat flux is mainly driven by the change of soil moisture (Hou et al., 2012), which has an important control on ET in water-limited regions (Jung et al., 2010).We note that ET from open water in river and floodplain inundation, which will lead to more changes in the surface heat fluxes (Dadson et al., 2010), are not included in the current two-way coupling implementation.A future study will include ELM-MOSART simulations with an active atmosphere model to investigate of the impacts of floodplain on land-atmosphere interaction.
The newly developed inundation scheme and land river two-way hydrological coupling are not without shortcomings and limitations.First, the new inundation scheme is not processbased but data driven that relies on accurate flood inundation satellite dataset for training.
Satellite inundation data products have considerable uncertainty associated with detection of small flooded areas (Prigent et al., 2007), presence of clouds (Policelli et al., 2017;Revilla-Romero et al., 2015), and densely vegetated areas (Wu et al., 2019).Second, given that inundation only occurs along the main river channels in MOSART, it is critical to accurately represent the river networks for inundation simulations.Thus, using a single main channel river network (Wu et al. (2011) at relatively coarse resolution (e.g.0.5 degree) would introduce additional source of uncertainty for inundation estimation in this study.Third, the current coupling scheme adds the floodplain infiltrated water within the entire coarse-scale grid cell, which can be unrealistic as the inundated area usually occupies only a very small fraction of the grid cell area.However, infiltration of the inundated water only within a soil column of floodplain is not supported by the current version of ELM, as hydrological processes are represented by a single soil column within each grid cell.A more realistic surface and subsurface interactions in two-way coupling scheme are warranted with the subgrid topographic land unit model setup that will be available in a future version of ELM (Tesfa et al., 2020).Lastly, ELM and MOSART separately estimate inland inundation (i.e., wetland inundation) and floodplain inundation with different inundation schemes.Therefore, it remains challenging to evaluate the total inundations with different components (e.g., land, river), as different inundations often occur concurrently.
In summary, we implemented a land river two-way hydrological coupling scheme in E3SM to understand changes in the land and river processes due to riverine inundation.Our results show considerable impacts of riverine inundation on land and river hydrological processes as well as partitioning of surface energy fluxes, which may have further impacts on the water and energy cycle through land-atmosphere interactions.River and land hydrological processes could be more resilient to climate change as two-way land-river interactions in the floodplains tend to reduce the variability of hydrological processes at global scale, but future investigations are needed using fully coupled E3SM simulations.Lastly, this study provides the necessary first step for representing thermal, sediments, salinity, nutrients exchanges between river and land to better understand the impacts of floodplain inundation on the biogeochemical cycle.

Figure 1 .
Figure 1.(a) The 2-D elevation profile of a grid cell (2.75°, 55.25°) in the Amazon basin used in Luo et al. (2017) derived from the cumulative density function of 90m HydroSHEDS DEM shown in (b).The blue dashed line in (a) denotes the elevation of water needed to inundate 40% of the grid cell.(c) Elevation profiles along the transect shown by the black line between two red cross signs in (b), and the blue shaded areas represent areas below the inundated elevation (blue dashed line) in (b).
Figure2.Conceptual examples of excess height ( ) on floodplain in log-linear inundation scheme.Subplot (a) represents a channel with depressions nearby, and subplot (b) represents a channel with smooth floodplain.The left panels show the reality, while the right panels show how the inundation volume and excess height are modelled in the log-linear inundation scheme.The blue, red, and green solid lines denote the water levels in different flooding scenarios, which result in inundation fraction  !",8 ,  !",6 , and  !",6 , respectively. !",8 ,  !",6 , and  !",5 are the associated inundation volume in the floodplain, and ℎ #,8 , ℎ #,6 , and ℎ #,5 are the corresponding excess height based on the log-linear inundation scheme.

h e Figure 3 .
Figure 3. River channel geometry for (a) river width and (b) river depth calibrated using observed streamflow.The color of the circles in subplot (a) shows the correlation coefficient between the best calibrated simulated streamflow and the observed streamflow from GSIM at each gauge for the validation period.The location of the circle denotes the GSIM gauge location, and the red line delineates the corresponding basin boundary.

Figure 4 .
Figure 4.The relationship between the floodplain excess height and floodplain inundation volume.Subplot (a) and (b) show the relationships from two example grid cells located at (3.75°, 104.75°) and (3.75°, 51.75°), respectively.Excess height is estimated from the inundation volume and inundation fraction using Eq (8) for the new inundation scheme (solid line), and the elevation profile for the default inundation scheme (dashed line).The red squares, blue circles, and green triangles represent the minimum, mean and maximum inundation fraction from the GIEMS dataset, respectively.

Figure 5 .
Figure 5.Comparison of observed and simulated river streamflow seasonality at 15 selected major basins.Streamflow simulated using one-way and two-way coupling are shown in solid red and dashed blue lines.

Figure 7 .
Figure 7. Simulated inundation fraction for Mackenzie (top row), Mississippi (middle row), and Amazon (bottom row) basins.(a) Monthly basin average fractional inundation for the validation period (2003 -2007).(b), (c), and (d) show the spatial distribution of mean inundation fraction for GIEMS, default inundation scheme, and new inundation scheme, respectively.

Figure 10 .
Figure 10.(a) Annual time series of total runoff, surface runoff changes, subsurface runoff changes, and total runoff changes between the ELM-MOSART-LLR-2way simulation (R2way) and the ELM-MOSART-LLR-1way simulation (R1way) averaged over the affected cells at global scale.Annual total runoff is from the ELM-MOSART-LLR-1way simulation.Subplot (b) shows relationship between floodplain infiltration volume and change of total runoff, and (c) shows the relationship between floodplain infiltration volume and fractional inundation from ELM-MOSART-LLR-2way simulation. is the corresponding correlation coefficient.

Figure 11 .
Figure 11.Relationship between the average annual total runoff and (a) change of average annual infiltration, and (b) change of average annual total runoff between the ELM-MOSART-LLR-2way simulation (I2way and R2way) and the ELM-MOSART-LLR-1way simulation (I1way and R1way).Only grid cells that are impacted by the inundation are presented in the density scatter plot.

Figure 12 .
Figure 12.Impacts of land river two-way coupling on (a) soil temperature at 10 cm, (b) growing season length, (c) latent heat flux, and (d) sensible heat flux in ELM.Subplots (a) and (b) show the maximum changes of soil temperature and growing season length from 30 years for each grid cell, respectively, between ELM-MOSART-LLR-2way and ELM-MOSART-LLR-1way.For subplots (c) and (d), the ratio of 30-year means between the ELM-MOSART-LLR-2way simulation and the ELM-MOSART-LLR-1way simulation is used.

Figure 13 .
Figure 13.Relationship between the annual soil temperature at 10cm and the changes of annual soil temperature at 10cm between the ELM-MOSART-LLR-2way simulation and the ELM-MOSART-LLR-1way simulation.Only grid cells that are impacted by the inundation are presented in the density scatter plot.
Lawrence et al.

Table 2 . Evaluation of streamflow over 15 selected major basins.
Note:  is correlation coefficient, NSE represents Nash-Sutcliffe efficiency, and NRMSE is the normalized root mean square error.