Journal of Hydrology 570 (2019) 598–611

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Drywell infiltration and hydraulic properties in heterogeneous soil profiles
Salini Sasidharan
a
b
c

a,b,⁎

b

a

c

, Scott A. Bradford , Jiří Šimůnek , Stephen R. Kraemer

T

Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USA
United States Department of Agriculture, Agricultural Research Service, U. S. Salinity Laboratory, Riverside, CA 92507, USA
Systems Exposure Division, National Exposure Research Laboratory, Office of Research and Development, US Environmental Protection Agency, Los Angeles, CA, USA

ARTICLE INFOABSTRACTThis manuscript was handled by C. Corradini,
Editor-in-Chief, with the assistance of Philip
Brunner, Associate Editor

Drywells are increasingly used to capture stormwater runoff for surface infiltration and aquifer recharge, butlittle research has examined the role of ubiquitous subsurface heterogeneity in hydraulic properties on drywell
performance. Numerical experiments were therefore conducted using the HYDRUS (2D/3D) software to systematically study the influence of subsurface heterogeneity on drywell infiltration. Subsurface heterogeneity wasdescribed deterministically by defining soil layers or lenses, or by generating stochastic realizations of soil hydraulic properties with selected variance (σ) and horizontal (X) and vertical (Z) correlation lengths. The infiltration rate increased when a high permeability layer/lens was located at the bottom of the drywell, and hadlarger vertical and especially horizontal dimensions. Furthermore, the average cumulative infiltration (I) for 100
stochastic realizations of a given subsurface heterogeneity increased with σ and X, but decreased with Z. This
indicates that the presence of many highly permeable, laterally extending lenses provides a larger surface areafor enhanced infiltration than the presence of isolated, highly permeable lenses. The ability to inversely determine soil hydraulic properties from numerical drywell infiltration results was also investigated. The hydraulicproperties and the lateral extension of a highly permeable lens could be accurately determined for certain
idealized situations (e.g., simple layered profiles) using constant head tests. However, variability in soil hydraulic properties could not be accurately determined for systems that exhibited more realistic stochastic heterogeneity. In this case, the heterogeneous profile could be replaced with an equivalent homogeneous profileand values of an effective isotropic saturated conductivity (K s ) and the shape parameter in the soil water retention function ( ) could be inversely determined. The average value of K s for 100 stochastic realizations
showed a similar dependency to I on σ, X, and Z. Whereas, the average value of had large confidence interval
for soil heterogeneity parameters and played a secondary role in drywell infiltration. This research provides
valuable insight on the selection of site, design, installation, and long-term performance of a drywell.

Keywords:
Drywell
Infiltration
Stochastic simulation
Inverse optimization
Effective hydraulic propertiesHYDRUS (2D/3D)

1. Introduction
Freshwater resources determine the balance between food security,
public health, industrial growth, economic development, environmental sustainability, and the survival of ecosystems across the world(Rosegrant et al., 2002). However, water scarcity occurs in many partsof the world. About 1.1 billion people worldwide lack access to safe
water and at least two-thirds of the global population, over 4.0 billion
people, live under conditions of severe water scarcity at least 1 montheach year (Kummu et al., 2016). Groundwater decline, low or nonexistent river flow, and water pollution are some of the major signs ofwater scarcity (Postel, 2000). Although the global hydrologic cycle
produces more freshwater than is needed to sustain the current global
population, > 69% of this resource is not available to humans due to⁎rapid runoff or distribution in geographically isolated regions (Postel,2000). Growing food to meet the needs of the global population requires extraction of a large volume of water and irrigated agriculture
accounts for two-thirds of all water withdrawals from natural water
resources (Shiklomanov, 2000). The unsustainable pumping ofgroundwater exceeds the renewable water limits in many regions(Postel, 2000; Postel et al., 1996). Therefore, there is an urgency to findalternative sources of freshwater, to invest in water infrastructure, todevelop innovative technologies to recycle, reuse, and preserve stormwater, and to improve water management.
Several countries around the world are practicing artificial ormanaged aquifer recharge using engineered systems to collect surface
water and enhance its infiltration and subsequent movement to aquifers
to augment groundwater resources. Infiltration and managed aquifer

Corresponding author at: USDA, ARS, Salinity Laboratory, Riverside, CA 92507, United States.
E-mail address: salinis@ucr.edu (S. Sasidharan).

https://doi.org/10.1016/j.jhydrol.2018.12.073
Received 20 September 2018; Received in revised form 21 December 2018; Accepted 26 December 2018
Available online 11 January 20190022-1694/ © 2019 Elsevier B.V. All rights reserved.

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

