Interactive comment on “Coupling a groundwater model with a land surface model to improve water and energy cycle simulation” by W. Tian et al

RC2: Improve English description of the manuscript. For example, Line 26-27, P 10946. ‘Although a coupled model is developed in this paper, there are many energy and water processes that are not considered in our model, such as surface water processes’ should be changed into ‘Although a coupled model was developed in this study, some the energy and water processes, such as surface water processes, are not considered in the model’.


Introduction
Water movement and energy transfer in the soil-plantatmosphere continuum are the main processes on the land surface, and the two processes strongly interact.Land surface models (LSMs) are often used to model these physical processes.However, almost all LSMs are 1-D vertical models because the initial aim of these models was to provide a land surface condition for atmospheric models, such as general circulation models (GCMs) and regional climate models.Therefore, these models generally cannot simulate subsurface lateral water movement.However, many studies have indicated that lateral water movement can significantly affect land surface water and energy processes (Holt et al., 2006;Maxwell et al., 2007;Kollet and Maxwell, 2008;Maxwell et al., 2011;Soylu et al., 2011).
Many groundwater models (GWMs), such as the MODFLOW-HYDRUS (Twarakavi et al., 2008) and ParFlow (Kollet and Maxwell, 2006) models, are based on hydrodynamic mechanisms.These models describe the water balance processes and physical mechanism of the subsurface water movement in both saturated and unsaturated zones, but they usually do not explicitly describe the physical mechanism of evapotranspiration (ET) processes.ET is an integration process that includes water, energy and biological processes.The latter two processes are usually not included in GWMs.
In MODFLOW (McDonald and Harbaugh, 1988), for example, ET is calculated using a linear empirical function of the groundwater table (GWT).Although this approach makes the groundwater modeling system self-closing, it oversimplifies the ET simulation, leading to simulation error.
Based on this comparison of the two types of models, coupling LSMs and GWMs could eliminate their separate disadvantages and make the simulation of energy and water cycles in the Earth's surface more accurate on mechanism.

W. Tian et al.: Coupling a groundwater model with a land surface model to improve water simulation
In recent years, many studies have been conducted on the coupling of GWMs with LSMs.Gutowski et al. (2002) developed the Coupled Land-Atmosphere Simulation Program (CLASP) to study the coupled aquifer, land surface and atmospheric hydrological cycle.This model considered the groundwater as a reservoir.York et al. (2002) improved this model in CLASP II by integrating the soil vegetation zone into the USGS groundwater flow model, MODFLOW.Liang et al. (2003) developed a 1-D dynamic groundwater parameterization and implemented it in a three-layer variable infiltration capacity (VIC-3L) model to simulate surface and groundwater interaction dynamics for LSMs.Gedney and Cox (2003) coupled the Hadley Centre Atmospheric Climate Model (HadAM3) with the Met Office Surface Exchange Scheme (MOSES), in which the local water table depth was used to estimate the saturated fraction of the soil layers and to improve the runoff and global wetland area simulation.Tian et al. (2006) implemented a subsurface runoff scheme with a variable water table in Community Land Model 2 (CLM2).Yeh and Eltahir (2005) incorporated a lumped unconfined aquifer model into the Land Surface Transfer Scheme (LSX) to address the deficiency of the simplified representation of subsurface hydrological processes in current land surface parameterization schemes.In the model, groundwater was modeled as a nonlinear reservoir that received the recharge from the overlying soils and discharged the runoff into streams.Niu et al. (2007) developed a simple groundwater model (SIMGM) by representing groundwater recharge and discharge processes, which were added as a single integration element below the soil of the CLM3 land surface model.Fan et al. (2007) coupled a simple 2-D groundwater flow model with the VIC model to estimate the equilibrium water table, which is the result of long-term climatic and geologic effects.Maxwell and Miller (2005) coupled the Common Land Model with a variably saturated GWM (ParFlow) to create a single-column model to improve the groundwater representation in land surface schemes.Kollet and Maxwell (2008) then improved this model and developed an integrated, distributed watershed model based on the coupling of ParFlow and the Common Land Model (PF.CLM).Yuan et al. (2008) coupled a groundwater model with the Biosphere-Atmosphere Transfer Scheme (BATS) and the RegCM3 regional climate model to investigate the local and remote effects of water table dynamics on the regional climate.Maxwell et al. (2011) coupled ParFlow with the Noah LSM, which is a land surface module of the Weather Research and Forecasting (WRF) model, to study the impact of groundwater on weather.These studies have shown that groundwater flow and land surface processes are closely related and that coupled models can simulate complex processes more realistically than uncoupled models.
The GWMs used in these coupled models can be classified into two categories: empirical, lumped GWMs (Liang et al., 2003;Gedney and Cox, 2003;Yeh and Eltahir, 2005;Niu et al., 2007) and physically based distributed GWMs (Gutowski et al., 2002;York et al., 2002;Fan et al., 2007;Kollet and Maxwell, 2008;Maxwell et al., 2011).Of the two categories, the lumped GWMs require less data and are more efficient, but they can only provide a groundwater impact factor for the LSMs and cannot truly describe groundwater flow.Although distributed dynamic GWMs increase the complexity of the simulation, they have the advantage that the groundwater flow mechanism is included in the coupled model.Consequently, effects such as the interaction between the groundwater dynamics and the land-surface energy and water processes can be described in the coupled model.Additionally, the water balance can be maintained when a groundwater model is used instead of the groundwater level in the simulation.
In this paper, a coupled model of groundwater flow with a simple biosphere (GWSiB) is developed by coupling a 3-D variably saturated dynamic GWM (AquiferFlow) (Wang, 2007;Wang et al., 2010) with a typical land surface model, the Simple Biosphere Model version 2 (SiB2) (Sellers et al., 1996b).In this coupled-model system, the coupling mode of physically based distributed GWMs is used in the GWSiB model, and the land surface parameterization scheme and the GWM used are different from other studies.Furthermore, flexible temporal and spatial coupling schemes are used to increase the applicability of the model.
The paper is arranged as follows.In Sect.2, the models to be coupled are briefly introduced, and the coupling method used in the GWSiB model is described in detail.In Sect.3, a sensitivity analysis of the model is performed, and the model is validated using data measured from two sites.The model is then tested in a regional simulation of the middle reach of the Heihe River basin in northwestern China in Sect. 4. The discussion and conclusions are presented in Sect. 5.

