Effects of downstream environmental flow release on enhancing the groundwater recharge and restoring the groundwater/surface-water connectivity in Yongding River, Beijing, China

The Yongding River (Beijing, China) was dry most times of the year, and groundwater storage was severely depleted. To address this issue, a river rehabilitation project was initiated. A downstream environmental flow release (EFR) project from upstream reservoirs has been implemented since 2019. This study evaluated the impact of EFR by constructing transient groundwater-flow and numerical tracer transport models to simulate the hydrogeological responses to the water release events in 2019–2020. The study identified two factors that significantly influence the river leakage rate, which are operational factors (i.e., water release rate and duration) and physical factors (i.e., hydraulic properties of the riverbed, regional hydraulic gradients, and groundwater depth) that determine the maximum water availability for groundwater recharge and maximum infiltration capacity, respectively. Predictive modelling was performed to assess the long-term effects of the proposed EFR scheme from 2021 to 2050, which showed that groundwater levels along the river will increase by 10–20 m by 2050. Groundwater storage is expected to be largely recovered and groundwater/surface-water connectivity in the middle reach of the river will be restored. This restoration will not only maintain the environmental flow for the benefit of ecosystems but also enhance groundwater recharge, promoting sustainable groundwater development in the region. Overall, this study provides valuable insights into the effectiveness of the proposed EFR scheme in achieving sustainable groundwater development in the region.


Introduction
Flow regimes in rivers globally have been altered by anthropogenic activities, including the regulation by large hydraulic infrastructures, changes in land use, and overabstraction of water (Harwood et al. 2017).Reductions in river flow result in the degradation of ecosystems globally (Bunn and Arthington 2002;Döll et al. 2009).In order to determine flows necessary to maintain riverine habitats and ecosystem services, the science of environmental flow assessment has been established (Tharme 2003).Internationally, environmental flows are defined as "the quantity, timing, and quality of freshwater flows and levels necessary to sustain aquatic ecosystems which, in turn, support human cultures, economies, sustainable livelihoods, and well-being" (Harwood et al. 2017).Rehabilitation projects to restore drying rivers have been implemented to maintain environmental flows, restore the aquatic ecosystem and biodiversity (Amoros et al. 2005;Rolls et al. 2012;Torabi Haghighi et al. 2018), reduce water supply deficit (Ramírez-Hernández et al. 2017;Salem et al. 2020;Tang et al. 2018), improve river water quality (Tang et al. 2018), and enhance groundwater recharge (Kennedy et al. 2017;Tosline et al. 2012).
Yongding River catchment, in North China, also faces problems of drying up of the main river courses and continuous decline of groundwater levels due to the alteration of its natural flow regime (Jiang et al. 2014).Before the 1950s, flood events occurred frequently and had catastrophic impacts on the agricultural activities in Beijing Plain; hence, to mitigate the flood risk, many reservoirs were constructed between the 1950s to 1970s.As a consequence, the river discharge was significantly reduced and the river channel downstream dried up.Water pollution was also exacerbated due to increasing discharge of industrial and domestic effluents.Large pits resulting from illicit sand mining were exposed in the downstream riverbed, which heavily altered the natural river morphology and destroyed riparian vegetation.
Moreover, the groundwater system has also been influenced by the alteration of the Yongding River flows.Yongding River was the main source of water supply to the city of Beijing during 1950Beijing during -1970.Ever since the drying up and quality degradation of the river, groundwater has been used as an alternative source for the water supply.As a result, the groundwater level in the Yongding River area has declined dramatically since 1980.The lowering of the groundwater level also disconnected the hydraulic connectivity between the groundwater and surface water, further aggravating the drying of the Yongding riverbed downstream (Zhang and Zhang 2017).
To resolve the aforementioned problems in the Yongding River, a river rehabilitation project has been initiated, which aims to restore the riverine ecosystem and maintain environmental flows (Shao et al. 2021).Since the completion of the Yellow River water diversion project in the upstream, regulated environmental flow releases (EFR) from the upstream reservoirs has been performed.The released water greatly increased the surface runoff in the Yongding River (Sun et al. 2021).
Most of the previous EFR-related studies focused mainly on surface hydrological processes.Groundwater responses to the EFR were neglected and lack detailed assessment; however, groundwater plays an important role in sustaining the streamflow, especially during dry periods of the year (Lu et al. 2021;McMahon and Nathan 2021).The disconnection of the hydraulic connectivity of the surface water and groundwater has affected the hyporheic zone (Hayashi and Rosenberry 2002); thus, an integrated investigation of surface water and groundwater connectivity is necessary to the design of the EFR operation scheme (Song et al. 2020).
Numerical simulation is an effective tool for analysis of the EFR operation scheme to quantify the surface and groundwater interaction.The MODFLOW program is a commonly used model for simulating groundwater flow systems, which provides packages to simulate the groundwater/ surface-water interaction at different levels (Harbaugh 2005;Hughes et al. 2012;Niswonger and Prudic 2005).The results of the groundwater model can provide quantitative information for designing the EFR operation scheme.
In this study, multiscale transient groundwater models were constructed focusing on the western suburb area of the Beijing Plain to assess the impact of EFR on the groundwater system and surface flows of the downstream Yongding River in the Beijing Plain area.The objectives of the study were to (1) calibrate a refined transient groundwater flow model using the observations of three water release events in 2019 to 2020; (2) estimate enhanced groundwater recharge through river leakage and recovery of groundwater levels and storage; (3) delineate sections of the river where the connectivity of groundwater and surface water is restored; (4) determine influencing factors on the river leakage; and (5) predict the long-term effect of the EFR operation on the groundwater sustainability in the area.

