Water Research 247 (2023) 120781

Contents lists available at ScienceDirect

Water Research
journal homepage: www.elsevier.com/locate/watres

Estimating the impact of vadose zone heterogeneity on agricultural
managed aquifer recharge: A combined experimental and modeling study
Tiantian Zhou a, Elad Levintal b, Giuseppe Brunetti c, Spencer Jordan a, Thomas Harter a,
Isaya Kisekka a, Jiří Šimůnek d, Helen E. Dahlke a, *a

Department of Land, Air and Water Resources, University of California, Davis, CA, 95616, United States of America
Zuckerberg Institute for Water Research, Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev, Sde Boker campus, 84990, Israel
Department of Civil Engineering, University of Calabria, Rende, Italyd
Department of Environmental Sciences, University of California, Riverside, CA, 92521, United Statesbc

A R T I C L E I N F O

A B S T R A C T

Keywords:
HYDRUS-1D
Dual porosity model
Preferential flow
Subsurface runoff
Bromide transport
Travel time

Agricultural managed aquifer recharge (Ag-MAR) is a promising approach to replenish groundwater resources
using flood water and cropland as spreading grounds. However, site selection, particularly the layering of
sediment deposits in the subsurface, can greatly influence Ag-MAR efficacy as it controls water flow and solutetransport in the vadose zone. In this study, we use the HYDRUS-1D software to simulate water flow and solutetransport from the land surface to the groundwater table in three vadose zone profiles (LS, MS, HS) characterized
by differing fractions of sand (44 %, 47 %, and 64 %). For each profile, the single- and dual-porosity models (i.e.,
considering or not nonequilibrium water flow and solute transport) were calibrated using observed surface
ponding, soil water content, and KBr breakthrough data. Water flow and bromide transport in the profile with thelowest sand fraction (LS) were best captured using the model that considered both preferential flow andnonequilibrium bromide transport. Water flow and bromide transport in the profile with the highest sandfraction (HS) was best simulated with the model that considered preferential flow and equilibrium bromidetransport. Uniform water flow and nonequilibrium bromide transport provided the best fit for the third profile(MS). The degree of preferential flow was highest in the profile with the largest sand fraction (HS), which also
showed the largest flow velocities compared to the profiles with lower sand amounts (LS and MS). Preferential
flow did not significantly impact the overall water balance (within 3 %), but caused a significant decrease invadose zone travel times (bromide) by up to 38 %, relative to a single-porosity model fit. Recharge efficiency
varied between 88 % and 90 %, while the average travel times from the soil surface to groundwater varied up to
119 % (from 3.6 to 7.9 days) between the three sites. This study demonstrates that similar recharge efficiency canbe achieved at sites with differing soil texture profiles, but subsurface heterogeneity can substantially affectcontaminant transport processes and their travel times.

1. Introduction(Sprenger et al., 2017; Wendt et al., 2021).
Agricultural managed aquifer recharge (Ag-MAR) is a fairly newpractice that recharges groundwater by transferring excess surface water
onto agricultural lands (Kocis and Dahlke, 2017). This practice can playa dual role in increasing agricultural water supply by replenishinggroundwater and improving groundwater quality by diluting pollutant
concentrations in groundwater (Bali et al., 2023). Ag-MAR is also lesscostly because it utilizes the existing irrigation infrastructure and farm
fields as spreading sites without the additional construction of rechargebasins or wells (Kourakos et al., 2019). Because of these advantages andGroundwater contributes about 25 % of irrigation water use, 50 % ofdomestic water use, and 40 % of industrial water use (Siebert et al.,2010; WWDR, 2022). In semi (arid) regions, groundwater withdrawal
rates may exceed aquifer recharge rates, which causes groundwater
depletion and associated environmental issues (Scanlon et al., 2012;
Siebert et al., 2010)(Marwaha et al., 2021). Increasing groundwater
recharge through managed aquifer management (MAR) is a promising
approach for the sustainable development of groundwater resources* Corresponding author.
E-mail address: hdahlke@ucdavis.edu (H.E. Dahlke).
https://doi.org/10.1016/j.watres.2023.120781
Received 26 April 2023; Received in revised form 9 October 2023; Accepted 22 October 2023
Available online 23 October 2023
0043-1354/© 2023 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

T. Zhou et al.

Water Research 247 (2023) 120781

