Gravitational and capillary soil moisture dynamics for distributed hydrologic models

Distributed and continuous catchment models are used to simulate water and energy balance and fluxes across varied topography and landscape. The landscape is discretized into computational plan elements at resolutions of 10–10 m, and soil moisture is the hydrologic state variable. At the local scale, the vertical soil moisture dynamics link hydrologic fluxes and provide continuity in time. In catchment models these local-scale processes are modeled using 1-D soil columns that are discretized into layers that are usually 10–10 m in thickness. This creates a mismatch between the horizontal and vertical scales. For applications across large domains and in ensemble mode, this treatment can be a limiting factor due to its high computational demand. This study compares continuous multi-year simulations of soil moisture at the local scale using (i) a 1pixel version of a distributed catchment hydrologic model and (ii) a benchmark detailed soil water physics solver. The distributed model uses a single soil layer with a novel dualpore structure and employs linear parameterization of infiltration and some other fluxes. The detailed solver uses multiple soil layers and employs nonlinear soil physics relations to model flow in unsaturated soils. Using two sites with different climates (semiarid and sub-humid), it is shown that the efficient parameterization in the distributed model captures the essential dynamics of the detailed solver.


Introduction
Soil moisture controls the partitioning of rainfall into infiltration and runoff, and it controls land surface temperature through its effect on the partitioning of available energy into sensible and latent heat fluxes.It is the hydrologic state variable, together with land temperature, in models of surface water and energy balance.The state dynamics are affected by hydrometeorological forcing of precipitation, radiation, and atmospheric evaporative demand.Furthermore, topography, land use, and soil properties across the landscape, affect soil moisture temporal evolution (Western and Grayson, 2000;Lawrence and Hornberger, 2007;Vereecken et al., 2007;Ivanov et al., 2010;Liu et al., 2012;Beven and Germann, 2013).
There are diverse methods for measuring soil moisture e.g dielectric-and heat dissipation-based approaches.The suitability of a certain method or system depends largely on the desired scale, accuracy, and resolution, both in space and in time.Unfortunately, all current observing systems have their shortcomings.For instance, in situ sensors can provide high accuracy and fine temporal resolution but at limited spatial footprint, sampling campaigns can provide better spatial resolution and coverage but at low sampling frequency and duration, while space-borne remote sensing platforms provide global spatial coverage for surface soil moisture sensing but at coarse spatial resolution and with infrequent revisits.
Numerical hydrologic models fill some of the shortcomings of observations.Incoming radiation and precipitation are used in conjunction with water and energy balance mod-Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Castillo et al.: Gravitational and capillary soil moisture dynamics els to simulate the evolution of soil moisture in the vadose zone and estimate the water and energy fluxes across the landscape.Harter and Hopmans (2004) describes how hydrologic models have traditionally been used by two largely disconnected groups: the watershed hydrologists (and recently also climate modelers), who deal with macro-processes, and the soil physicists, who study soil properties and states at the laboratory or local to plot scales.Watershed hydrologists have traditionally used lumped or semi-distributed models that treat the vadose zone as a 0-D black box.The computational time step is usually hourly, daily, or even longer.Two examples of heritage models used by watershed hydrologists are the semi-distributed models TOPMODEL (Beven and Kirby, 1979) and SAC-SMA (Burnash et al., 1973), both of which have been demonstrated as highly capable of simulating streamflow.Meanwhile, soil physicists, who have detailed measurements of soil properties and states at the local to plot scales, model unsaturated flow by discretizing the hydrologically active soil column into layers that are usually 10 −3 to 10 −1 m in thickness, using the Richards equation (RE) which can be written as where K is the hydraulic conductivity, ψ is pressure head, z is elevation with respect to a datum, θ is soil moisture, and t is time.For stability, this nonlinear partial differential equation is solved using sub-hourly time steps.Over the years, the modeling efforts of the two disciplines have started to converge, as manifested by the emergence of physically based distributed hydrologic models (DHMs).These models discretize the landscape in computational elements that are 10 1 to 10 3 m in the horizontal.Adopting the practice in soil physics, many DHMs employ RE and discretize the hydrologically active soil layer into vertical layers that are 10 −3 to 10 −1 m thick.Some DHMs that use RE include MIKE-SHE (Refsgaard and Storm, 1995) and ParFlow (Ashby and Falgout, 1996), which use grids for horizontal discretization, and PIHM (Qu and Duffy, 2007) and TRIBS (Ivanov et al., 2004), which use triangulated irregular network (TIN) as horizontal elements (see Table 1).More DHMs are discussed by Smith et al. (2004Smith et al. ( , 2012) ) under the context of the Distributed Model Intercomparison Project.There are example studies that demonstrate the advantages of DHMs over lumped and semi-distributed model (Bartholomes and Todini, 2005;Castelli et al., 2009;Smith et al., 2004;Vieux et al., 2004).Although promising, the use of DHMs has its own challenges and criticisms which include (i) the need for a high number of inputs that often should have fine spatiotemporal resolutions, (ii) the use of many parameters which makes the calibration process tedious and raises the concern on equifinality (Beven, 2006), and (iii) the high computational requirement (Smith et al., 2004(Smith et al., , 2012)).
The hydrologically active soil mantle is a thin layer draped over the landscape, and it serves as the intermediate water storage connecting the surface above and the deeper soil layers or groundwater aquifer below.Because of the scale mismatch between the vertical and horizontal discretization of DHMs (millimeters to centimeters in the vertical soil column vs. tens to hundreds of meters in the horizontal), DHMs often treat flow dynamics in the soil as 1-D, i.e., lateral subsurface flow is considered negligible.Exceptions include MIKE-SHE and ParFlow, which can be set up to solve the full 3-D RE.This treatment is, however, very computationally intensive, as demonstrated by Kollet et al. (2010): they utilized 16 384 processors to achieve reasonable run time for ParFlow simulations of a basin on the order of 10 3 km 2 at fine spatial resolution (10 0 to 10 1 m in the horizontal and 10 −2 to 10 −1 m in the vertical).The high computational demand and significantly increased number of parameters to calibrate and state variables to initialize can be limiting factors for applications across large domains and in ensemble mode.
Nonetheless, models based on RE are useful when the vertical profile of soil moisture is desired, especially when the soil column has complex layer sequences or the soil properties are not vertically homogeneous, which is common in real life.Also as mentioned, RE-based models are perhaps appropriate for hillslope, plot, and other small-scale applications, especially when information about the vertical soil structure is available.
However, it is questionable whether RE is an appropriate physical model for watershed and large-scale applications (Beven, 1995;Harter and Hopmans, 2004;Beven and Germann, 2013).Also, using this equation for plan elements that are in the order of 10 1 -10 3 m implicitly assumes that the vertical dynamics of soil moisture at the local scale is scale invariant (up to the limit of the plan element area).Contrarily, field measurements show that soil hydraulic conductivity and pore properties related to the soil retention curve (of ψ) vary significantly both in the horizontal and vertical (Gelhar et al., 1992;Rubin, 2003;Zhang et al., 2004).Furthermore, the review paper of Beven and Germann (2013) argues that the use of RE to model field soil should not be considered physics-based but rather a convenient conceptual approximation.They highlighted the importance of macropores and suggested the use of soil structure with at least two flow pathways.Models that use such structure are the 1-D model of Gerke and van Genuchten (1993), the 1-D model MACRO (Larsbo et al., 2005), and the 1-D or 2-D/3-D model HY-DRUS (Šimůnek and van Genuchten, 2008).In these three models, the soil column is composed of a macropore and a matric compartment, with the water flow in the matric compartment still solved using RE.The inclusion of macropore pathways is dependent on available direct and indirect information on their density and connectivity across the basin.The matric compartment still needs to be characterized in distributed models.
The aim of this study is to test a parsimonious and computationally efficient representation of the near-surface unsaturated zone processes including mass balance and control on We focus on a novel dual-pore parameterized approach.The pore space is divided into gravity and capillary components that each control a different set of hydrologic fluxes; the two are also connected.The partitioning allows the capture of two different time scales in the local-scale soil moisture processes while remaining efficient for applications in distributed hydrologic models (possibly even in ensemble mode).
The novel dual-pore parameterized approach tested is a 1pixel version of the Modello Bilancio Idrologico DIstributo e Continuo (MOBIDIC), a raster-distributed catchment hydrologic model that solves mass and energy balance simultaneously.Table 1 lists the features of MOBIDIC and compares it with some of the hydrologic models that have been mentioned.A key feature of MOBIDIC is its use of a single layer of soil with dual compartments -one for gravitational water and another for capillary-bound water.This representation accounts for both fast and slow processes.At the same time, it makes the model computationally efficient and it reduces the number of state variables in the overall dynamic modeling system.
Division of hillslope soil water into storage that drains under gravitational force and storage that is held under capillary action has been used in diverse applications.The concept of field capacity -variably defined as it may be (drainage after 3 days or water content at a given potential) -has been used in agronomy and irrigation applications.Gravitational water can be considered stored water in the soil above its field capacity.Gravitational water contributes to lateral ex-change and vertical percolation fluxes.It also can fill smaller pores that hold water under capillary action.Capillary water is stored water below the field capacity and can be defined to be limited to water above the residual content.Plant roots and evaporation in general can remove capillary water.Thus gravitational and capillary water dynamics affect different hydrologic fluxes.More recently, Brooks et al. (2009) used water isotope data in a humid catchment field experiment to also distinguish between "tightly bound water" that is used by trees and mobile water that participates in "translatory flow" and enters streams.The conceptualization of the soil matrix into a dual-pore structure with each storage affecting different hydrologic fluxes has been further suggested as a general framework for characterizing hydrologic and ecohydrological response (McDonnell, 2014).White and Toumi (2012) modified a land-surface model to adopt the tightly bound and mobile water parameterization that also each affect different hydrologic fluxes.
In this study we use the gravitational and capillary dualpore approach to modeling soil moisture dynamics in a distributed hydrologic model.We test the fidelity of this approach to local processes by comparing its soil moisture dynamics with that resulting from a benchmark numerical model that solves the vertical heat and moisture dynamics using detailed physics including (Eq. 1).In addition and consequently, since most of the previous applications of MOBIDIC assessed its performance based mainly on streamflow which is an area-integrated flux, this study also tests whether MO-BIDIC is capable of correctly simulating the dynamics of soil moisture, soil temperature, and evapotranspiration (ET).
Although several studies, e.g., Romano et al. (2011), have shown that single bucket-type models can also capture the temporal dynamics of depth-averaged soil moisture, a comparison is not made with this simpler model because we recognize the experimental evidences (e.g., Brooks et al., 2009) and follow the recommendation of, e.g., Beven and Germann (2013) as discussed above to use a model with a dual-pore soil structure.Moreover, although bucket-type models that use single soil moisture state with piece-wise defined functions, e.g., using different dynamics when soil is below or above field capacity, are quite similar to the approach of MO-BIDIC, there are some advantages of explicit representation of gravity and capillary water: processes acting separately on the dual reservoirs can occur simultaneously but not necessarily with predefined relative magnitude.
The paper begins with a description of the catchment hydrologic model MOBIDIC and a description of the 1-pixel version used in this particular study.This is followed by an overview of the selected benchmark model, the established 1-D Simultaneous Heat and Water (SHAW).Then the correspondence between SHAW and MOBIDIC variables, the measures of model performance, and the two study sites are described.