Ecological rehabilitation projects in the Yongding River catchment
Management of the Yongding River plays a vital role in controlling floods, securing the water supply, and preserving the riverine and riparian ecosystem of the city of Beijing (Jiang et al. 2014).Figure 1 shows a map of the Yongding River catchment.The river flows through five provinces including Inner Mongolia, Shanxi, Hebei, Beijing and Tianjin.The total catchment area is 47,000 km 2 (Fig. 1a).The length of the main stream is 369 km, with a 2.85‰ average slope of the rivered.In total, there are 267 reservoirs constructed along the river.The largest three reservoirs are Cetian, Youyi, and Guanting.The total capacity of these three reservoirs is about 3 billion m 3 , and their combined surface area is 43,402 km 2 .The annual average discharge above the Guanting Reservoir is around 2.08 billion m 3 .With a temperate continental climate, the average annual precipitation varies between 400 and 650 mm in different parts of the catchment.However, the temporal distribution of the rainfall and surface runoff is extremely uneven.Surface runoff during the flood season (July-August) accounts for about 60% of the total runoff throughout the year.The largest annual inflow to Guanting Reservoir was recorded as 3.06 billion m 3 in 1939 and the minimum was only 374 million m 3 in 1972.
Starting from the Guanting Reservoir, the total length of the river channel in Beijing is 189 km and can be characterized into three parts (Fig. 1b).The mountain part accounts for 108.5 km from the Guanting Reservoir to the Sanjiadian sluice, which meanders with a steep gradient.The upper plain part starts from the Sanjiadian sluice and continues to the Lugou Bridge sluice, with a length of 18.4 km.The river channel is wide and smooth with good permeability.The average width of the river channel is 300-500 m.The lower part of the river in the Beijing Plain is around 60 km in length, with an elevated riverbed and lower permeability.The upper part of the alluvial plain has relatively simple geology.The Quaternary deposits consist mainly of gravel and sand with high permeability.The average thickness of the unconfined aquifer varies from 30 to 50 m.This aquifer is underlain by bedrock characterized as a prevailing layer with significantly lower hydraulic conductivities.Figure 1d, e shows two cross sections indicating the hydrogeological conditions of the Beijing Plain and Yongding River.
In 2009, Beijing Water Authority launched the project "Yongding River Green Ecological Corridor", aiming to rehabilitate the riverine ecosystem and maintain the environmental flow of the river (BWA 2013).The water quality of the Guanting Reservoir was improved by constructing wetlands.Riparian vegetation destroyed by the sand mining in the Beijing Plain area was restored, which helped to prevent sandstorm hazards and soil erosion of the riverbed.Large sand mining pits were converted into four lakes: Mencheng Lake, Lianshi Lake, Xiaoyue Lake, and Wanping Lake (Fig. 1c).To prevent large river leakage, liners were placed on the middle part of the lake bottom, which help to maintain a minimal water level in the lake (Hu et al. 2019).Reclaimed water from wastewater treatment plants has been discharged to the Yongding river channel to maintain river flow since 2010.The estimated river leakage to the aquifer in 2010 was only 16,500 m 3 , which is negligible (Hu et al. 2018).Thus, most of the river channel remained dry, and only the lakes were filled with water.In 2017, the north route of the Yellow River Diversion Project was launched.Water from the Yellow River was diverted to several reservoirs along the Yongding River, which increased the water availability to implement the EFR (Liu et al. 2016).To better recover the riverine ecosystem in the plain part of the Yongding River, a number of wetlands were constructed including Yuanbo Lake (Fig. 1c).The total storage capacity of these lakes is ~21 million m 3 .
In 2019, the large-scale EFR from the upstream reservoirs to the Yongding River channel started on March 16th and continued for 93 days until June 16th.The total discharge from the Sanjiadian Sluice was 132 million m 3 (Ma et al. 2020).In 2020, there were two release events of shorter duration in spring and autumn.Before the flood season, 142 million m 3 of water were released during 32 days from April 21st (Zhao et al. 2020), while later on October 14th, 51 million m 3 of water were released in 22 days (Wu et al. 2021).It has been reported that with release events during 2019-2020, the riverine ecosystem in Yongding River has significantly improved (Shao et al. 2021).Besides, groundwater levels near the Yongding River channel recovered after the release events.Figure 2 shows the long-term historical groundwater level changes in an observation well near the Yongding River alluvial fan (location of the observation well Fig. 1c).Due to the continuous drought from 1999 to 2010 and groundwater overexploitation, the groundwater level decreased by more than 15 m.Groundwater levels started to recover slightly after 2015 because of the reduced groundwater abstraction when the transferred water from the south-north water diversion replaced partial groundwater abstraction in the area.The groundwater level has recovered significantly with the operation of the EFR since 2019.

