Variations of Groundwater Dynamics in Alluvial Aquifers with Reclaimed Water Restoring the Overlying River, Beijing, China

Some of the rivers in northern China are dried, and reclaimed water (RW) is used to restore these degraded river ecosystems, during which the RW could recharge the aquifer by river bank infiltration. From 2007 to 2018, 2.78 × 108 m3 of RW has been replenished to the dried Chaobai River (Shunyi reach), Beijing, China, which is located on the edge of one depression cone in groundwater caused by groundwater over-pumping. The groundwater hydrodynamic variations and the flow path of the RW were identified by eight-year hydrological, hydrochemical, and stable isotopic data, together with multivariate statistical analysis. The RW infiltration drastically impacts the groundwater dynamics with a spatiotemporal variation. The 30-m depth groundwater levels at Perennial intake reach increased quickly around 3 m after 2007, which indicated that they were dominated by RW infiltration. Other 30-m depth groundwater levels were controlled by precipitation recharge from 2007 to 2011, showing significant seasonal variations. In 2012, with more RW transferred to the river, the hydrodynamic impact of the RW on 30-m depth aquifer expanded downstream. However, the 50-m and 80-m depth groundwater levels showed decreasing trend with seasonal variations, due to groundwater pumping. The 30-m depth aquifer was mainly recharged by RW, being evidenced by the enriched δ2H and δ18O. The depleted δ2H and δ18O of the 50-m and 80-m depth groundwater indicated that they were dominated by regional groundwater with meteoric origin. The heterogenous properties of the multi-layer alluvial aquifer offer the preferential flow path for RW transport in the aquifers. The proportion of the RW in the aquifers decreases with depth that was calculated by the chloride conservative mixing model. The increased lateral hydraulic gradient (0.43%) contributes to the RW transport in the 30-m depth aquifer. RW usage changed 30-m depth groundwater type from Ca·Mg-HCO3 to Na·Ca·Mg-HCO3·Cl. RW preferentially recharged the 50-m and 80-m depth aquifers by vertical leakage.


Introduction
In the water shortage areas, climate change and anthropogenic activities have led to a declined water table, decreased runoff, and drying of rivers with a degraded ecosystem [1]. RW generally refers to the tertiary treated wastewater that is produced by Reclaimed Treatment Water Plant (RTWP) [2]. RW, an important alternative water resource with advantages of stability and controllability, is increasingly essential in addressing the watersupply limitations [3]. RW has been widely used to replenish urban landscape pounds for aesthetic and recreational purposes and to feed nature rivers for restoring the degraded This is the first study to use long-term (before and after RW recharge, eight years), multi-layer (30-m, 50-m, and 80-m depth aquifers) hydrochemical and hydrological data, as well as stable isotopes to discuss the impact of the RW utilization on the groundwater dynamic system. The objectives of this study are (i) to characterize the hydrodynamic variation in the aquifers, (ii) to illustrate the transport of RW in multi-layer aquifers, and (iii) to estimate the percentage of RW in the groundwater system. It can improve the understanding of the subsurface hydrological processes with RW utilization in the alluvial plain, Beijing, and then assist the restoration of a seriously degraded river ecosystem, and the sustainable groundwater management in other similar areas globally.

Location and Hydrogeological Setting
The Shunyi reach of the Chaobai River is located in the northeast of Beijing, China ( Figure 1). It is characterized by a continental monsoon climate, with four distinct seasons. Based on a 59-year record in a local meteorological station (Figure 1), the average annual air temperature is 11.2 °C, with a mean relative humidity of 58.6%, the pan evaporation is 1690 mm, and the mean annual precipitation is 643 mm. A total of 70% precipitation falls in July-September [37].  The study area belongs to the piedmont alluvial area of the Chaobai River. The elevation of the mountainous area in upstream of the Chaobai River plain is generally from 400 m to 1500 m. Archean metamorphite, the middle-upper Proterozoic dolomite, Cretaceous rhyolite, Cambrian, and Ordovician limestone outcrop in the mountainous area [39]. The thickness of the unconsolidated Quaternary deposits ranges between 300 and 400 m in the study area [39]. The bedrocks underlying the Quaternary deposits belong to the middle-upper Proterozoic, the Cambrian, and the Ordovician group [39]. There are three aquifers within 80 m depth of the Quaternary stratum. From north to south, the lithology of the aquifer changes from gravel to fine sand, and to silty clay with poor permeability (Figure 1). The first 30-m depth aquifer is mainly composed of gravel and fine sand, with a buried depth of 0-30 m. The second 50-m depth aquifer is interlayered by gravel and sand within 30-50 m depth. The third 80-m depth aquifer with a depth of 50-80 m is mainly dominated by fine sand (Figure 2). The hydraulic conductivity (K) of 50-m and 80-m depth aquifers are 1.2 × 10 −3 -4.6 × 10 −3 m/s and 6.9 × 10 −4 -1.2×10 −3 m/s, respectively [40]. Between the 30-m and 50-m depth aquifers, there is a discontinuous clay layer with a thickness of more than 10 m, while there is a continuous clay layer that is distributed between the 50-m and 80-m depth aquifers. The extremely heterogenous aquifer medium could lead to complex groundwater flow systems, which may also control the RW transport in the aquifers.  Figure 1a, showing the water sampling locations and groundwater monitoring cross-sections. The time of reclaimed water transfer in different reaches are shown in brackets. AM: Anti-seepage reach monitoring cross-section; PM: Perennial intake reach monitoring cross-section; IM: Intermittent intake reach monitoring cross-section. Groundwater level contours and depression cone in groundwater (enclosed by 0 m groundwater level) in Nov. 2011 are referenced from Li et al. [38].
The study area belongs to the piedmont alluvial area of the Chaobai River. The elevation of the mountainous area in upstream of the Chaobai River plain is generally from 400 m to 1500 m. Archean metamorphite, the middle-upper Proterozoic dolomite, Cretaceous rhyolite, Cambrian, and Ordovician limestone outcrop in the mountainous area [39]. The thickness of the unconsolidated Quaternary deposits ranges between 300 and 400 m in the study area [39]. The bedrocks underlying the Quaternary deposits belong to the middle-upper Proterozoic, the Cambrian, and the Ordovician group [39]. There are three aquifers within 80 m depth of the Quaternary stratum. From north to south, the lithology of the aquifer changes from gravel to fine sand, and to silty clay with poor permeability ( Figure 1). The first 30-m depth aquifer is mainly composed of gravel and fine sand, with a buried depth of 0-30 m. The second 50-m depth aquifer is interlayered by gravel and sand within 30-50 m depth. The third 80-m depth aquifer with a depth of 50-80 m is mainly dominated by fine sand (Figure 2). The hydraulic conductivity (K) of 50-m and 80-m depth aquifers are 1.2 × 10 −3 -4.6 × 10 −3 m/s and 6.9 × 10 −4 -1.2×10 −3 m/s, respectively [40]. Between the 30-m and 50-m depth aquifers, there is a discontinuous clay layer with a thickness of more than 10 m, while there is a continuous clay layer that is distributed between the 50-m and 80-m depth aquifers. The extremely heterogenous aquifer medium could lead to complex groundwater flow systems, which may also control the RW transport in the aquifers.  Figure 1 along the Chaobai River, showing the screened intervals of the monitoring wells. The groundwater levels are shown for monitoring results in September 2007 [15]. The Chaobai River was dried before the reclaimed water transfer, so the surface water was not shown.
The confined aquifers of the Chaobai River with a depth of less than 120 m are the most important drinking water source for Beijing [39,41]. There is a well field 4 km off the north of the study area ( Figure 1). The groundwater was also pumped for agricultural activities [39]. In 2000, the pumping yield along the Chaobai River was about 35% of the city's total amount of groundwater exploitation [39]. Beijing's primary depression cone of  Figure 1 along the Chaobai River, showing the screened intervals of the monitoring wells. The groundwater levels are shown for monitoring results in September 2007 [15]. The Chaobai River was dried before the reclaimed water transfer, so the surface water was not shown.
The confined aquifers of the Chaobai River with a depth of less than 120 m are the most important drinking water source for Beijing [39,41]. There is a well field 4 km off the north of the study area ( Figure 1). The groundwater was also pumped for agricultural activities [39]. In 2000, the pumping yield along the Chaobai River was about 35% of the city's total amount of groundwater exploitation [39]. Beijing's primary depression cone of groundwater level has been observed in the north of the study area ( Figure 1) [42]. Twenty years ago, the groundwater flow direction was consistent with the terrain, namely from north to south [43]. However, it has been modified by groundwater over-exploitation [38]. There is an industrial park for furniture and automobile manufacture in the center of the study area, which may be one potential contaminant source ( Figure 1).