Overview
The MOBIDIC is a physically based and raster-distributed catchment hydrology model that solves mass and energy balance simultaneously.It was developed by Castelli et al. (2009) for basin-scale catchment modeling.This study introduces some modifications to the original parameterization.MOBIDIC uses a single layer for each plan element or soil unit.To account for the different roles of gravity and capillary forces in moving and storing soil water, each soil unit has dual compartments: a gravity reservoir composed of large pores that drain under gravity and a capillary reservoir composed of smaller pores that do not drain under gravity and hold water under capillary action.This representation gives the model computational parsimony.
MOBIDIC is composed of several MATLAB ™ subroutines.Preprocessing of topographic and geomorphologic model inputs, e.g., pit-filling of digital elevation model, determination of flow directions, computation of flow accumulation, and delineation of the river network and the basin boundary, is done in ArcGIS ™ using the Hydrology Toolbox.Other required model inputs are land cover and soil maps, which are in turn used to derive parameters such as albedo, turbulent heat exchange coefficient (neutral), canopy interception capacity, and soil hydraulic properties.The model can output time series of streamflow at any point along the river network, and the hydrologic fluxes (e.g., infiltration, runoff, and ET) and states (e.g., soil temperature and water content of the soil capillary and gravity reservoirs) at any point in the basin.More details about MOBIDIC can be found in Campo et al. (2006) and Castelli et al. (2009).

