Water Research 267 (2024) 122442

Contents lists available at ScienceDirect

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

The impact of managed aquifer recharge on the fate and transport ofpesticides in agricultural soils
Tiantian Zhou a,* , Nels Ruud b , Jiří Šimůnek c , Giuseppe Brunetti d , Elad Levintal e,
Cristina Prieto García a , Helen E. Dahlke a,*a

Department of Land, Air and Water Resources, University of California, Davis, CA, 95616, USA
Department of Pesticide Regulation, California Environmental Protection Agency, Sacramento, CA, USA
Department of Environmental Sciences, University of California, Riverside, CA, 92521, USAd
Department of Civil Engineering, University of Calabria, Rende, Italye
Zuckerberg Institute for Water Research, The Jacob Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev, Sde Boker campus, 84990, Israelbc

A R T I C L E I N F O

A B S T R A C T

Keywords:
HYDRUS
Particle swarm optimization
Preferential flow
Pesticide leaching
Capillary barrier

Groundwater aquifers worldwide experience unsustainable depletion, compounded by population growth, eco­nomic development, and climate forcing. Managed aquifer recharge provides one tool to alleviate flood risk andreplenish groundwater. However, concerns grow that intentional flooding of farmland for groundwater recharge,
a practice known as Ag-MAR, may increase the leaching of pesticides and other chemicals into groundwater. This
study employs a physically based unsaturated flow model to determine the fate and transport of residues of four
pesticide in three vadose zone profiles characterized by differing fractions of sand (41 %, 61 %, and 84 %) in
California’s Central Valley. Here, we show that the complex heterogeneity of alternating coarse and fine-grainhydrogeologic units controls the transit times of pesticides and their adsorption and degradation rates. Unsat­
urated zones that contain a higher fraction of sand are more prone to support preferential flow, higher recharge
rates (+8 %), and faster (42 %) water flow and pesticide transport, more flooding-induced pesticide leaching(about 22 %), as well as more salt leaching correlating with increased risks of groundwater contamination.
Interestingly, considering preferential flow predicted higher degradation and retention rates despite shorter
travel times, attributed to the trapping of pesticides in immobile zones where they degrade more effectively. The
findings underscore the importance of considering soil texture and structure in Ag-MAR practices to minimize
environmental risks while enhancing groundwater recharge. The study also highlights that selecting less mobile
pesticides can reduce leaching risks in sandy areas.

1. Introduction
California is leading the United States agricultural production in over
77 different products including fruit, vegetable and specialty crops suchas tomatoes (70 %), bell peppers (38 % of bell peppers), or carrots (85%) (Lucier and Jerardo, 2006). Imidacloprid, thiamethoxam, chloran­
traniliprole, and methoxyfenozide are common pesticides used to cultivate
many of these crops. Pesticides have played a substantial role in
improving crop yields, but also bear the risk of polluting freshwater
resources with potentially harmful impacts on human health whenleached into drinking water (Tudi et al., 2021). Diffuse pesticideleaching from agricultural fields is the dominant source of pesticides instreams and groundwater.

Pesticide management practices for runoff or groundwater protec­
tion areas have been adopted in the US since the 1980; mandating forexample minimum irrigation amounts (6 – 25 mm) and banded appli­
cations to reduce surface runoff or leaching risk to groundwater (PCPA,
1985). The main processes controlling the fate of pesticides in theenvironment are degradation, retention, volatilization, drift, atmo­
spheric dispersion, runoff, and leaching (Mottes et al., 2014). Once
pesticides enter the soil or vadose zone, their fate is controlled byadsorption to clay particles or organic matter, transformation anddegradation processes, soil water flow, and crop uptake (Köhne et al.,2009).
While pesticide fate and transport has been studied at soil column toplot to catchment scales for decades, understanding their fate in natural* Corresponding authors.
E-mail addresses: tzhou035@ucr.edu (T. Zhou), hdahlke@ucdavis.edu (H.E. Dahlke).
https://doi.org/10.1016/j.watres.2024.122442
Received 21 June 2024; Received in revised form 10 September 2024; Accepted 12 September 2024
Available online 14 September 20240043-1354/© 2024 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/).

T. Zhou et al.

Water Research 267 (2024) 122442

or structured soils has been a particular challenge since structured soils
can preferentially leach pesticides or other contaminants via nonequilibrium preferential flow paths (e.g. macropore flow) before they
can degrade or be adsorbed. Most agricultural soils are structured soils
that contain various soil horizons with contrasting soil hydraulic or
chemical properties. In addition, agricultural soils receiving pesticide
applications often undergo a variety of natural processes (e.g. wettingdrying cycles, shrink-swell behavior) and management practices (e.g.
tillage, plowing) that can facilitate preferential flow.
Because structured soils are prone to preferential flow, pesticidetransport models (e.g. MACRO, HYDRUS, PEARL, DAISY, CRACK-NP,
RZWQM, PESTLIS, SIMULAT) that explicitly consider preferential flow
have been shown to outperform models that do not (Köhne et al., 2009;
Scorza and Boesten, 2005). For example, Holback et al. (2022) used theagrohydrological model DAISY to simulate the transport of bentazone
and imidacloprid to drainpipes in a cracking clay field. After incorpo­
rating preferential flow features such as biopores and cracks into themodel, it simulated satisfactorily water flow and pesticide leaching todrain tiles (Holbak et al., 2022). Imig et al. (2023b) observed thetransport of four herbicides (metolachlor, terbuthylazine, prosulfuron,and nicosulfuron) in two lysimeters filled with a sandy gravel and clayey
sandy silt over 4.5 years and simulated transport in HYDRUS-1D. They
found that pesticide transport could be adequately described using
nonequilibrium (i.e. dual-porosity) models, where each model had adifferently sized mobile and immobile zones. For the clayey sandy silt
the relatively small mobile zone was dominating solute transport,
leading to higher solute concentration in the column drainage (Imiget al., 2023b). Dusek et al. (2015) simulated the fate and transport of five
pesticides (atrazine, imazaquin, sulfometuron methyl, S-metolachlor,and imidacloprid) in an undisturbed tropical Oxisol soil column exper­
iment by considering (or not) chemical nonequilibrium in HYDRUS-1D.
They showed that atrazine, sulfometuron methyl, and S-metolachlor
were better described considering kinetic sorption, while the other two
pesticides could be sufficiently characterized using the chemical equi­
librium model (Dusek et al., 2015). Sidoli et al. (2016) investigated thefate and transport of metolachlor and its two metabolites incolumn-leaching experiments in glaciofluvial soils, considering both
physical and chemical nonequilibrium in HYDRUS-1D. They identified astrong influence of chemical nonequilibrium on the fate of metolachlorbut not the metabolites (Sidoli et al., 2016). Köhne et al. (2006) applied
multiple physical and nonequilibrium models in HYDRUS-1D to simu­
late the fate and transport of isoproturon, terburazine, and bromide(Br− ) in aggregated loam columns subjected to multiple irrigation cy­
cles. They showed that simultaneous consideration of preferential flow
and kinetic adsorption (using the dual-permeability model with two-site
sorption) provided the best model performance (Köhne et al., 2006).
Many studies to date have focused on pesticide leaching in agricul­
tural soils under controlled conditions (e.g. soil column or lysimeter
experiments) or natural precipitation or irrigation regimes to evaluatesoil management and pesticide management practices for different soil
types. These studies have resulted in pesticide best management prac­
tices such as pesticides with medium or high mobility should not beapplied to saturated soils or prior to heavy irrigation (Waskom, 1995).
However, not many studies have focused on the fate and transport ofpesticides under large water applications (> 150 mm [6 inches]), incoarse textured soils (e.g. soils that support large percolation rates),considering preferential flow and kinetic transformation processes(Pang et al., 2000). The study of pesticide fate under large water ap­
plications is of particular interest, as many groundwater-dependent
agricultural regions experience growing groundwater overdraft (Perezet al., 2024), requiring managed aquifer recharge of depleted aquifers byusing agricultural lands for both, agricultural production and managed
aquifer recharge (MAR) (Dahlke et al., 2018; Levintal et al., 2023a,
2023b). Agricultural managed aquifer recharge (Ag-MAR), a method
where groundwater is actively replenished by spreading flood flows onto
agricultural lands, may exacerbate the leaching of pesticide residues

