Coupling a three-dimensional subsurface flow and transport model with a land surface model to simulate stream-aquifer-land interactions (CP v1.0)

Abstract. A fully coupled three-dimensional surface and subsurface land model is developed and applied to a site along the Columbia River to simulate three-way interactions among river water, groundwater, and land surface processes. The model features the coupling of the Community Land Model version 4.5 (CLM4.5) and a massively parallel multiphysics reactive transport model (PFLOTRAN). The coupled model, named CP v1.0, is applied to a 400 m × 400 m study domain instrumented with groundwater monitoring wells along the Columbia River shoreline. CP v1.0 simulations are performed at three spatial resolutions (i.e., 2, 10, and 20 m) over a 5-year period to evaluate the impact of hydroclimatic conditions and spatial resolution on simulated variables. Results show that the coupled model is capable of simulating groundwater–river-water interactions driven by river stage variability along managed river reaches, which are of global significance as a result of over 30 000 dams constructed worldwide during the past half-century. Our numerical experiments suggest that the land-surface energy partitioning is strongly modulated by groundwater–river-water interactions through expanding the periodically inundated fraction of the riparian zone, and enhancing moisture availability in the vadose zone via capillary rise in response to the river stage change. Meanwhile, CLM4.5 fails to capture the key hydrologic process (i.e., groundwater–river-water exchange) at the site, and consequently simulates drastically different water and energy budgets. Furthermore, spatial resolution is found to significantly impact the accuracy of estimated the mass exchange rates at the boundaries of the aquifer, and it becomes critical when surface and subsurface become more tightly coupled with groundwater table within 6 to 7 meters below the surface. Inclusion of lateral subsurface flow influenced both the surface energy budget and subsurface transport processes as a result of river-water intrusion into the subsurface in response to an elevated river stage that increased soil moisture for evapotranspiration and suppressed available energy for sensible heat in the warm season. The coupled model developed in this study can be used for improving mechanistic understanding of ecosystem functioning and biogeochemical cycling along river corridors under historical and future hydroclimatic changes. The dataset presented in this study can also serve as a good benchmarking case for testing other integrated models.