Use of groundwater flow model to compute river leakage
To understand the response of the groundwater system to the EFR and the effectiveness of the enhanced groundwater recharge in restoring groundwater depletion, a multiscale 3D transient groundwater flow model was constructed and calibrated using MODFLOW-2005 (Harbaugh 2005).The refined local model for the Yongding River area was constructed based on the regional model of the Beijing Plain developed in the earlier study.The regional groundwater model was constructed to simulate the monthly groundwater level changes and compute the groundwater storage change during 1995-2018 in the Beijing Plain (Liu et al. 2022).The model only covers the Beijing Plain area (Fig. 1b).The west and north boundaries lie along the contact between the hard rocks of the mountain areas and the alluvial deposits of the plain area.The effects of hard rock aquifers in the mountain areas were simulated by setting a lateral inflow boundary along the boundary.The model consisted of nine model layers including five aquifers and four aquitards.Details of the regional model are provided in the electronic supplementary material ESM1.Table 1 summarizes the boundary conditions, parameter setting, sources and sinks of the regional model and specifies the package used in MODFLOW.
In this study, the regional model grid of 1,000 m × 1,000 m was refined with a local grid of 100 m × 100 m in Yongding River area (Fig. 3a). Figure 3b shows the refined local model grid and lakes and wetlands along the Yongding River.In total, 19 columns and 28 rows of the regional model were refined into a local model grid of 190 columns and 290 rows, accounting for 7.6% of the entire model domain.In the Yongding River alluvial fan, only a single layer of gravel was present (Fig. 1d, e) and was simulated as an unconfined aquifer with the top model layer.The model calibration period was set as three years from 2018 to 2020 when EFR was conducted.The stress period was set as daily during the EFR events and monthly for the rest of the time to save computation time.The starting head was taken from the computed groundwater heads from the regional model in January 2018.
River (RIV) package of MODFLOW was used to simulate the leakage from lakes and wetlands.Due to the channelization and lining of the river, most of the surface water leakage occurs at the constructed lakes and wetlands along the river channel.There are in total five lakes, two wetlands and a surface water reservoir, which have been conceptualized into 17 polygons based on their hydraulic properties (Fig. 3b).Data for the RIV package include river bottom elevation (H b ), the riverbed conductance (C RIV ) and river head stage (H River ), which are listed in Table 2.
The river stage of each polygon was initially estimated from the total released water reported at the two hydrological stations, and the infiltration area of each polygon was slightly adjusted based on the simulated results.River Table 1 General information on the regional groundwater flow model conductance was assigned an initial value based on the hydrogeological conditions and calibrated by matching the computed groundwater level with the observed levels from the 46 observation wells (Fig. 4) in the vicinity.Afterwards, the leakage rates for all lakes and wetlands were computed and the contribution of the enhanced groundwater recharge for the recovery of the aquifer system was analysed.