recharge can be achieved using surface infiltration (e.g., infiltration
basins), direct injection (e.g., aquifer storage, transfer, and recovery),and vadose zone infiltration (e.g., drywells) (Bouwer, 2002; Bouwer
et al., 1999; Dillon, 2005; Edwards et al., 2016; Sasidharan et al., 2018)techniques. In recent years, drywells have gained a lot of attention in
the United States and across the world. Drywells are vadose zone infiltration wells that end (at least 1–3 m) above the water table(LACDPW, 2014). There are several potential advantages associatedwith vadose zone infiltration over surface infiltration and direct injection wells. This includes a minimal evaporative loss of water, a small
installation area, a large ponding depth, and removal of contaminantsin the vadose zone (Edwards et al., 2016; Sasidharan et al., 2018).
There are currently > 75,000 drywells in the states of Washington,
Oregon, California, and Arizona of the United States (Edwards et al.,2016) and their number is increasing annually. However, not manystudies have evaluated the long-term performance of drywells and theirenvironmental impact.
The purpose of a drywell is to transfer water from a poorly drained
or low permeable surface capture area to an aquifer for subsurface
water storage. The current installation practice involves drilling adrywell at a site and inspecting soil borings or profiles for the presenceof layers with adequate permeability to ensure drainage for design
storm events (Watt and Marsalek, 2013). However, the lateral extension
and vertical continuity of these permeable layers are not determined
during well installation. It is logical to anticipate that the position,
lateral extension, thickness, and spatial variability of a permeable
layer/lens adjacent to a drywell will play a critical role in determining
the infiltration behavior. Indeed, an understanding of the soil hydraulicproperties has been reported to be a very important factor for the
successful design, installation, execution, and long-term performance ofa drywell (Edwards et al., 2016). In our previous study, numerical andfield-scale experiments were conducted to characterize the drywell infiltration behavior. HYDRUS (2D/3D) computer software was modified
to simulate the transient boundary conditions of the Maxwell Type IV
drywell geometry. Effective soil hydraulic properties of two drywellsites across California were estimated via inverse optimization of the
field-scale falling head data (Sasidharan et al., 2018). However, no
studies have systematically investigated the role of subsurface heterogeneity on drywell performance.
Previous studies have attempted to estimate the effective saturated
hydraulic conductivity of an unsaturated soil using a constant head
borehole test. For example, analytical solutions have been developed todetermine the hydraulic properties of a homogenous soil domain using
constant head data and numerical simulations have been conducted toevaluate the accuracy of these analytical solutions (Reynolds and Elrick,
1985; Xiang et al., 1997). However, estimates of the hydraulic conductivity were often unrealistic due to the assumption of soil homogeneity in these analytical solutions (Reynolds and Elrick, 1985) andinherent soil heterogeneity. Later, Xiang et al. (1997) proposed ananalytical approach to analyze the flow contribution of each layer tothe total flow and numerically evaluated the approach and its assumptions. However, this test neglected unsaturated flow and was only
suitable for determining the hydraulic conductivity of layered media inthick unsaturated zones. Drywells are used not only to recharge deep
groundwater aquifers, but also as a flood control system in paved urban
settings where a shallow aquifer may be present. Alternatively, numerical solutions that account for both unsaturated and saturated flow
can be used to inversely estimate the unsaturated soil hydraulic parameters of thin as well as thick soil layers/lenses/fractures that areneeded to accurately evaluate the performance of a drywell.
The above approaches considered idealized homogeneous orlayered soil systems. However, most field soils are highly heterogeneous, and the hydraulic conductivity may change over a short distance. The heterogeneity of a field site can be estimated using ground
topography, soil logs, and other geophysical techniques. However,these methods are very expensive and time-consuming and it is unlikelythat such detailed information can be obtained for all drywell sites.
Alternatively, numerical solutions and stochastic approaches, such as
the scaling methodology, can be used to study drywell infiltration. The
scaling approach was first introduced by Miller and Miller (1956) who
derived scaling relationships between the pressure head and the hydraulic conductivity using the microscopic laws for capillary pressure
forces and viscous flow based on the similarity of the pore space(Vereecken et al., 2007). The Miller and Miller (1956) similitude approach has been widely used in the past to analyze flow and transport
processes in heterogeneous unsaturated soil systems (Hammel andRoth, 1998; Roth, 1995; Roth and Hammel, 1996; Tseng and Jury,
1994). In the last few decades, several other approaches such as stochastic perturbation methods, the scaleway approach, the stream-tube
approach, the aggregation concept, inverse modeling approaches, anddata fusion approaches were developed for upscaling of soil water
processes (Vereecken et al., 2007). To the best of our knowledge, nopublished studies have investigated the influence of stochastic vadosezone heterogeneity on the drywell infiltration behavior and the determination of effective unsaturated soil hydraulic properties.
The objective of this study was to investigate the influence of subsurface heterogeneity on drywell infiltration. The HYDRUS (2D/3D)
software was used to directly simulate cumulative infiltration volumes
for selected drywell geometries and soil heterogeneities under constantor falling head conditions. Subsurface heterogeneity was described inthis model deterministically by defining soil layers or lenses, or by
generating stochastic realizations of soil hydraulic properties with selected variance and correlation lengths. The numerically generated data
were then used in inverse optimizations to determine the hydraulicproperties and the lateral extension of individual layers or lenses, or todetermine soil hydraulic properties of an equivalent homogeneous
profile. The influence of stochastic subsurface heterogeneity parameters(e.g., the variance, horizontal, and vertical correlation length) on cumulative drywell infiltration and equivalent homogeneous profile values of the saturated isotropic hydraulic conductivity (Ks ) and the soil
water retention function (van Genuchten, 1980) shape factor ( ) were
determined.
2. Materials and method
2.1. Mathematical model
The HYDRUS (2D/3D) software package (Šimůnek et al., 2018;Šimůnek et al., 2016) was used to simulate drywell infiltration andredistribution in the vadose zone, by numerically solving the 2-dimensional axisymmetric form of Richards equation. The unsaturated
soil hydraulic functions were those of van Genuchten (1980) and
Mualem (1976). Unless otherwise noted, default values of soil hydraulicparameters (e.g., the saturated soil hydraulic conductivity, Ks ; theshape parameter in the soil water retention function, ; the residual soil
water content, r ; the saturated soil water content, s ; the pore-size
distribution parameter in the soil water retention function, n; and thetortuosity parameter in the hydraulic conductivity function, l) weretaken from the HYDRUS (2D/3D) Soil Catalog (Šimůnek et al., 2011)(Table 1). The van Genuchten hydraulic parameters in the soil catalogof the HYDRUS (2D/3D) software for selected soil textures were takenfrom (Carsel and Parrish, 1988).
Fig. 1 shows a schematic of an example simulation flow domain,
initial conditions, and boundary conditions. This hypothetical 2-dimensional axisymmetrical domain was adapted from the engineering
design and the soil log information collected from the Fort Irwin drywell site presented in our previous study (Sasidharan et al., 2018). The
initial condition was specified in terms of the soil water pressure head h(x, z) and was set to a constant pressure head of −0.5 m for the entire
flow domain. The nodes at the bottom boundary were assigned a freedrainage boundary condition. A no-flux boundary condition was assigned to the upper, right, and lower left sides of the flow domain
599

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

dimensional triangular finite element mesh using the MESHGEN tool
available within HYDRUS (2D/3D) (Šejna et al., 2014). The mesh wasrefined at the left part of the domain where infiltration from the drywell
was simulated. To reduce the mass balance error, the finite element
mesh was adjusted such that the size of elements was smaller (0.05 m)
near the drywell and the grid size was gradually increased with theradial distance from the drywell, with a maximum element size of
0.75 m. In numerical experiments with multiple soil layers (Fig. 8), thefinite element mesh was further refined (finite elements between 0.05
and 0.25 m near the drywell) with a maximum element size of 0.75 m toobtain a smaller mass balance error.

Table 1
The residual soil water content ( r ), the saturated soil water content ( s ), theshape parameter ( ), the pore-size distribution parameter (n), the saturated
isotropic hydraulic conductivity (K s ) and the tortuosity parameter (l) for different textural classes obtained from the HYDRUS (2D/3D) Soil catalog(Šimůnek et al., 2011).
Soil

r [−]

s [−]

Sandy clay loam (SCL)
Sandy clay (SC)Sand (S)Sandy loam (SL)Silt loam (SiL)

0.1
0.1
0.045
0.065
0.067

0.39
0.38
0.43
0.41
0.45

[m−1]

n [−]

Ks [m min−1]

l [−]

5.9
2.7
14.5
7.5
2.0

1.48
1.23
2.68
1.89
1.41

2.18 × 10−4
2.0 × 10−5
4.95 × 10−3
7.37 × 10−4
7.5 × 10−5

0.5
0.5
0.5
0.5
0.5

2.2. Numerical experiments(Fig. 1). The complex geometry and the dynamics of water flow in andout of the drywell were implemented as a Reservoir Boundary Condition which is described in detail by Šimůnek et al. (2018) and
Sasidharan et al. (2018). The Reservoir Boundary Condition was specified along the left and bottom boundary of the drywell (Fig. 1). Somesimulations employed a constant head boundary condition for 200 min
at the surface of the drywell or at textural interfaces. Others consideredthat the drywell was initially filled to its surface, and then allowed forthe water level to fall over time to a selected depth as a result of infiltration into the soil profile.
The size of the hypothetical simulation domain was adjusted depending on the considered soil heterogeneity. Simulations that examined the location, lateral extension, thickness of a sandy lens systems(Fig. 2) employed a 3 m wide and 35 m deep domain, the layered
system simulations (Figs. 6–8) employed a 5 m wide and 35 m deep
domain and 5 m and 20 m wide domains for 1000 and 5000 min simulation (Fig. 8D and E), simulations that examined lateral extensionof lenses (Fig. 9) considered a 15 m wide and 50 m deep domain, andstochastic simulations (Figs. 3–5, 10–13) considered a 12 m wide by
40 m deep domain. The simulation domain was discretized into a two-