Abstract.
A fully coupled three-dimensional surface and subsurface land model is developed and applied to a site along the Columbia River to simulate three-way interactions among river water, groundwater, and land surface processes. The model features the coupling of the Community Land Model version 4.5 (CLM4.5) and a massively parallel multiphysics reactive transport model (PFLOTRAN). The coupled model, named CP v1.0, is applied to a 400 m × 400 m study domain instrumented with groundwater monitoring wells along the Columbia River shoreline. CP v1.0 simulations are performed at three spatial resolutions (i.e., 2, 10, and 20 m) over a 5-year period to evaluate the impact of hydroclimatic conditions and spatial resolution on simulated variables. Results show that the coupled model is capable of simulating groundwater-river-water interactions driven by river stage variability along managed river reaches, which are of global significance as a result of over 30 000 dams constructed worldwide during the past half-century. Our numerical experiments suggest that the land-surface energy partitioning is strongly modulated by groundwater-river-water interactions through expanding the periodically inundated fraction of the riparian zone, and enhancing moisture availability in the vadose zone via capillary rise in response to the river stage change. Meanwhile, CLM4.5 fails to capture the key hydrologic process (i.e., groundwater-river-water exchange) at the site, and consequently simulates drastically different water and energy budgets. Furthermore, spatial resolution is found to significantly impact the accuracy of estimated the mass exchange rates at the boundaries of the aquifer, and it becomes critical when surface and subsurface become more tightly coupled with groundwater table within 6 to 7 meters below the surface. Inclusion of lateral subsurface flow influenced both the surface energy budget and subsurface transport processes as a result of river-water intrusion into the subsurface in response to an elevated river stage that increased soil moisture for evapotranspiration and suppressed available energy for sensible heat in the warm season. The coupled model developed in this study can be used for improving mechanistic understanding of ecosystem functioning and biogeochemical cycling along river corridors under historical and future hydroclimatic changes. The dataset presented in this study can also serve as a good benchmarking case for testing other integrated models.
al., 2015), cloud formation (Rahman et al., 2015), and climate (Leung et al., 2011;Taylor et al., 2013). Lateral subsurface processes are fundamentally important on multiple spatial scales, including hill-slope scales (McNamara et al., 2005;Zhang et al., 2011), basin scales in semiarid and arid climates where regional aquifers sustain baseflows in rivers (Schaller and Fan, 2009) and wetlands (Fan and Miguez-Macho, 2011). However, some current-generation land surface models (LSMs) routinely omit explicit lateral subsurface processes (Clark et al., 2015;Kollet and Maxwell, 2008;Nir et al., 2014), while others include them (described below). Observational and modeling studies suggest that groundwater forms an environmental gradient in soil moisture availability by redistributing water that could profoundly shape critical zone evolution on continental to global scales Taylor et al., 2013). The mismatch between observed and simulated evapotranspiration by current LSMs could be explained by the absence of lateral groundwater flow (Maxwell and Condon, 2016).
It has been increasingly recognized that rivers, despite their small aerial extent on the landscape, play important roles in watershed functioning through their connections with groundwater aquifers and riparian zones . The interactions between groundwater and river water prolong physical storage and enhance reactive processing that alter water chemistry, downstream transport of materials and energy, and biogenic gas emissions (Fischer et al., 2005;Harvey and Gooseff, 2015). The Earth system modeling community recognizes such a gap in existing Earth system models and calls for improved representation of biophysical and biogeochemical processes within the terrestrial-aquatic interface (Gaillardet et al., 2014).
Over the past decade, much effort has been expended to include groundwater in LSMs. Groundwater is important to water and energy budgets such as evapotranspiration (ET), latent heat (LH), and sensible heat (SH), but also to biogeochemical processes such as gross primary production, heterotrophic respiration, and nutrient cycling. The lateral convergence of water along the landscape and two-way groundwater-surface-water exchange are identified as the most relevant subsurface processes to large-scale Earth system functioning (see review by Clark et al., 2015). However, the choice of processes, the approaches to represent multiscale structures and heterogeneities, the data and computational demands, etc. all vary greatly among the research groups, even those working on the same land models.
Most of the LSMs reviewed by Clark et al. (2015) do not explicitly account for stream-aquifer-land interactions. For example, the Community Land Model version 4.5 allows for reinfiltration of flooded waters in a highly parameterized way without explicitly linking to groundwater dynamics, and therefore only one-way flow from the aquifer to the stream is simulated (Oleson et al., 2013). The Land-Ecosystem-Atmosphere Feedback model treats river elevation as part of the two-dimensional vertically integrated groundwater flow equation and allows river and floodwater to infiltrate through sediments in the flood plain (Miguez-Macho and Fan, 2012).
In contrast, the fully integrated models, being a small subset of LSMs, explicitly represent the two-way exchange between groundwater aquifers and their adjacent rivers in a spatially resolved fashion. Such models couple a completely integrated hydrology model with a land surface model, so that the surface-water recharge to groundwater by infiltration or intrusion and base flow discharge from groundwater to surface waters can be estimated in a more mechanistic way.
Examples of the integrated models include (1) the coupling between the Common Land Model (CoLM) and a variably saturated groundwater model (ParFlow) (Maxwell and Miller, 2005); (2) the Penn State Integrated Hydrologic Model (PIHM) (Shi et al., 2013); (3) the coupling between the Process-based Adaptive Watershed Simulator (PAWS) and CLM4.5 (Ji et al., 2015;Pau et al., 2016;Riley and Shen, 2014); (4) the coupling between the CATchment HYdrology (CATHY) model and the Noah model with multiple parameterization schemes (Noah-MP) ; and (5) the coupling between CLM3.5 and ParFlow through the Ocean Atmosphere Sea Ice Soil external coupler (OASIS3) in the Terrestrial Systems Modeling Platform (TerrSysMP) (Shrestha et al., 2014;Gebler et al., 2017). The integrated models eliminate the need for parameterizing lateral groundwater flow and allow the interconnected groundwater-surface-water systems to evolve dynamically based on the governing equations and the properties of the physical system. Although such models often require robust numerical solvers on high-performance computing (HPC) facilities to achieve high-resolution, large-extent simulations , they have been increasingly applied for hydrologic prediction and environmental understanding. However, as a result of differences in physical process representations and numerical solution approaches in terms of (1) the coupling between the variably saturated groundwater and surface water flow, (2) representation of surface water flow, and (3) implementation of subsurface heterogeneity in the existing integrated models, significant discrepancies exist in their results when the models were applied to highly nonlinear problems with heterogeneity and complex water table dynamics, while many of the models show good agreement for simpler test cases where traditional runoff generation mechanisms (i.e., saturation and infiltration excess runoff) apply Maxwell et al., 2014).
The developments of the integrated models have enabled scientific explorations of interactions and feedback mechanisms in the aquifer-soil-vegetation-atmosphere continuum using a holistic and physically based approach (Shrestha et al., 2014;Gilbert et al., 2017). Compared to simulations of regional climate models coupled to traditional LSMs, such a physically based approach shows less sensitivity to uncertainty in the subsurface hydraulic characteristics that could propagate from deep subsurface to free troposphere (Keune et al., 2016), while other physical representations (e.g., Geosci. Model Dev., 10, 2017 www.geosci-model-dev.net/10/4539/2017/ parameterizations in evaporation and transpiration, atmospheric boundary layer schemes) could have significant effects on the simulations as well . Therefore, it is of great scientific interest to further develop the integrated models and benchmarks to achieve improved understanding of complex interactions in the fully coupled Earth system. Motivated by the great potentials of using an integrated model to explore Earth system dynamics, the objective of this study is three-fold. First, we aim to document the development of a coupled land surface and subsurface model as a first step toward a new integrated model, featuring the twoway coupling between two highly scalable and state-of-theart open-source codes: CLM4.5 (Oleson et al., 2013) and a reactive transport model PFLOTRAN (Lichtner et al., 2015). The coupled model mechanistically represents the two-way exchange of water and solute mass between aquifers and river, as well as land-atmosphere exchange of water and energy. The coupled model is therefore named as CP v1.0 hereafter. We note that in recent years, efforts have been made to implement carbon-nitrogen decomposition, nitrification, denitrification, and plant uptake from CLM4.5 in the form of a reaction network solved by PFLOTRAN to enable the coupling of biogeochemical processes between the two models (Tang et al., 2016). In addition, although PAWS is coupled to the same version of CLM (i.e., CLM4.5) (Ji et al., 2015;Pau et al., 2016), PFLOTRAN resolves the subsurface in a three-dimensional fashion, while PAWS approximates the three-dimensional Richards equation by dividing the subsurface into an unsaturated domain represented by the one-dimensional Richards Equation coupled with threedimensional saturated groundwater flow equation for subsurface flow, by assuming that there is no horizontal flow in an unsaturated portion of soil, and that lateral flux in saturated portion is evenly distributed.
Second, we describe a numerically challenging benchmarking case for verifying coupled land surface and subsurface models, featuring a highly dynamic river boundary condition determined by dam-induced river stage variations (Hauer et al., 2017), representative of managed river reaches that are of global significance as a result of dam constructions in the past few decades (Zhou et al., 2016). Third, we assess the effects of spatial resolution and projected hydroclimatic changes on simulated land surface fluxes and exchange of groundwater and river water using the coupled model and datasets from the benchmarking case. In Sect. 2, we describe the component models and our coupling strategy. In Sect. 3, we describe an application of the model to a field site along the Hanford reach of the Columbia River, where the subsurface properties are well characterized and long-term monitoring of river stage, groundwater table, and exchange of groundwater and river water exist. In Sect. 4, we assess the effects of spatial resolution and hydroclimatic conditions to simulated fluxes and state variables. In Sect. 5, conclusions and future work are discussed. CLM4.5 (Oleson et al., 2013) is the land component of the Community Earth System Model version 1 (CESM1) (Hurrell et al., 2013), a fully coupled numerical simulator of the Earth system consisting of atmospheric, ocean, ice, land surface, carbon cycle, and other components. It has been applied successfully to explore interactions among water, energy, carbon, and biogeochemical cycling on local to global scales (Leng et al., 2016b;Xu et al., 2016) and proven to be highly scalable on leading HPC facilities such as the US Department of Energy (USDOE)'s National Energy Research Scientific Computing Center (NERSC). The model includes parameterizations of terrestrial hydrological processes including interception, throughfall, canopy drip, snow accumulation and melt, water transfer between snow layers, infiltration, evaporation, surface runoff, subsurface drainage, and redistribution within the soil column, as well as groundwater discharge and recharge to simulate changes in canopy water, surface water, snow water, soil water, and soil ice, and water in the unconfined aquifer (Oleson et al., 2013). Precipitation is either intercepted by the canopy, falls directly to the snow-soil surface (throughfall), or drips off the vegetation (canopy drip). Water input at the land surface, the sum of liquid precipitation reaching the ground and meltwater from snow, is partitioned into surface runoff, surface water storage, and infiltration into the soil. Two sets of runoff generation parameterizations, including formulations for saturation and infiltration excess runoff and baseflow, are implemented into the model: the TOPMODEL-based runoff generation formulations (Beven and Kirkby, 1979;Niu et al., 2005Niu et al., , 2007 and the variable infiltration capacity (VIC)-based runoff generation formulations (Lei et al., 2014;Liang et al., 1994;Wood et al., 1992). Surface water storage and outflow in and from wetlands and small subgrid-scale water bodies are parameterized as functions of fine-spatial-scale elevation variations called microtopography. Soil water is predicted from a multilayer model based on the one-dimensional Richards equation, with boundary conditions and source-sink terms specified as infiltration, surface and subsurface runoff, gradient diffusion, gravity, canopy transpiration through root extraction, and interactions with groundwater. A groundwater component is added in the form of an unconfined aquifer lying below the soil column following Niu et al. (2007). The model computes surface energy fluxes following the Monin-Obukhov similarity theory using formulations in Zeng et al. (1998), which updates the calculation of boundary resistance to account for understory turbulence, sparse and dense canopies, and surface litter layers (Sakaguchi and Zeng, 2009;Zeng et al., 2005;Zeng and Wang, 2007). Water and energy budgets are conserved at every modeling step.

