Representing Intrahillslope Lateral Subsurface Flow in the Community Land Model

The concept of using representative hillslopes to simulate hydrologically similar areas of a catchment has been incorporated in many hydrologic models but few Earth system models. Here we describe a configuration of the Community Land Model version 5 in which each grid cell is decomposed into one or more multicolumn hillslopes. Within each hillslope, the intercolumn connectivity is specified, and the lateral saturated subsurface flow from each column is passed to its downslope neighbor. We first apply the model to simulate a headwater catchment and assess the results against runoff and evapotranspiration flux measurements. By redistributing soil water within the catchment, the model is able to reproduce the observed difference between evapotranspiration in the upland and lowland portions of the catchment. Next, global simulations based on hypothetical hillslope geomorphic parameters are used to show the model's sensitivity to differences in hillslope shape and discretization. Differences in evapotranspiration between upland and lowland hillslope columns are found to be largest in arid and semiarid regions, while humid tropical and high‐latitude regions show limited evapotranspiration increases in lowlands relative to uplands.


Introduction
The Community Earth System Model (CESM) simulates the evolution of the physical, chemical, and biological components of Earth's environment (Hurrell et al., 2013), and typical simulations have global spatial coverage spanning decades-to centuries-long time periods. With current computational resources, CESM is generally operated at spatial resolutions between 25 to 200 km (Bacmeister et al., 2018;Kay et al., 2015). Even at the lower end of this spatial range (e.g., 25 km), significant subgrid-scale heterogeneity exists in atmospheric, oceanic, and terrestrial conditions. Over land, such heterogeneity is represented within the Community Land Model version 5 (CLM) through the use of a tiling approach, in which the land surface is decomposed into tiles having different properties Oleson et al., 2013). In this approach, each grid cell may be comprised of subareas having different land use and land cover and different vegetation characteristics. However, topographic gradients within a grid cell are not explicitly considered, and the tiles are not hydrologically connected.
As reviewed by Fan et al. (2019), hydrological connectivity and the movement of water along topographic gradients are primary controls of the distribution of soil water within the landscape and therefore strongly linked to observed patterns of vegetation. For example, the convergence of water in topographic lows can support vegetation during the dry season, in contrast to upland areas that dry out more rapidly and experience a reduction in plant productivity. The spatial patterns of vegetation in turn influence climate both locally and globally (Ahlström et al., 2015;Bonan, 2008;Duveiller et al., 2018).
The hillslope is perhaps the fundamental spatial scale at which lateral redistribution of water occurs. At these scales (tens to thousands of meters), elevation differences between ridges and valleys cause downslope water movement and the accumulation of water in the lower areas of the landscape. However, explicitly resolving these spatial scales globally is not computationally feasible at present. Recent studies have used the concept of hydrologic similarity to show that it may be possible to capture the salient features of hillslope-scale hydrologic processes using reduced complexity models (Newman et al., 2014). In particular, representative hillslopes provide a means of explicitly including intrahillslope processes driving lateral water movement (Ajami et al., 2016;Chaney et al., 2018;Subin et al., 2014). In this formulation, topographic data are used to decompose a larger area (e.g., a grid cell or catchment) into hillslopes, which are then analyzed to identify one or more representative hillslopes that capture the average structure of "typical" hillslopes. The representative hillslope is divided into connected elements (e.g., toeslope, midslope, and upslope), and the land model is then run for each element, passing water between elements according to the interelement hydraulic gradient. The surface fluxes of energy, moisture, and carbon at the larger scale are then calculated as an area-weighted average of the fluxes from the hillslope elements.
In this study, we incorporate the representative hillslope concept into CLM. The hillslope configuration is first used to simulate the water balance of a well-monitored first-order catchment in which the dominant vegetation types cluster in spatially distinct areas within the catchment. The catchment is represented as a single hillslope comprised of four spatial units, and the model is then used to relate the observed vegetation spatial distribution to differences in soil water conditions among the hillslope units. By varying the intrahillslope connectivity, the impact of the lateral exchange of water between hillslope units on soil moisture is quantified.
After validating the CLM hillslope configuration at the catchment scale, we next examine how interactions between climate and hillslope form influence lateral water exchange and the distribution of soil water conditions within a hillslope. Two hypothetical hillslope profiles are contrasted, and two methods of discretizing the hillslope are compared to explore the sensitivity of the soil moisture simulation to structural differences in the CLM hillslope configuration. Finally, the impact of the differences in soil water state on evapotranspiration are described.