2.2.1. Deterministic heterogeneity – Highly permeable lensNumerical experiments were initially conducted to understand the effects of the position (Z = −8, −10.5, −20.5, and −29 m), the lateral extension (width (W) = 0.09, 0.39, 0.89, and 1.9 m and thickness (H) of
0.91 m), and the thickness (H = 0.91, 1.91, 3.91, 5.91 m and W = 0.09 m)of a highly permeable lens (i.e., a sandy lens surrounded by a sandy clay
loam soil) on drywell infiltration behavior for 1481 min. A homogeneous(sandy clay loam) flow domain (W = 3 m × H = 35 m) with an inclusionof a sandy lens was considered in these simulations. The two-dimensional(2D) axisymmetrical flow domain, drywell geometry, initial values of thepressure head at the beginning of the simulation, time-variable boundary
conditions, and observed (OBS) data (Fig. 2) were adapted from the Fort
Irwin drywell falling head experiment presented in Sasidharan et al.,(2018). In brief, a fire hydrant was used to flood the detention pond at theFort Irwin drywell. A pressure transducer connected to a datalogger was
lowered into the drywell and the water depth in the drywell was measured
every minute for 1481 min. The in-situ soil hydraulic properties (Ks and )were determined by the inverse parameter optimization of the observed
falling head data for an equivalent sandy clay loam soil profile using

Fig. 1. The geometry, the boundary conditions, and the soil profile for the heterogeneous (sandy clay loam and sand) 2D-axisymmetrical flow domain (with a highly
permeable sandy lens) (A) the corresponding results are presented in Fig. 2. The geometry and the boundary conditions for the sandy loam soil profile used instochastic simulations (B) the corresponding results are presented in Figs. 3, 5, 10, and 11. The corresponding drywell geometry is presented in our previous study(Sasidharan et al., 2018).
600

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

Fig. 2. Flow domains, pressure head profiles (observed data from the falling head experiment at the Fort Irwin drywell (OBS – Black solid line) and simulated data
(Model – Red solid line)), and temporal changes in the water level in the drywell (HDW ) for continuous sandy lenses at different depths (Z = −8, −10.5, −20.5, and−29 m) (A), for sandy lenses with different lateral extensions (W = 0.09, 0.39, 0.89, and 1.9 m), and (B), for sandy lenses with different vertical extensions(thickness) (H = 0.91, 1.91, 3.91, 5.91 m) (C). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of thisarticle.)HYDRUS (2D/3D) software. The fitted values of Ks and for sandy clay
loam were 2.25 × 10−6 m min−1 and 2.63 m−1, respectively, and were
employed in this study to be consistent with the drywell at Fort Irwin, CA(Sasidharan et al., 2018).
Additional drywell infiltration experiments (W = 15 m × H = 50 m)
were conducted to see if the lateral extension (W = 0, 2.5, 5, 7.5, or 10 m
and H = 0.5 m) of the sandy lens in a sandy clay domain could be accurately identified. In this case, a constant head was maintained at the surfaceof the drywell, and the cumulative infiltration volume (I) and infiltration
flux (i) were collected for 5000 min.

HYDRUS (2D/3D) Soil Catalog (Ks HYDRUS and HYDRUS, respectively,
Table 2). The flow domain was kept the same in both direct and indirect
simulations, while the hydraulic parameters for non-fitted layers were
kept constant. In falling head experiments, the final pressure headprofile after one falling head experiment (half-full well, Fig. 6) was
imported as the initial condition for the next simulation (full-well,Fig. 6), and the parameter estimation procedure was repeated for the
next layer (Fig. 6). On the other hand, for the constant head simulations, the initial condition was always kept as a constant pressure headof −0.5 m for the entire flow domain for both direct and inverse simulations (Figs. 7 and 8).

2.2.2. Deterministic heterogeneity – Layered soil profile
Numerical experiments were then conducted to determine the effective hydraulic properties of a layered soil profile. A simulation domain of W = 5 m × H = 35 m was used and divided into two or moresoil layers. Various layers of the flow domain were assigned different
soil hydraulic properties as described in Table 2. Sequential numericalexperiments were then conducted to determine infiltration in thelayered systems, and to inversely estimate the hydraulic properties fromsoil layers at the bottom to the top of the drywell using falling (from thedrywell inlet to the textural interface) or constant head (200 min at theinterface of a soil layer) conditions. The cumulative infiltration volumes(I) from the HYDRUS (2D/3D) direct simulations were determined overtime. Values of Ks and for the bottom soil type were then inversely
estimated by fitting to this I data using the parameter estimation procedure in HYDRUS (2D/3D) that is based on the Levenberg-Marquardt
nonlinear parameter optimization method (Marquardt, 1963). Table 2
provides a summary of numerical experiments carried out for various
two- and multi-layered soil systems, including: the sequence of different
soil layers and soil hydraulic parameters used in direct simulations.
Direct simulations employed Ks and values for each soil type from the

2.2.3. Stochastic heterogeneity
The HYDRUS (2D/3D) computer software has an option to generatestochastic (random) distributions of hydraulic conductivity and pressure head scaling factors (αK and αh, respectively) using the MillerMiller similitude approximation (Miller and Miller, 1956). Using thisoption, a specific Ks value can be assigned to the flow domain (or a partof it) that is then multiplied in each node by a scaling factor αK. The
Miller-Miller scaling procedure requires that αK = αh−2. The stochastic
distribution requires three input parameters: the standard deviation ( )of log10(αK), and its correlation length in the lateral (the X-correlation
length) and vertical (the Z-correlation length) directions. The value of
determines the extent of variations in the scaling factors, with higher
values leading to higher variations in scaling factors. On the other hand,the correlation length is a measure of the distance in a specific directionthat scaling factors are related. A high value of the X-correlation lengthmeans that the scaling factor maintains similar values for a greater
horizontal distance. The same applies for the Z-correlation length in thevertical direction (Rassam et al., 2003). In this study, the following
combinations of the scaling factors with variable ( = 0, 0.25, 0.5,
601

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

Fig. 3. Heterogeneous sandy loam profiles for the hydraulic conductivity scaling factor (SF-K s ) generated for = 0, 0.25, 0.50, and 1, X = 1 m, and Z = 0.1 m (A).
The corresponding pressure head (h) (B) and soil water content ( ) (C) profiles at 100 min for a constant head experiment.

1;X = 1 m; Z = 0.1 m), a variable X-correlation length ( = 1, with
X = 0.1, 1, 10 and Z = 0.1; and = 1, with X = 1, 10 andZ = 2), anda variable Z-correlation length ( = 1, with X = 1 and Z = 0.1, 1, 2)
were used. Sandy loam was used in the flow domain(W = 12 m × H = 40 m). Forward simulations were conducted for
100 min by maintaining a constant head in the drywell. The I output
was collected from each simulation and a total of 100 simulations weredone for each combination of stochastic parameters. The mean cumulative infiltration volume (μI), 95% Confidence Interval-Upper Limit(95% CI-UL), and 95% Confidence Interval-Lower Limit (95% CI-LL)
were calculated.
Accurate determination of soil hydraulic properties of a field site is
the key to estimating infiltration from a drywell. A numerical modelthat can determine the hydraulic properties based on the stratigraphic
information of a site can help to identify the best location for the drywell installation. Therefore, the possibility of developing a correlation
between the type of the stochastic parameter combination used in theforward simulation, the resultant cumulative infiltration volume, andthe effective hydraulic properties of the stochastic field was investigated. Inverse optimization of the I output data obtained fromeach direct simulation (a total of 100 simulations per a stochastic
parameter combination) was conducted. Only the hydraulic parameters
Ks and were optimized while the remaining hydraulic parameters ( r ,s , n, and l) were kept constant at default values from the HYDRUS (2D/
3D) Soil Catalog (Table 1) for sandy loam. IBM Cloud Computing Virtual Servers (IBM, USA) with four computing instances with48 × 2.00 GHz core and 192 GB RAM were used to perform 800 inverse
simulations. Python scripts were used for importing and exporting of