PFLOTRAN
PFLOTRAN is a massively parallel multiphysics simulator (Hammond et al., 2014) developed and distributed under an open-source GNU LGPL license and is freely available through Bitbucket (https://bitbucket.org/pflotran/pflotran). It solves a system of generally nonlinear partial differential equations (PDEs) describing multiphase, multicomponent, and multiscale reactive flow and transport in porous materials. The PDEs are spatially discretized using a finite-volume technique, and the backward Euler scheme is used for implicit time discretization. It has been widely used for simulating subsurface multiphase flow and reactive biogeochemical transport processes (Chen et al., , 2012Hammond and Lichtner, 2010;Hammond et al., 2011;Kumar et al., 2016;Lichtner and Hammond, 2012;Liu et al., 2016;Pau et al., 2014).
PFLOTRAN is written in object-oriented Fortran 2003/2008 and relies on the PETSc framework (Balay et al., 2015) to provide the underlying parallel data structures and solvers for scalable high-performance computing. PFLOTRAN uses domain decomposition and MPI libraries for parallelization. PFLOTRAN has been run on problems composed of over 3 billion degrees of freedom with up to 262 144 processors, but it is more commonly employed on problems with millions to tens of millions of degrees of freedom utilizing hundreds to thousands of processors. Although PFLOTRAN is designed for massively parallel computation, the same code base can be run on a single processor without recompiling, which may limit problem size based on available memory.
In this study, PFLOTRAN is used to simulate single-phase variably saturated flow and solute transport in the subsurface. Single-phase variably saturated flow is based on the Richards equation with the form with liquid density ρ, porosity φ, and saturation s. The Darcy velocity, q, is given by with liquid pressure p, viscosity µ, acceleration of gravity g, intrinsic permeability k, relative permeability k r , and elevation above a given datum z. Conservative solute transport in the liquid phase is based on the advection-dispersion equation with solute concentration C and hydrodynamic dispersion coefficient D. The discrete system of nonlinear PDEs for flow and transport are solved using the Newton-Raphson method.

Model coupling
In this study, CLM4.5's one-dimensional models for flow in unsaturated (Zeng and Decker, 2009) and saturated (Niu et al., 2007) zones are replaced by PFLOTRAN's RICHRADS mode to simulate unsaturated-saturated flow within the three-dimensional subsurface domain. Although PFLOTRAN is also capable of simulating coupled flow and thermal processes in the subsurface, including explicit representation of liquid-ice phase (Karra et al., 2014) as well as soil nutrient cycles (Hammond and Lichtner, 2010;Zachara et al., 2016;Tang et al., 2016), those processes are not coupled between the two models in this study. A schematic representation of the coupling between CLM4.5 and PFLO-TRAN is shown in Fig. 1. A model coupling interface based on PETSc data structures was developed to couple the two models, and the interface includes some key design features of the CESM coupler (Craig et al., 2012). The model coupling interface allows each model grid to have a different spatial resolution and domain decomposition across multiple processors. While CLM4.5 uses a round-robin decomposition approach, PFLOTRAN employs domain decomposition via PETSc (Fig. 1a). Interpolation of gridded data from one model onto the grids of the other is done through sparse matrix vector multiplication. As a preprocessing step, sparse weight matrices for interpolating data between the two models are saved as mapping files. Analogous to the CESM coupler, the mapping files are saved in a format similar to the mapping files produced by the ESMF_RegridWeightGen (https://www.earthsystemcog.org/ projects/regridweightgen). ESMF regridding tools provide multiple interpolation methods (conservative, bilinear, and nearest neighbor) to generate the sparse weight matrix. In this work, we have used a conservative remapping method to interpolate data between CLM and PFLOTRAN. During model initialization, the model coupling interface first collectively reads all required sparse matrices. Next, the model coupling interface reassembles local sparse matrices after accounting for domain decomposition of each model ( Fig. 1b and c). For a given time step, CLM4.5 first computes infiltration, evaporation, and transpiration within the domain and then sends the data to the model coupling interface. The model coupling interface for each processor receives relevant CLM data vector from all other processors, interpolates data from CLM's grid onto PFLOTRAN's grid via a local sparse matrix vector multiplication, and saves the resulting vector in PFLOTRAN's data structures as prescribed flow conditions (Fig. 1b). PFLOTRAN evolves the subsurface states over the given time step length. The updated soil moisture data simulated by PFLOTRAN are then provided back to the model coupling interface, which interpolates data from PFLOTRAN's grid onto CLM's grid (Fig. 1c). The interpolated data are saved in CLM4.5's data structure and used for simulating land water-and energy-budget terms in the next step. Figure  3 Site description and model configuration

The Hanford site and the 300 Area
The Hanford Reach is a stretch of the lower Columbia River extending approximately 55 km from the Priest Rapids hydroelectric dam to the outskirts of Richland, Washington, USA ( Fig. 3a) (Tiffan et al., 2002). The Columbia River above Priest Rapids Dam drains primarily mountainous regions in Canada, Idaho, Montana, and Washington, over which spatio-temporal distributions of precipitation and snowmelt modulate the timing and magnitude of river flows (Elsner et al., 2010;Hamlet and Lettenmaier, 1999). The Columbia River is highly regulated by dams for power generation, and river stage and discharge along the Hanford Reach displays significant variation on multiple timescales. Strong seasonal variations occur, with the greatest discharge (up to 12 000 m 3 s −1 ) occurring from May through July due to snowmelt, with less discharge (> 1700 m 3 s −1 ) and lower flows occurring in the fall and winter (Hamlet and Lettenmaier, 1999;Waichler et al., 2005). Significant variation in discharge also occurs on a daily or hourly basis due to power generation, with fluctuations in river stage of up to 2 m within a 6-24 h period being common (Tiffan et al., 2002). The Hanford site features an unconfined aquifer developed in Miocene-Pliocene fluvial and lacustrine sediments of the Ringold Formation, overlain by Pleistocene flood gravels of the Hanford Formation (Thorne et al., 2006) that is in hydrologic continuity with the Columbia River. The Hanford Formation gravel and sand, deposited by glacial outburst floods at the end of the Pleistocene (Bjornstad, 2007), has a high average hydraulic conductivity at ∼ 3100 m day −1 (Williams et al., 2008). The fluvial deposits of the Ringold Forma-tion have much lower hydraulic conductivity than those in the Hanford Formation but are still relatively conductive at 36 m day −1 (Williams et al., 2008). Fine-grained lacustrine Ringold silt has a much lower estimated hydraulic conductivity of 1 m day −1 . The hydraulic conductivity of recent alluvium lining the river channel is low relative to the Hanford Formation, which tends to dampen the response of water table elevation in wells near the river when changes occur in river stage (Hammond et al., 2011;Williams et al., 2008). Overall, the Columbia River through the Hanford Reach is a prime example of a hyporheic corridor with an extensive floodplain aquifer. It is consequently an ideal alluvial system for evaluating the capability of the coupled model in simulating stream-aquifer-land interactions.
The region is situated in a cold desert climate with temperatures, precipitation, and winds that are greatly affected by the presence of mountain barriers. The Cascade Range to the west creates a strong rain shadow effect by forming a barrier to moist air moving from the Pacific Ocean, while the Rocky Mountains and ranges to the north protect it from the more severe cold polar air masses and winter storms moving south across Canada. Meteorological data are collected by the Hanford Meteorological Monitoring Network (http://www.hanford.gov/page.cfm/hms), which collects me-teorological data representative of the general climatic conditions for the Hanford site.
A segment of the hyporheic corridor in the Hanford 300 Area was chosen to evaluate the model's capability in simulating river-aquifer-land interactions. Located at the downstream end of the Hanford Reach, the impact of dam operations on river stage is relatively damped, exhibiting a typical variation of ∼ 0.5 m within a day and 2-3 m in a year. The study domain covers an area of 400 m × 400 m along the Columbia River shoreline (Fig. 3b). Aquifer sediments in the 300 Area are coarse-grained and highly permeable Hammond and Lichtner, 2010). Coupled with dynamic river stage variations, the resulting system is characterized by stage-driven intrusion and retreat of river water into the adjacent unconfined aquifer system. During high-stage spring runoff events, river water has been detected in monitoring wells nearly 400 m from the shoreline (Williams et al., 2008). During baseline, low-stage conditions (October-February), the Columbia River is a gaining stream, and the aquifer pore space is occupied by groundwater.
The study domain is instrumented with groundwater monitoring wells (Fig. 3b) and a river gaging station that records water table elevations. A vegetation survey in 2015 was conducted to provide aerial coverages of grassland, shrubland, and riparian trees in the domain (Fig. 3b). A high-resolution topography and bathymetry dataset at 1 m resolution was assembled from multiple surveys by Coleman et al. (2010). The data layers originated from deep-water bathymetric boat surveys, terrestrial light detection and ranging (lidar) surveys, and special hydrographic lidar surveys penetrating through water to collect both topographic and bathymetric elevation data.

Model configuration, numerical experiments, and analyses
To assess the effect of spatial resolution on simulated variables such as latent heat, sensible heat, water table depth, and river water in the domain, we configured CP v1.0 simulations at three horizontal spatial resolutions: 2, 10, and 20 m over the 400 m × 400 m domain. For comparison purposes, we also configured a 2 m resolution CP v1.0 vertical-only simulation (i.e., S v2 m ) in which lateral transfers of flow and solutes in the subsurface are disabled. Due to the lack of observations of water and energy fluxes from the land surface, in this study we treat the 2 m resolution CP v1.0 as the baseline and compare simulation results at other resolutions to it. New hydrologic regimes are projected to emerge over the Pacific Northwest as early as the 2030s due to increases in winter precipitation and earlier snowmelt in response to future warming (Leng et al., 2016a). Therefore, we expect that spring and early summer river discharge along the reach might increase in the future. To evaluate how land surface-subsurface coupling might be modulated hydroclimatic conditions, we designed additional numerical experi- The PFLOTRAN subsurface domain, also terrainfollowing and extending from soil surface (including riverbed) to 32 m below the surface, was discretized using a structured approach with rectangular grids. For the 2, 10, and 20 m resolution simulations, each mesh element was 2 m × 2 m, 10 m × 10 m, and 20 m × 20 m in the horizontal direction and 0.5 m in the vertical direction, giving 2.56 × 10 6 , 99.2 × 10 3 , and 2.48 × 10 3 control volumes in total. The domain contained two materials with contrasting hydraulic conductivities: Hanford and Ringold (Fig. 4). Note that only the soil moisture and soil hydraulic properties within the top 3.8 m are transferred from PFLOTRAN to CLM4.5 to allow simulations of infiltration, evaporation, and transpiration in the next time step, as the CLM4.5 subsurface domain is limited to 3.8 m and cannot currently be easily modified. The hydrogeological properties of the Hanford and Ringold materials (Table 2) were taken from Williams et al. (2008). The unsaturated hydraulic conductivity in PFLTORAN simulations was computed using the Van Genuchten water retention function (van Genuchten, 1980) and the Burdine permeability relationship (Burdine, 1953).
We applied time-varying pressure boundary conditions to PFLOTRAN's subsurface domain at the northern, western, and southern boundaries. The transient boundary conditions were derived using kriging-based interpolations of hourly water table elevation measurements in wells inside and beyond the model domain, following the approach used by Chen et al. (2013). Transient head boundary conditions were applied at the eastern boundary, with water table elevations from the river gaging station and the gradient along the river estimated using water elevations simulated by a onedimensional hydraulic model along the reach, the Modular Aquatic Simulation System 1D (MASS1) (Waichler et al., 2005), with a Nash-Sutcliffe coefficient (Nash and Sutcliffe, 1970) of 0.99 in the simulation period (figure not shown). The river stage simulated by MASS1 was also used to fill river stage measurement gaps caused by instrument failures.
Geosci. Model Dev., 10, 4539-4562, 2017 www.geosci-model-dev.net/10/4539/2017/   A conductance value of 10 −12 m was applied to the eastern shoreline boundary to mimic the damping effect of lowpermeability material on the river bed (Hammond and Lichtner, 2010). A no-flow boundary condition was specified at the bottom of the domain to represent the basalt underlying the Ringold formation. Vegetation types (Fig. 3b) were converted to corresponding CLM4.5 plant functional types (PFTs) and bare soil ( Fig. 5). At each resolution, fractional area coverages of PFTs and bare soil are determined based on the base map and written into the surface dataset as CLM4.5 inputs (Figs. 5, S1, and S2 in the Supplement). The CLM4.5 domain is terrain-following by treating the land surface as the top of the subsurface domain, which is hydrologically active to a depth of 3.8 m. The topography of the domain is retrieved from the 1 m topography and bathymetry dataset (Coleman et al., 2010) based on the North American Vertical Datum of 1988 (NAD88) and resampled to each resolution (Fig. S3).
The simulations were driven by hourly meteorological forcing from the Hanford meteorological stations and hourly river stage from the gaging station over the period of 2009-2015. Precipitation, wind speed, air temperature, and relative humidity were taken from the 300 Area meteorological station (46.578 • N, 119.726 • W) located ∼ 1.5 km from the modeling domain. Other meteorological variables, such as downward shortwave and longwave radiation, were obtained from the Hanford Meteorological station (46.563 • N, 119.599 • W) located in the center of the Hanford site. The first 2 years of simulations (i.e., 2009 and 2010) were discarded as the spin-up period, so that 2011-2015 is treated as the simulation period in the analyses.
Among the hydroclimatic forcing variables (e.g., river stage, surface air temperature, incoming shortwave radiation, and total precipitation), river stage displayed the greatest interannual variability (Fig. 6). During the study period, high river stages occurred in early summer of 2011 and 2012 due to the melt of above-average winter snow packs in the upstream drainage basin, typical flow conditions occurred in 2013 and 2014, while 2015 was a year with low upstream snow accumulation. Meanwhile, the meteorological variables, especially temperature and shortwave radiation, do not show much interannual variability or changes in trends, while precipitation in late spring (i.e., May) of 2012 is higher than that in the other years, coincident with the high river stage in 2012. In the "elevated" experiments (i.e., S E2 m , S E10 m , and S E20 m ), the observed river stage (meters based on NAD88) was increased by 5 m at each hourly time step to mimic a perturbed hydroclimatic condition in response to future warming.
To evaluate effects of river water and groundwater exchanges on land surface energy partitioning, we separated the study domain for the 2 m simulations with lateral water B Figure 5. Plant function types at 2 m resolution as inputs for CLM4.5. exchange (i.e., S 2 m and S E2 m ) into two subdomains based on 2 m topography (shown in Fig. S3a): (a) the inland domain where the surface elevation is higher than 110 m; and (b) the riparian zone where the surface elevation is less than or equal to 110 m. In addition to the latent heat flux, the evaporative fraction, defined as the ratio of the latent heat flux to the sum of latent and sensible heat fluxes, was calculated over the subdomains for both observed and elevated conditions at a daily time step for all days with significant energy inputs (i.e., when net radiation is greater than 50 W m 2 ). The evaporative fraction is an indicator of the type of surface, as summarized in literature (Lewis, 1995): it is typically less than 1 over surfaces with abundant water supplies; ranges between 0.75 and 0.9, between 0.5 and 0.7, and between 0.15 and 0.3 for tropical rainforests, temperate forests and grasslands, and semiarid landscapes, respectively; and approaches 0 over deserts.
To better quantify the spatio-temporal dynamics of stream-aquifer interactions, a conservative tracer with a mole fraction of one was applied at the river boundary to track the flux of river water and its total mass in the subsurface domain. While a constant concentration was maintained at the river (i.e., eastern) boundary, the tracer was allowed to be transported out of the northern, western, and southern boundaries. Water infiltrating at the upper boundary based on CLM4.5 simulations was set to be tracer free, while a zero-flux tracer boundary condition was applied at the lower boundary. The initial flow condition was a hydrostatic pressure distribution based on the water table, as interpolated from the same set of wells that were used to create the transient lateral flow boundary conditions at the northern, western, and southern boundaries. The initial conservative tracer concentration was set to be zero for all mesh elements in the domain. The simulations were started on 1 January 2009 and the first 2 years were discarded as the spin-up period in the analysis. The mass of tracers in the domain and the fluxes of tracers across the boundary allow us to quantitatively understand how river water is retained and transported in the subsurface domain.
In addition to the CP v1.0 simulation, a standalone CLM4.5 simulation was also configured and performed (i.e., CLM 2 m in Table 1). CLM 2 m shared the same subsurface properties and initial conditions as the CLM4.5 setup in S 2 m and S v2 m where CP v1.0 were used. However, we  Figure 6. Hydrometeorological drivers in the study period: (a) monthly mean river stage, (b) monthly total precipitation, (c) monthly mean surface air temperature, and (d) monthly mean incoming shortwave radiation.
note that CLM 2 m simulations are not directly comparable to other simulations listed in Table 1 for the following reasons: (1) the CLM4.5 simulates subsurface hydrologic processes only down to 3.8 m below the surface, while the CP v1.0 subsurface domain extends to ∼ 30 m below the surface; (2) as discussed in Sect. 2.1, CLM4.5 uses TOPMODEL-based parameterizations to simulate surface and subsurface runoffs, as well as mean groundwater table depth, using formulations derived from catchment hydrology that are only applicable at coarser resolutions; and (3) the key hydrologic processes (i.e., the exchange of river water and groundwater at the eastern boundary and lateral transfer of water at all other boundaries) that affect the hydrologic budget of the system are missing from CLM4.5. Therefore, the simulation was performed to characterize how physical parameterizations from one scale (i.e., catchment scale) affect simulations on another scale (i.e., field scale) where those physical parameterizations may not apply.

Model evaluation
For the three-dimensional numerical experiments driven by the observed river stage time series (i.e., S 2 m , S 10 m , S 20 m ), CP v1.0 simulated soil water pressure was converted to water table depth and compared against observed values at selected wells that were distributed throughout the domain and of variable distances from the river (Figs. 7, S4 and Table 3). The model performed very well in simulating the temporal dynamics of the water table at all resolutions. The rootmean-square errors were 0.028, 0.028, and 0.023 m at 2, 10, and 20 m resolutions, respectively. The corresponding Nash-Sutcliffe coefficients were 0.998, 0.998, and 0.999. It was surprising that the performance metrics at 20 m resolution matches the observations better than those at finer resolutions, but the differences were marginal given the close match between the model simulation results and observations. River stage was clearly the dominant driving factor for water table fluctuations at the inland wells. In addition, errors in water and tracer budget conservations and surface energy conservation for each time step in S 2 m are shown in Fig. S5a, b, and c, respectively. The errors are sufficiently small when compared to the magnitudes of the related fluxes to ensure faithful simulations in CP v1.0. These results indicated that the coupled model was capable of simulating dynamic stream-aquifer interactions in the near-shore groundwater aquifer that experiences pressure changes induced by river stage variations on subdaily timescales.

Effect of stream-aquifer interactions on land surface energy partitioning
Next we evaluated the role of water table fluctuations on land surface variables, including latent heat (LH) and sensible heat (SH) fluxes. The site is characterized by an approximate 10 m vadose zone, and surface fluxes and groundwater dynamics are typically decoupled , especially over the inland portion of the domain covered by shallow-rooted PFTs and with higher surface elevations. However, river discharge and water table elevation displayed large seasonal and interannual variability in the study period. Therefore, we selected the month of June in each year to assess potential land surface-groundwater coupling because it is the month of peak river stage, while energy input is high and relatively constant across the years (Fig. 8a).
In June 2011 and 2012, high river stages push the groundwater table to ∼ 108 m (or ∼ 6 m below the land surface). Groundwater at that elevation can affect land surface water and energy exchanges with the atmosphere. The shrubs, including the patch of Basin big sagebrush and the mixture of  rabbitbrush and bunchgrass on the slope close to the river, are able to tap into the elevated water table with their deeper roots. In the inland portion of the domain, capillary supply was most evident in high-water years (i.e., 2011 and 2012), remains influential in normal years (i.e., 2013 and 2014), and is essentially disabled in low-water years (i.e., 2015). The lateral discharge of shallow groundwater to the river led to a band of negative difference in LH between S 2 m and S v2 m at the river boundary when the stage was low due to a decrease in rooting zone soil moisture for evapotranspiration by the riparian trees (Fig. 8b). This pattern was most evident in June 2015. Such a mechanism decreases in high-water and normal years because of more frequent inundation of the river bank and groundwater gradient reversal. Driven by elevated river stages, land surface energy partitioning in S E2 m (Figs. 9 and 10) was significantly shifted from that in S 2 m (Fig. 8a) through two mechanisms: (1) expanding the periodically inundated fraction of the riparian zone (i.e., surface elevation ≤ 110 m) and (2) enhancing moisture availability in the vadose zone in the inland domain (i.e., surface elevation > 110 m) through capillary rise. Both mechanisms led to general increases in simulated vadosezone moisture availability and therefore higher latent heat fluxes compared to the simulations driven by the observed condition. For the inland domain, the evaporative fraction clearly displayed an increasing trend as the groundwater ta- ble level becomes shallower, a trend which is consistent between the simulations (Fig. 10c). The daily evaporative fractions for the inland domain stayed well below 0.2 when the water table levels are less than 112 m, suggesting decoupled surface-subsurface conditions in a typical semiarid environment. When water table levels increased to be above 112 m, the evaporative fraction increases to ∼ 0.2, indicating that the surface and subsurface processes become more strongly coupled because of improved water availability for evapotranspiration, especially in the elevated simulation (i.e., S E2 m ). The evaporative fraction in the riparian zone remained close to 1.0, suggesting strong influences of the river and the role of deeper rooted plant types (e.g., riparian trees and shrubs) in modulating the energy partitioning (Fig. 10d) of riparian zones in the semiarid to arid environments. To confirm the above findings, the liquid saturation (unitless) and mass of river water (mol) in the domain from S 2 m and S E2 m on 30 June each year are plotted along a transect perpendicular to the river (y = 200 m) in Figs. 11 and S6, and across a x-y plane at an elevation of 107 m in Figs. S7 and S8, respectively. Driven by the pressure introduced by elevated river stages, river water not only intruded further toward or even across the western boundary in high-water years, but also led to a shallower water table and increased liquid saturation in the vadose zone due to capillary rise across the domain. In fact, liquid saturation in the shallow vadose zone could increase from 0.1-0.2 in S 2 m to 0.3-0.4 in S E2 m on these days because of river-water intrusion. The river-water tracer could show up in the near-surface vadose zone at a distance of ∼ 400 m from the river (Fig. S6). Interestingly, by comparing the spatial distributions of riverwater tracer in the low-water year (i.e., 2015) between the "observed" and "elevated" scenarios, the presence of river water in the domain was much less in the elevated scenario in terms of its spatial coverage (Figs. 11 and S6). This pattern suggests that after a number of years of enhanced riverwater intrusion into the domain, the hydraulic gradient between groundwater and river water could be reversed, so that groundwater discharging might be expected more frequently in low-water years in a prolonged elevated scenario.
The responses of LH and evaporative fraction (Figs. 9 and 10) indicated that a tight coupling among stream, aquifer, and land surface processes occurred in the elevated scenario, which could become realistic in 1 to 2 decades for the study site, or for other sites along the Hanford reach characterized by lower elevations under the current condition. As discussed in Sect. 2.1 and 3.2., the hydrologic parameterizations in the default CLM4.5 model are based on conceptual and physical understandings from watershed hydrology that do not apply on the scale of our study site, where the exchange of river water and groundwater dominates the hydrologic budget of the system. Nevertheless, a comparison between CLM4.5 and CP v1.0 helps characterize how scale inconsistencies in physical representations affect the simulations. Figure 12 shows comparisons of key components in the hydrologic budget between the two models. The simulated mean water table elevation of the domain from CLM4.5 ranges between 74 and 80 m (i.e., 35-40 m below the surface), while the observed water table elevation ranges between 104 and 108 m (i.e., 5-10 m below the surface), and was accurately reproduced by S 2 m (Fig. 12a). By using physics derived for the larger scale, CLM4.5 could not capture subsurface river-water and groundwater exchanges, and consequently cannot accurately simulate groundwater table dynamics for our study domain.
At this semiarid field site, the groundwater and riverwater exchanges represented in S 2 m recharges the unconfined aquifer, and hence maintains sufficient soil water availability in the top 3.8 m of the soil column, while the lack of groundwater and river-water interactions in CLM 2 m leads to overall declining soil water content with seasonal variability as a result of percolation of winter rainwater (Fig. 12b). The difference in soil moisture availability propagates to evapotranspiration (ET) and its components (Fig. 12c-f). Simulated summer ET in CLM 2 m shows a high-frequency signal in response to rainfall pulses through ground evaporation. Transpiration simulated by CLM 2 m is determined by soil water availability in the soil column. In the spring and early summer of 2011 and 2013, transpiration from CLM 2 m is close to that from S 2 m given sufficient soil water. For other periods, CLM 2 m simulates significant lower transpiration rates compared to S 2 m .
Simulated latent heat fluxes in June for the period of 2011-2015 from CLM 2 m and their differences from those in S 2 m are also illustrated in Fig. S9a and b. Evidently, the hydrologic gradient from river to inland is missing as CLM4.5 lacks the capability of capturing the river stage dynamics at such a resolution (in Fig. S9a). Instead, although initiated from the same initial condition as S 2 m on 1 January 2009 as discussed in the spin-up procedure in Sect. 3.2, soil moisture at the grid cells inundated or periodically inundated by the river is soon depleted through ET, surface runoff, or baseflow. However, latent heat from the inland domain is generally higher in CLM 2 m than in S 2 m due to ground evaporation in response to rainfall pulses. In short, CLM4.5 fails to capture the dynamics of groundwater and river-water exchanges. These biases propagate to simulated water and energy fluxes, which could have large impacts on boundary layer evolution, convection, and cloud formation in coupled land-atmosphere studies.

Effect of spatial resolution
To apply the model to large-scale simulations or over a long time period, it is important to assess how the model performs at coarser resolution, as the 2 m simulations are computationally expensive. Here, we use the 2 m simulations (i.e., S 2 m and S E2 m ) as benchmarks for this assessment. That is, S 2 m and S E2 m simulated variables are treated as the "truth" for "observed" and "elevated" river stage scenarios, and outputs from other simulations are compared to them to verify their performance. In the previous section, we showed that simu- lated water table levels from the model were virtually identical to observations. In this section, we further quantify biases of other variables of interest from the high-fidelity 2 m simulations.
The domain-averaged daily surface energy fluxes from S 2 m show clear seasonal patterns, which are consistent in terms of their magnitudes and timing, reflecting mean climate conditions at the site (Fig. S10). Driven by elevated river stages, latent heat from S E2 m is consistently higher than that from S 2 m . The mean latent heat and sensible heat fluxes simulated by S 2 m were 14.1 and 38.7 W m −2 over this pe-riod, compared to by 18.50 and 35.75 W m −2 in S E2 m . Figure 13 shows deviations of simulated LH and SH in the 20 and 10 m simulations from the corresponding 2 m simulations. The deviations of both LH and SH were small across all the simulations driven by the observed river stage when surface and subsurface were decoupled. In the elevated simulations (i.e., S E10 m and S E20 m ) when surface and subsurface processes are more tightly coupled, errors in surface fluxes became significant in the coarse-resolution simulations when compared to S E2 m . For example, the relative errors in LH (Table 4)  tively, compared to S 2 m , but grew as large as 33.84 and 33.19 % for S E20 m and S E10 m , respectively, when compared to S E2 m . The 10 m simulations outperformed the 20 m simulations under both scenarios but the magnitudes of errors were comparable. However, notably the vertical-only simulation (S v2 m ) has a small error of 5.67 % in LH compared to S 2 m , indicating that lateral flow is less important when the water table is deep.
To better understand how water in the river and the aquifer was connected, we also quantified the biases of subsurface state variables and fluxes including total water mass and tracer amount, as well as exchange rates of water and tracer at four boundaries of the subsurface domain using a similar approach (Figs. S11 and 14). Compared to the magnitude of total water mass in the domain (averaged 919.45 × 10 6 and 1020.19 × 10 6 kg in S 2 m and S E2 m ), errors introduced by coarsening the resolution were very small under the observed river stage condition (0.04 % for S 20 m and 0.03 % for S 10 m ) and grew to 9.85 % for S E20 m and 9.87 % for S E10 m in terms of total water mass in the domain (Table 5). However, for total tracer in the domain (averaged 142.07 × 10 6 and 172.46 × 10 6 mol in S 2 m and S E2 m ) as a result of transport of river water in lateral and normal directions to the river, resolution clearly makes a difference under both observed condition and elevated scenarios (relative errors of 5.44 % for S 10 m , 10.40 % for S 20 m , and 22.0 % for both S E10 m and S E20 m ). The magnitude of computed mass exchange rates at the four boundaries (Fig. S11) indicates that a coarse resolution promotes larger river-water fluxes and groundwater exchanges, especially during the period of spring river stage increase under the elevated scenario. This forcing contributes to a significant bias in total tracer amount by the end of the simulation. The exchange rates at the other three boundaries follow the same pattern but with smaller magnitudes, especially for the west boundary that requires a significant gradient high enough to push river water further inland. Figure 14. Deviations of total water mass, tracer, and exchange rates of water and tracer at four boundaries from those simulated by S 2 m (for S 10 m and S 20 m ) and by S E2 m (for S E10 m and S E20 m ).
The results of simulations at three different resolutions indicated the following: (1) the partitioning of the land surface energy budget is mainly controlled by near-surface moisture -spatial resolution did not seem to be a significant factor in the computation of surface energy fluxes when the water table was deep at the semiarid site; (2) if the surface and subsurface are tightly coupled as in the elevated river stage simulations, resolution becomes an important factor to consider for credible simulations of the surface fluxes, as the land surface, subsurface, and riverine processes are expected to be more connected and coupled; (3) regardless of whether a tight coupling between the surface and subsurface occurs, if mass exchange rates and associated biogeochemical reactions in the aquifer are of interest, a higher resolution is desired close to the river shoreline to minimize terrain errors.

Discussion and future work
A coupled three-dimensional surface and subsurface land model was developed and applied to a site along the Columbia River to simulate interactions among river water, groundwater, and land surface processes. The model features the coupling of the open-source and state-of-the-art models portable on HPCs, the multiphysics reactive transport model Geosci. Model Dev., 10, 4539-4562, 2017 www.geosci-model-dev.net/10/4539/2017/ By applying the coupled model to a field site along the Columbia River shoreline driven by highly dynamic river boundary conditions resulting from upstream dam operations, we demonstrated that the model can be used to advance mechanistic understanding of stream-aquifer-land interactions surrounding near-shore alluvial aquifers that experience pressure changes induced by river stage variations along managed river reaches, which are of global significance as a result of over 30 000 dams constructed worldwide during the past half-century. The land surface, subsurface, and riverine processes along such managed river corridors are expected to be more strongly coupled under projected hydroclimatic regimes as a result of increases in winter precipitation and early snowmelt. The dataset presented in this study can serve as a good benchmarking case for testing other coupled models for their applications to such systems. More data need to be collected to facilitate the application and validation of the model to a larger domain for understanding the contribution of near-shore hydrologic exchange to water retention, biogeochemical cycling, and ecosystem functions along the river corridors.
By comparing simulations from the coupled model (CPv1.0) to that from CLM4.5, we demonstrated that the catchment-scale physics imbedded in CLM4.5 does not apply on the field scale. By misrepresenting, or not including, key hydrologic processes on the scale of interest, CLM4.5 fails to capture groundwater table dynamics, which could propagate to water and energy budgets and have profound impacts on boundary layer, convection, and cloud formation in coupled land-atmosphere studies. Our finding is consistent with results from other recent studies in which integrated surface and subsurface models were compared to standalone land surface models (Fang et al., 2017;Niu et al., 2017).
By benchmarking the coarser-resolution simulations at 20 and 10 m against the 2 m simulations, we find that resolution is not a significant factor for surface flux simulations when the water table is deep. However, resolution becomes important when the surface and subsurface processes are tightly coupled, and for accurately estimating the rate of mass exchange at the riverine boundaries, which can affect the calculation of biogeochemical processes involved in carbon and nitrogen cycles.
Our numerical experiments suggested that riverine, land surface, and subsurface processes could become more tightly coupled through two mechanisms in the near-shore environments: (1) expanding the periodically inundated fraction of the riparian zone and (2) enhancing moisture availability in the vadose zone in the inland domain through capillary rise. Both mechanisms can lead to increases in vadose-zone moisture availability and higher evapotranspiration rates. The latter is critical for understanding ecosystem functioning, biogeochemical cycling, and land-atmosphere interactions along river corridors in arid and semiarid regions that are expected to experience new hydroclimatic regimes in a changing climate. However, these systems have been poorly accounted for in current-generation Earth system models and therefore require more attention in future studies.
We acknowledge that there are a number of limitations of this study that need to be addressed in future studies.
1. In order to understand the stream-aquifer-land interactions with a focus on groundwater and river-water interactions along a river corridor situated in a semiarid climate, the river boundary conditions were prescribed using observations with gaps filled by a onedimensional hydrodynamics model. Future versions of the CP model need to incorporate two-way interactions between stream and aquifer by developing a surface flow component and testing the new implementation against standard benchmark cases Maxwell et al., 2014).
2. We note that CLM estimates the surface heat and moisture fluxes using the Monin-Obukhov similarity theory (Sect. 2.1), which is only valid when the surface layer depth z z 0 , where z 0 is the aerodynamic roughness length. As reviewed by Basu and Lacser (2017), it is highly recommended that z > 50 z 0 , which should be proportional to the horizontal grid spacing to guaran-tee the validity of the Monin-Obukhov similarity theory (Arnqvist and Bergström, 2015). In our simulations, the majority of the Hanford 300 Area domain is covered by bare soil (z 0 = 0.01 m), grass (z 0 = 0.013 m), shrubs (z 0 = 0.026-0.043 m), and riparian trees (varies across the seasons, z 0 = 0.008 m when LAI = 2 in the summer and z 0 = 1.4 when LAI = 0 in the winter). Therefore, a 2 m resolution is sufficiently coarse under most conditions except for the grid cells covered by riparian trees in the winter. Nevertheless, the wintertime latent heat and sensible heat fluxes are nearly zero due to extremely low energy inputs. Therefore, the 2 m simulations supported by the dense groundwater monitoring network at the site provide a valid benchmark for the coarser-resolution simulations. For future applications of the coupled model, caution should be taken to evaluate the site condition for the validity of model parameterizations.
3. We used the simulated surface energy fluxes from S 2 m to verify coarser-resolution simulations. The simulated surface energy flux needs to be validated against eddy covariance tower observations, which are not available yet at the site. Nevertheless, we have made initial efforts to install eddy covariance systems at the site (see description in Sect. 3.1 of Gao et al., 2017), but the processing of the flux data is still preliminary. We will report flux observations and validations of the surface energy budget simulations in future studies.
4. Even when observed fluxes are available for validation, the model structural problems associated with ET parameterizations in CLM4.5 need to be addressed for reasonable simulations of the ET components, especially for the study site. That is, it has been well documented that ET simulated by CLM4.5 and CLM4 could be enhanced when vegetation is removed. This ET enhancement over bare soil has been documented as a counterintuitive bias for most unsaturated soils in CLM4 and CLM4.5 simulations (Lawrence et al., 2012;Tang and Riley, 2013a). Tang and Riley (2013a) explored a few potential causes for this likely bias (e.g., soil resistance, litter layer resistance, and numerical time step). They found the implementation of a physically based soil resistance lowered the bias slightly, but concluded that the bias remained (Tang and Riley, 2013b). Meanwhile, in studying ET over semiarid regions, Swenson and Lawrence (2014) proposed another soil resistance formulation to fix this excessive soil evaporation problem within CLM4.5. While their modification improved the simulated terrestrial water storage anomaly and ET when compared to GRACE data and FLUXNET-MTE data, respectively, the empirical nature of the soil resistance proposed could have underestimated the soil resistance variability when compared to other estimates (Tang and Riley, 2013b).
Code and data availability. CLM4.5 (Oleson et al., 2013) is an open-source software released as part of the Community Earth System Model (CESM) version 1.2 (http://www.cesm.ucar.edu/models/cesm1.2). The version of CLM4.5 used in CP v1.0 is a branch from the CLM developer's repository. Its functionality is scientifically consistent with descriptions in Oleson et al. (2013) with source codes refactored for a modular code design. Additional minor code modifications were added by the authors to support coupling with PFLOTRAN (Lichtner et al., 2015). Permission from the CESM Land Model Working Group has been obtained to release this CLM4. The Supplement related to this article is available online at https://doi.org/10.5194/gmd-10-4539-2017supplement.