Numerical modeling of changes in groundwater storage and nitrate load in the unconfined aquifer near a river receiving reclaimed water

Reclaimed water (RW) has been widely used as an alternative water resource to recharge rivers in mega-city Beijing. At the same time, the RW also recharges the ambient aquifers through riverbank filtration and modifies the subsurface hydrodynamic system and hydrochemical characteristics. To assess the impact of RW recharge on the unconfined groundwater system, we conducted a 3D groundwater flow and solute transport model based on 10 years of sequenced groundwater monitoring data to analyze the changes of the groundwater table, Cl− loads, and NO3-N loads in the shallow aquifer after RW recharge to the river channel. The results show that the groundwater table around the river channel elevated by about 3–4 m quickly after RW recharge from Dec. 2007 to Dec. 2009, and then remained stable due to the continuous RW infiltration. However, the unconfined groundwater storage still declined overall from 2007 to 2014 due to groundwater exploitation. The storage began to recover after groundwater extraction reduction, rising from 3.76 × 108 m3 at the end of 2014 to 3.85 × 108 m3 at the end of 2017. Cl− concentrations varied from 5–75 mg/L before RW recharge to 50–130 mg/L in 2 years (2007–2009), and then remained stable. The zones of the unconfined groundwater quality affected by RW infiltration increased from 11.7 km2 in 2008 to 26.7 km2 in 2017. Cl− loads in the zone increased from 1.8 × 103 t in 2008 to 3.8 × 103 t in 2017, while NO3-N loads decreased from 29.8 t in 2008 to 11.9 t in 2017 annually. We determined the maximum area of the unconfined groundwater quality affected by RW, and groundwater outside this area not affected by RW recharge keeps its original state. The RW recharge to the river channel in the study area is beneficial to increase the groundwater table and unconfined groundwater storage with lesser environmental impacts.