input and output data, and for further data analysis. Mean values of Ks ,, 95% CI-UL, and 95% CI-LL were calculated from a total of > 80
inverse simulations.
3. Results and discussion
3.1. A direct (forward) numerical simulations
3.1.1. Deterministic heterogeneity – Highly permeable lens
The current drywell installation practice involves drilling at a potential drywell site until a permeable layer is detected. The permeability of this layer is assumed to infiltrate water efficiently and guarantee the successful performance of the drywell over a specified time
duration. Numerical experiments were conducted to investigate thevalidity of this assumption by simulating drywell infiltration behavior
for different high permeability lens locations, lateral extensions, andthicknesses. Fig. 2 presents a comparison of the observed data obtained
from the Fort Irwin drywell field scale falling head experiment(Sasidharan et al., 2018) and simulated falling heads for heterogeneous(lenses) domains. Fig. 2A shows the influence of the vertical position of
a 50 cm thick continuous sandy lens on the simulated matrix potential(h ) profile in the soil adjacent to the drywell after 1481 min, andchanges in the drywell head (HDW ) over time. The drywell infiltrates
water and HDW falls at a faster rate when the sandy lens occurs at a
deeper vertical position because of a larger pressure head at the lensboundary. Fig. 2B and C show the effect of a sandy lens’s lateral extension and thickness, respectively, on changes in simulated soil matricpotential and the drywell HDW over time. Increasing the lateral602

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

drywell can infiltrate water at a much faster rate than the same sizedlens at the top of a drywell. Furthermore, the cumulative infiltration
volume from a drywell increases with an increasing lateral extension
and a thickness of a highly permeable lens. However, a thin lateral lensat the bottom of a drywell infiltrated a larger amount of water than anisolated highly permeable lens (a thick, vertically extended lens) due toits greater infiltration surface area. These simulation results clearly
indicate that the drywell infiltration performance is sensitive to thehigh permeability lens location, lateral extension, and thickness. Highpermeability lens that occurs deeper in the soil profile and has greater
lateral and vertical extensions are preferred for infiltration.
3.1.2. Stochastic heterogeneity
The above simulations considered highly idealized scenarios,
whereas the natural subsurface is inherently highly heterogeneous.
Stochastic simulations were therefore conducted to better understandthe effect of random heterogeneities represented by stochastic fields ofscaling factors on the cumulative infiltration volume from a drywell.
Fig. 3A shows a sandy loam soil profile with increasing heterogeneityfrom left to right (i.e., = 0, 0.25, 0.5, 1; X = 1 m; Z = 0.1 m). The
scaling factor for Ks (SF-Ks ) shows 8 orders of magnitude variations(0.0001–10000) for = 1, which indicates a highly heterogeneous soil
profile. Fig. 3B and C show the corresponding h and profiles, respectively, at 100 min of simulation for a constant head experiment (afull reservoir). Similar information is presented for the homogeneous
domain ( = 0) as a comparison for results from the heterogeneous
domains. The homogeneous wetting front reached a distance of 2.8 m inboth horizontal and vertical directions from the bottom of the drywell.
In contrast, the wetting front in the heterogeneous profiles extended
laterally 3.0, 4.0, and 5.6 m and vertically 2.6, 2.7, and 5.6 m when thevalue of increased to 0.25, 0.50, and 1, respectively.
The same stochastic parameters will produce different realizationsof heterogeneity and drywell infiltration behavior. Multiple stochastic
realizations were therefore conducted to determine the mean

Fig. 4. The mean cumulative infiltration volume ( µI ) (A) and the mean standard deviation of the cumulative infiltration volume ( I ) (B) as a function of thetotal number of simulations for all considered combinations of stochastic
parameters investigated in this study for a sandy loam flow domain (W = 12,
H = 40) for 100 min.

extension and thickness of the sandy lens also increases the water infiltration and the rate of decrease in HDW . The numerical results demonstrated that a highly permeable lens/layer at the bottom of a

Fig. 5. The mean cumulative infiltration volume ( µI ) (solid line) and the 95% Confidence Interval Upper (95%-CI-UL) and Lower (95%-CI-LL) Limits (dotted lines), asa function of ( = 0.25, 0.5, 1) for X = 1 and Z = 0.1 (A), as a function of X ( X = 0.1, 1, 10) for = 1 and Z = 0.1 (B), as a function of X ( X = 1, 10) for = 1
and Z = 2 (C), as a function of Z (Z = 0.1, 1, 2) for = 1 and X = 1 (D) for direct stochastic simulations for a sandy loam soil. The HomogeneousI (red broken line)
value for the homogeneous sandy loam flow domain with K s and from the HYDRUS (2D/3D) soil catalog (Table 1). (For interpretation of the references to color inthis figure legend, the reader is referred to the web version of this article.)603

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

Fig. 6. The two-layered flow domain with sandy clay above sandy loam (SLB + SCT) (A) and sandy loam above sandy clay (SCB + SLT) (D). The pressure head (h )profiles at 200 min for a half-full well (B, E) and (C, F) 400 min for a full well for falling head experiments. Note that the Reservoir 2 was initially half-full and thewater level was then raised at 200 min to the top of Reservoir 2.

cumulative infiltration volume ( µI ) and the mean standard deviation of
the cumulative infiltration volume ( I ). Fig. 4A and B show µI and I ,respectively, as a function of the total number of simulation realizations
for all different combinations of stochastic parameters. The values of µI
and I stabilize when the total number of simulations was > 75. To be
extra conservative, 100 realizations of stochastic heterogeneity parameters were conducted to determine values of µI and I in this study. Itshould be mentioned that the values of µI reflect the net effects of theposition, lateral extension, and thickness of high permeability lenses.
The value of I, µI , and its 95% confidence interval as a function of(=0.25, 0.5, and 1) when X = 1 and Z = 0.1 m is shown in Fig. 5A. For
these stochastic parameters, µI increased with and exhibited narrow

95% confidence intervals. An increase in increases the variability in
the Ks scaling factor and leads to the formation of more permeable
layers/regions in the soil profile and, consequently, enhanced infiltration from a drywell.
Multiple realizations of drywell infiltration were also conducted forother stochastic parameter combinations. Fig. 5 also presents the calculated values of µI and its 95% confidence interval as a function ofvarious X and Z correlation lengths when = 1. Fig. 5B indicates thatvalues of µI increased with X and showed narrow 95% confidence intervals when Z = 0.1. The value of µI also increased with X whenZ = 2(Fig. 5C), but this increase was more rapid, and the 95% confidence
intervals broadened with X in comparison to Fig. 5B. This increase in µI
604

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

Fig. 7. The two-layered flow domain with sandy clay above sandy loam (SLB + SCT) (A, B) and sandy loam above sandy clay (SCB + SLT) (E, G). The pressure head(h ) profiles at 200 min for constant head BC with a half-full (B, D) and full Reservoir 2 (C, H). The constant head boundary condition is highlighted with a red line.
(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)with X is due to the formation of highly permeable lateral lenses andthat enhances the rapid movement of the wetting front in the horizontal
direction, and facilities infiltration of water over a larger area. The
infiltration area is greater when Z was smaller, and this produced larger
values of µI (Fig. 5B in comparison to Fig. 5C). Similarly, Fig. 5D showsthat µI decreased and exhibited very narrow 95% confidence intervals
when Z increased. This decrease in µI with Z is attributed to the formation of highly permeable isolated lenses and ponded regions withlimited infiltration area. The cumulative infiltration volume averaged
from 100 stochastic realizations increased with increasing and a lateral correlation length, but decreased with an increasing vertical correlation length. Collectively these results indicate that the lateral extent
of highly permeable thin layers with a great vertical connectivity enhances rapid infiltration of water from a drywell to the surroundingsoil.

