Received: 1 October 2019

Accepted: 29 March 2020

DOI: 10.1002/vzj2.20027
Vadose Zone Journal

ORIGINAL RESEARCH ARTICLES

Scaling factors in HYDRUS to simulate a reduction in hydraulicconductivity during infiltration from recharge wells andinfiltration basins
Jana Glass1Jiří Šimůnek2

1 Research Group INOWAS, Dep. of

Hydrosciences, Technische Univ. Dresden,
Dresden, Germany
2 Dep. of Environmental Sciences, Univ. of

California Riverside, Riverside, CA, 92521,
USA
Correspondence
Jana Glass, Research Group INOWAS, Dep.
of Hydrosciences, Technische Univ. Dresden,
Dresden, Germany.
Email: jana.glass@tu-dresden.de
Funding information
Bundesministerium für Bildung und
Forschung, Grant/Award Number:
01LN1311A

Catalin Stefan1

Abstract
Managed aquifer recharge (MAR) represents a promising technique to cope withincreasing water stress worldwide. However, it can be challenging to operate MAR
systems, especially concerning clogging. Clogging reduces the infiltration capacity
and system efficiency of a MAR facility. Processes that cause clogging are difficult toquantify and assess, and their simulations in MAR schemes have so far been limited.
The variably saturated water flow model HYDRUS-2D was therefore modified toinclude time-variable hydraulic conductivities to more realistically represent clogging
at the infiltration interface of infiltration basins and recharge wells. An exponentialfunction with a time-variable scaling factor was implemented into HYDRUS to vary
the soil hydraulic conductivity over time during simulations. The new approach
was tested, in combination with the reservoir boundary condition, by simulating
two-dimensional cross-sections of two three-dimensional laboratory experiments
representing recharge from infiltration basins and injection wells. With the help of
the time-variable scaling factor, the increasing ponding depth in both experiments
due to progressive clogging was reproduced. Hypothetical simulations with various
well configurations, soils, and injection rates indicate that clogging influences theinfiltration volumes and subsurface infiltration area, which should be considered during the planning of MAR systems. The approach can be used to evaluate the resulting
infiltration capacity numerically and to help with the operation and design of MAR
facilities. However, further research is required to fully understand and integrate theprocesses that lead to clogging at MAR facilities into numerical simulations.

1

Abbreviations: FE, finite element; HDPE, high-density polyethylene; HIE,
horizontal infiltration extent; MAR, managed aquifer recharge; 2D,
two-dimensional.

I N T RO D U C T I O NManaged aquifer recharge (MAR) represents a viable technique to cope with worldwide increasing water demand. It
is mainly applied to increase the local groundwater storage,
thereby overcoming the temporal imbalance between local
water demand and restricted water availability (Bouwer, 2002;

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original
work is properly cited.
© 2020 The Authors. Vadose Zone Journal published by Wiley Periodicals, Inc. on behalf of Soil Science Society of America
Vadose Zone J. 2020;19:e20027.
https://doi.org/10.1002/vzj2.20027

wileyonlinelibrary.com/journal/vzj2

1 of 19

2 of 19

Vadose Zone Journal

Dillon, 2005). Managed aquifer recharge schemes are built tofulfill the following objectives: to secure irrigation or drinking
water supplies throughout the year, to mitigate seawater intrusion in coastal aquifers, to improve the water quality throughsoil aquifer treatment, and to prevent land subsidence (Dillon,
2005; Stefan & Ansems, 2018), among other objectives. A
wide range of MAR methods exists, including well, shaft, andborehole recharge, spreading practices, induced bank filtration, in-channel modifications, and runoff harvesting (Hannappel, Scheibler, Huber, & Sprenger, 2014). Dry (recharge)
wells are frequently used to recharge stormwater throughthe unsaturated zone (Edwards, Harter, Fogg, Washburn, &Hamad, 2016). The suitability of a MAR method depends ona variety of factors, including the site-specific topography,
land use, hydrology, hydrogeology, the quality of the source
and ambient groundwater, and the intended use of recovered
water (Dillon, Pavelic, Page, Beringen, & Ward, 2009).
Managed aquifer recharge plays a vital role in the development of water supply systems in Europe, Australia, and theUnited States and contributes substantially to drinking waterproduction (Sprenger et al., 2017; Stefan & Ansems, 2018).
Despite its financial and ecological benefits, the contribution of MAR to safe water supply at the global scale isstill somewhat limited, as indicated by a MAR growth rateof 0.5 km3 yr−1 between 2005 and 2015, which did not keep
pace with the worldwide increase of groundwater use of
53 km3 yr−1 between 2000 and 2010 (Dillon et al., 2018). On
the global scale, MAR currently covers ∼1% of groundwater
extraction, which equals an estimated 10 km3 yr−1 (Dillon
et al., 2018). The lack of data on MAR technological costs,
hydrogeological site-specific characterization, and operational uncertainties, in connection with the absence of propermonitoring and evaluation methods, raises questions about
the impact of MAR on the affected ecosystems, reduces thelevel of public trust, and hinders the optimal operational
management (Dillon et al., 2018). Hence, new or enhanced
methods are required, which can help to advance the safeand reliable planning, management, and operation of MAR(Dillon et al., 2018).
Clogging mainly occurs at the infiltration interface of MAR
schemes and involves various chemical, physical, mechanical,
and biological processes. Whereas physical clogging mainly
encompasses the accumulation of suspended solids, biological clogging involves bacteria and algae growth and accumulation, as well as the microbiological production of polysaccharides (Martin, 2013; Pérez-Paricio, 2001). Mechanicalclogging may be caused by air entrainment, whereas chemical
clogging involves various types of geochemical reactions
such as mineral precipitation and dissolution, as well as ion
exchange and adsorption (Martin, 2013). Besides the achievable recovery efficiency, clogging represents the most criticalissue during the MAR operation, since it can greatly reduce
system efficiency, as it leads to reduced infiltration rates over

GLASS ET AL.

Core Ideas• We implemented time-variable scaling factors into
HYDRUS.
• Scaling factors reduce the hydraulic conductivity
over time to represent clogging.
• Hydraulic conductivity reduction in two laboratory experiments was well replicated.
• The proposed approach improves simulations of
two MAR schemes exposed to clogging.

time. Because each MAR facility experiences some clogging
during its operational lifetime, the evaluation of the clogging potential and appropriate remediation and mitigation
measures through engineering design and operational management are required (Martin, 2013). A variety of processes
including dissolution and precipitation of minerals, such as
calcite, Fe oxides, and hydroxides (Anderson, Dewhurst,
Jones, & Baxter, 2006; Gutiérrez-Ojeda, Martínez-Morales,& Ortiz-Flores, 2007; Page et al., 2014; Vanderzalm, Page,
Barry, & Dillon, 2013), the accumulation of suspended
solids (Pavelic et al., 2007; Rinck-Pfeiffer, 2000), as well as
algae and microbial growth (Page et al., 2014; Pavelic et al.,2007; Thompson, Stotler, Macpherson, & Liu, 2015), wereidentified to cause clogging at MAR schemes. Althoughlaboratory experiments (Rinck-Pfeiffer, 2000; Wood et al.,2006), and to some extent field experiments, have been used
to analyze clogging mechanisms (Pavelic et al., 2007), only
a few simulation studies have attempted to represent theseprocesses (Majumdar, Sekhar, Sridharan, & Mishra, 2008;
Masciopinto, 2013; Pérez-Paricio, 2001; Rinck-Pfeiffer et al.,2013; Ringleb, Sallwey, & Stefan, 2016).
The use of unsaturated water flow models is required forrealistic simulations of clogging processes at the infiltrationinterface of MAR facilities that infiltrate water into the unsaturated zone from infiltration basins and dry wells. However,the use of unsaturated zone models to assess MAR schemes,compared with the widespread application of saturated flow
models, has been so far somewhat limited (Sasidharan, Bradford, Šimůnek, DeJong, & Kraemer, 2018; Sasidharan, Bradford, Šimůnek, & Kraemer, 2019, 2020), whereas clogging
has been mostly neglected (Sallwey, Glass, & Stefan, 2018),
which leads to the overestimation of infiltration and recharge
rates. For example, Sasidharan et al. (2018, 2019, 2020) haverecently used HYDRUS-2D (Radcliffe and Šimůnek, 2010)
in combination with the reservoir boundary to assess experimental data collected at two dry wells, as well as the effectsof heterogeneity on drywell infiltration and groundwater
recharge, without implicitly accounting for clogging.

