Development of Land‐River Two‐Way Hydrologic Coupling for Floodplain Inundation in the Energy Exascale Earth System Model

Floodplain inundation links river and land systems through significant water, sediment, and nutrient exchanges. However, these two‐way interactions between land and river are currently missing in most Earth System Models. In this study, we introduced the two‐way hydrological coupling between the land component, E3SM Land Model, and the river component, Model for Scale Adaptive River Transport, in Energy Exascale Earth System Model (E3SM) to study the impacts of floodplain inundation on land and river processes. We calibrated the river channel geometry and developed a new data‐driven inundation scheme to improve the simulation of inundation dynamics in E3SM. The new inundation scheme captures 96% of the spatial variation of inundation area in a satellite inundation product at global scale, in contrast with 7% when the default inundation scheme of E3SM was used. Global simulations including the new inundation scheme performed at 0.5°×0.5° $0.5{}^{\circ}\times \mathrm{0.5}{}^{\circ}$ resolution with and without land‐river two‐way coupling was used to quantify the impact of coupling. Comparisons show that two‐way coupling modifies the water and energy cycle in 20% of the global land cells. Specifically, riverine inundation is reduced by two‐way coupling, but inland inundation is intensified. Wetter periods are more impacted by the two‐way coupling at the global scale, while regions with different climates exhibit different sensitivities. The two‐way exchange of water between the land and river components of E3SM provides the foundation for enabling two‐way coupling of land‐river sediment and biogeochemical fluxes. These capabilities will be used to improve understanding of the interactions between water and biogeochemical cycles and their response to human perturbations.