from the root zone to groundwater because of the large amounts of water(e.g. 0.5 to over 10 m of recharge) applied for MAR (Bachand et al.,2014; Guo et al., 2023; Levintal et al., 2023a; Murphy et al., 2021;
Waterhouse et al., 2020). Detailed knowledge of the fate and transportof pesticides in the vadose zone of agricultural soils when flooded for
MAR is needed to adequately inform the placing of suitable Ag-MAR
locations. Evidence that these systems or conditions require more
research can be gleaned from a few studies that observed early arrivalsof pesticides at the groundwater table in response to strong precipitation
or irrigation events soon after pesticides were applied in coarser, moreheterogeneous soils (Czapar et al., 1992; Stagnitti et al., 1995).
Addressing this gap helps us assess the feasibility of implementing
Ag-MAR while avoiding (or minimizing) groundwater pollution.
In this study, we combine field observations and non-equilibrium
flow modeling in HYDRUS-1D to describe the fate and transport offour pesticides (imidacloprid, thiamethoxam, chlorantraniliprole, andmethoxyfenozide) and bromide in the vadose zone of three Ag-MAR fieldsites in the Central Valley, California, USA. Our main research objectives
are 1) to capture and compare the fate and transport of the four pesti­
cides in three agricultural soils in response to one large (1.2 m) water
application using water balance and mass balance approaches; 2) eval­uate the role of vadose zone heterogeneity (e.g. soil hydraulic andchemical properties) on the occurrence of non-equilibrium preferential
flow at the three sites and its impact on the water and contaminant (e.g.,
pesticides, Br− ) transit times; and 3) analyze sensitivity of solutetransport and reaction parameters in HYDRUS, their calibration using
global optimization methods, and how they related to dominant pro­
cesses and factors governing the fate and transport of pesticides.
2. Methods
2.1. Study site and soil texture profiles
The Ag-MAR experiment was conducted at the Terranova Ranch(36◦ 34′27″N 120◦ 05′39″W, 50 m), located southwest of Fresno, CA, USA(Fig. 1a), within the Kings River basin, which is underlain by a pre­dominantly sand and gravel aquifer bound by the Corcoran Clay at a
depth of 140 m, a thick aquitard that spans the western half of the Kings
River basin. The depth of the groundwater table was approximately 70
m at the time of data acquisition (February 2021).
The soil at the site is a Traver fine sandy loam (fine-loamy, mixed,
thermic Natric Haploxeralf). Laboratory sediment analysis of undis­
turbed soil samples of the top 2.5 m of each soil profile (P1, P2, and P3)
showed progressively increasing fractions of sand from P1 (41 %) to P3(61 %) and further to P2 (84 %) (Fig. 7 and Table S1), with a cemented
duripan at around 1 m depth at P1 and P3 (Bachand et al., 2014).
2.2. Application of flood water, bromide tracer, and pesticides
A 32,376 m2 (8 acres) recharge plot (Fig. 1b) was flooded in February
16–24, 2021, using pumped groundwater. A total of 38,774.74 m3 ofwater (1.2 m in depth) was applied at a flow rate of ~3.35 m3/min. The
start times of flooding at each soil profile location were deduced fromthe time of water content increase, which was 12:50 on Feb. 16 at P3,
14:00 on Feb. 17 at P2, and 8:00 on Feb. 18 at P1. The end times offlooding at each soil profile were deduced from the endpoints ofdecreasing levels in surface ponding, which were 12:40 on Feb. 24 at P3,
19:50 on Feb. 24 at P2, and 19:50 on Feb. 24 at P1.
At each profile, 541 g of Br− (806 g of KBr) dissolved in 100 L ofwater were applied. The area of the tracer application was 2.5 × 7 m2
within each profile (Fig. 1c). The KBr solution was applied on February
15th, 2021, at P1 and P2, on February 16th, 2021, at P3, between
7:30–10am. Br− was thus applied at a concentration of 5410 ppm at an
irrigation rate of 0.00381 cm/min.
Table S2 shows the history of pesticide applications and their rates
for the four pesticides studied here (imidacloprid, thiamethoxam,2

T. Zhou et al.

Water Research 267 (2024) 122442

Fig. 1. A schematic of the study site in California, USA (a), the recharge plot, three soil profiles P1, P2, P3, and the potassium bromide (KBr) application location ineach profile (b), sensors and suction cups in each profile (c), and sampling details (d).

chlorantraniliprole, and methoxyfenozide). Table S3 provides the basic
physical and chemical properties of these pesticides.

