Two-way coupling between the sub-grid land surface and river networks in Earth system models

Over the past decade, there has been appreciable progress towards modeling the water, energy, and carbon cycles at field-scales (10–100 m) over continental to global extents. One such approach, named HydroBlocks, accomplishes this task 10 while maintaining computational efficiency via sub-grid tiles, or Hydrologic Response Units (HRUs), learned via a hierarchical clustering approach from available global high-resolution environmental data. However, until now, there has yet to be a macroscale river routing approach that is able to leverage HydroBlocks’ approach to sub-grid heterogeneity, thus limiting the added value of field-scale land surface modeling in Earth System Models (e.g., riparian zone dynamics, irrigation from surface water, and interactive floodplains). This paper introduces a novel dynamic river routing scheme in HydroBlocks that is 15 intertwined with the modeled field-scale land surface heterogeneity. The primary features of the routing scheme include: 1) the fine-scale river network of each macroscale grid cell’s is derived from very high resolution (<100 m) DEMs; 2) the inlet/outlet reaches of each macroscale grid cell are linked to assemble the continental river networks; 3) the river dynamics are solved at a reach-level via the Kinematic wave assumption of the Saint-Venant equations; 4) a two-way coupling is established between each sub-grid tile and the river network. To implement and test the novel approach, a 1.0-degree bounding 20 box surrounding the Atmospheric Radiation and Measurement (ARM) Southern Great Plains (SGP) site in Northern Oklahoma (United States) is used. The results show: 1) the implementation of the two-way coupling between the land surface and the river network leads to appreciable differences in the simulated spatial heterogeneity of the surface energy balance; 2) a limited number of tiles (~300 per 0.25-degree cell) are required to approximate the fully distributed simulation adequately; 3) the surface energy balance partitioning is sensitive to the river routing model parameters. The resulting routing scheme provides 25 an effective and efficient path forward to enable a two-way coupling between the high-resolution river networks and existing tiling schemes within Earth system models. https://doi.org/10.5194/gmd-2020-291 Preprint. Discussion started: 28 October 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Recent years have seen a renewed effort to improve the representation of land surface heterogeneity in Earth System Models 30 (ESMs) (Chaney et al., 2018;Clark et al., 2015;Fan et al., 2018;Newman et al., 2014). This effort is driven, in part, by the limitations of existing upscaling parameterizations to adequately simulate the multi-scale interactions between the water, energy, and carbon cycles and must instead be modeled more explicitly at higher spatial resolutions (Bierkens et al., 2014;Wood et al., 2011). Although the desired goal might eventually be a modeling approach that models every meter-scale grid cell over the globe, for the foreseeable future, the significant computational barriers will continue to limit the feasibility of this 35 approach. This is especially true given the need for large ensemble frameworks to handle the resulting immense structural and parameter uncertainties that emerge in models when moving to very high spatial resolutions over large spatial extents (Beven et al., 2015). As a result, there is a persistent need to use reduced-order modeling approaches that are able to converge on the desired meter-scale modeling while minimizing computational expense.
One approach to sub-grid heterogeneity widely adopted within ESMs over the past decades is the use of "tiling" schemes, 40 which effectively partition the sub-grid into representative clusters (e.g., forests, grasslands, etc.). This concept is similar to the concept of hydrologic response units in hydrologic models. Since their origin in the late 1980s and early 1990s (Avissar and Pielke, 1989;Koster and Suarez, 1992), tiling schemes have moved to include spatial variability in land use, soil type, elevation gradients, and management practices, among others. Chaney et al., 2016 andNewman et al., 2014 take this a step further by designing a mechanism to derive these tiles by clustering available spatial environmental datasets of the drivers of 45 landscape heterogeneity. Furthermore, over the past decade, there has been an effort to adapt these tiling schemes to enable hillslope-scale processes (e.g., lateral flow) via Darcy flow along topographic gradients within a grid cell (Chaney et al., 2018;Clark et al., 2015;Fan et al., 2018;Subin et al., 2014). These additions are making it possible for ESMs to simulate the role of local lateral flow along topographic gradients which has implications for modeling the riparian ecosystems and surface energy balance partitioning. However, the two-way interconnectivity between the modeled hillslopes and their respective channels 50 has yet to be explored.
In parallel to the development of sub-grid tiling schemes, there has been significant innovation in river routing schemes within ESMs in recent years. For a comprehensive review of the advances in macroscale routing schemes the reader is referred to (Shaad and Di Baldassarre, 2018); a brief overview of these advances is provided here. Early developments in these schemes focused on the horizontal transfer time and integration of runoff over the land surface, allowing for the validation of early Land 55 Surface Model (LSM) outputs against streamflow observations. This initial work is embodied in the Impulse Response Function Unit Hydrograph (IRF-UH) approach described in (Lohmann et al., 1996), which represented the process as linear and time-invariant, both in grid cells and between grid cells in the river network. Alternative approaches were proposed, focusing instead on storage-based schemes to represent the spatial distribution of flow. The Total Runoff Integration Pathways (TRIP) model presented by (Oki and Sud, 1998) was the first to use a simplified form of the 1D Kinematic Wave routing equation with a constant velocity for water movement between grid cells, an approach that has since been implemented in many LSMs over varying scales of global river networks (Pappenberger et al., 2010). With recent computational and data advances, these basic models have been expanded to better utilize the outputs of general ESMs for global predictions beyond just discharge. For example, the recent MizuRoute model (Mizukami et al., 2016) builds upon the approaches described above.
MizuRoute combines both approaches to model the flow from hillslope to channels through a gamma-distribution based IRF-65 UH, while the main channel routing is performed using either a kinematic wave approximation or a linear IRF-UH method, enabling the use of higher resolution terrain and network data with relatively coarse ESM outputs. Alternative approaches for using high resolution vector river networks have been proposed. The Routing Application for Parallel Computation of Discharge model (RAPID) is one such model; it leverages advances in high performance computing to create a highly efficient vector river routing scheme that uses a matrix-based version of the Muskingum-Cunge method (David et al., 2011(David et al., , 2016. 70 Recent macroscale routing model development has sought to make use of advances in data and computational resources to further enable integrations with ESMs. The Model for Scale Adaptive River Transport (MOSART; Li et al., 2013), a recent example of this, divides a macroscale grid cell into hillslopes, tributaries, and a main channel and solves the flow of water through the river network via the Kinematic wave approximation of the St. Venant equations. This decoupling of the ESM grid from that of the routing model (network), allows for differences in scale to more accurately capture flow and channel 75 dynamics at any point within a given basin. Other models have attempted to move beyond the more simplistic representations of flow hydrodynamics to better represent flood inundation. A recent example of this is the Catchment-based Macro-scale Floodplain model (CaMa-Flood) (Yamazaki et al., 2011), which discretizes the macroscale grid cells into unit catchments and explicitly parameterizes the topography of the floodplain for each unit. CaMA-Flood then simulates open channel flow via the local inertial equation, a simplified form of the full 1-D Saint Venant equation, allowing for an estimation of water depth and 80 area of inundation in ESM grid cells.
Although there have been significant advances in river routing and sub-grid heterogeneity schemes over the past decade, there has yet to be a concerted effort to couple these two concepts in ESMs. This persistent lack of interconnectivity between the modeled river network and the land surface leads to: 1) oversimplified sub-grid river networks (e.g., tributaries are mostly ignored) (Jones and Schmidt, 2017;Rice, 2017;Swanson and Meyer, 2014); 2) a lack of interaction between the sub-grid land 85 surface tiles and the sub-grid river networks (e.g., simulated flooding of the Nile River doesn't recharge the land surface) (Bisht et al., 2017;Helton et al., 2012;Shen et al., 2016); 3) the water moving through the sub-grid river network does not influence the surface energy partitioning (i.e., macroscale schemes continue to mostly treat river networks as pipes) (Bisht et al., 2017;Sheng et al., 2017;Zampieri et al., 2011); 4) sub-grid irrigation and water management schemes are spatially agnostic and rely almost exclusively on the main channel (Pokhrel et al., 2015;Shaad and Di Baldassarre, 2018;Voisin et al., 2017). 90 All these deficiencies compound to illustrate a persistent weakness in Earth system models by ignoring critical processes that are known to play an important role in both the natural and engineered hydrologic systems (Fan et al., 2018;Pokhrel et al., 2016). With the recent advances in both river routing and sub-grid tiling schemes, the timing is right to enable two-way interactions between the sub-grid tiling schemes and river networks in ESMs. This study provides a path to address this persistent deficiency in ESMs. This is accomplished by implementing a reach-95 based routing model in the HydroBlocks land surface model (Chaney et al., 2016) and enabling a two-way coupling with the modeled hydrologic response units. The primary features of the novel river routing scheme include: 1) each macroscale grid cell is assigned its own field-scale river network delineated from DEMs; 2) the fine-scale inlet/outlet reaches of the macroscale grid cells are linked to assemble the continental river networks (and to ensure conservation of mass); 3) river dynamics are solved at the reach-level via an implicit solution of the Kinematic wave simplification of the St. Venant equations; 4) a two 100 way coupling is established between each cell's sub-grid tiles and the river network. The scheme is implemented over a 1.0degree bounding box around the Southern Great Plains site in Northern Oklahoma in the United States. Furthermore, a series of experiments are performed to understand: 1) the sensitivity of the land surface to the two-way coupling with the sub-grid river network; 2) the number of tiles that are required to adequately represent the sub-grid heterogeneity; 3) the impact of the uncertainty in the routing scheme parameters in the macroscale response. 105

