the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Development and evaluation of a variably saturated flow model in the global E3SM Land Model (ELM) version 1.0
William J. Riley
Glenn E. Hammond
David M. Lorenzetti
Improving globalscale model representations of nearsurface soil moisture and groundwater hydrology is important for accurately simulating terrestrial processes and predicting climate change effects on water resources. Most existing land surface models, including the default E3SM Land Model (ELMv0), which we modify here, routinely employ different formulations for water transport in the vadose and phreatic zones. Clark et al. (2015) identified a variably saturated Richards equation flow model as an important capability for improving simulation of coupled soil moisture and shallow groundwater dynamics. In this work, we developed the Variably Saturated Flow Model (VSFM) in ELMv1 to unify the treatment of soil hydrologic processes in the unsaturated and saturated zones. VSFM was tested on three benchmark problems and results were evaluated against observations and an existing benchmark model (PFLOTRAN). The ELMv1VSFM's subsurface drainage parameter, f_{d}, was calibrated to match an observationally constrained and spatially explicit global water table depth (WTD) product. Optimal spatially explicit f_{d} values were obtained for 79 % of global 1.9^{∘} × 2.5^{∘} grid cells, while the remaining 21 % of global grid cells had predicted WTD deeper than the observationally constrained estimate. Comparison with predictions using the default f_{d} value demonstrated that calibration significantly improved predictions, primarily by allowing much deeper WTDs. Model evaluation using the International Land Model Benchmarking package (ILAMB) showed that improvements in WTD predictions did not degrade model skill for any other metrics. We evaluated the computational performance of the VSFM model and found that the model is about 30 % more expensive than the default ELMv0 with an optimal processor layout. The modular software design of VSFM not only provides flexibility to configure the model for a range of problem setups but also allows for building the model independently of the ELM code, thus enabling straightforward testing of the model's physics against other models.
 Article
(6711 KB)  Fulltext XML

Supplement
(172 KB)  BibTeX
 EndNote