3.2. Indirect (Inverse) numerical simulations
Drywell infiltration experiments provide an opportunity to inversely
determine equivalent homogeneous hydraulic properties for heterogeneous soil profiles (Sasidharan et al., 2018). Below we numericallyinvestigate the ability to determine soil hydraulic properties of individual layers and lenses, as well as the potential limitations of such an
approach in more complex, heterogeneous systems.
3.2.1. Deterministic heterogeneity – Layered soil profile
Table 2 provides a summary of numerical experiments carried out
for various two- and multi-layered soil systems, hydraulic parameters,
and the statistical parameters for the inverse simulations. Initial simulations considered the effect of a simple two-layered system on drywell
infiltration, with a sandy clay layer on the top (SCT) and a sandy loam
layer on the bottom (SLB), using falling head (from the surface of thedrywell to the textural interface) or constant head (water level at the605

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

Fig. 8. The flow domain with seven different soil-type layers and boundary conditions (A). The pressure head (h) (B, D) and water content ( ) (C, E) profiles for aconstant head boundary (Drywell Reservoir) experiment at 200 min (B, C) and at 1000 and 5000 min (flow domain W = 20 m H = 35 m) (D, E).

textural interface) conditions in the drywell. Similar simulations werealso conducted with the soil layers reversed; e.g., sandy loam layer on
top (SLT) and sandy clay layer on the bottom (SCB). Table 2 indicatesthat a very good agreement was obtained between the direct and inversely simulated I data with a coefficient of determination(R2) > 0.99. Furthermore, the fitted values of Ks (Ks fit) and Ks HYDRUS
were always consistent. Similarly, the fitted values of ( fit) andHYDRUS showed good agreement for constant head experiments, but fit
for SL was underpredicted during the SCT + SLB experiment with thefalling head drywell condition. An explanation for this behavior can beobtained from the simulated profiles of h for these layered simulations(Figs. 6 and 7). A distinct wetting front occurs in both SC and SL layers,
as well as between them, during the constant head experiments due tothe continuous flow of water (Fig. 7). In contrast, only a relatively small
volume of water was available for infiltration during the falling headexperiments and this produced a gradual change in the pressure headsof both soil layers, especially when the highly permeable SL was below
a less conductive SC layer (Fig. 6).
606

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

Fig. 9. The cumulative infiltration volume (I ) (A) and the infiltration flux (i ) (B) for the sandy clay domain with a sandy lens with lateral extensions (X) of 0, 2.5, 5,
7.5, and 10 m and the vertical extension (H) of 0.5 m. Water content profiles for the sandy clay domain with a 0.5 × 10 m sandy lens at 100, 500, 800, and 2000 min
for a constant head boundary (Drywell Reservoir) experiment (C).

Fig. 10. The mean fitted values of the saturated isotropic hydraulic conductivity(Ks ) (solid line) and the 95% Confidence Interval Upper (95%-CI-UL) and Lower(95%-CI-LL) Limits (dotted lines), as a function of ( = 0.25, 0.5, 1) for X = 1 and Z = 0.1 (A), as a function of X ( X = 0.1, 1, 10) for = 1 andZ = 0.1 (B), as afunction of X ( X = 1, 10) for = 1 and Z = 2 (C), as a function of Z (Z = 0.1, 1, 2) for = 1 and X = 1 (D) for inverse stochastic simulations for a sandy loam soil.
The HomogeneousKs (red broken line) value for sandy loam soil from the HYDRUS (2D/3D) soil catalog (7.37 × 10−4 m min−1) (Table 1). (For interpretation of thereferences to color in this figure legend, the reader is referred to the web version of this article.)607

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

Fig. 11. The mean fitted values of the shape parameter ( ) (solid line) and the 95% Confidence Interval Upper (95%-CI-UL) and Lower (95%-CI-LL) Limits (dottedlines), as a function of ( = 0.25, 0.5, 1) for X = 1 and Z = 0.1 (A), as a function of X ( X = 0.1, 1, 10) for = 1 and Z = 0.1 (B), as a function of X ( X = 1, 10) for= 1 and Z = 2 (C), as a function of Z (Z = 0.1, 1, 2) for = 1 and X = 1 (D) for inverse stochastic simulations for a sandy loam soil. The Homogeneous (redbroken line) for sandy loam from the HYDRUS (2D/3D) soil catalog (7.5 m−1) (Table 1). (For interpretation of the references to color in this figure legend, the reader
is referred to the web version of this article.)

Additional numerical experiments were conducted to see if hydraulic properties for a soil profile consisting of 7 soil layers could beinversely determined from constant head infiltration studies at the topof each layer. From top to bottom the soil layers are SiL, S, SC, SCL, SL,
SiL, and SC (Fig. 8A); where S, SCL, and SiL denote sand, sandy clay

loam, and silty loam, respectively. Similar to the two-layered experiments, I data was well described by the model (R2 > 0.99) and valuesof Ks fit and fit were consistent with Ks HYDRUS and HYDRUS, respectively.
A very distinct pattern in h and profiles was observed in Fig. 8B and C,
respectively, which reflects the Ks value of each layer (Table 2). For

Fig. 12. A comparison of mean fitted values of Ks (A) and α (B) obtained by > 80 inverse optimization of the cumulative infiltration volume collected from eightstochastic realizations of different parameter combinations presented in this study. Also shown are the default Ks and α values of sandy loam from the HYDRUS (2D/
3D) soil catalog used in the homogeneous domain.
608

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

has a complex influence on the drywell infiltration capacity which
depends on its location and the hydraulic properties of the neighboringlayers.
Several field-scale approaches have been developed to determine
soil hydraulic properties of boreholes. The atmospheric pressure wave
test measures the bulk air permeability of the vadose zone sediments asa function of time lag and subsurface pressure attenuation between thesurface and subsurface (Weeks, 1978). However, this method only
provides information about the vertical permeability rather than thehorizontal permeability and is not suitable for soil with a high moisture
content or near-surface soil layers (Weeks, 1978). Similarly, gamma
and neutron logs of boreholes only provide information on the vertical
distribution of soil stratigraphy and water-filled porosity (Hodges et al.,2012). Xiang et al. (1997) presented a multistep (constant head) permeameter test to determine Ks value for a layered soil using the
LAYERK model. LAYERK is an analytic code that was designed to determine the field-saturated hydraulic conductivity and the matric flux
potential for individual layers using input parameters including layer
thickness, the water level in the borehole, the radius of the borehole,and flow rate. However, this approach is only suitable for estimating
the Ks value for thick unsaturated zones where the ratio of the water
depth (H ) and the borehole radius (a ) was H ≥ 20 (Jain et al., 2006;

Fig. 13. The cumulative infiltration volume (I ) for a homogeneous sandy loam
flow domain with K s and values obtained from the HYDRUS (2D/3D) soil
catalog (K s = 7.36 × 10−4 m min−1 and = 7.5 m−1) (black broken lineHomogeneous), and for a stochastic domain ( X = 10; Z = 0.1; = 1) withinversely optimized parameter values of average K s = 3.94 × 10−3 m min−1
and = 6.17 m−1 (blue solid line-Stochastic), Upper-Confidence-Interval (UCI)values of K s = 3.68 × 10−3 m min−1 and = 5.40 m−1 (magenta broken lineUCI),
andLower-Confidence-Interval(LCI)
valuesofK s = 4.20 × 10−3 m min−1 and = 6.94 m−1 (magenta broken line-UCI). The
K s and values were presented in Fig. 12. (For interpretation of the references
to color in this figure legend, the reader is referred to the web version of thisarticle.)a

Xiang and Chen, 1993; Xiang et al., 1997). In contrast, a drywell’s H
avalue can range from approximately 10 to 60. Also, in the absence of an
independent determination of α, there are some doubts about thepracticality of using borehole permeameter test to measure Ks (Archer
et al., 2014; Philip, 1985). Moreover, drilling of an additional borehole
to perform the above techniques might not be practical at every drywellsite because they are expensive, energy intensive, time-consuming, require a dry drilling method (a difficult technique), and need vast expertise. On the other hand, the presented numerical modeling approach
using HYDRUS (2D/3D) software can be used to successfully design
experimental methodologies to determine effective soil hydraulicproperties of a field site regardless of the depth and radius of a drywell.
However, we acknowledge that validation of the model using site
specific infiltration experiment is essential but goes beyond the scope of
this study.