Mass and energy balance
A schematic diagram of MOBIDIC's mass balance for a typical soil unit (on a hillslope) is shown in Fig. 1 where θ sat , θ fld , and θ res are the volumetric soil moisture [-] at saturation, field capacity, and residual content, respectively.The parameters θ sat , θ fld , and θ res are initialized based on soil texture type and using typical values reported by Rawls et al. (1982).The water holding capacity of the plant reservoir, W p,max , and the surface reservoir, W s,max , are parameterized based on topography, land cover, and land use.
Within each computational time step, dt [T], the hydrologic fluxes [L T −1 ] linking elements across the landscape include precipitation P , infiltration-excess runoff R H , partialarea (saturation from below) runoff R D , return flow R R , total runoff R T , infiltration I , absorption Q as from W g to W c , percolation Q per , lateral subsurface flow Q L , capillary rise Q cap , and evapotranspiration E. These water fluxes can be limited by the available water to be transported, the allowable transport rate, or the available receiving storage.
The water content states of the four reservoirs evolve according to Eqs. ( 4) to (7): The terms on the right-hand side of Eqs. ( 4) to ( 7) are the hydrologic fluxes described earlier.For each soil moisture storage unit, the allowable rate of infiltration I , absorption Q as from W g to W c , percolation Q per , and lateral subsurface flow Q L are formulated according to Eqs. ( 8) to ( 13): where K s is the soil saturated hydraulic conductivity [L T −1 ]; κ, γ and β are rate coefficients [1/T]; and W gu is simply an updated value of W g introduced for conciseness of Eqs. ( 12) and ( 13): The subscripts "up" and "down" denote incoming flow from upstream cell(s) and outgoing flow to downstream cell, respectively.The lateral subsurface flow Q L,down and total surface runoff R T,down to downstream cell are calculated as where φ cha is the channelization parameter [0-1], i.e., the fraction of the flow from upstream cell(s) that is routed directly to the downstream cell.Typically, using the full version of MOBIDIC in distributed catchment modeling, the fluxes from upstream cell(s) are calculated through flow routing.However, since a 1-pixel version of MOBIDIC is used in this research, the fluxes from upstream cell(s) are calculated as follows: which means that the outgoing flows to downstream cell computed at time step t = i, multiplied by a parameter [0-1], become the incoming flow from upstream cell(s) at the next time step, t = i + 1.
Infiltration fills the gravity storage at a rate limited by K s .Absorption flux Q as draws water from gravity storage into available capillary storage.The parameter κ is a linear rate coefficient.The water in gravity storage is lost to percolation or to lateral subsurface flow.Both are again characterized by linear rate coefficients γ and β. κ, γ , and β are dimensionless parameters with values from 0 to 1.For fine soil texture, typically κ is close to 1 since the capillary reservoir is filled first before any substantial filling of the gravity reservoir.Meanwhile, based on a comparison of Eq. ( 12) with the analytic percolation equation of Eagleson (1978), a good initialization of γ is K s /W g,max .
The conceptualization of soil water storage as gravity and capillary storage and the flux relations (see Eqs. 9 to 13) constitute the core of the simplified modeling system.Infiltration fills the larger pores increasing gravity storage.Water is moved from the gravity storage into the smaller capillary storage pores.Losses to the groundwater and lateral flow are only from gravity storage.Simple linear rate constants characterize the time scales of these exchanges.This simple representation is based on physical considerations and they result in a parsimonious and computationally efficient modeling approach.
The soil capillary water storage unit can also receive water from capillary rise from shallow water table.There are a number of available capillary rise models, e.g., Gardner (1958), Eagleson (1978), and Bogaart et al. (2008).They vary primarily based on their parameterization of K s and the soil matric potential ψ [L] as function of soil moisture.The capillary rise model of Salvucci (1993), shown in Eq. ( 19), was chosen because, unlike other models, it allows direct calculation of the capillary rise Q cap [L T −1 ] as a function of ψ and the mean distance of the unsaturated soil layer from the water table d w [L]: where ψ 1 [L] is the bubbling pressure, and n [-] is the product of the Brooks-Corey pore-size distribution index and pore-size disconnectedness index.Brooks and Corey (1964) is used to compute ψ: The effective soil saturation S [-] is computed as The evapotranspiration has three components: E 1 is evaporation from canopy retention, E 2 is evaporation from free surface water surfaces, and E 3 is evapotranspiration from the A. Castillo et al.: Gravitational and capillary soil moisture dynamics soil: Equation ( 25) has the form of an S-curve and was chosen because it mimics the nonlinear behavior of actual ET as a function of potential evapotranspiration (PET) and relative soil saturation S. It uses a single parameter ξ .S is multiplied by 10 for convenience such that ξ takes on non-negative integer values (suggested value: 2 or 3).
Except during a precipitation event and the subsequent draining period, most of the fluxes are inactive.During dry conditions, the only significant fluxes are ET 3 and The PET is determined through surface energy balance under potential (energy-limited) conditions as where ρ w is density of water, L v is the latent heat of vaporization, R n is net incoming radiation, H is sensible heat flux, and G is heat flux into the soil.Upon calculation of actual evaporation through Eqs. ( 22) to ( 25), the energy balance is solved again to update the surface temperature state.The turbulent fluxes are computed according to Eqs. ( 27) and ( 28), where ρ a is the density of air, C a is heat capacity of air, C H is turbulent heat exchange coefficient, and U is wind speed; T a and q a are the temperature and specific humidity of air, respectively; T s and q s are the temperature and specific humidity of the surface (soil and vegetation continuum), respectively.
The unknown surface temperature T s and soil heat flux G are estimated using the heat diffusion equation where ρ s is the density, C s is heat capacity, k is thermal conductivity, T is temperature of soil, and t is time.Equation (29) is integrated forward in time using a parsimonious three-point vertical discretization: where z d and z y are the damping depths of daily and yearly heatwaves, respectively.The lower boundary condition is a constant temperature (T constant ) roughly equal to the annual mean air temperature.The upper boundary condition is