GLASS ET AL.

Vadose Zone Journal

The main objective of this manuscript is to describe theimplementation and validation of a simplified representationof clogging into the variably saturated flow model HYDRUS
to improve its ability to simulate processes occurring at the
infiltration interface of MAR facilities, especially when an
unsaturated zone component is present (Glaß, 2019). For this,a time-variable scaling factor in the form of an exponentialfunction was implemented into HYDRUS to reduce thehydraulic conductivity with time during simulations. In
combination with the reservoir boundary condition, which
was recently developed to account for ponded infiltration
at MAR facilities (Šimůnek, Šejna, & van Genuchten,2018), the approach was tested using two three-dimensional
laboratory experiments representing MAR via recharge wells
and infiltration basins. This manuscript focuses on numericalsimulations of two laboratory experiments, integrating atime-variable scaling factor to validate the simple representation of clogging at MAR facilities. In addition, hypotheticalsimulations with varying soils, well configurations, andinfiltration rates were conducted to quantify the effects of
the reduction of hydraulic conductivity on the infiltration
volumes and the extent of the saturated infiltration areaaround vadose zone wells.
The focus of this manuscript is on the infiltration process
in MAR systems rather than on groundwater recharge. The
complex temporal relationship between these two processes
depends mainly on the thickness of the vadose zone between
the infiltration system (e.g., recharge wells, irrigation basins)
and the ground water table, the heterogeneity of the vadosezone, the volume of infiltrated water, and other factors. This
relationship between irrigation and groundwater recharge isa subject of the study by Sasidharan, Bradford, Šimůnek, andKraemer (2020) and is beyond the scope of this manuscript.

2

M AT E R I A L S A N D M E T H O D S

2.1 Simulation of flow in the unsaturated
zone
Variable-saturated water flow was simulated using the
HYDRUS-2D model (Šimůnek, van Genuchten, & Šejna,2016). The program uses the Garlekin-type finite element(FE) scheme to numerically solve the Richards equation forvariable-saturated water flow:
∂
∂θ
=
∂𝑡
∂𝑥𝑖

[

(
)]
A ∂ℎ
A
+ 𝐾𝑖𝑥
𝐾 𝐾𝑖𝑗
−𝑆
2
∂𝑥𝑖

(1)

where θ is the volumetric water content [–], h is the pressure
head [L], S is a sink term [T−1 ], xi (i = 1,2) are spatial coordinates [L] where x1 is the horizontal coordinate and x2 is thevertical coordinate, Kij A are components of a dimensionless

3 of 19