example, a large lateral extension of the wetting front is apparent in thehighly permeable S layer, followed by SL, SCL, SiL, and SC. Note,however, that the lateral extension of the wetting front in a specificlayer is also affected by its location in the profile and neighboring two
soil layers. For example, when the sand (Ks = 4.95 × 10−3 m min−1)
layer is on top of the SC layer (Ks = 2.0 × 10−5 m min−1), water extents laterally in the SC layer despite of its low permeability. The interface between the SC and S layers provides a large infiltration area forthe sandy layer. On the other hand, when the SCL(Ks = 2.18 × 10−4 m min−1) layer is present on top of the SL(Ks = 7.37 × 10−4 m min−1) layer, which has a similar permeability,
the two layers facilitates the vertical movement of water over time
(Fig. 8D and E). In contrast to Fig. 2, these results demonstrate that thepresence of a high permeable layer in a more heterogeneous soil profile

3.2.2. Deterministic heterogeneity – Highly permeable lens
The infiltration capacity of a drywell is not only determined by thepresence of a highly permeable lens but also by its lateral extent(Fig. 2). However, the inverse optimization technique discussed above
only considered layered soil profiles surrounding the drywell, but not

Table 2The summary of various numerical experiments with different positions of soil layers (SL – Sandy loam, SC – Sandy clay, SiL – Silt loam, SCL – Sandy clay loam, S –Sand; subscript B and T for bottom and top, respectively). The saturated isotropic hydraulic conductivity (K s ) from the HYDRUS (2D/3D) Soil Catalog (K s HYDRUS) and
fitted (K s fit), the standard or coefficient of K s (SE Ks ), the shape parameter from the HYDRUS (2D/3D) Soil Catalog ( HYDRUS) and fitted ( fit), the standard error
coefficient of (SE ), the coefficient of determination between the direct simulated and inversely optimized cumulative infiltration volume data (R2), and the massbalance error (M . B. Err ).
Experiment

Soil layers

Soil
[−]

K s HYDRUS[m min−1]

K s fit[m min−1]

SE K s
[−]

[m−1]

[m−1]

SE
[−]

R2
[−]

M. B. Err[%]

I (falling head experiments)

SLB + SCT

SL
SC
SC
SL

7.37 × 10−4
2.0 × 10-5
2.0 × 10-5
7.37 × 10-4

7.37 × 10−4
1.91 × 10-5
2.01 × 10-5
7.10 × 10-4

9.71 × 10-8
2.90 × 10-7
1.63 × 10-8
7.43 × 10-6

7.5
2.7
2.7
7.5

7.49
2.67
2.70
5.33

4.92 × 10-3
3.38 × 10-2
2.70 × 10-3
7.43 × 10-2

1
0.99
0.99
0.99

0.15
0.86
0.0037
0.02

SL
SC
SC
SL

7.37 × 10−4
2.0 × 10−5
2.0 × 10−5
7.37 × 10−4

7.37 × 10−4
1.98 × 10−5
2.00 × 10−5
7.37 × 10−4

4.75 × 10−7
2.39 × 10−7
3.65 × 10−88.90 × 10−7

7.5
2.7
2.7
7.5

7.47
2.71
2.70
7.51

5.28 × 10−2
7.51 × 10−22.40 × 10−2
8.99 × 10−2

1
1
1
1

0.007
0.01
0.00710.0099

SiL
SL
SCL
SC
S

7.50 × 10−5
7.37 × 10−4
2.18 × 10−4
2.0 × 10−5
4.95 × 10−3

7.50 × 10−5
7.37 × 10−4
2.33 × 10−4
1.33 × 10−4
4.43 × 10−3

4.36 × 10−8
6.64 × 10−7
1.12 × 10−6
4.85 × 10−7
6.07 × 10−6

2.0
7.5
5.9
2.7
14.5

2.0
7.46
6.66
1.32
14.43

8.0 × 10−3
0.11
0.22
0.06
0.48

0.99
0.99
0.99
0.99
1

0.11
0.03
0.02
0.02
0.03

SCB + SLTII (constant head experiments)

SLB + SCTSCB + SLTIII (constant head experiments)

SiLB
SiLB + SLTSiLB + SL + SCLTSiLB + SL + SCL + SCT
SiLB + SL + SCL + SC + ST

609

HYDRUS

fit

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

the lateral extension of a lens. Additional, numerical experiments were
therefore conducted to investigate the ability to identify the lateralextension of a sandy lens that embedded in a sandy clay domain. Thesesimulations considered a constant head boundary condition in thedrywell reservoir and different horizontal extensions (0, 2.5, 5, 7.5, and10 m) of a 0.5 m thick sandy lens. Fig. 9A and B show the cumulative
infiltration volumes (I ) and infiltration rates (i ) as a function of time inthese simulations. Note that I increases with the length of the sandylens in Fig. 9A. The value of i also decreases precipitously when the
sandy lens gets fully saturated (Fig. 9B). For example, Fig. 9C (at800 min) shows the distribution after 800 min when the lateral extension of the sandy layer is 10 m. In this case, the sandy lens has become completely saturated (at ∼761 min), while the corresponding
infiltration rate (blue line) in Fig. 9B shows a sudden decrease. However, it should be mentioned that it was much more difficult to discern
the influence of the horizontal length of lens on I and i when propertiesof sandy clay loam were assigned to the lens (rather than sand), when aconstant head boundary was assigned at the top of the sandy lens (rather than assuming a full reservoir), and for greater subsurface heterogeneity in soil hydraulic properties.