The CLM
The CLM version 5  is the land component of the CESM (Hurrell et al., 2013). CLM simulates the partitioning of moisture and energy from the atmosphere, the redistribution of moisture and energy within the land surface, and the export of fresh water and heat to the oceans. Biogeophysical processes represented in CLM include solar and longwave radiation interactions with the vegetation canopy and soil, turbulent fluxes of momentum, heat, and moisture from canopy and soil, heat and moisture transfer through soil and snow, and stomatal physiology and photosynthesis. Some of the hydrological processes included are interception of precipitation by the vegetation canopy, throughfall, infiltration into soils, surface and subsurface runoff, snow and soil moisture evolution, evaporation from soil, and vegetation and transpiration Oleson et al., 2013). Spatial land surface heterogeneity is represented as a nested subgrid hierarchy in which grid cells are composed of multiple land units; the CLM land units are vegetated, lake, urban, glacier, and crop. Each land unit can have a different number of columns, and each column can have multiple patches each with a specific plant functional type (PFT) or crop functional type.
The subsurface in CLM is vertically discretized into a series of layers whose number and thicknesses are prescribed. This vertical structure is the same for all columns of a land unit type. By default, CLM5 uses a 25-layer soil column spanning a depth of 50 m. While all layers are thermodynamically active, not all layers are hydrologically active. The hydrologically active region of the soil column ranges from 0.4-to 8.5-m depth and varies spatially. Layers below the hydrologically active region have thermal properties representative of bedrock and zero hydraulic conductivity. Within the hydrologically active region, the vertical redistribution of soil moisture is governed by Darcy's law. When saturated layers exist within the soil column, a water table is diagnosed from the vertical soil moisture profile. Surface runoff is parameterized as a function of water table depth. Subsurface lateral runoff (also referred to as baseflow) is parameterized as a function of the thickness of the saturated zone.
Surface and subsurface runoff generated by CLM are passed to the Model for Scale Adaptive River Transport (MOSART, Li et al., 2013) that routes the runoff across the land surface. MOSART simulates the movement of water across the landscape using a hierarchy of spatial units: the hillslope, the tributary stream network, and the main river channel. Surface runoff is routed across the hillslope, after which, it is combined with subsurface runoff and passed to the tributary spatial unit. The tributary spatial unit uses a lumped routing scheme having a transport capacity equal to the total channel capacity of the tributaries within the grid cell. The main channel receives water from both the tributary spatial unit and main channel discharge from upstream grid cells and discharges water downstream. MOSART uses Manning's equation to calculate the flow velocity. The hillslope and tributary spatial units utilize a kinematic formulation in which the friction slope is approximated by the topographic slope. The main channel flow velocity can be calculated using the kinematic or diffusive wave approximation.

Representing Hillslope-Scale Hydrological Processes in CLM
CLM is typically operated using a single column to represent the hydrologic state of the vegetated land unit. To represent hillslopes in CLM, we introduce multiple columns within the vegetated land unit. A hillslope consists of one or more hydrologically connected columns, and a land unit may contain more than one hillslope. The number of hillslopes and columns in each land unit is prescribed by an input data set and may vary spatially. The input data set also prescribes the connectivity between columns by indicating each column's downslope neighbor. Multiple columns may be connected to the same downslope column. Surface runoff is passed from each column to the routing model MOSART, where it is combined and routed as overland flow across the MOSART hillslope spatial unit. Within CLM, subsurface lateral flow is passed between neighboring columns and may occur in the upslope direction if the hydraulic gradient so indicates. At least one column is designated as the "lowland" column, and this column does not have a downslope neighbor but instead discharges its subsurface runoff directly to the MOSART tributary spatial unit. The depth of water in the tributary channel, which provides the boundary condition for the lowland column, is calculated from the volume of water stored in the tributary subnetwork by dividing the volume by the channel width and length, which are time invariant input parameters. Discharge from the lowland column to the tributary network occurs when a saturated zone exists in the soil profile and the elevation of the upper surface of the saturated zone (i.e., the water table) is greater than the elevation of the tributary water surface or the channel bottom if the tributary is dry. If the depth of water in the tributary channel is greater than 0 and the hydraulic gradient between the channel and the lowland column is positive towards the lowland column, water will flow from the channel into the lowland column. When multiple hillslopes are defined, they do not interact directly but indirect communication occurs via the shared boundary condition determined by the tributary channel depth. Intergrid cell lateral flow is not enabled. Figure 1 shows some examples of different hillslope configurations. The simplest approach is to use a single column, which is analogous to standard CLM operation. Alternatively, a hillslope can be modeled as a group of columns connected in series. In this approach, which might be applied to a higher-order stream reach, lateral flow originating in an upland column must pass through all intermediate columns to reach the stream. A hillslope may also be divided into multiple upland columns that are connected in parallel to a single lowland column. This approach might represent a headwater catchment in which strong spatial gradients in meteorological forcing occur. Finally, serial and parallel connectivities can be combined to represent more complicated hillslope flow paths. Column connectivity is determined from an input data set and remains fixed throughout the simulation.