The simulation of floodplain inundation in the coarse-scale (∼100 km) 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 river transport 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 sub-kilometer topography data (Figure 1b) by assuming that inundation always occurs from lower to higher elevation. Although some studies have demonstrated satisfactory performance of the sub-grid topography (SGT) schemes in capturing basin-averaged inundation dynamics 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., 10 m) DEM, aims to capture the role of river network . 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 coarse-scale 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 of land and river for floodplain inundation impact water and energy cycles 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?
In this study, we focus on the coupling between land and river in offline coupled simulations, while ET over the floodplain associated with the coupling between river and atmosphere will be considered in future development.
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 10.1029/2021MS002772 3 of 24 availability and global coverage (C. Huang et al., 2014;Di Baldassarre et al., 2011;Papa et al., 2010;Pekel et al., 2016;Prigent et al., 2001Prigent et al., , 2007Schroeder 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 data sets, 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 (H. 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, fp , is given by: where ex denotes the volume of water that is in excess of the main channel capacity and represents the relationship between ex and fp derived from the 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 for Floodplain Inundation
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 ( inf l,max ) is used to estimated floodplain infiltration in the two-way coupling scheme, which is formulated as: where sat is the fractional saturated area, Θice is an ice impedance factor, and sat 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 inf l,fp = min where fp represents floodplain inundation volume, inf l,max is the maximum soil infiltration capacity described in Equation 2, is the available capacity in the soil for infiltration in the floodplain, Δ is the time step, sat is the saturated soil moisture in the topsoil layer, and cur is the soil moisture. The coupling scheme does not allow floodwater to infiltrate when the soil is saturated. The use of fp in Equation 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 data set, which has been dynamically downscaled and bias-corrected based on the reanalysis data (Compo et al., 2011). Lawrence et al. (2019) found better model performance in GSWP3v1-based simulations than simulations driven by other atmospheric forcing data sets. The time step for ELM and MOSART is 30 and 60 min, respectively, and the model coupling time step is 180 min. MOSART uses sub-cycling and the local time step size is chosen to ensure numerical stability. The default 0.5 • × 0.5 • ELM surface data set was used in this study. The topographic parameters (i.e., flow direction, river length, slope, etc.) of MOSART were generated using the Dominant River Tracing 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 90 m-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 data set (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 data set was developed for a long-term global flood analysis called Global Reach-level Flood Reanalysis (GRFR; Yang et al., 2021) using the Variable Infiltration Capacity Model LSM that is calibrated and bias-corrected (Lin et al., 2019) against machine-learning derived global runoff characteristics (Beck et al., 2015). The GRFR runoff data set 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., 2001Prigent et al., , 2007Prigent et al., , 2012, respectively. The GSIM data set 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 data set 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 (Lehner & Döll, 2004) and the Monthly Irrigated and Rainfed Crop Areas (MIRCA2000) products (Portmann et al., 2010). Since the river channel area is included in GIEMS, inundation fraction refers to the sum of floodplain inundation fraction and river fraction. The modified GIEMS data set was upscaled to 0.5 • × 0.5 • for model calibration and evaluation. Monthly GIEMS from 1993 to 2002 was selected for model training, and an independent period of 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 . 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 data set 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 and are curve fitting parameters. The 95% confidence intervals for the parameters and 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 and 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 where and denote the observed and simulated streamflow in the ith month, is the averaged observed streamflow, and is the total number of months used for model calibration. The values of and 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 ( = 7.2 and = 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 and their corresponding impacts in the two-way coupled land-river model. Thus, a novel inundation scheme is implemented in E3SM to simulate floodplain inundation dynamics more accurately than the default inundation scheme of E3SM. Specifically, after logarithmic transformation, the simulated total channel volume exhibits a linear relationship with the GIEMS inundation fraction (Figure 2a). At basin scale, 65% and 50% of areas in the Mississippi and Amazon basins show good correlation (>0.5) between the simulated volume and observed inundation area, but only 10% of the Mackenzie basin has correlation coefficient larger than 0.5 ( Figure 2b). Overall, 30% of global grid cells with an averaged inundation fraction larger than 0.01 have correlation coefficients larger than 0.5 (red line in Figure 2b), although the relationship for less inundated areas (inundation fraction <0.01) is relatively weak (blue line in Figure 2b). Therefore, we developed a new inundation scheme that estimates the floodplain inundation fraction ( fp) using the log-linear regression (LLR) when the total volume exceeds the channel capacity ( ch) and the relationship is given by where ex represents the excess volume in the river channel, and and are tunable parameters. Parameter represents the sensitivity of inundation area to excess volume, such that for larger value of , larger area can be inundated given the same ex . Note that the log-linear volume-area relationship for the floodplain emerged from data fitting. Although the physical justification of the relationship is unclear and should be explored in the future, this study adopts the empirical relationship and compares simulations with the new and default inundation schemes. MOSART with the new inundation scheme is hereafter referred as MOSART-LLR. Once the inundation fraction is calculated from Equation 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 3): where ch is the channel fractional area, is the river width, is the river length, and represents the grid cell area. Next, the floodplain inundation volume ( fp) can be separated from the excess volume with following equation: And the channel volume ( ch ) is updated by: The procedures to calibrate the slope and intercept parameters and in the log-linear relationship of Equation 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 ith iteration. 3. Estimate the i + 1-th LLR parameter values using the total volume from the ith ELM-MOSART-LLR-2way simulation and the GIEMS inundation fraction. 4. Perform ELM-MOSART-LLR-2way with the i + 1-th LLR parameter values. If the averaged NRMSE of the global inundation fraction between the i + 1-th and the ith simulations is greater than 0.01, repeat step 2 to step 4.
The GIEMS inundation data during the training period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002) were used in the above procedures to fit the log-linear relationship, and an independent period of GIEMS (2003)(2004)(2005)(2006)(2007) were used to validate the fitted scheme as described in the following section. The LLR parameter values obtained after the above procedures were also used in the ELM-MOSART-LLR-1way simulation to exclude the impacts of parameter values when evaluating differences between the two-way coupled and one-way coupled simulations.

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 4a), 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 4a) and depth (Figure 4b) shows similar spatial variability as compared to other global studies 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 S1 in Supporting Information S1). 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 4).

Log-Linear Inundation Scheme
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 fitted parameters and are larger over flat areas along the river, where floodplain inundation area is more sensitive to inundation volume ( Figure S2 in Supporting Information S1). 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 5). 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 5), implying a more realistic relationship. The excess height in Figure 4 is related to the inundation fraction through Equation 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 5b). 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 3a). 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 3a). The reasonable magnitude of floodplain volume and the nonmonotonic relationship between the volume and excess height suggest that the connectivity may be implicitly included when the new inundation scheme is trained against observed data. 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.