The Chaobai River Restoration Project with RW
Since 2000, the Shunyi reach of the Chaobai river has been dried, with few aquatic organisms surviving. In addition, many residents are living near the Chaobai River. In order to restore the riparian eco-environment and to recover the river landscape for aesthetic and recreational purposes, the Beijing Water Authority implemented a Water Transfer Project in 2007. The Wenyu River water, which is mainly composed of wastewater treatment plant effluent, is further treated by Reclaimed Treatment Water Plant (RTWP) to produce RW. Subsequently, the RW was transferred to the Jian River by underground pipeline (Figure 1). The treating processes of the RTWP consist of ozone pre-oxidation, membrane bioreactor, chemical phosphorus removal, ultraviolet disinfection, constructed wetlands, and oxidation ponds [44]. The concentrations of total nitrogen (TN), total phosphorous (TP), NH 4 + -N, and chemical oxygen demand (CODcr) in the RW are 15, 0.2, 1, and 20 mg/L, respectively. Besides TN, the water quality of the RW achieved the GB/T18921-2002 for Class III criteria [45]. Because of the temperature dependence of water treatment capacity, the monthly transferred volume from May to October (wet season) is 2 × 10 6 -4 × 10 6 m 3 , and that from November to April next year (dry season) is 0.5 × 10 6 -2 × 10 6 m 3 .
Since October 2007, RW started to restore the Chaobai River. The river water is managed by sluices and dams, and it has been divided into several reaches. The reach between S01 and S08 is feed by RW all year round, called the First-stage perennial intake reach ( Figure 1). The RW is transferred to reach between S12 and S05 between June and October every year. From November to May next year, the river becomes dry by evaporation and infiltration. Hence, the reach between S12 and S05 is called intermittent intake reach. The RW was also transferred to upstream of S12 since April 2018, where the river bed was conducted by the anti-seepage treatment with the compacted clay layer, called anti-seepage reach. In 2012, the transferred volume of the RW from RTWP increased by approximately 1.2 × 10 7 m 3 , and the RW was expanded from S08 to S10. The reach between S08 and S10 is also fed by RW all year round, called second-stage perennial intake reach. From 2007 to 2018, the total transferred volume of the RW is 2.78 × 10 8 m 3 (Table S1), and it forms a 7.13 × 10 6 m 2 landscape water area. The compacted clay layer of the anti-seepage reach may prevent the RW infiltration, while the natural channels in other reaches may enhance this process ( Figure 2).