Study Domain
A 1.0-degree box in central northern Oklahoma and southern Kansas in the United States is used to implement and test the new routing scheme (see Figure 1). The region is generally flat with a shallow decreasing gradient in terrain, precipitation, and NDVI from west to east. The overall climate is generally dry with cold winters and hot summers. The Salt Fork Arkansas river, 110 Chikaskia, and to a lesser extent the Arkansas river traverse the region; all rivers in the domain eventually flow into the Arkansas river. The vegetation throughout the domain is primarily croplands with forested regions found along the riparian areas. There are also small urban areas dispersed throughout the region with the town of Enid in the southeast corner being the largest. This domain contains the Atmospheric Radiation and Measurement (ARM) Southern Great Plains central facility (SGP CF) among other ARM SGP facilities. As the largest climate facility in the world, ARM SGP collects a wealth of data on land-115 atmosphere interactions and atmospheric processes that are used frequently to evaluate and improve sub-grid atmospheric processes in ESMs. Understanding the role that the interconnected land surface and river network play in the partitioning of the surface energy balance (and its role in land-atmosphere interactions) is a primary motivation for using this domain in this study.

Land surface model: HydroBlocks 120
HydroBlocks is a field-scale resolving land surface model (Chaney et al., 2016) that accounts for the water, energy, and carbon balance to solve land surface processes at high spatial and temporal resolutions. HydroBlocks leverages the repeating patterns that exist over the landscape (i.e., the spatial organization) by clustering areas of assumed similar hydrologic behaviour into HRUs. The simulation of these HRUs and their spatial interactions allows the modeling of hydrological, geophysical, and biophysical processes at the field-scale (30 m) over regional to continental extents (Chaney et al., 2016b;Vergopolan et al., 2020). The core of HydroBlocks is the Noah-MP vertical land surface scheme (Niu et al., 2011). HydroBlocks applies Noah-MP in an HRU framework to explicitly represent the spatial heterogeneity of surface processes down to field scale. At each time step, the land surface scheme updates the hydrological states at each HRU; and the HRUs dynamically interact laterally via subsurface flow. In the original HydroBlocks, subsurface flow between HRUs was modeled via a subsurface kinematic wave. In this new implementation, following the approach described in Chaney et al. (2018), the subsurface flows have been 130 updated by computing the Darcy flux between adjacent HRUs at each subsurface level. These terms were included as divergences within the vertical one-dimensional solution of Richards' equation in Noah-MP. This allows for the flow of water between HRUs to also be driven by capillarity and not just the predefined topographic gradient.

