Effects of Decaying Hydraulic Conductivity on the Groundwater Flow Processes in a Managed Aquifer Recharge Area in an Alluvial Fan

: Groundwater artiﬁcial recharge and medium characteristics represent the major factors in controlling the groundwater ﬂow processes in managed aquifer recharge areas. According to the depositional features of alluvial fans, an analogous homogeneous phreatic sand tank aquifer and the corresponding inhomogeneous scale numerical models were established to investigate the groundwater ﬂow under the combined inﬂuence of artiﬁcial recharge (human activities) and decaying hydraulic conductivity (medium characteristics). In this study, groundwater ﬂow through a managed aquifer recharge area in an alluvial fan was analyzed under the conditions of decaying hydraulic conductivity (K) with depth or length from apex to apron. The results showed that groundwater ﬂow processes induced by artiﬁcial recharge were signiﬁcantly controlled by the increasing decay exponents of K. The decaying K with depth or length in alluvial fan areas expanded the degree of inﬂuence of artiﬁcial recharge on groundwater ﬂow. With the increase of decay exponents, the ﬂow directions gradually changed from a horizontal to vertical direction. Groundwater age and spatial variability could also be increased by the increasing decay exponents. The residence time distributions (RTDs) of ambient groundwater and artiﬁcially recharged water exhibited logarithmic, exponential, and power law behavior. Penetration depth and travel times of ambient groundwater ﬂow could be affected by artiﬁcial recharge and decay exponents. Furthermore, with the increase of decay exponents, the thickness of the artiﬁcially recharged water lens and travel times of artiﬁcially recharged water were increased. These ﬁndings have important implications for the performance of managed aquifer recharge in alluvial fan areas as well as the importance of considering the gradual decrease of K with depth and length.


Introduction
Alluvial fans are becoming an important area of interest for groundwater artificial recharge [1][2][3][4][5] given their substantial thickness, high porosity, huge storage capacity [6], and abundant groundwater resources such as those in Northwest China [7] and Hueco Bolson in the United States and Mexico [8]. In most of the alluvial fans, excessive groundwater withdrawal to meet increasing demand for water resources because of population increase and economic growth has resulted in water table decline, land subsidence, and saltwater intrusion [9][10][11]. Groundwater artificial recharge is an important method of

Laboratory Experiments
An acrylic glass tank measuring 2.05 (length) × 0.5 (height) × 0.16 m (width) was carried out to simulate a two-dimensional cross-section of a homogeneous phreatic aquifer and used for laboratory experiments (Figure 1a). The general setup of the laboratory experiments was adopted from Atlabachew et al. [45] and Wu et al. [15]. However, the sand tank was only used for the homogeneous scenario here. Computer-controlled pumps were connected to three injection wells (IW 1-3). In total, 19 high-sensitivity deflated pressure transducers (PTs) (HM20-1-A1-F2-W2, by German HELM) were buried directly into the back face of the sand tank to measure hydraulic pressure. The saturated hydraulic conductivity of the sand was 40 m/d as measured using the constant-head method, and the bulk porosity was found to be 0.3 using the oven-drying method (Chinese National Standard). Detailed descriptions of the filter sands and the methods employed for the determination of their particle size can be found in Wu et al. [15].

Laboratory Experiments
An acrylic glass tank measuring 2.05 (length) × 0.5 (height) × 0.16 m (width) was carried out to simulate a two-dimensional cross-section of a homogeneous phreatic aquifer and used for laboratory experiments (Figure 1a). The general setup of the laboratory experiments was adopted from Atlabachew et al. [45] and Wu et al. [15]. However, the sand tank was only used for the homogeneous scenario here. Computer-controlled pumps were connected to three injection wells (IW 1-3). In total, 19 high-sensitivity deflated pressure transducers (PTs) (HM20-1-A1-F2-W2, by German HELM) were buried directly into the back face of the sand tank to measure hydraulic pressure. The saturated hydraulic conductivity of the sand was 40 m/d as measured using the constant-head method, and the bulk porosity was found to be 0.3 using the oven-drying method (Chinese National Standard). Detailed descriptions of the filter sands and the methods employed for the determination of their particle size can be found in Wu et al. [15].
It is important to point out that the sand tank experiments acted as the foundation of the following numerical model. That is, the sand tank experiments had neutral conditions with non-decaying hydraulic conductivity. The sand tank experiments could also show phenomena that actually occurred in the artificially recharged aquifer, for example, the flow field above the water table.