Column Geomorphic Parameters
Each column within the hillslope is described with six geomorphic parameters: the column area A, the mean column elevation h, the width of the column at its downslope interface w, the mean distance of the column from the channel d, the mean slope of the column , and the column's aspect (defined with respect to north) . Figure 2 (left panel) shows these parameters for an individual column.
The geomorphic parameters can be determined from terrain analysis of a digital elevation model and are passed to the model via an input data file. Figure 2 (right panel) shows an example of a single hillslope comprised of four columns connected in series with a stream reach. In addition to having independent geomorphic parameters, each column's soil thickness may vary independently.

Subsurface Lateral Flow
It is assumed that saturated subsurface flow in soils along the hillslope profile is described by Darcy's law.
where q is the moisture flux (m 3 ·m −2 ·s −1 ), K s is the saturated hydraulic conductivity (m/s), Ψ is the hydraulic head (m), and x is the distance from the channel reach (m). In unconfined saturated flow, the hydraulic head is just the water table height, h.
Volumetric discharge Q (m 3 s −1 ) across a plane normal to x is obtained by multiplying by the saturated area.
where w is the width of the column (m) and D s is the thickness of the saturated zone (m). When the elevation gradient is much larger than the change in saturated thickness along the flow path, the head gradient can be approximated by the topographic gradient. With this assumption, called the kinematic wave approximation (Beven, 1981), equation (2) becomes where is the topographic slope. When the kinematic wave approximation is enabled, the interaction between the stream and the lowland column is unidirectional, that is, no flow may occur from the stream to the lowland column. Whether to calculate discharge using equations (2) or (3) is specified when a CLM simulation is configured and is the same for all grid cells for the duration of the simulation.
The soil water states in CLM are updated using an operator splitting approach (Clark & Kavetski, 2010), in which the Richards equation is solved first using a zero vertical flux lower boundary condition at the soil-bedrock interface, followed by the calculation of the change in storage in any soil layers that are saturated due to lateral subsurface flow. In the hillslope configuration of CLM, subsurface flow is passed to the column's downhill neighbor; in contrast, subsurface flow in the standard CLM configuration is passed directly to the river routing model. Thus, a column's net subsurface drainage (i.e., the difference between lateral subsurface inflows and outflows) can be both positive and negative in the hillslope configuration, with the latter case indicating moisture convergence in the column.

Reynolds Mountain East Catchment Application
Using one or more representative hillslopes to approximate the hydrologic and thermodynamic state of a large region such as a model grid cell is a strategy for improving the fidelity of the model while minimizing the computational cost (Ajami et al., 2016;Clark et al., 2015;Chaney et al., 2018;Newman et al., 2014). An additional feature of this model configuration is the ability to assess the simulation by comparing to observations made at the spatial scale of individual catchments. To demonstrate the ability of CLM to perform multicolumn, hillslope-scale simulations, we apply the model to simulate a single, well-monitored catchment.

Catchment Description
Reynolds Mountain East (RME) catchment is located in the headwaters of the Reynolds Creek Experimental Watershed in Idaho, United States (Flerchinger et al., 2010;Seyfried et al., 2009). The RME catchment has an area of 38 ha, with vegetation comprised of approximately 68% shrubland (mainly sagebrush) and 32% forest (mainly aspen; Marks et al., 2002). About 70% of mean annual precipitation occurs as snow, and the spatial distribution of the snowpack varies due to prevailing wind patterns and sheltering by vegetation and topographic features . The catchment experiences dry summer conditions, with July and August typically receiving the lowest precipitation amounts of the year (Hanson, 2001). The location of the stream channel with respect to the catchment elevation is shown in Figure 3, along with the topographic slope.

Site Measurements
Hourly meteorological data (downwelling solar radiation, temperature, humidity, wind speed, and precipitation rate; Reba et al., 2011) for the period 1984 to 2016 (water year) are available. These data are used as atmospheric boundary conditions for the CLM simulations. Validation data include snow water equivalent (SWE) surveys , catchment discharge at a weir located at the northern outlet (Pierson et al., 2001), and three eddy covariance (EC) sites (Flerchinger et al., 2010). The EC sites were Note. Units are meters, except for A that has units of 10 4 m 2 and , which is unitless. Column 1 is the lowland column. Subsurface lateral flow from each of the upland columns (Columns 2-4) is passed to Column 1. Abbreviation: RME = Reynolds Mountain East.