Hierarchical generation of hydrologic response units
The HRU generation scheme in the original HydroBlocks (Chaney et al., 2016) was not sufficient to adequately capture the 135 dynamics in riparian zones (e.g., runoff), which led to the development of the hierarchical multivariate clustering (HMC) approach introduced in Chaney et al. (2018). To properly define the riparian zone dynamics and to enable adequate coupling with the river network, this paper builds on HMC. The hierarchical scheme consists of: 1) Characteristic basins -The basins are first delineated from a 30-meter DEM. To avoid basins being split between macroscale grid cells, the fixed regular grid is adjusted such that basins are assigned to the macroscale grid cell where it 140 has the majority of its area (see Figure 2A for an example). These basins are then clustered using K-means to obtain a set of k characteristic basins. Figure 2B shows an example of this division for k = 10 for one of the 0.25 degree grid cells in the sub-domain. These characteristic basins were identified using latitude, longitude, flow accumulation area, and the natural logarithm of the flow accumulations area as feature predictors. Note that the concept of a characteristic basin is only used for the non-routing component of HydroBlocks. After clustering the basins, to maximize similarity among the 145 basins that belong to a given characteristic basin, the empirical cumulative distribution function (CDF) of the 30-meter height above nearest drainage (HAND; Nobre et al., 2011) values within each basin is matched to the characteristic basin's HAND empirical CDF. The characteristic basin's HAND empirical CDF is set to be the average of all of its basins' empirical CDFs (for a background on CDF matching see Reichle and Koster, (2004)).
2) Height bands -The 30-meter HAND values of each characteristic basin are discretized into height bands. First, all channel 150 grid cells within a given characteristic basin are grouped into one height band. Then, the non-channel grid cells are discretized into additional height bands by binning HAND values. The binning involves creating groups of HAND values that have an areal coverage n (user-defined) times larger than its adjacent lower height band. For example, the height band immediately above the channel height band will have an areal coverage n times larger than the channel height band. The next height band will have an areal coverage n 2 larger than that of the channel and so forth. Figure 2C shows the 155 discretization into height bands of the basins of one of the macroscale grid cells in the domain. An additional maxhb parameter is used to ensure the number of height bands per characteristic basin does not exceed this user-defined threshold (maxhb is fixed at 100 in this study).
3) Intra-band clusters -Each height band within each characteristic basin is further divided into clusters (i.e., HRUs) to represent intra-band heterogeneity of land use, soils, and elevation, among others. This is accomplished by iterating per 160 height band and clustering all the corresponding 30-meter grid cells into intra-band clusters. Given that the areal coverage of each height band can vary; the number of intra-band clusters per height band is set to be proportional to their fractional coverage of the characteristic basin. The user-defined p parameter is the average number of clusters per height band. Figure 2D shows an example of this division for p = 5 using latitude, longitude, land cover, and clay as covariate features.