Model development
The critical aspect of the coupling of the GWM and the LSM is the soil moisture movement in the vadose zone, which is simulated in both models.If the vadose zone water movement in the two models can be linked and integrated, a water cycle process can be completely simulated.In our study, a GWM called AquiferFlow that can simulate water movement in saturated and unsaturated zones (Wang, 2007;Wang et al., 2010) and an LSM model, SiB2, which is used to simulate the water and energy movement above the ground surface and the soil moisture movement in up to three layers of the subsurface, were chosen as the two models to be coupled.The two models are coupled by replacing the three-layer soil moisture simulation in the SiB2 code by the unsaturated zone water movement simulation of AquiferFlow.The principle of these two models and their coupling scheme are presented in detail in the following.AquiferFlow (Wang, 2007;Wang et al., 2010) is a typical numerical GWM in which rectangular grid cells are used to divide the simulation domain and the finite difference method (FDM) is used to solve groundwater movement equations in saturated and unsaturated zones.The main feature of AquiferFlow is that the conceptual model is similar to traditional saturated GWMs, but the Richards equation is fully incorporated to handle unsaturated flow, as follows: where x, y and z are coordinates (m) in which z is oriented positively upward, H is the hydraulic head (m) H = Z + ψ, where ψ is the water suction potential of the unsaturated soil (m), and Z is the relative hight hydraulic head according to the the z-coordinate (m); K xx , K yy and K zz are the saturated hydraulic conductivity (m s −1 ) in the x, y and z directions, respectively; K r is the relative hydraulic conductivity, defined as the unsaturated conductivity divided by the saturated conductivity; S s is the extended specific storage and can be used in the saturated and unsaturated zones (m −1 ); ε is the source and sink term (s −1 ); and t is the time (s).Compared with the traditional governing equation for saturated groundwater, Eq. ( 1) introduces new parameters: the relative hydraulic conductivity (K r ) and the extension of the specific storage (S s ) for the calculation of the water movement in an unsaturated zone.The relative hydraulic conductivity is a function of the pressure head.In AquiferFlow, Gardner and Fireman's method (1958) is used to relate K r to the soil moisture potential ψ: where C k is the attenuation coefficient of permeability (m −1 ), and ψ s is the saturation moisture potential (m).The specific storage is defined using different forms in the saturated (ψ ≥ ψ s ) and unsaturated (ψ < ψ s ) conditions.The specific storage depends on the compressibility of the porous media and the water in the saturated zone and is a function of the soil volumetric water content (θ) (m 3 m −3 ) and the soil moisture potential (ψ) in the unsaturated zone.
where C s is the specific moisture capacity (m −1 ), ρ w is the water density (kg m −3 ), α is the coefficient of soil compressibility (m s 2 kg −1 ), β is the coefficient of groundwater compressibility (m s 2 kg −1 ), φ is the porosity, and g is the gravitational acceleration (m s −2 ).
The relationship between the soil moisture potential ψ and the moisture content θ in the aquifer media is described by the suction curve, θ (ψ).In AquiferFlow, the suction curve can be created using the commonly used van Genuchten (1980) equation or a simple exponential equation such as where S e is the effective saturation, θ r is the residual saturation, and C w is the attenuation coefficient of the soil moisture (m −1 ), which is an empirical parameter.
The equations of AquiferFlow are based on the principles of groundwater dynamics; consequently, AquiferFlow can describe water movement not only in the saturated zone but also in the unsaturated zone (Wang et al., 2010).However, the ET simulation in AquiferFlow, as in most existing GWMs, is treated as a sink term for the groundwater system, and the ET calculation depends only on the soil moisture of the top soil layer and the potential ET of the location.The ET simulation in AquiferFlow is empirical and is unable to represent a complex physical process.In summary, the governing equations of GWMs are derived from water conservation and do not include the simulation of energy and biological processes.The water phase-change process of evaporation and the vegetation root-uptake process of transpiration are usually simplified in the parameterization scheme of GWMs.Therefore, a model that can simulate the energy cycle and the biological processes, such as an LSM, need to be coupled with the GWM to overcome this weakness.