Sampling and Analytical Procedures
Seven groundwater monitoring sections with 35 monitoring well groups (with 30 m, 50 m, and 80 m depth) and 12 surface water monitoring sections form a monitoring network in the study area ( Figure 1 and Table S2). The waters samples were collected at the Intermittent intake reach cross monitoring section 1 and 2 (IM1 and IM2), and at the Perennial intake reach monitoring cross-section 1, 2, 3, and 4 (PM1, PM2, PM3, and PM4), in order to reveal the impact of the RW recharging duration to the river channel on the groundwater hydrodynamics. The water samples were also collected at the Anti-seepage reach monitoring cross-section 1 (AM1) to identify the impact of riverbed permeability on RW transport in the aquifers. In each monitoring section, the groundwater samples were collected from different distances to the river bank in three-depth aquifers (Table S2). Additionally, the seven groundwater monitoring sections are in different sites of the groundwater flow systems with a hydraulic gradient ( Figure 2). Hence, the spatiotemporal variations of groundwater hydrodynamics and the RW transport in the aquifers can be identified. The groundwater levels of the wells were recorded bimonthly or monthly from Sept. 2007 to Dec. 2016. The groundwater levels are reported by elevation in meters above sea level (m.a.s.l.). The Arc Map version 10.2 [46] can be used for conducting groundwater level contours by the Kriging spatial interpolation technique [47].
The groundwater samples that are collected before RW transfer (Sept. 2007) are considered as the groundwater background. In spring, there is negligible precipitation or snow in the study area. To interpret the causes of groundwater dynamics variations and reveal the RW transport in the aquifers, seven sampling campaigns were conducted in the summer after the RW transfer (May 2009, May 2010, May 2013, Mar. 2015, May 2016, Mar. 2017, and May 2018). Groundwater samples were sampled after three to five well volumes of water was pumped out in order to ensure the collection of new groundwater. "U", "SC", and "DC" as prefixes of the sample labels represent groundwater from 30-m, 50-m, and 80-m depth aquifers, respectively. A total of 566 water samples were collected for cation and anion analysis, including 507 groundwater samples and 59 surface water samples (Table S3). The water samples in May 2018 were conducted for δ 2 H and δ 18 O analysis to identify the impact of RW on the aquifers.
Water samples for cation and anion analysis were filtered by 0.45 µm filters during the sampling campaigns. The pH values were measured in situ via a pH meter (WM22EP, TOA-DKK, Japan). The alkalinity was determined within 24 h by titration with 0.02 mol/L hydrochloric acid with methyl orange as an indicator. The hydrochemistry compositions of the water samples were performed in the physical and chemical analysis laboratory of the Institute of Geographic Sciences and Natural Resources Research (IGSNRR), Chinese Academy of Sciences (CAS). Cations were measured by using inductively coupled plasma optical emission spectrometry (ICP-OES) (PerkinElmer Optima 5300, DV, USA). Approximately 96% of the samples have charge balance errors of less than 5%, and 4% of the samples have charge balance errors of less than 10%. The detection limits of potassium permanganate index (COD Mn ) and total hardness (TH) are 0.05 mg/L and 1 mg/L, respectively. The isotopic analysis was performed by a liquid water isotope analyzer (LGR, USA) at the Key Laboratory of Water Cycle and Related Land Surface Process of IGSNRR, CAS, with precisions of ±1‰ and ±0.1‰ for δ 2 H and δ 18 O (1σ), respectively. The results are expressed in part per thousand deviations (‰) relative to the Vienna Standard Mean Ocean Water (V-SMOW).

Mixing Model
Chloride is still the most cost-effective and conservative tracer, which can be used to illustrate the groundwater dynamics [27]. The mixing calculations were conducted by the Cl − concentrations of the samples with a conservative two-end member mixing system [48]. The ratio of RW in each sample is expressed as f RW , while using the formula: where C S , C BG , and C RW refer to the Cl − concentration in the sample, background groundwater, and reclaimed water, respectively. The maximum of Cl − concentration of the reclaimed water outfall (S01) (114.0 mg/L) was selected as the endmember of reclaimed water (C RW ). The Cl − concentration of each groundwater sample before RW recharge to the Chaobai River (samples in 2007) was selected as the background groundwater endmember (C BG ) for each monitoring well.

Hierarchical Cluster Analysis
Hierarchical cluster analysis (HCA) is an efficient way to identify groups of groundwater samples based on their similar hydrogeochemical components, which has been widely used for identifying the groundwater flow path [23,49,50]. In this study, it is used to investigate the relationships between the hydrochemical datasets of the sampling sites in both 2007 and 2018. Na + , K + , Ca 2+ , Mg 2+ , Cl − , SO 4 2− , HCO 3 − , pH, COD Mn , and TH are considered for the HCA. Firstly, all parameters were checked for the skewness, which show positive or negative values. Hence, log-transformation was performed for these parameters (except for pH) to obtain the normal distribution. In addition, the standardization of log-transformed parameters was also performed to eliminate the magnitude of different parameters [50]. Finally, a new value of each parameter with zero mean and a range of −3 to 3 standard deviations was used as the HCA input. In this study, the Euclidian distance, along with the Ward's method for linkage, was conducted to generate clusters. This approach can give the most distinctive groups [50]. These procedures were performed by the statistical software SPSS 22. The mean value of the parameters in each cluster can be calculated, and the Piper and Stiff Diagrams plot the hydrochemical compositions of the clusters.

The 30-m Depth Aquifer
Long-term variations in the groundwater levels can be a good indicator for the interactions between surface water and groundwater [51]. In Sept. 2007, the groundwater levels of the 30-m depth aquifer decreased from 28 m in the southwest to 0 m in the northeast, with a hydraulic gradient of 0.43% (Figure 3a). In May 2018, they decreased from 30 m in the south to 2 m in the north, with a hydraulic gradient of 0.43% (Figure 3d). Generally, the 30-m depth groundwater levels increased by 0-3 m after the Chaobai River was restored by the RW for 11 years.   Hence, the seasonal variation of groundwater level at IM2 before 2012 was controlled by rainfall recharge. Subsequently, it was dominated by the RW infiltration, with the transferred volume of RW increased by 1.2 × 10 7 m 3 in 2012. At anti-seepage reach, the groundwater level at U01 increased from 3.3 to 8.2 m from Jun. 2008 to Sept. 2008, with obviously seasonal variation before 2012. After 2012, the groundwater level at U01 showed weak seasonal variation, and then slowly decreased to 3.5 m in Dec. 2016. Hence, the groundwater level at U01 was controlled by groundwater pumping, precipitation recharge, and RW infiltration. The groundwater levels of the 30-m depth groundwater decrease with increasing distance from the Chaobai River in the same monitoring sections, except for U06 and U05 ( Figure S1). It indicates that a long-term standing groundwater mound formed along the river. The losing regime of river replenished with transferred water was widely identified [20,51,52], showing the decreased effect of the river on groundwater with increasing distance to the river. The groundwater level at DC01 showed a similar variation trend and values with that of SC01, indicating that it was also controlled by groundwater pumping and RW infiltration. In the monitoring period, the groundwater levels at DC22 and DC15 varied from −8.9 to 1.4 m and −9.6 to 1.1 m, respectively. They varied almost simultaneously with time and they were characterized by obvious seasonal variations. The seasonal variations of groundwater levels in different sections increase with depth ( Figure 4). It indicates that the 80-m depth groundwater levels are controlled by groundwater pumping for agricultural irrigation concentrated between March and July, other than RW or rainfall recharge by surface infiltration.
were characterized by obvious seasonal variations. The seasonal variations of groundwater levels in different sections increase with depth ( Figure 4). It indicates that the 80-m depth groundwater levels are controlled by groundwater pumping for agricultural irrigation concentrated between March and July, other than RW or rainfall recharge by surface infiltration. In brief, most of the three-depth groundwater samples are hydrodynamically affected by RW infiltration, as mentioned above. The study area is characterized by multi-layer alluvial aquifers with extreme heterogeneity, resulting in a complex groundwater dynamic system [53]. Hydrogeological, long-term hydrochemical, and isotopic data, together with multivariate statistics, are together used to interpret the groundwater dynamics variations, and to identify the RW transport in the aquifer as follows.