Tracking the movement of leakage water by MT3DMS
The risk of groundwater pollution from the EFR implementation should be addressed since the surface-water quality has a larger variation (Alam et al. 2021).Tracking of the infiltrated water and the mixing process provides information to delineate the area at risk of pollution.The movement of the infiltrated water can be tracked by a particle tracking technique.In this study, the MT3DMS program (Zheng and Wang 1999) was used to track the movement of leakage water from the river during the EFR events.A numerical tracer, with a concentration of 100 mg/L, was assigned in river polygons, while the background tracer concentration was set to 0 mg/L.Effective porosity was specified based on the soil type properties for each parameter zone.The advective transport was run to simulate tracer spreading along the river, using the Method of Characteristics (MOC) technique.
The computed tracer concentration delineates the area of spreading leakage water and the percentage of leakage water mixing with native groundwater.

Prediction of long-term effects on sustainable groundwater development
The EFR is implemented as a measure to maintain the environmental flow and sustain the riverine ecosystem of Yongding River.The positive side effect is enhanced groundwater recharge; thus, the long-term effects of the EFR on the groundwater system should be assessed.The main aim of the prediction model is to quantify the contribution of the river leakage from the Yongding River channel to the groundwater storage recovery, and to investigate the potential future surface-water/groundwater interaction due to the enhanced groundwater recharge.The prediction model uses a monthly stress period and covers the period from 2021 to 2050.The initial condition was set as the computed groundwater heads in December 2020 from the current model.Groundwater abstraction remains the same as that of 2020 for future periods.Additionally, the simulation also includes the EFR in the Yongding River, which are modelled by the RIV package.The release scheme in 2019 is used in the prediction models for all the years.Other hydrologic sources and sinks of the model also remain the same as in 2020 for future stress periods.The groundwater recharge from precipitation was derived from the projected monthly precipitation of the Representative Concentration Pathways (RCP) 4.5 climate scenario of the CCLM5-0-2 Regional Climate Model.Variance-based bias correction was applied to the monthly precipitation output of the CCLM5-0-2 model before conducting further study.Since the lateral inflow from the mountains to Beijing Plain is also dependent on precipitation, it will be influenced by future precipitation variation.Therefore, the prediction model also considers the variation of lateral inflow.Water balance results from previous models indicated that the recharge from the lateral inflow contributes ~10% of the total inflow and 20% of the precipitation infiltration.Accordingly, the input for future lateral inflow was calculated as 20% of the projected precipitation with monthly variations.The results of the prediction model offer insight into the effects of the EFR project on the groundwater system over a long-term period.