Discharge
The simulated streamflow is very sensitive to the river channel geometry for the selected major basins ( Figure  S3 in Supporting Information S1). 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 6) 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 6). 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, Danube, Volga, Murray-Darling, Irrawaddy, and Congo basin indicate larger biases in the simulated streamflow. Further, the magnitude of annual averaged simulated streamflow is significantly  Table 2), suggesting the biases are from ELM as routing only affects the shape of hydrograph. Those biases can result from uncertainty in the ELM runoff generation processes because of (a) a lack of accounting for water management (Voisin et al., 2013); (b) atmosphere forcing uncertainty (H.-Y. Li et al., 2015); (c) surface and subsurface runoff parameters uncertainty (M. Huang et al., 2013); and (d) poor representation of snowmelt dynamics (Toure et al., 2018). The evaluation of the simulated streamflow at daily scale can be found in Table S1 in Supporting Information S1, and the performance is similar to the evaluation at monthly scale.
ELM-MOSART-LLR-1way simulated streamflow seasonality at the 15 major basins is essentially identical to that of ELM-MOSART-LLR-2way (Figure 6), suggesting that two-way coupling does not impact the long-term streamflow seasonality of large basins. The minor changes of streamflow seasonality could be caused by the floodplain infiltration being offset by the increase of runoff in the two-way simulation (see more discussion in Section 5.1), and lack of ET from river. 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 in Supporting Information S1), 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 7). The ELM-MOSART-LLR-2way simulated inundation explains about 96% of the spatial variation of GIEMS during the validation period ( 2 = 0.96 in Figure 7c inset), which is substantially superior to that of the default inundation scheme used in ELM-MOSART-SGT-1way ( 2 = 0.07 in Figure 7b 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 in Supporting Information S1). 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 in Supporting Information S1).
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 7b). The model biases are reduced with the new inundation scheme (Figure 7c), though some regions show a slightly overestimation of the inundation. Overall, the averaged global inundated Note. The metrics in the first row and second row are for the evaluation of ELM-MOSART-LLR-2way and ELM-MOSART-LLR-1way simulation, respectively. is correlation coefficient, NSE represents Nash-Sutcliffe efficiency, NRMSE is the normalized root mean square error, and Annual ratio is the ratio of annual averaged streamflow between simulation and observation. All the metrics are estimated with the whole monthly time series of available period. , 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 shown in Table 3).
Although the default inundation scheme shows a better correlation with the GIEMS data at Mackenzie (Table 3), 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 8c). 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 8c). 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 8d). We note the inundation fraction presented in Figure 8 includes the river areas, therefore, the model with both SGT and LLR inundation schemes yield non-zero inundation during the winter months for Mackenzie. The GIEMS data set shows zero inundation for the winter months of Mackenzie, when the river is frozen and covered with snow (Papa et al., 2010). The evaluation of basin averaged inundation for all the other major basins can be found in Table 3, which shows the performance of inundation simulation is much improved with the new inundation scheme.

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 9a), representing the driver for the changes of other processes.
Since the high latitude of the Northern Hemisphere has broader inundation extents (Figure 7a) due to wetter soil ( Figure S7 in Supporting Information S1), 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., fp > 0.01) distributed above 40°N based on GIEMS inundation data set. The increased infiltration leads to an increase in soil moisture ( Figure 9b) and a shallower water table (Figure 9c) through soil water movements over the affected areas. The higher soil moisture in the two-way coupled simulations causes larger surface runoff (Figure 9d) 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 9e). Additionally, the surface water fraction increases with the two-way coupling (Figure 9f) 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 & Winter, 1993).
In the two-way coupled simulation, the river loses water to land through the floodplain infiltration, but it also receives more runoff (Figures 9d and 9e), 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 10a), while simultaneously increasing the minimum daily discharge (Figure 10b). The reduction in daily streamflow variability implies that floodplain plays a critical role in mitigating extreme events like flooding. For example, the 30-year averaged annual maximum floodplain inundation is lowered in the two-way coupling simulation (Figure 10c). The averaged floodplain inundation also decreases, except in some areas over the Northern high latitude (Figure 10d). On average, the global total floodplain inundation area decreased by 3.3% during our simulation period due to two-way coupling. Table 4 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 two-way coupling are greater during years with higher annual total runoff (Figure 11a), 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 11b) due to larger inundation during wetter years (Figure 11c). Indeed, the more inundation water infiltrated into the soil, the more runoff will be generated attributed to a shallower water Note. The metrics in first row and second row are for the evaluation of new inundation scheme with two-way coupling and default inundation scheme, respectively. and Nash-Sutcliffe efficiency (NSE) were estimated with the whole time series during the validation period (2003)(2004)(2005)(2006)(2007), and the spatial correlation is estimated with averaged inundation fraction during the validation period.