Hydrochemical and Stable Isotopic Compositions of Water Samples
The time-series Cl − concentrations of surface water and groundwater are plotted in Figure 5 to identify the impact of RW on groundwater hydrochemistry (Table S3) (Table S3).  In brief, most of the three-depth groundwater samples are hydrodynamically affected by RW infiltration, as mentioned above. The study area is characterized by multilayer alluvial aquifers with extreme heterogeneity, resulting in a complex groundwater dynamic system [53]. Hydrogeological, long-term hydrochemical, and isotopic data, together with multivariate statistics, are together used to interpret the groundwater dynamics variations, and to identify the RW transport in the aquifer as follows.

Hydrochemical and Stable Isotopic Compositions of Water Samples
The time-series Cl − concentrations of surface water and groundwater are plotted in Figure 5 to identify the impact of RW on groundwater hydrochemistry (Table S3) (Table S3).  Table S4). In 2007, the groundwater samples from different aquifers in the same sections have similar hydrochemical features. The 80-  In 2018, the groundwater samples at AM1 show different hydrochemical compositions from those in 2007 (Figure 6a). When compared with those in 2007, the 50-m and 80-m depth groundwater samples in 2018 at IM1 have higher Cl − and Na + contents, while the 30-m depth groundwater samples have higher HCO 3 − contents (Figure 6b). Most of the 30-m depth groundwater at IM2 and PM1, and the 50-m depth groundwater at IM2 in 2018, are similar to that of the surface water. However, they have obviously different hydrochemical compositions from those in 2007 (Figure 6c,d). In contrast, the 50-m depth groundwater at PM1 and 80-m depth groundwater at IM2 and PM1 have similar hydrochemical compositions to those in 2007.
The δ 2 H and δ 18 O of different water bodies can be applied in order to identify the groundwater origin and potential subsurface hydrological processes in an alluvial aquifer, especially where surface water has undergone extreme evaporation [30,54]. To identify the impact of RW on the three-depth groundwater, the δ 2 H and δ 18 O of water samples in May 2018 are plotted on Figure 7. Figure 6. Stiff diagrams of hydrochemical data for the surface water and groundwater at (a) Antiseepage reach monitoring cross-section 1, (b) Intermittent intake reach monitoring cross-section 1, (c) Intermittent intake reach monitoring cross-section 2, and (d) Perennial intake reach monitoring cross-section 1. Blue, green, and yellow stiff diagrams are plotted for surface water sampled in May 2018, groundwater sampled in September 2007 and May 2018, respectively. The symbols of stratum are the same as in Figure 2.
In 2018, the groundwater samples at AM1 show different hydrochemical compositions from those in 2007 (Figure 6a). When compared with those in 2007, the 50-m and 80m depth groundwater samples in 2018 at IM1 have higher Cl − and Na + contents, while the 30-m depth groundwater samples have higher HCO3 − contents (Figure 6b). Most of the 30m depth groundwater at IM2 and PM1, and the 50-m depth groundwater at IM2 in 2018, are similar to that of the surface water. However, they have obviously different hydrochemical compositions from those in 2007 (Figure 6c,d). In contrast, the 50-m depth groundwater at PM1 and 80-m depth groundwater at IM2 and PM1 have similar hydrochemical compositions to those in 2007.
The δ 2 H and δ 18 O of different water bodies can be applied in order to identify the groundwater origin and potential subsurface hydrological processes in an alluvial aquifer, especially where surface water has undergone extreme evaporation [30,54]. To identify the impact of RW on the three-depth groundwater, the δ 2 H and δ 18 [55]; LMWL: local meteoric water line [56]; and, GWL: groundwater line.