uniform and can be described using the Richards andadvection-dispersion equations, respectively. DPM divides the soil porespace into mobile and immobile regions (i.e., considering preferential
flow/transport). Water flow or solute transport occurs only in the mobile
region, as described by the Richards and advection-dispersion equa­tions, respectively. At the same time, there can be water/solute transfer
between these two regions. The model setup, including input data, ini­
tial/boundary conditions and governing equations, is shown in Fig. 2
and Table S4. More details about the governing equations can be foundin Text S1 (Supporting Information).
The 250 cm soil profile was divided into five modeling layers ranging
from 0 to 46 cm, 47–76 cm, 77–122 cm, 123–182, and 183–250 cm bygrouping the original soil texture of the profile (Table S1) and associ­ating each model layer with one sensor depth (sensor installed at 0.2,
0.6, 1, 1.75 and 2.5 m). The simulation period was 104 days long, from0:00 on Dec. 17, 2020, to 24:00 on Mar. 31, 2021, which included preflooding, flooding (Feb. 16~Feb. 24, 2021), and post-flooding periods.
The spatial discretization was 1 cm throughout the soil profile, while thetemporal discretization was variable, with a minimum time step of 0.01
min.
The initial soil pressure head profile was obtained from measured
soil water contents at four soil depths (20, 60, 100, and 250 cm) and soil
water retention curves of typical soil texture classes (Radcliffe andŠimůnek, 2018) and then linearly interpolated between any two mea­
surement depths. The initial Br− concentrations were set to zerothroughout the soil profile, while the initial soil water pesticide con­centrations were prescribed based on field measurements.
For water flow, the upper boundary condition (BC) was set to an
atmospheric flux (with a maximum allowed surface water layer of 50

2.3. Field monitoring and data collection
The meteorological data, including precipitation (P) and potential
evaporation (Ep), were obtained from station 2 (Five Points) of the CIMIS
network (https://cimis.water.ca.gov/).
Sensors were installed at depths of 0.2, 0.6, 1.0, and 2.5 m at eachsoil profile. Sensors measuring ponding depth (CS-451, Campbell Sci­
entific, Logan, UT, USA), soil water content, electrical conductivity (EC),and soil temperature (TEROS-12, METER Inc., San Francisco, CA, USA),soil O2 (Figaro KE-25, Figaro Sensors, Rolling Meadows, IL, USA), andoxidation–reduction potential (ORP; built in house) were logged at a 10minute time interval (Fig. 1c-1d).
Breakthrough curve data for the KBr tracer and pesticides werecollected using high-volume dual-chamber suction cups (Model 1920
F1L24-B02M2, SoilMoisture Inc., Goleta, CA, USA), which were
installed at depths of 20, 60, 100, 175, and 250 cm and sampled at 7:00,
11:00, 15:00, 19:00, and 23:00 every day during flooding. The suction
cups were installed approximately 50 cm from one another within ahorizontal distance of 3.5 m from the sensor profile (Fig. 1c-1d).
2.4. HYDRUS-1D model setup
Water flow and the transport of potassium bromide (KBr) and pes­
ticides in the unsaturated zone were simulated using the single and dualporosity models (SPM and DPM) of the HYDRUS-1D software (Šimůnek
et al., 2016). SPM assumes that flow and transport processes in soil are3

T. Zhou et al.

Water Research 267 (2024) 122442

Fig. 2. Conceptual model setups for the (a) single-porosity and (b) dual-porosity models (SPM and DPM, respectively). Note that “W” and “S" represent water flow
and solute transport, respectively. The explanations of each variable are shown in Supplemental materials S1.

cm, i.e., the height of the berms). The potential water flux across the soil
surface is the difference between daily values of potential evaporation,
Ep, and precipitation (P) or flood irrigation (F). The lower BC at the soil
depth of 250 cm was set to free drainage because of the low water table
at the site (75 m). For solute transport, the upper BC was prescribed as asolute flux (i.e., a Cauchy BC), with bromide concentrations and irri­
gation fluxes during the bromide application as inputs. The model thenautomatically adjusts surface bromide concentrations depending on thedepth of the surface water level, evaporation, precipitation, irrigation
fluxes, and associated bromide concentrations. The lower BC was pre­scribed as a zero-concentration gradient (i.e., a Neumann BC when only
a convective solute flux occurs).

2.5. Parameter optimization and sensitivity analysis
In this study, we adopted a two-step optimization. Since bromidetransport is largely unaffected by adsorption/desorption and degrada­tion, measured breakthrough curves (BTCs) of bromide were used firstto estimate the soil hydraulic and basic (e.g., dispersivity) solutetransport parameters for both the single and dual-porosity models.
In the single-porosity models, θr was not optimized. Instead, thedefault values of corresponding soil textures were adopted first and thenmanually adjusted to obtain a better model fit. Therefore, five param­
eters (θs , α, n, Ks , and λ) were optimized for each layer (25 parameters intotal). The ranges of other parameters (Table S5) were prescribed as 1/
4

T. Zhou et al.

Water Research 267 (2024) 122442

5th and 6/5th of the minimum and maximum of their initial values forall soil textures in each modeling layer, obtained from the Rosetta
module in HYDRUS-1D, based on measured soil textural data (Table S1).
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; Šimůnek et al., 2001). Therefore, eight parameters (θmo,s , α,n, Ks , θim,r , θim,s , ωw , and λ) were optimized for each layer (40 parameters
in total). The ranges of α, n, Ks , and λ (Table S7) were the same as thosein the single-porosity models. The upper boundary of θmo,s and θim,s were
the same as ranges of θs in the single-porosity models, while the lower
boundary was set as 0.1 and 0.01, respectively. The ranges of θim,r andωw (Table S7) were set to 0~0.1 (Imig et al., 2023b) and 0~7e-3 min− 1(Brunetti et al., 2016), respectively.
In the second step, the parameters θs , α, n, Ks , and λ in the singleporosity model or θmo,s , α, n, Ks , θim,r , θim,s , ωw , and λ in the dualporosity model obtained during the first step were fixed, while thepesticide reaction parameters (i.e., adsorption and degradation) were
optimized based on measured pesticide BTCs. In this study, chemical
nonequilibrium was not considered. Therefore, two parameters (Kd andμ) were optimized for each layer (10 parameters in total) in the singleporosity model, while three parameters (Kd , μ, and ωs ) were optimized
for each layer (15 parameters in total) in the dual-porosity model.
Measurements of degradation or sorption parameters specific to thestudy site were unavailable. The ranges of Kd (Tables S6 and S8) wereprescribed as 0 and 6/5 times the maximum of Koc *TOC (i.e., the organiccarbon-water partition coefficient * the total organic carbon content) forall soil textures in each modeling layer (Tables S1 and S3). The ranges ofμ (Tables S6 and S8) were prescribed as 0 and 6/5 times the maximum of
ln2/t1/2 for all soil textures in each modeling layer (Table S3). The
ranges of ωs (Table S8) were set to 0~7e-3 min− 1 (Imig et al., 2023a).
The parameter optimization was performed using the particle swarm
optimization (PSO) algorithm using the PySwarm Library in Python. A
swarm of candidate solutions is moved around in the search space in the
PSO according to a few equations. The movement of the particles isguided by their optimal positions and that of the whole swarm. Once
improved positions are discovered, they guide the swarm’s movement.
This process is repeated until the global optimal position that all parti­
cles tend to follow is found (Shi and Eberhart, 1998). More details about
PSO can be found in (Brunetti et al., 2016) (Brunetti et al., 2017, 2022;
Zhou et al., 2022).
The PSO parameters (cognitive parameter c1 =− 0.267; social
parameter c2 =3.395; inertia-weight w=− 0.444) (Brunetti et al., 2016)
were used in this study. The number of swarm particles and iterations is200; thus, each parameter optimization had 40,000 runs. The Python
script produces an input parameter space, overwrites the input param­
eter file, and runs the executable module of HYDRUS-1D. For each PSO
run for the water flow and bromide transport parameters, theKling-Gupta efficiency (KGE) indices for the dynamics of surface pond­
ing levels (KGE_sp), soil water contents (KGE_wc), and bromide con­centrations(KGE_Br),andtheirweighted
average(KGE_avg=0.7*KGE_wc+0.2*KGE_sp+0.1*KGE_Br, where the weights
were determined based on the degree of trust in these three types ofdata) were calculated. For each PSO run for the pesticide transport andreaction parameters, the Kling-Gupta efficiency (KGE) indices forpesticide concentrations (i.e., KGE_Imidacloprid, KGE_Thiamethoxam,
KGE_Chlorantraniliprole, KGE_Methoxyfenozide) were calculated. The
KGE index compares the correlation coefficient (r), the ratio of meanvalues (β), and the ratio of variances (γ) between simulated andobserved data (Knoben et al., 2019). The value of the KGE index is al­ways smaller or equal to 1. The higher the KGE value, the better the fit
between the simulated and observed values.
KGE = 1 −[]0.5
(1 − r)2 + (1 − β)2 + (1 − γ)2