Introduction
The increasing demands for water resources, driven by a combination of population growth, socio-economic development, rapid urbanization, and climate change have led to water scarcity and quality deterioration (USEPA and USAID 2012;WWAP 2019). In response to the ongoing enormous stresses and challenges, reclaimed water (RW) treated from wastewater has been increasingly being used as an alternative water resource for a variety of purposes (Asano et al. 2007;Ahuja 2015;Angelakis et al. 2018). RW is utilized for agriculture, landscape irrigation, industrial reuse, municipal water, groundwater recharge, and even drinking water in many countries (De Gisi et al. 2017;Tortajada and van Rensburg 2020), such as the USA (Hess and Collins 2019), EU (Angelakis and Gikas 2014), Australia (Barry et al. 2017), China (Zhu and Dou 2018), etc. It provides a costeffective and ecologically sustainable solution to the problem of water shortages.
In China, RW utilization was 2.17 × 10 9 m 3 in 2009, accounting for 0.36% of the total water supply amount (Zhu and Dou 2018). Then, it was up to 8.74 × 10 9 m 3 in 2019, accounting for 1.5% of the total water supply (MWR 2021). Beijing, as a mega-city with serious water shortages characterized by the drying of some rivers with a degraded ecosystem, is the biggest RW user in China (Zhang et al. 2014;MWR 2021). Beijing's RW utilization was 6 × 10 8 m 3 in 2008, accounting for 17% of the total water supply amount, while it was up to 2 × 10 9 m 3 in 2020, accounting for approximately 30% of the total water supply amount (Beijing Water Authority 2020). Over 90% of the RW was primarily used for purposes of recreation and ecological conservation (e.g. water supply for recharging rivers or lakes) (Beijing Water Authority 2019). But compared with other natural water types, RW usually contains higher levels of salts, nitrogen, phosphorus, heavy metals, organic, and biological contaminants, which are controlled by water treatment technologies (Deng et al. 2019;Valhondo et al. 2020). As a result, groundwater may have potential health risks and environmental impacts as RW seeps into the aquifer through river leakage (Asano and Cotruvo 2004;Bekele et al. 2018).
RW recharge changes the interaction of surface-and ground-water and related hydraulic features (Guo et al. 2019;Yuan et al. 2020). Hydrochemical components of groundwater could be modified by the RW mixing with native groundwater along the flow path (Yu et al. 2017;Gilabert-Alarcón et al. 2018;Daesslé et al. 2020). Metalloids (arsenic) (Fakhreddine et al. 2021) and metals (Davranche and Bollinger 2000;Patterson et al. 2010) in sediments may be released into groundwater during fluid-rock interaction, bringing risks and challenges for water resources management. It has been demonstrated that some contaminants' concentrations decreased during the RW infiltration process (e.g., filtration, adsorption, redox reaction, and degradation), such as nutrients (Singer and Brown 2018), total organic carbon (Bekele et al. 2011), heavy metals (Zhang et al. 2020), trace organic compounds Ding et al. 2020), and pathogens (Page et al. 2010). About 85% of the nitrate was observed to be removed from the riverbed sediments in a river recharged by RW in Beijing (Pan et al. 2018). However, some studies indicate that agricultural irrigation with RW has resulted in nitrate pollution of groundwater Candela et al. 2007). The variational attenuation of NO 3 -N depends on the utilization methods and the hydrogeological conditions of the recharge fields.
The Shunyi reaches of the Chaobai River in Beijing have been dried since 1999, and RW has been utilized for recharging the river since December 2007 to restore the water ecology of the dry river. Many studies have been carried out on characterizing variations in physical and chemical parameters of groundwater after RW recharge in the study area, with the most significant changes in the unconfined aquifer. The groundwater table increased by ~ 3 m after RW supply to the Chaobai River from 2007 to 2018 (He et al. 2021). Column experimental results indicated that cation exchange took place between K + in the RW and Ca 2+ in the riverbed medium during RW infiltration (Liu et al. 2018). The hydrochemical type of the unconfined groundwater has been modified from HCO 3 -Ca·Mg in 2007 to HCO 3 ·Cl-Na·Ca·Mg in 2018. The latter is similar to the RW type (HCO 3 ·Cl-Na·Ca and HCO 3 ·SO 4 ·Cl-Na·Ca) Jiang et al. 2020;He et al. 2021). It has been observed that the salinity of the unconfined groundwater increased and the total hardness decreased (Liu et al. 2018). In the column experiments and field monitoring, nitrogen   (Liu et al. 2018;Li et al. 2019) and some endocrinedisrupting chemicals (Li et al. 2013;Ma et al. 2015;Wang et al. 2019) were observed to be removed during the infiltration process with different removal rates. Nitrogen, especially nitrate, is the typical contaminant in the RW, river, and unconfined groundwater in this area (Yang et al. 2016;Li et al. 2019). However, few studies have focused on quantifying the effects of RW recharge on the groundwater system.
The objectives of this study are to (i) conduct a 3D groundwater flow and solute (Cl − ) transport model, (ii) calculate the changes of unconfined groundwater storage, Cl − loads, and NO 3 -N loads, and (iii) identify the area of groundwater quality affected by the RW recharge and illustrate the environmental impacts in the study area. It can improve understanding of the effects of RW recharge to a river channel on the groundwater system and provide a scientific reference for the safe use of RW and water resources management in other similar regions globally.

Study area
The study area is northeast of Beijing, between longitudes 116°36′ to 116°46′ E and latitudes 40°02′ to 40°12′ N (Fig. 1a). Elevation ranges between 22 and 54 m.a.s.l. The area is in the northern temperate region, with a typical continental monsoon climate. The annual average temperature is ~ 11.5 °C. The mean annual precipitation is about 610 mm, chiefly from June to September. The average annual evaporation based on pan evaporation experiments is about 1108 mm.
The Chaobai River channel was dry before 2007. The Beijing Municipal Water Diversion Project was implemented to restore the dry river channel by the Beijing Water Authority since 2007. The Wenyu River water was pumped and treated by a membrane bioreactor (MBR) to produce RW. Then, the RW was transferred into the Jian River, and finally flowed into the Chaobai River (Fig. 1a). The transferred volume of RW is 2.9 × 10 7 m 3 in 2017, and the cumulative transferred volume is 2.3 × 10 8 m 3 from 2007 to 2017. RW has replenished the Jian River and Shunyi reach of the Chaobai River since October 2007. The Jian River is an artificial channel with a natural river bottom and sloping brick sides. It is about 4 km long, 50-90 m wide. The Chaobai River in the study area is a natural channel. It is about 15 km long, 200-400 m wide. Dams were set up to control the river water level and separate reaches of the river. The Jian River and Chaobai River hold approximately 5.1 × 10 5 m 3 and 8.91 × 10 6 m 3 of water under the control of a sluice and four dams, respectively. The section from the Earth dam to the Henan rubber dam (Fig. 1a) is a perennial water-receiving river channel. The section between the Xiangyang sluice and the Earth dam is replenished only in May and October as an intermittent water-receiving river channel. Moreover, the section from the Henan rubber dam to the Suzhuang rubber dam was perennially replenished since 2012.