River routing scheme 165
The river routing scheme implemented within HydroBlocks is a reach-based implicit solution of the St. Venant equations with the Kinematic wave assumption. The conservation of mass across a 1-dimensional reach is given by: !# is the change in the cross-sectional area of the flow with respect to time ) is the change in discharge across an infinitesimally small reach length (note that = where is the velocity of the flow), and S is a source/sink term of the 170 cross-sectional area. Following the Kinematic wave assumption, uniform flow is assumed (i.e., ( = ) ), and thus u can be estimated at each time step for a given reach via Manning's equation: Where is Manning's coefficient, * is the hydraulic radius, and ( is the slope of a given reach. Note that the hydraulic radius is given by where is the wetted perimeter. For this implementation, a compound channel composed of a 175 rectangular channel and a symmetric floodplain is assumed (see Figure 3). The segmented conveyance method is used where the conveyance of both channel, c, and floodplain, f, components are summed and the effective velocity for the compound channel is computed.
As shown in Figure 3, each reach's cross-sectional profile is composed of a rectangular channel and a symmetric floodplain 180 that is learned from the reach's discretized height bands (computed in Section 2.3). At each reach, the cross-sectional area of https://doi.org/10.5194/gmd-2020-291 Preprint. Discussion started: 28 October 2020 c Author(s) 2020. CC BY 4.0 License. the flow and the derived cross-sectional profile are used to determine the wetted perimeters / and ) and the cross-sectional areas / and ) at each time step.
To solve the St. Venant equation with the Kinematic wave assumption at each time step, a fully implicit first-order finite volume upwind scheme in space and backward Euler in time is used: 185 where the upstream 3 3 contribution to a given reach is the sum of the discharge of all the reaches that are immediately upstream of a given reach ( 3 3 = ∑ 3 3 ). In other words, the upstream internal boundary of reach i is the sum of the discharge of all reaches j that flow into reach i (Jacovkis and Tabak, 1996). Note that is the length of a given reach in the network and will rarely be the same across reaches. S is a source/sink term that can account for inflow from runoff as well as 190 water abstractions from irrigation and recharge of the land surface. This leads to a system of non-linear equations (per macroscale grid cell) to be solved per time step. Given the non-linearities in the hydraulic radius calculation as well as the need to iteratively exchange network boundary conditions between macroscale grid cells, an iterative time step is implemented.
In the current implementation, the Picard iteration approach is used to attempt convergence at each time step.
Similar to Mizukami et al. (2016), an impulse response function is used to route the runoff produced at the HRU level 195 in HydroBlocks to its corresponding channel. HydroBlocks' HRU delineation is first used to assemble a histogram of travel times (constant velocity of the flow per grid cell) of all 30-meter grid cells that belong to a given HRU to their closest channel.
The travel times consider a constant velocity of the flow per grid cell (fixed to 0.1 m/s in this study). The path of water flow in the landscape is computed via d8 flow direction. At each time step, the convolution of this histogram (~ unit hydrograph) with the HRU runoff is used to compute the runoff from that time step that reaches the channel at each future time step. 200

Assembling the stream network and reach profiles
For each macroscale grid cell, the river network is derived at a 30-meter spatial resolution following the method described in Chaney et al. (2018). For each river reach, the channel width and bankfull depth are computed using the functional relationships derived for the Contiguous United States described in (Bieger et al., 2015). To minimize inconsistencies between reaches that belong to the same characteristic basin, the computed channel width and bankfull depths are averaged across all corresponding 205 reaches. The cross-sectional profile for each reach (not characteristic basin) is learned by first extracting all the fine-scale grid cells in a given basin; the profile follows the height band discretization from Section 2.3. The precomputed channel length is used to compute each section of the profile's width, given the corresponding height band's area within that reach. Because the spatial resolution of the DEM (30 meters) is much larger than many of the computed channel widths of the delineated streams (~1 meter), the additional fractions of the channels' collocated grid cell that are not part of the channel are assigned to the 210 height band right above the channel (i.e., the lowest element of the floodplain). Figure 3 shows an example of the computed profile for reach 122 in grid cell 7 in the domain. To avoid "step changes", the shown linear interpolation between the hand values of each height band is used when computing the cross-sectional area ) and wetted perimeter ) of a given channel's floodplain (as shown in Figure 3). A different cross-sectional profile is made per reach in a given macroscale grid cell. It is important to note that the routing scheme is run on the complete DEM-derived sub-grid network of river reaches and not the 215 characteristic basins (defined in Section 2.3).