Response of the groundwater flow system to the EFR events
Figure 5 shows scatter plots of computed groundwater head versus observed head for the stress periods of the two EFR events in 2020-May and November.Most of the computed heads fit the observations well.The coefficients of determination of the two stress periods are 0.874 and 0.878, respectively.The average residual errors of the two stress periods are -0.42 and -0.71 m, while the root mean squared error (RMSE) values are 3.59 and 2.87 m, respectively.In general, the simulated groundwater heads are slightly overestimated at the upper reach and underestimated at the lower reach.However, the flow model was able to capture most of the peak groundwater levels in response to each EFR event.Figure 6 displays the time series of computed and observed groundwater heads in 12 observation wells (highlighted in red in Fig. 4) located in the upper, middle, and lower parts of the river.Three groundwater level peaks can be captured at most of the observation wells, which correspond to the three EFR events.The trends of the computed groundwater-level change are mostly in line with the observed trends, which verifies the reliability of the simulation results.
Figure 7 shows maps of the groundwater level contours before and after each EFR event.The groundwater flow direction follows mainly with the topographical gradient from the west suburb region towards the south and east.In general, the patterns of the contour are similar for the three EFR events.The groundwater levels beneath the river channel increased rapidly in response to river leakage during the EFR events; however, the magnitude and the area of groundwater-level rise are slightly different.
During the EFR event in spring 2019, two groundwater mounds were formed near Mencheng Lake (P3-P5) and Lianshi Lake (P7-P8) in the upper reach (Fig. 7b).The largest groundwater level rise reached 25 m (Fig. 8a).However, no significant groundwater level rise was found in the lower river reaches.After the rainy season, without receiving the river leakage from the EFR and with less natural groundwater recharge, the water mound under the riverbed dispersed and the regional groundwater level declined.During the EFR event in spring 2020 also, the groundwater level rose significantly (Fig. 8b); the maximum increase reached 14 m at Yuanbo Lake (P9-P10).Compared to the year 2019, groundwater level increase was found also in Wanping Lake (P15).Three groundwater mounds were formed at Mencheng lake (P3-P5), Yuanbo Lake (P9-P10) and Wanping Lake (P15), which gradually dispersed after the end of the EFR event.After the EFR event in the fall of 2020, the rise of groundwater level occurred only in the upper and middle reaches of the river with a smaller magnitude compared with the 2020 spring event (Fig. 8c).

Enhanced groundwater recharge and recovery of the groundwater storage
Table 3 summarises the annual groundwater budget of the local model area for 3 years.In the Yongding River area (523.6 km 2 ), natural recharge from precipitation and irrigation return flow contributes about one-third of the total inflow.Enhanced recharge from river leakage accounts for ~40-45%, and the rest comes from lateral inflow from the mountain boundary and other parts of the plain.Outflow components of the Yongding River area include groundwater abstraction, lateral flow to the downstream area of the Beijing Plain, evapotranspiration, and discharge to the river channel.Groundwater abstraction is still the major discharge, accounting for more than 75% of the total outflow in the region.The lateral flow to the other parts of the plain contributes to 20% of the total outflow; however, due to the trend of increasing groundwater levels, groundwater evapotranspiration and discharge of groundwater to the river channel increased in 2019 and 2020.
Compared to the situation in 2018, when there was no regulated EFR, the recharge through the riverbed leakage increased by 50 and 32% in 2019 and 2020, respectively.The recharge patterns of these three modelled years are also different.Figure 9 shows the change in infiltration rate through the riverbed with time.In 2018, most of the river leakage occurred between July and August from the stormwater released during heavy rain events.However, due to the EFR before the flood seasons in 2019 and 2020, the groundwater recharge maintained a relatively high rate during the EFR period.In 2019, 74.1% of the total river leakage (117.1 million m 3 ) came from the EFR period.In 2020 this percentage increased to 86.5% (120.4 million m 3 ).
The effect of the EFR on recovering the aquifer storage was significant.Figure 10