Hydrogeology conditions
The study area is in the lowland of the Chaobai alluvialdiluvial fan where Quaternary deposits are widely distributed (Li et al. 2013). The hydrogeological profile (Fig. 1b) shows the stratigraphic distribution. From north to south, sediment grains gradually change from gravel to fine sand and to silty clay, from coarse to fine. The vertical profile is characterized by four aquifers (gravel and sand) and three aquitards (silty clay). The first aquifer (layer I) is unconfined and mainly composed of gravel and fine sand, with a buried depth of 0-30 m. The second aquifer (layer III) consists of gravel and sand with a buried depth of 30-50 m. The third aquifer (layer V) and the fifth aquifer (layer VII) are mainly dominated by fine sand, with a buried depth of 50-80 m and below 80 m, respectively.
The groundwater is mainly recharged by rainfall, river leakage, lateral inflow, and irrigation seepage through the upper aquifer. Groundwater discharge is mainly via the extraction of groundwater below 80-m depth. The average depth of the groundwater table in the Xiangyang sluice (G01, G02, and G03) was 29.7 m in 2017, and the average depth for other observation wells in Fig. 1a was 4.4 m in 2017. The groundwater generally flows from southwest to northeast, which is the opposite of the terrain due to excessive groundwater exploitation (Zheng et al. 2015).

Data sources
Details on input data are summarized in Table 1. From December 2007 to December 2017, all data were processed into a monthly dataset, analyzed for the purpose of developing a conceptual model, and finally converted into a numerical model for groundwater flow and solute transport simulation.

Groundwater flow and solute transport modeling
Numerical simulations of groundwater flow and solute transport in this study area were performed using MOD-FLOW (McDonald and Harbaugh 1988) and MT3DMS The model was conducted with a total simulation time of 241 months, as stress periods. The model domain was discretized into 150 columns and 170 rows with a grid size of 100 m × 100 m. The total thickness of the aquifer system is 200 m including four aquifers and three aquitards (Fig. 1b). The aquifer system was divided into seven model layers to represent each aquifer and aquitard. The top surface elevation was determined using the Digital Elevation Model (DEM). The bottom of each layer was generated by interpolating the boreholes data. Further analysis and discussion in this study focus on the unconfined aquifer (first layer) because this aquifer has been significantly affected by RW relative to other aquifers. According to the monitoring data on the groundwater levels and the map of the groundwater level contour (He et al. 2021), the no-flow boundary was applied to the northern boundary AB, which is almost perpendicular to the groundwater level contour. And the general head boundary applied to the eastern boundary BC and western boundary AD. The south boundary CD was characterized by a constant head boundary. The top boundary was the free surface of the unconfined aquifer. The lowest boundary is the bottom of the seventh layer, most of which are limestone and dolomite with low permeability (4.6 × 10 −11 -1.2 × 10 −9 of hydraulic conductivity) (Louis 1972). This boundary is generalized to the confining boundary.
Water fluxes across the top boundary include recharge from the precipitation infiltration, leakage recharge from the river channel, and discharge through the leakage of the unconfined aquifer and evaporation. Recharge from precipitation is calculated as the product of the monthly values of precipitation in the period 2007-2017 and the local precipitation infiltration coefficient (Zheng 2012). The river leakage recharge, including the Chaobai River and Jian River, was calculated through water balance in the river channel: where Q R is the river leakage recharge, Q RW is RW transfer volume, Q P is precipitation recharge, Q E is evaporation from the water surface, Q T is the transferred river water to upstream, and the unit is m 3 /a. Due to RW restoring the river channel and forming a large surface water area (approximately 5.63 km 2 ), river leakage recharge was regarded as planar recharge, the same as precipitation infiltration. The temporal recharge rate in the upper layer can be estimated as a combination of rainfall and river leakage. Groundwater is mainly pumped from the aquifers below 80-m depth. The well fields were simulated (1) as pumping wells (Zheng 2012), and their discharge rates were assigned on the fifth and seventh layers. The amount of groundwater evaporation was very small due to the low groundwater table. When the groundwater table depth in the North China Plain is greater than 3 m, evaporation of the unconfined groundwater disappears (Huo 2015). Thus, the evaporation was simulated in the zones where the groundwater table depth is shallower than 3 m. The observed groundwater levels in the unconfined aquifer and the confined aquifers at the end of 2007 were interpolated as the initial groundwater levels.
The solute transport model was developed based on the groundwater flow model. Cl − as a conservative solute has been generally used to trace and evaluate the impact of RW recharge on subsurface flow systems (Yu et al. 2017;Zhang and Yu 2021). Cl − concentrations of RW are usually higher than that of natural water.
The Cl − concentration in the unconfined groundwater increased after RW restoring the river channel. Besides RW, industrial pollution in the study area is also the source of Cl − inputting into groundwater. But the latter didn't affect the unconfined groundwater surrounding the observation wells except for G32 (He et al. 2021). The high Cl − concentrations in the unconfined groundwater along the river are mainly affected by the RW infiltration. Thus, river leakage was assumed to be the only source of high concentrations of Cl − in the model. The river was treated as a constant concentration boundary condition with Cl − concentration of 60-120 mg/L according to the observed values. The observed data at the end of 2007 was interpolated to be the initial Cl − concentrations. The Cl − concentrations in rainfall and lateral flow can be set as 2 mg/L (Gao et al. 2015).
The study area is divided into various parameter zones with different layers. The initial values of hydraulic conductivity (K), storage coefficient (S s ) specific yield (S y ), porosity (n), and dispersion coefficient (α) can be obtained from the previous work done by Zheng et al. (2015) in the study area.