Two-way interaction: Coupling of the HRUs and reaches
This section explains how the sub-grid HRUs interact each time step with the routing scheme's river reaches. After the river routing scheme (i.e., Kinematic wave) is updated for a given time step ( Figure 4A), each reach's derived cross-section profile (see Section 2.5) is used to determine the inundation height over each of the height bands of that reach. This is done by 220 effectively pouring the volume of water contained within the reach and calculating the inundation heights that correspond to a flat surface at the top of the river. To assemble the inundation heights for a given characteristic basin, the inundation heights per height band are averaged across all basins that belong to a given characteristic basin ( Figure 5B); these computed inundation heights are then equally distributed to their corresponding HRUs ( Figure 5C). The computed inundation height per HRU is then added to the following time step's update of Noah-MP as a constant flux at the soil surface ( Figure 5D). Through 225 this approach, this water is not only able to infiltrate the soil in the HRU via Noah-MP, but it is also able to move laterally to another HRU via HydroBlocks' modeling of subsurface flow (see Section 2.2 for more details). The total column of runoff produced during the time step at the given HRU is then set as the HRU's new inundation height. The difference between the old and new inundation heights are computed ( Figure 5E) and averaged up to the height band level ( Figure 5F). To determine the changes of inundation height at each reach (not characteristic basin), the computed differences in inundation heights at the 230 characteristic basin level are scaled to the reach level. This is accomplished by multiplying the change in inundation height value by the ratio between the basin's original inundation height per height band and its corresponding characteristic basin inundation height. These differences are then aggregated into a difference in cross-sectional area Δ per basin ( Figure 5G). Finally, this difference is divided by the model time step and added to the reach's source/sink term in the routing scheme ( Figure 5H). Note that the source/sink term also includes runoff from the land surface that is not inundated; this water arrives 235 at the reach via the method described in Section 2.4. After the update of the cross-sectional area per reach via the river routing scheme, the two-coupling process starts anew ( Figure 5A). 2. Sensitivity to the two-way coupling -To evaluate the impact of the two-way interaction between the modelled rivers and the land surface, two model simulations are performed. For both simulations, the HRU generation scheme is the same as that used for the baseline experiment. The first one uses the default approach of not allowing the routing scheme to interact 260 with the land surface while the second enables the two-way connectivity scheme introduced in Section 2.6.

Model experiments
3. Convergence -Similar to the approach used in Chaney et al. (2018), nine different model experiments are run to evaluate how the HRU configuration parameters impact the model simulations. Simulations A-C focus on inter-basin heterogeneity by increasing k from 1 to 5 to 10, while setting p = 1 and n = 1000. Simulations D-F focus on the discretization of the height bands by decreasing n from 5 to 3 to 2, while setting p = 1 and k = 10. Finally, simulations G-I focus on the role of 265 intra-band heterogeneity by increasing the average number of intra-band clusters from 2 to 3 to 5, while setting n = 2 and p = 5.
4. Sensitivity to routing model parameters -To determine a first order understanding of the sensitivity of the updated HydroBlocks to the routing model parameters, a 25-member parameter ensemble is run. More specifically, the routing model parameters that are used include Mannings coefficients for the channel and floodplains ( and ), as well as the 270 channel width ( ), and the bankfull depth ( ). Given that each reach has its own parameter values, to minimize complexity, a set of scalar multiplier parameters are used to uniformly scale the precomputed parameters of all reaches.

Exploratory simulation
As an initial baseline simulation, the updated HydroBlocks model is run between 2015 and 2017 at an hourly resolution over the study domain using the HRU generation scheme parameters defined in Section 2.7. To provide a general understanding of the simulated land surface, routing scheme, and floodplain dynamics, this section focuses on a simulated inundation event that occurs during the first half of August. Figure  To further investigate the impact that this flooding event has on the land surface, Figure 6 shows the mapped simulated 290 inundation height, root zone soil moisture, and latent heat flux over the same four central 0.25-degree grid cells on August 5 th 14:00 UTC, August 11 th 14:00 UTC, and August 14 th 14:00 UTC. The chosen time steps coincide with a dry period before the simulated event, right after the rain event, and a few days after the event. A brief description of the simulated states and fluxes for the three days is shown below.
August 5 th , 2017 14:00 UTC -All rivers within the domain are within their banks with relatively low stage height 295 within the channels. Not surprisingly, the root zone soil moisture is very low throughout the entire area, with the notable exceptions being the Salt Fork Arkansas River (west to east) and the Chikaksia River (northeast) and lakes throughout the region. The tributaries throughout the region also have higher soil moisture due mostly to recharge from redistribution of runoff from the land surface via the river network. A similar story is evident in the mapped latent heat flux with the only noticeable difference being the higher latent heat fluxes in the south of the domain. 300 August 11 th , 2017 14:00 UTC -This time step falls within the flooding event. Several tributaries of the Salt Fork Arkansas River are flooding, and the stage height in the Salt Fork Arkansas River is appreciably higher than that only 6 six days before. The signal of the flooding is immediately apparent in the root zone soil moisture where the channel and adjacent HRUs are close to or at saturation. These differences are not as noticeable in the latent heat flux most likely since the entire area is at or close to field capacity due to the widespread rain event and thus the evapotranspiration is not constrained by soil 305 moisture.
August 15 th , 2017 14:00 UTC -The floodwaters have receded, and all rivers are again within their banks, although the stage height of the Salt Fork Arkansas River is still higher when compared to August 8 th . Except for the channel HRUs, root zone soil moisture has decreased when compared to August 11 th . However, the imprint of the flood event is still evident in the riparian zones, which to a lesser extent can also be seen as identified in the latent heat flux. 310 The simulated event in Figure 5 and 6 provides a first order understanding of how all the pieces of the implemented land surface and routing scheme come together to directly impact surface fluxes, soil moisture, and flooding over the domain.