The SHAW model
The SHAW models the transfer of heat, water, and solute within a 1-D vertical profile composed of multi-layered and multi-species plant cover, snow layer, dead plant residue layer, and multi-layered soil.It was first developed by Flerchinger and Saxton (1989)  In SHAW, a soil column is discretized into computational nodes.The fluxes between nodes are solved using implicit finite difference.The required inputs include general site information (e.g., location, elevation, aspect); parameters for soil, snow, and vegetation; meteorological forcings (precipitation, air temperature, total solar radiation, wind speed, and relative humidity); lower boundary conditions; and initial states for soil moisture and temperature.Optional inputs are time series of water sources or sinks and time series of vegetation parameters.The latter, which include canopy height, biomass, leaf diameter, leaf area index, and effective root depth, are specified in this study.

Correspondence between SHAW and MOBIDIC variables
In order to compare the soil moisture dynamics of SHAW and MOBIDIC, the parameters used in both models were set as consistently as possible.For example, the surface albedo is the same in both models.Also, the soil water content at saturation of MOBIDIC and the corresponding depth-averaged value of SHAW are the same.SHAW and MOBIDIC output different state variables.SHAW gives the volumetric soil moisture θ i [-] at each soil node i, while MOBIDIC gives the equivalent water depth W [L] stored as capillary and gravity water for its single soil layer.To allow comparison, the results of the two models were converted to depth-averaged soil moisture θ [−] averaged over MOBIDIC's soil depth d.Note that typically, as done in this study, SHAW's total soil depth is more than the depth of the hydrologically active soil layer.Let the superscripts "O", "S", and "M" denote observed, SHAWsimulated, and MOBIDIC-simulated variables, respectively.For SHAW, θ S is the depth-weighted average of the θ i values: where d is the sum of the thickness of each soil layers, d i [L], For MOBIDIC, θ M is the sum of the equivalent depth [L] of water stored in the capillary reservoir W M c , the gravity reservoir W M g , and the time-invariant residual water content W M r , normalized by d: The soil moisture can also be expressed as equivalent depth: Moreover, in order to compare with MOBIDIC's partitioning of soil moisture into gravity water and capillary-bound water, the total water content simulated by SHAW for the ith soil layer is partitioned into gravity water W S g,i and capillary water W S c,i .Water in excess of the field capacity is considered gravitational storage water, while water between residual water content and field capacity is considered capillary bound.
By summing over the same soil depth d, the corresponding total water stored in the gravity and capillary reservoirs simulated by SHAW are obtained by