Model Configuration
The RME catchment was decomposed into a lowland (riparian) column connected in parallel to three upland columns (Figure 3, right panel) based on topography, vegetation distribution, and wind-sheltering regime described in Winstral and Marks (2002). The standard tiling approach used by CLM allows for multiple PFTs to exist within a grid cell, but there is no spatial relationship between the PFTs. Because the columns within a hillslope do have an explicit spatial relationship, it is possible to account for the covariance between vegetation type and location within the hillslope. Based on field surveys  indicating that the primary vegetation types were species of sagebrush and aspen trees, PFTs were specified to be broad-leaf deciduous tree in Column 1 and broad-leaf deciduous shrub in the three upland columns. Soil  thicknesses of all columns were specified to be 2 m based on Seyfried et al. (2009). All other model parameters were default values used by CLM. The hillslope physiographic parameters, derived from analysis of a 10-m digital elevation model, are given in Table 1.
In addition to facilitating intercolumn lateral subsurface flow, the CLM hillslope configuration can be used to apply spatially variable meteorological inputs and surface properties within a grid cell. Because of exposure or shelter by topography and vegetation, wind-exposed areas in RME become scoured and sheltered areas develop snow drifts, leading to approximately 50% lower snow in exposed areas and 150-350% greater snow in areas where drifting occurs . To account for the effects of spatially varying wind exposure on snow deposition, snowfall was redistributed within the catchment keeping the total snow amount fixed. The amount of snowfall in Column 2 was doubled, snowfall in Column 4 was halved, and Columns 1 and 3 received unmodified snowfall rates. Rainfall rates were assumed to be uniform across the catchment. Figure 4 shows time series of SWE (top panel) for the period 1988 to 2008. The observations are shown in black, while the CLM grid cell average values are shown in red. Because the snowfall was explicitly modified across the hillslope, the snowpack varies between columns; the intercolumn spread of the CLM SWE is indicated by gray shading. Averaged over the time period, observed and CLM (grid cell average) annual maximum SWEs are 497 mm and 468 mm, respectively, with a root-mean-square difference of 74 mm or about 15%. Stream discharge for the period 1988 to 2008 is shown in the bottom panel of Figure 4. The timing and magnitude of the CLM stream discharge agrees well the observed discharge; however, CLM shows more discharge events during winter than observed. Annual mean discharge is 5.8 × 10 −3 and 6.0 × 10 −3 m 3 /s for the observed and CLM time series, respectively, with a root-mean-square difference of 7.8 × 10 −4 m 3 /s or about 13%. Column values represent net runoff, while the grid cell value represents outflow from the hillslope. Bottom panel: warm season (May-September) accumulated ET. "Os" indicates observed ET from the "Sagebrush" site, and "Oa" indicates observed ET from the "Aspen" site. Figure 5 compares evapotranspiration from the two EC sites in RME to evapotranspiration (ET) simulated by CLM. Latent heat flux measurements are converted to ET by dividing by the latent heat of vaporization (λ v = 2.5 MJ/kg). The three upper panels of Figure 5 compare the measured ET flux from the "Sagebrush" site to the three upland columns' ET during the years 2004-2007. The simulated ET fluxes generally reproduce the observed fluxes, with annual peak fluxes between 100 and 150 mm/month. The effect of the different dates of snowmelt for the three upland columns can be seen in the timing of the ET time series. Column 2, which represents the snow-drifted portion of the catchment and therefore has the highest peak SWE and latest melt date, shows a delayed onset of ET. Column 4, which represents the wind-scoured portion of the catchment, exhibits an earlier and more gradual increase in ET during the growing season. The remaining column received the mean grid cell snowfall amounts, and its behavior lies between the other two upland columns. The bottom panel of Figure 5 compares the "Aspen" site ET measurements to the lowland column ET for the years 2007 and 2008. Like the measured flux, the annual cycle of the CLM ET flux for the lowland column is larger in magnitude than the upland column, peaking between 200 and 250 mm/month. Figure 6 summarizes the comparison of the RME measurements of SWE, stream discharge, and ET flux to the CLM-simulated quantities during the period 2004-2008. The top panel of Figure 6 compares each column's maximum SWE and the grid cell mean maximum SWE to the observed value. The middle panel compares the net runoff of each column (i.e., outflows minus inflows) to the grid cell mean outflow and the observed stream discharge. While the CLM grid cell mean outflow will always be positive, the net runoff of the individual columns will be negative when runoff convergence occurs (i.e., inflow greater than outflow). The bottom panel compares the integrated warm season (May to September) evapotranspiration of each column to the EC measurements from the "Sagebrush" and "Aspen" sites.

Model Assessment
The total ET from the lowland column (Column 1) is greater than the total ET of any of the upland columns due to its higher summer soil water conditions. The left panels of Figure 7 show vertical profiles of CLM soil moisture as a function of time for the years 2004-2008. Because of its relatively low snowpack, the wind-scoured column (Column 4) receives less infiltration and is therefore drier than the other columns; saturated conditions are not present in Column 4 during this time period. Columns 2 and 3 are wetter than Column 4, with saturated conditions forming during the spring snowmelt season. Despite receiving earlier and less snowmelt than the wind-drifted column (Column 2), the lowland column has higher soil moisture conditions later in the summer season due to the additional water gained from subsurface lateral flow from the upland columns. The right panels of Figure 7 show the net lateral flow from each column. Because of its lack of saturated conditions, Column 4 produces no subsurface flow. Columns 2 and 3 both produce flow, which is passed to Column 1. The net lateral flow of Column 1 is sometimes positive but mainly shows negative values, indicating that the outflow from the lowland column is less than the inflow. This convergence of water increases the soil moisture that supplies the greater evapotranspiration flux from the lowland column.