models) or the length of the modeled hydrograph was shorter than the
total simulation period (149, 760 min), the run was considered nonconvergent. The run was then terminated, and a large negative value(− 1E+7) was assigned to the objective function. Only the parameters
leading to the maximum KGE_avg were retained as optimized
parameters.
The Sobol’ global sensitivity analysis was conducted using theSensitivity Analysis Library (SALib) in Python to identify the mostinfluential transport and reaction parameters of pesticides (Text S2).
3. Results
3.1. Model parameters and performance
The HYDRUS-1D single-porosity (SP) and dual-porosity (DP) flow
and transport models were first fit to the observed volumetric watercontents, ponding depths, and KBr tracer breakthrough data before
pesticide adsorption and degradation parameters were optimized. Witha few exceptions, the optimized soil hydraulic parameters (Tables S4-S5)
fell within the typical ranges of other studies (Text S1). Although all
three Ag-MAR sites were classified as Traver fine sandy loam, the P1 andP3 sites exhibited a clear cemented duripan at 77–122 cm depth(Bachand et al., 2014), which resulted in lower saturated hydraulicconductivities (Ks ) at these depths (about 0.022 cm/min). In general, Ks
values were largest at P2, the sandiest profile of the three (84 % sand),and lowest at P1 (41 % sand), while P3 (61 % sand) was in between the
two. Optimized adsorption and degradation coefficients for the pesti­
cides (Tables S4-S6) were largest at P1 and lowest at P2, following thepattern of their total organic matter contents (Fig. 7). In addition,
chlorantraniliprole and methoxyfenozide had relatively higher adsorption
but lower degradation coefficients than imidacloprid and thiamethoxam,
indicating that they were less mobile and more persistent in theenvironment.
The observed surface ponding levels (Fig. 3) quickly increased totheir maximum (about 11 cm, 12 cm, and 21 cm in profiles P1, P2, andP3, respectively) because of continuous flooding. After that, the ponding
level decreased when water application stopped. The water contents at
all depths (Fig. 4) exhibited increasing trends during the flooding period
and decreased during the post-flooding period.
The bromide BTCs at all depths (Fig. 5) first showed increased con­centrations, followed by decreased concentrations with time. The P2
profile overall had the lowest concentrations, the fastest arrival rates ofpeak concentrations after the flooding, and the least pronounced tailingof BTCs. In contrast, the P1 profile had the highest concentrations, theslowest arrival times of peak concentrations after the flooding started,and the most pronounced tailing of BTCs.
Most pesticide BTCs (Fig. 6 and Figs. S1-S3) showed reduced con­centrations in the surface layers due to leaching from flooding, except
for methoxyfenozide, which increased in concentration at 20 cm due tothe pesticide transfer from the topmost layer. However, changes inpesticide concentrations within the deeper layers varied a lot, fromrising (mostly at P1), falling (mostly at P2), or initially rising before
falling (mostly at P2 and P3). These patterns indicate different arrival
rates of peak pesticide concentrations during flooding: slowest at P1,
fastest at P2, and medium at P3.
In the single-porosity models (SPM) fitted to the BTCs, simulated
wetting fronts progressively arrived later at increasing depths than thoseobserved, especially at P2 and P3. Similar to wetting fronts, simulated
bromide breakthroughs arrived later than those observed at depths at P2
and P3. The simulated bromide BTCs also had larger tailings than thoseobserved. The main mismatch in pesticide modeling between simula­tions and observations was the late arrival of BTCs. These may suggest
the presence of preferential flow/transport at P2 and P3 due to theoccurrence of mobile and immobile zones and the solute mass transfer
between them. Due to fast water movement during flooding, mixing
with immobile water is limited. As a result, solute mixing occurs mainly(22)