Tracking of the leakage water and the mixing process
The result of the solute transport simulation of a numerical tracer by the MT3DMS program reveals the movement of the infiltrated water and the mixing of the infiltrated water with the native phreatic groundwater.Before the EFR events, the background concentration of the numerical tracer was 0 mg/L (Fig. 11a), thus, with the 100 mg/L concentration of the infiltrated water, the simulated concentration represents  the mixing percentage of the infiltrated water and ambient groundwater.As shown in Fig. 11b, after the first EFR in 2019, the groundwater beneath the river channel at the upper stream reach was mostly replaced by the infiltrated water; however, with less infiltration at the lower stream reach, the mixing of the infiltrated and natural groundwater was also less.The infiltrated water spread to the aquifer driven by advective transport.The spreading area of the infiltrated water was limited to a 500-m band on each side of the Yongding River channel.Before the spring EFR in 2020 (Fig. 11c), the infiltrated water from the last year's EFR spread further into the aquifer.And with other source water (e.g., stormwater, precipitation infiltration) entering the aquifer, the percentage of infiltrated water from the EFR event reduced significantly under the riverbed.The same mixing process occurred repetitively during the two EFR events in the 2020 spring and fall, resulting in an alternately distributed high mixing and low mixing zone of the infiltrated water and the ambient groundwater at the upper reach and middle part of the EFR site (Fig. 11d-f).The lower part of the EFR site remained at a relatively low mixing percentage throughout the EFR period, due to the low leakage rate.According to the mass balance result of the MT3DMS model, at the end of the year 2020, only 0.43% of the infiltrated water left the aquifer through groundwater abstraction and 0.59% returned to the riverbed, which is negligible.In all, 99% of the infiltrated water stayed in the aquifer as storage without spreading far from the Yongding River channel.At the end of the year 2020, the extent of the >10% mixing zone was 44.3 km 2 .Table 4 summarizes the average travel times of the infiltrated water for arrival at different distances, which were computed as time lags between the breakthrough curves by cross-correlation.More details can be found in Figure S2.1 of ESM2.These travel times are indicative of the spreading of the infiltrated water to various distances from the river.The actual travel times from the river to these distances are longer since the travel time of infiltration through the vadose zone was not considered in the model.In general, the estimated groundwater travel times at the upper and middle reaches are similar.However, the estimated groundwater travel time is much longer in the lower reach.This spatial difference in the groundwater travel time can be explained by (1) the difference in the hydraulic characteristics of the aquifer, and (2) the difference in the spatial hydraulic gradient.With higher permeability and larger hydraulic gradient, the groundwater travel time is much faster at the upper and middle reaches, which also corresponds to the difference in mixing percentage.
At present, the mixing of infiltrated water and native groundwater occurs only in a limited area, as there is no large-scale groundwater abstraction in the vicinity.As a result, the mixing occurs only as a result of the natural hydraulic gradient; however, it is expected that with the prolonged operation of the EFR in the future, the mixing zone will expand to a larger area, and more infiltrated water will be captured by the groundwater abstraction.

Factors controlling the infiltration capacity of the riverbed
It is observed that during the three EFR events, the infiltration rate from river leakage also differs.Table 5 summarizes the estimated surface-water budget of the three EFR events according to the available information.By comparing the differences of each flow component, the controlling factors that influence the groundwater recharge efficiency can be identified.
Firstly, operational factors, including the water release duration and water release rate, determine the water availability of the groundwater recharge.Most of the released water has been recharged to the aquifer as groundwater storage; however, by comparing the two spring EFR events in 2019 and 2020, surface-water outflow to the downstream accounts for only 5% of the total amount of released water in 2019, while this percentage increased to 25% in 2020.Hence with a longer water release duration and smaller flow rate, outflow to the downstream was reduced; moreover, a longer release duration allows the released water to be detained longer in the recharge lakes, which prolongs the infiltration process.
Secondly, physical factors, including hydraulic properties of the riverbed and the aquifer, the regional hydraulic gradient and regional groundwater depth, control the maximum infiltration capacity.These factors also vary spatially in different river reaches.
Figure 12 shows the infiltration rate of nine recharge lakes located in the upper, middle, and lower reaches.The mean infiltration rate in the three EFR events was calculated as 0.10, 0.16, and 0.11 m/day, respectively.A high infiltration rate was found at the upper reach lakes (P1-P9).During the three EFR events, the infiltration rates were more than 0.5 m/day, and the maximum rate could be up to 0.8 m/ day.An intermediate infiltration rate of around 0.1-0.5 m/ day was found in the middle reach.And the infiltration rate decreased significantly along the river channel.At the lower reach, the infiltration rate dropped below 0.1 m/day.Only in Lake P14, a 0.2 m/day infiltration rate was found during the 2020 spring release event.The infiltration rate of the same lake also differs in each EFR event.Generally, the highest infiltration rate was detected during the 2020 spring EFR event in most of the lakes.A constant infiltration rate can be maintained in most of the lakes at the upper and middle reaches; however, in some lakes in the lower reach, e.g., P10, P15, and P17, the infiltration rate was decreasing during the EFR events in 2020 due to the rise of the ambient groundwater level.
As a typical alluvial fan-plain hydrogeological environment, hydraulic conductivities were higher at the alluvial fan and gradually decreased along the alluvial plain, which is the main reason for the differences in the infiltration rate in the upper, middle and lower reaches.
The regional hydraulic gradients along the river determine how fast the infiltrated water can be spread into the aquifer.As can be seen from the results of the mixing process in Fig. 11, the recharged water was not able to travel far due to the mild hydraulic gradient in the vicinity.
Besides, groundwater depth also determines the available underground space for infiltration.Figure 13 shows the groundwater depth before and after the EFR in 2019 and 2020.Compared with the groundwater depth at the lower reach, the groundwater depth is much deeper in the upper reach lakes.Surface water was hydraulically disconnected from the groundwater.Leakage through the riverbed was mainly driven by gravity, which allowed river leakage at a constant infiltration rate during each EFR event; however, due to the groundwater level increase, groundwater depth became much shallower before the EFR in 2020.Especially at the middle and lower reaches, groundwater levels near lakes P10, P15, P16 and P17 rose higher than the riverbed elevation, so the infiltration rate started dropping after the beginning of the 2020 spring EFR.Hence, it can be concluded that the permeability of the riverbed deposits, the groundwater depth and the hydraulic gradients near the streambed are three important physical factors that control the infiltration capacity during the EFR.