Test sites
The comparison is performed using two sites with contrasting climatic regimes.The first site is the "Lucky Hills" catchment in Walnut Gulch experimental watershed, Arizona.The climate is semiarid with two-thirds of the annual  precipitation occurring during the North American monsoon from July to September (Goodrich et al., 2008;USDA-ARS, 2007).The site has a mild topography with deep groundwater table.The vegetation is dominated by shrubs (creosote bush or Larrea tridentata) with sparse grass (USDA-ARS, 2007).The soil is gravelly sand and loam.Meteorological data and measurements of soil moisture and temperature are available from the USDA-ARS Southwest Watershed Research source (http://www.ars.usda.gov/main/site_main.htm?modecode=20-22-10-00).Soil moisture is measured at seven depths (5, 15, 30, 50, 75, 100, and 200 cm).For consistency, the SHAW model is set up with nine soil nodes with the two extra nodes located at 0 and 300 cm.A subset of the calibrated soil parameters of the SHAW model for this site is shown in Table 2.The soil composition (percent sand, silt, and clay) is based on measurements, while the rest is manually calibrated but constrained to be within the typical range of values for given soil texture recommended by Rawls et al. (1982).
The second site is the USDA Soil Climate Analysis Network (SCAN) station "Mayday" in Yazoo, west central Mississippi (32 • 52 N, 90 • 31 W, elevation 33 m a.s.l.).Located on the Mississippi Delta, this site is characterized by thick clayey alluvial soil, flat topography, shallow groundwater table, and agricultural land use.Its humid subtropical climate is significantly influenced by the warm and moist air often originating from the Gulf of Mexico.In contrast to Site 1, precipitation here is almost evenly distributed throughout the year.Hourly meteorological data and measurements of soil moisture and soil temperature are available from the SCAN source (http://www.wcc.nrcs.usda.gov/scan).Soil moisture and temperature are measured at five depths (5, 10, 20, 50, and 100 cm).The SHAW model was setup with eight soil nodes with the three extra nodes located at 0, 75, and 150 cm.

Calibration
The periods simulated for both sites comprise a 1-year warmup period, 2-year calibration period, and 1-year validation period.The use of a warm-up period greatly reduces possible errors that can arise from incorrect model initialization.
First, the SHAW model is calibrated.Next, the albedo and depth-averaged saturated water content of the calibrated SHAW model are copied to the 1-pixel MOBIDIC model.Since for both sites the soil moisture observed at depths greater than 50 cm are quite stable, the soil depth chosen for comparison between SHAW and MOBIDIC is the top 50 cm.This is also the soil depth d used for MOBIDIC.With d and θ M sat fixed, the remaining parameters to be calibrated to set MOBIDIC's W c,max and W g,max are θ M fld and θ M res (recall Eqs. 2 and 3).Once the MOBIDIC model is calibrated, the values of θ M fld and θ M res are used to calculate SHAW's θ S fld,i and θ S res,i for each layer, such that z = 0-50 cm, θ S fld = θ M fld , and θ S res = θ M res .Another set of simulations is performed with the calibrated MOBIDIC model but using d = 30 cm.The results are also compared against observed and SHAW-simulated values averaged over this depth.Table 3 lists the calibrated soil properties for both the SHAW and MOBIDIC models.For the SHAW layers, the calibrated θ sat,i ranges from 0.19 to 0.21.Although low, these values are expected because the site is very gravelly and rocky.Also listed are the calibrated capacities of the dual-pore of MOBIDIC and the corresponding soil water contents at saturation, field capacity, and residual content for the top 50 cm of soil.
To guide the manual calibration, several objective and qualitative checks were performed.The Pearson correlation coefficient R and the absolute value bias B are used as objective measures of goodness of fit.R measures the phase relationship or the match in timing between the modeled and observed values.Its main drawback is that a model which systematically over-or underpredicts the data can still have R close to unity.This drawback is addressed by also computing the absolute bias.
However, the objective of the calibration was not simply to get the best value of the objective metrics.Emphasis was also given to the realism of the model.For instance, param- eters such as θ sat and θ fld were constrained based on literature values.Moreover, the time series of SHAW-simulated soil moisture at various depths were also plotted and visually compared against observations.For MOBIDIC, the hourly time series and annual total of fluxes, e.g., of ET, were qualitatively checked and compared against reported values.

Site 1 -Lucky Hills, Arizona
The soil moisture simulated by SHAW for the Lucky Hills site at various depths is plotted alongside observed values in Fig. 2. The magnitude range and temporal dynamics of θ for all seven nodes are in close agreement especially near the surface.Particularly, SHAW reproduced the sharp difference between the drier and more dynamic top four soil nodes (z = 5, 15, 30, 50 cm) and the wetter and less dynamic bottom three nodes (z = 75, 100, 200 cm).Notice also that during precipitation events, the top four layers become wetter than the deeper layers, a process called "profile inversion".This particular phenomenon cannot be resolved in singlelayer models such as MOBIDIC.
Next, the modeled W S and W M for z = 0 to 50 cm are plotted alongside observed values in Fig. 3.Both SHAW (R = 0.89, B = 0.018) and MOBIDIC (R = 0.88, B = 0.023) accurately reproduced the observations for the 2-year calibration period.More importantly, the performance of MOBIDIC in capturing the magnitude range and temporal dynamics of soil moisture is comparable to that of SHAW. Figure 3 also shows the time series of observed precipitation and the MOBIDIC-simulated ET.High ET occurs around Julian days 200-300, with a maximum of about 5 mm day −1 .For the rest of the year, ET rarely exceeds 0.5 mm day −1 .These are realistic values.
To illustrate the adequacy of the dual-pore soil structure of MOBIDIC, the W c and W g simulated by MOBIDIC are plotted against the corresponding values derived from the outputs of SHAW (see Fig. 4a and b).As shown, the two models are in general agreement indicating that the magnitude range and temporal dynamics of MOBIDIC's W c and W g have correspondence in SHAW.Two plots are used to highlight the difference in time scale between the capillary-bound and gravity water.Gravity storage is filled during rain storms and it is emptied rapidly.In contrast, capillary-bound water has multiday time scale in its dynamics with its recession lasting for months.
Using the SHAW and MOBIDIC models calibrated for the top 50 cm, the performance metrics were also evaluated for z = 0-30 cm.Table 4 summarizes the results for Site 1.The degradation of model performance in the validation period is minimal.Actually, the performance even improved for soil moisture in the validation period for the top 30 cm.

Site 2 -Mayday, Mississippi
In contrast to Site 1, this site is sub-humid.Figure 5 shows the soil moisture simulated by SHAW (lines) vs. observations (points).The soil moisture generally increases and becomes more stable with depth indicating the presence of a shallow water table.The soil node at z = 50 cm remained practically saturated for the entire simulated period.Overall, the SHAW-simulated θ S i at various depths resembles the magnitude range and temporal dynamics of the observations.However, θ S i does not dip down as low as the observations e.g., during days 170-270, probably because during this dry period the effect of plant transpiration through root suction is not correctly captured by the SHAW model of the site.
Figure 6a plots the time series of observed precipitation and MOBIDIC-simulated ET.After precipitation wetting events, the evapotranspiration rate can be as high as about 12 mm day −1 .During the rest of the year, ET is normally 1-3 mm day −1 .
The two objective measures of goodness-of-fit are evaluated using only the equivalent depth of water stored in the top 50 cm of soil.For the 2-year calibration period, SHAW performed well (R = 0.78, B = 0.005) while MOBIDIC performed slightly better (R = 0.86, B = 0.001); see Fig. 6b.For the validation period, both models significantly underestimate θ .As shown in Fig. 6b, the soil column remained saturated during almost the entire validation period, whereas SHAW and MOBIDIC naturally predicted the recession of θ due to ET and drainage.A possible reason for the discrepancy is irrigation in upstream areas, which causes significant lateral subsurface flow, raises the groundwater table, and is not properly accounted for in the two models applied without upstream conditions.
As expected of a site with shallow groundwater table, clayey soil, and sub-humid climate, the soil capillary reservoir remains full during non-drought years, i.e., the soil remains near or above field capacity.The fluctuation of the total soil moisture at this site is associated only with the soil grav-A.Castillo et al.: Gravitational and capillary soil moisture dynamics ity reservoir.Figure 6c shows that the MOBIDIC-simulated W S g and the equivalent values W S c derived from SHAW track one another in both magnitude range and dynamics.Again, this indicates that MOBIDIC's dual-pore soil has behavioral correspondence in the RE-based SHAW model.
The values of the performance metrics for soil moisture and temperature for Site 2 are summarized in Table 5.Similar to the findings in Site 1, the results here show that MO-BIDIC's simple dual-pore storage model captures the essential local-scale soil moisture dynamics that is comparable to those simulated with a solver like SHAW.Furthermore, the two models performed relatively better in Site 1 than in Site 2 because the former is well-represented by an independent vertical soil column, whereas in the latter, lateral subsurface fluxes and groundwater interactions are important.

Summary
The local scale (referring to vertical discretization of the soil column) in distributed hydrologic models is often modeled using grids with millimeters to centimeters spacing.This is required for the stable and correct solution of vertical soil moisture dynamics based on the Richards equation.This local-scale treatment is embedded in distributed models with lateral gridding with a scale of tens to hundreds of meters.The distributed models are applied across entire basins.The desired applications to larger domains and in ensemble mode is limited by (1) the computational demand of the detailed treatment of local-scale processes and (2) the number of model states that need to be initialized.
In this study we compared the effective performances of two distinct approaches to the characterization of the local scale.In the detailed approach, a numerical solver of the Richards equation for the vertical soil moisture dynamics (coupled to heat flow) is used.In the simpler and computationally efficient and parsimonious conceptual approach, a dual-pore characterization of a single soil unit is used.The various hydrologic fluxes act on the two reservoirs in different ways.Also, an exchange flux links the two pore storages.This conceptual approach is based on physical reasoning and is embedded in the MOBIDIC distributed hydrologic model.
The soil moisture state variables simulated by the two models are compared to field observations.The comparisons are made at two sites with contrasting climate (semiarid and sub-humid).The parameters that can be linked between the two models are constrained to be consistent.The calibrated models are then compared with each other and the observations.At each of the two sites, the magnitude range and temporal dynamics of the gravity storage water and the capillary storage water are comparable.This result is the basis for using the simplified local-scale characterization to largedomain and ensemble distributed hydrologic model applications.
Macropore pathways and vertical structure in the soil column that are associated with horizons and parent geology cannot be resolved in the dual-pore conceptual approach.The application of models like MOBIDIC is justified where there is limited or no information on the soil vertical stratification and macropores connectivity.Finally, the role of roots cannot be captured or represented in both detailed and simplified conceptual approaches.Extensive field observations are required before an approach capturing these complications can be designed and implemented.

Figure 1 .
Figure 1.Schematic diagram of MOBIDIC's mass balance at each soil unit.

Figure 4 .Figure 5 .
Figure 4.The SHAW-and MOBIDIC-simulated equivalent depth [cm] of water stored in the soil (a) capillary reservoir and (b) gravity reservoir for the top 50 cm of Site 1.

Figure 6 .
Figure 6.Site 2: (a) observed precipitation [mm day −1 ] and MOBIDIC-simulated ET [mm day −1 ]; (b) observed, SHAW-, and MOBIDICsimulated equivalent depth [mm] of soil water stored in the top 50 cm; and (c) MOBIDIC-and SHAW-simulated equivalent depth [mm] of water in the capillary reservoir of the top 50 cm of soil.

Table 2 .
Rawls et al. (1982)perties of the SHAW model of Site 1. b and ψ e are the Campbell pore-size distribution index and air-entry potential, respectively.zwasfixed; percent sand, silt, and clay is based on observations, while the rest is manually calibrated but constrained by typical values reported byRawls et al. (1982).

Table 3 .
Summary of modeled soil depths and calibrated soil water capacities.D is the total soil depth modeled in SHAW while d is the soil depth modeled in MOBIDIC and the depth used for comparison.

Table 4 .
Performance of the SHAW and MOBIDIC models of Site 1 for the calibration period(years 2007 and 2008)and validation period (year 2009).

Table 5 .
Performance of the SHAW and MOBIDIC models of Site 2 for the calibration period (water years 2006 and 2007) and validation period (water year 2008).