Numerical Simulation
Exploration into the roles of hydrodynamics (surface and subsurface) based on an actual research area requires detailed geomorphological and hydrogeological data, which are rarely available [46], so a synthetic model is always needed for solving such problems. In this study, a theoretical groundwater flow model generalized from an alluvial fan was constructed to explore the role of the artificial recharge and decaying K in controlling the It is important to point out that the sand tank experiments acted as the foundation of the following numerical model. That is, the sand tank experiments had neutral conditions with non-decaying hydraulic conductivity. The sand tank experiments could also show phenomena that actually occurred in the artificially recharged aquifer, for example, the flow field above the water table.

Numerical Simulation
Exploration into the roles of hydrodynamics (surface and subsurface) based on an actual research area requires detailed geomorphological and hydrogeological data, which Water 2021, 13, 1649 4 of 20 are rarely available [46], so a synthetic model is always needed for solving such problems. In this study, a theoretical groundwater flow model generalized from an alluvial fan was constructed to explore the role of the artificial recharge and decaying K in controlling the groundwater flow by using the finite element model FEFLOW, which is a variably saturated groundwater flow modeling program. The homogeneous numerical model was constructed according to the sand tank and adjected slightly. Then, the numerical models with decaying K were expanded based on the homogeneous one. The finite element mesh and boundary conditions for the numerical solution are shown in Figure 1b.
The basic equation describing the groundwater flow simulated in the FEFLOW code is shown below.
where s is the saturation of fluid in the void space " [1], S 0 is the specific storage coefficient  [1], Ø is the buoyancy coefficient [1], and e is the gravitational unit vector [1]. The nonlinear depth-dependent K (z) was described by using the exponential decay model [35,36]. Assuming a locally isotropic condition, the depth-dependent K (z) is given by: where K 0 is the K at the ground surface. The K 0 adopted in the numerical model was the same as that measured in the sand tank. A 1 is the decay exponent that dictates the rate of decrease in K with depth, z s (x) is a function of the ground surface elevation, and z is the elevation of the aquifer bottom.
There is a lack of relevant research on the variation law of K with the length from apex to apron zone, which approximately follows an exponential law in the Tailan river basin ( Figure 2). The main purpose of this study was to investigate how the decaying hydraulic conductivity affects the groundwater flow; a further precise study on the variation law was not conducted here, and the length-dependent K(x) was assumed to be: where A 2 is the decay exponent that dictates the rate of decrease in K with length, and x u (z) is a function of the vertical boundary on the upper reaches side. groundwater flow by using the finite element model FEFLOW, which is a variably saturated groundwater flow modeling program. The homogeneous numerical model was constructed according to the sand tank and adjected slightly. Then, the numerical models with decaying K were expanded based on the homogeneous one. The finite element mesh and boundary conditions for the numerical solution are shown in Figure 1b. The basic equation describing the groundwater flow simulated in the FEFLOW code is shown below.
where s is the saturation of fluid in the void space ε [1], S 0 is the specific storage coeffi- , ε is the porosity (void space) [1], ∇ is the Nabla (vector) operator [L −1 ], q is the Darcy velocity of fluid [L T −1 ], Q is the bulk source/sink term of flow [L 3 T −1 ], QEOB is the correction sink/source term of the extended Oberbeck-Boussinesq approximation [T −1 ], k r is the relative permeability [1], K is the tensor of hydraulic conductivity [LT −1 ], f μ is the viscosity relation function [1], is the buoyancy coefficient [1], and e is the gravitational unit vector [1]. The nonlinear depth-dependent K (z) was described by using the exponential decay model [35,36]. Assuming a locally isotropic condition, the depth-dependent K (z) is given by: where K0 is the K at the ground surface. The K0 adopted in the numerical model was the same as that measured in the sand tank. A1 is the decay exponent that dictates the rate of decrease in K with depth, zs(x) is a function of the ground surface elevation, and z is the elevation of the aquifer bottom.
There is a lack of relevant research on the variation law of K with the length from apex to apron zone, which approximately follows an exponential law in the Tailan river basin ( Figure 2). The main purpose of this study was to investigate how the decaying hydraulic conductivity affects the groundwater flow; a further precise study on the variation law was not conducted here, and the length-dependent K(x) was assumed to be: where A2 is the decay exponent that dictates the rate of decrease in K with length, and xu(z) is a function of the vertical boundary on the upper reaches side.