Sensitivity to the two-way coupling
Two different model experiments were run to further investigate the role that including a two-way coupling between the routing scheme and the land surface has on the modeled states and fluxes. The first one called "uncoupled" uses the default approach 315 (i.e., not allowing the routing scheme to interact with the land surface). The second simulation is called "coupled" and it enables the connectivity scheme introduced in Section 2.6. Figure 7 shows the difference in simulated annual mean sensible heat flux, latent heat flux, root zone soil moisture, and land surface temperature over the four central macroscale grid cells in the domain. For each variable, the results are shown for both the uncoupled (i.e., the routing scheme does not interact with the land surface) and the coupled simulation. These 320 simulations would suggest that over this region, the local lateral flow of subsurface flow is only one of the contributors to sustaining the riparian zones (as represented in HydroBlocks). The redistribution of water via surface flow, flooding, and recharge also plays a key role in maintaining these ecosystems. The comparison between the annual mean makes this more apparent in the root zone soil moisture, which then appears in the latent heat flux, sensible heat flux, and skin temperature. Figure 8 investigates the role of this coupling further by plotting the time series of the spatial mean and spatial standard 325 deviation of sensible heat flux, inundation height, latent heat flux, land surface temperature, and root zone soil moisture. In general, the results show how the coupling can lead to differences in both the spatial means and spatial standard deviations.
The differences are more appreciable for the spatial standard deviation, as would be expected from the spatial maps from Figure 7. While the relative change in the spatial mean between the two is on the order of 0.1-2 %, the relative change in the spatial standard deviation can be around 20-50% for variables like latent heat flux and sensible heat flux. The differences are 330 more extreme during the summer when the region is drier and, thus, a recharge mechanism for the riparian zones will play a much more important role (as seen in Section 3.1). Overall, the results appear to show that adding a two-way coupling between the land surface and routing scheme will increase the macroscale evaporative fraction (decrease Bowen ratio).

Convergence
To determine how well the reduced order model representation in HydroBlocks can represent the heterogeneity of the domain, 335 a convergence analysis is performed over the study domain. More specifically, the model is run for increasingly complex HRU configurations. Simulations A-C focus on inter-basin heterogeneity by increasing k from 1 to 5 to 10 while setting p = 1 and n = 1000. Simulations D-F focus on the discretization of the height bands by decreasing n from 5 to 3 to 2 while setting p = 1 https://doi.org/10.5194/gmd-2020-291 Preprint. Discussion started: 28 October 2020 c Author(s) 2020. CC BY 4.0 License. and k = 10. Finally, simulations G-I focus on the role of intra-band heterogeneity by increasing the average number of intraband clusters from 2 to 3 to 5 while setting n = 2 and p = 5. 340 Figure 9 illustrates how the increase in heterogeneity complexity impacts the fine scale simulated sensible heat flux.
More specifically, it shows the annual mean sensible heat flux over the domain. The explanation of the results below is split into understanding the role of each parameter: A. Increasing the number of characteristic basins (increasing k; Figure 9A-C) -Initially, there are two HRUs represented per 0.25-degree grid cell; one channel and one "floodplain". In this scenario, upon recharge of the floodplain, the characteristic 345 inundation heights are effectively an area-weighted average between all of the inundation heights of all reaches (remember that the flow through each reach is still resolved in the routing scheme). Increasing the number of characteristic basins from 5 to 10 starts to show differences in the connectivity between the rivers and the land surface. Increasing the number of characteristic basins leads to a separation between main channels and tributaries and, thus, a distinction between their interaction with their respective riparian zones. Increasing the number of characteristic basins also leads to an increase in 350 sub-grid heterogeneity due to the ability to represent spatial differences in land cover and soil properties. For example, with 10 characteristic basins the southwest grid cell is able to start to represent the town of Enid.
B. Increasing the number of height bands (decreasing n; Figure 9D-F) -As the n parameter decreases (ratio between the area of a height band and its adjacent height band below it), the number of height band increases. This enables a finer discretization at the channel/floodplain interface, which in turn leads to more realistic cross-section profiles per reach and 355 thus improved floodplain dynamics. When the number of height bands is too low, the inundation height that should correspond only to the region immediately adjacent to the channel is instead evenly distributed to a much larger area upslope, thus, diffusing its influence on recharge of the riparian area.
C. Increasing the number of intra-band clusters (increasing p; Figure 9G-I) -The increase of intra-band clusters leads to a substantial increase in the heterogeneity within each 0.25-degree macroscale grid cell. The most noticeable difference is 360 the emergence of the urban areas and the country roads and interstates (although these are not clearly visible at 1.0 degree).
The other changes are driven by the separation between crops, grasslands, forests, and bare soil as well as soil properties.
Note how, in this case, the boundaries between the macroscale grid cells effectively disappear. Figure 10 formalizes the convergence analysis by showing how the temporal mean of the spatial mean and standard deviation of the plotted 30-meter maps per 0.25 degree varies as a function of the number of HRUs. For all cases, the largest 365 differences in spatial means are across the different macroscale grid cells, which illustrates the controling role of climate and the model land characteristics in the grid cell mean values for the majority of the variables. This can be seen also in Figure 9, where the west to east gradient in precipitation/temperature/vegetation is readily apparent in the modeled sensible heat flux.
For each case, the largest differences in the spatial mean occur when increasing the number of characteristic basins. For the spatial standard deviation (which can be interpreted as a metric of heterogeneity), the changes are more abrupt, with the largest 370 changes occurring when increasing the number of characteristic basins. For all grid cells in the domain, the convergence is relatively quick (although one could argue that root zone soil moisture has not yet converged). It is encouraging that all macroscale grid cells follow similar paths and that the use of 300-350 HRUs appears to be sufficient to adequately model the fine-scale features while maintaining computational efficiency-the 16 interconnected cells take 5 minutes per year of simulation on 16 cores. This setup coincides with k = 10, n = 2, and p = 5, which is the reasoning behind its use throughout 375 the paper. It should be noted that this convergence analysis is non-exhaustive since it only looks at the predefined HMC parameter path configuration (see Section 2.7). For example, the role of the number of characteristic basins might be different if one starts by increasing the number of intra-band clusters instead of the number of characteristic basins. Work is ongoing to find the optimal path configuration to minimize the number of HRUs even further.