Table 3 Evaluation of Basin Averaged Inundation Fraction for the Validation Period
Over 15 Selected Major Basins 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 in Supporting Information S1), regional differences in the changes of surface and subsurface runoff due to different climate characteristics are noticeable ( Figure S9 in Supporting Information S1). 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 in Supporting Information S1). 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.
Two-way coupling reduces daily streamflow variability but increases the interannual runoff variability globally (Figures 10 and 11). The effects of two-way coupling may also differ spatially. Here we find that areas with lower annual runoff (e.g., <500 [mm∕yr]) have larger changes in infiltration due to two-way coupling (Figure 12a), resulting in larger changes in total runoff (Figure 12b). 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 (Equation 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 13a). 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 14). Specifically, floodplain inundation results in warmer soil for grid cells with annual soil temperature lower than about 8 • C ; 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 in Supporting Information S1), the GSL can be shorter (Figure 13b). 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 15a). However, the GSL of warmer areas is not affected by two-way coupling as only the hot months are cooled (Figure 15b). The partitioning of net radiation is modified by two-way coupling as well, with a higher latent heat flux ( Figure 13c) and lower sensible heat flux (Figure 13d) 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 13c) because ET is not water limited for these regions (Brum et al., 2018;Xu et al., 2019).

Conclusion and Discussion
In this study, we developed two-way hydrological coupling between the land and river components of E3SM for modeling floodplain inundation. 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 captures 96% of the spatial variance of a satellite-based observational data set. We calibrated river geometry parameters for MOSART and parameters for the new inundation scheme. ELM-MOSART coupled simulations were performed 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 impacts river processes 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  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°C or less than −0.1°C are counted as affected. 1 • C and −1 • C are used as criteria for significant increase and decrease for the soil temperature, respectively.

Table 4
Percentage of Global Cells That Are Affected by Two-Way Coupling 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 the 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 <500 [mm∕yr] ) exhibit higher sensitivity to the two-way land-river 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 two-way coupling than wetter regions, where the inundation infiltration may be constrained by the lower infiltration capacity of the soil (Equation 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. Figure 11. (a) Annual time series of total runoff, surface runoff changes, subsurface runoff changes, and total runoff changes between the ELM-MOSART-LLR-2way simulation (R 2way ) and the ELM-MOSART-LLR-1way simulation (R 1way ) 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 panel (c) shows the relationship between floodplain infiltration volume and fractional inundation from ELM-MOSART-LLR-2way simulation. is the corresponding correlation coefficient.
The simulated energy cycle is also modified by the land river two-way hydrological coupling. Since only hydrological exchange is implemented between land and river in our implementation, the 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). The evaporation from open water in river and floodplain inundated, which represents the vertical coupling between river/floodplain with the atmosphere, is not included due to the limitation of the current model structure of E3SM. Specifically, the river model does not occupy a fraction of a land grid cell in E3SM and does not exchange fluxes of water and energy with the atmosphere model. However, we found the evaporation from river and floodplain water represents an insignificant part in the energy cycle in the selected basins (Text S1 in Supporting Information S1). Consequently, the simulated changes in the energy Figure 14. Relationship between the annual soil temperature at 10 cm and the changes of annual soil temperature at 10 cm 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. cycle are driven by the lateral water exchanges between land the river. While previous studies argued that evaporation from floodplain is significant for water and energy balances in arid regions such as the Niger delta (Dadson et al., 2010;Getirana et al., 2021) or some high latitude regions , floodplain evaporation could have been overestimated due to model structure uncertainty. For example, floodplain evaporation may be sensitive to interactions between river and atmosphere and the river thermal processes, which could affect the water temperature, but such processes were not explicitly represented in previous studies. Future studies will include floodplain evaporation through three-way hydrologic and thermal coupling of the atmosphere, land, and river components to further investigate the impacts of river and floodplain on the coupled system.
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 process-based but data driven that relies on accurate flood inundation satellite data set for training. However, 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). Another caveat of the new inundation scheme is that its parameters are fitted with GSWP3 in this study, and additional calibration/ fitting is needed when another atmospheric forcing is used. Since the GIEMS inundation data is only available for a short period under current climate condition, the trained inundation scheme may not capture the extreme inundation in the historical and future periods accurately. The sensitivity of the log-linear relationship to forcing data should be investigated for more general applications in ESMs. 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°) 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).
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 land-river two-way 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.