Scenario Definition
Three scenarios (scenarios A-C) were conducted to mimic the groundwater flow induced by artificial recharge and decaying K. Scenario A was conducted in the sand tank aquifer and corresponding numerical model and represents the condition of artificial recharge in a homogeneous phreatic aquifer. The water levels were controlled at 43 and 25 cm above the inner base of the tank at the entrance (left) and exit (right) water chambers, respectively. The flow fields were measured in the steady-state condition.
Scenario B was conducted using a numerical model, based on the homogeneous one corresponding to sand tank experiment, with artificial recharge at different locations (infiltration basin and injection wells at different depths of the lab-scale aquifer) and decaying K with depth. Various decay exponents with depth (A 1 ) have been reported in the literature, mostly ranging from 0.003 to 0.5 [20,33,35,36]. Here, we set the decay exponents as 0.01, 0.05, and 0.1 m −1 . The corresponding K fields and their relationships to depth are shown in Figures 3 and 4, respectively. The artificial recharge rate through the infiltration basin was set to 5 and 0.2 m/d in the injection wells. The water level in the right and left water chambers was the same as that in scenario A.
Three scenarios (scenarios A-C) were conducted to mimic the groundwater flow induced by artificial recharge and decaying K. Scenario A was conducted in the sand tank aquifer and corresponding numerical model and represents the condition of artificial recharge in a homogeneous phreatic aquifer. The water levels were controlled at 43 and 25 cm above the inner base of the tank at the entrance (left) and exit (right) water chambers, respectively. The flow fields were measured in the steady-state condition.
Scenario B was conducted using a numerical model, based on the homogeneous one corresponding to sand tank experiment, with artificial recharge at different locations (infiltration basin and injection wells at different depths of the lab-scale aquifer) and decaying K with depth. Various decay exponents with depth (A1) have been reported in the literature, mostly ranging from 0.003 to 0.5 [20,33,35,36]. Here, we set the decay exponents as 0.01, 0.05, and 0.1 m −1 . The corresponding K fields and their relationships to depth are shown in Figures 3 and 4, respectively. The artificial recharge rate through the infiltration basin was set to 5 and 0.2 m/d in the injection wells. The water level in the right and left water chambers was the same as that in scenario A.
Scenario C was also conducted with a numerical model, with artificial recharge at different locations and decaying K with length from apex to apron zone. Various decay exponents along the length from the apex to apron are not reported. For the main purpose of investigating the response of artificially recharged groundwater flow to the decay of K along the length direction, the A2 was assumed as 0.005, 0.01, and 0.02 m −1 here. The corresponding K fields and their relationships with length are shown in Figures 5 and 6, respectively. The artificial recharge rate and water level in the water chambers were the same as those in scenario B.  Scenario C was also conducted with a numerical model, with artificial recharge at different locations and decaying K with length from apex to apron zone. Various decay exponents along the length from the apex to apron are not reported. For the main purpose of investigating the response of artificially recharged groundwater flow to the decay of K along the length direction, the A 2 was assumed as 0.005, 0.01, and 0.02 m −1 here. The corresponding K fields and their relationships with length are shown in Figures 5 and 6,     The flow fields became increasingly complex when driven by the same artificial recharge rate with the increase of decay exponents. That is, the flow directions changed from a horizontal to vertical direction, especially near the artificially recharged areas. The antidromic flow area, where the groundwater flows in the opposite direction of the ambient groundwater flow, became larger with the increasing decay exponents, which can be more obviously seen in Figure 8f.

Groundwater Flow Patterns
Stagnation points, where groundwater flow converges from and parts in opposite directions, appeared at different locations when driven by the same artificial recharge rate and when the decay exponents were set to different values. Locations of stagnation points were not obvious when the decay exponents were small and had a small injection rate. Left side figures in Figures 7 and 8 show that the locations of the stagnation points produced by the infiltration basin (SP1) moved to the upstream direction with the increase of decay exponents. Those produced by injection well IW3 (SP2) moved to the upstream direction and the deeper layer of the aquifer in the same conditions.  The flow fields became increasingly complex when driven by the same artificial recharge rate with the increase of decay exponents. That is, the flow directions changed from a horizontal to vertical direction, especially near the artificially recharged areas. The antidromic flow area, where the groundwater flows in the opposite direction of the ambient groundwater flow, became larger with the increasing decay exponents, which can be more obviously seen in Figure 8f.