Kinematic Wave and No-Lateral Flow Simulations
The higher soil moisture values in the lowland column relative to the upland columns shown in Figure 7 are the result of both subsurface lateral flow convergence and drainage from the stream to the lowland column. To examine the relative impact of these processes, we contrast the Darcy flow simulation with a simulation using the kinematic wave approximation and a simulation in which lateral subsurface flow from every column is passed directly to the stream rather than to another column. The latter behavior is analogous to the manner in which runoff is treated in the standard (single column) CLM.
When the kinematic wave approximation (equation 3) is used, the hydraulic gradient depends only on topographic slope and not on the water table gradient between columns. Because of this, lateral flow between columns may occur only in the downslope direction, and water in the stream channel may not infiltrate back to the lowland column as the lowland column's water table declines. In the no-lateral flow simulation, moisture convergence between columns is eliminated such that the soil moisture of the lowland column is only supplied by local infiltration. Figure 8 compares the lowland column soil moisture for the three simulations. In the no-lateral flow simulation, soil moisture is on average the lowest of the three simulations. The soil profile experiences saturated or near-saturated conditions briefly in the spring but dries relatively rapidly. The kinematic wave simulation shows wetter conditions for a longer period of time, indicating that the additional inputs of water from upland columns are a significant source of moisture for the lowland column. The Darcy flow simulation has somewhat wetter conditions than the kinematic wave simulation because the stream acts to buffer surface runoff within the catchment, which may reinfiltrate and maintain the relatively high soil moisture values in the lowland column.
The differences in soil moisture values during the summer shown in Figure 8 are reflected in the respective ET fluxes. Figure 9 compares the mean annual cycle of ET for the three simulations. The Darcy flow simulation increases more rapidly in the spring and maintains the highest values during summer. The ET fluxes in the kinematic wave simulation are slightly lower but also extend longer into the summer than do the ET fluxes from the no-lateral flow simulation. During spring, when infiltration is relatively high due to snowmelt, the no-lateral flow simulation has ET fluxes that are almost as high as those of the kinematic wave simulation but because it depends entirely on local moisture inputs, the no-lateral flow ET fluxes decline more rapidly during summer. Fan et al. (2019) postulated that the impact of lateral groundwater redistribution on evapotranspiration fluxes varies regionally due to the interaction of terrain and climate. For example, in regions experiencing seasonal drought, upland water surpluses during the wet season may subsidize lowland water deficits during the dry season. Conversely, in arid regions where no surpluses are generated or in wet regions where no deficits occur, the presence or impact of lateral flow may be minimal. In this section, we examine the effect of including intrahillslope lateral subsurface flow in CLM under a variety of climate conditions. Global simulations were performed using hypothetical hillslope profiles. All hillslopes are specified to have uniform channel-to-ridge distances and widths of 500 m. Two profiles were selected: one linear and one nonlinear. At each grid cell, the hillslope profile was scaled by the standard deviation of the topographic elevation such that the height difference between the valley and ridge was greater in areas of the world with greater topographic variability. The height at any point on the hillslope for on the linear hillslope profile is given by

Global Simulations
where h max is the elevation of the top of the hillslope relative to the stream channel (m), x is the distance from the stream channel reach (m), and L is the total distance from the stream to the ridge (m).
The height at any point on the hillslope based on the nonlinear hillslope profile is given by where and k and n are shape parameters; in this study, k = 5 and n = 1.6. Figure 10a shows the linear hillslope profile scaled to a maximum height of 50 m and Figure 10d shows the nonlinear profile, similarly scaled.
The manner in which a hillslope is discretized into columns can impact the simulation. Of particular significance is the elevation of the lowland column (i.e., the column that interacts with the stream). The vertical reference point is the stream bankfull height D bankfull ; and the stream bed thus is at a depth D bankfull below z = 0, while the surface of the lowland column is at a height h 1 above z = 0. If the thickness of the soil profile is D soil , then the base of the soil profile is z base = h 1 − D soil . If z base is greater than 0, water movement from the channel into the lowland column cannot take place and stream-soil water interactions are unidirectional.
To illustrate the effect of hillslope discretization choice, two binning strategies are contrasted. The first uses equal height bins to divide the hillslope into columns, and the second uses a fixed lower height bin spanning 0-2 m, while the upper three bins are 25%, 50%, and 25% of the remaining height range, respectively. Making the first height bin relatively narrow ensures that the soil depth of the lowland column overlaps with the stream channel and therefore allows possible bidirectional stream-soil water interactions. Figures 10b and  10e show the results of applying equal height bins to the profiles in Figures 10a and 10d. For the linear profile, this results in columns of equal width in the across-channel direction and therefore also of equal area. For the nonlinear profile, the areas are unequal. Figures 10c and 10f show the results of applying the varying height bins to the profiles in Figures 10a and 10d. In this case, the areas of the columns vary across the hillslope for both profiles.