If a HYDRUS-1D run was not finished within a prescribed time (i.e.,
60 s for the single-porosity models and 600 s for the dual-porosity
5

T. Zhou et al.

Water Research 267 (2024) 122442

Fig. 3. Observed and simulated (using the single-porosity [SPM] and dual-porosity [DPM] models) surface ponding water levels in the three soil profiles P1, P2, andP3 (left to right). The blue-shaded areas indicate the flooding period.

Fig. 4. Observed and simulated (using the single-porosity [SPM] and dual-porosity [DPM] models) soil water contents at different depths (20, 60, 100, and 250 cm;
top to bottom) in the three soil profiles P1, P2, and P3 (left to right). The blue-shaded areas indicate the flooding period.

within the available mobile region (Imig et al., 2023b), leading to higher
observed solute concentrations and less tailing than simulated by SPM at
P2 and P3. This also highlights the potential utility of applying aweighted KGE (Gupta et al., 2009; Lamontagne et al., 2020), which al­lows for specifying the relative importance of the three KGE components(i.e., correlation, variability bias, and mean bias between observed andsimulated values) to better capture solute arrivals (Table 1).
The dual-porosity models (DPM) improved the model performancein simulating surface ponding levels, water contents, and bromide BTCs
at P1, but worsened it for pesticide BTCs. DPM improved the model
performance in simulating surface ponding levels and bromide, imida­
cloprid, and chlorantraniliprole BTCs at P2. However, it worsened themodel performance in simulating soil water contents and the thiame­thoxam and methoxyfenozide BTCs. DPM improved the model perfor­
mance at P3 in simulating surface ponding levels and bromide,
imidacloprid, and methoxyfenozide BTCs. However, it worsened it for soil
water contents and thiamethoxam and chlorantraniliprole BTCs. In otherwords, the model performance in simulating the bromide and pesticides
BTCs cannot be improved simultaneously, and it can even be worsened
when switching from SPM to DPM. This may be because the solutetransport parameters were optimized for bromide BTCs and thus may
not apply to pesticide BTCs.
Overall, SPM and DPM could both capture the general trends ofobserved surface ponding levels, soil water contents, and bromide/
pesticide BTCs (Table 1). Köhne et al. (2009) compared different

pesticide transport models and concluded that a model gives satisfactory
predictions if the ratio between measured and modeled concentrations is<3–5 times (Köhne et al., 2009). The ratio was generally <3 in thisstudy, indicating that the fit was fairly good. The model performance forimidacloprid, thiamethoxam, and chlorantraniliprole was much better thanfor methoxyfenozide. This may be related to the omission of potential
chemical nonequilibrium in this study since rapid water flow during
intensive flooding made it more difficult to reach equilibrium adsorption(Dusek et al., 2015).
Surface ponding level
Water content
BromideImidacloprid
3.2. Water mass balances and travel time of flooding water
The water balance and recharge amounts estimated with the SPMs
and DPMs showed distinct differences with variations of up to 6 %(Fig. 7, Table 2). Groundwater recharge was largest at P2 (SPM: 89.7 %,DPM: 90.8 %), smallest at P1 (82.6 %, 83.4 %), and intermediate (butclose to P1) at P3 (83.6 %, 83.6 %). The groundwater recharge varied by
up to 8 % between P2 and the other two profiles.
Considering preferential flow by DPMs resulted in lower bromide
travel times by up to 23 % and higher flow velocities by up to 31 %compared to SPMs. As expected, the sandiest profile P2 had the highest
6

T. Zhou et al.

Water Research 267 (2024) 122442

Fig. 5. Observed and simulated (using the single-porosity [SPM] and dual-porosity [DPM] models) bromide concentrations at different depths (20, 60, 100, and 250
cm; top to bottom) in the three soil profiles P1, P2, and P3 (left to right). The blue-shaded areas indicate the flooding period.

Fig. 6. Observed and simulated (using the single-porosity [SPM] and dual-porosity [DPM] models) imidacloprid concentrations at different depths (20, 60, 100, and250 cm; top to bottom) in the three soil profiles P1, P2, and P3 (left to right). The blue-shaded areas indicate the flooding period.

flow velocity, followed by P1, while P3 had the lowest. Accordingly,
bromide transport velocities between the three profiles differed by about
42 %, ranging between 16.54 and 91.84 cm/day in the 2.5 m nearsurface unsaturated zone (Fig. 7, Table 3).

thiamethoxam had significantly higher leaching rates, averaging 22.3 %more than chlorantraniliprole and methoxyfenozide (Fig. 7, Table 4).
With respect to degradation (Dp ), P1 had the highest degradation
rate, averaging 31.9 %, while P2 and P3 were at similar levels (27.7 %and 26.6 %, respectively). The DPM model predicted about 2.7 % moredegradation on average than the SPM, with P1 showing a 4.8 % higher
degradation than the other sites. Imidacloprid and thiamethoxam alsodegraded more, with an average of 19.1 % higher than chlorantraniliprole
and methoxyfenozide (Fig. 7, Table 4).
In terms of retention (Sp,final ), P1 and P3 retained the most (averaging
56.1 % and 57.9 %, respectively), while P2 retained the least (43.2 %).
Retention was about 3.8 % higher in the DPMs than in the SPMs, with P1

3.3. Pesticide mass balance
Among the sites, P2 showed the highest leaching efficiency (Lp ) withan average of 37.5 %, followed by P3 at 19.4 %, and P1 at 12.2 %. The
DPM model predicted slightly higher leaching (1.8 % on average)
compared to the SPM model, with P2 experiencing approximately 21.6
% more leaching than the other sites on average. Imidacloprid and7

T. Zhou et al.

Water Research 267 (2024) 122442

root zone but gains in the deep vadose zone (Fig. 7, Table 4).
Overall, P1 and P3 showed more retention or degradation, while P2
exhibited greater leaching. These findings are supported by the fact thatthe adsorption and degradation coefficients were largest at P1, followedby P3, while they were lowest at P2, as discussed in Section 3.1.
Furthermore, the DPMs generally showed more leaching, degradation,and retention than the SPMs

Table 1Model performance (Kling-Gupta Efficiency) of the single-porosity [SPM] anddual-porosity [DPM] models to simulate surface ponding levels, soil watercontents, and concentrations of bromide and four pesticides (imidacloprid, thia­
methoxam, chlorantraniliprole, and methoxyfenozide) in the three soil profiles (P1,
P2, and P3).
Variable
Surface ponding level
Water content
BromideImidacloprid
Thiamethoxam
Chlorantraniliprole
Methoxyfenozide

P1

P2

P3

SPM

DPM

SPM

DPM

SPM

DPM

0.67
0.87
0.44
0.78
0.67
0.81
0.65

0.76
0.89
0.38
0.73
0.59
0.56
0.58

0.02
0.94
− 0.69− 0.08
0.68
0.24
0.2

0.1
0.93
− 0.38
0.7
0.66
0.31
0.14

0.44
0.83
0.27
0.7
0.7
0.89
0.81

0.6
0.8
0.3
0.73
0.65
0.75
0.85

3.4. Model sensitivity analysis
The Sobol’ global sensitivity analysis was conducted to identify themost influential parameters controlling the transport and reaction ofpesticides within the SPMs (Table 5 and Text S2). The reaction param­
eters (especially the adsorption coefficient (Kd ) of Layer 1 (the surface
layer) were the most influential parameters in simulating chloran­
traniliprole at all three profiles and imidacloprid at P2. This finding un­
derscores the critical role that organic matter in the surface layer playsin the adsorption of pesticides, which, in turn, affects their movement
through the soil (Fig. 7). In addition to the surface layer’s role, thedegradation parameter (μ) of Layer 3 (the restrictive layer where thecemented duripan occurs) was the most impactful parameter in simu­
lating imidacloprid, thiamethoxam, and methoxyfenozide at P1 and P3.
This is likely because the presence of a restrictive layer at P1 and P3
increased the residence time of the infiltrating water and the potential

and P3 showing an average of 8.2 % more retention than P2. Chloran­
traniliprole and methoxyfenozide exhibited greater retention in the soil,
averaging 33.9 % more than imidacloprid and thiamethoxam (Fig. 7,
Table 4).
Regarding changes in pesticide storage (ΔSp,RZ and ΔSp,DVZ ), imida­
cloprid and thiamethoxam in general lost significant amounts of pesti­
cides from both the root zone and deep vadose zone, whilechlorantraniliprole and methoxyfenozide typically showed losses in theFig. 7. Conceptual model of flow and pesticide transport processes, water and pesticide mass balances, and bromide travel times during Ag-MAR in the top 2.5 m of
the three soil profiles (P1, P2, and P3). F: Flooding; E: Evaporation; H: horizontal flow caused by a capillary barrier; D: Deep drainage. The two numbers to the rightof these terms are water mass amounts (in cm) calculated using the single-porosity (SPM) and dual-porosity (DPM) models. The bottom bar plots show adsorption,
degradation, and leaching amounts (in ppb.cm) for the four pesticides, including imidacloprid (IMCP), thiamethoxam (TMTX), chlorantraniliprole (CRNP), andmethoxyfenozide (MTFZ), at the end of the simulations. TOC and MIM represent the total organic carbon content (%) and Mobile-Immobile zones, respectively.
8

T. Zhou et al.

Water Research 267 (2024) 122442

Table 2
Water mass balance components for the three soil profiles P1, P2, and P3 calculated using the SPMs and DPMs.
Term

P1

P2

SPM

DPM

cm
P+I
E
D
ΔSRZ
ΔSDVZ
GR

128.4
21.7
99.1
4.6
6.9
106.0

%

Difference

SPM

%

(%)

cm

​

​

cm
128.4
19.5
94.0
4.4
13.1
107.1

​

16.9
77.2
3.6
5.4
82.6

15.2
73.2
3.4
10.2
83.4

− 1.7
− 4.0
− 0.2
4.8
0.9

128.6
11.3
119.7
5.7
− 4.3
115.4

P3
DPM
%

cm

8.8
93.1
4.4
− 3.3
89.7

128.6
9.0
114.2
7.6
2.6
116.8

Difference

SPM

%

(%)

cm

​

​

7.0
88.8
5.9
2.0
90.8

− 1.8
− 4.3
1.5
5.3
1.1

DPM
%

128.4
14.3
98.3
10.6
9.0
107.3
​

11.1
76.6
8.3
7.0
83.6

Difference

cm
128.4
14.1
105.6
11.1
1.7
107.3

%

(%)

​

​

11.0
82.2
8.6
1.3
83.6

− 0.1
5.6
0.3
− 5.7
0

P: precipitation, F: flood irrigation, E: evaporation, D: drainage, ΔS: storage change in the root zone 0~100 cm (ΔSRZ ) and deep vadose zone 100–250 cm (ΔSDVZ ), GR:
groundwater recharge including D and ΔSDVZ (since water flow is considered one-dimensional, deep drainage below the root zone will eventually recharge
groundwater) (de Vries and Simmers, 2002).

Table 3
Travel times and average velocities of bromide (calculated by the peak displacement method) from the soil surface to different soil depths at three soil profiles P1, P2,
and P3 calculated using the SPMs and DPMs.
Term
Travel time (day)Flow velocity (cm/day)

Depth (cm)
20
60
100
175
250
20
60
100
175
250

P1

P2

P3

SPM

DPM

Relative difference (%)

SPM

DPM

Relative difference (%)

SPM

DPM

Relative difference (%)

1.13
1.79
2.67
3.29
5.63
17.71
33.54
37.46
53.2
44.41

1.08
1.63
2.29
3.63
4.5
18.53
36.83
43.68
48.22
55.57

− 4.4
− 8.9
− 14.2
10.3− 20.1
4.6
9.8
16.6
− 9.4
25.1

1.21
1.83
2.96
4.33
5.17
16.54
32.8
33.79
40.42
48.36

1.21
1.71
2.54
3.5
3.96
16.54
35.1
39.38
50.01
63.14

0
− 6.6
− 14.2− 19.2− 23.4
0
7
16.5
23.7
30.6

0.22
1.13
2.47
4.38
5.93
91.84
53.2
40.52
39.97
42.17

0.22
1.13
2.13
3.72
4.88
91.84
53.2
47
47.07
51.25

0
0
− 13.8− 15.1− 17.7
0
0
16
17.8
21.5

Note that the peak displacement method estimates travel times from the time lag between signals in the input (flooding) and output (different soil depths) of bromide
BTCs (Zhou et al., 2023).

for degradation. This finding again validates the conclusion that P1 andP3 show more retention or degradation, while P2 show more leaching,as discussed in Section 3.3.

4.2. Impacts of lateral flow on pesticide fate and transport
In addition to preferential flow, lateral flow also played a crucial rolein pesticide transport, especially in profiles P1 and P3. The low watercontent (about 0.04 cm3/cm3 at P1 and 0.05 cm3/cm3 at P3, Fig. 4) inthe deeper soil layers indicate unsaturated soil conditions. These un­saturated conditions below the duripan and notably higher pesticideconcentrations at 100 cm than in the other layers (Fig. 6 and Figs. S1-S3)
in the P1 and P3 profiles could potentially be related to the restrictive
layer forming a capillary barrier that promotes lateral flow within the
low-conductivity layer to the overlying coarser-textured soil layers(Fig. 7). This lateral flow led to the perching of water and pesticides at
specific depths (Zhou et al., 2023), particularly affecting less mobile
pesticides like chlorantraniliprole and methoxyfenozide which exhibited
higher degradation and retention than other two pesticides (Table 4).
This is also supported by the sensitivity analysis results which indicatethat the degradation coefficient of the restrictive layer is a key factor in
the behavior of less mobile pesticides like chlorantraniliprole andmethoxyfenozide in profiles P1 and P3 (Table 5).

4. Discussion
4.1. Impacts of preferential flow on pesticide fate and transport
Preferential flow paths were more pronounced in the sandier pro­files, such as P2 and P3, as indicated by the earlier arrival of observed
water and solute fronts, higher peak concentrations, and narrower
bromide and pesticide BTCs than simulated (e.g., for imidacloprid andchlorantraniliprole at P2 and thiamethoxam at P3) (Figs. 5-6 and Figs. S1S3). Preferential flow led to faster water flow (Table 3), more ground­water recharge (Table 2) and more pesticide leaching (Table 4) in theunsaturated zone at P2 and P3, especially for more mobile pesticides likeimidacloprid and thiamethoxam. The possible mechanisms that caused
preferential were elaborated in similar Ag-MAR research (Zhou et al.,2023).
Interestingly, considering preferential flow using the dual-porosity
models (DPMs) predicted higher degradation and retention rates eventhough they showed shorter travel times. This is because DPMs accountfor the unique interaction between mobile and immobile regions in thesoil. In DPMs, pesticides can move quickly through the mobile regions,but a significant portion gets temporarily trapped in the immobile re­gions, where the water is largely stagnant. This trapping allows pesti­
cides to stay in contact with the soil for longer periods, which enhances
both their retention and degradation (Ray et al., 2004).

4.3. Differences in pesticide behavior
Our results reveal distinct differences in the behavior of the four
pesticides studied—imidacloprid, thiamethoxam, chlorantraniliprole, andmethoxyfenozide—highlighting the importance of both chemical prop­
erties and soil conditions. We showed that Imidacloprid and thiame­thoxam were more mobile (more prone to leaching as seen in Table 4)
than chlorantraniliprole, which is consistent with previous research(Tingle et al., 2003). Because of this, imidacloprid and thiamethoxam
were significantly affected by preferential flow, particularly in thesandier P2 profile (discussed in Section 4.1). Conversely,9

T. Zhou et al.

Table 4
Solute mass balance components for the four pesticides at the three soil profiles P1, P2, and P3 calculated using the SPMs and DPMs.
Pesticide

Term

P1

P2

SPM
ppb⋅cm10

Imidacloprid
​
​
​
​
​
Thiamethoxam
​
​
​
​
​
Chlorantraniliprole
​
​
​
​
​
Methoxyfenozide
​
​
​
​
​

Sp,initLp
Dp
Sp,final
ΔSp,RZ
ΔSp,DVZ
Sp,initLp
Dp
Sp,final
ΔSp,RZ
ΔSp,DVZ
Sp,initLp
Dp
Sp,final
ΔSp,RZ
ΔSp,DVZ
Sp,initLp
Dp
Sp,final
ΔSp,RZ
ΔSp,DVZ

280.2
9.5
180.9
89.0
− 120.5− 70.71052.8
576.0
69.0
409.6
− 232.0− 411.1
4394.08.8
1073.43304.1− 1035.3− 54.7
440.0
7.4
200.4
232.1
− 103.0− 104.9

DPM
%
​

3.4
64.6
31.8
​
​
​

54.7
6.6
38.9
​
​
​

0.2
24.4
75.2
​
​

​

​

​

1.7
45.5
52.8

ppb⋅cm
363.6
10.8
161.6
189.7
− 118.6− 55.31527.8471.0
17.3
1032.7− 309.7− 185.5
2599.623.1
730.7
1844.7− 681.1− 73.8
336.4
9.8
137.1
198.5
− 69.3− 68.5%Difference
%

​

​

​

​
​

​
​

​

​

​
​

​
​

​

​

​
​

​
​

​

​

​
​

​

3.0
44.4
52.5

30.8
1.1
67.6

− 0.4
− 20.2
20.7

0.9
28.1
71.0

− 23.9
− 5.5
28.7

2.9
40.8
59.0

0.7
3.7
− 4.2

1.2
− 4.7
6.2

P3

SPM
ppb⋅cm
31.6
18.8
0
12.8
− 18.2
− 0.6
592.8
164.0
382.7
45.0
− 136.1− 411.7
1689.4
362.4
448.2
879.2
− 861.2
51.0
170.4
29.1
0
141.3
− 32.6
3.5

DPM
%
​

59.5
0.0
40.5
​
​
​

27.7
64.6
7.6
​
​
​

21.5
26.5
52.0
​
​

​

​

​

17.1
0.0
82.9

ppb⋅cm
33.1
23.2
17.0
14.5
− 16.8
− 1.9
421.7
132.0
244.2
35.6
− 92.6− 293.5
1054.9
581.0
206.8
289.9
− 244.0
479.0
120.7
19.3
1.7
99.9
− 24.6
3.8%

Difference
%

​

​

​

​
​

​
​

​

​

​
​

​
​

​

​

​
​

​
​

​

​

​
​

​

71.4
51.4
43.8

31.3
57.9
8.4

15.1
54.7
6.1

55.1
19.6
27.5

3.6
− 6.7
0.8

16.0
1.4
82.8

33.6
− 6.9
70.3

− 1.1
1.4
− 0.1

SPM
ppb⋅cm
147.2
42.9
51.0
53.4
− 102.1
8.3
874.9
328.5
340.7
206.3
− 303.4− 365.1
702.9
65.7
39.1
598.4
− 181.0
76.4
326.4
13.2
55.1
258.2
− 75.8
7.6

DPM
%
​

29.1
34.6
36.3
​
​
​

37.5
38.9
23.6
​
​
​

9.3
5.6
85.1
​
​

​

​

​

4.0
16.9
79.1

ppb⋅cm117.68
55.8
43.1
50.3
79.0
11.6
1643.9
234.6
963.4
517.4
182.7
− 666.5
561.6
58.1
39.8
460.6
− 162.6
61.7
344.8
12.0
48.7
285.4
− 65.2
5.7%

Difference
%

​

​

​

​
​

​
​

​

​

​
​

​
​

​

​

​
​

​
​

​

​

​
​

​

47.4
36.6
42.7

14.3
58.6
31.5

18.3
2.0
6.4

10.4
7.1
82.0

− 23.3
19.7
7.9

3.5
14.1
82.8

1.1
1.5
− 3.1

− 0.5
− 2.8
3.7

Note that Sinit and Sfinal are the initial and final pesticide storages in the soil profile, respectively, Lp is the pesticide mass leached through drainage, Dp is the pesticide mass degraded through chemical or biologicalreactions, and ΔSp is pesticide storage change in the root zone 0~100 cm (ΔSp,RZ ) and deep vadose zone 100–250 cm (ΔSp,DVZ ).
Water Research 267 (2024) 122442

T. Zhou et al.

Water Research 267 (2024) 122442

Table 5The most influential soil layers and reaction parameters (adsorption coefficient Kd ; degradation coefficient μ) for modeling (using the single-porosity [SPM]) of four
pesticides (imidacloprid, thiamethoxam, chlorantraniliprole, methoxyfenozide) at the three soil profiles (P1, P2, and P3).

chlorantraniliprole and methoxyfenozide, which are less mobile, weremore influenced by soil heterogeneity and lateral flow, leading togreater retention and degradation in profiles with finer textures andmore organic matter, such as P1 and P3 (discussed in Section 4.2). The
DPMs’ ability to capture these nuanced behaviors underscores its value
in simulating the real-world dynamics of pesticide fate and transport in
Ag-MAR systems.

research should aim to expand the scope of pesticides studied andexplore the use of three-dimensional models to capture more detailed
spatial variations in Ag-MAR systems.
5. Conclusions
This study emphasized the significant role of soil texture and het­
erogeneity in influencing the fate and transports of pesticides in thevadose zone under agricultural managed aquifer recharge (Ag-MAR)
through detailed field experiments and numerical modeling using both
single and dual-porosity models in HYDRUS-1D. Based on these waterflow and pesticide transport information, the suitability of implement­
ing Ag-MAR was evaluated by combining other biogeochemical in­
dicators including oxygen, oxidation–reduction potential, and electricalconductivity levels in soil water.
The sandier profile (P2) demonstrated more pronounced preferential
flow, faster water flow, higher groundwater recharge efficiency, but also
more soil pesticide and salt leaching correlating with increased risks ofgroundwater contamination. Therefore, choosing less mobile pesticides(e.g., chlorantraniliprole and methoxyfenozide) and extending the timeinterval between the last pesticide application and the Ag-MAR event
are important for reducing groundwater pollution risks. In contrast,
profiles with finer textures (P1) showed less preferential flow, slower
water flow, lower groundwater recharge efficiency, but raised concerns
about the accumulation of pesticides or salts and elongated anaerobic
conditions in upper soil layers due to capillary barriers, which could
severely impact soil structure, biochemical reactions and crop growth.
Our results underscore the necessity of considering soil heterogene­
ity and implementing site-specific management practices in the designand operation of Ag-MAR systems to maximize groundwater recharge
while minimizing environmental risks. Furthermore, our study calls forthe integration of more advanced models that can adequately capture
the complex processes of pesticide fate and transports, such as chemicaland physical nonequilibrium processes. On the other hand, the relative
magnitude of preferential flow/transport and lateral flow remains to be
distinguished.

4.4. Implications for Ag-MAR practice and future research
The three profiles had the same land use and hydroclimatologicalconditions, yet they varied in the subsurface hydrogeology. Conse­
quently, this study offers valuable insights into the impacts of soil
texture and heterogeneity on groundwater recharge and pesticideleaching that may be encountered when implementing Ag-MAR.
In terms of groundwater recharge and pesticide residue leaching, P2(the soil profile with the highest sand content) tended to facilitate
preferential flow, leading to faster water flow and pesticide transport (byabout 42 %) and increased recharge rates (by approximately 8 %). These
conditions also resulted in greater pesticide leaching due to flooding (byabout 22 %). Therefore, Ag-MAR should be implemented with caution insandy soils. To reduce this risk, it is recommended to give the soil time toabsorb and break down these chemicals before any planned rechareactivities. Choosing the right type of pesticides also matters. Imidaclo­prid and thiamethoxam, are more likely to leach into groundwater,
especially in sandy soils. In contrast, pesticides like chlorantraniliprole
and methoxyfenozide, which are less mobile and break down more easily,pose a lower risk.
For oxygen (O₂) levels (Figs. S7-S9) during Ag-MAR flooding, there
were noticeable sharp drops especially for shallow layers at P1 and P3,
which validated the water accumulation at shallow layers because of therestrictive layers mentioned in Section 4.2. These significant reductionsof O₂ levels indicated the development of anaerobic conditions, which
could potentially restrict crop growth. Correspondingly, P1 and P3 also
had relatively lower oxidation–reduction potential (ORP) levels than P2
during Ag-MAR (Figs. S7-S9).
There were significant electrical conductivity (EC) spikes (evenhigher than irrigation water) at P3 during Ag-MAR, which highlighted afaster soil salt leaching and higher groundwater salinization risk. These
ions could originate from the mobilization and release of fertilizer or
pesticide residue in the topsoil. P2 and P1 also presented EC spikes, but
in general lower than irrigation water and to a lesser extent than P3,
suggesting some issues with salt leaching and groundwater salinization
but potentially more manageable.
Overall, each profile illustrated different challenges and opportu­
nities for groundwater recharge, environmental protection, and agri­
cultural productivity. Tailoring management practices to the specificconditions reflected in these profiles can balance these aspects.
While this study provides valuable insights, it is important torecognize its limitations. The reliance on HYDRUS-1D models, thoughrobust, does not capture the full three-dimensional complexity of waterand solute movement in the field. Additionally, the study focused on alimited set of pesticides, which may not represent the full range of be­
haviors exhibited by other chemicals under similar conditions. Future

CRediT authorship contribution statement
Tiantian Zhou: Writing – original draft, Visualization, Methodol­
ogy, Investigation, Formal analysis. Nels Ruud: Validation, Resources,
Methodology, Conceptualization. Jiří Šimůnek: Writing – review &editing, Supervision, Software, Methodology. Giuseppe Brunetti:
Writing – review & editing, Validation, Supervision, Methodology. Elad
Levintal: Supervision, Methodology, Investigation, Data curation,
Conceptualization. Cristina Prieto García: Methodology, Investigation,
Formal analysis, Data curation. Helen E. Dahlke: Writing – review &editing, Supervision, Resources, Project administration, Funding
acquisition, Conceptualization.
Declaration of competing interest
The authors declare that they have no known competing financial
interests/personal relationships that could have appeared to influence11

T. Zhou et al.

Water Research 267 (2024) 122442

the work reported in this paper.

Knoben, W.J.M., Freer, J.E., Woods, R.A., 2019. Technical note: inherent benchmark or
not? Comparing nash-sutcliffe and kling-gupta efficiency scores. Hydrol. Earth Syst.
Sci. 23 (10), 4323–4331. https://doi.org/10.5194/hess-23-4323-2019.
Köhne, J.M., Köhne, S., Šimůnek, J., 2006. Multi-process herbicide transport instructured soil columns:: experiments and model analysis. J. Contam. Hydrol. 85(1–2), 1–32. https://doi.org/10.1016/j.jconhyd.2006.01.001.
Köhne, J.M., Köhne, S., Šimůnek, J., 2009. A review of model applications for structured
soils: b) Pesticide transport. J. Contam. Hydrol. 104 (1–4), 36–60. https://doi.org/
10.1016/j.jconhyd.2008.10.003.
Lamontagne, J.R., Barber, C.A., Vogel, R.M., 2020. Improved estimators of model
performance efficiency for skewed hydrologic data. Water Resour. Res. 56 (9), 25.
https://doi.org/10.1029/2020wr027101.
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.
Lucier, G., Jerardo, A., 2006. Vegetables and Melons outlook. VGS-317. US Departmentof Agriculture, Economic Research Service.
Mottes, C., Lesueur-Jannoyer, M., Le Bail, M., Malézieux, E., 2014. Pesticide transfer
models in crop and watershed systems: a review. Agron. Sustain. Dev. 34 (1),229–250. https://doi.org/10.1007/s13593-013-0176-3.
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.
Pang, L.P., Close, M.E., Watt, J.P.C., Vincent, K.W., 2000. Simulation of picloram,atrazine, and simazine leaching through two New Zealand soils and into
groundwater using HYDRUS-2D. J. Contam. Hydrol. 44 (1), 19–46. https://doi.org/
10.1016/s0169-7722(00)00091-7.
PCPA, 1985.
Perez, N., Singh, V., Ringler, C., Xie, H., Zhu, T.J., Sutanudjaja, E.H., Villholth, K.G.,
2024. Ending groundwater overdraft without affecting food security. Nat. Sustain.
14. https://doi.org/10.1038/s41893-024-01376-w.
Radcliffe, D.E., Šimůnek, J., 2018. Soil Physics With HYDRUS: Modeling andApplications. CRC Press.
Ray, C., Vogel, T., Dusek, J., 2004. Modeling depth-variant and domain-specific sorption
and biodegradation in dual-permeability media. J. Contam. Hydrol. 70 (1–2), 63–87.
https://doi.org/10.1016/j.jconhyd.2003.08.009.
Scorza, R.P., Boesten, J., 2005. Simulation of pesticide leaching in a cracking clay soil
with the PEARL model. Pest Manag. Sci. 61 (5), 432–448. https://doi.org/10.1002/
ps.1004.
Shi, Y.H., Eberhart, R., 1998. A modified particle swarm optimizer. In: IEEE International
Conference on Evolutionary Computation. Anchorage, Ak. Ieee, pp. 69–73. https://doi.org/10.1109/icec.1998.699146.
Sidoli, P., Lassabatere, L., Angulo-Jaramillo, R., Baran, N., 2016. Experimental andmodeling of the unsaturated transports of S-metolachlor and its metabolites inglaciofluvial vadose zone solids. J. Contam. Hydrol. 190, 1–14. https://doi.org/
10.1016/j.jconhyd.2016.04.001.
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.
Stagnitti, F., Parlange, J.-Y., Steenhuis, T.S., Boll, J., Pivetz, B., Barry, D., 1995.
Transport of moisture and solutes in the unsaturated zone by preferential flow.
Environmental Hydrology. Springer, pp. 193–224.
Tingle, C.C.D., Rother, J.A., Dewhurst, C.F., Lauer, S., King, W.J., 2003. Fipronil:
environmental fate, ecotoxicology, and human health concerns. In: Ware, G.W. (Ed.),
Reviews of Environmental Contamination and Toxicology: Continuation of Residue
Reviews. Springer New York, New York, NY, pp. 1–66. https://doi.org/10.1007/
978-1-4899-7283-5_1.
Tudi, M., Ruan, H.D., Wang, L., Lyu, J., Sadler, R., Connell, D., Chu, C., Phung, D.T.,
2021. Agriculture development, pesticide application and its impact on theenvironment. Int. J. Environ. Res. Public Health 18 (3), 23. https://doi.org/
10.3390/ijerph18031112.
Waskom, R., (1995), Best management practices for agricultural pesticide use.
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.
Zhou, T., Levintal, E., Brunetti, G., Jordan, S., Harter, T., Kisekka, I., Šimůnek, J.,
Dahlke, H.E., 2023. Estimating the impact of vadose zone heterogeneity on
agricultural managed aquifer recharge: a combined experimental and modelingstudy. Water Res. 247, 120781. https://doi.org/10.1016/j.watres.2023.120781.
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 612, 128100. https://doi.org/10.1016/j.
jhydrol.2022.128100.

Data availability
Data will be made available on request.

Acknowledgments
We acknowledge Mr. Don Cameron from Terranova Ranch for his
support in collecting original data. This work was partially funded by
USDA NIFA (Award # 2021–68012–35914), USDA NRCS (grant noNR203A750023C017), the Gordon and Betty Moore Foundation (grantno 7975), and the California Department of Pesticide Regulation(Agreement 17-C0099).
Supplementary materials
Supplementary material associated with this article can be found, inthe online version, at doi:10.1016/j.watres.2024.122442.
References
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.
Brunetti, G., Šimůnek, J., Piro, P., 2016. A comprehensive numerical analysis of thehydraulic behavior of a permeable pavement. J. Hydrol 540, 1146–1161. https://doi.org/10.1016/j.jhydrol.2016.07.030.
Brunetti, G., Šimůnek, J., Turco, M., Piro, P., 2017. On the use of surrogate-basedmodeling for the numerical analysis of low impact development techniques. J Hydrol
548, 263–277. https://doi.org/10.1016/j.jhydrol.2017.03.013.
Brunetti, G., Stumpp, C., Šimůnek, J., 2022. Balancing exploitation and exploration: anovel hybrid global-local optimization strategy for hydrological model calibration.
Environ. Modell. Softw. 150, 11. https://doi.org/10.1016/j.envsoft.2022.105341.
Czapar, G.F., Horton, R., Fawcett, R.S., 1992. Herbicide and tracer movement in soilcolumns containing an artificial macropore. J. Environ. Qual. 21 (1), 110–115.
https://doi.org/10.2134/jeq1992.00472425002100010016x.
Dahlke, H.E., LaHue, G.T., Mautner, M.R.L., Murphy, N.P., Patterson, N.K.,
Waterhouse, H., Yang, F.F., Foglia, L., 2018. Managed aquifer recharge as a tool toenhance sustainable groundwater management in california: examples from field
and modeling studies. In: Friesen, J., RodriguezSinobas, L. (Eds.), Advanced Tools
for Integrated Water Resources Management. Academic Press Ltd-Elsevier Science
Ltd, London, pp. 215–275. https://doi.org/10.1016/bs.apmp.2018.07.003.
Advances in Chemical Pollution Environmental Management and Protection.
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.
Dusek, J., Dohnal, M., Snehota, M., Sobotkova, M., Ray, C., Vogel, T., 2015. Transport ofbromide and pesticides through an undisturbed soil column: a modeling study with
global optimization analysis. J. Contam. Hydrol. 175, 1–16. https://doi.org/
10.1016/j.jconhyd.2015.02.002.
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.
Gupta, H.V., Kling, H., Yilmaz, K.K., Martinez, G.F., 2009. Decomposition of the meansquared error and NSE performance criteria: implications for improving hydrologicalmodelling. J. Hydrol 377 (1–2), 80–91. https://doi.org/10.1016/j.
jhydrol.2009.08.003.
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 313 (3–4), 257–273. https://doi.org/
10.1016/j.jhydrol.2005.03.035.
Holbak, M., Abrahamsen, P., Diamantopoulos, E., 2022. Modeling preferential waterflow and pesticide leaching to drainpipes: the effect of drain-connecting and matrixterminating biopores. Water Resour. Res. 58 (7), 21. https://doi.org/10.1029/
2021wr031608.
Imig, A., Augustin, L., Groh, J., Putz, T., Elsner, M., Einsiedl, F., Rein, A., 2023a. Fate ofherbicides in cropped lysimeters: 2. leaching of four maize herbicides considering
different processes. Vadose. Zone J. 22 (5), 14. https://doi.org/10.1002/vzj2.20275.
Imig, A., Augustin, L., Groh, J., Pütz, T., Zhou, T., Einsiedl, F., Rein, A., 2023b. Fate ofherbicides in cropped lysimeters: 1. Influence of different processes and model
structure on vadose zone flow. Vadose Zone J. e20265.

12