Groundwater Dynamics and Hydrochemistry before RW Transfer
The areas of groundwater recharge and discharge can be identified by the groundwater level contours characteristics before the RW replenished to the Chaobai River [57]. It was obvious that the groundwater level contours in the 30-m depth aquifer become denser in the north, indicating quick groundwater head loss (Figure 3a). The widely distributed gravel sediments in north with well hydraulic conductivity (4.0 × 10 −4 m/s) [15] should not be the case. Hence, it was accounted for the shorter distance to the well field where underwent intensively groundwater abstraction [39] (Figure 1). Additionally, the northward groundwater flow in the three-depth aquifers was opposite with the terrain, indicating the drastic effect of the groundwater over-exploitation in the well field ( Figure 3). The dense level contours in the south of the 50-m depth aquifer might be controlled by the silty clay with weak hydraulic conductivity (0.7 × 10 −4 m/s) [15] (Figures 2 and 3b).
In 2007, most of the groundwater samples in the same layer are clustered together based on the hydrochemical data ( Figure 8). Namely, the clusters A2, A3, and A1 represent the 30-m, 50-m, and 80-m depth groundwater samples, respectively. It shows the trend of decreasing TDS and total hardness (TH) with depth ( Figure 8a and Table 1). COD Mn is a common indicator of organic pollution in water [58]. Cluster A2 shows higher COD Mn (2.1 mg/L), than clusters A1 and A3 (1.3 mg/L), which indicates the potential pollution input into the 30-m depth aquifer.
The 50-m depth groundwater and most 30-m depth groundwater at PM1, PM3, and PM4 with the water type of Ca·Mg-HCO 3 show Cl − concentrations lower than 30 mg/L ( Figure 5). It reflects the local groundwater with a shorter flow path around the recharge zone of the aquifers [37,53]. The water types of 80-m depth groundwater are featured by Ca·Na·Mg-HCO 3 ( Figure S2) and they maintain the hydrochemical composition of regional groundwater with a longer flow path in the deeper layer [39].
Cl − can be selected as a conservative tracer to identify the pollutant source of groundwater before RW recharge to the channel. In 2007, the Cl − concentrations of groundwater decreased with depth ( Figure 5). The U28 has the lowest Cl − concentration (9.7 mg/L) in the 30-m depth aquifer, which was located in upstream of the groundwater system ( Figure 3). Hence, it can be regarded as the background of the 30-m depth groundwater. The Cl − concentrations of most 30-m depth groundwater at AM1, IM1, IM2, and PM2 are higher than 40 mg/L, which is four-fold of that of the U28. The higher Cl − concentrations of these groundwater samples indicate that anthropogenic activities influence the 30-m depth aquifer. Firstly, PM2 is characterized by the highest Cl − concentration (mean 104.2 mg/L) of the 30-m depth aquifer at all monitoring cross-sections. Before the RW recharges to the Jian River, there is a municipal dump around this section. The leaching of municipal dump may be responsible for the higher Cl − concentration at PM2 in the 30-m depth aquifer.
The microbial source tracking that was conducted by Zhang et al. [59] and the distribution of antibiotics conducted by Chen et al. [60] can further confirm the contribution of domestic sewage to the Chaobai River. Secondly, U32 is located within 50 m away from an industrial sewage outlet. It also has the highest Cl − concentration (141.0 mg/L) (Figure 1), which indicates the local industrial contaminant input from furniture and automobile manufacture into the 30-m depth aquifer (Figure 1). It may also explain the higher Cl − concentrations of the 30-m depth groundwater at IM2 located in the west of the industrial park. The distribution of polychlorinated biphenyls in the sediment samples that were measured in the previous study can provide evidence of the point-source pollution from local factories [61]. The 30-m depth groundwater samples at AM1 are featured by higher Cl − (40.3 mg/L), K + (13.9 mg/L), and NO 3 -(13.6 mg/L) concentrations (Table S4), demonstrating the agricultural pollution to the 30-m depth aquifer at AM1 [62,63]. The 50-m depth groundwater and most 30-m depth groundwater at PM1, PM3, and PM4 with the water type of Ca·Mg-HCO3 show Cl − concentrations lower than 30 mg/L ( Figure 5). It reflects the local groundwater with a shorter flow path around the recharge zone of the aquifers [37,53]. The water types of 80-m depth groundwater are featured by Ca·Na·Mg-HCO3 ( Figure S2) and they maintain the hydrochemical composition of regional groundwater with a longer flow path in the deeper layer [39]. Cl − can be selected as a conservative tracer to identify the pollutant source of groundwater before RW recharge to the channel. In 2007, the Cl − concentrations of groundwater decreased with depth ( Figure 5). The U28 has the lowest Cl − concentration (9.7 mg/L) in

The Impact of the RW on Different Aquifers
The variety of recharge sources of the aquifers could have different responses to the water transfer project [19]. Hence, the origin of groundwater in different aquifers is identified to recognize the impact of RW on groundwater. The obviously different stable isotopic compositions of three-depth groundwater indicate their different origins (Figure 7 and Figure S3). Most of the 30-m depth groundwater samples are plotted below and far away from LMWL, showing similar δ 18 O and δ 2 H values with surface water (Figure 7). It indicates that the isotopically enriched surface water may recharge into the 30-m depth aquifer. In addition, a few 30-m depth groundwater samples, such as U05, U12, and U32, are close to LMWL, which suggests that the 30-m depth groundwater far away from the river channel is dominated by the precipitation. However, the d-excess values of the 30-m depth groundwater (mean 1.6‰) are lower than that of surface water (mean 2.0‰) in May 2018, which indicates the directed evaporation of some 30-m depth groundwater samples in the riparian zone. Their groundwater level depths are around 2 m, which is lower than the phreatic limit evaporation depth (4 m) in north China [64]. The irrigation by pumping groundwater and the associated return flow around the anti-seepage reach could also be responsible for the lower d-excess values of the 30-m depth groundwater [63]. It is also verified by the higher K + (13.9 mg/L) and NO 3 − (13.6 mg/L) of the 30-m depth groundwater at AM1 (Table S4), which may result from the fertilizer input that accompanies the irrigation return flow [62][63][64].
The δ 18 O of 80-m depth groundwater (mean −8.3‰) is more depleted than the weighted annual local precipitation (−7.2‰) [65] and surface water (mean −6.1‰). It suggests that the recharge of the 80-m depth aquifer could take place during the cooler climatic conditions, such as in the high mountains or the cold periods of the Late Pleistocene [66]. However, the residence time of groundwater within 100 m depth in the study area was 15-35 years, as identified by 3 H [53]. Additionally, the δ 18 O of paleo-water in the North China Plain ranges from −11.6 to −9.5‰, which is more depleted than that of the 80-m depth aquifer [67]. It can be inferred that the 80-m depth groundwater may be dominated by the lateral groundwater flow, with modern precipitation being recharged around the north or northwest mountain areas [39] (Figure 1). However, the 80-m depth groundwater samples around the depression cone, such as DC01, DC02, DC04, and DC07, are characterized by enriched δ 18 O and δ 2 H, as well as lower d-excess values (Figure 7 and Figure S3). It could be explained by the vertical mixing of the 80-m depth groundwater with 30-m depth groundwater, which is recharged by the isotope-enriched surface water [53], accelerated by the increased vertical hydraulic gradient.
In Figure 7 and Figure S3, the 50-m depth groundwater was distributed between the 30-m and 80-m depth groundwater. The stable isotope compositions of the 50-m depth groundwater (δ 2 H −58‰) are similar to that of the 30-m depth groundwater (δ 2 H −56‰), indicating the close hydrological link between the 30-m and 50-m depth aquifers. The 50-m depth groundwater may also be recharged by surface water. Generally, the groundwater is characterized by a decreased renew rate with depth (15-35a), which is affected by the depression cone in groundwater [53].