Model Configuration
The global CLM simulations were forced with atmospheric boundary conditions for the period 1991-2010 from the Global Soil Wetness Project 3 (GSWP3). GSWP3 is a 3-hourly 0.5 • global forcing product  that was developed for the third phase of GSWP. It is based on the Twentieth Century Reanalysis version 2 performed with the National Centers for Environmental Prediction (NCEP) model (Compo et al., 2011). The reanalysis was dynamically downscaled to T248 (0.5 • ) resolution and bias corrected for temperature, precipitation, longwave radiation, and shortwave radiation. A wind-induced undercatch correction was also applied. The spatial resolution of the simulations is 1.25 • longitude × 0.9 • latitude. A constant soil thickness of 8 m was specified globally.
For the global simulations performed for this study, two PFTs are specified for each grid cell: a lowland PFT and an upland PFT. All upland columns are assigned the same PFT for this study, but this is not required. Figure 11 shows the distribution of PFT for the lowland and upland columns. Larger PFTs (i.e., trees) are typically associated with the lowland column, while smaller PFTs (i.e., grasses) are generally associated with the upland columns.
Because the PFTs are distributed across a hillslope based on column type (lowland vs. upland), the PFT area will also depend on the hillslope delineation method. The elevation of the top of the hillslope varies spatially according to where elevation is the standard deviation of the elevation within each grid cell calculated from a high-resolution digital elevation model; values of h max used by the global simulations are shown in Figure 12. Figure 13 compares the mean annual water table depth for two CLM simulations. The first (linear, equal height bins ["LIN_EQ"]) uses the linear hillslope profile discretized using equal height bins (Figure 10b). The second (linear, varying height bins ["LIN_VAR"]) uses the same linear hillslope profile but is discretized using the varying height bins shown in Figure 10c. The top row of Figure 13 shows the mean annual water table depths of the highest column (columns are numbered from 1 = lowland to 4 = upland) for the LIN_EQ (left panel) and LIN_VAR (right panel) simulations. Note that when no saturated layers exist, CLM assigns the water table a value equal to the depth of the bottom of the soil profile, but this does indicate that saturated conditions exist. Saturated conditions are indicated only when the water table depth is less than the soil thickness depth. Both maps in the top row of Figure 13 show shallower water table depths in humid areas of low relief and deeper water tables in regions that are dry or have high topographic variability.