hydraulic conductivity anisotropy tensor KA , K is the unsaturated hydraulic conductivity function given by the productof the relative hydraulic conductivity Kr and the saturated
hydraulic conductivity Ks [L T−1 ], and t is time [T].
Although the spatial discretization is carried out using the
FE method to account for irregular geometries and required
local FE refinements, the time is discretized into finite time
steps by applying a fully implicit finite difference scheme.
The relationship between pressure heads and water contents
is described using the soil water retention curve, which canbe parameterized according to van Genuchten (1980):{
θ (ℎ) =θ −θsr
θr + 1+(−α|ℎ|)𝑛 𝑚, ψ < 0
[
]
θs ,
ψ>0

(2)

where θs is the saturated water content [–], θr is the residual water content [–], and α [L−1 ], n [–], and m [–] areempirical fitting parameters (Šimůnek, Šejna, Saito, Sakai, &van Genuchten, 2013). The van Genuchten–Mualem equation(van Genuchten, 1980) is used to describe the relationshipbetween unsaturated hydraulic conductivity and soil saturation:
[
(
)𝑚 ]2
( )
1∕𝑚
𝐾 𝑆e = 𝐾s 𝑆e𝑙 1 − 1 − 𝑆e(3)where Se is the effective saturation [–], and l is the pore connectivity parameter, often assumed to be 0.5 for most soils(Mualem, 1976). With the help of Equation 3, the unsaturated
hydraulic conductivity function can be described using soil
water retention parameters.
The reservoir boundary condition is applied to simulate infiltration during laboratory infiltration experiments.
This boundary condition calculates the water level in theinfiltration unit, depending on the infiltration area andinfiltration rates (Šimůnek et al., 2018). The utilization of thereservoir boundary condition restricts the simulation to a twodimensional (2D) cross-section of the infiltration structure
or its 2D radially symmetric representation. A seepage faceboundary condition cannot be defined elsewhere in the transport domain while the reservoir boundary condition is used.
Temporal changes in the hydraulic conductivity in
HYDRUS-2D are considered using a time-dependent scaling
factor that decreases exponentially with time:
𝐾s,𝑡 = 𝐾s,initial 𝑒−λ𝑡(4)where Ks,t is the saturated hydraulic conductivity at time t [LT−1 ], Ks,initial is the initial saturated hydraulic conductivity [LT−1 ], and λ is the hydraulic conductivity reduction parameter(λ ≥ 0), also referred to as the clogging factor [T−1 ]. Similar
equations have been applied for almost 60 yr to describe theimpact of clogging on surface systems (Pérez-Paricio, 2001).

4 of 19

Vadose Zone Journal

GLASS ET AL.

F I G U R E 1 An experimental laboratory tank setup representing an infiltration basin. The figure also shows the locations of various measuring
devices, such as tensiometers (2–6), time domain reflectometry (TDR) probes (2–6), and a pressure sensor recording the water level in the infiltration
basin (1)F I G U R E 2 The particle size distribution of the experimental sand

2.2
2.2.1

Infiltration basin experiment
Experimental setup

The laboratory experiment to simulate MAR via an infiltration basin consisted of a three-dimensional, rectangularly
shaped, stainless steel tank (1.5 × 1.0 × 1.0 m, L × W
× H). An infiltration basin (0.44 × 0.27 × 0.06 m) was
installed in the center of its surface (Figure 1). The bottom of the tank consisted of a stainless-steel grid, which
was overlaid by a 14.5-cm filter layer with increasing

grain sizes (from sand to gravel) to avoid leaching of finematerial. Above the sand filter, the main experimental
soil was installed, which consists of 84 cm of sand (thebulk density of 1.6 g cm−3 ). The particle size distribution
is displayed in Figure 2. The tank was equipped withtensiometers and TDR sensors to measure time-variable
pressure heads and water contents during infiltration cycles,
respectively, as well as with suction cups and electricalconductivity sensors to monitor the water quality (Figure 1).
The sensors were installed at various positions at three different depths (0.15, 0.25, and 0.66 m). The water level in the

GLASS ET AL.

Vadose Zone Journal

5 of 19F I G U R E 3 The numerical model setup representing the infiltration basin, and showing the finite element grid, boundary conditions, andobservation points

TABLE 1

Soil hydraulic parameters used in simulations of waterflow in the laboratory infiltration basin tankSoilmaterial

Description

𝛉r

𝛉s

𝛂
m−1

n

Ks
cm min−11

Sand (exposed toclogging)

0.066 0.292 0.236 1.91 0.21

2

Sand (replaced at the 0.066 0.292 0.236 1.91 3.15
beginning of eachscenario)3

Sand (not replaced)

0.121 0.228 0.224 1.66 1.47

4

Sand filter

0.045 0.300 0.300 1.51 3.26

5

Coarse cobble filter

0.000 0.300 0.050 1.35 50.00

Note. θr , residual water content; θs , saturated water content; α and n, empiricalfitting parameters; Ks , saturated hydraulic conductivity.

infiltration basin was measured using a pressure sensor. The
tank was installed in a fully automatic climate tank where thetemperature was controlled at 17 ◦ C and the relative humidity
was controlled at 70% (Fichtner, Barquero, Sallwey, & Stefan,2019). A variety of scenarios with different wet/dry ratios(a water application followed by water redistribution) andinfiltration rates have been conducted (Fichtner et al., 2019).
To demonstrate the applicability of the time-variable scaling factor to reproduce clogging in laboratory measurements,
the numerical analysis focused on one experimental scenario(unpublished raw data provided by T. Fichtner, Technische
Universität Dresden, 2018) with 24 h of infiltration, followedby 144 h of redistribution and a hydraulic loading rate of
300 m yr−1 . Infiltration water (0.54 L min−1 ) was takenfrom the Elbe River with natural water quality fluctuations
during the experiment. The concentration of total suspended
solids in river water fluctuated between 7 and 40 mg L−1 ,whereas the dissolved organic C concentration stayed almost
constant at 30 mg L−1 (5 mg L−1 in the Elbe River water and
25 mg L−1 added). At the beginning of the scenario, the sandunderneath the infiltration basin to a depth of 26 cm was
replaced by new sand with the same initial soil properties, as
indicated in Figure 1.
Tracer breakthrough curves using NaCl were used todetermine the reduction of hydraulic conductivity with time
during the laboratory experiment. The tracer tests were6 of 19GLASS ET AL.

Vadose Zone Journal

F I G U R E 4 Laboratory tank setup simulating infiltration well recharge (after Bonilla Valverde, 2018)

performed by adding 1 g L−1 NaCl solution to the infiltrationwater to determine the changing median flow velocity(Fichtner et al., 2019). Electrical conductivity monitored atobservation points (for locations see Figure 1) was used todetermine the median arrival time of the tracer (t50 ), fromwhich the hydraulic conductivity was calculated. For further
information, including detailed results and interpretation of
the laboratory experiments, the reader is referred to Fichtner
et al. (2019).

2.2.2 Model setup
A vertical cross-section of the laboratory tank that includes
the infiltration basin was set up in HYDRUS-2D (Figure 3).
Due to assumed symmetry, only half of the transport domain
was simulated. The reservoir boundary condition of the furrow type (Šimůnek et al., 2018) was assigned to the infiltration basin with a bottom width of 22 cm and a side angleof 45◦ . As the use of the reservoir boundary condition prevents the use of the seepage face boundary condition elsewhere in the model domain, the lower boundary was definedas free drainage. This approach is regarded as feasible, since
the tank outflow is not considered in the calibration and does
not influence the measurements at the six observation points

TABLE 2

Soil hydraulic parameters for the laboratory infiltration
well simulations (after Bonilla Valverde, 2018)
Soil

𝛉r

𝛉s

𝛂

n

m−1

Ks
m min−1

Sand exposed to clogging

0.02

0.36

28

1.94

0.186

Sand

0.02

0.36

28

1.94

0.372

Geotextile

0.02

0.36

28

1.94

7.200

Note. θr , residual water content; θs , saturated water content; α and n, empiricalfitting parameters; Ks , saturated hydraulic conductivity.

defined in the HYDRUS model (Figure 3). An Observation
Point 1 was defined at the bottom of the infiltration basin torepresent the observation of the water level in the infiltration
basin. Five additional observation points were defined in the
HYDRUS model that match the locations of sensors in thelaboratory tank (Figure 3) to compare measured and simulated water contents and pressure heads. By using a targetedFE size of 3.1 cm, a mesh refinement of 2.5 cm along soilmaterial boundaries, and a mesh refinement of 1 cm along theinfiltration basin and the left boundary between the infiltration
basin and Observation Point 5, the 2D triangular FE mesh was
generated. The final FE mesh had a total number of 2,656 FE
nodes (Figure 3). The simulation time was 50,400 min (35 d).
The model was divided into five different soil materials,
including the cobble filter at the bottom of the tank, theGLASS ET AL.

Vadose Zone Journal

7 of 19F I G U R E 5 A numerical model setup representing laboratory well injection experiments, and showing assigned boundary conditions, anobservation point for calibration, well radii (rw ) and the distribution of soil materials

sand filter, and three types of experimental sands (Figure 3,
Table 1). The first type was the sand, which does not get
replaced (Soil Material 3), the second type was the sand,which gets replaced at the beginning of each scenario (SoilMaterial 2), and the third type was the sand exposed toclogging in the first 2.5–3.0 cm below the infiltration basin(Soil Material 1). The laboratory experiments showed thatthe thickness of the clogging layer always ranged between 2
and 3 cm. The initial soil hydraulic parameters were definedusing the results of the HYPROP soil laboratory analysis
for Soil Materials 1–3 and subsequently calibrated using
HYDRUS-2D inverse simulations (Table 1). As observation
data, tensiometer measurements at five observation points(2–6, Figure 3) were used. The soil hydraulic parameters forSoil Materials 4 and 5 were defined using texture-related
values, as simulations showed that the model results were
not sensitive to these parameters. For Soil Material 1, thetime-variable scaling factor λ for the hydraulic conductivity

reduction was inversely calibrated using water level measurements in the infiltration basin. The soil hydraulic parametersof Soil Materials 1 and 2 were identical, with the exceptionof the hydraulic conductivity of the surface soil layer (SoilMaterial 1), which had to be reduced by a factor of 15
compared with the replaced soil (Soil Material 2) due to thehigh initial water level during the first infiltration cycle.

2.3

Recharge well experiment2.3.1 Experimental setup
The experimental laboratory setup is based on the work of
Bonilla Valverde (2018) (Figure 4), who tested the influenceof the screen length on the injection rate in an aquifer–wellsystem representing an isotropic, homogeneous, uniform soil,
and an infinite aquifer. The aquifer–well system consisted of afiberglass cylinder 1.2 m high with an internal radius of 0.5 m

8 of 19GLASS ET AL.

Vadose Zone Journal

F I G U R E 6 Schematic of the model domain for hypotheticalsimulations to study the influence of the hydraulic conductivity
reduction on the extent of the subsurface infiltration area for domains
with different heights, widths, well depths, well radii (rw ), and initialwell water levels (hw )TABLE 3

Soil hydraulic parameters used in hypotheticalsimulations to quantify the extent of the infiltration area
Soil

𝛉r

𝛉s

𝛂

n

m−1

Ks

from the bottom of the tank were active, and a high-headbuildup in the well could be observed. The first part of theexperiment involved the operation of the aquifer–well system
as a closed system, using pure tap water where outflow waterflows back to the water storage tank. This ensured that thetank reached equilibrium conditions when the water table in
the well (22.4 cm) and outflow remained constant. This partof the experiment was used to validate the general numericalmodel setup (see Section 2.3.2). The second part of theexperiment involved physical clogging of the laboratory tank.
For this, Friedländer blue clay with a concentration of 200 mg
L−1 was added to the water storage tank. The Friedländer
blue clay has the following particle size distribution, which
was determined by the sieve and sedimentation analysis
according to DIN18312: <2.0 μm = 43%, 2.0– 6.3 μm =22.75%, 6.3–20.0 μm = 15.25%, 20.0–63.0 μm = 5.5%,
63.0–200.0 μm = 6.5%, and >200 μm = 7% (Schüler, 2016).
Pumps, which mixed water in the water storage tank and theconstant head tank, ensured the suspension of clay. Outflow
pipes were diverted to a water discharge tank so that thequality of inflow water was kept constant during the second
part of the experiment. A pressure logger positioned at thebottom of the well was used to measure the water level inthe well automatically every 5 min, which was corrected forchanges in the air pressure, which were measured using abarometer. Water levels in piezometers away from the well
were measured manually during the experiment.

m min−1

Loam

0.078

0.43

3.6

1.56

0.0104

Sandy loam

0.065

0.41

7.5

1.89

0.0442Note. θr , residual water content; θs , saturated water content; α and n, empiricalfitting parameters; Ks , saturated hydraulic conductivity.

and various piezometers at the bottom of the tank to measure
the water head (Figure 4). The fully penetrating smalldiameter, high-density polyethylene (HDPE) well screen witha radius of 12.7 mm in the center of the aquifer tank was used
to inject water. A constant head tank ensured a constant infiltration rate throughout the experiment. A series of four side
perforations (diameter = 0.05 m) positioned at a 90◦ angle
at five depths served as an outflow system and ensured radial
flow. A geotextile was installed between the tank wall and theaquifer material as a drainage system to maximize radial flow.
The lowest outflow level was located 0.1 m above the tankbottom (Figure 4). The aquifer–well system was filled with amedium to coarse sand with a saturated hydraulic conductivity Ks of 0.372 m min−1 and a bulk density of 1.5 g cm−3 .
The buildup of the pressure head in the well due toclogging of the laboratory tank was examined using the fully
open screen length, where infiltration took place along thewhole well wall. The infiltration rate was kept at 900 ml
min−1 so that only the lowest outflow perforations at 0.1 m

2.3.2 Model setup
A two-dimensional axisymmetrical model with a total size of
0.5 × 1.1. m was set up in HYDRUS-2D to simulate the wellinjection experiment (Figure 5). A previous numerical studyof the laboratory aquifer–well system by Bonilla Valverde(2018) used the HYDRUS-3D model and approximated thewell with the help of a highly permeable soil with a constanthead. In contrast with that, the well reservoir boundary
condition (Šimůnek et al., 2018) was considered in this study
at the inflow boundary along the well casing. This has theadvantage that infiltration along the well wall is dependent on
the head in the well. The well had a radius of 0.0127 m, withan initial water level in the well of 0.23 m and a recharge rateof 0.0009 m3 min−1 (Figure 5). For inverse simulations, one
observation point was defined at the well bottom. The tankoutflow was defined as a constant head boundary condition.
During the clogging experiment, the head at the outflow
boundary stayed constant (0.127 m), as measured in thelaboratory experiment. The model was spatially discretized
into a 2D triangular FE mesh with a targeted FE size of
0.03 m and the mesh refinement of 0.01 m along the reservoir
and the constant head boundaries, resulting in a total numberof 4,212 nodes (Figure 5). The mass balance error in all

GLASS ET AL.

Vadose Zone Journal

9 of 19F I G U R E 7 Simulated (sim) and observed (obs) (a) pressure heads and (b) water contents at three observation points defined in Figure 3 (R2 =.96, weighted RMSE = 0.0083)

simulated scenarios was always considerably below 1%,
which is considered acceptable and assures the quality of the
FE mesh (Brunetti, Šimůnek, Turco, & Piro, 2017).
The soil hydraulic parameters for the experimental sandand the geotextile were taken from Bonilla Valverde (2018)(Table 2). To consider clogging close to the well casing,
the experimental sand was divided into two parts. The firstpart represents the sand close to the well casing where thehydraulic conductivity changes over time, and the clogging
factor is activated. For this soil type, the initial saturated
hydraulic conductivity was halved compared with Bonilla
Valverde (2018) due to the resistance of the HDPE well screen
to replicate the initial well water level. The second part covers
the rest of the tank, where the hydraulic conductivity is kept
constant during the simulations (Figure 5). To represent
progressive clogging along the well casing, the clogging area
was defined as a cone with a decreasing thickness from thebottom (8.5 cm) to the top of the well casing (2 cm). Thisis necessary only for Soil Material 1, in which the clogging
factor λ and, as a consequence, the temporal reduction of
the hydraulic conductivity was considered. However, thisapproach has the disadvantage that clogging of the soil starts
immediately from the beginning of the simulation in the

F I G U R E 8 Simulated vs. observed water levels in the infiltration
basin as a function of time (hydraulic conductivity reduction [λ] = 1.9
× 10−5 min−1 )entire defined clogging area and is not dependent on the waterlevel in the well in the numerical model. The clogging factor
λ is calibrated using the measured water head in the well.

2.4

Hypothetical experiments

In this section, the influence of the hydraulic conductivity
reduction over time on the infiltration volumes and thehorizontal extent of the saturated infiltration area around10 of 19

GLASS ET AL.

Vadose Zone Journal

F I G U R E 9 Comparison of normalized soil hydraulicconductivities (Kt /Kstart ) of the upper soil layer determined using tracer
breakthrough curves (K experiment) and calibrated with the help of thewater level measurements in the infiltration basin (K simulation,
hydraulic conductivity reduction [λ] = 1.9 × 10−5 min−1 )F I G U R E 10 Results of the laboratory well injection experiment:(a) the water level in the well as a function of time measured manuallyand using the pressure logger, and (b) cross-sections of manually
measured water heads at the beginning and end of the clogging
experimentvadose zone recharge wells for various well system setups
was examined. Depending on the scenario, the model domain
had a height of 5–90 m and a width of 15–30 m (Figure 6).
The simulation time for each scenario was 360 h (i.e., 15 d).
A reservoir boundary condition was specified in the wellwith initial well water levels of 0.2–1 m. A free drainage
boundary condition was assigned at the bottom of the model

domain. The dependence of the infiltration volume (over theduration of the simulation; i.e., 360 h) and the size of theinfiltration area on the extent of clogging was analyzed by
performing simulations with a clogging factor λ between
0 and 0.1 h−1 . Various well configurations (a well depth
between 2 and 30 m, a well radius between 0.25 and 1.2 m),
injection rates (defining the inflow into the well) between
0.006 and 1 m3 h−1 , and two different soils (Table 3) were
simulated to evaluate how the infiltration volume and thehorizontal infiltration extent (HIE) depends on the hydraulicconductivity change over time. The extent of the HIE wasdetermined by finding the farthest point of the zero pressure
head isoline from the infiltration well after 360 hours of
infiltration. The scenarios were chosen to represent a rangeof dry well setups from small-scale recharge systems aspresented by Händel et al. (2016) and Händel, Liu, Dietrich,
Liedl, and Butler (2014) up to large-scale systems as used,
for example, in Tucson, AZ (Edwards et al., 2016).

3

RESULTS AND DISCUSSIONS

3.1

Infiltration basin

The simulation of the laboratory experiment with cyclic
infiltration of Elbe River water and subsequent drying phases
focused on the replication of the water level in the infiltration
basin. During both wetting and drying phases of the experiment, time series of pressure heads and water contents werereproduced by the numerical model reasonably well (Figure 7). The overall fit between the simulated and observed
measurements is good (R2 = .96, weighted RMSE = 0.008),
despite differences in the maxima and minima of measured
and simulated water contents and pressure heads. The R2
and RMSE statistical measures are calculated for weighted
deviations between measured and simulated water contents
and pressure heads, whereas measurement variances are used
as weights to minimize differences in weighting between
different data types (Šimůnek et al., 2013). Deviations could
be caused by differences in the placement of measurement
devices within the laboratory tank (various positions of watercontent and pressure head sensors at the same depth), by local
soil heterogeneities due to the installation of measurement
devices, or by nonhomogeneous packing of the laboratory
tank, as well as by preferential flow. The replacement of thesoil material at the beginning of the scenario underneath theinfiltration basin could have influenced the tensiometer andwater content measurements in the tank.
For Soil Material 1, which is located in the first 2–3 cm
underneath the infiltration basin, the clogging factor λ was
activated to reproduce the reduction in the soil hydraulicconductivity and the subsequent water level rise in theinfiltration basin. To adequately reproduce the initial water

Sandy loam

Loam

Loam

Loam

Sandy loam

Sandy loam

Sandy loam

Sandy loam

Loam

Loam

Sandy loam

Sandy loam

Loam

Loam

Loam

Sandy loam

Sandy loam

Loam

Loam

Loam

Sandy loam

Sandy loam

Loam

Loam

Loam

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

Sandy loam

Sandy loam

6

7

Loam

Sandy loam

4

5

Sandy loam

Loam

2

Sandy loam

1

3

Soil

20

m

0.250

0.250

0.250

0.250

0.250

0.750

0.750

0.750

0.750

0.750

0.500

0.500

0.500

0.500

0.500

0.250

0.250

0.250

0.250

0.500

0.500

0.500

0.250

0.250

0.250

0.250

0.250

0.250

0.250

0.075

0.075

0.075

Well radius

0.060

0.030

0.010

0.100

0.060

0.100

0.030

0.010

0.100

0.500

0.100

0.050

0.010

0.100

0.500

0.030

0.010

0.100

0.030

0.100

0.010

0.010

0.030

0.010

0.100

0.050

0.010

0.006

0.006

0.006

0.050

0.010

m3 h−1

Injection rate

23.69

11.63

4.25

43.74

24.83

36.39

11.05

4.17

42.02

178.34

37.61

18.68

3.95

42.05

179.96

11.41

4.07

42.13

13.72

42.10

3.70

4.06

11.54

4.05

43.52

22.00

3.65

2.21

2.44

2.44

18.72

3.61

23.39

11.56

4.31

43.41

23.26

35.63

10.66

3.43

41.70

179.36

37.41

18.35

3.96

41.92

180.68

11.29

4.14

42.21

12.78

41.90

3.63

4.03

11.46

4.08

43.39

22.14

3.61

2.17

2.48

2.63

19.01

3.59

m3

23.50

11.42

4.14

42.20

21.48

35.42

10.31

3.35

42.19

177.28

37.09

18.19

3.60

43.47

179.69

11.24

3.92

42.52

11.58

43.42

3.57

3.85

11.27

4.11

43.63

21.84

3.60

2.17

2.13

2.24

19.66

3.59

4.11

1.90

0.46

1.60

0.82

6.12

1.31

0.11

1.11

7.90

6.55

3.03

0.22

1.32

8.31

1.90

0.47

1.60

0.31

1.35

0.00

0.25

1.83

0.47

1.61

0.67

0.00

0.00

0.16

0.33

0.66

0.01

4.88

2.74

1.29

2.39

1.36

6.98

1.85

1.02

2.24

9.70

7.18

3.77

1.10

2.19

9.54

2.44

1.11

2.24

0.89

2.30

0.33

1.02

2.46

0.99

2.18

1.23

0.35

0.28

0.65

0.68

1.01

0.30

𝛌 = 0.01 h−1Well water level
𝛌 = 0.1 h−1

𝛌 = 0 h−1

𝛌 = 0.01 h−1

Cumulative infiltration
𝛌 = 0 h−1

5.11

2.83

1.40

2.54

1.49

7.37

2.80

1.26

2.39

9.76

7.42

3.97

1.19

2.36

9.64

2.59

1.22

2.35

0.98

2.44

0.42

1.08

2.58

1.11

2.26

1.30

0.39

0.33

0.68

0.75

1.14

0.34

𝛌 = 0.1 h−1

HIE

4.76

2.50

0.81

2.03

1.12

6.63

4.01

0.87

2.31

9.95

6.67

2.49

0.76

2.30

9.99

2.47

0.81

2.29

0.63

2.32

0.00

0.75

2.50

0.82

2.03

1.08

0.00

0.00

0.44

0.43

0.64

0.08

m

𝛌 = 0 h−1

4.60

2.45

1.08

2.07

1.07

6.42

3.78

1.46

2.66

9.95

6.52

2.12

1.25

2.36

9.95

2.38

1.03

2.30

0.85

2.45

0.53

0.91

2.52

1.03

2.04

1.26

0.32

0.28

0.54

0.51

0.66

0.15

𝛌 = 0.01 h−1(Continues)

4.78

2.57

1.14

2.12

1.13

6.60

4.05

1.67

2.79

10.06

6.64

3.01

1.32

2.50

10.04

2.50

1.11

2.33

0.89

2.56

0.70

1.13

2.60

1.12

2.08

1.30

0.42

0.35

0.65

0.52

0.73

0.18

𝛌 = 0.1 h−1

Vadose Zone Journal

20

20

20

20

10

10

10

10

10

10

10

10

10

10

10

10

10

10

5

5

5

5

5

5

2

2

2

2

2

2

2

Well depth

The horizontal infiltration extent (HIE) after 360 h of infiltration for different well setups and different scaling (clogging) factors, λ (0, 0.01, and 0.1 h−1 )No.

TABLE 4GLASS ET AL.
11 of 19

Sandy loam

Sandy loam

Sandy loam

55

56

57

Loam

Loam

53

54

Sandy loam

Loam

51

52

Sandy loam

Sandy loam

49

50

Loam

Loam

47

48

Sandy loam

Sandy loam

45

46

Loam

Sandy loam

43

44

Loam

Loam

41

42

Sandy loam

Sandy loam

39

40

Loam

Loam

37

38

Loam

Loam

35

36

Sandy loam

Sandy loam

33

30

30

30

30

30

30

30

30

30

30

30

30

30

30

30

30

30

20

20

20

20

20

20

20

20

Well depth(Continued)

34

Soil

No.

TABLE 4

1.200

1.200

1.200

1.200

1.200

1.200

0.750

0.750

0.750

0.750

0.750

0.500

0.500

0.500

0.500

0.500

0.500

0.750

0.750

0.750

0.750

0.500

0.500

0.500

0.500

Well radius

1.000

0.500

0.100

0.500

0.100

0.050

0.800

0.300

0.100

0.300

0.030

0.900

0.300

0.100

0.250

0.100

0.030

0.300

0.100

0.100

0.030

0.050

0.010

0.300

0.100

Injection rate

352.56180.3342.38

105.53

35.81

18.96

291.07116.22

44.27

113.31

11.74

328.57117.42

44.27

97.32

38.28

11.59

116.29

44.04

37.44

11.70

23.69

4.78

115.9843.76

348.24176.2341.39

103.41

34.38

17.41

292.20115.81

43.48

112.3710.99

327.50116.89

43.61

97.02

37.40

11.12

115.85

42.73

36.70

11.18

18.04

4.07

115.4143.03

348.59177.0239.36

102.03

33.36

15.93

288.03116.7441.75

111.859.88

328.31116.25

42.36

96.98

37.72

11.07

115.28

35.12

36.28

11.07

17.73

3.87

114.77

35.11

14.95

6.97

0.68

16.77

5.54

2.32

12.27

4.81

1.08

17.87

1.31

14.56

4.93

1.33

15.67

6.47

1.56

4.84

1.07

6.20

1.32

3.07

0.29

5.11

1.31

18.21

9.85

3.38

18.82

6.89

3.89

14.12

6.43

2.80

19.22

2.99

15.28

6.18

2.55

16.71

7.66

2.85

6.03

2.11

7.17

2.53

3.82

1.10

6.17

2.36

𝛌 = 0.01 h−1Well water level
𝛌 = 0.1 h−1

𝛌 = 0 h−1

𝛌 = 0.01 h−1

Cumulative infiltration
𝛌 = 0 h−1

19.08

10.49

3.85

19.68

7.74

4.75

14.19

6.89

3.15

20.06

3.39

15.90

6.66

3.01

17.48

8.29

3.40

6.35

2.17

7.66

3.02

4.27

1.49

6.46

2.19

𝛌 = 0.1 h−1

HIE

19.46

10.84

2.26

12.42

6.82

4.20

15.56

6.64

2.06

12.68

2.53

17.22

6.67

2.08

11.58

6.82

2.50

6.70

2.05

6.77

2.55

4.04

0.86

6.64

2.03

𝛌 = 0 h−1

19.34

10.16

3.51

12.13

6.39

3.91

15.32

6.39

2.61

12.40

2.74

16.90

6.41

2.28

11.40

6.50

2.47

6.57

2.23

6.51

2.60

3.72

1.11

6.55

2.19

𝛌 = 0.01 h−1

18.60

10.88

3.89

12.48

6.94

4.63

15.50

6.79

2.84

12.57

3.10

17.17

6.54

2.51

11.61

6.86

2.92

6.69

2.19

6.82

3.07

4.03

1.39

6.66

1.96

𝛌 = 0.1 h−1

12 of 19
Vadose Zone Journal

GLASS ET AL.

GLASS ET AL.

Vadose Zone Journal

13 of 19F I G U R E 11 Simulated pressure heads in the aquifer–well system after 320 min when the effect of clogging on the hydraulic conductivity(a) was not and (b) was considered (using the time-variable scaling factor with hydraulic conductivity reduction [λ] = 0.029 min−1 ) in the simulationlevel in the infiltration basin, the initial saturated hydraulicconductivity of Soil Material 1 was defined as 15 times
lower than for Soil Material 2. The low initial saturated soil
hydraulic conductivity of Soil Material 1 could be due to anincreased infiltration impedance of the soil surface caused by,for example, the air pressure, the low water content in the soilat the beginning of the experiment, air entrapment, or water
repellency.
The best fit between simulated and observed water levelsin the infiltration basin was reached by setting the clogging
factor λ to 1.9 × 10−5 min−1 for the first soil layer (Figure 8).
The reduction of the saturated hydraulic conductivity isapproximated by an exponential function, which results in acontinuous increase in the water level in the infiltration basinin subsequent infiltration cycles in the numerical simulations.
In contrast with that, the water level measured in the laboratory experiment did not increase continuously with infiltration
cycles and even slightly decreased from Infiltration Cycle

2 to 3 (Figure 8). The air in soil pores may have caused adecrease in the infiltration capacity during Infiltration Cycle
2. As the water level in the infiltration basin dropped during
Infiltration Cycle 3 by ∼1 cm, the air seemed to have been
displaced by infiltrating water. This is supported by an abrupt
water level drop during the third infiltration cycle.
Temporal reductions in the hydraulic conductivity determined using tracer breakthrough curves and the calibration of
the clogging factor in HYDRUS were of the same magnitude(Figure 9). At the end of the experiment, the hydraulicconductivity was reduced by 60% in the upper soil layer
due to clogging. Although the reduction in the saturated
hydraulic conductivity was the same as that determined bymeasurements at the end of the simulation period, the model
underestimated the hydraulic conductivity reduction during
the first 20 d of the experiment (Figure 8). The measured
water level in the infiltration basin is higher in the secondand early part of the third infiltration cycle compared with

14 of 19

GLASS ET AL.

Vadose Zone Journal

F I G U R E 12 Water levels in the well as a function of time measured during the laboratory experiment and simulated while either considering(with hydraulic conductivity reduction [λ] = 0.029 min−1 ) or not considering clogging

the simulation, which is in agreement with the measured
hydraulic conductivity reduction determined by tracer tests.

3.2

Recharge well

The well water level in the laboratory aquifer–well systemincreased from 22 to 103 cm during 320 min of infiltration
with a suspension with blue clay (Figure 10a). The water head
only rose in the well and its close vicinity. The water head
10 cm away from the well was not influenced by physicalclogging (Figure 10b).
The numerical simulation of the laboratory experimentthat considered the time-variable scaling factor resulted in amore realistic pressure head distribution in the model domain(Figure 11). The water level in the well and the pressure headclose to the well were higher when clogging was accounted
for in the simulation using the time-variable clogging factor.
In agreement with the laboratory experiments, clogging
influenced the pressure head only in the close vicinity of thewell (Figure 10b, 11b).
By considering the time-variable scaling factor in the
numerical simulations, the clogging-induced water level risein the infiltration well can be approximated well (Figure 12,
R2 = .98, RMSE = 0.024). Small fluctuations of the waterlevel in the well, especially visible in the simulation without clogging and at the beginning of the simulation withclogging, were due to the used iteration criteria (pressure
head and water content tolerances). Especially during thefirst part of the experiment, the simulated water level rise in

the well is in good agreement with the measured water levelrise. During the second part of the experiment, observed
and simulated water levels start deviating, likely due tothe exponential equation used for the time-variable scaling
factor in combination with the defined clogging area in theupper part of the well not describing exactly the conductivity
reduction over time due to clogging.

3.3

Hypothetical simulations

Table 4 lists all simulated well configurations and summarizes various simulated results, such as the final water levelsin the well, cumulative infiltration volumes, and the final
HIEs for different clogging factors (λ = 0, 0.01, and 0.1 h−1 ),and an increase in HIE when clogging is or is not considered(Figures 13 and 14). Different configurations involve twodifferent soils (sandy loam and loam), five different welldepths (2, 5 10, 20, and 30 m), different injection rates (0.006,
0.01, 0.03, 0.05, 0.06, 0.1, 0.3, 0.5, 0.8, 0.9, and 1.0 m3 h−1
for sandy loam, and 0.01, 0.03, 0.05, 0.1, 0.25, 0.3, and 0.5
m3 h−1 for loam), and different initial water levels in thewell (0.2 and 1.0 m). There are very complex relationships
between all of these variables due to the nonlinearity offlow in the vadose zone and interactions of various fluxes(injection and infiltration) in the well.
Table 4 reports only cumulative infiltrations from thewell rather than groundwater recharge. To be able to reportgroundwater recharge, additional assumptions about the
thickness of the vadose zone between the infiltration system

GLASS ET AL.

Vadose Zone Journal

15 of 19F I G U R E 14 Change in the horizontal extent of infiltration (HIE,%) after 360 h of infiltration as a function of the clogging factor (whenhydraulic conductivity reduction [λ] changes from 0 to 0.01 h−1 or
0.1 h−1 ) and (a) the injection rate (m3 h−1 ) and (b) final water level inthe well (m)F I G U R E 13 Infiltration volumes (m3 ) after 360 h of infiltration
as a function of the clogging factor (hydraulic conductivity reduction[λ] = 0, 0.01, and 0.1 h−1 ) and (a) the injection rate (m3 h−1 ), (b) the
final water level in the well (m), and (c) the well radius (m)and the groundwater table, and the homogeneity or heterogeneity of the vadose zone, would have to be made. For
the system analyzed in these hypothetical simulations, it ispossible to assume that there are only minimal losses due toevaporation (due to infiltration occurring at depths far fromthe soil surface) and that infiltrating water will eventuallyreach the groundwater. There are, however, complex tempo-

ral relationships between well infiltration and groundwater
recharge (depending mainly on the thickness of the vadosezone and its heterogeneity), and these relationships were thesubject of a study by Sasidharan et al. (2020) and are beyond
the scope of this manuscript.
Interestingly, there are only relatively minor differences incumulative infiltration volumes due to clogging for particularwell configurations (Table 4). Increased clogging (higher λ)leads to a reduction in the infiltration capacity for a given
water level in the well. This reduced infiltration capacity (andreduced infiltration) results in an increase in the water volume(and thus the water level) in the well. This increase in thewater level in the well then increases the infiltration area andcapacity of the well and compensates for the effects of clogging. The water level in the well thus automatically responds
to the reduction in the infiltration capacity due to clogging
and compensates for it. This compensation can occur inrecharge wells where an increase in the water level leads to a

16 of 19

GLASS ET AL.

Vadose Zone Journal

corresponding increase in the infiltration area. This would not
happen in infiltration basins, where an increase in the waterlevel leads to only a minimal increase in the infiltration area.
Figure 13 shows cumulative infiltration volumes after 360 hof infiltration as a function of the clogging factor (λ = 0.01 h−1and λ = 0.1 h−1 ), the injection rate, the final water level in thewell, and the well radius. As discussed above, the infiltration
volume is not significantly dependent on the clogging factor
and the injection rate. The same injection rates result in thesame infiltration volumes for different clogging factors (Figure 13a) by the system automatically adjusting the water levelin the well (Table 4, Figure 13b). Figure 13b shows that higher
water levels in the well (see the rightward trend from blue
symbols [no clogging] to green and red symbols [significant
clogging]) are needed to infiltrate the same water volumes forthe same well configurations. Finally, Figure 13c shows, as
expected, that larger well radii can infiltrate larger water volumes independently of the clogging factor, and that the cumulative infiltration volume depends mainly on the injectionrate.
The HIE varies depending on the well depth and radius,
the injection rate, and the reduction clogging factor (Table 4).
Figure 14 shows how much HIE (in %) changes after 360 h
for different injection rates and final water tables in the wellwhen the clogging factor increases from 0 (no clogging)to 0.01 h−1 (moderate clogging) or 0.1 h−1 (significant
clogging). For low injection rates and diameters, HIE (whichis relatively small) increases when the reduction in the soil
hydraulic conductivity over time in the vicinity of the wellis considered (Figure 14a). This increase is insignificant forlarger injection rates. Similarly, low water levels in the wellresult in an increase in the HIE when the reduction in thehydraulic conductivity over time is considered (Figure 14b).
Again, this increase is relatively insignificant for higher
water levels in the well. Larger clogging factors (λ) generally
result in a greater expansion of the HIE. Due to clogging, thewater level in the well and the infiltration surface increases,
leading to a broader spread of the infiltration water in thehorizontal direction. However, this effect is important only
for low injection rates and low water levels in the well andcan be neglected for larger injection rates and higher watertables.
To summarize, the effect of the reduction of the hydraulicconductivity in the vicinity of the recharge well on theinfiltration volume is relatively small. This effect is similarly
small on HIE for larger injection rates and is significant
only for relatively small injection rates when HIE is alsoquite small. This relatively small effect of clogging on
infiltration from recharge wells, which can compensate theeffect of clogging by increasing the water level in the welland thus the infiltration area, is unlikely to be found for
infiltration basins, where such compensation is unlikely
to occur.

4

CONC LU SI ON S

The simulation of laboratory experiments using modified
HYDRUS-2D demonstrates the successful implementationof the time-variable scaling factor to represent hydraulicconductivity changes over time, which is especially relevant for the simulation of MAR schemes that are prone toclogging.
Cyclic infiltration of Elbe River water in an infiltration basin and increasing clogging of the surface layer,
mainly caused by physical and biological processes, were
well reproduced by the HYDRUS-2D simulation, which
considered the time-variable scaling factor in combination
with the “furrow”-type reservoir boundary condition. Water
content and pressure head changes during the wetting
and drying phases of the experiment could be, in general,
well replicated by the model. The water level rise in theinfiltration basin during subsequent infiltration cycles was
reproduced, and the hydraulic conductivity reduction inducedby the calibrated clogging factor λ corresponded well to themeasured hydraulic conductivity reduction in the laboratory
experiment. The numerical simulation of the aquifer–wellexperiment was similarly improved by the inclusion of thescaling factor reducing the soil hydraulic conductivity as
the increase of the water level in the well due to physicalclogging was accurately reproduced when considering thesoil hydraulic conductivity reduction over time.
By simulating two laboratory experiments representing
basin and well recharge, it was demonstrated that the representation of clogging in HYDRUS-2D using the time-variablescaling factor is feasible. The simulations indicate that the
time-variable scaling factor, especially in combination withthe reservoir boundary condition, can improve simulationsof clogging in MAR schemes by HYDRUS. The current
setup is only valid if the hydraulic conductivity reduction
is in the form of an exponential function. Although thismay not always be the case, the literature indicates thatthe exponential function is frequently used to approximate
clogging (Pérez-Paricio, 2001). One also must consider thatthe numerical model that uses the scaling factor to reduce thehydraulic conductivity is only valid for a limited time period,
which in general covers the time between maintenance intervals. At some point, either the infiltration structure overflows,or the soil material gets fully clogged (Ks approaches 0),
which usually leads to the secession of infiltration and theapplication of maintenance measures such as scraping of theinfiltration surface in the infiltration basins or backflushingof injection wells.
Hypothetical simulations indicate that clogging has arelatively small effect on infiltration volumes from rechargewells, which can compensate the clogging effect by adjusting(increasing) the water level in the well, and thus the wellinfiltration area. Clogging can influence the HIE in the

GLASS ET AL.

subsurface of vadose zone wells, depending on the well configuration and operation, only for relatively small injection
volumes. As a consequence, clogging is likely to influence
the infiltration capacity of certain MAR systems (e.g.,
infiltration basins) more than that of others (e.g., rechargewells), and this should therefore be considered during thedesign and planning process.
As recognized by Martin (2013), each MAR system
undergoes some kind of clogging during its operational
lifetime, and thus clogging represents an important issue
determining the success or failure of a system. In this way, thepresented approach can be applied to numerically evaluate
the resulting infiltration capacity, and it can furthermore help
to improve the design and operation of MAR facilities withregard to clogging. The scaling factor reducing the hydraulicconductivity over time only represents an initial, very simple
approximation to simulate clogging in a MAR system.
The approach can be further enhanced by including,for example, multiple infiltration structures to study theinteraction among infiltration structures, various clogging
factors to account for different clogging rates in various areasof the simulation domain, or the dependence of clogging
on the water content or the water level in the infiltration
structure. If needed, the equation used for the hydraulicconductivity reduction (Equation 4) can be expanded toinclude an asymptotic term, so that Ks approaches a certain
value over time instead of 0, as described in Pérez-Paricio(2001). Furthermore, the implementation of other functions
to represent clogging could make the developed approach
more flexible. Also, clogging should not only be related tochanges in hydraulic conductivity over time but also otherrelevant soil hydraulic parameters such as soil porosity or
the van Genuchten soil retention parameters. This could beachieved by including, for example, the Kozeny–Carman
relationship in the simulations, as demonstrated by Xie
et al. (2015). To further enhance simulation capabilities,
clogging processes could be simulated not only as a functionof the soil type, but also the quality of infiltrating water.
For that, the application of more complex reactive transportmodels such as HP1 (Jacques, Šimůnek, Mallants, & van
Genuchten, 2018), MIN3P (Mayer, Frind, & Blowes, 2002),
or PHREEQC (Parkhust & Appelo, 1999) may be necessary.
ACKNOW LEDGMENTS
This study was supported by the German Federal Ministryof Education and Research (BMBF), Grant no. 01LN1311A(Junior Research Group “INOWAS”). The authors would
like to thank Thomas Fichtner (Junior Research Group“INOWAS”, Technische Universität Dresden) for providing the measurements for the laboratory infiltration basin
experiment (including tensiometer measurements, watercontent, the water level in the infiltration basin, and tracer
breakthrough curves). Also, the authors would like to thankVadose Zone Journal

17 of 19

José Bonilla Valverde (Instituto Costarricense de Acueductosy Alcantarillados, San José, Costa Rica) for the construction and initial setup of the aquifer–well system, and Fritz
Kalwa (Institute of Groundwater Management, Technische
Universität Dresden, Germany) for support in conducting thelaboratory aquifer–well experiment.
CO N F L I C T O F I N T E R E ST
The authors declare no conflict of interest.
O RC I D
https://orcid.org/0000-0003-4517-8091
Jana Glass
https://orcid.org/0000-0001-8250-2749
Catalin Stefan
REFERENCES
Anderson, M., Dewhurst, R., Jones, M., & Baxter, K. (2006). Characterisation of turbidity and well clogging processes in a double porosity
Chalk aquifer during the South London Artificial Recharge Scheme
trials. In Recharge systems for protecting and enhancing groundwater
resources: Proceedings of the 5th International Symposium on Management of Aquifer Recharge ISMAR 5 (pp. 593–598). Paris: United
Nations Education, Scientific and Cultural Organization.
Bonilla Valverde, J. P. (2018). Managed aquifer recharge assessment
to overcome water scarcity during the dry season in Costa Rica
(Doctoral dissertation). Dresden, Germany: Technische Universität
Dresden.
Bouwer, H. (2002). Artificial recharge of groundwater: Hydrogeology
and engineering. Hydrogeology Journal, 10, 121–142. https://doi.
org/10.1007/s10040-001-0182-4
Brunetti, G., Šimůnek, J., Turco, M., & Piro, P. (2017). On the use of
surrogate-based modeling for the numerical analysis of low impact
development techniques. Journal of Hydrology, 548, 263–277.
https://doi.org/10.1016/j.jhydrol.2017.03.013
Dillon, P. (2005). Future management of aquifer recharge. Hydrogeology
Journal, 13, 313–316. https://doi.org/10.1007/s10040-004-0413-6
Dillon, P., Stuyfzand, P., Grischek, T., Lluria, M., Pyne, R. D. G., Jain, R.
C., … Sapiano, M. (2018). Sixty years of global progress in managed
aquifer recharge. Hydrogeology Journal, 27, 1–30. https://doi.org/10.
1007/s10040-018-1841-z
Dillon, P. J., Pavelic, P., Page, D., Beringen, H., & Ward, J. (2009). Managed aquifer recharge: An introduction. Canberra: National Water
Commission.
Edwards, E. C., Harter, T., Fogg, G. E., Washburn, B., & Hamad,
H. (2016). Assessing the effectiveness of drywells as tools forstormwater management and aquifer recharge and their groundwater contamination potential. Journal of Hydrology, 539, 539–553.
https://doi.org/10.1016/j.jhydrol.2016.05.059
Fichtner, T., Barquero, F., Sallwey, J., & Stefan, C. (2019). Assessing
managed aquifer recharge processes under three physical model concepts. Water, 11(1). https://doi.org/10.3390/w11010107
Glaß, J. (2019). New advances in the assessment of managed aquifer
recharge through modelling (Doctoral dissertation). Dresden, Germany: Technische Universität Dresden.
Gutiérrez-Ojeda, C., Martínez-Morales, M., & Ortiz-Flores, G. (2007).
Artificial recharge of groundwater in a coal mining area of northeast
Mexico. In P. Fox (Ed.), Management of aquifer recharge for sustainability: Proceedings of the 6th International Symposium on Managed

18 of 19

Vadose Zone Journal

Artificial Recharge of Groundwater, ISMAR6 (pp. 74–83). Phoenix,
AZ: Acacia Publishing.
Händel, F., Liu, G., Dietrich, P., Liedl, R., & Butler, J. J. (2014). Numerical assessment of ASR recharge using small-diameter wells and surface basins. Journal of Hydrology, 517, 54–63. https://doi.org/10.
1016/j.jhydrol.2014.05.003
Händel, F., Liu, G., Fank, J., Friedl, F., Liedl, R., & Dietrich, P. (2016).
Assessment of small-diameter shallow wells for managed aquifer
recharge at a site in southern Styria, Austria. Hydrogeology Journal,
24, 2079–2091. https://doi.org/10.1007/s10040-016-1442-7
Hannappel, S., Scheibler, F., Huber, A., & Sprenger, C. (2014). Recommendations for further data generation. DEMEAU. Retrieved fromhttps://demeau-fp7.eu/system/files/results/M11%20catalogue%
20of%20european%20MAR%20applications%20final.pdf
Jacques, D., Šimůnek, J., Mallants, D., & van Genuchten, M. Th. (2018).
The HPx software for multicomponent reactive transport during
variably-saturated flow: Recent developments and applications. Journal of Hydrology and Hydromechanics, 66, 211–226. https://doi.org/
10.1515/johh-2017-0049
Majumdar, P. K., Sekhar, M., Sridharan, K., & Mishra, G. C. (2008).
Numerical simulation of groundwater flow with gradually increasing
heterogeneity due to clogging. Journal of Irrigation and Drainage
Engineering, 134, 400–404. https://doi.org/10.1061/(ASCE)07339437(2008)134:3(400)Martin, R. editor. (2013). Clogging issues associated with managed
aquifer recharge methods. International Association of Hydrogeologists Commission on Managing Aquifer Recharge. Retrieved fromhttps://recharge.iah.org/files/2015/03/Clogging_Monograph.pdf
Masciopinto, C. (2013). Management of aquifer recharge in Lebanon
by removing seawater intrusion from coastal aquifers. Journal ofEnvironmental Management, 130, 306–312. https://doi.org/10.1016/
j.jenvman.2013.08.021
Mayer, K. U., Frind, E. O., & Blowes, D. W. (2002). Multicomponent
reactive transport modeling in variably saturated porous media using
a generalized formulation for kinetically controlled reactions: Reactive transport modeling in variably saturated media. Water Resources
Research, 38, 13–21. https://doi.org/10.1029/2001WR000862
Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 12,
513–522. https://doi.org/10.1029/WR012i003p00513
Page, D., Vanderzalm, J., Miotliński, K., Barry, K., Dillon, P., Lawrie,K., & Brodie, R. S. (2014). Determining treatment requirements forturbid river water to avoid clogging of aquifer storage and recovery
wells in siliceous alluvium. Water Research, 66, 99–110. https://doi.
org/10.1016/j.watres.2014.08.018
Parkhust, D. L., & Appelo, C. A. J. (1999). PHREEQC (Version 2): A
computer program for speciation, batch-reaction, one-dimensionaltransport, and inverse geochemical calculations. Denver, CO: USGS.
Pavelic, P., Dillon, P. J., Barry, K. E., Vanderzalm, J. L., Correll, R. L.,& Rinck-Pfeiffer, S. M. (2007). Water quality effects on clogging
rates during reclaimed water ASR in a carbonate aquifer. Journal of
Hydrology, 334, 1–16. https://doi.org/10.1016/j.jhydrol.2006.08.009
Pérez-Paricio, A. (2001). Integrated modeling of clogging processes in
artificial groundwater recharge (Doctoral dissertation). Barcelona,
Spain: Technical University of Catalonia.
Radcliffe, D. E., & Šimůnek, J. (2010). Soil physics with HYDRUS: Modeling and applications. Boca Raton, FL: CRC Press.
Rinck-Pfeiffer, S. (2000). Interrelationships between biological, chemical, and physical processes as an analog to clogging in aquifer stor-

GLASS ET AL.

age and recovery (ASR) wells. Water Research, 34, 2110–2118.
https://doi.org/10.1016/S0043-1354(99)00356-5
Rinck-Pfeiffer, S., Dillon, P., Ragusa, S., Hutson, J., Fallowfield, H. et al.
(2013). Reclaimed water for aquifer storage and recovery: A column
study of well clogging. In R. Martin (Ed.), Clogging issues associated with managed aquifer recharge methods (pp. 26–33). Australia:
International Association of Hydrogeologists Commission on Managing Aquifer Recharge.
Ringleb, J., Sallwey, J., & Stefan, C. (2016). Assessment of managed aquifer recharge through modeling: A review. Water, 8(12).
https://doi.org/10.3390/w8120579
Sallwey, J., Glass, J., & Stefan, C. (2018). Utilizing unsaturated soil
zone models for assessing managed aquifer recharge. Sustainable
Water Resources Management, 4, 383–397. https://doi.org/10.1007/
s40899-018-0214-z
Sasidharan, S., Bradford, S. A., Šimůnek, J., DeJong, B., & Kraemer,
S. R. (2018). Evaluating drywells for stormwater management andenhanced aquifer recharge. Advances in Water Resources, 116, 167–
177. https://doi.org/10.1016/j.advwatres.2018.04.003
Sasidharan, S., Bradford, S. A., Šimůnek, J., & Kraemer, S. R. (2019).
Drywell infiltration and hydraulic properties in heterogeneous soil
profiles. Journal of Hydrology, 570, 598–611. https://doi.org/10.
1016/j.jhydrol.2018.12.073
Sasidharan, S., Bradford, S. A., Šimůnek, J., & Kraemer, S. R. (2020).
Groundwater recharge from drywells under constant head conditions. Journal of Hydrology, 583, 124569. https://doi.org/10.1016/
j.jhydrol.2020.124569
Schüler, D. (2016). Entwicklung eines Bewertungsschemas zur Anwendbarkeit von Versickerungsbrunnen im Vergleich zu Infiltrationsbecken (Master’s thesis). Dresden, Germany: Technische Universität
Dresden.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., & van Genuchten, M. Th.
(2013). The HYDRUS-1D software package for simulating the movement of water, heat, and multiple solutes in variably saturated media.
Riverside: University of California Riverside.
Šimůnek, J., Šejna, M., & van Genuchten, M. Th. (2018). New featuresof version 3 of the HYDRUS (2D/3D) computer software package.
Journal of Hydrology and Hydromechanics, 66, 133–142. https://doi.
org/10.1515/johh-2017-0050
Šimůnek, J., van Genuchten, M. Th., & Šejna, M. (2016). Recent developments and applications of the HYDRUS computer software packages. Vadose Zone Journal, 15(7). https://doi.org/10.2136/vzj2016.
04.0033
Sprenger, C., Hartog, N., Hernández, M., Vilanova, E., Grützmacher, G.,
Scheibler, F., & Hannappel, S. (2017). Inventory of managed aquifer
recharge sites in Europe: Historical development, current situationand perspectives. Hydrogeology Journal, 25, 1909–1922. https://doi.
org/10.1007/s10040-017-1554-8
Stefan, C., & Ansems, N. (2018). Web-based global inventory ofmanaged aquifer recharge applications. Sustainable Water Resources
Management, 4, 153–162. https://doi.org/10.1007/s40899-0170212-6
Thompson, A. R., Stotler, R. L., Macpherson, G. L., & Liu, G. (2015).
Laboratory study of low-flow rates on clogging processes forapplication to small-diameter injection wells. Water Resources
Management, 29, 5171–5184. https://doi.org/10.1007/s11269015-1111-y
van Genuchten, M. Th. (1980). A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Science Society

GLASS ET AL.

of America Journal, 44, 892–898. https://doi.org/10.2136/sssaj1980.
03615995004400050002x
Vanderzalm, J. L., Page, D. W., Barry, K. E., & Dillon, P. J. (2013).
Application of a probabilistic modelling approach for evaluation of
nitrogen, phosphorus and organic carbon removal efficiency during four successive cycles of aquifer storage and recovery (ASR)
in an anoxic carbonate aquifer. Water Research, 47, 2177–2189.
https://doi.org/10.1016/j.watres.2013.01.038
Wood, M., Le Gal La Salle, C., Dillon, P., Hutson, J., Fallowfield, H. et al.
2006. Effect of grain size on biological clogging in porous media
leached with synthetic nutrient solutions. In Recharge systems forprotecting and enhancing groundwater resources: Proceedings of the
5th International Symposium on Management of Aquifer Recharge
ISMAR 5 (pp. 630–635). Paris: United Nations Education, Scientificand Cultural Organization.

Vadose Zone Journal

19 of 19

Xie, M., Mayer, K. U., Claret, F., Alt-Epping, P., Jacques, D., Steefel,C., … Simunek, J. (2015). Implementation and evaluation of
permeability-porosity and tortuosity-porosity relationships linked tomineral dissolution-precipitation. Computational Geosciences, 19,655–671. https://doi.org/10.1007/s10596-014-9458-3

How to cite this article: Glass J, Šimůnek J, Stefan
C. Scaling factors in HYDRUS to simulate a reduction
in hydraulic conductivity during infiltration fromrecharge wells and infiltration basins. Vadose Zone J.
2020;19:e20027. https://doi.org/10.1002/vzj2.20027