Model calibration and validation
From December 2007 to December 2015, the model was calibrated by the trial-and-error method based on the groundwater table and the Cl − concentrations. Four typical observation wells (G01, G15, G22, and G30) distributed along the Chaobai River channel from upstream to downstream were selected to represent the groundwater table around the perennial and intermittent water-receiving river channel. The parameters were adjusted to minimize the difference between the observed and simulated outputs in the four typical observation wells. On account of the unconfined aquifer being most significantly affected by RW, the first layer is the focus of research. The validation period ranges from January 2016 to December 2017. The root mean squared error (RMSE) (Eq. (2)) and the coefficient of determination (R 2 ) (Eq. (3)) were used to evaluate the goodness-of-fit of the calibrated model.
where n is the number of observation points, O i is the i th observed data, S i is the i th simulated result, O m is the mean value of observed data, and S m is the mean value of simulated result.

Scenario design
Six modeling scenarios were designed based on different conditions to predict groundwater changes and analyze the influences of precipitation, groundwater extraction, and RW recharge (Table 2). A baseline scenario was designed to compare the modeling results of different scenarios. Then, the groundwater dynamics in the next 10 years after 2017 were predicted using the calibrated numerical model. Different annual and seasonal precipitation level changes would affect the groundwater system (Schoenheinz and Grischek 2011). However, it is difficult to predict the precipitation in the next 10 years based on the existing meteorological data. Therefore, it is necessary to assume multiple datasets reflecting different precipitation trends for prediction. Based on the monthly precipitation data of the Miyun Station from 1956 to 2020, 1993, 2005, and 1982 were selected to represent the dry, normal, and wet year, respectively. Their annual precipitation was 454.8 mm, 574.6 mm, and 710.9 mm, respectively. Accordingly, three precipitation scenarios (Scenario 1-3) were considered to repeat the monthly precipitation series of the three representative years. For groundwater exploitation, Beijing has taken a series of measures (SCPRC 2014). The extraction volume of groundwater in well field A has been halved since 2014. Therefore, two different groundwater exploitation scenarios were set up to discuss the impact of groundwater exploitation on the RW infiltration in the study area. Based on the extraction volume (2.7 × 10 7 m 3 /a) of well field A in 2017, the groundwater extraction volume of well field A in the next 10 years (2018-2027) was designed to decrease by 25% or increase by 25% in Scenario 4 and Scenario 5, respectively. Other well fields remain unchanged.
RW infiltration is an important recharge source for groundwater around the river channel. Scenario 6 was set up as no RW recharge in the next 10 years (2018)(2019)(2020)(2021)(2022)(2023)(2024)(2025)(2026)(2027). We compared the difference with or without RW recharge to discuss the impact of RW recharge on surrounding groundwater.