The 30-m Depth Aquifer
The hydrochemical compositions of the 30-m depth groundwater samples at IM2 and PM1 in 2018 are very similar to that of surface water, with the increased Cl − , Na + and decreased Ca 2+ , HCO 3 − concentrations when compared with those in 2007 (Figure 6c,d). The 30-m depth groundwater samples at IM2 and PM1 are clustered into B4 and B3, showing similarities with RW ( Figure 8b). Additionally, the Cl − concentrations of most 30-m depth groundwater samples increase with time at IM2 and PM1 (Figure 5c,d). It indicates that the hydrochemistry of the 30-m depth groundwater samples at IM2 and PM1 are obviously affected by the RW. The Cl − concentrations of most 30-m depth groundwater samples at PM2, PM3, and PM4 present the synchronous variations with that of surface water from 2015 to 2018 (Figure 5f,g), revealing the hydrochemical impact of the RW on these 30-m depth groundwater samples.
The Cl − concentrations of most 30-m depth groundwater samples at AM1 increase with time, reflecting the influence of the RW (Figure 5a). However, the hydrochemical compositions of 30-m depth groundwater samples at AM1 are different from that of surface water (S12) (Figure 6a), indicating that the 30-m depth aquifer is dominated by the regional groundwater and slightly influenced by the RW. It is in line with the 30-m depth groundwater samples at AM1 being clustered into B1, showing large dissimilarities with RW (Figure 8b). At IM1, the 30-m depth groundwater far away from the river channel (U06 and U12) in 2018 shows higher Ca 2+ and HCO 3 − concentrations than that close to the river channel (U08) (Figure 6b). U06 and U12 also have more depleted isotopic compositions than that of surface water (Figure 7). Hence, the RW only locally affects the groundwater close to the river channel (U08) at IM1, which has similar hydrochemical compositions with the surface water (S12).
There are two potential RW sources in the 30-m depth aquifer at AM1, including the vertical infiltration of RW through the riverbed at AM1 and the lateral groundwater flow from the Intermittent and Perennial intake reaches. When considering the anti-sewage treatment of the riverbed at AM1, the former source of RW is excluded. It is in line with the different hydrochemical compositions of 30-m depth groundwater and surface water at AM1 in 2018 (Figure 6a). The similar variation trend of the groundwater levels at U15 and U01 reveals the lateral hydraulic connectivity of the 30-m depth aquifer between AM1 and IM2 (Figure 4a,b). Hence, the later source could contribute to the RW in the 30-m depth aquifer at AM1.
The groundwater pumping led to the increased lateral hydraulic gradient (0.43%) in 30-m depth aquifer before reclaimed water transfer (Figure 3a). It does not change (0.43%) in 2018 (Figure 3d), even though the RW has been recharged the aquifer for 11 years through the riverbed infiltration. The increased lateral hydraulic gradient could enhance the transport of the RW in the aquifers. Zhu et al. also found that the groundwater pumping greatly determines the water head differences between the river water and groundwater, which subsequently impacts the mixing of river water and groundwater significantly [68]. In the Yellow River, water was transferred to a river, the larger longitudinal hydraulic gradient in the underlying aquifer (0.46%) confines the hydrochemical impact of the transferred water to the riparian zone [20]. It can be concluded that the groundwater dynamic is one of the decisive factors for the RW transport in the aquifers.

The 50-m Depth Aquifer
The Cl − concentrations in most of 50-m depth groundwater samples increase with time, indicating the potential impact of RW ( Figure 5). However, only the 50-m depth groundwater samples at IM2 are clustered into B3, showing similarities with RW ( Figure 8b). Additionally, the groundwater flow in the 50-m depth aquifer between PM1 and IM2 inversed from northward in 2007 to southward in 2018 (Figure 3b,e). This may be explained by the vertical leakage from the 30-m depth aquifer to the 50-m depth aquifer.
The downward leakage from the 30-m to 50-m depth aquifer can be accelerated by the groundwater pumping in the latter. It resulted in the 30-m depth groundwater levels being much higher than those in the 50-m depth aquifer. In the center and coastal of North China Plain, the deep groundwater pumping has resulted in the leakage from the 30-m depth to the confined aquifer, which accounted for 70% of total withdrawals from the deep aquifer [69,70]. Additionally, geologic heterogeneity could strongly affect the movement of water and solutes through the subsurface [71]. The interconnectedness of the sand body could control the groundwater flow paths in the multiple-aquifer system [72]. During the managed aquifer recharge (MAR), the interconnected coarse texture also resulted in the rapid recharge to the aquifer [30]. The vertical leakage from the 30-m to 50-m depth aquifer could be speculated when considering the thin clay layer (4 m) between the 30-m and 50-m depth aquifers at IM2 (Figure 6d) and the increased vertical hydraulic gradient by the confined groundwater pumping.
At PM1, the hydrochemical compositions of the 50-m depth groundwater in 2018 are similar to those in 2007 (Figure 6d). They are included in cluster B1, showing different hydrochemical features with the RW (Figure 8b). The increasing trend of Cl − concentrations in groundwater of the 50-m depth aquifer is slower than that of the 30-m depth aquifer (Figure 5d). Therefore, the 50-m depth groundwater at PM1 could mix with the RW, but it might be dominated by the regional groundwater. The groundwater levels of U22 and SC22 exhibit obviously different variations (Figure 4c), demonstrating their weak hydraulic connectivity, due to the 15-m thick clay layer between the 30-m and 50-m depth aquifers (Figure 6d). However, similar variations in the groundwater levels at SC15 and SC22 indicate the well lateral hydraulic connectivity of the 50-m depth aquifer between IM2 and PM1. Thus, the RW in the 50-m depth aquifer at PM1 is laterally migrated from IM2 under a hydraulic gradient (Figure 3e), other than by vertical leakage from the 30-m depth aquifer at PM1.
The groundwater levels of U01 and SC01 are characterized by different variation trends (Figure 4a). Additionally, the variations of Cl − concentrations in the 30-m and 50-m depth aquifers at AM1 and IM1 (Figure 5a,b) reflect the weak hydraulic connectivity between the two aquifers. At AM1 and IM1, the Cl − concentrations in the 50-m depth aquifer are higher than that in the 30-m depth aquifer. Hence, similar to PM1, the RW from IM2 also laterally recharge to the 50-m depth aquifer at AM1 and IM1. In general, the RW vertically recharge to the 50-m depth at the intermittent intake reach by leakage, due to the heterogenous properties of the alluvial aquifer. Subsequently, RW gradually migrates to other sites in a horizontal direction by groundwater advection under a hydraulic gradient.