Long-term effect on groundwater sustainability
In the Yongding River area, most of the EFR water ultimately percolates to the underlying aquifer of the river channel, which is considered as "surface-water loss" from the perspective of surface-water management.For the purpose of restoring the natural river flow regime, surfacewater losses should be avoided as much as possible so that the environmental flow can be sustained; however, from an integrated water management viewpoint, managed surface-water infiltration is an extra source of groundwater recharge replenishing the shallow aquifer.Thus, from the perspective of sustainable groundwater management, this enhanced groundwater recharge from the EFR increases groundwater storage, which is also a strategic resource to achieve water supply security.The concept of this innovative option for water storage has been noticed in recent years but is still under investigation (Richter et al. 2012); therefore, when it comes to the design of EFR schemes, it is meaningful to consider both the restoration of the environmental flows and the response of the groundwater system.Previous studies show that only a very small portion of the total released water could arrive at the river downstream, and most of the water will contribute to groundwater recharge (Kennedy et al. 2017;Tosline et al. 2012).Long-term EFR operation is an extra source of groundwater recharge that will gradually replenish the depleted groundwater storage and support sustainable groundwater development.
The prediction model constructed in this study predicts the groundwater levels and water budget change up to 2050.In general, groundwater heads are predicted to increase by 10-20 m in different locations with distinct seasonal variations in the next 30 years.Long-term predicted groundwater level from 2021 to 2050 can be found in Figure S2.2 of ESM2.The high water level occurs during the EFR period annually and falls back during the dry season.The inter-annual variations can also be detected due to the variation of predicted groundwater recharge derived from the projected precipitation.Compared with the upper reach, the inter-annual fluctuation in the lower reach is relatively smaller.
With the implementation of the current EFR plan over the next 30 years, the estimated river leakage will vary from 99.9 to 112.6 million m 3 /year.The occurrences of low river leakage are mostly in wet years with larger natural groundwater recharge from precipitation.When the aquifer receives more natural recharge, the groundwater level near the river channel will rise beyond the river bottom elevation and the river leakage rate will be reduced; hence, with the groundwater level rise, the infiltration through river leakage will be less.The mean river leakage is predicted to be 106.1 million m 3 /year from 2020 to 2050, which is larger than the recharge from precipitation infiltration.Moreover, this extra source of groundwater recharge will mitigate the groundwater storage depletion accumulated over the last two decades.Groundwater storage change in 2020-2050 will fluctuate between 0.87 and -0.42 million m 3 under the influence of the precipitation infiltration variation.The groundwater balance will reach a new equilibrium state, which can accomplish healthy groundwater development in this region.A detailed groundwater balance chart can be found in Figure S2.3 of ESM2.