Estimation of Cl − and NO 3 -N loads in the unconfined aquifer
The pollutant load can be calculated from the groundwater storage and pollutant concentration for a given area, reflecting the changes in groundwater pollution under the combined effect of groundwater quantity and quality (Su et al. 2019;Kang et al. 2021). In this study, the concentration gradient was not considered. The formula is as follows: Table 2 Scenario design for the groundwater model in the study area a "2007-2017 series" means the monthly precipitation from 2018 to 2027 repeats the monthly precipitation series from 2008 to 2017. "1993 (dry year)" means the monthly precipitation from 2018 to 2027 repeats the monthly precipitation in 1993. And so on, for 2005 (normal year) and 1982 (wet year). b "2017" means the groundwater extraction volume in all well fields from 2018 to 2027 is the same as the groundwater extraction volume in 2017. "25% reduction (increase)" means the groundwater extraction volume of well field A from 2018 to 2027 is decreased (increased) by 25% compared with 2017. c "Yes" means RW recharge from 2018 to 2027 remains unchanged, the same as in 2017. "No" means there is no RW recharge from 2018 to 2027. where P refers to the Cl − or NO 3 -N loads (t), C i is the solute concentration of the groundwater in the i th cell (mg/L), and S i is groundwater storage in the i th cell (m 3 ).
The groundwater table and Cl − concentration of each cell can be obtained from the modeling results. The NO 3 -N concentration of each cell can be gained by interpolating the observed NO 3 -N concentrations in the discrete grids of the model.

Groundwater model calibration and validation
The fitted results for the simulated groundwater table and Cl − concentrations in four typical observation wells of the unconfined aquifer were presented in Fig. 2. During the calibration and validation periods, the RMSE value was 0.35 m, and the R 2 value was 0.99 for the groundwater table. The RMSE value was 11.9 mg/L, and the R 2 value was 0.68 for Cl − concentrations. It indicates that the model can be used to predict the variations of the groundwater table and simulate Cl − transport in the aquifer. The calibrated parameters of the layers referenced from Zheng et al. (2015) are shown in Table S1.

Changes in the groundwater table
The simulated results of the groundwater level at the end of 2008, 2009, 2012, and 2017 after RW restoring the river channel were shown in Fig. 3. The area with the lowest groundwater table is located in the north of the Xiangyang sluice. The groundwater tables were still falling during the 10 years (2007-2017) after RW recharge, along with groundwater exploitation. The groundwater table in the south of the Xiangyang sluice increased quickly, reaching around 3 m after the RW recharge for two years   [2007][2008][2009]. The closer to the river, the more the groundwater table rises due to the RW infiltration.

Groundwater budget analysis
The groundwater budget of the unconfined aquifer was summarized in Table 3. The total recharge during the simulation period (Dec. 2007-Dec. 2017) was less than the discharge, with a reduction in groundwater volume of 1.30 × 10 8 m 3 . The main recharge sources are precipitation infiltration (7.57 × 10 8 m 3 , 84.6%), lateral inflow (0.78 × 10 8 m 3 , 8.7%), and RW infiltration (0.59 × 10 8 m 3 , 6.7%). For discharge items, the leakage of the unconfined aquifer was 9.97 × 10 8 m 3 (97.4%), which was greater than the total recharge. The groundwater table of the unconfined aquifer was higher than that of the first and second confined aquifers. An increase in the vertical hydraulic gradient due to the deep groundwater pumping enhanced the RW infiltration into the deep aquifers. In contrast, the lateral outflow and evaporation were 2.3 × 10 7 m 3 (2.2%) and 4 × 10 6 m 3 (0.4%). They accounted for a small part of the total discharge volume. Therefore, the discharge leakage of the unconfined aquifer driven by the deep groundwater withdrawal was responsible for the negative balance of the unconfined groundwater storage.

Variations of the Cl − concentrations
The contour maps of the simulated Cl − concentrations at the end of the 1st, 2nd, 5th, and 10th years after RW restoring the river channel were shown in Fig. 4. The Cl − concentrations in the unconfined groundwater increased from 5-75 mg/L to 50-130 mg/L from the end of 2007 to 2010, and then remained stable (Fig. 2d). The groundwater Cl − concentrations in this area decrease with the increasing distance from the river channel. Both precipitation infiltration and subsurface lateral flow played a dilutive effect, the Cl − concentration in the areas away from the river channel was gradually decreased.