Groundwater, which accounts for 30 % of freshwater reserves globally, is a vital human water resource. It is estimated that groundwater provides 20–30 % of global freshwater withdrawals (Petra, 2009; Zektser and Evertt, 2004), and that irrigation accounts for ∼70 % of these withdrawals (Siebert et al., 2010). Climate change is expected to impact the quality and quantity of groundwater in the future (Alley, 2001). As temporal variability in precipitation and surface water increases in the future due to climate change, reliance on groundwater as a source of fresh water for domestic, agriculture, and industrial use is expected to increase (Taylor et al., 2013).
Local environmental conditions modulate the impact of rainfall changes on groundwater resources. For example, highintensity precipitation in humid areas may lead to a decrease in groundwater recharge (due to higher surface runoff), while arid regions are expected to see gains in groundwater storage (as infiltrating water quickly travels deep into the ground before it can be lost to the atmosphere; Kundzewicz and Doli, 2009). Although global climate models predict changes in precipitation over the next century (Marvel et al., 2017), few global models that participated in the recent Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2012) were able to represent global groundwater dynamics accurately (e.g., Swenson and Lawrence, 2014)
Modeling studies have also investigated impacts, at watershed to global scales, on future groundwater resources associated with landuse (LU) and landcover (LC) change (Dams et al., 2008) and ground water pumping (Ferguson and Maxwell, 2012; Leng et al., 2015). Dams et al. (2008) predicted that LU changes would result in a small mean decrease in subsurface recharge and large spatial and temporal variability in groundwater depth for the Kleine Nete basin in Belgium. Ferguson and Maxwell (2012) concluded that groundwaterfed irrigation impacts on water exchanges with the atmosphere and groundwater resources can be comparable to those from a 2.5 ^{∘}C increase in air temperature for the Little Washita basin in Oklahoma, USA. By performing global simulations of climate change scenarios using CLM4, Leng et al. (2015) concluded that the water source (i.e., surface or groundwater) used for irrigation depletes the corresponding water source while increasing the storage of the other water source. Recently, Leng et al. (2017) showed that the irrigation method (drip, sprinkler, or flood) has impacts on water balances and water use efficiency in global simulations.
Groundwater models are critical for developing understanding of groundwater systems and predicting impacts of climate (Green et al., 2011). Kollet and Maxwell (2008) identified critical zones, i.e., regions within the watershed with water table depths between 1 and 5 m, where the influence of groundwater dynamics was largest on surface energy budgets. Numerical studies have demonstrated impacts of groundwater dynamics on several key Earth system processes, including soil moisture (Chen and Hu, 2004; Liang et al., 2003; Salvucci and Entekhabi, 1995; Yeh and Eltahir, 2005), runoff generation (Levine and Salvucci, 1999; Maxwell and Miller, 2005; Salvucci and Entekhabi, 1995; Shen et al., 2013), surface energy budgets (Alkhaier et al., 2012; Niu et al., 2017; Rihani et al., 2010; Soylu et al., 2011), land–atmosphere interactions (Anyah et al., 2008; Jiang et al., 2009; Leung et al., 2011; Yuan et al., 2008), vegetation dynamics (Banks et al., 2011; Chen et al., 2010), and soil biogeochemistry (Lohse et al., 2009; Pacific et al., 2011).
Recognizing the importance of groundwater systems on terrestrial processes, groundwater models of varying complexity have been implemented in land surface models (LSMs) in recent years. Groundwater models in current LSMs can be classified into four categories based on their governing equations. Type1 models assume a quasisteady state equilibrium of the soil moisture profile above the water table (Hilberts et al., 2005; Koster et al., 2000; Walko et al., 2000). Type2 models use a θbased (where θ is the water volume content) Richards equation in the unsaturated zone coupled with a lumped unconfined aquifer model in the saturated zone. Examples of onedimensional Type2 models include Liang et al. (2003), Yeh and Eltahir (2005), Niu et al. (2007), and Zeng and Decker (2009). Examples of quasithreedimensional Type2 models are York et al. (2002), Fan et al. (2007), MiguezMacho et al. (2007), and Shen et al. (2013). Type3 models include a threedimensional representation of subsurface flow based on the variably saturated Richards equation (Maxwell and Miller, 2005; Tian et al., 2012). Type3 models employ a unified treatment of hydrologic processes in the vadose and phreatic zones but lump changes associated with water density and unconfined aquifer porosity into a specific storage term. The fourth class (Type4) of subsurface flow and reactive transport models (e.g., PFLOTRAN, Hammond and Lichtner, 2010; TOUGH2, Pruess et al., 1999; and STOMP, White and STOMP, 2000) combine a water equation of state (EoS) and soil compressibility with the variably saturated Richards equation. Type4 models have not been routinely coupled with LSMs to address climate change relevant research questions. Clark et al. (2015) summarized that most LSMs use different physics formulations for representing hydrologic processes in saturated and unsaturated zones. Additionally, Clark et al. (2015) identified the incorporation of variably saturated hydrologic flow models (i.e., Type3 and Type4 models) in LSMs as a key opportunity for future model development that is expected to improve simulation of coupled soil moisture and shallow groundwater dynamics.
The Energy, Exascale, Earth System Model (E3SM) is a new Earth system modeling project sponsored by the US Department of Energy (E3SM Project, 2018). The E3SM model started from the Community Earth System Model (CESM) version 1_3_beta10 (Oleson et al., 2013). Specifically, the initial version (v0) of the E3SM Land Model (ELM) was based off the Community Land Model's (CLM's) tag 4_5_71. ELMv0 uses a Type2 subsurface hydrology model based on Zeng and Decker (2009). In this work, we developed in ELMv1 the Type4 Variably Saturated Flow Model (VSFM) to provide a unified treatment of soil hydrologic processes within the unsaturated and saturated zones. The VSFM formulation is based on the isothermal, single phase, variably saturated (RICHARDS) flow model within PFLOTRAN (Hammond and Lichtner, 2010). While PFLOTRAN is a massively parallel, threedimensional subsurface model, the VSFM is a serial, onedimensional model that is appropriate for climatescale applications.
This paper is organized into several sections: (1) a brief review of the ELMv0 subsurface hydrology model; (2) an overview of the VSFM formulation integrated in ELMv1; (3) an application of the new model formulation to three benchmark problems; (4) development of a subsurface drainage parameterization necessary to predict global water table depths (WTDs) comparable to recently released observationally constrained estimates; (5) comparison of ELMv1 global simulations with the default subsurface hydrology model and VSFM against multiple observations using the International Land Model Benchmarking package (ILAMB; Hoffman et al., 2017); and (6) a summary of major findings.
2.1 Current model formulation
Water flow in the unsaturated zone is often described by the θbased Richards equation:
where θ (m^{3} of water m^{−3} of soil) is the volumetric soil water content, t (s) is time, q (m s^{−1}) is the Darcy water flux, and Q (m^{3} of water m^{−3} of soil s^{−1}) is a soil moisture sink term. The Darcy flux, q, is given by
where K (m s^{−1}) is the hydraulic conductivity, z (m) is height above some datum in the soil column and ψ (m) is the soil matric potential. The hydraulic conductivity and soil matric potential are modeled as nonlinear functions of volumetric soil moisture following Clapp and Hornberger (1978):
where K_{sat} (m s^{−1}) is saturated hydraulic conductivity, ψ_{sat} (m) is saturated soil matric potential, B is a linear function of percentage clay and organic content (Oleson et al., 2013), and Θ_{ice} is the ice impedance factor (Swenson et al., 2012). ELMv0 uses the modified form of the Richards equation of Zeng and Decker (2009) that computes Darcy flux as
where C is a constant hydraulic potential above the water table, z_{∇}, given as
where ψ_{E} (m) is the equilibrium soil matric potential, z (m) is height above a reference datum, θ_{E} (m^{3} m^{−3}) is volumetric soil water content at equilibrium soil matric potential, and z_{∇} (m) is the height of the water table above a reference datum. ELMv0 uses a cellcentered finite volume spatial discretization and backward Euler implicit time integration. By default, ELMv0's vertical discretization of a soil column yields 15 soil layers of exponentially varying soil thicknesses that reach a depth of 42.1 m Only the first 10 soil layers (or top 3.8 m of each soil column) are hydrologically active, while thermal processes are resolved for all 15 soil layers. The nonlinear Darcy flux is linearized using Taylor series expansion and the resulting tridiagonal system of equations is solved by LU factorization.
Flow in the saturated zone is modeled as an unconfined aquifer below the soil column based on the work of Niu et al. (2007). Exchange of water between the soil column and unconfined aquifer depends on the location of the water table. When the water table is below the last hydrologically active soil layer in the column, a recharge flux from the last soil layer replenishes the unconfined aquifer. A zeroflux boundary condition is applied to the last hydrologically active soil layer when the water table is within the soil column. The unconfined aquifer is drained by a flux computed based on the SIMTOP scheme of Niu et al. (2007) with modifications to account for frozen soils (Oleson et al., 2013).
2.2 New VSFM model formulation
In the VSFM formulation integrated in ELMv1, we use the mass conservative form of the variably saturated subsurface flow equation (Farthing et al., 2003; Hammond and Lichtner, 2010; Kees and Miller, 2002):
where ϕ (m^{3} m^{−3}) is the soil porosity, s_{w} () is saturation, ρ (kg m^{−3}) is water density, q (m s^{−1}) is the Darcy velocity, and Q (kg m^{−3} s^{−1}) is a water sink. We restrict our model formulation to a onedimensional system and the flow velocity is defined by Darcy's law:
where k (m^{2}) is intrinsic permeability, k_{r} () is relative permeability, μ (Pa s) is viscosity of water, P (Pa) is pressure, g (m s^{−2}) is the acceleration due to gravity, and z (m) is elevation above some datum in the soil column.
In order to close the system, a constitutive relationship is used to express saturation and relative permeability as a function of soil matric pressure. Analytic water retention curves (WRCs) are used to model effective saturation (s_{e})
where s_{w} is saturation and s_{r} is residual saturation. We have implemented Brooks and Corey (1964; Eq. 10) and van Genuchten (1980; Eq. 11) WRCs:
where P_{c} (Pa) is the capillary pressure, ${P}_{\mathrm{c}}^{\mathrm{0}}$ (Pa) is the air entry pressure, α (Pa^{−1}) is inverse of the air entry pressure, λ () is the exponent in the Brooks and Corey parameterization, and n () and m () are the van Genuchten parameters. The capillary pressure is computed as ${P}_{\mathrm{c}}=P{P}_{\mathrm{ref}}$ where P_{ref} is ${P}_{\mathrm{c}}^{\mathrm{0}}$ for Brooks and Corey WRC and typically the atmospheric pressure (=101 325 Pa) is used for the van Genuchten WRC. In addition, a smooth approximation of Eqs. (10) and (11) was developed to facilitate convergence of the nonlinear solver (Appendix A). Relative soil permeability was modeled using the Mualem (1976) formulation:
where $m=\mathrm{1}\mathrm{1}/n$. Lastly, we used an EoS for water density, ρ, that is a nonlinear function of liquid pressure, P, and liquid temperature, T, given by Tanaka et al. (2001):
where
The sink of water due to transpiration from a given plant functional type (PFT) is vertically distributed over the soil column based on area and root fractions of the PFT. The top soil layer has an additional flux associated with balance of infiltration and soil evaporation. The subsurface drainage flux is applied proportionally to all soil layers below the water table. Details on the computation of water sinks are given in Oleson et al. (2013). Unlike the default subsurface hydrology model, the VSFM is applied over the full soil depth (in the default model, 15 soil layers). The VSFM model replaces both the θbased Richards equation and the unconfined aquifer of the default model and uses a zeroflux lower boundary condition. In the VSFM model, water table depth is diagnosed based on the vertical soil liquid pressure profile. Like the default model, drainage flux is computed based on the modified SIMTOP approach and is vertically distributed over the soil layers below the water table.
2.2.1 Discrete equations
We use a cellcentered finite volume discretization to decompose the spatial domain, Ω, into N nonoverlapping control volumes, Ω_{n}, such that $\mathrm{\Omega}={\cup}_{n=\mathrm{1}}^{N}{\mathrm{\Omega}}_{i}$ and Γ_{n} represents the boundary of the nth control volume. Applying a finite volume integral to Eq. (7) and the divergence theorem yields
The discretized form of the lefthandside term and first term on the righthand side of Eq. (14) are approximated as
where ${\mathbf{A}}_{nn{}^{\prime}}$ (m^{2}) is the common face area between the nth and n^{′}th control volumes. After substituting Eqs. (15) and (16) in Eq. (14), the resulting ordinary differential equation for the variably saturated flow model is
We perform temporal integration of Eq. (17) using the backwardEuler scheme:
Rearranging terms of Eq. (18) results in a nonlinear equation for the unknown pressure at time step t+1 as
In this work, we find the solution to the nonlinear system of equations given by Eq. (19) using Newton's method via the Scalable Nonlinear Equations Solver (SNES) within the Portable, Extensible Toolkit for Scientific Computing (PETSc) library (Balay et al., 2016). PETSc provides a suite of data structures and routines for the scalable solution of partial differential equations. VSFM uses the composable data management (DMComposite) provided by PETSc (Brown et al., 2012), which enables the potential future application of the model to solve tightly coupled, multicomponent, multiphysics processes as discussed in Sect. 3.4. A smooth approximation of the Brooks and Corey (1964; SBC) water retention curve was developed to facilitate faster convergence of the nonlinear solver (Appendix A). ELMv0 code for subsurface hydrologic processes only supports two vertical mesh configurations and a single set of boundary and source–sink conditions. The monolithic ELMv0 code does not allow for building of code for individual process representations independent of ELMv0 code, thus precluding easy testing of the model against analytical solutions or simulation results from other models. The modular software design of VSFM overcomes ELMv0's software limitation by allowing VSFM code to be built independently of the ELM code. This flexibility of VSFM's build system allows for testing of the VSFM physics in isolation without any influence from the rest of ELM's physics formulations. Additionally, VSFM can be easily configured for a wide range of benchmark problems with different spatial grid resolutions, material properties, boundary conditions, and source–sink forcings.
2.3 VSFM singlecolumn evaluation
We tested the VSFM with three idealized 1dimensional test problems. First, the widely studied problem for the 1D Richards equation of infiltration in dry soil by Celia et al. (1990) was used. The problem setup consists of a 1.0 m long soil column with a uniform initial pressure of −10.0 m (=3535.5 Pa). The time invariant boundary conditions applied at the top and bottom of the soil column are −0.75 m (=93989.1 Pa) and −10.0 m (=3535.5 Pa), respectively. The soil properties for this test are given in Table 1. A vertical discretization of 0.01 m is used in this simulation.
Second, we simulated transient onedimensional vertical infiltration in a twolayered soil system as described in Srivastava and Yeh (1991). The domain consisted of a 2 m tall soil column divided equally into two soil types. Except for soil intrinsic permeability, all other soil properties of the two soil types are the same. The bottom soil is 10 times less permeable than the top (Table 1). Unlike Srivastava and Yeh (1991), who used exponential functions of soil liquid pressure to compute hydraulic conductivity and soil saturation, we used Mualem (1976) and van Genuchten (1980) constitutive relationships. Since our choice of constitutive relationships for this setup resulted in absence of an analytical solution, we compared VSFM simulations against PFLOTRAN results. The domain was discretized in 200 control volumes of equal soil thickness. Two scenarios, wetting and drying, were modeled to test the robustness of the VSFM solver robustness. Initial conditions for each scenario included a time invariant boundary condition of 0 m ($=\mathrm{1.01325}\times {\mathrm{10}}^{\mathrm{5}}$ Pa) for the lowest control volume and a constant flux of 0.9 and 0.1 cm h^{−1} at the soil surface for wetting and drying scenarios, respectively.
Third, we compare VSFM and PFLOTRAN predictions for soil under variably saturated conditions. The 1dimensional 1 m deep soil column was discretized in 100 equal thickness control volumes. A hydrostatic initial condition was applied such that the water table is 0.5 m below the soil surface. A time invariant flux of $\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{5}}$ m s^{−1} is applied at the surface, while the lowest control volume has a boundary condition corresponding to the initial pressure value at the lowest soil layer. The soil properties used in this test are the same as those used in the first evaluation.
2.4 Global simulations and groundwater depth analysis
We performed global simulations with ELMv1VSFM at a spatial resolution of 1.9^{∘} (latitude) × 2.5^{∘} (longitude) with a 30 (min) time step for 200 years, including a 180year spinup and the last 20 years for analysis. The simulations were driven by CRUNCEP meteorological forcing from 1991 to 2010 (Piao et al., 2012) and configured to use prescribed satellite phenology.
For evaluation and calibration, we used the Fan et al. (2013) global ∼1 km horizontal resolution WTD dataset (hereafter F2013 dataset), which is based on a combination of observations and hydrologic modeling. We aggregated the dataset to the ELMv1VSFM spatial resolution. ELMVSFM's default vertical soil discretization uses 15 soil layers to a depth of ∼42 m, with an exponentially varying soil thickness. However, ∼13 % of F2013 land grid cells have a water table deeper than 42 m. We therefore modified ELMv1VSFM to extend the soil column to a depth of 150 m with 59 soil layers; the first nine soil layer thicknesses were the same as described in Oleson et al. (2013) and the remaining layers (10–59) were set to a thickness of 3 m.
2.5 Estimation of the subsurface drainage parameterization
In the VSFM formulation, the dominant control on longterm ground water depth is the subsurface drainage flux, q_{d} (kg m^{−2} s^{−1}), which is calculated based on water table depth, z_{∇} (m), Niu et al. (2005):
where q_{d,max} (kg m^{−2} s^{−1}) is the maximum drainage flux that depends on grid cell slope and f_{d} (m^{−1}) is an empirically derived parameter. The subsurface drainage flux formulation of Niu et al. (2005) is similar to the TOPMODEL formulation (Beven and Kirkby, 1979) and assumes the water table is parallel to the soil surface. While Sivapalan et al. (1987) derived q_{d,max} as a function of lateral hydraulic anisotropy, hydraulic conductivity, topographic index, and decay factor controlling vertical saturated hydraulic conductivity, Niu et al. (2005) defined q_{d,max} as a single calibration parameter. ELMv0 uses f_{d}=2.5 m^{−1} as a global constant and estimates maximum drainage flux when WTD is at the surface as ${q}_{\mathrm{d},\mathrm{max}}=\mathrm{10}\mathrm{sin}\left(\mathit{\beta}\right)$ kg m^{−2} s^{−1}, where β (radians) is the mean grid cell topographic slope. Of the two parameters, f_{d} and q_{d,max}, available for model calibration, we choose to calibrate f_{d} because the uncertainty analysis by Hou et al. (2012) identified it as the most significant hydrologic parameter in CLM4. To improve on the f_{d} parameter values, we performed an ensemble of global simulations with f_{d} values of 0.1, 0.2, 0.5, 1.0, 2.5, 5.0, 10.0, and 20 m^{−1}. Each ensemble simulation was run for 200 years to ensure an equilibrium solution, and the last 20 years were used for analysis. A nonlinear functional relationship between f_{d} and “WTD” was developed for each grid cell and then the F2013 dataset was used to estimate an optimal f_{d} for each grid cell.
2.6 Global ELMVSFM evaluation
With the optimal f_{d} values, we ran a ELMVSFM simulation using the protocol described above. We then used the ILAMB package to evaluate the ELMv1VSFM predictions of surface energy budget, terrestrial water storage anomalies (TWSAs), and river discharge (Collier et al., 2018; Hoffman et al., 2017). ILAMB evaluates model prediction bias, RMSE, and seasonal and diurnal phasing against multiple observations of energy, water, and carbon cycles at insitu, regional, and global scales. Since ELMVSFM simulations in this study did not include an active carbon cycle, we used the following ILAMB benchmarks for water and energy cycles: (i) latent and surface energy fluxes using sitelevel measurements from FLUXNET (Lasslop et al., 2010) and globally from FLUXNETMTE (Jung et al., 2009); (ii) TWSA from the Gravity Recovery And Climate Experiment (GRACE) observations (Kim et al., 2009); and (iii) stream flow for the 50 largest global river basins (Dai and Trenberth, 2002). We applied ILAMB benchmarks for ELMv1VSFM simulations with default and calibrated f_{d} to ensure improvements in WTD predictions did not degrade model skill for other processes.
3.1 VSFM singlecolumn evaluation
For the 1D Richards equation infiltration in dry soil comparison, we evaluated the solutions at 24 h against those published by Celia et al. (1990; Fig. 1). The VSFM solver accurately represented the sharp wetting front over time, where soil hydraulic properties change dramatically due to nonlinearity in the soil water retention curve.
For the model evaluation of infiltration and drying in layered soil, the results of the VSFM and PFLOTRAN are essentially identical. In both models and scenarios, the higher permeability top soil responds rapidly to changes in the top boundary condition and the wetting and drying fronts progressively travel through the less permeable soil layer until soil liquid pressure in the entire column reaches a new steady state by about 100 h (Fig. 2).
We also evaluated the VSFM predicted water table dynamics against PFLOTRAN predictions from an initial condition of saturated soil below 0.5 m depth. The simulated water table rises to 0.3 m depth by 1 day and reaches the surface by 2 days, and the VSFM and PFLOTRAN predictions are essentially identical (Fig. 3). Soil properties, spatial discretization, and time step used for the three singlecolumn problems are summarized in Table 1 These three evaluation simulations demonstrate the VSFM accurately represents soil moisture dynamics under conditions relevant to ESMscale prediction.
3.2 Subsurface drainage parameterization estimation
The simulated nonlinear WTD–f_{d} relationship is a result of the subsurface drainage parameterization flux given by Eq. (20; Fig. 4a and b). For $\mathrm{0.1}\le {f}_{\mathrm{d}}\le \mathrm{1}$, the slope of the WTD–f_{d} relationship for all grid cells is log–log linear with a slope of $\mathrm{1.0}\pm \mathrm{0.1}$. The log–log linear relationship breaks down for f_{d}>1, where the drainage flux becomes much smaller than infiltration and evapotranspiration (Fig. 4c and d). Thus, at larger f_{d}, the steady state z_{∇} becomes independent of f_{d} and is determined by the balance of infiltration and evapotranspiration.
For 79 % of the global grid cells, the ensemble range of simulated WTD spanned the F2013 dataset. The optimal value of f_{d} for each of these grid cells was obtained by linear interpolation in the log–log space (e.g., Fig. 4a). For the remaining 21 % of grid cells where the shallowest simulated WTD across the range of f_{d} was deeper than that in the F2013 dataset, the optimal f_{d} value was chosen as the one that resulted in the lowest absolute WTD error (e.g., Fig. 4b). At large f_{d} values, the drainage flux has negligible effects on WTD, yet simulated WTD is not sufficiently shallow to match the F2013 observations, which indicates that either evapotranspiration is too large or infiltration is too small. There was no difference in the mean percentage of sand and clay content between grids cells with and without an optimal f_{d} value. The optimal f_{d} has a global average of 1.60 m${}^{\mathrm{1}}\pm \mathrm{2.68}$ m^{−1} and 72 % of global grid cells have an optimal f_{d} value lower than the global average (Fig. 5).
3.3 Global simulation evaluation
The ELMv1VSFM predictions are much closer to the F2013 dataset (Fig. 6a) using optimal globally distributed f_{d} values (Fig. 6c) compared to the default f_{d} value (Fig. 6b). The significant reduction in WTD bias (model – observation) is mostly due to improvement in the model's ability to accurately predict deep WTD using optimal f_{d} values. In the simulation using optimal globally distributed f_{d} values, all grid cells with WTD bias >3.7 m were those for which an optimal f_{d} was not found. The mean global bias, RMSE, and R^{2} values improved in the new ELMv1VSFM compared to the default model (Table 2). The 79 % of global grid cells for which an optimal f_{d} value was estimated had significantly better water table prediction with a bias, RMSE, and R^{2} of −0.04, 0.67, and 0.99 m, respectively, as compared to the remaining 21 % of global grid cells that had a bias, RMSE, and R^{2} of −9.82, 18.08, and 0.31 m, respectively. The simulated annual WTD range, which we define to be the difference between maximum and minimum WTD in 1 year, has a spatial mean and standard deviation of 0.32 and 0.58 m, respectively, using optimal f_{d} values (Fig. 7a). The annual WTD range decreased by 0.24 m for the 79 % of the grid cells for which an optimal f_{d} value was estimated (Fig. 7b).
Globally averaged WTD in ELMv1VSFM simulations with default f_{d} and optimal f_{d} values were 10.5 and 20.1 m, respectively. Accurate prediction of deep WTD in the simulation with optimal f_{d} caused very small differences in nearsurface soil moisture (Fig. 8). The 79 % of grid cells with an optimal f_{d} value had deeper globally averaged WTDs than when using the default f_{d} value (24.3 vs. 8.6 m). For these 79 % of grid cells, the WTD was originally deep enough to not impact nearsurface conditions (Kollet and Maxwell, 2008); therefore, further lowering of WTD led to negligible changes in nearsurface hydrological conditions.
The ILAMB package (Hoffman et al., 2017) provides a comprehensive evaluation of predictions of carbon cycle states and fluxes, hydrology, surface energy budgets, and functional relationships by comparison to a wide range of observations. We used ILAMB to evaluate the hydrologic and surface energy budget predictions from the new ELMv1VSFM model (Table 3). Optimal f_{d} values had inconsequential impacts on simulated surface energy fluxes at sitelevel and global scales. Optimal f_{d} values led to improvement in prediction of deep WTD (with a mean value of 24.3 m) for grid cells that had an average WTD of 8.7 m in the simulation using default f_{d} values. Thus, negligible differences in surface energy fluxes between the two simulations are consistent with the findings of Kollet and Maxwell (2008), who identified decoupling of groundwater dynamics and surface processes at a WTD of ∼10 m. There were slight changes in the bias and RMSE for predicted TWSA, but the ILAMB score remained unchanged. The TWSA amplitude is lower for the simulation with optimal f_{d} values, consistent with the associated decrease in annual WTD range. ELM's skill in simulating runoff for the 50 largest global watersheds remained unchanged. Two additional 10yearlong simulations were performed to investigate the sensitivity of VSFM simulated WTD to spatial and temporal discretization. Results show that simulated WTD is insensitive to temporal discretization, and has small sensitivity to vertical spatial resolution. See the Supplement for details regarding the setup and analysis of results from the two additional simulations.
Finally, we evaluated the computational costs of implementing VSFM in ELM and compared them to the default model. We performed 5yearlong simulations for default and VSFM using 96, 192, 384, 768, and 1536 cores on the Edison supercomputer at the National Energy Research Scientific Computing Center. Using an optimal processor layout, we found that ELMv1VSFM is ∼30 % more expensive than the default ELMv1 model (the Supplement, Fig. S1). We note that the relative computational cost of the land model in a fully coupled global model simulation is generally very low. Dennis et al. (2012) reported computational cost of the land model to be less than 1 % in ultrahighresolution CESM simulations. We therefore believe that the additional benefits associated with the VSFM formulation are well justified by this modest increase in computational cost. In particular, VSFM allows for a greater variety of mesh configurations and boundary conditions, and can accurately simulate WTD for the ∼13 % of global grid cells that have a water table deeper than 42 (m) (Fan et al., 2013).
3.4 Caveats and future work
The significant improvement in WTD prediction using optimal f_{d} values demonstrates VSFM's capabilities to model hydrologic processes using a unified physics formulation for unsaturated–saturated zones. However, several caveats remain due to uncertainties in model structure, model parameterizations, and climate forcing data.
In this study, we assumed a spatially homogeneous depth to bedrock (DTB) of 150 m. Recently, Brunke et al. (2016) incorporated a global ∼1 km dataset of soil thickness and sedimentary deposits (Pelletier et al., 2016) in CLM4.5 to study the impacts of soil thickness spatial heterogeneity on simulated hydrological and thermal processes. While inclusion of heterogeneous DTB in CLM4.5 added more realism to the simulation setup, no significant changes in simulated hydrologic and energy fluxes were reported by Brunke et al. (2016). Presently, work is ongoing in the E3SM project to include variable DTB within ELM and future simulations will examine the impact of those changes on VSFM's prediction of WTD. Our use of the “satellite phenology” mode, which prescribes transient leaf area index profiles for each PFT in the grid cell, ignored the likely influence of water cycle dynamics and nutrient constraints on the C cycle (Ghimire et al., 2016; Zhu et al., 2016). Estimation of soil hydraulic properties based on soil texture data is critical for accurate LSM predictions (Gutmann and Small, 2005) and this study does not account for uncertainty in soil hydraulic properties.
Lateral water redistribution impacts soil moisture dynamics (Bernhardt et al., 2012), biogeochemical processes in the root zone (Grant et al., 2015), distribution of vegetation structure (Hwang et al., 2012), and land–atmosphere interactions (Chen and Kumar, 2001; Rihani et al., 2010). The ELMv1VSFM developed in this study does not include lateral water redistribution between soil columns and only simulates vertical water transport. Lateral subsurface processes can be included in LSMs via a range of numerical discretization approaches of varying complexity, e.g., adding lateral water as source/sink terms in the 1D model, implementing an operator split approach to solve vertical and lateral processes in a noniterative approach (Ji et al., 2017), or solving a fully coupled 3D model (Bisht et al., 2017, 2018; Kollet and Maxwell, 2008). Additionally, lateral transport of water can be implemented in LSMs at a subgrid level (Milly et al., 2014) or grid cell level (MiguezMacho et al., 2007). The current implementation of VSFM is such that each processor solves the variably saturated Richards equation for all independent soil columns as one single problem. Thus, extension of VSFM to solve the tightly coupled 3D Richards equation on each processor locally while accounting for lateral transport of water within grid cells and among grid cells is straightforward. The current VSFM implementation can also be easily extended to account for subsurface transport of water among grid cells that are distributed across multiple processors by modeling lateral flow as source/sink terms in the 1D model. Tradeoffs between approaches to represent lateral processes and computational costs need to be carefully studied before developing quasi or fully threedimensional LSMs (Clark et al., 2015).
Transport of water across multiple components of the soilplantatmosphere continuum (SPAC) has been identified as a critical process in understanding the impact of climate warming on the global carbon cycle (McDowell and Allen, 2015). Several SPAC models have been developed by the ecohydrology community and applied to study site level processes (Amenu and Kumar, 2008; Bohrer et al., 2005; Manoli et al., 2014; Sperry et al., 1998), yet implementation of SPAC models in global LSMs is limited (Clark et al., 2015). Similarly, current generation LSMs routinely ignore advective heat transport within the subsurface, which has been shown to be important in highlatitude environments by multiple field and modeling studies (Bense et al., 2012; Frampton et al., 2011; Grant et al., 2017; Kane et al., 2001). The use of PETSc's DMComposite in VSFM provides flexibility for solving a tightly coupled multicomponent problem (e.g., transport of water through the soilplantatmosphere continuum) and multiphysics problem (e.g., fully coupled conservation of mass and energy equations in the subsurface). DMComposite allows for an easy assembly of a tightly coupled multiphysics problem from individual physics formulations (Brown et al., 2012).
Starting from the climatescale land model ELMv0, we incorporated a unified physics formulation to represent soil moisture and groundwater dynamics that are solved using PETSc. Application of VSFM to three benchmarks problems demonstrated its robustness to simulate subsurface hydrologic processes in coupled unsaturated and saturated zones. Ensemble global simulations at 1.9^{∘} × 2.5^{∘} were performed for 200 years to obtain spatially heterogeneous estimates of the subsurface drainage parameter, f_{d}, that minimized mismatches between predicted and observed WTDs. In order to simulate the deepest water table reported in the Fan et al. (2013) dataset, we used 59 vertical soil layers that reached a depth of 150 m.
An optimal f_{d} was obtained for 79 % of the grids cells in the domain. For the remaining 21 % of grid cells, simulated WTD always remained deeper than observed. Calibration of f_{d} significantly improved global WTD prediction by reducing bias and RMSE and increasing R^{2}. Grids without an optimal f_{d} were the largest contributor to error in WTD predication. ILAMB benchmarks on simulations with default and optimal f_{d} showed negligible changes to surface energy fluxes, TWSA, and runoff. ILAMB metrics ensured that model skill was not adversely impacted for all other processes when optimal f_{d} values were used to improve WTD prediction.
The Brooks and Corey (1964) water retention curve of Eq. (10) has a discontinuous derivative at ${P}_{\mathrm{c}}={P}_{\mathrm{c}}^{\mathrm{0}}$. Figure A1 illustrates an example. To improve convergence of the nonlinear solver at small capillary pressures, the smoothed Brooks–Corey function introduces a cubic polynomial, B(P_{c}), in the neighborhood of ${P}_{\mathrm{c}}^{\mathrm{0}}$.
where the breakpoints P_{u} and P_{s} satisfy ${P}_{\mathrm{u}}<{P}_{\mathrm{c}}^{\mathrm{0}}<{P}_{\mathrm{s}}\le \mathrm{0}$. The smoothing polynomial
introduces four more parameters, whose values follow from continuity. In particular, matching the saturated region requires $B\left({P}_{\mathrm{s}}\right)={b}_{\mathrm{0}}=\mathrm{1}$, and a continuous derivative at P_{c}=P_{s} requires $B{}^{\prime}\left({P}_{\mathrm{s}}\right)={b}_{\mathrm{1}}=\mathrm{0}$. Similarly, matching the value and derivative at P_{c}=P_{u} requires
where $\mathrm{\Delta}={P}_{\mathrm{u}}{P}_{\mathrm{s}}$. Note that ${P}_{\mathrm{u}}\le \mathrm{\Delta}<\mathrm{0}$.
In practice, setting P_{u} too close to ${P}_{\mathrm{c}}^{\mathrm{0}}$ can produce an unwanted local maximum in the cubic smoothing regime, resulting in s_{e}>1. Avoiding this condition requires that B(P_{c}) increase monotonically from P_{c}=P_{u}, where $B{}^{\prime}\left({P}_{\mathrm{c}}\right)>\mathrm{0}$, to P_{c}=P_{s}, where $B{}^{\prime}\left({P}_{\mathrm{c}}\right)=\mathrm{0}$. Thus a satisfactory pair of breakpoints ensures
throughout ${P}_{\mathrm{u}}\le {P}_{\mathrm{c}}<{P}_{\mathrm{s}}$.
Let ${P}_{\mathrm{c}}^{\ast}$ denote a local extremum of B, so that $B{}^{\prime}\left({P}_{\mathrm{c}}^{\ast}\right)=\mathrm{0}$. If ${P}_{\mathrm{c}}^{\ast}\ne {P}_{\mathrm{s}}$, it follows ${P}_{\mathrm{c}}^{\ast}{P}_{\mathrm{s}}=\mathrm{2}{b}_{\mathrm{2}}/\left(\mathrm{3}{b}_{\mathrm{3}}\right)$. Rewriting Eq. (22), $B{}^{\prime}\left({P}_{\mathrm{c}}\right)=\left({P}_{\mathrm{c}}{P}_{\mathrm{s}}\right)\mathrm{3}{b}_{\mathrm{3}}\left({P}_{\mathrm{c}}{P}_{\mathrm{c}}^{\ast}\right)$ shows that $B{}^{\prime}\left({P}_{\mathrm{u}}\right)>\mathrm{0}$ requires either: (1) b_{3}<0 and ${P}_{\mathrm{c}}^{\ast}<{P}_{\mathrm{u}}$; or (2) b_{3}>0 and ${P}_{\mathrm{c}}^{\ast}>{P}_{\mathrm{u}}$. The first possibility places ${P}_{\mathrm{c}}^{\ast}$ outside the cubic smoothing regime, and so does not constrain the choice of P_{u} or P_{s}. The second possibility allows for an unwanted local extremum at ${P}_{\mathrm{u}}<{P}_{\mathrm{c}}^{\ast}<{P}_{\mathrm{s}}$. In this case, b_{3}>0 implies b_{2}<0 (since ${P}_{\mathrm{c}}^{\ast}<{P}_{\mathrm{s}}\le \mathrm{0})$. Then since $B{}^{\prime \prime}\left({P}_{\mathrm{c}}^{\ast}\right)=\mathrm{2}{b}_{\mathrm{2}}$, the local extremum is a maximum, resulting in ${s}_{\mathrm{e}}\left({P}_{\mathrm{c}}^{\ast}\right)>\mathrm{1}$.
Given a breakpoint P_{s}, one strategy for choosing P_{u} is to guess a value, then check whether the resulting b_{2} and b_{3} produces ${P}_{\mathrm{u}}<{P}_{\mathrm{c}}^{\ast}<{P}_{\mathrm{s}}$. If so, P_{u} should be made more negative. An alternative strategy is to choose P_{u} in order the guarantee acceptable values for b_{2} and b_{3}. One convenient choice forces b_{2}=0. Another picks P_{u} in order to force b_{3}=0. Both of these reductions (1) ensure B(P_{c}) has a positive slope throughout the smoothing interval, (2) slightly reduce the computation cost of finding s_{e}(P_{c}) for P_{c} on the smoothing interval, and (3) significantly reduce the computational cost of inverting the model, in order to find P_{c} as a function of s_{e}.
As shown in Fig. A1, the two reductions differ mainly in that setting b_{2}=0 seems to produce narrower smoothing regions (probably due to the fact that this choice gives zero curvature at P_{c}=P_{s}, while b_{3}=0 yields a negative second derivative there). However, we have not verified this observation analytically.
Both reductions require solving a nonlinear expression, either Eqs. (A3) or (A4), for P_{u}. While details are beyond the scope of this paper, we note that we have used a bracketed Newton–Raphson method. The search switches to bisection when Newton–Raphson would jump outside the bounds established by previous iterations, and by the requirement ${P}_{\mathrm{u}}<{P}_{\mathrm{c}}^{\mathrm{0}}$. In any event, since the result of this calculation may be cached for use throughout the simulation, it need not be particularly efficient.
The residual equation for the VSFM formulation at t+1 time level for nth control volume is given by
where ϕ (mm^{3} mm^{3}) is the soil porosity, s_{w} () is saturation, ρ (kg m^{−3}) is water density, ${}_{nn{}^{\prime}}$ (m s^{−1}) is the Darcy flow velocity between nth and n^{′}th control volumes, ${A}_{nn{}^{\prime}}$ (m s^{−1}) is the interface area between nth and n^{′}th control volumes, and Q (kg m^{−3} s^{−1}) is a sink of water. The Darcy velocity is computed as
where κ (m^{−2}) is intrinsic permeability, κ_{r} () is relative permeability, μ (Pa s) is viscosity of water, P (Pa) is pressure, g (m s^{−2}) is the acceleration due to gravity, d_{n} (m) and ${d}_{n{}^{\prime}}$ (m) are the distances between centroid of the nth and n^{′}th control volume to the common interface between the two control volumes, ${\mathbf{d}}_{nn{}^{\prime}}$ is a distance vector joining centroid of nth and n^{′}th control volume, and ${\mathbf{n}}_{nn{}^{\prime}}$ is a unit normal vector joining centroid of the nth and n^{′}th control volume.
The density at the interface of control volume, ${\mathit{\rho}}_{nn{}^{\prime}}$, is computed as the inverse distance weighted average by
where ω_{n} and ${\mathit{\omega}}_{n{}^{\prime}}$ are given by
The first term on the righthand side of Eq. (27) is computed as the product of distanceweighted harmonic average of intrinsic permeability, ${k}_{nn{}^{\prime}}$, and upwinding of ${k}_{r}/\mathit{\mu}(=\mathit{\lambda})$ as
where
By substituting Eqs. (28), (29), and (30) in Eq. (27), we obtain
The discretized equations of VSFM leads to a system of nonlinear equations given by ${\mathbf{R}}^{t+\mathrm{1}}\left({\mathbf{P}}^{t+\mathrm{1}}\right)=\mathrm{0}$, which are solved using Newton's method using the Portable, Extensible Toolkit for Scientific Computing (PETSc) library. The algorithm of Newton's method requires solution of the following linear problem
where ${\mathbf{J}}^{t+\mathrm{1},k}\left({\mathbf{P}}^{t+\mathrm{1},k}\right)$ is the Jacobian matrix. In VSFM, the Jacobian matrix is computed analytically. The contribution to the diagonal and offdiagonal entry of the Jacobian matrix from nth residual equations are given by
The derivative of the accumulation term in J_{nn} is computed as
The derivative of flux between nth and n^{′}th control volume with respect to pressure of each control volume is given as
Lastly, the derivative of the Darcy velocity between the nth and n^{′}th control volume with respect to pressure of each control volume is given as
VSFM uses a twostage check to determine an acceptable numerical solution