played a secondary role in determining the cumulative infiltration
volume and long-term performance of a drywell. A large Ks value indicates the presence of a highly conductive soil that leads to the faster
infiltration. Whereas, a large
value indicates a uniform and sharpwetting front and a small value indicates a nonuniform and dispersed
wetting front (Carsel and Parrish, 1988; Elrick et al., 1989; Rawls et al.,
1982; van Genuchten et al., 1991). Additional numerical experiments
are needed to understand this topic better, which was not further investigated in this study.
Fig. 13 shows a comparison of I simulated for a homogeneous sandy
loam soil profile with Ks and values obtained from the HYDRUS (2D/
= 7.5 m−1)
3D) Soil catalog (Ks = 7.37 × 10−4 m min−1 and(Table 1), inversely optimized average value of Ks = 3.94 × 10−3 m
min−1 and = 6.17 m−1, Upper-Confidence-Interval values of
Ks = 3.68 × 10−3 m min−1 and = 5.40 m−1, and Lower-ConfidenceInterval values of Ks = 4.20 × 10−3 m min−1 and = 6.94 m−1 for theheterogeneous domain with the following stochastic parameters: = 1,
X = 10, and Z = 0.1. Average values of Ks and I increased by 435%
and 396% for the stochastic domain. Note that all these simulations
were conducted for a homogeneous flow domain with a sandy loam soilfor 5760 min (96 h) by keeping a constant head at the surface of thedrywell.
Replacing the heterogeneous domain with an equivalent homogeneous domain as presented in this study can accurately predict longterm cumulative infiltration (96 h), which is relevant for the Low
Impact Development Standard for drywells (LACDPW, 2014). Therefore, a detailed characterization of heterogeneity in soil hydraulicproperties, which is very time consuming and expensive, may not be
required when water quantity is the only criterion for evaluating drywell performance. This research demonstrates that falling or constanthead drywell infiltration experiments can be used in conjunction with
HYDRUS (2D/3D) numerical modeling to inversely estimate effective
hydraulic properties for a heterogeneous drywell site that can be used
to accurately assess its infiltration performance. A constant head experiment can accurately provide the average hydraulic properties of adrywell site. However, a constant head experiment may not account forhydraulic properties (position, lateral extension, and thickness of highpermeable lenses) of sites presented in Fig. 2 that are highly relevant
factors affecting the cumulative infiltration volume in a falling headexperiment. The hydraulic properties determined using a falling headexperiment may thus underestimate or overestimate the cumulative
infiltration volume compared to the hydraulic properties obtained froma constant head experiment. Therefore, this research demonstrates thata constant head drywell infiltration experiment coupled with HYDRUS(2D/3D) inverse optimization of hydraulic parameters of each layer of amulti-layered soil system (Table 2 and Fig. 8) is the best method toaccurately estimate the hydraulic property of a drywell site and to assess the long-term performance of a drywell.

3.2.3. Stochastic heterogeneity
Although the above results demonstrate that the soil hydraulicproperties of individual soil layers can be accurately determined during
constant head infiltration experiments, the noted limitations indicatethat it is unlikely to accurately determine the lateral extension of theselayers/lenses under heterogeneous field conditions. The heterogeneous
stochastic simulations discussed in Figs. 3 and 5 were therefore analyzed in terms of equivalent homogeneous profile values of Ks andthat were inversely fitted to I for each stochastic realization. Averagevalues of Ks and , and their 95% confidence intervals, were determined from 100 realizations using a given set of stochastic parameters. Fig. 10 presents average values of Ks as a function of the stochastic parameters (Fig. 10A with X = 1 and Z = 0.1), X (Fig. 10B
with = 1 and Z = 0.1, and Fig. 10C with = 1 and Z = 2), and Z(Fig. 10D with = 1 and X = 1) and the Ks value of sandy loam fromthe HYDRUS (2D/3D) Soil catalog. The average values of Ks were foundto be highly sensitive to different stochastic parameters and exhibitedvery narrow 95% confidence intervals. In particular, average values of
Ks increase with and X , but decrease with Z. Similar to I in Figs. 3 and5, this can be attributed to increasing amounts of high permeability
and X , andlenses with greater lateral extensions with increasing
isolated high permeability zones with increasing Z. The average Ks
value for stochastic heterogeneous domains was always higher than thehomogeneous domain (Fig. 12A). A highly heterogeneous (σ = 1) domain with many high and low permeable lateral lenses connectedvertically via short highly permeable lenses (i.e., X = 10, Z = 0.1;
X = 10, Z = 2; X = 0.1, Z = 0.1; and X = 1, Z = 0.1) had the highest
average Ks value (Fig. 12A).
Fig. 11 presents average values of as a function of the stochastic
parameters (Fig. 11A with X = 1 and Z = 0.1), X (Fig. 11B with = 1
and Z = 0.1, and Fig. 11C with = 1 and Z = 2), and Z (Fig. 11D with
= 1 and X = 1) and the value of sandy loam from the HYDRUS (2D/
3D) Soil catalog. Average values ofranged between 2.44 and8.46 m−1 and had a relatively large confidence interval (Fig. 12B)
compared to the average values of Ks with a very tight confidence interval (Fig. 12A). Furthermore, a sensitivity analysis ofKs and wereconducted where one of the parameters (i.e. either Ks or ) was kept asconstant in each direct simulation and the cumulative infiltration volume was determined (Data is not shown). Average value of
Ks = 3.94 × 10−3, 8.19 × 10−4
and 2.44 × 10−4 m min−1 and,= 8.46, 7.5, and 2.44 m−1 were used (Fig. 12). The results from thesensitivity analysis demonstrated that the cumulative infiltration volume increases with increase in Ks and values at a given condition.
However, the observations from sensitivity analysis, Figs. 10–12 indicate that the average value of Ks plays a primary role and the value of

4. Conclusions
The numerical modeling results presented in this study demonstrate
the application of HYDRUS (2D/3D) in predicting the cumulative infiltration volume and in determining soil hydraulic properties adjacent
to a drywell. In particular, numerical experiments were conducted tobetter understand the influence of deterministic and stochastic heterogeneity in soil hydraulic properties on drywell infiltration behavior.
Results demonstrated that a highly permeable lens/layer at the bottomof a drywell can infiltrate water at a much faster rate than the same
sized lens at the top of a drywell. Furthermore, the cumulative infiltration volume from a drywell increases with an increasing lateralextension and thickness of a high permeable lens. Similarly, the average
cumulative infiltration volume from 100 stochastic realizations increased with increasing and lateral correlation length but decreased
with an increasing vertical correlation length.
The result from direct and indirect numerical experiments provided
610

Journal of Hydrology 570 (2019) 598–611

S. Sasidharan et al.

insight in developing field scale falling head and constant head experimental methods for drywell to determine the in-situ soil hydraulicproperties of multiple soil layers and high permeable lenses to characterize drywell infiltration. Numerical experiments demonstrated that
constant head experiments provide higher accuracy in inversely optimized hydraulic parameters than falling head experiments. However,
inverse determination of soil hydraulic properties for individual soil
lenses of varying lateral extensions may be impractical for highly heterogeneous systems. The hydraulic properties for stochastic experiments with heterogeneous soil distributions were therefore characterized in terms of an equivalent homogeneous soil profile. The averagevalue of Ks from multiple stochastic realizations with the same parameter sets had the greatest influence and played a primary role on
drywell cumulative infiltration. The value showed a larger confidence
interval and thus played a secondary role in the overall drywell performance. Similar to cumulative infiltration, the average value of Ks
increased with and the lateral correlation length, but decreased withan increasing vertical correlation length.
The equivalent homogeneous hydraulic properties obtained from aheterogeneous domain can be used to estimate the infiltration from adrywell. However, the hydraulic properties obtained from a fallinghead experiment may over or underestimate the cumulative infiltration
volume compared to the hydraulic properties obtained from a constanthead experiment. Therefore, a constant head experiment for multilayered soil system would be the best method to accurately determine
the drywell site hydraulic property and its long-term performance. This
research demonstrates that extensive soil characterization may not benecessary if water quantity is the only concern. However, if the drywell
performance also needs to account for water quality, including the fate
and transport of pathogenic microbial contaminants such as virus and/or bacteria, an accurate estimation of site heterogeneity is expected tobe very important. Ongoing research in our lab is currently investigating this further.

Hodges, M.K.V., Orr, S.M., Potter, K.E., LeMaitre, T., 2012. Construction diagrams, geophysical logs, and lithologic descriptions for boreholes, Idaho National Laboratory,
USGS 103, 105, 108, 131, 135, NRF-15, and NRF-16.
Jain, P., Powell, J., Townsend, T.G., Reinhart, D.R., 2006. Estimating the hydraulicconductivity of landfilled municipal solid waste using the borehole permeameter test.
J. Environ. Eng. 132 (6), 645–652.
Kummu, M., Guillaume, J.H.A., De Moel, H., Eisner, S., Flörke, M., Porkka, M., Siebert, S.,
Veldkamp, T.I.E., Ward, P.J., 2016. The world’s road to water scarcity: shortage andstress in the 20th century and pathways towards sustainability. Sci. Rep. 6, 38495.
https://doi.org/10.1038/srep38495.
LACDPW, 2014. Low Impact Development. County of Los Angeles Department of Public
Works Standards Manual (2014) 1-496.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear parameters. J. Soc. Ind. Appl. Math. 11 (2), 431–441. https://doi.org/10.1137/0111030.
Miller, E., Miller, R., 1956. Physical theory for capillary flow phenomena. J. Appl. Phys.
27 (4), 324–332. https://doi.org/10.1063/1.1722370.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 12 (3), 513–522. https://doi.org/10.1029/
WR012i003p00513.
Philip, J., 1985. Approximate analysis of the borehole permeameter in unsaturated soil.
Water Resour. Res. 21 (7), 1025–1033.
Postel, S.L., 2000. Entering an era of water scarcity: the challenges ahead. Ecol. Appl. 10(4), 941–948. https://doi.org/10.2307/2641009.
Postel, S.L., Daily, G.C., Ehrlich, P.R., 1996. Human appropriation of renewable fresh
water. Science 271 (5250), 785–788. https://doi.org/10.1126/science.271.5250.
785.
Rassam, D., Šimůnek, J., Van Genuchten, M.T., 2003. Modeling variably saturated flow
with HYDRUS-2D. ND Consult.
Rawls, W.J., Brakensiek, D., Saxtonn, K., 1982. Estimation of soil water properties. Trans.
ASAE 25 (5), 1316–1320.
Reynolds, W., Elrick, D., 1985. In situ measurement of field-saturated hydraulic conductivity, sorptivity, and the [alpha]-parameter using the Guelph permeameter. Soil
Sci. 140 (4), 292–302.
Rosegrant, M.W., Cai, X., Cline, S.A., 2002. World Water and food to 2025: dealing withscarcity. Int. Food Policy Res. Inst.
Roth, K., 1995. Steady state flow in an unsaturated, two-dimensional, macroscopicallyhomogeneous miller-similar medium. Water Resour. Res. 31 (9), 2127–2140.
Roth, K., Hammel, K., 1996. Transport of conservative chemical through an unsaturated
two-dimensional Miller-similar medium with steady-state flow. Water Resour. Res.
32 (6), 1653–1663. https://doi.org/10.1029/95WR00946.
Sasidharan, S., Bradford, S.A., Šimůnek, J., DeJong, B., Kraemer, S.R., 2018. Evaluating
drywells for stormwater management and enhanced aquifer recharge. Adv. Water
Resour. 116, 167–177. https://doi.org/10.1016/j.advwatres.2018.04.003.
Šejna, M., Šimůnek, J., van Genuchten, M.T., 2014. The HYDRUS Software Package forSimulating the Two- and Three-Dimensional Movement of Water, Heat, and Multiple
Solutes in Variably-Saturated Porous Media, User Manual, Version 2.04. PC-Progress,
Prague, Czech Republic. 1-307.
Shiklomanov, I.A., 2000. Appraisal and assessment of world water resources. Water Int.
25 (1), 11–32. https://doi.org/10.1080/02508060008686794.
Šimůnek, J., van Genuchten, M. Th., and Šejna, M., 2011. The HYDRUS Software Package
for Simulating Two- and Three-Dimensional Movement of Water, Heat, and Multiple
Solutes in Variably-Saturated Porous Media, Technical Manual, Version 2.0, PC
Progress, Prague, Czech Republic, pp. 258.
Šimůnek, J., Šejna, M., van Genuchten, M.T., 2018. New features of version 3 of the
HYDRUS (2D/3D) computer software package. J. Hydrol. Hydromech. 66 (2),
133–142. https://doi.org/10.1515/johh-2017-0050.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent developments and applications of the HYDRUS computer software packages. Vadose Zone J. 15 (7), 25. https://doi.org/10.2136/vzj2016.04.0033.
Tseng, P.H., Jury, W.A., 1994. Comparison of transfer function and deterministic modeling of area-averaged solute transport in a heterogeneous field. Water Resour. Res.
30 (7), 2051–2063. https://doi.org/10.1029/94WR00752.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898.
van Genuchten, M.T., Leij, F., Yates, S., 1991. The RETC code for quantifying the hydraulic functions of unsaturated soils. IAG-DW12933934, U.S. Salinity Laboratory, U.
S. Department of Agriculture, Agricultural Research Service, Riverside, California.
Vereecken, H., Kasteel, R., Vanderborght, J., Harter, T., 2007. Upscaling hydraulicproperties and soil water flow processes in heterogeneous soils. Vadose Zone J. 6 (1),
1–28. https://doi.org/10.2136/vzj2006.0055.
Watt, E., Marsalek, J., 2013. Critical review of the evolution of the design storm event
concept. Can. J. Civ. Eng. 40 (2), 105–113. https://doi.org/10.1139/cjce-2011-0594.
Weeks, E.P., 1978. Field determination of vertical permeability to air in the unsaturated
zone. Department of the Interior, Geological Survey.
Xiang, J., Chen, L., 1993. Users’ manual for LAYERK.VI-a code to evaluate the tield hydraulic conductivity layered media in unsaturated for zone. Bureau of Economic
Geology, The University of Texas at Austin.
Xiang, J., Scanlon, B., Mullican, W., Chen, L., Goldsmith, R., 1997. A multistep constanthead borehole test to determine field saturated hydraulic conductivity of layeredsoils. Adv. Water Resour. 20 (1), 45–57. https://doi.org/10.1016/S0309-1708(96)00017-6.

Acknowledgments
Funding for this research was provided by the U.S. Environmental
Protection Agency (US EPA) through an interagency agreement with
the United States Department of Agriculture (EPA DW-012-92465401;
ARS 60-2022-7-002). The views expressed in this article are those of theauthors and do not necessarily represent the views or policies of the
U.S. Environmental Protection Agency. Mention of commercial products does not constitute an endorsement.
References
Archer, N., Bonell, M., MacDonald, A., Coles, N., 2014. A constant head well permeameter
formula comparison: its significance in the estimation of field-saturated hydraulicconductivity in heterogeneous shallow soils. Hydrol. Res. 45 (6), 788–805.
Bouwer, H., 2002. Artificial recharge of groundwater: hydrogeology and engineering.
Hydrogeol. J. 10 (1), 121–142. https://doi.org/10.1007/s10040-001-0182-4.
Bouwer, H., Back, J.T., Oliver, J.M., 1999. Predicting infiltration and ground-watermounds for artificial recharge. J. Hydrol. Eng. 4 (4), 350–357. https://doi.org/10.
1061/(ASCE)1084-0699(1999) 4:4(350).
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24 (5), 755–769.
Dillon, P., 2005. Future management of aquifer recharge. Hydrogeol. J. 13 (1), 313–316.
https://doi.org/10.1007/s10040-004-0413-6.
Edwards, E.C., Harter, T., Fogg, G.E., Washburn, B., Hamad, H., 2016. Assessing the effectiveness of drywells as tools for stormwater management and aquifer recharge andtheir groundwater contamination potential. J. Hydrol. 539, 539–553. https://doi.
org/10.1016/j.jhydrol.2016.05.059.
Elrick, D.E., Reynolds, W.D., Tan, K.A., 1989. Hydraulic conductivity measurements inthe unsaturated zone using improved well analyses. Ground Water Monit. R 9 (3),184–193. https://doi.org/10.1111/j.1745-6592.1989.tb01162.x.
Hammel, K., Roth, K., 1998. Approximation of asymptotic dispersivity of conservative
solute in unsaturated heterogeneous media with steady-state flow. Water Resour. Res.
34 (4), 709–715. https://doi.org/10.1029/98WR00004.

611