SiB2
The simple biosphere model (SiB) (Sellers et al., 1986) is a typical land surface model that can be used to calculate the transfer of energy, mass, and momentum between the atmosphere and the vegetated surface of the Earth.Sellers et al. (1996b) produced a new version of the model, SiB2, by improving the hydrological sub-model, increasing the canopy photosynthesis/conductance model and introducing snowmelt process simulation into the SiB model.Since the release of SiB, the model has been verified and applied in many case studies (Sellers et al., 1989;Colello et al., 1998;Sen et al., 2000;Baker et al., 2003;Li and Koike, 2003;Gao et al., 2004;Vidale and Stockli, 2005).The study results show that the model can provide an adequate description of the energy and water processes above the ground surface.
The model structure of SiB2 is divided into five layers in the vertical direction.One vegetation layer and one ground layer lie above the ground surface.The three layers below the ground surface represent the surface soil layer, the root layer, and the deep soil layer; the characteristics of these layers are determined according to the soil properties of a study area.Ground vegetation is classified into nine types to represent W. Tian et al.: Coupling a groundwater model with a land surface model to improve water simulation various global vegetation conditions.The main inputs of the model include meteorological data, soil data, and the morphological, physiological, and biophysical parameters of the vegetation.To explain the model coupling, we only provide a detailed description of the soil moisture movement and evaporation processes; descriptions of the other processes can be found in the relevant literature (Sellers et al., 1986(Sellers et al., , 1989(Sellers et al., , 1996a,b),b).
In the SiB2 model, precipitation reaches the ground surface after the canopy interception, and some water infiltrates to the subsurface, limited by the local soil infiltration capacity.If the residual precipitation still exceeds the groundwater storage capacity, runoff is generated.This process can be expressed as where R 0 is the runoff (m s −1 ), P g is the precipitation reaching the ground surface (m s −1 ), Q 1 is the infiltrated water from the ground surface to the first soil layer (m s −1 ), and E gi is the evaporation from the water intercepted by the ground surface (m s −1 ).
As the water infiltrates into the soil, the water balance in the three soil layers is defined as where W i is the wetness in the i-th layer, W i = θ i /θ s , θ i is the volumetric soil moisture content in the i-th layer, and θ s is the volumetric soil moisture content for the saturation condition; E gs is the evaporation from the surface soil layer (m s −1 ); E ct is the vegetation canopy transpiration (m s −1 ); Q i,i+1 is the water exchange between the i and i + 1 layers (m s −1 ); Q 3 is the base flow that drains out from the bottom of the soil (m s −1 ); and D i is the thickness of the i-th layer (m).
The moisture movement between the soil layers is described by the Richards equation, and the unsaturated zone hydraulic conductivity K (m s −1 ) and the soil moisture potential ψ (m) are calculated using the scheme of Clapp and Hornberger (1978), which is a function of the soil wetness (W ).
where K s is the saturated hydraulic conductivity (m s −1 ), ψ s is the saturation moisture potential (m), and B is an empirical parameter.K s , ψ s and B are dependent on the soil characteristics and can be obtained from the empirical equations, which are defined as a function of the soil texture (Cosby et al., 1984;Yang et al., 2005).
The ET processes in the SiB2 model consist of four parts: vegetation canopy transpiration (E ct ), evaporation from the interception of the canopy (E ci ), evaporation from the interception of the ground (E gi ) and evaporation from the surface soil layer (E gs ).The formulas used in the calculation of the ET are similar with the electrical analog form; the ET, defined as the latent heat flux, is equal to the vapor pressure difference divided by the resistances between the different simulation points.The formulas are as follows: where λ is the latent heat of vaporization (J kg −1 ); λE is the latent heat deduced from the ET (W m −2 ); e * (T ) is the saturated vapor pressure at temperature T (Pa); e a is the canopy air space vapor pressure (Pa); ρ is the density of air (kg m −3 ); c p is the specific heat of air (J kg −1 K −1 ); γ is the psychrometric constant (Pa K −1 ); W c is the fractional wetted area of the canopy; W g is the fractional wetted area of the ground surface; g c is the canopy conductance (m s −1 ), a parameter associated with the biological processes and the vegetation growth environment (e.g., the water potential of the root zone and the temperature); r b is the bulk canopy boundary layer resistance (s m −1 ), which is a function of the wind speed, the temperature, the canopy structure and other factors and is considered an energy-related parameter; r d is the aerodynamic resistance between the ground and the canopy air space (s m −1 ) and is related to the wind speed, the ground surface roughness, and other factors; h soil is the relative humidity of the soil pore space; and r soil is the soil surface resistance (s m −1 ), representing the impact of the soil on the water vapor diffusion.
It can be seen from the above description that the ET simulation in the SiB2 model is based on a physical mechanism in which the impacts of the water, energy and biological processes are all considered.However, in the SiB2 model, the description of the water movement in the subsurface is relatively simple, and the soil water movement is limited to a shallow unsaturated zone; the GWT and lateral flows are not considered in the model.These simplifications in the water cycle will eventually cause errors in the energy cycle calculation and the biological process simulation.Coupling the LSM and the GWM provides a good approach to improving the accuracy of the simulated model.

Coupled model approach
As described above, the SiB2 and AquiferFlow models have strengths and weaknesses.Coupling the two models can overcome their weaknesses and more accurately simulate the water and energy cycles.In the study presented here, the two models were tightly coupled from the model codes, and a new model, named GWSiB, was developed.We will introduce the coupling mechanism in the following section.
In the coupled model, SiB2 simulates the energy balance, the vegetation root water uptake and the hydrologic processes above the ground surface, and AquiferFlow simulates water movement in the subsurface, including the saturated and unsaturated zones.Specifically, the SiB2 model is used to calculate the precipitation infiltration (Q 1 ), the moisture evaporation (E gs ) and the transpiration (E ct ) based on the energy balance and the water balance.The calculated results are used as the sinks and sources (ε) and are input into Aquifer-Flow to calculate the groundwater potential (ψ).The obtained water potential is then used to calculate the groundwater movement in the model grids.The new groundwater condition obtained is transferred back to SiB2 to complete the calculation cycle in one time step.A flowchart of the coupling procedure is illustrated in Fig. 1.
The coupling of the two models includes spatial and temporal coupling.First, we discuss coupling the two models in space.The SiB2 model, a vertical 1-D model, must be extended horizontally to match AquiferFlow, a 3-D model.In our study, the mesh of the coupled model uses the Aquifer-Flow scheme in the horizontal dimensions, which means that on every topmost cell of the AquiferFlow grid, a SiB2 simulation is built.In the vertical space, AquiferFlow has a more flexible layer structure than SiB2; consequently, the three subsurface layers in SiB2 are preserved, and the three top layers in AquiferFlow are set to be consistent with them.The infiltration and the soil evaporation are linked with the top layer of AquiferFlow, and the root zone uptake is linked with the second layer.
Although the surface runoff (R 0 ) and the base flow (Q 3 ) are calculated on a vertical column in SiB2, the surface water convergence between cells is not taken into account in the coupled model.This simplification will not cause a significant deviation when the model is used in the middle or lower reaches of an arid or semi-arid basin, because there are almost no flow confluence processes in these regions.However, if the model is used in the upper reach of a basin, the errors cannot be ignored.Wang et al. (2009) handled this problem by coupling SiB2 with a geomorphology-based hydrological model (GBHM).In our study, the model validation and tests were performed in the middle reach of the Heihe River basin where the surface runoff is not the key hydrological process; consequently, the coupled model can be used here.
We now discuss how to handle the temporal discretization of the coupled models.LSMs usually use a time step of 1 h or less, because the energy and mass variables simulated by the LSM, such as the soil surface temperature, vary rapidly and vary significantly from day to night.However, the groundwater head and flow vary more slowly; therefore, the time intervals used in GWMs are usually one day or longer.LSMs are more sensitive to the time resolution than GWMs and generally cannot accept time intervals greater than 1 h.If the time step in an LSM is one day, the temporal fluctuations that occur in hours will smooth out, generating significant calculation errors and even making the simulation meaningless.A shorter time interval will not significantly affect the simulation accuracy of GWMs but will significantly reduce their computational efficiency.
Considering the time steps used in the two models, two alternative time coupling schemes are implemented in the coupled model.One scheme is to use the AquiferFlow time step (which is the same as that used in SiB2), which is set to be either 1 h or 30 min.The second scheme is to adopt a time step of one day in AquiferFlow, which is the time step normally used in the GWMs, while using a time step of 1 h in SiB2.The fluxes that are accumulated over one day in SiB2 are then exchanged with AquiferFlow.The first scheme has a higher precision but requires more computation; it is thus suitable for theoretical analysis or small-scale simulation.The second scheme greatly improves the computational efficiency and achieves acceptable calculation accuracy.This scheme is more suitable for large area simulation.
Both AquiferFlow and SiB2 use Richards equation as their control equations for soil water movement, but they adopt different parameterization schemes to describe the relationship between the unsaturated hydraulic conductivity (K) and the soil moisture potential ( ).The Gardner and Fireman (1958) method is used in AquiferFlow, and the Clapp and Hornberger (1978) scheme is used in SiB2.The different schemes would make a huge difference in the calculation of Richards equation.Because the different parameter schemes can strongly affect the calculation results of Richards equation, a discontinuity would be created in the soil moisture at the two model communication times, especially when using the second scheme, in which the water accumulation is exchanged between the GWM and the LSM.To solve this problem, the Clapp and Hornberger (1978) soil moisture scheme used in SiB2 is introduced into the Aquifer-Flow model framework.In AquiferFlow, the relative permeability (K r ) and the effective saturation (S e ) are the two key parameters for soil moisture movement and content.These two parameters are defined as fractions in AquiferFlow and are used to control the moisture movement in the unsaturated zone by adjusting the saturated hydraulic conductivity (K s ) and the saturated soil water content (θ s ), which is approximately equal to the porosity (φ).The parameters used in Clapp and Hornberger's (1978) soil moisture scheme are also based on the saturated moisture potential (ψ s ) and the saturated hydraulic conductivity (K s ); these equations can thus be transformed to the AquiferFlow framework as Replacing the K r and S e parameters of the AquiferFlow model by Eqs. ( 10) and ( 11) makes the vadose zone parameters of AquiferFlow consistent with SiB2 and reduces the soil moisture discontinuities in the model at the time of coupling.
After the GWSiB model is built, a sensitivity test about the key parameters of the model is performed, and the model validation is conducted based on the measured data of two observation stations.The model is then applied to the middle reach of the Heihe River basin to test the applicability of the model on the regional scale.The following sections will describe the model validation in detail.

Model validation
Many processes, such as the water cycle, the energy cycle, and biochemical processes, are integrated into the GWSiB model.Because the ET is determined by the combined effects of these processes, the validation of the simulated ET can assess whether these processes are adequately simulated.In this paper, ET is used to validate the GWSiB model.The main feature of the GWSiB model is that the 3-D groundwater movement is added to the land surface model.This makes it possible to simulate the impact of the saturated groundwater level and the lateral groundwater flow on the land surface processes.To analyze the relationship between the groundwater and the ET, a sensitivity test is performed prior to the model validation.

Sensitivity test
A synthetic domain is used to perform the sensitivity analysis of the GWSiB model.In this domain, there are two rivers separated by 200 m, with a platform located between the two rivers.The altitude of the platform is 1500 m, and its soil textures are homogeneous and isotropic.The rivers can recharge the groundwater of the platform through lateral flow, so the river levels are generally representative of the groundwater level in the area.The synthetic domain is divided into an array of uniform vertical columns extending 3 columns wide parallel to the rivers and 10 columns long between the rivers; each column is 5 m wide and 20 m long (Fig. 2).In the vertical direction, the soil is divided into four layers representing the surface soil layer, the root layer, the deep soil layer and the phreatic aquifer layer, at 0.02, 0.48, 1.5, and 50 m, respectively.Consequently, the simulated domain forms a 120element grid (10 × 3 × 4).
The model structure allows the groundwater level in the model to be easily controlled by directly changing the water level of the two rivers.The model structure can also incorporate the characteristics of the lateral groundwater flow.
The forcing data used were measured at the Linze grassland station (LZG) of the Heihe River basin at 17:00 LT (local time) on 12 August 2008; these data represent typical moderate-radiation atmospheric conditions in this arid region.The forcing data are maintained constant throughout the simulation period in the sensitivity test.
The land cover of the simulated area is assumed to be grassland.The vegetation parameters used are the default parameters for the "short vegetation/C4 grassland" vegetation type, which is one of the nine types of vegetation derived from Sellers et al. (1996a) that are defined in SiB2.The leaf area index (LAI) used to represent vegetation growth is kept at a constant value of 2 during the simulation period.In the model, the two rivers are defined as fixed head boundaries.The other groundwater boundaries are defined using the no-flow condition.The river levels, the soil texture and the groundwater hydraulic parameters are specified according to the sensitivity tests described in detail in the following.The first time-coupling scheme described in Sect.2.3 is used in the sensitivity test, and the time step is set to 1 h.The simulation period is 1488 h.
Using the model, the impact of the groundwater depth on the ET is first analyzed.Five groundwater depths are simulated: 1.5, 3.0, 5.0, 8.0, and 10.0 m (corresponding to river levels of 1498.5, 1497.0,1495.0,1492.0, and 1490.0 m, respectively).In this experiment, the soil texture is set to 30 % clay, 30 % silt, and 40 % sand.The results from the cell located in the center of the platform are used for the analysis.The analysis results are shown in Fig. 3.
The results show that the ET decreases as the groundwater depth increases.There is a difference of up to 80 % in the ET between the 10 and 1.5 m groundwater depths in this experiment.This is because a lower groundwater level reduces the water supply from the saturated groundwater to the surface and the root soil layer by capillary action, and the reduction in the surface water limits ET.The simulated ET also decreases continuously with time, except in the case of the 1.5 m groundwater depth.We believe this is due to the velocity of groundwater flow, which includes the velocity of the lateral flow and the vertical flow.When the water lost in ET is greater than the water gained through groundwater recharge, the soil gradually dries, and ET decreases.However, there is a sufficient water supply in the case of the 1.5 m groundwater depth, so ET does not significantly decline during the simulation period.
This analysis shows that the groundwater flow can significantly affect ET.In the GWSiB model, the flow characteristics of the groundwater are determined by the soil texture and calculated using an empirical formula (Yang et al., 2005).Consequently, in the second experiment, the effect of different soil textures on the ET is analyzed.Three soil texture types, including sand-based soil (clay: 20 %, silt: 20 %, sand: 60 %), silt-based soil (clay: 20 %, silt: 60 %, sand: 20 %), and clay-based soil (clay: 60 %, silt: 20 %, sand: 20 %), are used in this experiment.The relevant parameters defined by Clapp and Hornberger (1978) are K s = 0.66 m day −1 , s = 0.12 m, and B = 6.09 for sandbased soil; K s = 0.16 m day −1 , s = 0.41 m, and B = 6.09 for silt-based soil; and K s = 0.16 m day −1 , s = 0.41 m, and B = 12.45 for clay-based soil.In this experiment, the groundwater level is set to 3 m.The results are shown in Fig. 4.
Figure 4 shows that the hydraulic characteristics of the groundwater have a significant effect on ET.The slope of the decrease in ET of the sand-based soil is less than the other soil types.This is because the sand-based soil has a higher hydraulic conductivity than the other soil types, and the water in the soil that is lost by ET can be recovered quickly.In contrast, the clay-based and silt-based soils have greater capillary action (i.e., these soil types have a higher saturation moisture potential).These soil types can provide more water to the ET process early in the simulation, but, as the soil moisture decreases, the water in the soil cannot be quickly replenished, and the rate of ET decreases rapidly.

Model validation
Based on the sensitivity analysis model, the GWSiB is validated at the Linze grassland station (LZG) and the National Observatory on Climatology at Zhangye (ZYNOC), which represent shallow and deep groundwater conditions, respectively.The measured data from the two stations, including atmospheric driving data, vegetation data, soil textures and groundwater levels, are input to the model as the true value.The vegetation parameters are calibrated for each station before the model validation is performed.The detailed process is as follows.

Validation of the model for the Linze grassland station (LZG)
The LZG is located in the middle reach of the Heihe River basin in the northwest of China.The longitude is 100.07 • E, and the latitude is 39.25 • N. The land cover in the LZG is mainly wetland, grassland and salinized meadow.An automatic meteorological station (AMS) built by the Watershed Allied Telemetry Experimental Research (WATER) project (Li et al., 2009) in the LZG was used for observations from 1 October 2007 to 27 October 2008.The AMS provides all the necessary atmospheric forcing data for our modeling study.Although there is no direct measurement of the latent heat at the LZG station, it can be obtained from the sensible heat by the energy balance equation: where λE is the latent heat (W m −2 ); λ is the heat of vaporization (J kg −1 ); E is the ET (m); R n is the net radiation (W m −2 ), equal to the difference of the downward radiation and the upward radiation, which can be obtained from the atmospheric forcing data; and H is the sensible heat.In the WATER experiment, a large-aperture scintillometer (LAS) flux system was used from 19 May 2008 to 27 August 2008 to obtain the sensible heat data for the model.G is the ground heat flux (W m −2 ) and is assumed to be proportional to the net radiation (Su, 2002): where c = 0.05 for a full vegetation canopy, s = 0.315 for bare soil, and f c is the fractional canopy coverage, which is set to 0.81 based on observations in the LZG.The latent heat of the LZG is calculated according to a variety of observational data and the energy balance equation; the ET is then deduced based on the latent heat.varies between 1.2 and 1.9 m during the simulation period.The vegetation type used in the LZG model is "C3 grassland", and the parameters from Sellers et al. (1996a) for this vegetation type are calibrated according to the ET measured in the LZG from 19 May to 1 July 2008.The parameters used in the model are listed in Table 1.In the GWSiB model, the LAI has often been used to characterize the vegetation growth process.The LAI data for the LZG were obtained from MODIS global LAI and FPAR (fraction of photosynthetically active radiation) products (MCD15A3) with a 4day time resolution and a 1-km spatial resolution.These data were revised according to observation and interpolated to the 1-h time resolution used in the model.The LAI varies from 2.1 to 3.4 during the simulation period.
The same model structure as in the sensitivity test is used in the model validation, including the same spatial structure and the same time step.The data described above are used as the input to the GWSiB model to simulate the energy and water cycles of the LZG from 19 May 2008 to 27 August 2008.In the simulation, the data from 19 May to 1 July are used for the calibration of the model parameters, and the data from 2 July to 27 August are used to validate the model of the LZG.For comparison, the SiB2 model is executed using the same conditions.The simulation results are shown in Fig. 5.
The simulation results show that the GWSiB simulation agrees well with the observed ET trends and magnitudes.However, the SiB2 simulation significantly underestimates the ET in the LZG.These findings can also be observed in the scatter plot (Fig. 6), which shows the relationship between the observed data and the results of the two models.The simulated results from the SiB2 model are far below the 1 : 1 regression line.This underestimation is confirmed by the statistical analysis.The mean value of the observed ET is 3.93 mm per day during the simulation period, and the mean value from the GWSiB model is 3.98 mm per day; the mean relative error (MRE) between them is 1.4 %.However, the mean value from the SiB2 model is 0.76 mm per day, and the MRE is 80.7 %.
The GWSiB model can provide a more realistic simulation than the SiB2 model, because the movement of the groundwater is taken into account in the GWSiB model but not in the SiB2 model.The lateral flow of groundwater causes groundwater accumulation in the shallow-water region and raises the groundwater level.The saturated groundwater supplies water to the soil near the ground surface by capillary action, leading to greater ET at the land surface in the LZG.Compared with the GWSiB model, the SiB2 model is a vertical 1-D model that cannot model the process of groundwater recharge by lateral water movement; consequently, the soil water moisture of the land surface is underestimated, and the lower soil moisture reduces the ET observed in the SiB2 model.

Validation of the model for the National Observatory on Climatology at Zhangye
The Fig. 8.A scatter plot comparing the GWSiB-simulated and the SiB2-simulated evapotranspiration with the measured evapotranspiration in the National Observatory on Climatology station at Zhangye.Sellers et al. (1996a) for this vegetation type are calibrated according to the ET measurements.The typical soil textures of the ZYNOC obtained from field measurements are 23 % clay, 30 % sand and 47 % silt.The groundwater level data come from a well approximately 2 km from ZYNOC.The variation of the groundwater depth is not significant during the simulation period; the depth varies from 25.4 to 26.2 m.The groundwater levels are used as fixed-head boundary conditions in the model.The LAI data are obtained from the MODIS products (MCD15A3), in a manner similar to the validation of the LZG.The LAI of the ZYNOC increases from 0.1 to 0.4 during the simulation period.
The setup of the model for the ZYNOC validation is the same as that used for the LZG validation.The simulation period is from 28 June 2008 to 22 August 2008.The ET data from 28 June to 21 July are used to calibrate the vegetation parameters of the model.The parameters that are held constant throughout the simulation period are listed in Table 1.The remainder of the simulation period is used to validate the model.The ET processes are simulated in the SiB2 model using the same conditions, except that the groundwater level is not considered in the SiB2 model.The simulation results from the two models are shown in Figs.7 and 8.
Figures 7 and 8 show that the two models produce similar results for the ZYNOC.Both models provide a good simulation of the magnitude and the variation of the ET process, and there is no significant difference between the results from the two models.The mean value of the observed ET is 0.72 mm per day during the simulated period, and the mean value from the GWSiB model is 0.74 mm per day; the MRE between the observation and the simulation is 3.0 %.The mean value from the SiB2 model is 0.69 mm per day, and the MRE is 4.7 %.
We believe that the groundwater level and the lateral flow of the groundwater have a small effect on the energy and water cycles on the land surface because the groundwater is deep (approximately 26 m below ground surface).In the thick vadose zone, water movement occurs primarily in the vertical dimension.Both the SiB2 and the GWSiB models can simulate this process adequately, and they provide similar and realistic results.The GWSiB and SiB2 results nevertheless show some minor differences; e.g., the variations in the simulated ET are greater in the SiB2 model than in the GWSiB model, and the SiB2 model tends to produce larger or smaller extreme values.We believe that the difference in the results can be attributed to the difference in the two model structures.

Regional test of the model
The validation of the GWSiB model in the shallow groundwater site (LZG) and the deep groundwater site (ZYNOC) demonstrates the adequacy of the model.However, the water cycle can usually only be understood comprehensively on a regional scale.Because the GWSiB model is fundamentally a 3-D model, it has the ability to simulate regional water and energy cycles.In the following, the GWSiB is applied in the middle reach of the Heihe River basin to test the simulation capabilities of the model in the region.

Study area
The middle reach of the Heihe River basin is an arid inland river basin located in northwestern China.In this basin, an integral groundwater cell, which is hydrogeologically described as the Zhangye basin, is selected as the study area.The latitude of the study area ranges from 38.7 • N to 39.8 • N, the longitude from 98.5 • E to 102 • E, and the total area is approximately 12 825 km 2 (Fig. 9).Irrigated agriculture, grassland/steppe, wetland, and Gobi are the main types of land cover in the study area.Some numerical groundwater simulation studies have been performed of this area, such as studies by Su (2005) and Wen et al. (2007), who used the FEFLOW model to simulate the groundwater movement in the area and forecast the groundwater trends.Hu et al. (2007) developed a 3-D GWM and used it to study the groundwater interaction with rivers and springs in the area.Ding et al. (2009) developed a 2-D numerical model to simulate the groundwater dynamics of the area.Zhou et al. (2011) built a GWM of the area to quantify the effects of land use and anthropogenic activities on the groundwater system.These groundwater modeling studies provide excellent references for this study.

Model settings
The study area is uniformly discretized into 79 rows and 32 columns horizontally, and each numerical cell has dimensions of 3 km × 3 km.In the vertical direction, the soil below the ground surface is divided into six layers.The upper three layers correspond to the soil layers of the SiB2 model and represent the surface soil, the soil root, and the deeper soil layers.The thickness of these soil layers is set to 0.02, 0.48, and 1.5 m, respectively, according to the average conditions of the soils in the middle reach of the Heihe River basin.The lower three layers are used to describe the hydrogeologic structure in the study area, representing the unconfined aquifer, the aquitard and the confined aquifer.The thicknesses of the lower three layers are determined by the interpretation of the logging data obtained from 108 boreholes in the region, as proposed by Zhou et al. (1990).The study area contains a total of 15 168 (79 × 32 × 6) cells, as shown in Fig. 9.
The topography of the study area is determined from 90m resolution digital elevation model (DEM) data obtained from the Shuttle Radar Topography Mission (SRTM) and upscaled to 3-km resolution.The initial GWT distribution in the study area is obtained by the interpolation of GWT measurements conducted in December 2003 from 36 observation wells in the study area.Any GWT positions not available from the measurements are determined from the relevant literature (Su, 2005;Hu et al., 2007;Wen et al., 2007;Zhou et al., 2011).The initial GWT data show that the main direction of the groundwater flow of the middle reach of the Heihe River basin is from south to north and is roughly consistent with the river flow direction.The GWT is lower in the north of the study area, where the groundwater discharges to the river.The boundary conditions and the saturated hydraulic conductivity (K s ) of the model were initially assigned values according to previous study results from the Heihe River basin (Su, 2005;Hu et al., 2007;Wen et al., 2007;Ding et al., 2009;Zhou et al., 2011).Later, these parameters were optimized through trial-and-error calculations using GWT data obtained in January 2008.Ultimately, the boundary conditions of the model are set as fixed-flow conditions: the southern boundary has 1.62 × 10 8 m 3 water inflow every year; the northern boundary has 0.37 × 10 8 m 3 a −1 water inflow; and the western and eastern boundaries have a total of 0.08 × 10 8 m 3 a −1 water inflow.These water inflows were allocated to each of the active cells of the boundary grid in the model.The saturated hydraulic conductivity (K s ) field in the study area was divided into 24 sub-regions, with values ranging from 0.5 m d −1 to 20 m d −1 .The distribution of the specific storage (S s ) was represented by 10 sub-regions, with values ranging from 0.003 m −1 to 0.17 m −1 .The water potential parameters of the unsaturated zone were determined according to the soil texture.The parameter values for the soil characteristics used in this study were obtained through the analysis of the Chinese dataset of the multi-layer soil-particle size distribution, which has a 1-km resolution (Shangguan et al., 2012).
The atmospheric data used in the model, including the incident solar radiation, incident longwave radiation, wind speed, air pressure, vapor pressure, air temperature, and precipitation, are taken from the Global Land Data Assimilation System (GLDAS) project (Rodell et al., 2004).The spatial resolution of the original data is 25 km, and the temporal resolution is 3 h.The data are interpolated to a spatial resolution of 3 km and a temporal resolution of 1 h to fit the resolution of the coupled model.The temporal interpolations of the data were performed using a statistical method provided by the Global Soil Wetness Project 2 (GSWP2) (for the precipitation data) and the cubic spline method (Dai et al., 2003) (for other data).The high-resolution meteorological interpolation model MicroMet (Liston and Elder, 2006) is used in the spatial interpolation of the data.These data are used in all of the model cells except those cells where an AMS is located.In those cells, the interpolated atmospheric data are replaced by the measured data.
The land cover data used in the model are derived from the Multi-source Integrated Chinese Land Cover (MICL Cover) data (Ran et al., 2012).The International Geosphere-Biosphere Programme (IGBP) land cover classification system is used in the MICL Cover data, and the land surface is classified into 17 types.In this study, the 17 types are grouped into nine types corresponding to the vegetation classifications of the SiB2 model (Sellers et al., 1996b).The parameters of each vegetation type defined by Sellers et al. (1996a) are calibrated according to the measured ET data from this study.
There is a considerable amount of farmland in the study area, and the energy and water cycles are affected by the irrigation.In this model, the irrigation is modeled by sink terms added to the groundwater system on the grid cells at which the vegetation type of cells is defined as "agricultural".The irrigation data used in the model come from the local administrative department of agriculture.
The time-dependent vegetation parameters used in the coupled model were obtained from satellite data.The level-4 combined (Terra and Aqua) MODIS global LAI and FPAR products (MCD15A3), which are provided every 4 days at a resolution of 1 km, are linearly interpolated to a temporal scale of 1 h and are resampled to a spatial resolution of 3 km for the coupled model.
Because of the scope of the study area and the quantity of cells, the computational effort required in the simulation is considerable.Therefore, the second time-coupling scheme is used in the case study; i.e., an hourly time step is used in the SiB2 model and a daily time step is used in AquiferFlow.At the start of each simulated day, the two models exchange fluxes.Because the groundwater flow equation is highly nonlinear, the convergence of the groundwater model needs careful adjustment of the slover.In this study, the model convergence was implemented through the adjustment of the relaxation factor, the number of iterations, and the convergence criterium, which were set to 1.3, 10 000, and 0.001 m, respectively.
The initial conditions of the model, such as the soil moisture, surface temperature, canopy temperature and other variables, are determined according to the general conditions in the middle reach of the Heihe River basin in winter.The initial conditions are specified in the coupled model for 1 January 2004.The model is run for 4 yr as a spin-up from 1 January 2004 to 1 January 2008 to ensure that the model is initially in equilibrium.The simulated values at the end of the spin-up period (1 January 2008) are treated as the initial conditions of the simulation period for the coupled model.

Analysis of the model results
Using the GWSiB model of the middle reach of the Heihe River basin, the energy budget and the water movements in this region were simulated from 1 January to 31 December 2008.We analyzed the results of the model to test the model on the regional scale.First, the simulated GWT data from 20 December 2008 are compared with the measured GWT data for the same day obtained from the interpolation of the data from the 36 observation wells.The results are shown in Fig. 10.
The simulation results and the measurements agree well (Fig. 10) except on the west side of the study area, where the groundwater depth is greater than 100 m and there are few groundwater observation wells.The initial GWT data for these regions are taken from the existing literature (Su, 2005;Wen et al., 2007), and the measured data for these regions are extrapolated from the 36 observation wells.We believe the uncertainty of the GWT data for these areas is the main reason for the simulation errors.Additionally, in the upper section of the Heihe River, the simulated GWT results are higher than the measured GWT data (e.g., at the 1450 m GWT) in Fig. 10.This is because the upper section of the Heihe River is considered to be the region of surface water recharge to the groundwater in the coupled model.The significant water infiltration in this region and the partial GWT increase are simulated in the model, but this change is not caught in the measurements because there are few observed wells in this region to provide data for the GWT interpolation.In general, the GWSiB simulation produces a good model of the groundwater conditions of the middle reach of the Heihe River basin, and the GWSiB model can provide a continuous and time-varying depiction of the impact of the groundwater on the land surface energy and water cycles.
The regional ET is analyzed next.In this test, the GWSiB model is used to calculate the hourly ET for the study area in 2008.Because there is no way to directly measure the regional ET, we compare the GWSiB simulation results with ET data obtained by remote sensing.Li et al. (2012) calculated ET from the NOAA/AVHRR satellite data of the Heihe River basin and obtained a regional ET with a 1-km resolution at 15:00 LT on 2 August 2008.The regional ET simulated by GWSiB for the same time is compared with these data in Fig. 11.The results of the model simulation and the remote sensing calculation have a very similar spatial distribution, and both show the ET distribution characteristics of the middle reach of the Heihe River basin.In this region, the banks of the Heihe River and the irrigated area have greater ET because the river recharges the groundwater through lateral flow and raises the groundwater levels of the banks while the irrigation increases the soil moisture of the irrigated area.The average ET in the study area determined from the model simulation is 0.24 mm, and the value calculated from the remote sensing data is 0.31 mm.There is a 23 % MRE between the two values.However, some details are not consistent between the two results.The main reason for the inconsistency, besides the error caused by the remote sensing terms, is the model uncertainties, including uncertainty in the meteorological driving data, the land cover data, the soil texture data, the groundwater data, and the grid structure of the model.
On the whole, the GWSiB model can simulate regional energy and water cycles, although there are some errors.In this test, the GWSiB model achieves an acceptable regional ET simulation of the middle reach of the Heihe River basin both in terms of absolute values as spatial distribution.

Discussion and conclusions
LSMs can describe the land surface energy and water cycles well, but the water movement in the subsurface is oversimplified, and the movement of saturated groundwater is ignored.Conversely, GWMs can describe the dynamic movement of subsurface water, but they cannot simulate the physical mechanisms of ET, which is an important component of the water cycle because this process involves the energy cycle and the biological processes.Coupling the two types of models can effectively overcome their respective shortcomings, and, by linking the energy and water cycles together, these processes can be simulated more comprehensively and potentially more accurately.
In this study, a 3-D dynamic groundwater model, Aquifer-Flow, and a typical land surface model, SiB2, are fully coupled.In the coupling scheme, infiltration, evaporation and transpiration, which are simulated by the SiB2 model, are used as inputs to the AquiferFlow model, and the soil moisture values calculated by AquiferFlow are used in SiB2.In the sensitivity analysis of the coupled model, the effects of the groundwater level and the hydraulic parameters of the groundwater on the energy and water cycles are analyzed.The excellent performance of the GWSiB model in the LZG and ZYNOC validation studies demonstrates that the coupled model has the proper structure.Additionally, the case study of the middle reach of the Heihe River basin demonstrates the ability of the GWSiB model to simulate regional-scale dynamics.
There are many factors that can affect the accuracy of the simulation of the ET in a region, such as the meteorological driving data, the groundwater conditions determined by the groundwater hydraulic parameters and boundary conditions, the vegetation type and parameters, the dynamic variation of the vegetation, the soil texture, and the resolution of the model grid.In the simulation of the middle reach of the Heihe River basin, after we adjust the hydraulic conductivities of the groundwater by the trial-and-error method, calibrate the vegetation parameters according to in situ observations, and select the appropriate land cover and soil texture databases, we achieve good simulation results, as demonstrated by the comparison with the remote sensing data.However, the regional simulation using the model is directly affected by the input data in addition to the physical structure of the model.We believe that the model errors will be gradually reduced as the accuracy of the regional data improves.
In our study, the ET is used as the validation variable because the ET is the key variable in the energy and water cycles, and it provides an integrated measure of the accuracy of the simulation of the energy cycle, the water movement and vegetation photosynthesis.However, the GWSiB model is not limited to the calculation of ET.The GWSiB model can describe the energy and water processes as a system, and the model can be widely used in the study of Earth science.
From our study, the following four conclusions can be obtained.
1.The groundwater depth can significantly affect the ET on the land surface; the ET increases as the groundwater depth decreases.Additionally, the groundwater hydraulic parameters and the soil structure can affect ET through the vertical and lateral movement of the groundwater.
2. In a shallow groundwater depth zone, the GWSiB model, which incorporates the groundwater movement, simulates the ET process on the land surface more accurately than the SiB2 model, in which the groundwater movement is not simulated.The ET will be underestimated if the groundwater movement is ignored in this region.
3. The interaction of the groundwater and the land surface processes is weak in zones of large depth to groundwater, and the subsurface water movement is dominated by vertical movement under these conditions.The GWSiB model produces results similar to the SiB2 model, and each of the models can simulate the ET process in this region well.
4. The GWSiB model can simulate regional energy and water cycles.The GWSiB simulation accurately models the middle reach of the Heihe River basin.The accuracy not only depends on the correctness of the model structure but is also directly affected by the model data.
In summary, the coupling of the groundwater and land surface models allows the land surface and subsurface processes to be simulated as a system, and the coupled model can depict the interaction of the groundwater and the energy and water movement on the land surface.This improves the simulation of the energy and water cycles.Although a coupled model was developed in this study, some the energy and water processes, such as surface water processes, water resource allocation, and soil freezing and thawing processes, are not considered in the model.Additionally, the validation of the GWSiB model is limited by the shortage of available data.Further validation and improvement of the coupled model will be the primary thrust of future studies.

Fig. 1 .
Fig. 1.The coupling between SiB2 and AquiferFlow in GWSiB.The bold solid lines with 2

Fig. 2 .Fig. 2 .
Fig. 2. The synthetic model structure used in the sensitivity analysis of the GWSiB.The results 3

Fig. 4 .
Fig. 4. Simulated ET results at different soil texture conditions.3 4 Fig. 4. Simulated ET results at different soil texture conditions.

Fig. 5 .Fig. 5 .
Fig. 5.The measured evapotranspiration, the GWSiB-simulated evapotranspiration, and the 3 SiB2-simulated evapotranspiration for the Linze grassland station from May 19, 2008 to August 4 27, 2008.The data in the shadowed part of the graph are used for the model parameter calibration, 5 while the other data are used for the model validation.6 7

8 Fig. 7 .
Fig. 7.The measured evapotranspiration, the GWSiB-simulated evapotranspiration, and the 3 SiB2-simulated evapotranspiration for the National Observatory on Climatology station at 4 Zhangye from June 28, 2008 to August 22, 2008.The data in the shadowed part of the graph are 5 used for the model parameter calibration, while the other data are used for the model validation.6 7 8 Fig. 8.A scatter plot comparing the GWSiB-simulated and the SiB2-simulated evapotranspiration 2 with the measured evapotranspiration in the National Observatory on Climatology station at 3 Zhangye.4

Fig. 9 .
Fig. 9.The location of the study area in the middle reach of the Heihe River basin and the model 2

7 Fig. 10 .
Fig. 10.The observed and simulated groundwater table (GWT) of the study area in December 3

Fig. 11 .Fig. 11 .
Fig. 11.Evapotranspiration calculated by remote sensing (a) and simulated by the GWSiB model 3 (b) in the middle reach of the Heihe River basin at 15:00 local time on August 8, 2008.The map (a) 4 is derived from Li et al. (2011).5 6

Table 1 .
The soil texture and groundwater depth data come from observations in the LZG.The soil is composed of 13.4 % clay, 31.6 % sand and 55 % silt.The groundwater depth The vegetation parameters for each station.