Parameter sensitivity 380
To assess the role that the routing module parameters' plays in the two-way coupling between the rivers and the land surface, a 25 member Latin Hypercube Sample is run. This sensitivity experiment allows to perform a preliminary (yet non-exhaustive) analysis of the role that the Mannings coefficients (channel and floodplain), the bankfull depth, and the channel width have on the simulated macroscale states and fluxes. Figure 11 shows the range of differences between each ensemble member and the baseline uncoupled simulation that was used in Section 3.2. These results focus exclusively on July of 2017 since it covers a 385 time period that shows a rapid wetting of the soil, followed by a prolonged drying period. The biggest difference between the uncoupled and coupled simulations happens primarily during the dry down period, where the latent heat spatial mean and spatial standard deviation are significantly enhanced with respect to the baseline simulation. It is important to note that for all cases, even though there is a clear signal in the sensitivity of the coupling to the parameters, the general role that the river plays in recharging the local area is consistent for most parameters sets. This indicates that the first order effect of river recharge 390 over this landscape is fairly insensitive to plausible model parameters. The parameter sensitivity plays an important role during the dry periods, with the flooding and recharge components playing a key role in these differences. As expected, the spread in the surface fluxes is most pronounced during the middle of the day since this is when available soil moisture will play the largest role in transpiration.

Coupling the river network and the land surface in Earth system models
The existing lack of interconnectivity between the modeled river network and the land surface has been a long-ignored process in ESMs. The results from this study aim to show a path forward while maintaining the computational efficiency of existing land surface schemes in ESMs. That being said, there are still a number of additional challenges that need to be addressed to make this approach suitable for ESMs: 400 Reducing the number of sub-grid reaches -Although the number of reaches per 0.25 degree grid cell in this study (~400 reaches per cell) remains computationally tractable, this will not always be the case. This is especially true in models where the number of sub-grid tiles is usually kept to under a few dozen; and thus, the routing component would consume a significant portion of the land model computation. Furthermore, the increase in the number of reaches for coarser macroscale grid cells would quickly add a significant computational burden (e.g., 2,500-20,000 reaches per 1.0-degree grid cell). Therefore, moving 405 forward, there remains a need to simplify the sub-grid networks. One approach being explored by the co-authors is to cluster the lower stream orders. For example, given that the first order reaches usually take up more than half of the reaches, if these were clustered and the numerical solver was adapted to handle this framework, this would reduce the number of reaches significantly. This approach could be applied to all tributaries while maintaining the sub-grid reaches of the main channels that traverse the grid cell. 410 Reducing the number of iterations -Another weakness of the implemented approach is the need to update the boundary conditions iteratively per time step to ensure convergence. This can lead to multiple data exchanges across the HPC network per time step. Although appropriate for an implementation where the number of land model tiles/HRUs is sufficient to still take up most of the computation, this could quickly become prohibitive for ESMs. To handle this challenge, one could run the routing scheme at a higher temporal resolution; however, this would not really solve the problem, as the number of iterations 415 would be replaced by more time steps. The most obvious solution is an adaptive time time-stepping approach. However, this would also cause synchronization problems over large domains as there would often be a region where due to a flood event a large number more of iterations would be necessary, thus slowing down the entire model update for that time step.