Groundwater Flow Patterns
Stagnation points, where groundwater flow converges from and parts in opposite directions, appeared at different locations when driven by the same artificial recharge rate and when the decay exponents were set to different values. Locations of stagnation points were not obvious when the decay exponents were small and had a small injection rate. Left side figures in Figures 7 and 8 show that the locations of the stagnation points produced by the infiltration basin (SP1) moved to the upstream direction with the increase of decay exponents. Those produced by injection well IW3 (SP2) moved to the upstream direction and the deeper layer of the aquifer in the same conditions.

Groundwater Age
Figures 9 and 10 show the groundwater age fields produced by the infiltration basin and injection well IW3 with variational decay exponents of K in different directions. It can be seen that the smaller artificial recharge rate (the artificial recharge rate in infiltration basin is 5 and 0.2 m/d in injection well IW3) has a limited effect on the groundwater age distribution when the decay exponents were small. The groundwater age always increased along the flow direction from upstream to downstream area with the decaying K. For example, groundwater age ranged from 0 to 2547.3 s (Figure 9a,b) when the depthdecay exponents were set to 0.01 m −1 . As shown in Figure 10a,b, age mainly ranged from 0 to 2711 s when the length-decay exponents were set to 0.005 m −1 . Furthermore, the increasing decay exponents can also increase the groundwater age and expand the degree of influence of the artificial recharge on groundwater age distribution. Taking left side figures in Figure 9 as an example, with the increase of depth-decay exponents, the spatial variability of groundwater age becomes increasingly obvious.

Groundwater Age
Figures 9 and 10 show the groundwater age fields produced by the infiltration basin and injection well IW3 with variational decay exponents of K in different directions. It can be seen that the smaller artificial recharge rate (the artificial recharge rate in infiltration basin is 5 and 0.2 m/d in injection well IW3) has a limited effect on the groundwater age distribution when the decay exponents were small. The groundwater age always increased along the flow direction from upstream to downstream area with the decaying K. For example, groundwater age ranged from 0 to 2547.3 s (Figure 9a,b) when the depthdecay exponents were set to 0.01 m −1 . As shown in Figure 10a,b, age mainly ranged from 0 to 2711 s when the length-decay exponents were set to 0.005 m −1 . Furthermore, the increasing decay exponents can also increase the groundwater age and expand the degree of influence of the artificial recharge on groundwater age distribution. Taking left side figures in Figure 9 as an example, with the increase of depth-decay exponents, the spatial variability of groundwater age becomes increasingly obvious. Figure 11 shows the groundwater age along the bottom of the aquifer induced by injection at different locations with variational depth-decay exponents, indicating little distinction for the age along the aquifer bottom when the decay exponent was small. On the contrary, the variation of groundwater age showed an obvious distinction when the decay exponent was larger. The results shown by Figure 11 could also demonstrate that the decaying K could expand the degree of influence of artificial recharge on the distributions of groundwater age.   Figure 11 shows the groundwater age along the bottom of the aquifer induced by injection at different locations with variational depth-decay exponents, indicating little distinction for the age along the aquifer bottom when the decay exponent was small. On the contrary, the variation of groundwater age showed an obvious distinction when the decay exponent was larger. The results shown by Figure 11 could also demonstrate that the decaying K could expand the degree of influence of artificial recharge on the distributions of groundwater age.

Residence Time Distributions
RTDs of the sand tank phreatic aquifer induced by artificial recharge and decaying hydraulic conductivity were considered. The results showed that not only the ambient groundwater but also the artificially recharged water exhibit different RTD laws with the decaying K in different directions and different decay exponents. Taking the RTDs driven Artificial recharged through injection well IW1 when decay exponent is 0.01 Artificial recharged through injection well IW3 when decay exponent is 0.01 Artificial recharged through injection well IW2 when decay exponent is 0.01

Residence Time Distributions
RTDs of the sand tank phreatic aquifer induced by artificial recharge and decaying hydraulic conductivity were considered. The results showed that not only the ambient groundwater but also the artificially recharged water exhibit different RTD laws with the decaying K in different directions and different decay exponents. Taking the RTDs driven Artificial recharged through injection well IW1 when decay exponent is 0.01 Artificial recharged through injection well IW3 when decay exponent is 0.01 Artificial recharged through injection well IW2 when decay exponent is 0.01 Figure 11. Groundwater age at the bottom of the aquifer driven by artificial recharge through injection wells at different locations and different decay exponents.