Variations of the unconfined groundwater storage
The unconfined groundwater storage was calculated by the simulated groundwater table each year at the end of December (Fig. 5). The storage decreased year by year from 4.83 × 10 8 m 3 at the end of 2007 to 3.88 × 10 8 m 3 at the end of 2011. Although the groundwater table around the river has risen after the RW recharged, the groundwater table in the north of Xiangyang sluice has been falling due to the continuous groundwater exploitation. Since 2012, the section between Henan rubber dam and Suzhuang rubber dam has been replenished by the RW. The storage increased to 4.1 × 10 8 m 3 by the end of 2012, then continued to decline and was up to 3.76 × 10 8 m 3 at the end of 2014 (Fig. 5).
In December 2014, the South-to-North Water Diversion Project began to transfer water to Beijing. The amount of groundwater extraction in the well field, which is located in the north of the Xiangyang sluice, began to be halved.  The declining trend of the confined groundwater level has slowed down. Moreover, the annual average leakage of the unconfined aquifer has decreased from 1.04 × 10 8 m 3 /a (2008-2013) to 8.81 × 10 7 m 3 /a (2014-2017). Due to the continuous RW infiltration and the groundwater extraction reduction, storage began to rise slowly since 2014. It rose from 3.76 × 10 8 m 3 at the end of 2014 to 3.85 × 10 8 m 3 at the end of 2017 (Fig. 5).

The zones of unconfined groundwater quality affected by the RW infiltration
The increase of Cl − concentrations reflects the impact of RW infiltration on groundwater. The zones where the Cl − concentration is higher than the initial concentration (5-75 mg/L) after RW recharge is defined as being affected by RW infiltration. We mapped the zones affected by RW for the unconfined aquifer accordingly (Fig. 6). The Cl − transport is controlled by both the groundwater flow rate and the concentration gradient. In the early stage of RW replenishment, the groundwater table elevated rapidly due to the RW infiltration, resulting in quick Cl − movement with the groundwater flow. Thus, the zones affected by RW are distributed on both sides of the river. Due to the unconfined groundwater flows from southwest to northeast, the affected zone of the left bank is more extensive than that at the right bank (Fig. 6). The affected zones gradually expanded with the passage of time and recharged channels (Table 4). But the annual increasing rate of the zones' area decreased year by year when the transferred volume of RW and the extent of the recharged channel were stable. It indicates that Cl − movement slowed down with the stable groundwater table around the river.

Cl − loads in the unconfined groundwater
The Cl − loads in the model area decreased annually, dropping from 1.7 × 10 4 t at the end of 2007 to 6.5 × 10 3 t at the end of 2017 (Fig. 5). In contrast, the Cl − loads in the affected zones rose from 1.8 × 10 3 t at the end of 2008 to 3.8 × 10 3 t at the end of 2017. There are rising trends both in the Cl − concentrations and Cl − loads of the unconfined groundwater in the affected zones. The zone can represent where the unconfined groundwater quality is affected by RW in terms of solute concentration and corresponding loads.

NO 3 -N loads in the affected zones
NO 3 -N is the typical contaminant in the RW (Pan et al. 2018;Li et al. 2019). Although NO 3 -N has been partially removed during the RW infiltration process, NO 3 -N from the RW still enters the groundwater ). The observed NO 3 -N concentrations of some observation wells occasionally exceed 10 mg/L (guideline for drinking water recommended by WHO (World Health Organization 2017). Furthermore, the change trends of groundwater NO 3 -N concentrations during monitoring periods were consistent with those of groundwater Cl − concentrations (Li 2020). It is one piece of evidence that RW influenced the concentration and distribution of NO 3 -N in the ambient groundwater.
Groundwater NO 3 -N loads in the affected zones were calculated at the end of each year (Fig. 5). It shows that the loads declined from 29.8 t at the end of 2008 to 11.9 t at the end of 2017. According to the results, RW infiltration did not increase the NO 3 -N loads in the unconfined groundwater. Furthermore, NO 3 -N attenuated in the affected zone, resulting in a reduction of the NO 3 -N loads.
Although the maximum NO 3 -N concentration in the RW is as high as 20.2 mg/L, the unconfined groundwater has  been subjected to denitrification during the RW infiltration, especially in the south of the Xiangyang sluice (Liu et al. 2018). Because denitrification hotspots are easily formed during surface water infiltration (Rivett et al. 2008;Ranalli and Macalady 2010;Trauth et al. 2018). The attenuation rate of NO 3 -N near the observation well G22 can reach 99.6% ). However, not all areas in the study area have a great denitrification environment. The NO 3 -N contents near the Xiangyang sluice are relatively high due to the thick gravel layer, which is not conducive to removing NO 3 -N (Xiong 2009). It resulted in a high level (48.8 t) of NO 3 -N loads at the end of 2011. The north channel of the Earth dam has been replenished only in May and October. Sediments can adsorb NH 4 -N during the wet period, which can be converted by nitrification reaction into NO 3 -N during the dry period (Böhlke et al. 2006). Therefore, the unconfined groundwater near the intermittent water-receiving channel appeared high NO 3 -N concentrations (~ 14.7 mg/L).