Including water management in HydroBlocks
One of the primary motivations for the implementation of a routing scheme coupled to the land surface in HydroBlocks is to 420 enable the implementation of water management in HydroBlocks (e.g., irrigation from surface water abstraction, reservoir operations, etc). Although additional work is necessary, the reach-based routing scheme provides the scaffolding to allow for a robust incorporation of water management within HydroBlocks. Reservoir implementation can now be readily implemented by fixing the output flows of the corresponding reach by the reservoirs' operation rules. The flooding component of the scheme will then enable the valley to fill-up and, thus, producing a first-order representation of the time-varying reservoir spatial 425 coverage. The reach based implicit method also allows for surface water abstraction where each HRU that has irrigation can draw water from either the most accessible reservoir (as defined by the closest reach that is categorized as a reservoir), or directly from the closest channel.

Improving the modeling of channel and floodplain dynamics
Although the routing scheme implemented in HydroBlocks is able to represent the role of floodplains in the compound channel 430 flow (see Section 2.4), there is still a lack of proper representation of the separation between channel and floodplain flow.
Moving forward, an alternative would be to adapt the extended Saint-Venant equations used in HEC-RAS (Brunner, 2010).
This approach solves the flow on both the channel and the floodplain by coupling two solutions of the Saint-Venant equations via flow between the channel and floodplain.
In addition, the implemented routing scheme uses the kinematic wave approximation of the Saint-Venant equations. 435 However, this assumption is not suitable for many flat regions (e.g., Mississippi, Amazon). The simplest path forward appears to be the implementation of a diffusive assumption of the Saint-Venant equations. The added value of this approach would extend beyond a more realistic movement of water through the channel. It would also be useful lakes and reservoirs are included in HydroBlocks as backwater effects will need to be represented.
Finally, the numerical scheme currently implemented is an implicit finite volume upwind scheme, which can cause 440 significant numerical diffusion at large time steps. Although this helps account for floodwave diffusion in the numerical solution, this is simply an artifact of the numerical method and not an actual represented process. Moving forward, replacing the current scheme with a higher-order implicit scheme should help address this problem.

Preparing continental land surface modeling for meter-scale elevation data
The United States Geographical Survey (USGS) is currently working to collect LIDAR data over the entire contiguous United 445 States (CONUS) to assemble a 1-meter elevation dataset over CONUS in the coming decade. This significant increase in the spatial detail and quality of elevation data will provide an improved representation of small channel networks and flood control infrastructure. However, the nature of fully distributed models, makes the use of these data mostly prohibitive in the short term, especially if we wish to account for the two-way coupling explored in this study. The clustering approaches explored using HydroBlocks provide a path forward to be able to quickly make use of these data; the HydroBlocks tiling approach seeks 450 to leverage hydrologic similarity to represent what matters over the landscape while minimizing redundancies.

Conclusions
The existing lack of interconnectivity between the modeled river network and the land surface in Earth system models leads to: 1) oversimplified sub-grid river networks (e.g., tributaries are mostly ignored); 2) no existing mechanism for sub-grid river networks to interact with the sub-grid land surface tiles; 3) water moving through the sub-grid river network does not influence 455 the surface energy partitioning; 4) sub-grid irrigation and water management schemes are spatially agnostic and rely mainly on the main channel. This study illustrates a path forward to address this persistent deficiency in Earth system models. This is accomplished by implementing a reach-based routing model in HydroBlocks and enabling a two-way coupling with the modeled hydrologic response units. The primary features of the novel river routing scheme include: 1) each macroscale grid cell is assigned its own field-scale river network delineated from DEMs; 2) the fine-scale inlet/outlet reaches of the macroscale 460 grid cells are linked to assemble the continental river networks; 3) river dynamics are solved at the reach-level via an implicit solution of the Kinematic wave simplification of the St. Venant equations; 4) a two-way coupling is established between each cell's sub-grid tiles and the river network. The experiments run over the study domain illustrate the sensitivity of the land surface to recharge due to floodplain dynamics and the role that it can have on surface energy partitioning. Furthermore, the         clustering algorithm. Simulations A-C increase k from 1 to 5 to 10 while setting p = 1 and n = 1000; simulations D-F decrease n from 5 to 3 to 2 while setting p = 1 and k = 10; simulations G-I increase the average number of intra-band clusters from 2 to 3 to 5 while setting n = 2 and p = 5. The average number of HRUs per 0.25-degree grid cell is shown as the title of each panel (the number can differ per grid cell).