vast agricultural lands available in many groundwater-dependent re­gions, a boom in Ag-MAR implementation has been witnessed in recent
years, especially in the USA and Europe (Levintal et al., 2023b). How­ever, there is uncertainty about how to identify suitable Ag-MAR loca­
tions that allow recharge of large water volumes but do not create
additional environmental issues (Alam et al., 2021). One of the mostproblematic issues is the leaching of nitrate, salts, pesticides, and mi­
crobes from surface water or farmlands to groundwater (Bachand et al.,2014; Guo et al., 2023; Levintal et al., 2023a; Murphy et al., 2021;
Waterhouse et al., 2020).
The prerequisite for addressing the above issues is an accurate un­derstanding of water flow and solute transport processes in the vadosezone during Ag-MAR and their influencing factors. The influencing
factors can be, in general, divided into three groups: a) hydrometeoro­
logical conditions, such as the intensity, amount, duration, and qualityof precipitation or the applied water and the evaporation rate; b)hydrogeological conditions, including the vadose zone thickness, soil
texture, and structure, soil chemical composition, and c) the type of land
use and its management (Dahlke et al., 2018; Perzan et al., 2023; Siebert
et al., 2010). Field experiments cannot always observe solute transport
processes and groundwater contamination risks at the spatial-temporalscales at which they occur or might introduce the contaminant (Sasid­haran et al., 2021; Wang et al., 2022, 2020). Numerical modeling canprovide hydrologic information at a high spatial-temporal resolution,
but the model calibration is often highly dependent on experimental
data and is subject to the equifinality problem (Zhou et al., 2023). In
contrast, field experiments combined with numerical models, such as
HYDRUS, are more popular tools to investigate the impacts of differenthydrometeorological and hydrogeological factors on groundwater
recharge quantity and quality and associated crop root zone status. For
example, Bali et al. (2023) used HYDRUS-2D to simulate the waterbalance of an alfalfa field in the San Joaquin Valley comparing different
irrigation treatments in the summer (full and deficit irrigation) andwinter flooding for groundwater recharge. Their results show that the
fully irrigated alfalfa stand that was flooded in winter for groundwater
recharge showed recharge efficiencies of 85 %, 89 % and 84 % in 2020,2021 and 2022, while those values were only 78 %, 79 % and 76 % for
the deficit irrigation treatment. Similarly, Ganot and Dahlke (2021) used
HYDRUS-1D to validate an analytical solution-based root zone residence
time model to develop guidance on the safe flooding duration for
Ag-MAR for different soil textures (Ganot and Dahlke, 2021). Bothstudies estimated water balance components for different treat­ments/scenarios using a daily time step. However, to understand theeffects of Ag-MAR on the fate and transport of solutes such as nitrate orsalts, high-resolution (e.g., 10 min) data and numerical modeling areneeded to fully quantify the complete recharge cycle of water andassociated solutes from the soil surface through the vadose zone to thegroundwater table.
In addition, preferential flow (a.k.a. nonequilibrium or nonuniform
water flow in which the moisture front can propagate quickly to sig­
nificant depths while bypassing a large part of the matrix pore-space)(Šimůnek et al., 2003) can lead to the earlier arrival of pathogens(Bradford et al., 2017), bacteria such as Escherichia coli (Arnaud et al.,2015), and other contaminants at the groundwater table. These adverse
effects may pose a higher risk under Ag-MAR because of intensive ap­
plications of fertilizers and pesticides on agricultural lands. Preferential
flow is often poorly characterized in field studies (Nimmo et al., 2021)
and rarely considered in current numerical modeling analyses of
Ag-MAR (Levintal et al., 2023b). Quantifying the degree of preferential
flow and how it affects recharge and contaminants’ transit times is oneof the top priorities for selecting appropriate Ag-MAR sites.
In this study, we simulate water flow and solute transport in thevadose zone of three Ag-MAR plots in California at a temporal resolutionof 10 min, using the HYDRUS-1D software, considering both single- anddual-porosity models. All three plots are located within the same 15 haagricultural field, thus allowing the assessment of how within-fieldsubsurface heterogeneity affects Ag-MAR vadose zone flow and trans­
port processes. The specific objectives are to (1) identify and quantify towhat degree preferential flow controls water flow and solute transport in
the vadose zone under continuous ponded conditions, (2) quantify theflow of the recharge water, (3) estimate transit times of recharge waterand contaminants from the land surface to the groundwater table using
multiple indicators (bromide tracer, water table dynamics, soil andgroundwater salt dynamics, and soil aeration dynamics), and (4) inte­
grate hydrological and biogeochemical field and simulation data toassess Ag-MAR suitability.
2. Materials and methods
2.1. Study area
The Ag-MAR experiments were conducted in the source areas ofthree groundwater wells (LS, MS, HS) at a fallowed almond orchard westof Modesto, California (37◦ 37″N, 121◦ 05′W, Fig. 1). The area has aMediterranean climate with mean annual precipitation of 330 mm, andpotential evaporation of about 1356 mm (https://cimis.water.ca.gov/).
Sediments in the vadose zone belong to the distal portions of the Tuo­lumne River alluvial fan, with channel sands as the base and heteroge­
neous floodplain sediments at the top. Groundwater flows from
northeast to southwest (Fig. 1a) (Gurevich et al., 2021). The day before
the flooding for Ag-MAR was started, groundwater table depths in thethree wells were 6.34, 6.41, and 6.06 m, respectively. Therefore, we
consider 7 m deep soil profiles in this study. While dominated by silty
clay and silty clay loam, the profiles show progressively increasing
fractions of sand from LS (44 %), MS (47 %) to HS (64 %) (Fig. 2). The
detailed soil particle distributions are shown in Table S1.
2.2. Field monitoring and data collection
Recharge plots shown in Fig. 1a were bermed with native soil (0.5 m
height, 1 m width). The field monitoring lasted from April 20 to June 14,2022. The plots were flooded with surface water for 28 days, from May 3(13:30) to May 31, 2022. Flooding was manually operated from May
3–5, 2022 to determine an irrigation schedule that would maintain aconstant head (ponded water) on the field. During this period, the irri­
gation rate was considered continuous and constant. From 12:00 on May
6 to 8:00 on May 10, water was applied at all three sites for 2-hour in­tervals: 12:00–14:00, 18:00–20:00, 00:00–02:00, and 6:00–8:00. From
12:00 on May 10 to 7:40 on May 31, water applications were reduced to
100-minute intervals: 12:00–13:40, 18:00–19:40, 00:00–01:40, and
6:00–7:40. The flooding stopped at 9:00 on May 31, 2022. Totals of
about 5750 m3, 6281 m3 and 5969 m3 of water were applied within aflooding area of 694 m2, 735 m2, and 690 m2 (i.e., 829.7, 855.8, and866.3 cm in water depth) at LS, MS, and HS, respectively. Irrigation rates
were calculated by dividing the total volume of applied water by theflooded surface area and the irrigation duration at corresponding stages(Table 1). As can be seen, between 5/6/2022 12:00 and 5/10/2022
11:43, water was applied at a higher rate. As a result, plots experienced aleak and some of the applied water spilled outside the berms overnight.
Because the amount of water that spilled is unknown, a scaling factor
was introduced to reduce the applied water amount on those days suchthat the simulated and observed ponding matched. As a result, scaled
water application amounts reduced from 0.0850 to 0.0572 cm/min at
LS, from 0.0820 to 0.0608 cm/min at MS, and from 0.0858 to 0.0613
cm/min at HS.
The bromide application area was about 7 m from each site’s cor­
responding monitoring well (Fig. 1b). Br- was applied up-gradient toensure it would move to the well area once it reached the groundwater.
At each site, 4.833 kg of KBr (i.e., 3.245 kg of Br− ) was dissolved inapproximately 100–110 gallons (380 L) of water and then applied overan area of 36.8 m2 (i.e., 9.2 m* 4 m) for 20 min (starting at 14:00, 12:40,
and 10:52 on May 2, 2022 in recharge plots LS, MS, and HS, respectively,2

T. Zhou et al.

Water Research 247 (2023) 120781

Fig. 1. Location of the study site in California, USA. Locations of the three recharge plots upgradient of the three groundwater wells LS, MS and HS (a), of
instrumentation and the potassium bromide (KBr) application site in each recharge plot (b), of sensors or suction cups at each profile (c). The blue arrow in (a) shows
the regional groundwater flow direction.

i.e., one day before flooding) using a custom made drip irrigation system
and a water tank. Br− was thus applied at a concentration of 8164.18
ppm at an irrigation rate of 0.0540 cm/min.
The meteorological data, including daily precipitation (P) and po­
tential evaporation (Ep) (Fig. S1) were obtained from the CIMIS website(https://cimis.water.ca.gov/) from station 71 – Modesto, located 7 km
west of the study site. The total potential evaporation during themonitoring period was 35.1 cm (Fig. S1).
The ponding level at the soil surface and the groundwater table depth
were measured using pressure transducers (CS451, Campbell Scientific,
UT, USA). Each recharge plot was instrumented with sensors at depths of20, 60, 100, 250, and 500 cm at profile LS, 20, 60, 100, 300, and 450 cm
at MS, and 20, 60, 100, 275, and 430 cm at HS measuring soil watercontent, electrical conductivity (EC), soil temperature (TEROS12, Meter,
WA, USA), O2 (KE-25, Figaro, Japan), and oxidation–reduction potential(ORP) using constructed platinum electrodes. Hydrological (ponding
level, soil moisture, soil temperature, groundwater table depth) andbiogeochemical (O2, ORP, EC) measurements were logged (using
CR1000 and CR800, Campbell Scientific, UT, USA) at a 10-min timeinterval. Breakthrough curve data, monitoring the transport of the KBr
tracer, were collected at minimum 4 h using suction cups (LT-DBL, Irr­
ometer, CA, USA) at the same depths as sensors.
The sensors or suction cups were installed along the long edge of therecharge plot, perpendicular to the groundwater flow direction (Fig. 1b).
Each sensor or suction cup was about 0.5 m apart from the other to avoid
interference. The distance between sensors and corresponding suction

cups was about 1.12 m (Fig. 1c). The temporal distributions of relevant
variables are shown in Figs. S2-S8.
2.3. HYDRUS-1D model setup
2.3.1. Initial and boundary conditions
Water flow and KBr transport in the unsaturated zone were simulatedusing the HYDRUS-1D software (Šimůnek et al., 2016). The
700-cm-thick soil profile was divided into five modeling layers,
including 0–33 cm, 34–66 cm, 67–200 cm, 201–400 cm, and 401–700
cm, in which the third through the fifth modeling layers represent
multiple physical soil layers (Fig. 2). The number of modeling layers was
smaller than the number of physical layers because not all physicallayers had a sensor collecting data. The simulation period was 55 days
long, from April 20 to June 14, 2022, which included pre-flooding,flooding, and post-flooding periods. Since field water contents obser­
vations were limited to five depths, the simulation period was extended
by considering a "spin-up" period (three months) before the flooding
experiment. The spatial discretization resolution was 1 cm throughout
the soil profile, while the temporal discretization resolution was vari­able, with a minimum time step of 0.01 min.
The initial soil pressure head profile was first converted frommeasured water contents at five soil depths by using soil water retention
curves of typical soil textures (Radcliffe and Šimůnek, 2018), and thenlinearly interpolated between any two measurement depths. The initial
Br− solute concentration was zero throughout the soil profile.
3

T. Zhou et al.

Water Research 247 (2023) 120781

Fig. 2. Soil texture in the study profiles.

For water flow, the upper boundary condition (BC) was set to an
atmospheric BC (with a maximum surface water layer of 50 cm, i.e., theheight of the berms). In this BC, when the soil surface is not flooded, thepotential water flux across the soil surface equals the difference between
daily values of potential evaporation, E0, and precipitation (P) or irri­
gation (I). When the soil surface is flooded, the boundary pressure headis equal to the water level at the surface, and the model calculates theinfiltration flux. Depending on the soil moisture status, this atmosphericBC may appear as a Neumann BC (when the surface pressure head iswithin a critical range) or a Dirichlet BC (when the surface pressure headexceeds these critical values). The lower BC was set to a variable

pressure head BC (i.e., Dirichlet BC), defined by the measured position of
the groundwater table.
For solute transport, the upper BC was prescribed as a solute flux BC(i.e., a Cauchy BC), with bromide concentrations and irrigation fluxes
during the bromide application as inputs. The model then automaticallyadjusts surface bromide concentrations depending on the thickness ofthe surface water level and evaporation/precipitation/irrigation fluxes
and associated bromide concentrations. The lower BC was prescribed asa zero concentration gradient (i.e., a Neumann BC when only aconvective solute flux occurs).

4

T. Zhou et al.

Water Research 247 (2023) 120781

Table 1
Irrigation schedule and rates.
LS

Flowmeter
reading (m3)

Floodingrate (cm/
min)

MS

Flowmeter
reading (m3)

Flooding rate
(cm/min)HS

Flowmeter reading(m3 )Flooding rate
(cm/min)5/3/2022
13:30
5/3/2022
17:06
5/6/2022
12:00
5/10/2022
11:43
5/26/2022
12:45
5/31/2022
7:40

00.021500.0231

0.02200.021599.0

0.0231108.0

0.0220708.5

0.0572805.2

0.0608683.1

0.06131841.20.0675

1962.90.07031819.20.07164839.30.06565270.80.06874981.10.07165750.106281.30

5/3/2022
13:30
5/3/2022
16:58
5/6/2022
12:00
5/10/2022
11:52
5/26/2022
12:04
5/31/2022
7:40

0

117.8

5/3/2022
13:30
5/3/2022
17:01
5/6/2022
12:00
5/10/2022
11:39
5/26/2022
12:36
5/31/2022
7:40

5969.40

2.3.2. Single-porosity model (SPM)
Vapor flow can be neglected for conditions considered in thisexample. The one-dimensional uniform soil water movement in
HYDRUS-1D can then be described using the Richards equation:[()]∂θ(h) ∂
∂h
=
K(h)
+1
(1)
∂z
∂t
∂z

2.3.3. Dual-porosity model (DPM)
The dual-porosity model divides the soil pore space into mobile andimmobile regions. Water flow occurs only in the mobile region,
described by the Richards equation, while water can also be stored but
does not flowin the immobile region. The governing equations for waterflow in the dual-porosity model are (Šimůnek et al., 2003):[()]∂θmo (h) ∂
∂h
(8)
=
K(h)
+ 1 − Γw∂z∂t∂z

where θ is the volumetric water content [L3L− 3], t is time [T], h is thewater pressure head [L], z is the spatial coordinate [L] (positive up­ward), and K is the hydraulic conductivity [LT− 1]. The soil water
retention and hydraulic conductivity functions are described using thevan Genuchten-Mualem (VGM) equations (Mualem, 1976; van Gen­uchten, 1980):⎧
θs − θr⎪⎨
θr +m h < 0[1 + |αh|n ]θ(h) =
(2)
⎪
⎩
h≥ 0
θs
[
(
)m ]2
1
k(h) = ks Sel 1 − 1 − Sem(4)m = 1 − 1/n (n > 1)(5)where ωw is the first-order rate coefficient for water transfer between thetwo regions [T − 1], and Se,mo and Se,im are effective saturations in the tworegions [-], respectively. Compared with the single-porosity model, thedual-porosity model additionally considers three parameters, including
the residual (θim,r ) and saturated (θim,s ) water contents in the immobile
region, and ωw . The other parameters are the same as in the singleporosity model and described by the VGM equations, except that now
they are referred to as θmo,r , θmo,s , α, n, and Ks .
The dual-porosity model assumes that solute transport is limited tothe mobile region, but there is a solute transfer between mobile andimmobile regions. The governing equations for solute transport in the
dual-porosity system are:[]∂(θmo Cmo ) ∂∂Cmo∂(qmo Cmo )(11)=
θmo Dmo
−
− Γs
∂z
∂t
∂z
∂z

where θr and θs are the residual and saturated water contents [L3L− 3],
respectively; Ks is the saturated hydraulic conductivity [LT− 1]; Se is theeffective saturation [-]; l is the pore connectivity parameter (about 0.5);
n is an empirical parameter related to the pore-size distribution [-], andα is an empirical parameter related to the inverse of the air-entry suction[L − 1].
The governing equation for solute transport is the advectiondispersion equation:
(
)
∂θC ∂
∂C
∂(qC)
=
θD
−
(6)
∂z
∂t
∂z
∂z
where C is the solute concentrations of soil water (ppm), q is the water
flux [LT− 1], and D is the effective dispersion coefficient of solute in soil
water [L2T− 1] given by:
D = λv +

D0 τ
θ

(9)

where θmo and θim are water contents in the mobile and immobile regions[L3L− 3], respectively. and Γw is the water transfer rate between the tworegions [T − 1], which can be described as (Gerke and van Genuchten,
1993; Šimůnek et al., 2003):()(10)Γw = ωw Se,mo − Se,im(3)

θ − θrθs − θ r

Se =∂θim (h)= Γw∂t∂(θim Cim )
= Γs
∂t

(12)

Γs = ωs (Cmo − Cim ) + Γw c∗(13)

where Dmo , Cmo , and qmo are the dispersion coefficient [L], solute con­centration (ppm), and water flux [LT− 1] in the mobile region [L2T− 1],
respectively, Cim is the solute concentration in the immobile region(ppm), Γs is the solute mass transfer between mobile and immobile re­gions [ML− 3T− 1], and ωs is the solute mass transfer coefficient between
the two regions [T − 1]. c∗ is the solute concentration that depends on the
direction of mass transfer and equals Cmo for Γw >0 and Cim for Γw <0. In
this study, we consider three cases of the dual-porosity model (Fig. 3b, c,
d).

(7)

where λ is the soil dispersivity [L], v is the pore-water velocity [LT− 1], D0
is the molecular diffusion coefficient [L2T− 1], which is about 1.584 cm2/
d for Br− (Isch et al., 2019; Köhne et al., 2004), and τ is the tortuosity
factor [-].

5

T. Zhou et al.

Water Research 247 (2023) 120781

Fig. 3. Conceptual schematics of different model setups (adapted from (Šimůnek and van Genuchten, 2008)).

2.4. Parameter estimation and model performance evaluation

Jacobian approximation of the Hessian matrices around the optima(Šimůnek and Hopmans, 2002). This enabled us to detect and discussparameters’ interaction and test the appropriateness of applying thesingle or dual-porosity models.

The Levenberg-Marquardt algorithm in HYDRUS-1D was used to
optimize soil hydraulic and solute transport parameters. This algorithm
aims to minimize the sum of squared weighted deviations (SSQtotal)
between observed and simulated surface ponding levels (SSQsp), soil
water contents (SSQwc), and bromide concentrations in soil water(SSQBr) (SSQtotal= SSQwc + SSQsp + SSQBr). The square of the correla­tion coefficient (R2), normalized-root-mean-square error (NRMSE), andKling-Gupta efficiency (KGE) were calculated to evaluate the model
performance. While the NRMSE index represents the average deviationof the residuals, R2 measures the linear relationship between simulated
and measured values, and KGE is a comprehensive indicator combining
correlation and bias. The lower the SSQ and NRMSE, and the higher theR2 and KGE, the better the fit between the simulated and observed
values.
In this study, the dispersivity λ was not optimized since preliminarymodel runs indicated that the model performance was not sensitive tothis parameter and only limited solute concentration data were avail­able. The dispersivity was instead assumed to be equal to 70 cm, i.e., 1/
10th of the total travel distance, representing a one-dimensional effec­
tive macrodispersivity (Gelhar et al., 1992).Soil hydraulic and othersolute transport parameters were optimized using measured surface
ponding levels, soil water contents, and bromide concentrations. During
optimizations, parameters were adjusted layer by layer from top tobottom.
In the single-porosity model, the residual water contents were not
optimized. Instead, the default values of corresponding soil textures
were adopted first and then manually adjusted for the model to fit themeasured data better. Therefore, four soil hydraulic parameters (θs , α, n,
Ks ) had to be optimized for each layer. The initial estimates of saturated
water contents were manually set based on the steady water contents
and common values for similar soil textures in Tables S2 (https://structx.com), while initial values of other parameters were obtained from the
Rosetta module in HYDRUS-1D, based on measured average soil particle
distribution data from this orchard (Table S1). The initial parameters forthe soil layer with multiple soil textures were prescribed as those fromthe dominant soil texture.
To reduce the number of optimized parameters in the dual-porosity
model, θm,r was set to zero, as done in many similar studies (Haws
et al., 2005; Imig et al., 2023; Šimůnek et al., 2001). Therefore, eight soil
hydraulic and solute transport parameters (θmo,s , α, n, Ks , θim,r , ωw , ωs )were optimized for each layer of profile LS (40 parameters in total), twosoil hydraulic and solute transport parameters (θim , ωs ) were optimized
for each layer of profile MS (10 parameters in total), and seven soil
hydraulic parameters (θmo,s , α, n, Ks , θim,r , ωiw,r , ωw ) were optimized foreach layer of profile HS (35 parameters in total). The previously opti­
mized parameters of the single-porosity model were used as the initialvalues of parameters in the mobile zone of the dual-porosity model. For
the immobile zone, the initial values θim,r were set the same as those in
the single porosity model, while initial values of θim,r were prescribed as
0.1. The initial values of ωw , ωs were obtained from the literature (Imiget al., 2023; Isch et al., 2019; Köhne et al., 2004).
The correlation matrices for each scenario were calculated using the

3. Results
3.1. Single-Porosity model: parameters and performance
The optimized parameters for the single-porosity model are shown in
Table 2. The soil retention and hydraulic conductivity curves for opti­
mized parameters are shown in Fig. S9. Overall, the optimized param­
eters were within typical values from literatures findings (Text S1). The
exceptions are that the saturated hydraulic conductivities (0.027–0.034
cm/min) for the silty clay layers (0~66 cm) were far higher than the
typical values (0.0062 cm/min in Table S4) according to (Clapp andHornberger, 1978; Li et al., 1976). The saturated water contents for thesilty clay loam and silt loam layers (67~200 cm) were 0.21–0.29
cm3/cm3, below typical values of 0.29–0.52 cm3/cm3 for these texturalclasses, according to Table S2.
The model performance is shown in Table 4 and Figs. 4-6. Overall,
the simulated values and trends of surface ponding levels, soil watercontents, and bromide concentrations matched the observations wellwhen SPM was used. The observed surface ponding levels quickly
increased to their maximum (about 13 cm, 17 cm, and 13 cm for profiles
LS, MS, and HS, respectively) because of intense and continuous irri­
gation at the beginning of the experiment. After that, the ponding level
decreased and remained relatively stable as irrigation became inter­mittent (about every 4 h) and potential evaporation increased (Fig. S1).
The water contents at all depths exhibited increasing trends during
the flooding period and decreased during the post-flooding period.
Simulated wetting fronts arrived later than those observed as depth
increased, especially at LS. The bromide concentrations at all depths firstincreased and then decreased with time. Similarly to wetting fronts,
simulated early breakthroughs at Profiles LS and HS arrived later thanthose observed as depth increased.
Table 2
Optimized parameters of the single-porosity model.
Site

Depth
(cm)

θr (cm3/
cm3)

θs (cm3/
cm3)

α

n (-)

Ks (cm/
min)

LS

0–33
34–66
67–200
201–400
401–700
0–33
34–66
67–200
201–400
401–700
0–33
34–66
67–200
201–400
401–700

0.105
0.115
0.125
0.045
0.151
0.105
0.080
0.095
0.085
0.091
0.050
0.105
0.080
0.025
0.040

0.374
0.301
0.291
0.230
0.528
0.374
0.284
0.210
0.329
0.428
0.305
0.284
0.241
0.210
0.210

0.016
0.012
0.007
0.013
0.061
0.009
0.010
0.006
0.003
0.088
0.011
0.009
0.006
0.013
0.061

1.747
1.305
1.215
3.800
2.000
1.831
1.468
1.598
1.084
3.235
1.232
1.276
1.300
2.943
3.099

0.031
0.029
0.015
9.000
0.200
0.028
0.027
0.015
0.417
0.030
0.034
0.017
0.017
6.000
2.000

MS

HS

6

(cm− 1)

T. Zhou et al.

Water Research 247 (2023) 120781

Fig. 4. Observed and simulated (using the single-porosity model [SPM] and dual-porosity model [DPM]) surface ponding water levels for the three soil profiles (LS,MS, HS). The blue shaded area indicates the flooding period. .

At LS, simulated breakthrough curves (BTCs) at the bottom (250 and
500 cm) showed a slower response to flooding than those observed. In
addition, the observed BTCs had much more significant tailings thanthose simulated. Observed BTCs also displayed secondary peaks at 100
and 250 cm. At HS, simulated BTCs always occurred slightly later thanthose observed, except at depths of 10, 20, and 275 cm. At MS, thearrival of the simulated BTCs matched well with those observed. How­ever, the model could not capture strong tailings at 300 and 450 cm
depths.

caused by diffusive mass transfer between mobile and immobile regions
and the fast bromide front arrival by nonuniform (preferential) waterflow (Isch et al., 2019; Köhne et al., 2004), and therefore, both ωw and ωswere optimized (Fig. 3d).
The optimized parameters for the three dual-porosity models areshown in Table 3. The model performance is shown in Table 4 andFigs. 4-6. DPM provided a slightly better fit to the observed data than
SPM. Notably, the use of the dual-porosity model (considering
nonequilibrium water flow and solute transport) resulted in a slightdecrease in NRMSEs and similar correlation coefficients (R2) for simu­
lated and observed BTCs at profiles LS and MS (Table 4). Although using
the dual-porosity model increased R2, it also increased NRMSEs for BTCs
at Profile HS. The cumulative water and bromide transfer from themobile to immobile zone (CumΓw and CumΓw , a higher value means
higher degree of nonequilibrium flow and solute transport) are shown in
Fig. S16. Therefore, Profile HS showed a higher propensity for prefer­ential flow than LS, while MS showed the least indication of preferentialflow.
The correlation matrices for parameters of different models areshown in Fig. S10-S15. When considering only the strong correlations(R>0.6), only α3 -α4 , Ks1 -Ks3 were negatively correlated when using SPM
at LS, while many more parameters were positively and negatively
correlated when using DPM. This suggests that DPM improved model
performance due to over-parameterization. However, both SPM andDPM structures could not capture the observed bromide BTCs well
(Table 4). At MS, both SPM and DPM showed a positive correlation
between Ks3 -α1 and negative correlations between Ks1 -α1 , Ks1 -α2 ,Ks5 -α2 , Ks3 -Ks1 . Both SPM and DPM performed well and similarly. This
indicates that SPM and DPM structures were equivalent and sufficient indescribing the observations, and it was unnecessary to apply DPM at MS.
At HS, for SPM, n1 -α3 , Ks3 -α3 , Ks3 -α4 were positively correlated, while
α2 -θs,4 , α4 -α2 , Ks1 -α3 , Ks3 -Ks2 , Ks1 -n1 , Ks2 -α4 were negatively correlated.
For DPM, only four parameter pairs were correlated, including Ks1 -α1(positive), and n1 -θs,1 , Ks3 -Ks2 , ωw3 -θim,r5 (negative), and model

3.2. Dual-Porosity model: parameters and performance
As discussed in Section 3.1, at profile MS, arrival times of thesimulated wetting and bromide fronts at depths of 300 and 450 cm
matched very well those observed, but simulated bromide leaching was
faster than observed when SPM was used (i.e., strong tailing of observed
BTCs). Therefore, the special case of the dual-porosity model, i.e., themobile-immobile model (MIM), was used. This model assumes thatwater flow is uniform (i.e., no preferential flow), the immobile watercontent is constant, water mass transfer between the two regions is zero,and solute is transported between the mobile and immobile regions bydiffusive exchange (Isch et al., 2019; Köhne et al., 2004). In other words,
the water transfer coefficient, ωw , was equal to zero, while ωs was
optimized (Fig. 3b).
At profile HS, no significant tailing was observed in the bromide
BTCs, and arrival times of bromide fronts simulated using SPM wereretarded compared to those observed. Thus, it was assumed that theearlier arrival of observed bromide fronts compared to those simulated
was caused only by nonuniform (preferential) water flow (Haws et al.,2005). In this case, ωs was set to zero while ωw was optimized (Fig. 3c).
At profile LS, significant tailing was observed in the bromide BTCs,
and the arrival of bromide fronts simulated using SPM was retarded
compared to those observed. The observed BTCs also displayed sec­
ondary peaks at 100 and 250 cm. In this case, bromide tailing was likely
Table 3
Optimized parameters of the dual-porosity model.
Site

LS

MS

HS

Depth (cm)θmo,r (cm3/cm3)

θmo,s (cm3/cm3)α(cm− 1)n (-)

Ks (cm/min)θim,r (cm3/
cm3)

θim,s (cm3/
cm3)

ωw

0~33
34~66
67~200201~400
401~700
0~33
34~66
67~200201~400
401~700
0~33
34~66
67~200201~400
401~700

0
0
0
0
0
0.105
0.080
0.095
0.085
0.091
0
0
0
0
0

0.254
0.181
0.176
0.103
0.228
0.374
0.284
0.210
0.329
0.428
0.205
0.134
0.134
0.160
0.110

0.012
0.012
0.007
0.014
0.061
0.009
0.010
0.006
0.003
0.088
0.011
0.008
0.006
0.014
0.067

1.487
1.303
1.231
3.096
2.100
1.831
1.468
1.598
1.084
3.235
1.207
1.324
1.473
2.956
2.849

0.033
0.029
0.015
16.033
0.050
0.028
0.027
0.015
0.417
0.030
0.035
0.017
0.016
6.515
2.077

0.077
0.115
0.112
0.048
0.150

0.100
0.130
0.130
0.110
0.450

4.100E-08
3.270E-07
5.870E-07
5.870E-07
2.670E-08

0.031
0.100
0.091
0.030
0.036

0.100
0.150
0.100
0.050
0.100

6.082E-04
2.978E-04
5.291E-03
3.633E-03
6.877E-03

7(min− 1)θim (cm3/
cm3)

6.116E-02
5.359E-02
2.710E-03
5.271E-02
8.007E-02

ωs (min− 1)
1.052E-03
6.458E-03
1.276E-05
2.890E-04
4.614E-06
7.000E-03
7.000E-03
4.616E-03
1.318E-03
3.528E-03

T. Zhou et al.

Water Research 247 (2023) 120781

Table 4The performance of the single-porosity model [SPM] and dual-porosity model [DPM] to simulate surface ponding levels, soil water contents, and bromide concen­
trations for the three soil profiles (LS, MS, HS).
Simulated variable

Indicator

LS
SPM

DPM

MS
SPM

DPM

HS
SPM

DPM

Surface ponding level

R2
NRMSE
KGE
R2
NRMSE
KGE
R2
NRMSE
KGE

0.435
0.276
0.545
0.852
0.201
0.913
0.238
1.150
− 0.046

0.460
0.263
0.544
0.873
0.190
0.903
0.237
1.139
0.047

0.624
0.287
0.656
0.764
0.220
0.849
0.569
1.406
0.446

0.624
0.287
0.656
0.764
0.220
0.849
0.569
1.393
0.458

0.632
0.279
0.334
0.908
0.181
0.939
0.646
1.856
0.412

0.634
0.248
0.337
0.907
0.185
0.935
0.791
1.933
0.592

Water content
Bromide concentration

performance improved a little compared to SPM. This emphasizes thenecessity of employing DPM at HS.

where Lj and Ks,j are the thickness and the saturated hydraulic con­ductivity of each layer (Tables 2~3), respectively.
At LS, MS, and HS, the KsTop,eff values were calculated as 0.018,
0.018, and 0.019 cm/min for SPM, and 0.018 cm/min for DPM. The
corresponding mean ponding depths measured during Ag-MAR were
6.29, 7.72, and 7.14 cm. The groundwater recharge rates were thus
estimated to be 0.019, 0.019, and 0.020 cm/min for SPM and 0.019 cm/
min for DPM. In other words, groundwater recharge was determined by(and close to) KsTop,eff , because gravity gradients prevailed over the
ponding gradients. These values were also consistent with water fluxes
at 200 cm (0.018, 0.019, and 0.020 cm/min at LS, MS, and HS,
respectively for both SPM and DPM) simulated by HYDRUS-1D(Fig. S17), validating their accuracy. Therefore, the soil water balance
was similar between the three soil profiles.

3.3. Water mass balance
The water balance calculation results obtained using the single- anddual-porosity models were very similar (differences were within 2 %), as
shown in Table 5 and Fig. 7b. Regarding the amount of groundwater
recharge as a fraction of the total surface water applied (a.k.a. recharge
efficiency, calculated as (D+ΔSDVZ )/(P + I)), profile HS yielded thelargest recharge, LS the smallest, while MS ranked between the two. All
three profiles showed a similar groundwater recharge efficiency (88 %−90 %). Overall differences in mass balance components between thethree profiles were not very large (within 3 %).
This is because the top three layers (0~200 cm) of all three soil
profiles quickly reached saturation during intensive Ag-MAR (Fig. 5). In
such a case, evaporation from the soil surface was close to potential
evaporation and thus was similar between different sites. Water move­
ment from the soil surface to the bottom of the third layer was only
driven by the ponding and gravity gradients. According to Darcy’s law,soil drainage rates below 200 cm (equivalent to groundwater recharge
since we included soil water storage in the deep vadose zone into
groundwater recharge as well) can be calculated as:()
LTop + SP
GR = − KsTop,eff ∗
(14)
LTop

3.4. Bromide travel time
The peak displacement method estimates travel times from the timelag between peaks in the measured input (irrigation water) and output(soil water at different depths) bromide BTCs (Zhou et al., 2021). The
travel times and average velocities of bromide front from the soil surface
to different soil depths calculated with the peak displacement method
are shown in Table 6. In general, the mean velocities of bromide front
increased as depth increased for both the single- and dual-porosity
models. Due to preferential flow, the travel times of bromide front
from the soil surface to different depths of the soil profiles decreased by
up to 38 %, while the transport velocities increased by up to 61 %,compared to the single-porosity model. Overall, travel times (flow ve­
locities) were longest (slowest) at MS, followed by LS, and shortest(fastest) at HS. The travel times from land surface to groundwater tablevaried from 3.6 to 7.9 days, resulting in an overall average transport
velocity difference between the three sites of up to 119 %. Travel times
and transport velocities inferred from water table dynamics, soil

where LTop is the total thickness of the top three soil layers (200 cm), SP
is the surface ponding depth, KsTop,eff is the effective saturated hydraulicconductivity of the top three layers as follows:∑3
j=1 LjKsTop,eff = ∑3 Lj(15)
j=1 Ks,j

Table 5
Water balance components for different soil profiles.
Term

P+F
L
R
E
D
ΔSRZ
ΔSDVZ
GR

LS
SPM

DPM

cm

%

cm

%

Relativedifference
%

829.7
63.9
0
23.9
687.1
10.7
43.8
730.9

7.7
0
2.9
82.8
1.3
5.3
88.1

829.7
63.9
0
24.9
693.2
11.8
35.7
728.9

7.7
0
3.0
83.5
1.4
4.3
87.9

0
0
0.1
0.7
0.0
− 1.0
− 0.2

MS
SPM

DPM

cm

%

cm

%

Relativedifference
%

855.8
51.5
0
24.7
729.5
20.3
28.7
758.2

6.0
0
2.9
85.2
2.4
3.4
88.6

855.8
51.5
0
24.7
729.5
9.1
39.9
769.4

6.0
0
2.9
85.2
1.1
4.7
89.9

0
0
0
0
− 1.3
1.3
1.3

HS
SPM

DPM

cm

%

cm

%

Relativedifference
%

866.3
53.0
0.0
24.4
738.2
10.0
42.6
780.8

6.1
0
2.8
85.2
1.2
4.9
90.1

866.3
53.0
0
24.3
739.3
10.1
42.2
781.5

6.1
0
2.8
85.3
1.2
4.9
90.2

0
0
0
0.1
0
0
0.1

P: precipitation, F: flood irrigation, L: water loss outside the berms, R: runoff, E: evaporation, D: drainage, ΔS: storage change in the root zone 0~150 cm (ΔSRZ ) anddeep vadose zone (ΔSDVZ ), GR: groundwater recharge including D and ΔSDVZ because water flow is considered to be one-dimensional and thus deep drainage below theroot zone will eventually recharge groundwater with a delay (de Vries and Simmers, 2002).

8

T. Zhou et al.

Water Research 247 (2023) 120781

Fig. 5. Observed and simulated (using the single-porosity model [SPM] and dual-porosity model [DPM]) soil water contents at different depths for the three soil
profiles (LS, MS, HS).
Table 6
Travel times and average velocities of bromide front from the soil surface to different soil depths.
Term

Depth (cm)

LS
SPM

DPM

Relative difference%

Travel time (day)

20
60
100
200
250 (300,275)
500 (450,
430)
700
20
60
100
200
250 (300,275)
500 (450,
430)
700

1.09
1.38
1.69
2.55
2.71

1.08
1.38
1.51
1.90
2.33

4.97

Bromide front velocity (cm/
day)

MS
SPM

DPM

Relative difference
%

HS
SPM

DPM

Relative difference
%

− 0.9
0.0
− 10.7− 25.5− 14.0

1.12
1.42
1.61
2.46
3.46

1.12
1.42
1.62
2.44
3.46

0.0
0.0
0.6
− 0.8
0.0

1.19
1.37
1.57
2.33
2.52

1.17
1.33
1.47
1.90
1.97

− 1.7
− 2.9
− 6.4
− 18.5− 21.83.09− 37.8

5.10

5.08

− 0.4

3.23

2.50

− 22.6

7.86
18.35
43.48
59.17
78.43
92.25

5.17
18.52
43.48
66.23
105.26107.30− 34.2
0.9
0.0
11.9
34.2
16.3

7.91
17.86
42.25
62.11
81.30
86.71

7.90
17.86
42.25
61.73
81.97
86.71

− 0.1
0.0
0.0
− 0.6
0.8
0.0

4.43
16.81
43.80
63.69
85.84
109.13

3.60
17.09
45.11
68.03
105.26139.59− 18.7
1.7
3.0
6.8
22.6
27.9

100.60161.81

60.8

88.24

88.58

0.4

133.13172.00

29.2

89.06

135.40

52.0

88.50

88.61

0.1

158.01

194.4423.1

aeration, and soil and groundwater salt leaching were also analyzed(Texts S2~S4), which were overall comparable to those in Table 6.
Darcy’s law calculates water fluxes through the entire cross-sectionalarea, but water flow occurs only in soil pores. Therefore, the pore watervelocity v (or bromide front velocity when considering only convective
bromide transport) is related to Darcy flux Jw by soil water content θ(Radcliffe and Šimůnek, 2018):
v = Jw /θ

and 0.022, 0.014, and 0.019 cm/min for DPM at LS, MS, and HS,
respectively. These values were overall consistent with those in Section
3.3, despite some differences probably due to dispersive or diffusive
bromide transport in this study.
4. Discussion
4.1. Impact of preferential flow on model performance(16)

The lag between observed and simulated wetting and/or bromide
fronts in the deeper profiles of LS and HS when using the single porosity
model (Figs. 5-6) clearly indicates the existence of preferential flow.
Preferential flow likely occurred due to the combined effects of intense
infiltration, dry climate (initially dry soil), soil texture heterogeneity,presence of macrofauna (e.g., earthworms), and active and decaying
crop roots. Continuous water ponding at the soil surface resulting from

During Ag-MAR, the thickness-weighted water contents at the top
three layers θTop were 0.303, 0.248, and 0.261 cm3/cm3 at LS, MS, andHS, respectively. The bromide front velocities at 200 cm v200 cm were
78.43, 81.30, and 85.84 cm/day for SPM, and 105.26, 81.97, and105.26 cm/day for DPM at LS, MS, and HS, respectively (Table 6). The
corresponding Darcy fluxes for groundwater recharge GR calculatedusing Eq. (16) were therefore 0.017, 0.014, and 0.016 cm/min for SPM,
9

T. Zhou et al.

Water Research 247 (2023) 120781

Fig. 6. Observed and simulated (using the single-porosity model [SPM] and dual-porosity model [DPM]) bromide concentrations at different depths for the three soil
profiles (LS, MS, HS).

the application of about 8–9 m of water in a month (Table 1) is more
likely to produce preferential flow than intermittent flooding fromnatural precipitation (Chen et al., 2002; Mitchell and van Genuchten,
1993; Selker et al., 1995). The semi-arid climate and particularly therain-free summer likely promoted the formation of macropores, espe­
cially in the top soil, which had a higher clay content and is therefore
more prone to shrink-swell dynamics that can create desiccation cracks(Jiang et al., 2010). In addition, the almond orchard was just fallowed,
with almond wood chips incorporated into the soil to about 50 cm
resulting in higher hydraulic conductivities during the experimental
period (Fig. 7a(1)). Research has shown that the infiltration rate may
increase in soils with decaying plant roots or with wormholes serving aspreferential flow paths (Fig. 7a(3) and 7a(4)) (Mitchell et al., 1995).
Such features could explain preferential flow in the top 0–66 cm of LS
and HS (Fig. 2). The presence of a coarser-textured soil layer (near thebottom of the profiles) overlain by a fine-textured soil layer may produce
funnel or fingered flow (Fig. 7a(5) and 7a(6)) (Council, 2001; Wanget al., 2018), resulting in preferential flow in the deep layers of profile LS
and HS (Figs. 5-6).
There were only small differences between the single- and dualporosity models in the simulated surface ponding levels and soil watercontents (see the model performances in Table 4). However, the dualporosity model (considering preferential flow) produced much betterfits for the bromide BTC values and trends observed at LS (Fig. 6).
However, both SPM and DPM could not capture the observed bromide
BTCs very well (Table 4). This may be associated with observation errors
since the bromide samples were not taken at the exact same locations as
the soil sensors (at a horizontal distance of 1.12 m). It might also indi­cate that some other hydrological processes may be occurring thatcannot be described by SPM and DPM (discussed in Section 4.2). At HS,
the dual-porosity model captured overall trends better (increased R2),
but it did not capture observed values (increased NRMSEs) as well as the
single porosity model. The dual-porosity model simulated much higher
peak values of the bromide BTCs than were observed at HS. This may berelated to the fact that the temporal resolution of the bromide samples
taken in the field was at minimum 4 h and hence may not have accu­
rately captured the real peak values that the dual-porosity model sug­gested. Finer spatial and temporal resolutions of field measurements
would better constrain model parameters and improve model

performance.
4.2. Impacts of flow dimensionality on model performance
Measured soil water contents showed abrupt decreases (before May16) and/or increases (after May 16) in all three profiles at depths of 60
and 100 cm that cannot be captured by the model (Fig. 5). There wasalso strong tailing in measured BTCs at shallow depths of LS (Fig. 6),
indicating significant initial tracer storage and slow re-release at theseshallow depths. In addition, measured steady water contents at layers
below 200 cm were very low, e.g., about 0.1 cm3/cm3 at 250 cm of LS,
and about 0.1 cm3/cm3 at 275 and 430 cm of HS, far smaller than theirsaturated water contents (Table 2). The O2 concentrations at these
depths were still high during flooding (Fig. S5a, S6a, S7a), indicating thesoil was mostly unsaturated at those depths.
As shown in Tables 2-3, all three profiles had a low conductivity layer(silty clay or silty clay loam, as shown in Fig. 2) at 67–200 cm depth(Ks=0.015–0.017 cm/min) followed by a higher conductivity layer(Ks=0.4–16 cm/min in LS, HS) at 201–400 cm depth. This lithologicalcombination may form the capillary barrier, where water flux to thedeeper profiles was limited and likely only occurred as finger flow, whilelateral flow within the low-conductivity layers (i.e., interflow H2 in
Fig. 7) likely dominated (Ho and Webb, 1998). As a result, water andsolutes could have been perched at these less permeable layers, allowing
less water and solutes to move downwards to deeper layers until acritical soil water potential is reached (Si et al., 2011), explaining thelow, unsaturated water contents at soil depths below 200 cm. Thisphenomenon was also studied and discussed in previous studies. For
example, (Botros et al., 2012; Harter and Yeh, 1996) demonstrated thelateral spreading of solutes due to the heterogeneous unsaturated zone
leads to extensive tailing in the observed breakthrough curves.
The three profiles had almond wood chips incorporated into the top
soil layer, following removal of a 20-year-old orchard, creating flow
pathways with much higher soil hydraulic conductivities within the top
30 cm (Fig. 7a(1)). This is likely the reason that estimated soil hydraulicconductivities for the silty clay layers were higher than their typicalrepresentative values (Tables 2-3 and Text S1). This much higher satu­
rated hydraulic conductivity of the top soil layer, compared to the un­derlying clay soils, likely promoted lateral water flow (i.e., interflow H1
10

T. Zhou et al.

Water Research 247 (2023) 120781

Fig. 7. Conceptual models of water flow and bromide transport (a), water and bromide mass balance (b), and bromide travel times (c) during Ag-MAR in the study
profiles. Blue polygons and orange circles represent water and bromide molecules, respectively. Blue and orange arrows represent water flow and bromide transport
directions, respectively. I: Bromide application; F: Flooding; L: Water loss outside the berms; Ea : Evaporation; H1: Horizontal flow through wood chips; H2: horizontal
flow caused by capillary barrier; D: Deep drainage; C: Capillary rise; G: Groundwater flow. (1)~(6) are possible preferential flow mechanisms caused by (1) woodchips, (2) soil repellency, (3) wormhole, (4) decayed roots, (5) impeding layers, and (6) layer boundaries. Γw and Γs represent the water and solute transfer terms in= 0, and Γs ∕=0 at LS (preferential flow and nonequilibrium bromide transport); Γw ∕= 0, and Γs =0 at HS (preferential
Eqs. (10) and 13, respectively. In this study, Γw ∕flow and equilibrium bromide transport); Γw = 0, and Γs ∕=0 at MS (uniform water flow and nonequilibrium bromide transport). SPM and DPM represent single anddual-porosity model in HYDRUS, respectively.

significantly differ from the site-specific real values. On the other hand,the simultaneous fitting of 4 (single-porosity model) or up to 8 param­
eters (dual-porosity model) for each soil layer is likely to result in nonunique and local optimal parameter sets (Hopmans and Šimůnek, 1999).
Global parameter optimization algorithms may help improve this aspect(Zhou et al., 2022). However, since a single HYDRUS model execution
requires up to 50 s in this study, global optimization may face a very
high computational cost. Alternatively, we can identify highly corre­
lated parameters by correlation matrices (Figs. S10-S15) and fix some ofthem in future parameter optimization, thus alleviating the computa­
tional burden.

in Fig. 7a). The relatively lower-than-expected saturated water content
of third layer, which may be the outcome of compaction from agricul­
tural operations, further contributes to increased lateral flow within the
top soil layers.
Owing to the occurrence of lateral flow, the wetted unsaturated zone
is likely broader than the size of the ponding basins. In the 1D effective
parameter model representation of this system, the effective, fitted
saturated water content for the silty clay loam/silt loam layer at
67~200 cm is therefore lower than their typical values (Text S1). The
fitted, high effective saturated hydraulic conductivities (about 0.4–16
cm/min) at 201~400 cm (Tables 2-3) allow for the speed up of thisexcess downward water flow, which resulted in a better fit in watercontents at these depths, but it also promoted bromide leaching and thus
caused lower simulated bromide concentrations than observed at the
deeper depths (e.g., 250 and 500 cm of LS; 300 and 450 cm of MS,
Fig. 6). Based on the analysis above, the conceptual models of water flow
and bromide transport in the three profiles can be deduced (Fig. 7a).

4.4. Suitability of implementing ag-mar
Because of their close proximity, the three profiles had the same land
use and hydroclimatological conditions but differed in subsurface hy­
drogeology. This study can therefore provide some insights into thefield-scale variability that one can expect when implementing Ag-MAR
at the field scale.
While HS provided the largest recharge efficiency compared to MS
and LS, recharge efficiency between sites varied only between 88 % and
90 %, because of similar effective saturated hydraulic conductivities atlayers above 200 cm as discussed in Sections 3.3 and 3.4. However, thedegree of preferential flow varied distinctly between the three profiles,
with HS showing the largest degree of preferential flow, and MS the leastindication of preferential flow, and LS being in between the two(Fig. S16). Similarly, travel times (flow velocities) were longest (slow­est) at MS, followed by LS, and were shortest (fastest) at HS (Table 6,
Fig. 7c). Bromide transport velocities differed by as much as 119 %between the three sites. This can also be verified by the dynamics ingroundwater table depth and EC (Fig. S8). For example, there wereabrupt decreases in groundwater depth and increases in EC after the
beginning of flooding at LS and HS, which reached steady rates after a
few days. In addition, the peaks in groundwater EC (about 2000 µS/cm
at LS and 3500 µS/cm at HS) were much higher than pore water ECs in

4.3. Other possible reasons behind model deficiencies
First, as discussed in Section 4.2, there were abrupt decreases inwater contents at 60 and 100 cm at LS, at 100 cm at MS, and at 60 and
100 cm at HS at the beginning of the flooding (before May 17) that could
not be captured by the model (Fig. 5). This could be due to some sub­surface heterogeneities (discussed in Section 4.2) or measurement errors
such as overshoot during saturation increase (Xiong, 2014). Water
repellency, which may develop in long-term dry soils and reduce the soil
infiltration capacity (and take some time to overcome), may be anotherexplanation for these differences (Fig. 7a(2)) (Doerr et al., 2000).
HYDRUS cannot simulate flow in repellent soils, which may result insimulation errors.
Second, detailed measurements of soil particle distributions andsaturated hydraulic conductivities were also unavailable. Therefore,
parameter optimization runs started with the average soil hydraulicparameters of different soil textures in this region (Table S1), which may
11T. Zhou et al.

Water Research 247 (2023) 120781

the bottom part (about 300~500 cm) of the soil profiles before flooding(about 1450~1650 µS/cm at LS, and 1050 µS/cm at HS (Fig. S5d,
Fig. S7d). This suggests that preferential flow transported salts from thesoil surface layers (with much higher pore-water ECs) to groundwater
while bypassing the soil matrix in the upper parts of the soil profiles. In
contrast, the groundwater table rise and EC dynamics were always more
subtle at MS, representing a slower flow rate in the form of piston flow.
The occurrence of preferential flow helps accelerate the timing of soil
salt leaching at HS and LS, but also poses a greater risk for microbes orcontaminants to be transported to groundwater since it reduces the time
available for chemical or pathogen immobilization or degradation(Willkommen et al., 2021). However, evaluating the amount of anon-conservative pollutant transported to groundwater because ofpreferential flow needs further modeling studies.
Overall, all three profiles were able to achieve similar groundwater
recharge efficiencies under the tested flooding regime for Ag-MAR(Table 1) considering their varying soil textures (Fig. 2). However, thewater flow and solute transport processes might be very distinct. The
suitability of implementing Ag-MAR depends on specific needs. In
practical applications, the vadose zone with higher sand contents (suchas HS) may imply more preferential flow (Sendros et al., 2020), which
promotes more focused soil salt leaching, while a vadose zone with moresilt contents (such as MS) would likely have a more muted contaminanttransport response.

Acknowledgement
This research was supported by the Gordon and Betty Moore Foun­dation (grant no. 7975), the USDA National Institute of Food andAgriculture (grant no. CA-D-LAW-2513-H), USDA National Institute ofFood and Agriculture (grant no. 2021–68012–35914), and USDA Nat­ural Resources Conservation Agreement NR183A750023C005. We also
acknowledge the editors and two anonymous reviewers for theirconstructive comments on this manuscript.
Supplementary materials
Supplementary material associated with this article can be found, inthe online version, at doi:10.1016/j.watres.2023.120781.
References
Alam, S., Borthakur, A., Ravi, S., Gebremichael, M., Mohanty, S.K., 2021. Managed
aquifer recharge implementation criteria to achieve water sustainability. Sci. Total
Environ. 768, 19. https://doi.org/10.1016/j.scitotenv.2021.144992.
Arnaud, E., Best, A., Parker, B.L., Aravena, R., Dunfield, K., 2015. Transport ofEscherichia coli through a thick vadose Zone. J. Environ. Qual. 44 (5), 1424–1434.
https://doi.org/10.2134/jeq2015.02.0067.
Bachand, P.A.M., Roy, S.B., Choperena, J., Cameron, D., Horwath, W.R., 2014.
Implications of using on-farm flood flow capture to recharge groundwater andmitigate flood risks along the Kings River, CA. Environ. Sci. Technol 48 (23),13601–13609. https://doi.org/10.1021/es501115c.
Bali, K.M., Mohamed, A.Z., Begna, S., Wang, D., Putnam, D., Dahlke, H.E., Eltarabily, M.
G., 2023. The use of HYDRUS-2D to simulate intermittent Agricultural Managed
Aquifer recharge (Ag-MAR) in Alfalfa in the San Joaquin valley. Agric. Water
Manage. 282, 15. https://doi.org/10.1016/j.agwat.2023.108296.
Botros, F.E., Onsoy, Y.S., Ginn, T.R., Harter, T., 2012. Richards equation-based modelingto estimate flow and nitrate transport in a deep alluvial vadose zone. Vadose Zone J.
11 (4), 16. https://doi.org/10.2136/vzj2011.0145.
Bradford, S.A., Leij, F.J., Schijven, J., Torkzaban, S., 2017. Critical role of preferential
flow in field-scale pathogen transport and retention. Vadose Zone J. 16 (4), 13.
https://doi.org/10.2136/vzj2016.12.0127.
Chen, C.C., Roseberg, R.J., Selker, J.S., 2002. Using microsprinkler irrigation to reduceleaching in a shrink/swell clay soil. Agric. Water Manage. 54 (2), 159–171. https://doi.org/10.1016/s0378-3774(01)00150-0.
Clapp, R.B., Hornberger, G.M., 1978. Empirical equations for some soil hydraulicproperties. Water Resour. Res. 14 (4), 601–604. https://doi.org/10.1029/
WR014i004p00601.
Council, N.R., 2001. Conceptual Models of Flow and Transport in the Fractured Vadose
Zone. National Academies Press.
Dahlke, H.E., LaHue, G.T., Mautner, M.R., Murphy, N.P., Patterson, N.K.,
Waterhouse, H., Yang, F., Foglia, L., 2018. Managed aquifer recharge as a tool toenhance sustainable groundwater management in California: examples from field
and modeling studies. Advances in Chemical pollution, Environmental Management
and Protection. Elsevier, pp. 215–275.
de Vries, J.J., Simmers, I., 2002. Groundwater recharge: an overview of processes andchallenges. Hydrogeol. J. 10 (1), 5–17. https://doi.org/10.1007/s10040-001-01717.
Doerr, S.H., Shakesby, R.A., Walsh, R.P.D., 2000. Soil water repellency: its causes,
characteristics and hydro-geomorphological significance. Earth-Sci. Rev. 51 (1–4),
33–65. https://doi.org/10.1016/s0012-8252(00)00011-8.
Ganot, Y., Dahlke, H.E., 2021. A model for estimating Ag-MAR flooding duration basedon crop tolerance, root depth, and soil texture data. Agric. Water Manage. 255, 12.
https://doi.org/10.1016/j.agwat.2021.107031.
Gelhar, L.W., Welty, C., Rehfeldt, K.R., 1992. A critical-review of data on field-scaledispersion in aquifers. Water Resour. Res. 28 (7), 1955–1974. https://doi.org/
10.1029/92wr00607.
Gerke, H.H., van Genuchten, M.T., 1993. Evaluation of a 1st-order water transfer termfor variably saturated dual-porosity flow models. Water Resour. Res. 29 (4),1225–1238. https://doi.org/10.1029/92wr02467.
Guo, Z., Fogg, G.E., Chen, K., Pauloo, R., Zheng, C., 2023. Sustainability of regional
groundwater quality in response to managed aquifer recharge. Water Resour Res 59(1), e2021WR031459.
Gurevich, H., Baram, S., Harter, T., 2021. Measuring nitrate leaching across the criticalzone at the field to farm scale. Vadose Zone J. 20 (2), 16. https://doi.org/10.1002/
vzj2.20094.
Harter, T., Yeh, T.C.J., 1996. Conditional stochastic analysis of solute transport inheterogeneous, variably saturated soils. Water. Resour. Res. 32 (6), 1597–1609.
https://doi.org/10.1029/96wr00503.
Haws, N.W., Rao, P.S.C., Šimůnek, J., Poyer, I.C., 2005. Single-porosity and dual-porosity
modeling of water flow and solute transport in subsurface-drained fields using
effective field-scale parameters. J. Hydrol. (Amst.) 313 (3–4), 257–273. https://doi.
org/10.1016/j.jhydrol.2005.03.035.
Ho, C.K., Webb, S.W., 1998. Capillary barrier performance in heterogeneous porous
media. Water Resour. Res. 34 (4), 603–609. https://doi.org/10.1029/98wr00217.

5. Conclusions
Our modeling results show that the dual-porosity models (consid­
ering preferential flow) can better fit the arrival times of bromide fronts
but cannot significantly improve the overall model performance. Pref­
erential flow occurred due to the combined effects of dry antecedent soil
moisture followed by flooding, dry climate, soil texture, and the incor­
poration of almond wood chips into the topsoil, etc.
Preferential flow did not significantly impact the water balancecalculations (within 2 %), but it decreased the travel times of bromide
from the soil surface to different depths of the soil profiles by up to 38 %,compared to the predictions provided by the single-porosity model.
In terms of groundwater recharge potential, HS showed a higher
efficiency than MS or LS, but the differences were relatively minor(within 2 %). LS showed the highest degree of preferential flow, fol­
lowed by HS and MS, and the overall average bromide transport ve­
locities differed by up to 119 %. In brief, similar recharge efficiency canbe achieved at sites with differing soil texture profiles but subsurface
heterogeneity can have substantial effects on salt and contaminanttransport dynamics, which should be considered when implementing
Ag-MAR.
The potential occurrence of lateral interflow is another importantreason behind the model deficiency and may lead to errors in the waterbalance calculation in our relatively small experimental plots. In addi­tion, we focused mainly on the effects of soil textural differences (i.e.,
between-lithofacies or large-scale heterogeneity) on Ag-MAR recharge
and neglected the impact of horizontal heterogeneities within lithofacies(small-scale heterogeneity). Future work should extend the current 1Dmodeling analysis to 2D/3D to get full insight into soil heterogeneity’s
impacts (especially within-lithofacies or small-scale heterogeneity) on
Ag-MAR recharge.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Data availability
Data will be made available on request.
12

T. Zhou et al.

Water Research 247 (2023) 120781
Selker, J., Cao, W.D., Roseberg, R., 1995. Use of ultra-low rate application devices toeliminate macropore flow during irrigation. In: 5th International Microirrigation
Congress - Microirrigation for a Changing World: Conserving Resources/Preserving
the Environment, 95. Amer Soc Agricultural Engineers, Orlando, Fl, pp. 54–59.
Sendros, A., Himi, M., Lovera, R., Rivero, L., Garcia-Artigas, R., Urruela, A., Casas, A.,
2020. Electrical resistivity tomography monitoring of two managed aquifer recharge
ponds in the alluvial aquifer of the Llobregat River (Barcelona, Spain). Near Surf.
Geophys. 18 (4), 353–368. https://doi.org/10.1002/nsg.12113.
Si, B., Dyck, M., Parkin, G.W., 2011. Flow and transport in layered soils PREFACE. Can. J.
Soil Sci. 91 (2), 127–132. https://doi.org/10.4141/cjss11501.
Siebert, S., Burke, J., Faures, J.M., Frenken, K., Hoogeveen, J., Doll, P., Portmann, F.T.,
2010. Groundwater use for irrigation - a global inventory. Hydrol. Earth Syst. Sci. 14(10), 1863–1880. https://doi.org/10.5194/hess-14-1863-2010.
Šimůnek, J., Hopmans, J.W., 2002. 1.7 parameter optimization and nonlinear fitting.
Method. Soil Anal.: Part 4 Phys. Methods 5, 139–157.
Šimůnek, J., Jarvis, N.J., van Genuchten, M.T., Gardenas, A., 2003. Review andcomparison of models for describing non-equilibrium and preferential flow andtransport in the vadose zone. J. Hydrol. (Amst.) 272 (1–4), 14–35.
Šimůnek, J., van Genuchten, M.T., 2008. Modeling nonequilibrium flow and transport
processes using HYDRUS. Vadose Zone J. 7 (2), 782–797. https://doi.org/10.2136/
vzj2007.0074.
Šimůnek, J., van Genuchten, M.T., Sejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15 (7), 25.
https://doi.org/10.2136/vzj2016.04.0033.
Šimůnek, J., Wendroth, O., Wypler, N., van Genuchten, M.T., 2001. Non-equilibrium
water flow characterized by means of upward infiltration experiments. Eur. J. Soil
Sci. 52 (1), 13–24. https://doi.org/10.1046/j.1365-2389.2001.00361.x.
Sprenger, C., Hartog, N., Hernandez, M., Vilanova, E., Grutzmacher, G., Scheibler, F.,
Hannappel, S., 2017. Inventory of managed aquifer recharge sites in Europe:
historical development, current situation and perspectives. Hydrogeol. J. 25 (6),1909–1922. https://doi.org/10.1007/s10040-017-1554-8.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. https://doi.
org/10.2136/sssaj1980.03615995004400050002x.
Wang, C.Z., Liu, G., McNew, C.P., Volkmann, T.H.M., Pangle, L., Troch, P.A., Lyon, S.W.,
Kim, M., Huo, Z.L., Dahlke, H.E., 2022. Simulation of experimental synthetic DNA
tracer transport through the vadose zone. Water Res. 223, 10. https://doi.org/
10.1016/j.watres.2022.119009.
Wang, C.Z., Wang, R.Y., Huo, Z.L., Xie, E., Dahlke, H.E., 2020. Colloid transport throughsoil and other porous media under transient flow conditions-a review. Wiley
Interdiscip. Rev.-Water 7 (4), 33. https://doi.org/10.1002/wat2.1439.
Wang, Y.C., Li, Y., Wang, X.F., Chau, H.W., 2018. Finger flow development in layeredwater-repellent soils. Vadose Zone J. 17 (1), 11. https://doi.org/10.2136/
vzj2017.09.0171.
Waterhouse, H., Bachand, S., Mountjoy, D., Choperena, J., Bachand, P.A.M., Dahlke, H.
E., Horwath, W.R., 2020. Agricultural managed aquifer recharge - water quality
factors to consider. Calif. Agric. 74 (3), 144–154. https://doi.org/10.3733/
ca.2020a0020.
Wendt, D.E., Van Loon, A.F., Scanlon, B.R., Hannah, D.M., 2021. Managed aquifer
recharge as a drought mitigation strategy in heavily-stressed aquifers. Environ. Res.
Lett. 16 (1), 13. https://doi.org/10.1088/1748-9326/abcfe1.
Willkommen, S., Lange, J., Ulrich, U., Pfannerstill, M., Fohrer, N., 2021. Field insights
into leaching and transformation of pesticides and fluorescent tracers in agriculturalsoil. Sci. Total Environ. 751, 12. https://doi.org/10.1016/j.scitotenv.2020.141658.
WWDR, 2022. Groundwater: making the invisible visible. The United Nations World
Water Development Report.
Xiong, Y.W., 2014. Flow of water in porous media with saturation overshoot: a review.
J. Hydrol. (Amst.) 510, 353–362. https://doi.org/10.1016/j.jhydrol.2013.12.043.
Zhou, T., Šimůnek, J., Braud, I., 2021. Adapting HYDRUS-1D to simulate the transport ofsoil water isotopes with evaporation fractionation. Environ. Modell. Softw. 143,
105118 https://doi.org/10.1016/j.envsoft.2021.105118.
Zhou, T., Šimůnek, J., Braud, I., Nasta, P., Brunetti, G., Liu, Y., 2022. The impact of
evaporation fractionation on the inverse estimation of soil hydraulic and isotope
transport parameters. J. Hydrol. (Amst.) 612, 128100. https://doi.org/10.1016/j.
jhydrol.2022.128100.
Zhou, T., Šimůnek, J., Nasta, P., Brunetti, G., Gaj, M., Neukum, C., Post, V., 2023. The
impact of soil tension on isotope fractionation, transport, and interpretations of theroot water uptake origin. Water Resour. Res. 59, e2022WR034023 https://doi.org/
10.1029/2022WR034023.

Hopmans, J., Šimůnek, J., 1999. Review of inverse estimation of soil hydraulicproperties. In: van Genuchten, M.Th., Leij, F.J., Wu, L. (Eds.), Characterization andMeasurement of the Hydraulic Properties of Unsaturated Porous Media. University ofCalifornia, Riverside, CA, pp. 643–659.
Imig, A., Augustin, L., Groh, J., Putz, T., Zhou, T., Einsiedl, F., Rein, A., 2023. Fate ofherbicides in cropped lysimeters: 1. Influence of different processes and model
structure on vadose zone flow. Vadose Zone J. 22 (5), 13. https://doi.org/10.1002/
vzj2.20265.
Isch, A., Montenach, D., Hammel, F., Ackerer, P., Coquet, Y., 2019. A comparative studyof water and bromide transport in a bare loam soil using lysimeters and field plots.
Water (Basel) 11 (6), 25. https://doi.org/10.3390/w11061199.
Jiang, S., Pang, L.P., Buchan, G.D., Šimůnek, J., Noonan, M.J., Close, M.E., 2010.
Modeling water flow and bacterial transport in undisturbed lysimeters underirrigations of dairy shed effluent and water using HYDRUS-1D. Water Res. 44 (4),1050–1061. https://doi.org/10.1016/j.watres.2009.08.039.
Kocis, T.N., Dahlke, H.E., 2017. Availability of high-magnitude streamflow forgroundwater banking in the Central Valley, California. Environ. Res. Lett. 12 (8), 13.
https://doi.org/10.1088/1748-9326/aa7b1b.
Köhne, J.M., Köhne, S., Mohanty, B.P., Šimůnek, J., 2004. Inverse mobile-immobile
modeling of transport during transient flow: effects of between-domain transfer andinitial water content. Vadose Zone J. 3 (4), 1309–1321. https://doi.org/10.2113/
3.4.1309.
Kourakos, G., Dahlke, H.E., Harter, T., 2019. Increasing Groundwater availability andseasonal base flow through agricultural managed aquifer recharge in an irrigatedbasin. Water Resour Res 55 (9), 7464–7492. https://doi.org/10.1029/
2018wr024019.
Levintal, E., Huang, L., García, C.P., Coyotl, A., Fidelibus, M.W., Horwath, W.R.,
Rodrigues, J.L.M., Dahlke, H.E., 2023a. Nitrogen fate during agricultural managed
aquifer recharge: linking plant response, hydrologic, and geochemical processes. Sci.
Total Environ. 864, 161206.
Levintal, E., Kniffin, M.L., Ganot, Y., Marwaha, N., Murphy, N.P., Dahlke, H.E., 2023b.
Agricultural managed aquifer recharge (Ag-MAR)-a method for sustainablegroundwater management: a review. Crit. Rev. Environ. Sci. Technol. 53 (3),291–314. https://doi.org/10.1080/10643389.2022.2050160.
Li, E., Shanholtz, V., Carson, E., 1976. Estimating Saturated Hydraulic Conductivity andCapillary Potential at the Wetting Front, Dep. of Agr. Eng. Va. Polytech. Inst. andState Univ., Blacksburg.
Marwaha, N., Kourakos, G., Levintal, E., Dahlke, H.E., 2021. Identifying agricultural
managed aquifer recharge locations to benefit drinking water supply in ruralcommunities. Water. Resour. Res. 57 (3), 24. https://doi.org/10.1029/
2020wr028811.
Mitchell, A.R., Ellsworth, T.R., Meek, B.D., 1995. Effect of root systems on preferential
flow in swelling soil. Commun. Soil Sci. Plant Anal. 26 (15–16), 2655–2666. https://doi.org/10.1080/00103629509369475.
Mitchell, A.R., van Genuchten, M.T., 1993. Flood irrigation of a cracked soil. Soil Sci.
Soc. Am. J. 57 (2), 490–497. https://doi.org/10.2136/
sssaj1993.03615995005700020032x.
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.
Murphy, N.P., Waterhouse, H., Dahlke, H.E., 2021. Influence of agricultural managed
aquifer recharge on nitrate transport: the role of soil texture and flooding frequency.
Vadose Zone J. 20 (5), 16. https://doi.org/10.1002/vzj2.20150.
Nimmo, J.R., Perkins, K.S., Plampin, M.R., Walvoord, M.A., Ebel, B.A., Mirus, B.B., 2021.
Rapid-response unsaturated zone hydrology: small-scale data, small-scale theory, big
problems. Front. Earth Sci. 9, 7. https://doi.org/10.3389/feart.2021.613564.
Perzan, Z., Osterman, G., Maher, K., 2023. Controls on flood managed aquifer recharge
through a heterogeneous vadose zone: hydrologic modeling at a site characterizedwith surface geophysics. Hydrol. Earth Syst. Sci. 27 (5), 969–990. https://doi.org/
10.5194/hess-27-969-2023.
Radcliffe, D.E., Šimůnek, J., 2018. Soil Physics with HYDRUS: Modeling andApplications. CRC Press.
Sasidharan, S., Bradford, S.A., Šimůnek, J., Kraemer, S.R., 2021. Virus transport fromdrywells under constant head conditions: a modeling study. Water Res. 197, 14.
https://doi.org/10.1016/j.watres.2021.117040.
Scanlon, B.R., Faunt, C.C., Longuevergne, L., Reedy, R.C., Alley, W.M., McGuire, V.L.,
McMahon, P.B., 2012. Groundwater depletion and sustainability of irrigation in the
US High plains and central valley. Proc. Natl. Acad. Sci. USA. 109 (24), 9320–9325.
https://doi.org/10.1073/pnas.1200311109.

13