Residence Time Distributions
RTDs of the sand tank phreatic aquifer induced by artificial recharge and decaying hydraulic conductivity were considered. The results showed that not only the ambient groundwater but also the artificially recharged water exhibit different RTD laws with the decaying K in different directions and different decay exponents. Taking the RTDs driven by injection at well IW1 as an example, the frequency distributions (FDs) to RTDs of artificially recharged water driven by variational decay exponents are shown in Figure 12, and those of ambient groundwater are shown in Figure 13. The results of depth-decay condition (Figure 12), with a smaller decay exponent (when the depth-decay exponent was set to 0.01 m −1 ), showed the fitted trend line for FDs to RTDs toward a quadratic polynomial law and the frequency of different residence time showed less distinction. With the increase of depth-decay exponents, when the depth-decay exponent was set to 0.05 m −1 , the fitted trend line for FDs to RTDs was divided into two steps: the early step toward a logarithmic behavior and the later step toward a power behavior. When the depth-decay exponent was set to a sufficiently large value (A 1 = 0.1 m −1 ), the two-step fitted trend line merged into one and toward a logarithmic behavior. Furthermore, similar results could be obtained with the K decay in the length direction. With the increase of the length-decay exponents, the fitted lines for FDs to RTDs exhibited a quadratic polynomial behavior and then was divided into two steps. Decay of K in different directions could affect the RTD behavior obviously.
by injection at well IW1 as an example, the frequency distributions (FDs) to RTDs of artificially recharged water driven by variational decay exponents are shown in Figure 12, and those of ambient groundwater are shown in Figure 13. The results of depth-decay condition (Figure 12), with a smaller decay exponent (when the depth-decay exponent was set to 0.01 m −1 ), showed the fitted trend line for FDs to RTDs toward a quadratic polynomial law and the frequency of different residence time showed less distinction. With the increase of depth-decay exponents, when the depth-decay exponent was set to 0.05 m −1 , the fitted trend line for FDs to RTDs was divided into two steps: the early step toward a logarithmic behavior and the later step toward a power behavior. When the depth-decay exponent was set to a sufficiently large value (A1 = 0.1 m −1 ), the two-step fitted trend line merged into one and toward a logarithmic behavior. Furthermore, similar results could be obtained with the K decay in the length direction. With the increase of the length-decay exponents, the fitted lines for FDs to RTDs exhibited a quadratic polynomial behavior and then was divided into two steps. Decay of K in different directions could affect the RTD behavior obviously.
FDs to RTDs of ambient groundwater are shown in Figure 13, with different decay direction and decay exponents of K when artificially recharged at injection well IW1. In contrast to the RTD behaviors of artificially recharged water, when the decay exponents were set to smaller values, the fitted trend lines for FDs to RTDs always showed two steps. However, the two-step lines were merged into one when the decay exponents were set to larger values. Most of the artificially recharged water RTD fitted trend lines showed logarithmic behavior and power behavior.

Ambient Groundwater Flow Paths
Ambient groundwater flow paths induced by artificial recharge and decaying K were considered here. Several particles were released at tracing points in the numerical model on the left side to trace the ambient groundwater flow. Taking those produced by the infiltration basin and injection well IW3 in the aquifer with decaying K as examples, the ambient groundwater flow paths are shown in Figures 14 and 15 with different decay directions. The left panel of Figure 14 shows the flow paths driven by the infiltration basin and variational decay exponents in the depth direction, indicating that the increasing decay exponents could increase the penetration depth and travel times of ambient groundwater. The penetration depth here means the depth of certain groundwater flow paths at certain sections. The right panel of Figure 14 shows the ambient groundwater flow paths induced by injection well IW3 and decaying K in the depth direction. With the increase of depth-decay exponents, the ambient groundwater flow paths were gradually fractionated by artificial recharge, meaning that the penetration depth of shallow layer's flow decreased with the increasing decay exponent, in contrast to that of the deep layer. The penetration depth of the ambient groundwater flow paths was decreased by the decaying K when the decay exponent was sufficiently large. Similar to the results shown by right panel of Figure 14, the ambient groundwater flow paths could also be found gradually fractionated by the screen of injection wells with the increase of the length-decay expo- Figure 13. Frequency distribution of RTDs of ambient groundwater with different decay direction and decay exponents of hydraulic conductivity. "D-0.01" means the depth-decay exponent was set to 0.01, and "L-0.005" means the length-decay exponent was set to 0.005. The red lines are fitted trend lines.
FDs to RTDs of ambient groundwater are shown in Figure 13, with different decay direction and decay exponents of K when artificially recharged at injection well IW1. In contrast to the RTD behaviors of artificially recharged water, when the decay exponents were set to smaller values, the fitted trend lines for FDs to RTDs always showed two steps. However, the two-step lines were merged into one when the decay exponents were set to larger values. Most of the artificially recharged water RTD fitted trend lines showed logarithmic behavior and power behavior.