Effect of different recharge and discharge conditions on the groundwater system
The difference between the modeling results of groundwater Cl − concentrations under different precipitation scenarios is less than 10 mg/L. The prediction results show that Cl − concentrations in the unconfined groundwater could be rapidly diluted without RW recharge, as RW is the only source of high concentrations of Cl − in the model. There is almost no difference in the modeling results in Cl − concentration under different groundwater extraction scenarios. Therefore, only the influence of different scenarios on the groundwater table will be discussed below.
Precipitation infiltration as the main source of recharge for unconfined groundwater affects the variations of groundwater levels. In different precipitation scenarios, the groundwater table in the wet scenario > baseline scenario > normal scenario > dry scenario (Table 5). The interannual variation of the groundwater table under the baseline scenario is the smallest, approximated as a state of balance. After 10 years, the groundwater table difference between the dry, normal, and wet scenarios and the baseline scenario is up to − 0.9 m, − 1.3 m, and + 2.9 m, respectively. Long-term low rainfall would cause the groundwater table to drop, and greater precipitation would result in an increased groundwater table, especially in the higher permeability areas.
The scenarios with baseline precipitation and RW recharge (Scenarios 4 and 5) demonstrated the negative effect of groundwater overexploitation. The results show that the reduction in groundwater withdrawal can recover the groundwater table (Table 5). The closer to the well field, the more significant the groundwater table changes. The groundwater table at G01 increases by 1.2 m under the reduction scenario and decreases by 0.96 m under the increase scenario in 2027 compared with 2017. It suggests that the groundwater table recovery is not easy to achieve in either the high precipitation scenario or the scenario with a small reduction in groundwater extraction. Limiting groundwater extraction and taking other measures need to be considered.
The scenario of baseline precipitation and groundwater extraction (scenario 6) was designed to demonstrate the influence of RW recharge. Compared with the baseline scenario, the groundwater table dropped when there was no RW recharge and then reached a relatively stable state. The difference in groundwater table between the baseline scenario and the no RW recharge scenario is about 2.5 m, 10.0 m, 4.5 m, and 4.2 m at G01, G15, G22, and G30, respectively (Table 5). Scenario 6 has the most significant effect on groundwater level changes. It suggests that continuous RW recharge contributes to a faster and more significant increase and maintenance of the groundwater table than high precipitation and reduced groundwater extraction.

Environmental implications for management of a river channel recharged by RW
RW has been used to restore the Chaobai River, while the unconfined aquifer near the river channel was unintentionally recharged. It can be considered as a riverbank filtration (RBF) system, a way of anthropogenic aquifer recharge (AAR) (Todd 1959;Morel-Seytoux 1985;Dillon 2005;NRMMC et al. 2009;Maliva 2020a). As a result, the RW infiltration affects the groundwater level and quality around the river channel, especially the unconfined groundwater.
RW infiltration has quickly replenished the unconfined groundwater, resulting in the groundwater table around the river channel rising rapidly by 3-4 m in the first two years  (Zheng et al. 2015;He et al. 2021). However, increased recharge via RW infiltration may not necessarily result in a corresponding increase in the groundwater storage volume (Maliva 2020b). Increasing recharge in isolation may not solve the problem of groundwater over-extraction (Gale et al. 2006;Foster and Garduño, 2013). According to the findings of this study of unconfined groundwater storage, recovering and maintaining the groundwater table and storage requires a combination of increased RW recharge and reduction extractions.
Water quality issues are usually the most concern in RW utilization. The infiltration process, driven by natural filtration and groundwater pumping, improves or degrades the recharged water quality. It depends on the quality and hydrochemical characteristics of RW and native groundwater, and the geochemical processes occurring during RW passes through the riverbed and underlying aquifer (Ray 2008;Stuyfzand 2011;Tyagi et al. 2013;Hu et al. 2016). In the study area, the subsurface hydrochemistry was modified by mixing the infiltrated RW with native groundwater and potential water-rock interactions (Yu 2013;Liu et al. 2018). The hydraulic travel time of RW infiltration into the 30-m depth was about 6.5 months . The hydrochemical type of the unconfined groundwater has changed from HCO 3 -Ca·Mg into HCO 3 ·Cl-Na·Ca, which is the RW type Jiang et al. 2020). The average proportion of the RW in the unconfined aquifers was about 53% from 2007 to 2018 (He et al. 2021). NO 3 -N in the RW was well attenuated with an attenuation rate of 99.6% during infiltration in G22 . Additionally, some emerging contaminants have been detected in the unconfined groundwater from this area, such as endocrine-disrupting compounds (Li et al. 2013;Ma et al. 2015) and antibiotics and antibiotic resistance genes (Zhang et al. 2019).
Although the RW recharge to a river channel may inevitably affect the groundwater quality of the underlying aquifers, the affected zones are limited. In the affected zones, the Cl − loads of the unconfined groundwater were increased and the NO 3 -N loads of the unconfined groundwater attenuated in the zones. But outside of these affected zones, the groundwater quality is characterized by the initial concentration of the aquifer.
RW utilization for recharging the river channel and the ambient groundwater system is usually a double-edged sword, which can bring the environment both benefits and harms. The main benefits are restoring the riverine ecosystem, increasing and maintaining the groundwater table and storage, and reducing groundwater pumping. At the same time, the disadvantage lies in the potential deterioration of groundwater quality, which depends on the inputs, attenuation of contaminants, and water-rock interaction in river-aquifer systems. Where conditions are hydro-geologically favorable, RW recharge to the river channel can be a valuable way of alleviating water shortages (Maliva 2020a). It is important and necessary to investigate and evaluate whether the field is suitable for RW recharge (Alam et al. 2021). Based on the numerical modeling, this study shows that apart from the affected zones, the RW infiltration has lesser environmental impacts on ambient groundwater quality. In future work, we should pay more attention to variations of water quality after RW recharge, and take necessary measures to manage RW recharge and ambient groundwater. It is an essential aspect of ensuring the safe use of RW. Additionally, as the unconfined groundwater table rises, a threshold of water level should be set to prevent soil salinization in the RW receiving area.

