Articles | Volume 11, issue 10
Model description paper
11 Oct 2018
Model description paper |  | 11 Oct 2018

Development and evaluation of a variably saturated flow model in the global E3SM Land Model (ELM) version 1.0

Gautam Bisht, William J. Riley, Glenn E. Hammond, and David M. Lorenzetti

Improving global-scale model representations of near-surface 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 ELMv1-VSFM's subsurface drainage parameter, fd, was calibrated to match an observationally constrained and spatially explicit global water table depth (WTD) product. Optimal spatially explicit fd 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 fd 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.

1 Introduction

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, high-intensity 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 land-use (LU) and land-cover (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 groundwater-fed 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. Type-1 models assume a quasi-steady state equilibrium of the soil moisture profile above the water table (Hilberts et al., 2005; Koster et al., 2000; Walko et al., 2000). Type-2 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 one-dimensional Type-2 models include Liang et al. (2003), Yeh and Eltahir (2005), Niu et al. (2007), and Zeng and Decker (2009). Examples of quasi-three-dimensional Type-2 models are York et al. (2002), Fan et al. (2007), Miguez-Macho et al. (2007), and Shen et al. (2013). Type-3 models include a three-dimensional representation of subsurface flow based on the variably saturated Richards equation (Maxwell and Miller, 2005; Tian et al., 2012). Type-3 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 (Type-4) 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. Type-4 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., Type-3 and Type-4 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 Type-2 subsurface hydrology model based on Zeng and Decker (2009). In this work, we developed in ELMv1 the Type-4 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, three-dimensional subsurface model, the VSFM is a serial, one-dimensional model that is appropriate for climate-scale 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 Methods

2.1 Current model formulation

Water flow in the unsaturated zone is often described by the θ-based Richards equation:

(1) θ t = - q - Q ,

where θ (m3 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 (m3 of water m−3 of soil s−1) is a soil moisture sink term. The Darcy flux, q, is given by

(2) q = - K ( ψ + z ) ,

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 non-linear functions of volumetric soil moisture following Clapp and Hornberger (1978):


where Ksat (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

(5) q = - K ( ψ + z - C ) ,

where C is a constant hydraulic potential above the water table, z, given as

(6) C = ψ E + z = ψ sat θ E ( z ) θ sat - B + z = ψ sat + z ,

where ψE (m) is the equilibrium soil matric potential, z (m) is height above a reference datum, θE (m3 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 cell-centered 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 tri-diagonal 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 zero-flux 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):

(7) ( ϕ s w ρ ) t = - ( ρ q ) - Q ,

where ϕ (m3 m−3) is the soil porosity, sw (-) 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 one-dimensional system and the flow velocity is defined by Darcy's law:

(8) q = - k k r μ ( P + ρ g z ) ,

where k (m2) is intrinsic permeability, kr (-) 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 (se)

(9) s e = s w - s r 1 - s r ,

where sw is saturation and sr is residual saturation. We have implemented Brooks and Corey (1964; Eq. 10) and van Genuchten (1980; Eq. 11) WRCs:


where Pc (Pa) is the capillary pressure, Pc0 (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 Pc=P-Pref where Pref is Pc0 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:

(12) κ r s e = s e 0.5 1 - 1 - s e 1 / m m if P < P ref , 1 if P P ref ,

where m=1-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):




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 zero-flux 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 cell-centered finite volume discretization to decompose the spatial domain, Ω, into N non-overlapping control volumes, Ωn, such that Ω=n=1NΩi and Γn represents the boundary of the nth control volume. Applying a finite volume integral to Eq. (7) and the divergence theorem yields

(14) t Ω n ϕ s w ρ d V = - Γ n ρ q d A - Ω n Q d V .

The discretized form of the left-hand-side term and first term on the right-hand side of Eq. (14) are approximated as


where Ann (m2) is the common face area between the nth and nth control volumes. After substituting Eqs. (15) and (16) in Eq. (14), the resulting ordinary differential equation for the variably saturated flow model is

(17) d d t ϕ s w ρ V n = - n ρ q n n A n n - Q n V n .

We perform temporal integration of Eq. (17) using the backward-Euler 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, multi-component, multi-physics 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.

Table 1Soil properties and discretization used in the three test problems described in Sect. 2.3.

Download Print Version | Download XLSX

2.3 VSFM single-column evaluation

We tested the VSFM with three idealized 1-dimensional test problems. First, the widely studied problem for the 1-D 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 one-dimensional vertical infiltration in a two-layered 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 (=1.01325×105 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 1-dimensional 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 2.5×10-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 ELMv1-VSFM at a spatial resolution of 1.9 (latitude) × 2.5 (longitude) with a 30 (min) time step for 200 years, including a 180-year spin-up 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 ELMv1-VSFM spatial resolution. ELM-VSFM'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 ELMv1-VSFM 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 long-term ground water depth is the subsurface drainage flux, qd (kg m−2 s−1), which is calculated based on water table depth, z (m), Niu et al. (2005):

(20) q d = q d , max exp - f d z ,

where qd,max (kg m−2 s−1) is the maximum drainage flux that depends on grid cell slope and fd (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 qd,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 qd,max as a single calibration parameter. ELMv0 uses fd=2.5 m−1 as a global constant and estimates maximum drainage flux when WTD is at the surface as qd,max=10sinβ kg m−2 s−1, where β (radians) is the mean grid cell topographic slope. Of the two parameters, fd and qd,max, available for model calibration, we choose to calibrate fd because the uncertainty analysis by Hou et al. (2012) identified it as the most significant hydrologic parameter in CLM4. To improve on the fd parameter values, we performed an ensemble of global simulations with fd 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 fd and “WTD” was developed for each grid cell and then the F2013 dataset was used to estimate an optimal fd for each grid cell.

2.6 Global ELM-VSFM evaluation

With the optimal fd values, we ran a ELM-VSFM simulation using the protocol described above. We then used the ILAMB package to evaluate the ELMv1-VSFM 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 in-situ, regional, and global scales. Since ELM-VSFM 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 site-level measurements from FLUXNET (Lasslop et al., 2010) and globally from FLUXNET-MTE (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 ELMv1-VSFM simulations with default and calibrated fd to ensure improvements in WTD predictions did not degrade model skill for other processes.

Figure 1Comparison of VSFM simulated pressure profile (blue line) against data (red square) reported in Celia et al. (1990) at time = 24 h for infiltration in a dry soil column; Initial pressure condition is shown by green line.


3 Results and discussion

3.1 VSFM single-column 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).

Figure 2Transient liquid pressure simulated for a two-layer soil system by VSFM (solid line) and PFLOTRAN (square) for (a) wetting and (b) drying scenarios.


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 single-column problems are summarized in Table 1 These three evaluation simulations demonstrate the VSFM accurately represents soil moisture dynamics under conditions relevant to ESM-scale prediction.

Figure 3Transient (a) liquid pressure and (b) soil saturation simulated by VSFM (solid line) and PFLOTRAN (square) for the water table dynamics test problem.


3.2 Subsurface drainage parameterization estimation

The simulated nonlinear WTD–fd relationship is a result of the subsurface drainage parameterization flux given by Eq. (20; Fig. 4a and b). For 0.1fd1, the slope of the WTD–fd relationship for all grid cells is log–log linear with a slope of -1.0±0.1. The log–log linear relationship breaks down for fd>1, where the drainage flux becomes much smaller than infiltration and evapotranspiration (Fig. 4c and d). Thus, at larger fd, the steady state z becomes independent of fd and is determined by the balance of infiltration and evapotranspiration.

Figure 4(a–b) The nonlinear relationship between simulated water table depth (WTD) and fd for two grid cells within ELM's global grid. WTD from the Fan et al. (2013) dataset and optimal fd for the two grid cells are shown with dashed red and dashed black lines, respectively. (c–d) The simulated drainage, evapotranspiration, and infiltration fluxes as functions of optimal fd for the two ELM grid cells.


Figure 5Global estimate of fd.


For 79 % of the global grid cells, the ensemble range of simulated WTD spanned the F2013 dataset. The optimal value of fd 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 fd was deeper than that in the F2013 dataset, the optimal fd value was chosen as the one that resulted in the lowest absolute WTD error (e.g., Fig. 4b). At large fd 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 fd value. The optimal fd has a global average of 1.60 m-1±2.68 m−1 and 72 % of global grid cells have an optimal fd value lower than the global average (Fig. 5).

Figure 6(a) Water table depth observation from Fan et al. (2013); (b) water table depth biases (= Model–Obs) from ELMv1-VSFM using default spatially homogeneous fd; and (c) water table depth biases from ELMv1-VSFM using spatially heterogeneous fd.


Figure 7(a) Annual range of water table depth for ELMv1-VSFM simulation with spatially heterogeneous estimates of fd and (b) difference in annual water table depth range between simulations with optimal and default fd.


Table 2Bias, root mean square error (RMSE), and correlation (R2) between simulated water table depth and Fan et al. (2013) data.

Download Print Version | Download XLSX

3.3 Global simulation evaluation

The ELMv1-VSFM predictions are much closer to the F2013 dataset (Fig. 6a) using optimal globally distributed fd values (Fig. 6c) compared to the default fd 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 fd values. In the simulation using optimal globally distributed fd values, all grid cells with WTD bias >3.7 m were those for which an optimal fd was not found. The mean global bias, RMSE, and R2 values improved in the new ELMv1-VSFM compared to the default model (Table 2). The 79 % of global grid cells for which an optimal fd value was estimated had significantly better water table prediction with a bias, RMSE, and R2 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 R2 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 fd values (Fig. 7a). The annual WTD range decreased by 0.24 m for the 79 % of the grid cells for which an optimal fd value was estimated (Fig. 7b).

Globally averaged WTD in ELMv1-VSFM simulations with default fd and optimal fd values were 10.5 and 20.1 m, respectively. Accurate prediction of deep WTD in the simulation with optimal fd caused very small differences in near-surface soil moisture (Fig. 8). The 79 % of grid cells with an optimal fd value had deeper globally averaged WTDs than when using the default fd value (24.3 vs. 8.6 m). For these 79 % of grid cells, the WTD was originally deep enough to not impact near-surface conditions (Kollet and Maxwell, 2008); therefore, further lowering of WTD led to negligible changes in near-surface hydrological conditions.

Figure 8Seasonal monthly mean soil moisture differences for top 10 cm between ELMv1-VSFM simulations with optimal and default fd values.


Table 3ILAMB benchmark scores for latent heat flux (LH), sensible heat flux 640 (SH), terrestrial water storage anomaly (TWSA), and surface runoff. The calculation of ILAMB metrics and scores are described at (last access: 4 October 2018).

Download Print Version | Download XLSX

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 ELMv1-VSFM model (Table 3). Optimal fd values had inconsequential impacts on simulated surface energy fluxes at site-level and global scales. Optimal fd 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 fd 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 fd 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 10-year-long 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 5-year-long 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 ELMv1-VSFM 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 ultrahigh-resolution 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 fd 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 ELMv1-VSFM 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 1-D 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 3-D 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 (Miguez-Macho 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 3-D 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 1-D model. Tradeoffs between approaches to represent lateral processes and computational costs need to be carefully studied before developing quasi- or fully three-dimensional LSMs (Clark et al., 2015).

Transport of water across multiple components of the soil-plant-atmosphere 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 high-latitude 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 multi-component problem (e.g., transport of water through the soil-plant-atmosphere continuum) and multi-physics problem (e.g., fully coupled conservation of mass and energy equations in the subsurface). DMComposite allows for an easy assembly of a tightly coupled multi-physics problem from individual physics formulations (Brown et al., 2012).

4 Summary and conclusion

Starting from the climate-scale 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, fd, 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 fd 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 fd significantly improved global WTD prediction by reducing bias and RMSE and increasing R2. Grids without an optimal fd were the largest contributor to error in WTD predication. ILAMB benchmarks on simulations with default and optimal fd 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 fd values were used to improve WTD prediction.

Appendix A: Smooth approximation of Brooks–Corey water retention curve

The Brooks and Corey (1964) water retention curve of Eq. (10) has a discontinuous derivative at Pc=Pc0. 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(Pc), in the neighborhood of Pc0.

(A1) s e = - α P c - λ if P c P u , B P c if P u < P c < P s , 1 if P s P s ,

where the breakpoints Pu and Ps satisfy Pu<Pc0<Ps0. The smoothing polynomial


introduces four more parameters, whose values follow from continuity. In particular, matching the saturated region requires BPs=b0=1, and a continuous derivative at Pc=Ps requires BPs=b1=0. Similarly, matching the value and derivative at Pc=Pu requires


where Δ=Pu-Ps. Note that PuΔ<0.

In practice, setting Pu too close to Pc0 can produce an unwanted local maximum in the cubic smoothing regime, resulting in se>1. Avoiding this condition requires that B(Pc) increase monotonically from Pc=Pu, where BPc>0, to Pc=Ps, where BPc=0. Thus a satisfactory pair of breakpoints ensures

(A5) B P c = P c - P s 2 b 2 + 3 b 3 P c - P s > 0

throughout PuPc<Ps.

Figure A1The Brooks–Corey water rendition curve for estimating liquid saturation, se, as a function of capillary pressure, Pc, shown in solid black line and smooth approximation of Brooks–Corey (SBC) are shown in dashed line.


Let Pc denote a local extremum of B, so that BPc=0. If PcPs, it follows Pc-Ps=-2b2/(3b3). Rewriting Eq. (22), BPc=Pc-Ps3b3Pc-Pc shows that BPu>0 requires either: (1) b3<0 and Pc<Pu; or (2) b3>0 and Pc>Pu. The first possibility places Pc outside the cubic smoothing regime, and so does not constrain the choice of Pu or Ps. The second possibility allows for an unwanted local extremum at Pu<Pc<Ps. In this case, b3>0 implies b2<0 (since Pc<Ps0). Then since BPc=-2b2, the local extremum is a maximum, resulting in sePc>1.

Given a breakpoint Ps, one strategy for choosing Pu is to guess a value, then check whether the resulting b2 and b3 produces Pu<Pc<Ps. If so, Pu should be made more negative. An alternative strategy is to choose Pu in order the guarantee acceptable values for b2 and b3. One convenient choice forces b2=0. Another picks Pu in order to force b3=0. Both of these reductions (1) ensure B(Pc) has a positive slope throughout the smoothing interval, (2) slightly reduce the computation cost of finding se(Pc) for Pc on the smoothing interval, and (3) significantly reduce the computational cost of inverting the model, in order to find Pc as a function of se.

As shown in Fig. A1, the two reductions differ mainly in that setting b2=0 seems to produce narrower smoothing regions (probably due to the fact that this choice gives zero curvature at Pc=Ps, while b3=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 Pu. 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 Pu<Pc0. In any event, since the result of this calculation may be cached for use throughout the simulation, it need not be particularly efficient.

Appendix B: Residual equation of VSFM formulation

The residual equation for the VSFM formulation at t+1 time level for nth control volume is given by


where ϕ (mm3 mm3) is the soil porosity, sw (-) is saturation, ρ (kg m−3) is water density, nn (m s−1) is the Darcy flow velocity between nth and nth control volumes, Ann (m s−1) is the interface area between nth and nth control volumes, and Q (kg m−3 s−1) is a sink of water. The Darcy velocity is computed as

(B2) q n n = - k k r μ n n P n - P n - ρ n n ( g d n n ) d n + d n n n n ,

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, dn (m) and dn (m) are the distances between centroid of the nth and nth control volume to the common interface between the two control volumes, dnn is a distance vector joining centroid of nth and nth control volume, and nnn is a unit normal vector joining centroid of the nth and nth control volume.

The density at the interface of control volume, ρnn, is computed as the inverse distance weighted average by

(B3) ρ n n = ω n ρ n + ω n ρ n ,

where ωn and ωn are given by

(B4) ω n = d n d n + d n = ( 1 - ω n ) .

The first term on the right-hand side of Eq. (27) is computed as the product of distance-weighted harmonic average of intrinsic permeability, knn, and upwinding of kr/μ(=λ) as

(B5) k k r μ n n = k n n k r μ n n = k n k n ( d n + d n ) k n d n + k n d n ) λ n n ,


(B6) λ n n = k r / μ n if q n n > 0 , k r / μ n otherwise .

By substituting Eqs. (28), (29), and (30) in Eq. (27), we obtain

Appendix C: Jacobian equation of VSFM formulation

The discretized equations of VSFM leads to a system of nonlinear equations given by Rt+1Pt+1=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

(C1) J t + 1 , k ( P t + 1 , k ) Δ P t + 1 , k = - R t + 1 , k ( P t + 1 , k ) ,

where Jt+1,k(Pt+1,k) is the Jacobian matrix. In VSFM, the Jacobian matrix is computed analytically. The contribution to the diagonal and off-diagonal entry of the Jacobian matrix from nth residual equations are given by


The derivative of the accumulation term in Jnn is computed as

(C4) ( ρ ϕ s w ) P n = ϕ s w ρ P n + ρ s w ϕ P n + ρ ϕ s w P n .

The derivative of flux between nth and nth control volume with respect to pressure of each control volume is given as


Lastly, the derivative of the Darcy velocity between the nth and nth control volume with respect to pressure of each control volume is given as

Appendix D: Numerical checks in VSFM

VSFM uses a two-stage check to determine an acceptable numerical solution

  • Stage-1: 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.

  • Stage-2: 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 re-enters Stage-1. 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.

Code availability

The stand-alone VSFM code is available at (last access: 25 September 2018) and the DOI for the code is (last access: 25 September 2018). Notes on how to run the VSFM for all benchmark problems and compare results against PFLOTRAN at (last access: 25 September 2018). The research was performed using E3SM v1.0 and the code is available at (last access: 25 September 2018).


The supplement related to this article is available online at:

Author contributions

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.

Competing interests

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. DE-AC02-05CH11231 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,, 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 soil-root moisture transport, Hydrol. Earth Syst. Sci., 12, 55–74,, 2008. 

Anyah, R. O., Weaver, C. P., Miguez-Macho, G., Fan, Y., and Robock, A.: Incorporating water table dynamics in climate modeling: 3. Simulated groundwater influence on coupled land-atmosphere variability, J. Geophys. Res.-Atmos., 113, D07103,, 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 LaboratoryANL-95/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,, 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,, 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 three-dimensional subsurface flow and transport model with a land surface model to simulate stream-aquifer-land interactions (CP v1.0), Geosci. Model Dev., 10, 4539–4562,, 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 ELM-3D v1.0, Geosci. Model Dev., 11, 61–76,, 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,, 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,, 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 mass-conservative 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., Keppel-Aleks, 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 land-use change and its impact on the groundwater system of the Kleine Nete catchment, Belgium, Hydrol. Earth Syst. Sci., 12, 1369–1385,, 2008. 

Dennis, J. M., Vertenstein, M., Worley, P. H., Mirin, A. A., Craig, A. P., Jacob, R., and Mickelson, S.: Computational performance of ultra-high-resolution 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:, last access: 23 April 2018. 

Fan, Y., Miguez-Macho, 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,, 2007. 

Fan, Y., Li, H., and Miguez-Macho, 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,, 2012. 

Frampton, A., Painter, S., Lyon, S. W., and Destouni, G.: Non-isothermal, three-phase simulations of near-surface 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 CO2 and CH4 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,, 2005. 

Hammond, G. E. and Lichtner, P. C.: Field-scale model for the natural attenuation of uranium at the Hanford 300 Area using high-performance computing, Water Resour. Res., 46, W09527,, 2010. 

Hilberts, A. G. J., Troch, P. A., and Paniconi, C.: Storage-dependent drainable porosity for complex hillslopes, Water Resour. Res., 41, W06001,, 2005. 

Hoffman, F. M., Koven, C. D., Keppel-Aleks, G., Lawrence, D. M., Riley, W. J., Randerson, J. T., Ahlstrom, A., Abramowitz, G., Baldocchi, D. D., Best, M. J., Bond-Lamberty, 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,, 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,, 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,, 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,, 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,, 2009. 

Kane, D. L., Hinkel, K. M., Goering, D. J., Hinzman, L. D., and Outcalt, S. I.: Non-conductive 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 two-phase 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,, 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,, 2008. 

Koster, R. D., Suarez, M. J., Ducharne, A., Stieglitz, M., and Kumar, P.: A catchment-based 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,, 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. 

Miguez-Macho, 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,, 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 Earth-System 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 TOPMODEL-based runoff parameterization (SIMTOP) for use in global climate models, J. Geophys. Res.-Atmos., 110, D21106,, 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,, 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., Riveros-Iregui, 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 global-scale assessment, Environ. Res. Lett., 4, 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,, 2012. 

Pruess, K., Oldenburg, C., and Moridis, G.: TOUGH2 User's Guide, Version 2.0, Lawrence Berkeley National Laboratory, Berkeley, CALBNL-43134, 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,, 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 subsurface-land 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,, 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 semi-arid grassland region, Hydrol. Earth Syst. Sci., 15, 787–806,, 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 one-dimensional, 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 layer-based soil resistance parameterization for the Community Land Model using GRACE and FLUXNET-MTE 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,, 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,, 2012. 

van Genuchten, M. T.: A Closed-form 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 Atmosphere-Biophysics-Hydrology 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,, 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. 

Short summary
Most existing global land surface models used to study impacts of climate change on water resources routinely use different models for near-surface unsaturated soil and the deeper groundwater table. We developed a model that uses a unified treatment of soil hydrologic processes throughout the entire soil column. Using a calibrated drainage parameter, the new model is able to correctly predict deep water table depth as reported in an observationally constrained global dataset.