Ambient Groundwater Flow Paths
Ambient groundwater flow paths induced by artificial recharge and decaying K were considered here. Several particles were released at tracing points in the numerical model on the left side to trace the ambient groundwater flow. Taking those produced by the infiltration basin and injection well IW3 in the aquifer with decaying K as examples, the ambient groundwater flow paths are shown in Figures 14 and 15 with different decay directions. The left panel of Figure 14 shows the flow paths driven by the infiltration basin and variational decay exponents in the depth direction, indicating that the increasing decay exponents could increase the penetration depth and travel times of ambient groundwater.
The penetration depth here means the depth of certain groundwater flow paths at certain sections. The right panel of Figure 14 shows the ambient groundwater flow paths induced by injection well IW3 and decaying K in the depth direction. With the increase of depthdecay exponents, the ambient groundwater flow paths were gradually fractionated by artificial recharge, meaning that the penetration depth of shallow layer's flow decreased with the increasing decay exponent, in contrast to that of the deep layer. The penetration depth of the ambient groundwater flow paths was decreased by the decaying K when the decay exponent was sufficiently large. Similar to the results shown by right panel of Figure 14, the ambient groundwater flow paths could also be found gradually fractionated by the screen of injection wells with the increase of the length-decay exponents when artificially recharged through injection well IW2. In contrast to those shown by Figure 14 and the left panel of Figure 15, the right panel of Figure 15 shows that the penetration depth of ambient groundwater flow paths decreased with the increase of the length-decay exponents when artificially recharged by injection well IW1. Furthermore, with the change of ambient groundwater flow paths, travel times could also be increased by the increasing decay exponents.
Water 2021, 13, x FOR PEER REVIEW 13 of 21 by Figure 14 and the left panel of Figure 15, the right panel of Figure 15 shows that the penetration depth of ambient groundwater flow paths decreased with the increase of the length-decay exponents when artificially recharged by injection well IW1. Furthermore, with the change of ambient groundwater flow paths, travel times could also be increased by the increasing decay exponents.

Artificially Recharged Water Lens
An artificially recharged water lens was formed in artificial recharge areas due to artificially induced recharge to the lens. The interface between artificially recharged water and ambient groundwater is a part of the boundary of such lenses. Artificially recharged water lenses produced by the infiltration basin and injection wells in the aquifer with decaying K are shown in Figures 16 and 17, in the form of a tracing-line cluster from infiltration basins or injection wells. Left panels of Figures 16 and 17 show the artificially recharged water lenses produced by the infiltration basins with different decay directions and decay exponents, indicating that the increasing decay exponents increased the thickness of the artificially recharged water lens and the travel times of artificially recharged water. The penetration depth of the interface of artificially recharged-ambient groundwater gradually increased with the increase of decay exponents and moved to the deep layers of the aquifer. Especially when the length-decay exponent was set to a sufficiently large value (A 2 = 0.002 m −1 , Figure 17e), the aquifer was almost filled with artificially recharged water, and the left recharge boundary became a discharge boundary. The right panels of Figure 16 show that the thickness of the artificially recharged water lens increased with the increasing depth-decay exponent, with the up-interface moving to the shallow layer and the down-interface to deep layer. Similar results can also be seen in right panel of Figure 17 when the decay of K was set in the length direction and artificially recharged through injection well IW1. With the increase of length-decay exponents, the interface of the artificially recharged-ambient groundwater interface gradually moved to the shallow layer, and the thickness of the artificially recharged water lens increased. Furthermore, with the increase of decay exponents, increased travel times could be seen in every scenario.

