Scaling factors in HYDRUS to simulate a reduction in hydraulic conductivity during infiltration from recharge wells and infiltration basins

Managed aquifer recharge (MAR) represents a promising technique to cope with increasing 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 to quantify and assess, and their simulations in MAR schemes have so far been limited. The variably saturated water flow model HYDRUS‐2D was therefore modified to include time‐variable hydraulic conductivities to more realistically represent clogging at the infiltration interface of infiltration basins and recharge wells. An exponential function 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 the infiltration 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 the processes that lead to clogging at MAR facilities into numerical simulations.

The use of unsaturated water flow models is required for realistic simulations of clogging processes at the infiltration interface 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, Šimůnek, DeJong, & Kraemer, 2018;Sasidharan, Bradford, Šimůnek, & Kraemer, 2019, 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. ( , 2019Sasidharan et al. ( , 2020 have recently used HYDRUS-2D (Radcliffe and Šimůnek, 2010) in combination with the reservoir boundary to assess experimental data collected at two dry wells, as well as the effects of heterogeneity on drywell infiltration and groundwater recharge, without implicitly accounting for clogging.
The main objective of this manuscript is to describe the implementation and validation of a simplified representation of 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 exponential function was implemented into HYDRUS to reduce the hydraulic conductivity with time during simulations. In combination with the reservoir boundary condition, which was recently developed to account for ponded infiltration at MAR facilities (Šimů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 numerical simulations of two laboratory experiments, integrating a time-variable scaling factor to validate the simple representation of clogging at MAR facilities. In addition, hypothetical simulations with varying soils, well configurations, and infiltration rates were conducted to quantify the effects of the reduction of hydraulic conductivity on the infiltration volumes and the extent of the saturated infiltration area around 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 vadose zone, the volume of infiltrated water, and other factors. This relationship between irrigation and groundwater recharge is a subject of the study by Sasidharan, Bradford, Šimůnek, and Kraemer (2020) and is beyond the scope of this manuscript.

Simulation of flow in the unsaturated zone
Variable-saturated water flow was simulated using the HYDRUS-2D model (Šimůnek, van Genuchten, & Šejna, 2016). The program uses the Garlekin-type finite element (FE) scheme to numerically solve the Richards equation for variable-saturated water flow: where θ is the volumetric water content [-], h is the pressure head [L], S is a sink term [T −1 ], x i (i = 1,2) are spatial coordinates [L] where x 1 is the horizontal coordinate and x 2 is the vertical coordinate, K ij A are components of a dimensionless hydraulic conductivity anisotropy tensor K A , K is the unsaturated hydraulic conductivity function given by the product of the relative hydraulic conductivity K r and the saturated hydraulic conductivity K s [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 can be parameterized according to van Genuchten (1980): where θ s is the saturated water content [-], θ r is the residual water content [-], and α [L −1 ], n [-], and m [-] are empirical fitting parameters (Šimůnek, Šejna, Saito, Sakai, & van Genuchten, 2013). The van Genuchten- Mualem equation (van Genuchten, 1980) is used to describe the relationship between unsaturated hydraulic conductivity and soil saturation: where S e 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 the infiltration unit, depending on the infiltration area and infiltration rates (Šimůnek et al., 2018). The utilization of the reservoir boundary condition restricts the simulation to a twodimensional (2D) cross-section of the infiltration structure or its 2D radially symmetric representation. A seepage face boundary 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: where K s,t is the saturated hydraulic conductivity at time t [L T −1 ], K s,initial is the initial saturated hydraulic conductivity [L T −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 the impact of clogging on surface systems (Pérez-Paricio, 2001).

Vadose Zone Journal
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

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 fine material. Above the sand filter, the main experimental soil was installed, which consists of 84 cm of sand (the bulk density of 1.6 g cm −3 ). The particle size distribution is displayed in Figure 2. The tank was equipped with tensiometers and TDR sensors to measure time-variable pressure heads and water contents during infiltration cycles, respectively, as well as with suction cups and electrical conductivity sensors to monitor the water quality ( Figure 1  Note. θ r , residual water content; θ s , saturated water content; α and n, empirical fitting parameters; K s , saturated hydraulic conductivity. infiltration basin was measured using a pressure sensor. The tank was installed in a fully automatic climate tank where the temperature 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) and infiltration 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, followed by 144 h of redistribution and a hydraulic loading rate of 300 m yr −1 . Infiltration water (0.54 L min −1 ) was taken from 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 sand underneath 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 to determine the reduction of hydraulic conductivity with time during the laboratory experiment. The tracer tests were 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 infiltration water to determine the changing median flow velocity (Fichtner et al., 2019). Electrical conductivity monitored at observation points (for locations see Figure 1) was used to determine the median arrival time of the tracer (t 50 ), from which 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).

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 (Šimůnek et al., 2018) was assigned to the infiltration basin with a bottom width of 22 cm and a side angle of 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 defined as 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 Note. θ r , residual water content; θ s , saturated water content; α and n, empirical fitting parameters; K s , saturated hydraulic conductivity.
defined in the HYDRUS model ( Figure 3). An Observation Point 1 was defined at the bottom of the infiltration basin to represent 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 the laboratory tank ( Figure 3) to compare measured and simulated water contents and pressure heads. By using a targeted FE size of 3.1 cm, a mesh refinement of 2.5 cm along soil material boundaries, and a mesh refinement of 1 cm along the infiltration 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 Vadose Zone Journal F I G U R E 5 A numerical model setup representing laboratory well injection experiments, and showing assigned boundary conditions, an observation point for calibration, well radii (r w ) 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 (Soil Material 2), and the third type was the sand exposed to clogging in the first 2.5-3.0 cm below the infiltration basin (Soil Material 1). The laboratory experiments showed that the thickness of the clogging layer always ranged between 2 and 3 cm. The initial soil hydraulic parameters were defined using 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 for Soil 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, the time-variable scaling factor λ for the hydraulic conductivity reduction was inversely calibrated using water level measurements in the infiltration basin. The soil hydraulic parameters of Soil Materials 1 and 2 were identical, with the exception of the hydraulic conductivity of the surface soil layer (Soil Material 1), which had to be reduced by a factor of 15 compared with the replaced soil (Soil Material 2) due to the high initial water level during the first infiltration cycle.

Experimental setup
The experimental laboratory setup is based on the work of Bonilla Valverde (2018) (Figure 4), who tested the influence of the screen length on the injection rate in an aquifer-well system representing an isotropic, homogeneous, uniform soil, and an infinite aquifer. The aquifer-well system consisted of a fiberglass cylinder 1.2 m high with an internal radius of 0.5 m F I G U R E 6 Schematic of the model domain for hypothetical simulations 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 (r w ), and initial well water levels (h w ) Note. θ r , residual water content; θ s , saturated water content; α and n, empirical fitting parameters; K s , saturated hydraulic conductivity.

T A B L E 3 Soil hydraulic parameters used in hypothetical simulations to quantify the extent of the infiltration area
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 with a 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 the aquifer material as a drainage system to maximize radial flow. The lowest outflow level was located 0.1 m above the tank bottom ( Figure 4). The aquifer-well system was filled with a medium to coarse sand with a saturated hydraulic conductivity K s 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 to clogging of the laboratory tank was examined using the fully open screen length, where infiltration took place along the whole well wall. The infiltration rate was kept at 900 ml min −1 so that only the lowest outflow perforations at 0.1 m from the bottom of the tank were active, and a high-head buildup in the well could be observed. The first part of the experiment involved the operation of the aquifer-well system as a closed system, using pure tap water where outflow water flows back to the water storage tank. This ensured that the tank reached equilibrium conditions when the water table in the well (22.4 cm) and outflow remained constant. This part of the experiment was used to validate the general numerical model setup (see Section 2.3.2). The second part of the experiment 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 the constant head tank, ensured the suspension of clay. Outflow pipes were diverted to a water discharge tank so that the quality of inflow water was kept constant during the second part of the experiment. A pressure logger positioned at the bottom of the well was used to measure the water level in the well automatically every 5 min, which was corrected for changes in the air pressure, which were measured using a barometer. Water levels in piezometers away from the well were measured manually during the experiment.

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 well injection experiment ( Figure 5). A previous numerical study of the laboratory aquifer-well system by Bonilla Valverde (2018) used the HYDRUS-3D model and approximated the well with the help of a highly permeable soil with a constant head. In contrast with that, the well reservoir boundary condition (Šimůnek et al., 2018) was considered in this study at the inflow boundary along the well casing. This has the advantage that infiltration along the well wall is dependent on the head in the well. The well had a radius of 0.0127 m, with an initial water level in the well of 0.23 m and a recharge rate of 0.0009 m 3 min −1 ( Figure 5). For inverse simulations, one observation point was defined at the well bottom. The tank outflow 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 the laboratory 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 number of 4,212 nodes ( Figure 5). The mass balance error in all Vadose Zone Journal F 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 (R 2 = .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, Šimůnek, Turco, & Piro, 2017).
The soil hydraulic parameters for the experimental sand and 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 first part represents the sand close to the well casing where the hydraulic 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 the bottom (8.5 cm) to the top of the well casing (2 cm). This is 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, this approach 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 water level in the well in the numerical model. The clogging factor λ is calibrated using the measured water head in the well.

Hypothetical experiments
In this section, the influence of the hydraulic conductivity reduction over time on the infiltration volumes and the horizontal extent of the saturated infiltration area around Vadose Zone Journal F I G U R E 9 Comparison of normalized soil hydraulic conductivities (K t /K start ) of the upper soil layer determined using tracer breakthrough curves (K experiment) and calibrated with the help of the water 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 manually and using the pressure logger, and (b) cross-sections of manually measured water heads at the beginning and end of the clogging experiment vadose 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 well with 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 the duration of the simulation; i.e., 360 h) and the size of the infiltration 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 m 3 h −1 , and two different soils (Table 3) were simulated to evaluate how the infiltration volume and the horizontal infiltration extent (HIE) depends on the hydraulic conductivity change over time. The extent of the HIE was determined 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 range of dry well setups from small-scale recharge systems as presented 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).

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 were reproduced by the numerical model reasonably well (Figure 7). The overall fit between the simulated and observed measurements is good (R 2 = .96, weighted RMSE = 0.008), despite differences in the maxima and minima of measured and simulated water contents and pressure heads. The R 2 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 (Šimůnek et al., 2013). Deviations could be caused by differences in the placement of measurement devices within the laboratory tank (various positions of water content 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 the soil material at the beginning of the scenario underneath the infiltration basin could have influenced the tensiometer and water 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 hydraulic conductivity and the subsequent water level rise in the infiltration basin. To adequately reproduce the initial water

F 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 simulation level in the infiltration basin, the initial saturated hydraulic conductivity 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 an increased infiltration impedance of the soil surface caused by, for example, the air pressure, the low water content in the soil at the beginning of the experiment, air entrapment, or water repellency.
The best fit between simulated and observed water levels in 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 is approximated by an exponential function, which results in a continuous increase in the water level in the infiltration basin in 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 a decrease 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 hydraulic conductivity 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 by measurements 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 second and early part of the third infiltration cycle compared with 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.

Recharge well
The well water level in the laboratory aquifer-well system increased 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 physical clogging (Figure 10b).
The numerical simulation of the laboratory experiment that considered the time-variable scaling factor resulted in a more realistic pressure head distribution in the model domain ( Figure 11). The water level in the well and the pressure head close 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 the well (Figure 10b, 11b).
By considering the time-variable scaling factor in the numerical simulations, the clogging-induced water level rise in the infiltration well can be approximated well (Figure 12, R 2 = .98, RMSE = 0.024). Small fluctuations of the water level in the well, especially visible in the simulation without clogging and at the beginning of the simulation with clogging, were due to the used iteration criteria (pressure head and water content tolerances). Especially during the first part of the experiment, the simulated water level rise in the well is in good agreement with the measured water level rise. During the second part of the experiment, observed and simulated water levels start deviating, likely due to the exponential equation used for the time-variable scaling factor in combination with the defined clogging area in the upper part of the well not describing exactly the conductivity reduction over time due to clogging. Table 4 lists all simulated well configurations and summarizes various simulated results, such as the final water levels in 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 two different soils (sandy loam and loam), five different well depths (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 m 3 h −1 for sandy loam, and 0.01, 0.03, 0.05, 0.1, 0.25, 0.3, and 0.5 m 3 h −1 for loam), and different initial water levels in the well (0.2 and 1.0 m). There are very complex relationships between all of these variables due to the nonlinearity of flow in the vadose zone and interactions of various fluxes (injection and infiltration) in the well. Table 4 reports only cumulative infiltrations from the well rather than groundwater recharge. To be able to report groundwater recharge, additional assumptions about the thickness of the vadose zone between the infiltration system 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 is possible to assume that there are only minimal losses due to evaporation (due to infiltration occurring at depths far from the soil surface) and that infiltrating water will eventually reach the groundwater. There are, however, complex tempo- ral relationships between well infiltration and groundwater recharge (depending mainly on the thickness of the vadose zone and its heterogeneity), and these relationships were the subject of a study by Sasidharan et al. (2020) and are beyond the scope of this manuscript.

Hypothetical simulations
Interestingly, there are only relatively minor differences in cumulative infiltration volumes due to clogging for particular well 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 (and reduced infiltration) results in an increase in the water volume (and thus the water level) in the well. This increase in the water level in the well then increases the infiltration area and capacity 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 in recharge wells where an increase in the water level leads to a corresponding increase in the infiltration area. This would not happen in infiltration basins, where an increase in the water level leads to only a minimal increase in the infiltration area. Figure 13 shows cumulative infiltration volumes after 360 h of infiltration as a function of the clogging factor (λ = 0.01 h −1 and λ = 0.1 h −1 ), the injection rate, the final water level in the well, 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 the same infiltration volumes for different clogging factors (Figure 13a) by the system automatically adjusting the water level in 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 for the 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 injection rate.
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 well when 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 (which is relatively small) increases when the reduction in the soil hydraulic conductivity over time in the vicinity of the well is considered (Figure 14a). This increase is insignificant for larger injection rates. Similarly, low water levels in the well result in an increase in the HIE when the reduction in the hydraulic 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, the water level in the well and the infiltration surface increases, leading to a broader spread of the infiltration water in the horizontal direction. However, this effect is important only for low injection rates and low water levels in the well and can be neglected for larger injection rates and higher water tables.
To summarize, the effect of the reduction of the hydraulic conductivity in the vicinity of the recharge well on the infiltration 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 also quite small. This relatively small effect of clogging on infiltration from recharge wells, which can compensate the effect of clogging by increasing the water level in the well and thus the infiltration area, is unlikely to be found for infiltration basins, where such compensation is unlikely to occur.

CONCLUSIONS
The simulation of laboratory experiments using modified HYDRUS-2D demonstrates the successful implementation of the time-variable scaling factor to represent hydraulic conductivity changes over time, which is especially relevant for the simulation of MAR schemes that are prone to clogging.
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 the infiltration basin during subsequent infiltration cycles was reproduced, and the hydraulic conductivity reduction induced by the calibrated clogging factor λ corresponded well to the measured hydraulic conductivity reduction in the laboratory experiment. The numerical simulation of the aquifer-well experiment was similarly improved by the inclusion of the scaling factor reducing the soil hydraulic conductivity as the increase of the water level in the well due to physical clogging was accurately reproduced when considering the soil 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-variable scaling factor is feasible. The simulations indicate that the time-variable scaling factor, especially in combination with the reservoir boundary condition, can improve simulations of 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 this may not always be the case, the literature indicates that the exponential function is frequently used to approximate clogging (Pérez-Paricio, 2001). One also must consider that the numerical model that uses the scaling factor to reduce the hydraulic 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 (K s approaches 0), which usually leads to the secession of infiltration and the application of maintenance measures such as scraping of the infiltration surface in the infiltration basins or backflushing of injection wells.
Hypothetical simulations indicate that clogging has a relatively small effect on infiltration volumes from recharge wells, which can compensate the clogging effect by adjusting (increasing) the water level in the well, and thus the well infiltration area. Clogging can influence the HIE in the 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., recharge wells), and this should therefore be considered during the design 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, the presented 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 with regard to clogging. The scaling factor reducing the hydraulic conductivity 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 the interaction among infiltration structures, various clogging factors to account for different clogging rates in various areas of 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 hydraulic conductivity reduction (Equation 4) can be expanded to include an asymptotic term, so that K s 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 to changes in hydraulic conductivity over time but also other relevant soil hydraulic parameters such as soil porosity or the van Genuchten soil retention parameters. This could be achieved 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 function of the soil type, but also the quality of infiltrating water. For that, the application of more complex reactive transport models such as HP1 (Jacques, Šimůnek, Mallants, & van Genuchten, 2018), MIN3P (Mayer, Frind, & Blowes, 2002), or PHREEQC (Parkhust & Appelo, 1999) may be necessary.

ACKNOWLEDGMENTS
This study was supported by the German Federal Ministry of 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, water content, the water level in the infiltration basin, and tracer breakthrough curves). Also, the authors would like to thank José Bonilla Valverde (Instituto Costarricense de Acueductos y 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 the laboratory aquifer-well experiment.