Restoration of groundwater/surface-water connectivity
Although it is predicted that the groundwater system in the Yongding River region will reach a new steady state in the next 30 years, the predicted result shows that the groundwater and surface-water system will remain disconnected.Groundwater levels will be higher than the riverbed bottom only in the middle reach so that the groundwater and river become connected hydraulically.Furthermore, the formation of a hyporheic zone beneath and adjacent to the streambed could provide soil water for the growth of riparian vegetation, which improves the groundwater-dependent ecosystem (Conant et al. 2019); however, the upper reach of the river channel will remain a permanently disconnected stream.
The restoration of the stream-aquifer connectivity can be achieved by several measures.Under the condition of adequate water availability from the upstream reservoirs, measures that increase the riverbed permeability, prolong the water release duration and increase the water release rate will enhance sufficient leakage and reconnect the surface water with the groundwater system.Thus, removing the liner that is currently installed at the river bottom to increase the riverbed permeability can be easily implemented compared with the other measures.The numerical model developed in this study can be used to evaluate the effect of the proposed measure by adapting the model settings accordingly.
This study tested a scenario assuming that the hydraulic conductance of all the lakes will be 1.5 higher than the current situation after removing the geomembrane liner at the river bottom, while the wetland remains the same.A small increase of the river conductance was used in consideration of possible clogging of the lake bottom during long-term operation.The configuration of the river stage and river bottom using the RIV package and other sources and sinks remains unchanged in the prediction model.By mapping out the area with higher groundwater level than the river bottom, the area with a connected stream can be delineated (details can be found in Figure S2.4 of ESM2).The results show the current extent of the connected stream, which only accounts for 22.2% of the total river channel and is mainly at the middle reach; however, with higher hydraulic conductance, the extent of the connected stream increases to 34% of the total river channel.
Moreover, during long-term operations, it is important to remain aware of the risk of clogging caused by the formation of a bioactive layer at the river bottom.The accumulation of this layer during low flow periods can impede the effectiveness of environmental flow release, by reducing water flow, and it can pose potential risks to groundwater quality.To address this issue, regular monitoring and maintenance of the river channel may be necessary to prevent the build-up of a bioactive layer.

Conclusions
To maintain the environmental flow and improve the riverine ecosystem, EFR has been implemented in the Yongding River area in Beijing, China.Released water from the upstream reservoirs replenishes river flow in the Beijing Plain.Transient groundwater flow and numerical tracer transport models were used to assess the observed EFR events in 2019-2020 and predict the effects of long-term operation of the EFR in 2021-2050.The results show that during the three EFR events from 2019 to 2020, groundwater levels increased significantly near the Yongding River channel.Groundwater storage increased by 97 million m 3 at the end of 2020.It was found that the river leakage rate was controlled by two main factors.Operational factors control the water availability for the groundwater recharge and physical factors determine the maximum infiltration capacity.With long-term operation of the proposed EFR scheme from 2021-2050, groundwater levels along the Yongding River are predicted to rise by 10-20 m.The enhanced river leakage will contribute significantly to groundwater sustainability in the region.The groundwater system will be able to reach a new equilibrium state in this region.The middle reach of the river channel will be reconnected with the groundwater hydraulically, which benefits the recovery of groundwater-dependent ecosystems.

Fig. 1
Fig. 1 a-c Maps of the Yongding River location and catchment at different scales and d-e two cross sections

Fig. 2
Fig. 2 Groundwater levelelevation, m above sea level (asl)-from January (01) 1995 to 2021 in the Yongding River alluvial fan Fig. 3 a Spatial discretization of the regional model and the local refinement.b Spatial discretization and the conceptualization of the lakes and wetland parks in Yongding River

Fig. 4
Fig. 4 Location of the observations near the Yongding River groundwater EFR site (Beijing Water Authority 2021)

Fig. 5
Fig. 5 Scatter plots of the computed versus observed groundwater heads (elevation, m asl) from the groundwater flow model results in a May and b November 2020

Fig. 7
Fig. 7 Groundwater level (elevation, m asl) contour maps before and after three EFR events in Yongding River: a the 2019 EFR, b the 2020 spring EFR and c the 2020 fall EFR shows the groundwater storage change from 2018 to 2020 in the Yongding River region.In 2018, groundwater storage was depleted during the dry season.Groundwater storage only increased during July and September; however, positive storage change occurred from March until October 2019 and most of the months in 2020.The cumulative storage change reached 97 million m 3 , which demonstrates the effectiveness of the EFR in recovering groundwater storage.

Fig. 10
Fig. 10 The monthly change of groundwater storage and the cumulative storage change in the Yongding River region

Fig. 12
Fig. 12 The infiltration rate of nine lakes in the a upper, b middle and c lower reaches from January (01) 2018-2020

Table 2
Data

Table 3
Annual groundwater water balance in Yongding River area

Table 4
Estimated groundwater travel time based on the cross-correlation analysis results