Artificially Recharged Water Lens
An artificially recharged water lens was formed in artificial recharge areas due to artificially induced recharge to the lens. The interface between artificially recharged water and ambient groundwater is a part of the boundary of such lenses. Artificially recharged water lenses produced by the infiltration basin and injection wells in the aquifer with decaying K are shown in Figures 16 and 17, in the form of a tracing-line cluster from infiltration basins or injection wells. Left panels of Figures 16 and 17 show the artificially recharged water lenses produced by the infiltration basins with different decay directions and decay exponents, indicating that the increasing decay exponents increased the thickness of the artificially recharged water lens and the travel times of artificially recharged water. The penetration depth of the interface of artificially recharged-ambient groundwater gradually increased with the increase of decay exponents and moved to the deep layers layer, and the thickness of the artificially recharged water lens increased. Furthermore, with the increase of decay exponents, increased travel times could be seen in every scenario.

Discussion
Alluvial plains are always deeply affected by human activities and richness of fresh groundwater [47]. Artificial recharge involves intensive activities to regulate water resources and represents a major factor in controlling the groundwater flow system on a regional scale [15]. At the same time, the spatial variation law of K in aquifers is funda-

Discussion
Alluvial plains are always deeply affected by human activities and richness of fresh groundwater [47]. Artificial recharge involves intensive activities to regulate water resources and represents a major factor in controlling the groundwater flow system on a regional scale [15]. At the same time, the spatial variation law of K in aquifers is fundamental to understanding the groundwater flow processes. Thus, based on the depositional features of the alluvial fan, lab-scale experiments and numerical simulations were performed in this study to investigate the groundwater flow features under the combined effects of artificial recharge and decaying K. Settings of artificial recharge modes were adopted from Wu et al. [15] and generalized from the Tailan river basin. Because surface water flowing through alluvial fans in arid and semiarid areas is often disconnected from groundwater, the infiltration basin was set to a prescribed flux (Neumann type) boundary [15,41,48].
The gradual decline in K with depth and length is a hallmark of piedmont alluvial fan hydrogeology [17,19], which influences features of subsurface flow such as the distribution of path lines, velocity, and residence time along these path lines. Depth-dependent K is a manifestation of regional-scale heterogeneity and was widely observed in the geological medium [20,[35][36][37]49]. Different relationships of K with depth have been fitted in different sites with different K range, for example, the power function [23], logistic function [23], and exponential function [33]. In our analysis, the aquifer was isotropic, and an exponentially decaying function was used to describe the decrease of K with depth for the similar K range to that of Ameli et al. [20]. There is a lack of data and results on the variation law of K with length from apex to apron area in alluvial fans, but the variation trend could be generalized based on the depositional features of Tailan alluvial fans. Thus, considering that the main purpose of this study was to investigate the groundwater flow induced by artificial recharge and the decaying K, here we assumed an exponentially decaying function to describe the gradual decrease of K with length from apex to apron area in the alluvial fan.
Previous research on groundwater flow under depth decay conditions was performed by Jiang et al. [35] and Ameli et al. [20]. Jiang's [35] results showed that the development of local versus regional flow systems is sensitive to the depth-decay exponent of K. With higher rates of decrease in K with depth, the penetration depth of local flow systems increases. That is, penetration depth of flow path from a certain recharge point decreased with the increase of depth-decay exponents. Similar to Jiang's [35] results, Ameli et al. [20] demonstrated that the increasing depth-decay exponents decreased the penetration depth of flow paths. In contrast with the previous work, numerical simulations conducted in this study showed that the penetration depth of flow paths from infiltration basin increased with the increase of depth-decay exponents. The opposite results were attributed to the differences of boundary conditions. In Jiang's [35] and Ameli's [20] studies, no flux boundaries were set to the left and right boundary in the numerical model, but they were constant head boundaries in this study. The main flow directions of ambient groundwater were horizontal when the constant head boundary was used in the model. With the increase of decay exponents, the horizontal flow velocities were lower and exerted less influence on the artificially recharged water flow, which was almost vertical. Thus, the resultant velocities were almost vertical, and the penetration depth was deeper.
The RTDs, the descriptions of residence times from recharge to discharge areas, determine the type and rate of many biological, biochemical, and geochemical processes occurring on the surface, near the surface, and in deep environments [34]. The most commonly used model used for RTDs is the exponential model [50], and some recent studies demonstrated that regional and multiscale flow systems exhibit a power law [34,50,51]. RTDs driven by human activities, such as artificial recharge and exploitation, exhibit different laws with the natural state. Most of them showed exponential law and logarithmic behavior [15]. According to Cardenas and Jiang's [34] results, the RTDs for both a simple basin with one flow system and a basin with nested local and regional systems gradually evolve to a power law RTD with exponentially depth-decreasing K and porosity. Here, the results from the numerical simulations demonstrated that the combined effect of artificial recharge and decaying K could lead to a different RTD behavior compared with the previous studies. That is, some of the sub-steps exhibit logarithmic law behavior or quadratic polynomial law behavior.
Notably, similar to the results presented by Wu et al. [15], the groundwater flow could pass across the water table, which is common in the simulated results and lab-scale experiments. This was attributed to the limitation of using the small-scale model. In the lab-scale model, the capillary zone thickness relative to the saturated zone thickness was extremely large in these simulations, which seems unreasonable when considering a field-scale case. The capillary zone, in agreement with Lehmann et al. [52], was 13 and 25 cm for packed sand in this study. Thus, the limitation did not affect the main purpose of this study (i.e., understanding how decaying hydraulic conductivity and artificial recharge affect the groundwater flow dynamics). As such, the process of groundwater flow in these experiments and simulations was plausible.
The differences in properties between ambient groundwater and artificially recharged water were not taken into account in this study, although, to the best of our knowledge, the artificial recharge accelerated the interaction of surface water and groundwater. Especially when artificially recharged through the injection wells, the physical and chemical properties of artificially recharged water are different from those of ambient groundwater due to the absence of vadose filtration of the aquifer and the corresponding physical and chemical reactions. For example, the temperature difference between surface and groundwater has been widely observed in many sites worldwide and used to trace the water flow [44,53,54]. Artificial recharge also increases the overall dissolved oxygen (DO) concentration in recharged areas, as it is commonly believed that DO concentration of groundwater is much lower than that of surface water [55]. Thus, accurate simulation of the artificially recharged water lenses and the artificially recharged-ambient groundwater interface was essential to developing adequate groundwater management and protection schemes. According to the results of this paper, decaying K increases the thickness of the artificially recharged water lens and changes the locations of the artificially rechargedambient groundwater interface compared with the homogeneous conditions. Due to economic development, population growth, urbanization, and irrigation-area increase, human activities are being viewed as the main driver of hydrogeological processes [49,56,57]. However, when coupled with the decaying K, which has been widely observed in geological medium both in the depth direction and length direction in alluvial fan areas, the groundwater flow processes caused by human activities increase. As shown by the results of this study, ambient groundwater and artificially recharged water flow paths, groundwater age distributions, locations of stagnation points, RTDs, and the locations of the artificially recharged-ambient groundwater interface were totally different from each other with the variation of decay direction and exponents. Thus, the decay of K should not be neglected when analyzing hydrological problems related to artificial regulated groundwater flow. In other words, based on the results of this study, neglecting the decaying K may cause an underestimation of the effect of human activities on groundwater systems.