Linear Hillslope Profile
All columns within a grid cell experience the same atmospheric boundary conditions but because of the lateral movement of water between columns, some lowland columns have higher water storage and therefore shallower water tables. The lower rows of Figure 13 show the differences in water table depths moving downslope from the highest column. Negative values indicate shallower water table depths in the downslope columns and thus wetter conditions. In the LIN_EQ simulation, water convergence from ridge to valley can be seen in some areas, such as southeastern Asia, eastern North America, and northern Europe. Few differences between columns can be seen in arid and semiarid parts of the world because saturated conditions rarely exist and therefore subsurface lateral flow is rarely generated. In Amazonia, differences in water table depth are not observed because the hillslopes in this region are already near saturation.
The right panels of Figure 13 show the intercolumn differences for LIN_VAR. While water convergence in humid areas is similar for the two simulations, the LIN_VAR simulation exhibits shallower water table depths in the lowland column in arid and semiarid regions. In this simulation, the use of a relatively small height bin (0-2 m) to define the lowland column results in a soil profile that extends to greater depth than the stream channel bottom depth and therefore allows bidirectional exchange of water between the soil column and the stream. During the early part of the the LIN_VAR simulation, this results in the lowland soil column gaining water at the expense of a reduction in stream discharge. When the water table has reached the elevation of the stream bottom, this exchange from stream to lowland column may cease or reverse.  Figure 14 shows the profiles for the LIN_EQ simulation. The profiles in this simulation are roughly the same for all columns, with some differences in the early part of the time period due to lateral subsurface flow convergence in the lower columns. Because of the equal height binning used in LIN_EQ, the elevation of the lowland column is 18.4 m, and the bottom of the soil profile is therefore at an elevation of 10.4 m. Because the bottom of the soil profile is higher than the stream channel bankfull elevation (i.e., 0 m), lateral subsurface flow is unidirectional and there can be no exchange from the stream to the lowland column.
The bottom panel of Figure 14 shows soil moisture profiles at the same location in the LIN_VAR simulation. The profiles of the three upland columns are similar to those in the LIN_EQ simulation, while the lowland  column is significantly wetter. The higher soil moisture values in the LIN_VAR lowland column are the result of exchange with the stream channel. Because the height bin of the lowest column is specified to be 0 to 2 m, the elevation of the surface of the lowland column is 1.0 m and, therefore, the elevation of the bottom of the soil profile is −7.0 m. The stream channel bottom elevation at this location is −2.0 m, which results in a negative elevation gradient from the stream to the lowland column. When unsaturated conditions exist in the soil layers of the lowland column that are below the elevation of the stream channel bottom, water flows from the stream to the lowland column. Figure 15 shows profiles of soil moisture for a hillslope in a more humid location, a grid cell located at 310 E longitude and 20 S latitude. This grid cell is located in southeastern Brazil at an elevation of 417 m. Annual precipitation during the simulation period averaged approximately 1,310 mm. The average diurnal temperature range during JJA was 17 • C to 32 • C while in DJF, the range was 23 • C to 33 • C. The profiles for the LIN_EQ simulation (top panel) show the impact of lateral flow transport on soil water, as the thickness of the saturated zone progressively increases in the downhill direction. In the LIN_VAR simulation(bottom panel), similar behavior can be seen, although the lowland column shows higher soil water conditions because its lower elevation (relative to the lowland column of the LIN_EQ simulation) results in the stream channel being closer to the surface of the lowland column. The stream acts to buffer the seasonal drawdown of soil moisture apparent in the LIN_EQ simulation, leading to a less variable water table in the LIN_VAR simulation at this location. Figure 16 compares the mean annual water table depth from two simulations that are discretized using varying height bins: the LIN_VAR simulation and a simulation using the nonlinear hillslope profile shown in Figure 10d ("NONLIN_VAR"). The results shown in the left column of Figure 16 are for the LIN_VAR simulation and are thus the same as those shown in the right column of Figure 13. The right column of Figure 16 shows the effect of changes in topographic slope along the hillslope. The nonlinear hillslope profile is relatively flat near the ridge, resulting in smaller topographic gradients relative to the linear profile, while the midslope is more steeply sloped, leading to relatively large topographic gradients. These differences are reflected in the Column 4 panels, which show shallower water table depths in the NONLIN_VAR simulation in more humid locations. In contrast, the Column 3 panels indicate that Column 3 in the NONLIN_VAR simulation is actually drier than its upslope neighbor (Column 4) despite receiving additional water inputs from Column 4. Figure 17 shows profiles of soil moisture for a hillslope in the grid cell located at 5 E longitude and 14 N latitude. This grid cell is located in the Sahel region of southern Niger at an elevation of 268 m. Annual precipitation during the simulation period averaged approximately 505 mm. The average diurnal temperature range during JJA was 25 • C to 37 • C while in DJF, the range was 18 • C to 37 • C. The profiles for the LIN_VAR simulation (top panel) show a monotonic increase in saturated thickness due to increasing moisture convergence in the downhill direction. The NONLIN_VAR simulation, in contrast, shows that Column 3 actually exhibits moisture divergence, as outflow is greater than inflow due to the variation in topographic slope across the hillslope.

Evapotranspiration
Spatial patterns of evapotranspiration are often related to the spatial distribution of soil moisture where the latter is limiting (Koster et al., 2004). The simulations of the RME catchment showed that the differences in soil moisture due to the use of the Darcy flow, kinematic wave approximation, and no-lateral flow model configurations influenced the predicted lowland column evapotranspiration (Figure 9). To examine  the relationship between soil moisture conditions and ET globally, we performed variants of the LIN_VAR simulation using the kinematic wave approximation and no-lateral flow model configurations. Figure 18 shows that the water table depths in the lowland column are much shallower in the Darcy flow simulation than in the no-lateral flow simulation for most regions of the world. Exceptions include extremely wet regions (e.g., Amazonia), some deserts, and permafrost regions. In humid regions, the kinematic wave approximation simulation also shows shallower water tables as a result of moisture convergence of subsurface lateral flow but in drier regions, little or no differences are apparent relative to the no-lateral flow simulation. In some locations, higher ET results from increased soil moisture availability ( Figure 18). In both the Darcy and kinematic wave approximation simulations, this is most apparent in the tropics and subtropics. The Darcy simulation also shows higher ET values in arid and semiarid due to the buffering effect of bidirectional groundwater-surface water interaction. At higher latitudes, areas having shallower water tables do not consistently show higher ET values, indicating that these areas are likely not water limited.