Stage1: at any temporal integration stage, the model attempts to solve the set of nonlinear equations given by Eq. (19) with a given time step. If the model fails to find a solution to the nonlinear equations with a given error tolerance settings, the time step is reduced by half and the model again attempts to solve the nonlinear problem. If the model fails to find a solution after a maximum number of time step cuts (currently 20), the model reports an error and stops execution. None of the simulations reported in this paper failed this check.

Stage2: after a numerical solution for the nonlinear problem is obtained, a mass balance error is calculated as the difference between input and output fluxes and change in mass over the integration time step. If the mass balance error exceeds 10–5 kg m^{−2}, the error tolerances for the nonlinear problem are tightened by a factor of 10 and the model reenters Stage1. If the model fails to find a solution with an acceptable mass balance error after 10 attempts of tightening error tolerances, the model reports an error and stops execution. None of the simulations reported in this paper failed this check.
The standalone VSFM code is available at https://github.com/MPPLSM/MPP (last access: 25 September 2018) and the DOI for the code is https://doi.org/10.5281/zenodo.1412949 (last access: 25 September 2018). Notes on how to run the VSFM for all benchmark problems and compare results against PFLOTRAN at https://bitbucket.org/gbisht/notesforgmd201844 (last access: 25 September 2018). The research was performed using E3SM v1.0 and the code is available at https://github.com/E3SMProject/E3SM (last access: 25 September 2018).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1140852018supplement.
GB developed and integrated VSFM in ELM v1.0. GB and WJR designed the study and prepared the manuscript. GEH provided a template for VSFM development. DML implemented smooth approximations of water retention curves in VSFM.
The authors declare that they have no conflict of interest.
This research was supported by the Director, Office of Science, Office of
Biological and Environmental Research of the US Department of Energy under
contract no. DEAC0205CH11231 as part of the Energy Exascale Earth System
Model (E3SM) programs.
Edited by: Jatin Kala
Reviewed by: Minki Hong and two anonymous referees
Alkhaier, F., Flerchinger, G. N., and Su, Z.: Shallow groundwater effect on land surface temperature and surface energy balance under bare soil conditions: modeling and description, Hydrol. Earth Syst. Sci., 16, 1817–1831, https://doi.org/10.5194/hess1618172012, 2012.
Alley, W. M.: Ground Water and Climate, Ground Water, 39, 161–161, 2001.
Amenu, G. G. and Kumar, P.: A model for hydraulic redistribution incorporating coupled soilroot moisture transport, Hydrol. Earth Syst. Sci., 12, 55–74, https://doi.org/10.5194/hess12552008, 2008.
Anyah, R. O., Weaver, C. P., MiguezMacho, G., Fan, Y., and Robock, A.: Incorporating water table dynamics in climate modeling: 3. Simulated groundwater influence on coupled landatmosphere variability, J. Geophys. Res.Atmos., 113, D07103, https://doi.org/10.1029/2007JD009087, 2008.
Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W. D., Kaushik, D., Knepley, M. G., McInnes, L. C., Rupp, K., Smith, B. F., Zampini, S., Zhang, H., and Zhang, H.: PETSc Users Manual, Argonne National LaboratoryANL95/11 – Revision, 3.7, 1–241 pp., 2016.
Banks, E. W., Brunner, P., and Simmons, C. T.: Vegetation controls on variably saturated processes between surface water and groundwater and their impact on the state of connection, Water Resour. Res., 47, W11517, https://doi.org/10.1029/2011WR010544, 2011.
Bense, V. F., Kooi, H., Ferguson, G., and Read, T.: Permafrost degradation as a control on hydrogeological regime shifts in a warming climate, J. Geophys. Res.Earth Surf., 117, F03036, https://doi.org/10.1029/2011JF002143, 2012.
Bernhardt, M., Schulz, K., Liston, G. E., and Zängl, G.: The influence of lateral snow redistribution processes on snow melt and sublimation in alpine regions, J. Hydrol., 424–425, 196–206, 2012.
Beven, K. J. and Kirkby, M. J.: A physically based, variable contributing area model of basin hydrology/Un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant, Hydrol. Sci. Bull., 24, 43–69, 1979.
Bisht, G., Huang, M., Zhou, T., Chen, X., Dai, H., Hammond, G. E., Riley, W. J., Downs, J. L., Liu, Y., and Zachara, J. M.: Coupling a threedimensional subsurface flow and transport model with a land surface model to simulate streamaquiferland interactions (CP v1.0), Geosci. Model Dev., 10, 4539–4562, https://doi.org/10.5194/gmd1045392017, 2017.
Bisht, G., Riley, W. J., Wainwright, H. M., Dafflon, B., Yuan, F., and Romanovsky, V. E.: Impacts of microtopographic snow redistribution and lateral subsurface processes on hydrologic and thermal states in an Arctic polygonal ground ecosystem: a case study using ELM3D v1.0, Geosci. Model Dev., 11, 61–76, https://doi.org/10.5194/gmd11612018, 2018.
Bohrer, G., Mourad, H., Laursen, T. A., Drewry, D., Avissar, R., Poggi, D., Oren, R., and Katul, G. G.: Finite element tree crown hydrodynamics model (FETCH) using porous media flow within branching elements: A new representation of tree hydrodynamics, Water Resour. Res., 41, W11404, https://doi.org/10.1029/2005WR004181, 2005.
Brooks, R. H. and Corey, A. T.: Hydraulic properties of porous media, Colorado State University, Fort Collins, CO, 1964.
Brown, J., Knepley, M. G., May, D. A., McInnes, L. C., and Smith, B.: Composable Linear Solvers for Multiphysics, 2012 11th International Symposium on Parallel and Distributed Computing, Munich, 55–62, https://doi.org/10.1109/ISPDC.2012.16, 2012.
Brunke, M. A., Broxton, P., Pelletier, J., Gochis, D., Hazenberg, P., Lawrence, D. M., Leung, L. R., Niu, G.Y., Troch, P. A., and Zeng, X.: Implementing and Evaluating Variable Soil Thickness in the Community Land Model, Version 4.5 (CLM4.5), J. Climate, 29, 3441–3461, 2016.
Celia, M. A., Bouloutas, E. T., and Zarba, R. L.: A general massconservative numerical solution for the unsaturated flow equation, Water Resour. Res., 26, 1483–1496, 1990.
Chen, J. and Kumar, P.: Topographic Influence on the Seasonal and Interannual Variation of Water and Energy Balance of Basins in North America, J. Climate, 14, 1989–2014, 2001.
Chen, Y., Chen, Y., Xu, C., Ye, Z., Li, Z., Zhu, C., and Ma, X.: Effects of ecological water conveyance on groundwater dynamics and riparian vegetation in the lower reaches of Tarim River, China, Hydrol. Process., 24, 170–177, 2010.
Chen, X. and Hu, Q.: Groundwater influences on soil moisture and surface evaporation, J. Hydrol., 297, 285–300, 2004.
Clapp, R. B. and Hornberger, G. M.: Empirical equations for some soil hydraulic properties, Water Resour. Res., 14, 601–604, 1978.
Clark, M. P., Fan, Y., Lawrence, D. M., Adam, J. C., Bolster, D., Gochis, D. J., Hooper, R. P., Kumar, M., Leung, L. R., Mackay, D. S., Maxwell, R. M., Shen, C., Swenson, S. C., and Zeng, X.: Improving the representation of hydrologic processes in Earth System Models, Water Resour. Res., 51, 5929–5956, 2015.
Collier, N., Hoffman, F. M., Lawwrence, D. M., KeppelAleks, G., Koven, C. D., Riley, W. J., Mu, M., and Randerson, J. T.: The International Land 1 Model Benchmarking (ILAMB) System: Design, Theory, and Implementation, J. Adv. Model. Earth Syst., in review, 2018.
Dai, A. and Trenberth, K. E.: Estimates of Freshwater Discharge from Continents: Latitudinal and Seasonal Variations, J. Hydrometeorol., 3, 660–687, 2002.
Dams, J., Woldeamlak, S. T., and Batelaan, O.: Predicting landuse change and its impact on the groundwater system of the Kleine Nete catchment, Belgium, Hydrol. Earth Syst. Sci., 12, 1369–1385, https://doi.org/10.5194/hess1213692008, 2008.
Dennis, J. M., Vertenstein, M., Worley, P. H., Mirin, A. A., Craig, A. P., Jacob, R., and Mickelson, S.: Computational performance of ultrahighresolution capability in the Community Earth System Model, The Int. J. High Perform. Comput. Appl., 26, 5–16, 2012.
E3SM Project: DOE, Energy Exascale Earth System Model, Computer Software, available at: https://github.com/E3SMProject/E3SM.git, last access: 23 April 2018.
Fan, Y., MiguezMacho, G., Weaver, C. P., Walko, R., and Robock, A.: Incorporating water table dynamics in climate modeling: 1. Water table observations and equilibrium water table simulations, J. Geophys. Res.Atmos., 112, D10125, https://doi.org/10.1029/2006JD008111, 2007.
Fan, Y., Li, H., and MiguezMacho, G.: Global Patterns of Groundwater Table Depth, Science, 339, 940–943, 2013.
Farthing, M. W., Kees, C. E., and Miller, C. T.: Mixed finite element methods and higher order temporal approximations for variably saturated groundwater flow, Adv. Water Resour., 26, 373–394, 2003.
Ferguson, I. M. and Maxwell, R. M.: Human impacts on terrestrial hydrology: climate change versus pumping and irrigation, Environ. Res. Lett., 7, 044022, https://doi.org/10.1088/17489326/7/4/044022, 2012.
Frampton, A., Painter, S., Lyon, S. W., and Destouni, G.: Nonisothermal, threephase simulations of nearsurface flows in a model permafrost system under seasonal variability and climate change, J. Hydrol., 403, 352–359, 2011.
Ghimire, B., Riley, W. J., Koven, C. D., Mu, M., and Randerson, J. T.: Representing leaf and root physiological traits in CLM improves global carbon and nitrogen cycling predictions, J. Adv. Model. Earth Syst., 8, 598–613, 2016.
Grant, R. F., Humphreys, E. R., and Lafleur, P. M.: Ecosystem CO_{2} and CH_{4} exchange in a mixed tundra and a fen within a hydrologically diverse Arctic landscape: 1. Modeling versus measurements, J. Geophys. Res.Biogeosci., 120, 1366–1387, 2015.
Grant, R. F., Mekonnen, Z. A., Riley, W. J., Wainwright, H. M., Graham, D., and Torn, M. S.: Mathematical Modelling of Arctic Polygonal Tundra with Ecosys: 1. Microtopography Determines How Active Layer Depths Respond to Changes in Temperature and Precipitation, J. Geophys. Res.Biogeosci., 122, 3161–3173, 2017.
Green, T. R., Taniguchi, M., Kooi, H., Gurdak, J. J., Allen, D. M., Hiscock, K. M., Treidel, H., and Aureli, A.: Beneath the surface of global change: Impacts of climate change on groundwater, J. Hydrol., 405, 532–560, 2011.
Gutmann, E. D. and Small, E. E.: The effect of soil hydraulic properties vs. soil texture in land surface models, Geophys. Res. Lett., 32, L02402, https://doi.org/10.1029/2004GL021843, 2005.
Hammond, G. E. and Lichtner, P. C.: Fieldscale model for the natural attenuation of uranium at the Hanford 300 Area using highperformance computing, Water Resour. Res., 46, W09527, https://doi.org/10.1029/2009WR008819, 2010.
Hilberts, A. G. J., Troch, P. A., and Paniconi, C.: Storagedependent drainable porosity for complex hillslopes, Water Resour. Res., 41, W06001, https://doi.org/10.1029/2004WR003725, 2005.
Hoffman, F. M., Koven, C. D., KeppelAleks, G., Lawrence, D. M., Riley, W. J., Randerson, J. T., Ahlstrom, A., Abramowitz, G., Baldocchi, D. D., Best, M. J., BondLamberty, B., Kauwe, M. G. D., Denning, A. S., Desai, A. R., Eyring, V., Fisher, J. B., Fisher, R. A., Gleckler, P. J., Huang, M., Hugelius, G., Jain, A. K., Kiang, N. Y., Kim, H., Koster, R. D., Kumar, S. V., Li, H., Luo, Y., Mao, J., McDowell, N. G., Mishra, U., Moorcroft, P. R., Pau, G. S. H., Ricciuto, D. M., Schaefer, K., Schwalm, C. R., Serbin, S. P., Shevliakova, E., Slater, A. G., Tang, J., Williams, M., Xia, J., Xu, C., Joseph, R., and Koch, D.: International Land Model Benchmarking (ILAMB) 2016 Workshop Report, U.S. Department of Energy, Office of Science, 159 pp., 2017.
Hou, Z., Huang, M., Leung, L. R., Lin, G., and Ricciuto, D. M.: Sensitivity of surface flux simulations to hydrologic parameters based on an uncertainty quantification framework applied to the Community Land Model, J. Geophys. Res.Atmos., 117, D15108, https://doi.org/10.1029/2012JD017521, 2012.
Hwang, T., Band, L. E., Vose, J. M., and Tague, C.: Ecosystem processes at the watershed scale: Hydrologic vegetation gradient as an indicator for lateral hydrologic connectivity of headwater catchments, Water Resour. Res., 48, W06514, https://doi.org/10.1029/2011WR011301, 2012.
Ji, P., Yuan, X., and Liang, X.Z.: Do Lateral Flows Matter for the Hyperresolution Land Surface Modeling?, J. Geophys. Res.Atmos., 12077–12092, https://doi.org/10.1002/2017JD027366, 2017.
Jiang, X., Niu, G.Y., and Yang, Z.L.: Impacts of vegetation and groundwater dynamics on warm season precipitation over the Central United States, J. Geophys. Res.Atmos., 114, D06109, https://doi.org/10.1029/2008JD010756, 2009.
Jung, M., Reichstein, M., and Bondeau, A.: Towards global empirical upscaling of FLUXNET eddy covariance observations: validation of a model tree ensemble approach using a biosphere model, Biogeosciences, 6, 2001–2013, https://doi.org/10.5194/bg620012009, 2009.
Kane, D. L., Hinkel, K. M., Goering, D. J., Hinzman, L. D., and Outcalt, S. I.: Nonconductive heat transfer associated with frozen soils, Global Planet. Change, 29, 275–292, 2001.
Kees, C. E. and Miller, C. T.: Higher order time integration methods for twophase flow, Adv. Water Resour., 25, 159–177, 2002.
Kim, H., Yeh, P. J. F., Oki, T., and Kanae, S.: Role of rivers in the seasonal variations of terrestrial water storage over global basins, Geophys. Res. Lett., 36, L17402, https://doi.org/10.1029/2009GL039006, 2009.
Kollet, S. J. and Maxwell, R. M.: Capturing the influence of groundwater dynamics on land surface processes using an integrated, distributed watershed model, Water Resour. Res., 44, W02402, https://doi.org/10.1029/2007WR006004, 2008.
Koster, R. D., Suarez, M. J., Ducharne, A., Stieglitz, M., and Kumar, P.: A catchmentbased approach to modeling land surface processes in a general circulation model: 1. Model structure, J. Geophys. Res.Atmos., 105, 24809–24822, 2000.
Kundzewicz, Z. W. and Doli, P.: Will groundwater ease freshwater stress under climate change?, Hydrol. Sci. J., 54, 665–675, 2009.
Lasslop, G., Reichstein, M., Papale, D., Richardson, A. D., Arneth, A., Barr, A., Stoy, P., and Wohlfahrt, G.: Separation of net ecosystem exchange into assimilation and respiration using a light response curve approach: critical issues and global evaluation, Global Change Biol., 16, 187–208, 2010.
Leng, G., Huang, M., Tang, Q., and Leung, L. R.: A modeling study of irrigation effects on global surface water and groundwater resources under a changing climate, J. Adv. Model. Earth Syst., 7, 1285–1304, 2015.
Leng, G., Leung, L. R., and Huang, M.: Significant impacts of irrigation water sources and methods on modeling irrigation effects in the ACMELand Model, J. Adv. Model. Earth Syst., 9, 1665–1683, 2017.
Leung, L. R., Huang, M., Qian, Y., and Liang, X.: Climate–soil–vegetation control on groundwater table dynamics and its feedbacks in a climate model, Clim. Dynam., 36, 57–81, 2011.
Levine, J. B. and Salvucci, G. D.: Equilibrium analysis of groundwater–vadose zone interactions and the resulting spatial distribution of hydrologic fluxes across a Canadian Prairie, Water Resour. Res., 35, 1369–1383, 1999.
Liang, X., Xie, Z., and Huang, M.: A new parameterization for surface and groundwater interactions and its impact on water budgets with the variable infiltration capacity (VIC) land surface model, J. Geophys. Res.Atmos., 108, 8613, https://doi.org/10.1029/2002JD003090, 2003.
Lohse, K. A., Brooks, P. D., McIntosh, J. C., Meixner, T., and Huxman3, T. E.: Interactions Between Biogeochemistry and Hydrologic Systems, Ann. Rev. Environ. Resour., 34, 65–96, 2009.
Manoli, G., Bonetti, S., Domec, J.C., Putti, M., Katul, G., and Marani, M.: Tree root systems competing for soil moisture in a 3D soil–plant model, Adv. Water Resour., 66, 32–42, 2014.
Marvel, K., Biasutti, M., Bonfils, C., Taylor, K. E., Kushnir, Y., and Cook, B. I.: Observed and Projected Changes to the Precipitation Annual Cycle, J. Climate, 30, 4983–4995, 2017.
Maxwell, R. M. and Miller, N. L.: Development of a Coupled Land Surface and Groundwater Model, J. Hydrometeorol., 6, 233–247, 2005.
McDowell, N. G. and Allen, C. D.: Darcy's law predicts widespread forest mortality under climate warming, Nat. Clim. Change, 5, 669–672, 2015.
MiguezMacho, G., Fan, Y., Weaver, C. P., Walko, R., and Robock, A.: Incorporating water table dynamics in climate modeling: 2. Formulation, validation, and soil moisture simulation, J. Geophys. Res.Atmos., 112, D13108, https://doi.org/10.1029/2006JD008112, 2007.
Milly, P. C. D., Malyshev, S. L., Shevliakova, E., Dunne, K. A., Findell, K. L., Gleeson, T., Liang, Z., Phillipps, P., Stouffer, R. J., and Swenson, S.: An Enhanced Model of Land Water and Energy for Global Hydrologic and EarthSystem Studies, J. Hydrometeorol., 15, 1739–1761, 2014.
Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522, 1976.
Niu, G.Y., Yang, Z.L., Dickinson, R. E., and Gulden, L. E.: A simple TOPMODELbased runoff parameterization (SIMTOP) for use in global climate models, J. Geophys. Res.Atmos., 110, D21106, https://doi.org/10.1029/2005JD006111, 2005.
Niu, G.Y., Yang, Z.L., Dickinson, R. E., Gulden, L. E., and Su, H.: Development of a simple groundwater model for use in climate models and evaluation with Gravity Recovery and Climate Experiment data, J. Geophys. Res.Atmos., 112, D07103, https://doi.org/10.1029/2006JD007522, 2007.
Niu, J., Shen, C., Chambers, J. Q., Melack, J. M., and Riley, W. J.: Interannual Variation in Hydrologic Budgets in an Amazonian Watershed with a Coupled Subsurface–Land Surface Process Model, J. Hydrometeorol., 18, 2597–2617, 2017.
Oleson, K. W., Lawrence, D. M., Bonan, G. B., Drewniak, B., Huang, M., Koven, C. D., Levis, S., Li, F., Riley, W. J., Subin, Z. M., Swenson, S. C., Thornton, P. E., Bozbiyik, A., Fisher, R., Kluzek, E., Lamarque, J.F., Lawrence, P. J., Leung, L. R., Lipscomb, W., Muszala, S., Ricciuto, D. M., Sacks, W., Sun, Y., Tang, J., and Yang, Z.L.: Technical Description of version 4.5 of the Community Land Model (CLM), National Center for Atmospheric Research, Boulder, CO, 422 pp., 2013.
Pacific, V. J., McGlynn, B. L., RiverosIregui, D. A., Welsch, D. L., and Epstein, H. E.: Landscape structure, groundwater dynamics, and soil water content influence soil respiration across riparian–hillslope transitions in the Tenderfoot Creek Experimental Forest, Montana, Hydrol. Process., 25, 811–827, 2011.
Pelletier, J. D., Broxton, P. D., Hazenberg, P., Zeng, X., Troch, P. A., Niu, G.Y., Williams, Z., Brunke, M. A., and Gochis, D.: A gridded global data set of soil, intact regolith, and sedimentary deposit thicknesses for regional and global land surface modeling, J. Adv. Model. Earth Syst., 8, 41–65, 2016.
Petra, D.: Vulnerability to the impact of climate change on renewable groundwater resources: a globalscale assessment, Environ. Res. Lett., 4, 035006, https://doi.org/10.1088/17489326/4/3/035006, 2009.
Piao, S. L., Ito, A., Li, S. G., Huang, Y., Ciais, P., Wang, X. H., Peng, S. S., Nan, H. J., Zhao, C., Ahlström, A., Andres, R. J., Chevallier, F., Fang, J. Y., Hartmann, J., Huntingford, C., Jeong, S., Levis, S., Levy, P. E., Li, J. S., Lomas, M. R., Mao, J. F., Mayorga, E., Mohammat, A., Muraoka, H., Peng, C. H., Peylin, P., Poulter, B., Shen, Z. H., Shi, X., Sitch, S., Tao, S., Tian, H. Q., Wu, X. P., Xu, M., Yu, G. R., Viovy, N., Zaehle, S., Zeng, N., and Zhu, B.: The carbon budget of terrestrial ecosystems in East Asia over the last two decades, Biogeosciences, 9, 3571–3586, https://doi.org/10.5194/bg935712012, 2012.
Pruess, K., Oldenburg, C., and Moridis, G.: TOUGH2 User's Guide, Version 2.0, Lawrence Berkeley National Laboratory, Berkeley, CALBNL43134, 1999.
Rihani, J. F., Maxwell, R. M., and Chow, F. K.: Coupling groundwater and land surface processes: Idealized simulations to identify effects of terrain and subsurface heterogeneity on land surface energy fluxes, Water Resour. Res., 46, W12523, https://doi.org/10.1029/2010WR009111, 2010.
Salvucci, G. D. and Entekhabi, D.: Hillslope and Climatic Controls on Hydrologic Fluxes, Water Resour. Res., 31, 1725–1739, 1995.
Shen, C., Niu, J., and Phanikumar, M. S.: Evaluating controls on coupled hydrologic and vegetation dynamics in a humid continental climate watershed using a subsurfaceland surface processes model, Water Resour. Res., 49, 2552–2572, 2013.
Siebert, S., Burke, J., Faures, J. M., Frenken, K., Hoogeveen, J., Döll, P., and Portmann, F. T.: Groundwater use for irrigation – a global inventory, Hydrol. Earth Syst. Sci., 14, 1863–1880, https://doi.org/10.5194/hess1418632010, 2010.
Sivapalan, M., Beven, K., and Wood, E. F.: On hydrologic similarity: 2. A scaled model of storm runoff production, Water Resour. Res., 23, 2266–2278, 1987.
Soylu, M. E., Istanbulluoglu, E., Lenters, J. D., and Wang, T.: Quantifying the impact of groundwater depth on evapotranspiration in a semiarid grassland region, Hydrol. Earth Syst. Sci., 15, 787–806, https://doi.org/10.5194/hess157872011, 2011.
Sperry, J. S., Adler, F. R., Campbell, G. S., and Comstock, J. P.: Limitation of plant water use by rhizosphere and xylem conductance: results from a model, Plant, Cell Environ., 21, 347–359, 1998.
Srivastava, R. and Yeh, T. C. J.: Analytical solutions for onedimensional, transient infiltration toward the water table in homogeneous and layered soils, Water Resour. Res., 27, 753–762, 1991.
Swenson, S. C. and Lawrence, D. M.: Assessing a dry surface layerbased soil resistance parameterization for the Community Land Model using GRACE and FLUXNETMTE data, J. Geophys. Res.Atmos., 119, 10299–210312, 2014.
Swenson, S. C., Lawrence, D. M., and Lee, H.: Improved simulation of the terrestrial hydrological cycle in permafrost regions by the Community Land Model, J. Adv. Model. Earth Syst., 4, M08002, https://doi.org/10.1029/2012MS000165, 2012.
Tanaka, M., Girard, G., Davis, R., Peuto, A., and Bignell, N.: Recommended table for the density of water between 0 ^{∘}C and 40 ^{∘}C based on recent experimental reports, Metrologia, 38, 301, 2001.
Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An Overview of CMIP5 and the Experiment Design, B. Am. Meteorol. Soc., 93, 485–498, 2012.
Taylor, R. G., Scanlon, B., Döll, P., Rodell, M., Van Beek, R., Wada, Y., Longuevergne, L., Leblanc, M., Famiglietti, J. S., and Edmunds, M.: Ground water and climate change, Nat. Clim. Change, 3, 322–329, 2013.
Tian, W., Li, X., Cheng, G.D., Wang, X.S., and Hu, B. X.: Coupling a groundwater model with a land surface model to improve water and energy cycle simulation, Hydrol. Earth Syst. Sci., 16, 4707–4723, https://doi.org/10.5194/hess1647072012, 2012.
van Genuchten, M. T.: A Closedform Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1, Soil Sci. Soc. Am. J., 44, 892–898, 1980.
Walko, R. L., Band, L. E., Baron, J., Kittel, T. G. F., Lammers, R., Lee, T. J., Ojima, D., Sr., R. A. P., Taylor, C., Tague, C., Tremback, C. J., and Vidale, P. L.: Coupled AtmosphereBiophysicsHydrology Models for Environmental Modeling, J. Appl. Meteorol., 39, 931–944, 2000.
White, M. and Stomp, O. M.: Subsurface transport over multiple phases; Version 2.0; Theory Guide, Pacific Northwest National Laboratory, 2000.
Yeh, P. J.F. and Eltahir, E. A. B.: Representation of Water Table Dynamics in a Land Surface Scheme. Part I: Model Development, J. Climate, 18, 1861–1880, 2005.
York, J. P., Person, M., Gutowski, W. J., and Winter, T. C.: Putting aquifers into atmospheric simulation models: an example from the Mill Creek Watershed, northeastern Kansas, Adv. Water Resour., 25, 221–238, 2002.
Yuan, X., Xie, Z., Zheng, J., Tian, X., and Yang, Z.: Effects of water table dynamics on regional climate: A case study over east Asian monsoon area, J. Geophys. Res.Atmos., 113, D21112, https://doi.org/10.1029/2008JD010180, 2008.
Zektser, I. S. and Evertt, L. G.: Groundwater resources of the world and their use, United Nations Educational, Scientific and Cultural Organization7, place de Fontenoy, 75352 Paris 07 SP, 2004.
Zeng, X. and Decker, M.: Improving the Numerical Solution of Soil Moisture–Based Richards Equation for Land Models with a Deep or Shallow Water Table, J. Hydrometeorol., 10, 308–319, 2009.
Zhu, Q., Riley, W. J., Tang, J., and Koven, C. D.: Multiple soil nutrient competition between plants, microbes, and mineral surfaces: model development, parameterization, and example applications in several tropical forests, Biogeosciences, 13, 341–363, 2016.
 Abstract
 Introduction
 Methods
 Results and discussion
 Summary and conclusion
 Appendix A: Smooth approximation of Brooks–Corey water retention curve
 Appendix B: Residual equation of VSFM formulation
 Appendix C: Jacobian equation of VSFM formulation
 Appendix D: Numerical checks in VSFM
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 References
 Supplement
 Abstract
 Introduction
 Methods
 Results and discussion
 Summary and conclusion
 Appendix A: Smooth approximation of Brooks–Corey water retention curve
 Appendix B: Residual equation of VSFM formulation
 Appendix C: Jacobian equation of VSFM formulation
 Appendix D: Numerical checks in VSFM
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 References
 Supplement