Conclusions
Since the operation of the Beijing Municipal Water Diversion Project in 2007, it has been more than 10 years since the implementation of RW supply to the Chaobai River, resulting in some environmental effects (e.g. groundwater level rise and water quality change). This study established a three-dimensional groundwater flow and solute transport model based on long time-sequenced groundwater monitoring data to evaluate the impact of RW recharge on the groundwater table and storage change, as well as Cl − and NO 3 -N loads in the unconfined aquifer. Seven scenarios are designed to reveal the influence of different precipitation, groundwater extraction, and RW recharge conditions on the unconfined groundwater system.
The results show that the groundwater table around the river channel rose by about 3-4 m from Dec. 2007 to Dec. 2009 after RW recharge and then remained stable. Groundwater exploitation caused the declining groundwater table near the well field in the north of the Xiangyang sluice, and the unconfined groundwater was in a negative equilibrium state from Dec. 2007 to Dec. 2017. The unconfined groundwater storage increased with the continuous recharge of RW and the reduced groundwater exploitation. Reducing the amount of groundwater exploitation by 25% can increase the surrounding groundwater table by ~ 1 m in 10 years (2018)(2019)(2020)(2021)(2022)(2023)(2024)(2025)(2026)(2027). Stop supplying RW, and the unconfined groundwater table around the river channel will drop by about 5 m in 10 years (2018)(2019)(2020)(2021)(2022)(2023)(2024)(2025)(2026)(2027).
The range of Cl − concentrations in the unconfined groundwater significantly increased from 5-75 mg/L before RW recharge to 50-130 mg/L due to the RW infiltration. The affected zones can be identified by the increase of Cl − concentration. The RW-affected area increased from 11.7 km 2 in 2008 to 26.7 km 2 in 2017. The affected zones tend to stabilize year by year, with a steady RW supply. The Cl − loads of the unconfined groundwater in the affected zones increased annually, but the NO 3 -N Environmental Science and Pollution Research (2022) 29:36100-36114 36111 loads decreased as a whole due to the NO 3 -N attenuated in the zones. In general, the river channel recharged by RW in the study area is beneficial for increasing and maintaining the groundwater table and the groundwater storage with lesser environmental impacts. To a large extent, the pressure on water demand in Beijing has been relieved. Long-term monitoring of surface and groundwater quality around the river channel should be implemented to keep a watchful eye on changes in groundwater quality and related environmental impacts. Improving the RW quality and optimizing replenishment plans through modeling are recommended for water management around the river channel recharged by RW. To minimize the negative environmental effects, it is necessary to regulate the water quality and quantity of RW in the long run to realize the optimal utilization of RW.