Discussion and Summary
The RME catchment provides a rare opportunity for assessing the performance of the CLM hillslope configuration and for testing some underlying assumptions. The presence of spatial gradients in topography, vegetation type, and water inputs to the soil, combined with the availability of meteorological data to provide model boundary conditions and stream discharge and latent heat flux measurements with which to validate model output, make it well suited to address questions related to landscape spatial heterogeneity.
In particular, the data from the two EC measurement instruments reveal large differences in evapotranspiration occurring over a relatively small distance (i.e., with respect to the spatial scale of an ESM grid cell).
To explicitly simulate the observed ET difference using a gridded, global model would require a spatial resolution orders of magnitude higher than what is currently feasible. An alternative approach explored here is to apply the concept of hydrologic similarity to define a spatial discretization that minimizes the number of spatial elements while still differentiating between areas of distinct hydrologic function. Newman et al. (2014) also used this general approach to study the RME catchment, but an important difference here is that the catchment subunits (i.e., the hillslope columns) are hydrologically connected. The results shown in Section 4.1 indicate that redistribution of soil water within the catchment influences the availability of moisture for evapotranspiration and that without lateral moisture fluxes and the resulting convergence of moisture, the relatively high evapotranspiration rates in the lower portion of the catchment could not be supported.
The results of Section 4.1 show that the representative hillslope methodology can reproduce key features of the spatial distribution of hydrologic states within a single catchment but to apply the methodology globally requires the characterization of geomorphic parameters from high-resolution global topographic data. Suitable analyses are not presently available, but the development of global data sets such as the Multierror-Removed Improved-Terrain digital elevation model (Yamazaki et al., 2019) in conjunction with terrain analysis software such as the Soil Moisture and Runoff simulation Toolkit (Ajami et al., 2016) or the hierarchical multivariate clustering approach (Chaney et al., 2018) will provide the necessary information in the near future.
In the absence of globally extensive input data, we used synthetic geomorphic parameters derived from hypothetical hillslope profiles to examine the sensitivity of simulated soil water states to hillslope form and model structural decisions. The regions in which the simulations showed the greatest differences resulting from model structural choices were mainly in arid and semiarid climates. In these areas, simulations in which bidirectional groundwater-stream interactions were either explicitly disabled (kinematic wave and no-lateral flow) or implicitly disabled due to height binning that vertically separated the lowland column from the stream channel (equal height bins) resulted in minimal differences in soil moisture between upland and lowland columns. In contrast, the simulations in which bidirectional stream interactions were enabled (i.e., using Darcy's law and a 0-2-m height bin to define the lowland column) displayed significantly wetter conditions in the lowland column. This result is due to the fact that saturated conditions rarely develop in these regions and therefore moisture convergence due to saturated subsurface flow was not the dominant process by which lowland columns received nonlocal moisture inputs. Instead, surface runoff that was buffered by the stream and allowed to reinfiltrate into the lowland column was the process by which the model maintained higher riparian soil moisture.
In humid parts of the world, lateral subsurface flow caused differences in soil moisture not just in the lowland column but across the hillslope. However, by contrasting the linear and nonlinear hillslope profiles, it can be seen that both positive and negative moisture convergence are possible; and therefore, the location of the driest part of the hillslope does not depend solely on elevation but also on elevation gradient. While not evaluated here, variations along the hillslope profile of other geomorphic parameters such as area and width may also lead to nonmonotonically increasing variations in soil moisture across a hillslope.
Vegetation provides a key link between the hydrologic cycle and the carbon cycle, and the capacity of plants to assimilate carbon is tightly coupled to their ability to transpire. Because the lowland column has a significantly shallower water table depth than the upland columns in most regions in the LIN_VAR simulation, evapotranspiration is also higher. However, not all locations show an increase in ET, which likely indicates that these areas are not water limited. In such areas, which may be seasonally or perennially inundated, higher water tables may actually inhibit transpiration (Fan et al., 2017), but this process is not currently represented in CLM.
In summary, representing multiscale spatial heterogeneity in an ESM is a challenging task. One approach to capturing this heterogeneity is to simply increase the spatial resolution utilizing a finer spatial grid. However, this approach is computationally expensive and may not significantly close the gap between the model resolution and the spatial resolution at which key processes occur, in this case, the hillslope scale. For example, increasing the resolution of CLM from 100 to 10 km leads to a hundredfold increase in computational expense but still would not capture the hillslope-scale spatial gradients in vegetation distribution nor would it permit simulation of lateral flow across the landscape without reconfiguring the model to include grid cell to grid cell communication, which adds additional computational complexity and computational costs.
Here we show that, for a modest increase in computational cost, CLM can be configured to model these processes at the relevant spatial scales. The representative hillslope concept has additional benefits, one of which is the ability to relate simulated variables directly to observations at the spatial scales at which they are measured. Another benefit of this modeling approach is that covariation between hillslope elements and meteorological forcing (e.g., wind-induced snow redistribution), soil type and thickness, and vegetation distribution can be explicitly represented.