Conclusions
Human activities and medium characteristics represent the major factors in controlling the groundwater flow system on a regional scale. Zonal heterogeneity and layered heterogeneity of K were applied in numerical models to investigate groundwater flow affected by human activities. The systematic decrease in K in a certain direction, especially with length, was seldom accounted for in studies on groundwater flow. Based on the lab-scale sand tank and depositional features of the alluvial fan, lab experiments and numerical simulations were performed in this study to investigate the variation of groundwater flow under the combined influences of artificial recharge and decaying K with the depth and length direc-tion. The lab-scale sand tank experiments and numerical simulations demonstrated that the groundwater flow features induced by artificial recharge were significantly controlled by the decaying K. The decay of K with depth or length in alluvial fan areas expanded the degree of influence of artificial recharge on ambient groundwater flow processes.
Our findings also have some practical implications for the performance of managed aquifer recharge. Based on these results, recommendations for the performance of MAR in alluvial fan aquifers could be made as follows: (a) A relatively low artificial recharge rate should be adopted in alluvial fan aquifers, so as not to affect the natural groundwater cycle significantly. (b) The screening of recovery wells should be much longer than that for homogeneous aquifers in order to capture more recharged water, as the thickness of the artificially recharged water lens is much larger in an alluvial fan aquifer with decaying hydraulic conductivity. (c) More attention should be paid to the water quality, as the chemical loading rates were strongly controlled by the groundwater flow paths, which change much more obviously in an alluvial fan aquifer. (d) A relatively low pumping rate should be adopted in managed alluvial fan aquifers, as the increased groundwater age due to decaying hydraulic conductivity results in a low cycle rate.