The 80-m Depth Aquifer
The 80-m depth groundwater samples in 2018 at most monitoring sections have similar hydrochemical compositions to those in 2007 ( Figure 6). They are all clustered into B1 and B2, showing minimal similarities with that of the RW (Figure 8b). Hence, the RW has less impact on the 80-m depth aquifer. It can be confirmed by the lowest endocrine disrupting chemicals detection frequency and concentrations observed in the 80-m depth groundwater [40]. The Cl − concentrations of most of the 80-m depth groundwater samples at AM1 and IM1 increase the time ( Figure 5). However, it is not the case at other monitoring cross-sections. The approximately same values of groundwater level in SC01 and DC01 indicate the well hydraulic connectivity between the 50-m and 80-m depth aquifers at AM1 (Figure 4a). It is in line with the relatively thin silty clay between them at AM1 (around 1.5 m) (Figure 6a). Hence, the RW may leak into the 80-m depth aquifer from the 50-m depth aquifer at AM1, which could be accelerated by the high permeability of the aquifer [19,53].

Proportion of Reclaimed Water in Groundwater
Chloride, as a conservative tracer, has been widely applied to trace the mixing of different water sources, but the leaching of Cl − from the soil could be another source to groundwater during the RW infiltration. The Cl − concentrations of topsoil (at 20-30 cm depth) of the study area in 2019 range from 7 to 53 mg/L (mean 25 mg/L) [73], which is significantly lower than that (114.0 mg/L) of the RW. Hence, chlorine can be selected to calculate the ratio of RW in groundwater (Figure 9). Besides RW, there are other sources of Cl − inputting into groundwater, such as industrial discharge. Only the groundwater from one well U32 that was obviously impacted by the industrial pollutant was sampled once before the RW transferred. It is difficult to identify the industry endmember. The Cl − concentration of groundwater higher than 60 mg/L in 2007 can be excluded for the calculation conservative mixing model for the obvious pollution in order to eliminate the uncertainty of industrial discharge.
In general, the 11-year average proportion of the RW in the 30-m, 50-m, and 80-m depth aquifers are 53%, 39%, and 15%, respectively. The proportion of the RW in the 30-m and 50-m depth groundwater increases fast than that in the 80-m depth groundwater, indicating the well hydraulic connectivity between the 30-m and 50-m depth aquifers with surface water. In the 30-m depth aquifer, the proportion of RW decreases from the perennial intake reach (PM1) (average 66%) to the intermittent intake reach (IM2) (average 55%), and to the anti-seepage reach (AM1) (average 34%). The higher proportion of RW in the 30-m depth groundwater at PM1 than that at IM2 could be explained by the larger transfer volume of RW in the perennial intake reach. The 50-m depth aquifer at IM2 has a relatively higher proportion of the RW (average 47%) than that at the PM1 (average 35%), which reflects the preferential recharge of the RW into the 50-m depth aquifer around IM2. In contrast, the proportion of the RW in the 80-m depth aquifer decreases from the AM1 (average 23%), to the IM2 (average 16%), and to the PM1 (average 4%), which suggests the leakage recharge of the 80-m depth aquifer around AM1. DC01 indicate the well hydraulic connectivity between the 50-m and 80-m depth aquifers at AM1 (Figure 4a). It is in line with the relatively thin silty clay between them at AM1 (around 1.5 m) (Figure 6a). Hence, the RW may leak into the 80-m depth aquifer from the 50-m depth aquifer at AM1, which could be accelerated by the high permeability of the aquifer [19,53].

Proportion of Reclaimed Water in Groundwater
Chloride, as a conservative tracer, has been widely applied to trace the mixing of different water sources, but the leaching of Cl − from the soil could be another source to groundwater during the RW infiltration. The Cl − concentrations of topsoil (at 20-30 cm depth) of the study area in 2019 range from 7 to 53 mg/L (mean 25 mg/L) [73], which is significantly lower than that (114.0 mg/L) of the RW. Hence, chlorine can be selected to calculate the ratio of RW in groundwater (Figure 9). Besides RW, there are other sources of Cl − inputting into groundwater, such as industrial discharge. Only the groundwater from one well U32 that was obviously impacted by the industrial pollutant was sampled once before the RW transferred. It is difficult to identify the industry endmember. The Cl − concentration of groundwater higher than 60 mg/L in 2007 can be excluded for the calculation conservative mixing model for the obvious pollution in order to eliminate the uncertainty of industrial discharge. In general, the 11-year average proportion of the RW in the 30-m, 50-m, and 80-m depth aquifers are 53%, 39%, and 15%, respectively. The proportion of the RW in the 30m and 50-m depth groundwater increases fast than that in the 80-m depth groundwater, indicating the well hydraulic connectivity between the 30-m and 50-m depth aquifers with surface water. In the 30-m depth aquifer, the proportion of RW decreases from the perennial intake reach (PM1) (average 66%) to the intermittent intake reach (IM2) (average 55%), and to the anti-seepage reach (AM1) (average 34%). The higher proportion of RW in the 30-m depth groundwater at PM1 than that at IM2 could be explained by the larger transfer volume of RW in the perennial intake reach. The 50-m depth aquifer at IM2 has a relatively higher proportion of the RW (average 47%) than that at the PM1 (average 35%), which reflects the preferential recharge of the RW into the 50-m depth aquifer around IM2. In contrast, the proportion of the RW in the 80-m depth aquifer decreases from the AM1

Conceptual Model for Groundwater Flow Systems Restored by RW
The RW has been transferred to the Chaobai River for the river ecosystem restoration, where the depression cone in groundwater level had been formed by the groundwater over-exploitation. The groundwater system in the alluvial plain is largely constrained by the pumping in addition to the heterogeneity of hydrogeological conditions, which determined the hydraulic and transport properties [68]. The higher proportion of the RW in the 50-m depth aquifer at the intermittent intake reach and in the 80-m depth aquifer at the anti-seepage reach indicated that the vertical leakage recharge controls the RW transport in confined aquifers.
A conceptual model of the groundwater flow system in the alluvial aquifers beneath the Chaobai River can be proposed, as shown in Figure 10. A total of 2.78 × 10 8 m 3 RW has been transferred to the nature dry channel since October 2007. RW is infiltrated to the aquifers at the Perennial intake reach all of the year, and at the intermittent intake reach during the RW periodically transferred into this reach. Subsequently, the RW migrates to the anti-seepage reach in the 30-m depth aquifer with a hydraulic gradient. The infiltrated RW forms a long-term standing groundwater mound along the river channel. The higher Cl − concentrations and isotope-enriched compositions of RW are recorded in the groundwater beneath. The direct evaporation of groundwater could be associated with irrigation return flow, resulting in the lower d-excess values of the 30-m depth groundwater than that of surface water. In addition, the industrial effluent with higher Cl − concentration (141.0 mg/L) locally affects the 30-m depth groundwater. The vertical leakage as a preferential recharge from the 30-m depth aquifer to 50-m depth aquifer at the intermittent intake reach controls the RW transport in the 50-m depth aquifer. It results in the relatively higher Cl − concentration (93.0 mg/L) and enriched δ 18 O (−7.2‰) of the 50-m depth groundwater at the intermittent intake reach in May 2018. Alarmingly, few RW also recharge the 80-m depth aquifer from 50-m depth aquifer at the anti-seepage reach. The confined groundwater can be recharged by the lateral groundwater flow from the mountain area [53]. The depression cone of the groundwater levels caused by groundwater pumping modifies the groundwater flow. The well field became a local discharge area, enhancing the RM transport. The declined groundwater levels would decrease the dilution capacity of groundwater to pollutants [36]. Hence, the northward transport of the RW with higher ionic components may threaten the groundwater safety of the well field.  (19) irrigation return flow. SW: surface water; RW: reclaimed water; B1, B2, B3, and B4 represent the clusters generated by HCA using hydrochemistry data (Figure 8).

Conclusions
The RW has been restored the dried Chaobai River since 2007. The groundwater dynamic variations and the transport of the RW in the multi-layer alluvial aquifers have been tracked with eight-year hydrological, hydrochemical and isotopic data, by the conservative mixing model and HCA methods. The major conclusions are as follows:  The impact of RW infiltration on groundwater dynamics shows significantly spatiotemporal variation. The 30-m depth aquifer at Perennial intake reach increased by 3 m in four months and then kept stable, which indicated that they were dominated by RW infiltration. However, the 30-m depth groundwater levels at intermittent intake and anti-seepage reaches were controlled by precipitation recharge before 2012, then they were dominated by RW infiltration with more RW transferred to the river. The 30-m depth groundwater levels decreased with increasing distance to the river, showing a decrease in the effect of the RW on groundwater. The 50-m and 80-m depth groundwater levels decreased by 6-9 m and 6-15 m, respectively, being dominated  (19) irrigation return flow. SW: surface water; RW: reclaimed water; B1, B2, B3, and B4 represent the clusters generated by HCA using hydrochemistry data ( Figure 8).

Conclusions
The RW has been restored the dried Chaobai River since 2007. The groundwater dynamic variations and the transport of the RW in the multi-layer alluvial aquifers have been tracked with eight-year hydrological, hydrochemical and isotopic data, by the conservative mixing model and HCA methods. The major conclusions are as follows:

•
The impact of RW infiltration on groundwater dynamics shows significantly spatiotemporal variation. The 30-m depth aquifer at Perennial intake reach increased by 3 m in four months and then kept stable, which indicated that they were dominated by RW infiltration. However, the 30-m depth groundwater levels at intermittent intake and anti-seepage reaches were controlled by precipitation recharge before 2012, then they were dominated by RW infiltration with more RW transferred to the river. The 30-m depth groundwater levels decreased with increasing distance to the river, showing a decrease in the effect of the RW on groundwater. The 50-m and 80-m depth groundwater levels decreased by 6-9 m and 6-15 m, respectively, being dominated by the groundwater pumping.

•
The RW has a significant impact on the water stable isotopic compositions of 30-m depth aquifer. However, regional groundwater with meteoric origin mainly recharges 50-m depth and 80-m depth groundwater. The RW usage significantly increases the Na + and Cl − concentrations in the groundwater. The groundwater types of the 30-m and 50-m depth aquifers change from Ca·Mg-HCO 3 in 2007 to Na·Ca·Mg-HCO 3 ·Cl and Ca·Na·Mg-HCO 3 in 2018, while that of the 80-m depth aquifer does not change (Ca·Na·Mg-HCO 3 and Na·Mg·Ca-HCO 3 ). The chloride conservative mixing model shows that the averaged proportion of the RW in 30 m, 50 m, and 80 m-depth aquifers are 53%, 39%, and 15%, respectively.

•
Our study confirms that the heterogenous properties of the multi-layer alluvial aquifer offer the preferential flow path for RW transport in the aquifers. The RW mainly infiltrates into the 30-m depth aquifer around the intermittent and perennial intake reaches. However, the RW recharges to the 50-m and 80-m depth aquifers by leakage in the intermittent intake reach and anti-seepage reach, respectively. This leads to the higher mixing ratio of the RW in the confined groundwater at intermittent intake reach and anti-seepage reach than that of the perennial intake reach, where there are more RW in the river channel.

•
The RW utilization can significantly alleviate the local water shortage. However, the increased hydraulic gradient between surface water and groundwater by groundwater pumping could enhance the RW transport in the aquifers. The shorter residence time of RW in the aquifers may restrain the removal rate of pollutants in the RW. It would pose a potential pollution to the groundwater. The restriction of groundwater pumping could decrease the hydraulic gradient. Hence, the RW transfer to the river channel and groundwater exploitation should be considered together for water source management in northern China and other similar areas around the world.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4 441/13/6/806/s1, Table S1: Recharge volume of reclaimed water transferred to the Chaobai River from 2007 to 2018; Table S2. Detailed information of seven monitoring sections; Table S3. Detailed information of seven monitoring sections; Table S3. Average chlorine concentration with a standard deviation of the water samples collected from Chaobai River alluvial plain from 2007 to 2018; Table  S4. Average hydrochemical and isotopic data of the water samples collected from the Chaobai River alluvial plain in September 2007 and May 2018; Figure S1. Hydrographs for well groups with screened sections in the 30-m depth aquifer at a) Anti-seepage reach monitoring cross-section1(AM1), b) Intermittent intake reach monitoring cross-section 1(IM1), c) Intermittent intake reach monitoring cross-section 2 (IM2), and d) Perennial intake reach monitoring cross-section 1(PM1); Figure S2 Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.