Development of a groundwater flow and reactive solute transport model in the Yongding River alluvial fan, China

Abstract

The Yongding River in the western suburbs of Beijing has been recharged with reclaimed water since 2010 for the purpose of ecological restoration. Where the reclaimed water is not well treated, it poses a danger to the aquifer underneath the river. To provide a reliable tool which could be used in future research to quantify the influence of reclaimed water in the Yongding River on the local groundwater environment, a transient groundwater flow and reactive solute transport model was developed using FEFLOW^TM in the middle-upper part of the Yongding River Alluvium Fan. The numerical model was calibrated against the observed groundwater levels and the concentrations of typical solutes from June 2009 to May 2010 and validated from June 2010 to December 2010. The average RMSE and R² of groundwater level at four observation wells are 0.48 m and 0.61, respectively. The reasonable agreement between observed and simulated results demonstrates that the developed model is reliable and capable of predicting the behavior of groundwater flow and typical contaminant transport with reactions. Water budget analysis indicates that the water storage in this aquifer had decreased by 43.76×10⁶ m³ from June 2009 to December 2010. The concentration distributions of typical solutes suggest that the middle and southern parts of the unconfined aquifer have been polluted by previous discharge of industrial and domestic sewage. The results underscore the necessity of predicting the groundwater response to reclaimed water being discharged into the Yongding River. The study established a coupled groundwater flow and reactive solute transport model in the middle-upper part of the Yongding River Alluvium Fan, one of the drinking water supply sites in Beijing city. The model would be used for risk assessment when reclaimed water was recharged into Yongding River.

Keywords

Yongding River alluvium fan / groundwater flow / reactive solute transport / FEFLOW / ecological restoration

Introduction

The Yongding River is the largest river flowing through the Beijing municipality, commonly known as the “mother river” of Beijing. The Yongding River started to dry up in 1980 and has been permanently dry since 2010 because of overexploitation in upstream areas (^{Peng et al., 2013; Jiang et al., 2014}). In order to improve the environmental quality and the socio-economic situation along the river, a project named “the Yongding River Green Ecological Corridor” was launched by the Beijing local government in 2009 (^{BWA, 2009a}). The proposed creation of lakes and wetlands were linked to the river which formed the Yongding River ecological system. Two wetlands were constructed near the Guanting reservoir, which is located in the canyon section of the Yongding River. Four lakes along the urban section of the Yongding River have been constructed and recharged with reclaimed water. Seepage-proof liners on the riverbed help to reduce river leakage and maintain the river water storage (Fig. 1). Reclaimed water commonly contains various kinds of solutes and may pose a risk to groundwater (through river leakage) if the sewage treatment does not meet the required standard, especially either over the longer term or at damaged sections of the bed liner. Thus, it is necessary to understand the mechanisms of groundwater flow and the fate of contaminants in the aquifer beneath the river, so as to provide early warnings on possible groundwater contamination disasters.

Numerical simulation is a useful tool to represent spatial and temporal variability of groundwater flow and transport. Several numerical models have been successfully constructed for simulating groundwater flow and to quantify the impact of human activities on groundwater in the Yongding River Alluvial Fan (YRAF) and areas of similar climatic conditions. For example, a transient groundwater flow model was developed to evaluate the applicability of artificial recharge of groundwater in YRAF using MODFLOW (^{Hao et al., 2014}). ^{Sun et al. (2011)} established a transient groundwater flow model for the Nankou area in the northwestern edge of the Beijing plain indicating unsustainable groundwater exploitation from 1998 to 2007. ^{Cao et al. (2013)} constructed a groundwater flow model for the entire North China Plain (NCP) to assess several management strategies for groundwater sustainability in the NCP. ^{Wang et al. (2008)} analyzed the water budget in NCP by MODFLOW integrated with MAPGIS. All the above mentioned studies evaluated the dynamic nature of groundwater systems and underscored the groundwater drawdown of these areas over the past few decades.

The section of the Yongding River flowing through the western suburb of Beijing is the middle-upper part of YRAF. These areas show high risk of groundwater pollution as the sediments of both vadose zones and aquifer systems are coarse-grained and thus susceptible to the contamination posed by human activities (^{Wang et al., 2012}; ^{Guo et al., 2014}). The YRAF has suffered from severe groundwater pollution according to the contamination assessment of the shallow aquifer in the Beijing plain (^{Guo et al., 2014}). Nitrogen (N) is not only a common type of contaminant in shallow aquifers on a world wide scale (^{Stuart and Lapworth, 2016; Collins et al., 2017}; ^{Soldatova et al., 2017}), but also a major water quality indicator in the Yongding River replenished by reclaimed water (^{Yu et al., 2017}). The hydrological processes, such as runoff, leaching, and groundwater recharge can introduce various contaminants generated by human activities in groundwater. To address the effects of different contaminant sources on groundwater systems, several solute transport models have been applied to several vicinities of the YRAF. Pollution ranges near a garbage dump field in the northern suburb of Beijing was determined by a simple one-dimensional groundwater solute transport model (^{Liu et al., 2005}). A three-dimensional (3D) transport model with streams and waste-irrigated fields as groundwater pollution sources has been set up to predict contamination transport pathways over time in the northwest of the Beijing plain (^{Sun et al., 2012}). Environment tracers were also examined to identify different sources of nitrate (NO₃⁻) from different land use types (^{Liu et al., 2014}; ^{He et al., 2016}; ^{Wang et al., 2017}). The above results showed that irrigation drainage and sewage were major inputs of NO₃⁻ in farmland and residential areas, respectively. To mitigate the likely impacts of human activities on groundwater quality, a sound understanding of not only the sources, but also the transport and transformation of contaminants in groundwater is necessary (^{Collins et al., 2017}). In situ monitoring results indicated that reclaimed water infiltration did not cause N contamination of groundwater in the Daxing District of Beijing, downstream of the YRAF, because of effective attenuation in a sand-clay alternating lithologic structure (^{Yin et al. 2016}). Reactive transport of three forms of N, i.e., NO₃⁻, nitrite (NO₂⁻), and ammonium (NH₄⁺) in the urban area of Beijing was simulated by a 3D groundwater model MODFLOW from 2004 to 2008, suggesting the decreasing attenuation potential of N along the groundwater flow direction (^{Dong, 2010}). The description of reactive transport of N in the middle and upper YRAF has not yet been reported.

As a first step of a comprehensive study, i.e., the impact of reclaimed water used as ecological water on the local groundwater environment, the main purpose of this study is to develop a numerical groundwater flow model combined with a reactive solute transport model of typical pollutants for the middle-upper part of YRAF. The flow and reactive transport model was calibrated under transient conditions and used to quantify the flow budget, flow pattern, and the concentration distributions of chloride (Cl⁻) and NO₃⁻. Both the calibration and validation used data for periods that preceded the project implementation. Water budget analysis was also conducted using the established model. The impact of river leakage on the local groundwater environment after the project construction will be predicted in a companion paper based on the model developed in the current study.

Materials and methods

Study site

The study site is located in the western suburb of the Beijing plain, a typically semi-arid region (Fig. 2). The total study area is 365 km². It is highly urbanized, where the dominant land use is residential with a small proportion of dry lands, croplands, forest lands, and grasslands (^{Zhu et al., 2013; Liu et al., 2014}; ^{Yu et al., 2015}). The focused urban section of the Yongding River is from Sanjiadian station (upstream) to Efang village (downstream). The study area includes four connected lakes, i.e., Mencheng Lake, Lianshi Lake, Xiaoyue Lake, and Wanping Lake along the Yongding River, with a total channel length of 37 km and channel area of 6.8 km².

The climate of the study area is characterized by a continental monsoon climate, with an average annual precipitation of 563 mm (1961–2010). 80% of the precipitation falls in the period of June to September. The annual average water surface evaporation is 1182 mm based on the records from 1999 to 2009 (^{Wang et al., 2010; Zhang et al., 2013}).

The study area is bounded on the west by the Western Hills of Beijing, receiving lateral inflows. The study area extends eastward to the city center. The aquifers mainly consist of Quaternary alluvial deposits with relatively coarse grains. The hydrogeological characteristics of the study area are shown in Fig. 3(a). According to the hydrogeological profile (Fig. 3(b)), the major part of the studied aquifer is a single layer of sand and gravel with hydraulic conductivity (K) higher than 200 m/d, and to the east, thinner layers of pebbles and gravels with K values ranging from 100 to 200 m/d. These small pebble and gravel layers are separated by one to two thin layers of clay. Besides, the groundwater table dropped down below the first layer of the aquitard in the study area after 2010 (Fig. 3(b)). The unconfined aquifer is thus generalized as a single layer composed of pebbles and gravels, which accepts precipitation recharge directly. The base of the Quaternary deposits is formed by impermeable sediments.

Groundwater within the alluvium is an important source for industry and domestic water supply. No. 3 and No. 4 pumping well groups are located in the middle part of study area as shown in Fig. 2. The Shijingshan pumping well group located in the northern part of the study area is the second largest groundwater supply field in Beijing. There are also some privately-owned wells located in the study area. The western suburb area of Beijing not only provides water for the increased local demand, but also for the water demand of the central districts, which had inevitably caused over-exploitation of local groundwater. In the middle and upper part of YRAF, the groundwater levels have dropped by 10–25 m from 1975 to 2010. A groundwater depression funnel has appeared in the northeast region of the study area. As of 2013, the shallow aquifers in western areas such as in the Fengtai district have almost dried up (^{Yu et al., 2015}). There was a rising trend of groundwater level in the short term, such as the years 1995‒1996, with artificial recharge activities and in 1998 with flooding (^{Liu, 2012}).

Groundwater modelling

A quasi-3D transient flow and reactive transport model was established using FEFLOW^TM.

Governing equations

The groundwater flow in the unconfined aquifer is controlled by the following equation (^{Bear, 1977}; ^{Feng et al., 2011}),

(1)\{\begin{array}{l}\mu {\displaystyle \frac{\partial H}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left(K_{x}h{\displaystyle \frac{\partial H}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(K_{y}h{\displaystyle \frac{\partial H}{\partial y}}\right)+\epsilon x,y\in \mathrm{\Omega },t\ge \mathrm{0}\\ H(x,y,t){|}_{t=0}={{\displaystyle H}}_{\mathrm{0}}x,y\in \mathrm{\Omega },t=\mathrm{0}\\ H(x,y,t){|}_{{\mathit{\Gamma }}_1}=H(x,y,t)x,y\in {{\displaystyle \Gamma }}_{\mathrm{1}},t\ge \mathrm{0}\\ {\displaystyle \frac{\partial H}{\partial n}}{|}_{{\mathit{\Gamma }}_2}=\mathrm{0}x,y\in {{\displaystyle \Gamma }}_{\mathrm{2}},t\ge \mathrm{0}\end{array},

where µ is the gravitational specific yield of the unconfined aquifer (-); H is the groundwater level (m); h is the distance between the bottom elevation of the phreatic aquifer and the phreatic free surface (m); K_x and K_y are the hydraulic conductivities along the x and y coordinate axes (m/d); e is a volumetric flux per volume representing sources and/or sinks of water (m/d); t is the time (d); H₀ is the initial water level (m); W is the seepage domain; G₁ is the first type boundary of the seepage domain; G₂ is the second type boundary of the seepage domain, and n is the direction of the outer normal line of the second boundary.

The reactive solute transport is described as follows,

(2) \{\begin{array}{l}{\displaystyle \frac{\partial C^j}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left(D_x{\displaystyle \frac{\partial C^j}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(D_y{\displaystyle \frac{\partial C^j}{\partial y}}\right)-{\displaystyle \frac{\partial v_x{C}^j}{\partial x}}-{\displaystyle \frac{\partial v_y{C}^j}{\partial y}}+r^{j}x,y\in \mathrm{\Omega ,}t\ge \mathrm{0}\\ C^j(x,y,t){|}_{t=0}=C_{\mathrm{0}}^{j}x,y\in \mathrm{\Omega ,}t=\mathrm{0}\\ C^j(x,y,t){|}_{{\mathit{\Gamma }}_1}=C^j(x,y,t)x,y\in {\Gamma }_{\mathrm{1}},t\ge \mathrm{0}\\ {\displaystyle \frac{\partial C^j}{\partial n}{|}_{{\mathit{\Gamma }}_2}}=\mathrm{0}x,y\in {\Gamma }_{\mathrm{2}},t\ge \mathrm{0}\end{array},

where C^jdenotes the concentration of species j in the aquifer (mg/L); D_x and D_y are the hydrodynamic dispersion coefficients along the x and y coordinate axes (m/d); v_x and v_y are the velocities along the x and y coordinate axes; r^jis the sources and/or sinks of solute j (mg/L/d), and C_{\mathrm{0}}^jis the initial concentration of species j (mg/L).

The hydrochemistry and groundwater contamination change gradually along the flow paths in the YRAF. The groundwater hydrochemistry evolved from HCO₃-Ca-Mg to HCO₃-Na, with complex types of HCO₃-Cl-Ca-Mg and HCO₃-SO₄-Ca-Mg-Na in the middle part of the YRAF (^{Lin, 2004}; ^{Zhai et al., 2013}). The typical contaminants include NO₃⁻, Cl⁻, Mg²⁺, and Ca²⁺ (Total Hardness, TH), and the other indicators (^{Lin, 2004}; ^{He et al., 2016}). Both conservative solute Cl⁻ and reactive solute N were considered in the reactive solute transport model. Three forms of N can be involved in different N transformations in groundwater. According to ^{Dong (2010)}, adsorption of NH₄⁺ and nitrification (transformation of NH₄⁺ to NO₃⁻) which are two major attenuation mechanism of NH₄⁺ can occur in the unconfined aquifer of the Beijing Plain. The two processes were expressed by Eqs. (3) and (4). Denitrification is the main attenuation mechanism of NO₃⁻ in groundwater. The evidence of denitrification also has been found in the study area (^{Cheng, 2007; Dong, 2010}). This process was expressed by Eq. (5). NO₂⁻ is an intermediate species in both denitrification and nitrification. As NO₂⁻ is labile and most NO₂⁻ concentrations in groundwater samples in 2010 were relatively low (<0.01 mg/L), NO₂⁻ is not included as a transform component in the simulation. Nitrification and denitrification are biological reactions and commonly simulated by the Monod model. In our study area, substrate concentrations (NH₄⁺ and NO₃⁻) were much smaller than the published Monod constant so the Monod model was simplified to first order kinetics (Sheibley et al., 2003; ^{Cui et al., 2014}). According to the soil column experiments we conducted in 2011 (Hu et al., 2015), the adsorption of NH₄⁺ by riverbed soils of the Yongding River satisfied the first-order kinetic, so the process was also included. The source and sink terms of NH₄⁺ to NO₃⁻ are shown in Eqs. (6) and (7),

(3){\mathrm{NH}}_{\mathrm{4}}^{+}\to {({\mathrm{NH}}_{\mathrm{4}}^{+})}_s,

(4){\mathrm{NH}}_{\mathrm{4}}^{+}\to {\mathrm{NO}}_{\mathrm{3}}^{-},

(5){\mathrm{NO}}_{\mathrm{3}}^{-}\to {\mathrm{N}}_{\mathrm{2}},

(6)r^{{{\displaystyle \mathrm{N}\mathrm{H}}}_{\mathrm{4}}^{+}}={\displaystyle \frac{\mathrm{d}{C}^{{{\displaystyle \mathrm{N}\mathrm{H}}}_{\mathrm{4}}^{+}}}{\mathrm{d}t}}=-{{\displaystyle k}}_{\mathrm{1}}{C}^{{{\displaystyle \mathrm{N}\mathrm{H}}}_{\mathrm{4}}^{+}}-{{\displaystyle k}}_{\mathrm{2}}{C}^{{{\displaystyle \mathrm{N}\mathrm{H}}}_{\mathrm{4}}^{+}},

(7)r^{{{\displaystyle \mathrm{N}\mathrm{O}}}_{\mathrm{3}}^{-}}={\displaystyle \frac{\mathrm{d}{C}^{{{\displaystyle \mathrm{N}\mathrm{O}}}_{\mathrm{3}}^{-}}}{\mathrm{d}t}}={{\displaystyle k}}_{\mathrm{2}}{C}^{{{\displaystyle \mathrm{N}\mathrm{H}}}_{\mathrm{4}}^{+}}-{{\displaystyle k}}_{\mathrm{3}}{C}^{{{\displaystyle \mathrm{N}\mathrm{O}}}_{\mathrm{3}}^{-}},

where the subscript “s” means the species in the adsorbed state, k₁ is the first-order adsorption rate constant of NH₄⁺ (s⁻¹), k₂ is the first-order nitrification rate constant (s⁻¹), and k₃ is the first-order denitrification rate constant (s⁻¹) of NO₃⁻.

Groundwater sources and sinks

The flow region is recharged by rainfall infiltration, river leakage, seepage of pipelines, and irrigation drainage. Of these, rainfall infiltration contributes the highest, which can be calculated in each rainfall event as the product of the equivalent thickness of rainfall and spatially varying recharge coefficients. The recharge rate had a decreasing trend along the groundwater flow direction from 0.4–1.44 m/yr in the northwest YRAF to 0.07–0.48 m/yr in the southeast (^{Zhai et al., 2013}). The recharge coefficients were regionalized according to the lithology of subareas (Fig. 4(a)). The relatively high recharge coefficient implies high permeability of the subarea. The values were assigned based on the study of ^{Xu (2006)}. Lateral inflow from the Western Hills and the northern boundary is also an important source of the groundwater system.

Each well field within the simulated area was regarded as a point sink with seasonally-varying pumping rates. The pumping rates were obtained from the annual total consumption from each main drinking water company. The phreatic water evaporation is negligible considering that the groundwater depth of the site is normally larger than 20 m (^{Zhai et al., 2013; Hao et al., 2014}).

The pollutants mainly come from domestic and industrial sewage, which leach from sewage pipes and landfill sites into the groundwater. Some of the sewage was directly discharged into the seepage wells, seepage pits, or abandoned wells (^{Lin, 2004; Sun et al., 2011}). The concentrations of Cl⁻, NH₄⁺, and NO₃⁻ from the above-mentioned point sources were assigned according to the reported values from literatures, 200 mg/L, 2 mg/L, and 30 mg/L (^{Dong, 2010}; ^{Sun et al., 2011}).

Parameter determination

Hydraulic conductivity, specific yield, and effective porosity of the aquifer were calibrated based on documented values of these parameters. The values of hydraulic conductivity were assigned based on the study of Xu (2006), with small adjustments by calibrations (Fig. 4(b)). The partition and value of specific yield are similar to that set by ^{Hao et al. (2014)}, as shown in Fig. 4(c). Longitudinal dispersivity of the study area was set as 5 m according to tracer tests conducted by the Beijing Institute of Hydrogeology and Engineering Geology (Lin et al., 2011). Transverse dispersivity has been taken as one-tenth of the longitudinal dispersivity, 0.5 m.

The first-order adsorption constant of NH₄⁺ was assigned as 1×10⁻⁹ s⁻¹. The values of nitrification and denitrification rate constants were obtained from the literature and set as 5.8×10⁻⁹ s⁻¹ and 0.5×10⁻⁹ s⁻¹, respectively (^{Dong, 2010}).

Boundary conditions

The west and northwest boundaries accept lateral inflow from the mountains and were defined as spatially varied Dirichlet boundaries (20–66 m) according to the survey of water level distribution in June 2009 (Fig. 3(a)). The Dirichlet boundary condition is flexible compared with the Neumann boundary condition for further scenario predictions. Considering the western boundary is closed to the Yongding River, the boundary could be influenced when simulating the scenario of reclaimed water being discharged into the river. The eastern and southern boundaries were assigned as no flow boundaries, as these boundaries nearly overlap the groundwater flow lines.

The typical pollutants can be introduced from lateral infiltration of groundwater. Concentrations of Cl⁻ and NO₃⁻ along the western and northern boundaries were assumed temporally constant during the simulation period. The western boundary of our model domain is the boundary between the mountain areas and plains, the inflow concentrations are relatively stable compared with the concentrations at boundaries of the plain areas. The northern boundary receives mass input largely caused by human activities. Considering the typical contaminants show a relatively stable level in recent years (^{Lin, 2004}), the boundary concentrations were thus set to be invariant with time. The other is zero-flux concentration boundaries. The Cl⁻ concentration varied spatially, ranging between 5 and 150 mg/L according to the initial concentration distribution, while NO₃⁻ concentration varied in a range of 0–20 mg/L.

Initial conditions

The spatial distributions of initial groundwater head and concentrations of Cl⁻, NH₄⁺, and NO₃⁻ were interpolated by ordinary Kriging based on measured values in June 2009 as shown in Fig. 5. Relatively high concentrations of Cl⁻ were shown in the middle part of the study area and high concentrations of NO₃⁻ were found in the southeastern part. The initial concentration of NH₄⁺ is quite low in the unconfined aquifer, as most of the NH₄⁺ was adsorbed in the vadose zone before reaching the groundwater (^{Dong, 2010}).

Spatial distribution of the studied aquifer based on DEM

The elevations of the ground surface and the bottom of unconfined aquifer in 2010 were interpolated by Kriging based on the digital elevation (DEM) data. The spatial distribution of the unconfined aquifer is shown in Fig. 6. The northern part of the aquifer, where the Shijingshan water plant is located, is much thicker than the other parts of the study area, suggesting a larger capacity for groundwater supply.

Numerical simulation scheme

Based on the observed data of groundwater levels, the calibration period was selected as June 2009 to May 2010 and the validation period from June 2010 to December 2010. Restricted by the limited data available on solute concentrations in the aquifer, only the concentrations of Cl⁻ and NO₃⁻ in October 2010 were used for calibration of the solute transport model.

During the simulation, the time step was automatically adjusted according to the convergence history in the previous time step. To minimize numerical oscillation, different numerical methods for solving the model were tried. A shock-capturing method was proved to have a good feasibility and applied in the model. The modeling area was discretized using triangular elements, with 26,248 finite elements and 26,700 nodes. A manual trial-and-error calibration of the groundwater model provided an adequate fit between observed and simulated groundwater levels and reproduced the spatial and temporal groundwater head variations. The same calibration method was also applied for fitting concentrations of Cl⁻ and NO₃⁻.

Groundwater model performance criteria

The simulated values were compared with the measured values using the root mean square error (RMSE) and the coefficient of determination (R²). These metrics are given by

(8) \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \sum _{i=\mathrm{1}}^N{(O_i-S_i)}^{\mathrm{2}}}},

where N is the number of observations, O_i are the observed values, S_i are the simulated values. A smaller RMSE indicates better model accuracy. The coefficient of determination R² is

(9) R^{\mathrm{2}}={\displaystyle \frac{{\left({\displaystyle \sum _{i=\mathrm{1}}^N(O_i-O_a)(S_i-S_a)}\right)}^{\mathrm{2}}}{{\displaystyle \sum _{i=\mathrm{1}}^N{(O_i-O_a)}^{\mathrm{2}}}{\displaystyle \sum _{i=\mathrm{1}}^N{(S_i-S_a)}^{\mathrm{2}}}}},

where O_a is the average value of the observed data; S_a is the average value of the simulated data. The closer R² is to 1, the more reliable the model.

Results

Calibration and validation of groundwater level

The model was calibrated and validated by three types of data: (i) monthly averaged groundwater heads at 4 observation wells from June 2009 to December 2010; (ii) the contours of observed groundwater heads; and (iii) observed Cl⁻ and NO₃⁻ concentrations. It can be seen from Fig.7(a) that the simulated groundwater levels from 4 observation wells reasonably match the measured values. The correlation coefficient R² is above 0.48, which indicates good agreement between the observed and the simulated groundwater levels (Table 1). The averaged RMSE and R² of groundwater heads at the 4 wells for the whole simulation period are 0.48 m and 0.61, respectively. The deficiencies of RMSE and R² can be attributed to aggregation and simplification of parameters for the groundwater flow system (^{Qin et al., 2013}).

Figures 7(b)–7(e) demonstrate that the simulated groundwater levels can reproduce the measured data reasonably. The difference between the observation and simulation is possibly caused by the discrepancy between the actual location of pumping wells and approximate location represented by triangle nodes in the model. In the simulation, the local pumping wells were grouped to one well and assumed to pump water at one model grid. There are discrepancies in the spatial scales between the simulated and observed groundwater levels (^{Qin et al., 2013}). Other reasons include the lack of knowledge of the other potential sources/sinks, and the spatial heterogeneity that failed to be captured in our model.

Both measured and simulated results show declining trends of groundwater heads in the four wells during the calibration stage, while slightly rising trends in some wells during the validation stage. The intensive groundwater pumping activities should be the main reason for the decline of groundwater levels. The rise of groundwater levels after the fall of 2010 is caused by the higher precipitation in 2010, which allowed more recharge into the unconfined aquifer (^{BWA, 2009b}, ²⁰¹⁰).

As shown in Fig. 8(a), the simulated distribution of groundwater heads agrees with the measured distribution. Groundwater flows generally towards the east and north. The groundwater also flows towards the depression funnel in the northwest area caused by overexploitation.

Calibration and validation of solute concentrations

The simulated distributions of Cl⁻ and NO₃⁻ concentrations in the study area are in reasonable agreement with the observed data (Figs. 8(b) and 8(c)). The Cl⁻ concentration in the middle and southeastern part is higher than the other parts of the study area, mainly caused by the leaching from landfill sites located in those areas, such as the Beitiantang and Yonghezhuang landfill sites in Fengtai district, and the Xihongmen landfill sites in Daxing district. Both the Cl⁻ and NO₃⁻ concentration distributions are significantly affected by the initial concentration distribution.

Generally, the numerical model can capture the spatially and temporally varied groundwater levels in the study area. The results of calibration and validation indicate that the groundwater model can be applied for predicting the groundwater system responses to different water management strategies of the YRAF.

Discussion

Groundwater budget

Understanding the groundwater budget, i.e., the quantity of each water balance term, is essential for sustainable water management. Table 2 shows the groundwater budgets during the model calibration and validation periods. In general, the total inflow is 3.3674×10⁸ m³ while total outflow is 3.805×10⁸ m³. Accordingly, there is a decrease of water storage, i.e., water storage deficit, of 0.4376×10⁸ m³, for the aquifer from June 2009 to December 2010. The decreasing trend of groundwater storage in the study area agrees with most other studies conducted in this area. For example, ^{Xu (2006)} used MODFLOW2000 to simulate the groundwater flow of the Beijing plain; the water storage deficit of the phreatic aquifer for the year 2000 is 3.61×10⁸ m³. In the study by ^{Hao et al. (2014)}, the groundwater storage deficit for the year 2010–2011 is 1.1975×10⁸ m³. All these results indicate overexploitation in this area and thus water management strategies are needed to meet the water demand.

Precipitation and lateral inflow are the two largest inflow components, accounting for nearly 60% of the total inflow. The lateral inflow includes the lateral recharge from the western mountain front and subsurface flow from the northern boundary in our study. Groundwater exploitation is a dominant component of groundwater discharge and the amount of groundwater withdrawal is similar to the result reported by ^{Hao et al. (2014)}. The pumping water is mainly for urban water use, including both domestic water consumption and industrial use, while the percentage of urban groundwater use has increased due to the economic development and population growth (^{Qin et al., 2013}).

Recharge percentage of precipitation α

Precipitation accounts for the largest proportion of water inflow of the studied groundwater system. The rainfall recharge in our study from June 2009 to May 2010 is 0.491×10⁸ m³, indicating the average recharge percentage of precipitation to be 27%. This value was obtained by the calculation below:

(10)\alpha ={\displaystyle \frac{RR}{P\cdot A}}={\displaystyle \frac{\mathrm{0.491}\times {\mathrm{10}}^{\mathrm{8}}}{\mathrm{500}\times {\mathrm{10}}^{-\mathrm{3}}\times \mathrm{3.65}\times {\mathrm{10}}^{\mathrm{8}}}}=\mathrm{0.27},

where RR is the rainfall recharge (m³), P is the precipitation from June 2009 to May 2010 (mm), and A is the area of their study site (m²).

In the study of ^{Hao et al. (2014)}, the rainfall recharge from 2010 to 2011 is 1.1471×10⁸ m³ and the average ratio of rainfall recharge to precipitation, α , would be approximately 0.59, which may be too high. The reason is that although the western suburb of Beijing is characterized by Quaternary sediments with relatively high permeability, nowadays the study area has largely become urbanized. Furthermore, excessive pumping of groundwater has created a thick vadose zone (>20 m, up to 60 m in the center of depression cone). This implies that a large proportion of precipitation will remain in the unsaturated zone instead of reaching the groundwater (^{Zhu et al., 2013}). This was also confirmed by the results of many studies of this area with a variety of estimation methods. For example, the precipitation infiltration coefficient ranged from 0.01 to 0.4 in the study area over the past few decades, according to the empirical analysis by ^{Meng et al. (2013)}. Rainfall infiltration accounted for 12% of the rainfall in 2007 based on the water balance analysis of the Beijing plain (^{Zhu et al., 2013}). It was also confirmed by an environmental tracer technique (using Cl⁻) that groundwater recharge accounted for 7.2%–13.9% of precipitation and irrigation in Shijiazhuang, which is also located in the NCP (^{Lin et al., 2013}).

Concentration distributions of typical contaminants

There is a high degree of spatial and temporal variation of Cl⁻ and NO₃⁻ concentrations in the study area. The observed and simulated concentration distributions of Cl⁻ and NO₃⁻ agree with the results reported by ^{Wang et al. (2012)} that the groundwater quality in western part of Beijing is very poor. Two factors primarily contribute to the results. One reason is that the aquifer shows a high level of intrinsic vulnerability (^{Wang et al., 2012}). Unconsolidated pebbles and gravels are mainly distributed in YRAF and high aquifer permeability of the alluvial fan allows contaminants to percolate easily through the vadose zone and reach the groundwater. This area is an important groundwater recharge zone in the Beijing Plain (^{Zhou et al., 2013}). On the other hand, the infiltrating contaminant loads increase year by year due to human activities taking place in this area over a long period of time. Typical contaminant indicators show an increasing trend of concentration from the 1970s until the 2010s. The concentration in the middle YRAF, such as in the Fengtai district, is higher than in the upper YRAF, such as at Shougang, a steel production corporation (Figs. 9(a) and 9(b)). Cl⁻ concentration in the sample point of Fengtai increased rapidly from 43 mg/L in 1975 to 168 mg/L in 2010. Cl⁻ had a relatively low concentration of 34 mg/L in the sample point of Shougang in 1975, and increased rapidly from 1980 to 1990, and almost remained unchanged as 77 mg/L during 1990–2000 and then rose up to 122 mg/L in 2010.

Groundwater in YRAF is enriched in nitrogen. In the upper and middle parts of YRAF, 58% of NO₃⁻ in the unconfined water is from manure, 41% from chemical fertilizer, and only 1% from precipitation based on nitrogen isotope (¹⁵N-NO₃⁻) analysis (^{He et al., 2016}). The result implies that NO₃⁻ contamination in groundwater is basically human induced. Animal wastewater and urban wastewater introduce a large proportion of NO₃⁻ into the study area. The contaminant load of NO₃⁻ increased year by year, as shown in the pumping well groups No. 3 and No. 4 (Fig. 9(c)). The deteriorating trend of groundwater quality seems to have been controlled in recent years.

The high level of background concentration of Cl⁻ and NO₃⁻ should be considered during simulation of local groundwater contaminant to reflect the real conditions of groundwater quality. Although we could only conduct a basic calibration of solute transport given the limited data available, the observed and simulated concentration distributions of Cl^‒ and NO₃⁻ are in reasonable agreement. This fact provides a foundation for the future prediction of the influence of reclaimed water used as ecological water of the Yongding River on the local groundwater environment.

Mechanism of nitrogen transforms

The N cycling in aquifers involves a complex set of processes and reactions. Besides the three primary transformations of N considered in the current study, the other N involved reactions can also occur in the local aquifer, such as cation exchange of NH₄⁺. The concentration of each species is comprehensively determined by a series of reactions. Nitrification produces acid which enables mineral dissolution in the aquifer. As a result, Mg²⁺, Ca²⁺, Na⁺, and other metal cations can be dissolved by water-rock interactions (^{Yin et al., 2016}). Meanwhile, the adsorbed cations can be released by introducing NH₄⁺ into groundwater, so the TH of the groundwater tends to increase. The increasing trend of TH is obvious in the unconfined water of YRAF as shown in Figs. 9(b) and 9(c). The TH concentration is higher in the middle part of YRAF, such as in Fengtai and Panjiamiao than the upper part of the study area. The rising groundwater hardness is linearly correlated with nitrate nitrogen in YRAF (Fig. 9(d)).

In our study, denitrification is assumed to occur widely in the model scope. However, denitrification can be inhibited when the Dissolved Oxygen (DO) concentration is high in shallow groundwater. The threshold DO concentration for denitrification to occur is 2 mg/L (^{Rivett et al., 2008}). According to the survey conducted by ^{He et al. (2016)}, DO concentration varied in the range of 6 to 7 mg/L in some areas, which means that denitrification can hardly occur in the aerobic environment. The results imply that simulated results might underestimate NO₃⁻ concentration in some parts of the aquifer and call attention to control nitrogen contaminant.

Conclusions

A groundwater flow and reactive solute transport model developed in this study permits the delineation of flow and contaminant transport in the western suburbs of Beijing. The model was calibrated and validated against a time series of groundwater head data from four observation wells, with the average RMSE and R² of 0.48 m and 0.61. Water budget analysis indicates that 60% of the inflow water enters the model area through precipitation and lateral inflow, while 87% of the outflow water left the hydrological system through pumping, with a small amount leaving the system through lateral outflows. In the simulation period from June 2009 to December 2010, a negative groundwater budget of 0.4376×10⁸ m³ was shown in this area, which reveals unsustainable water resource development. The high background concentrations of Cl⁻ and NO₃⁻ in the groundwater associated with industry development in the late 1990s suggests the necessity of taking into account the impact of human activities on the local groundwater environment. The average rainfall infiltration coefficient of 0.27 provided a reliable value for estimating the amount of rainfall infiltration in the study area.

To summarize, a groundwater reactive transport model of the study area was established for the first time. The evaluation of the response of the local groundwater environment after reclaimed water discharging into the Yongding River will be carried out based on the model.

References

[1] Bear J (1977). On the aquifer integrated balance equations. Adv Water Resour, 1(1): 15–23

[2] Beijing Water Authority (BWA) (2009a). Green Yongding River: the construction plan for an ecological corridor. Beijing, China (in Chinese)

[3] Beijing Water Authority (BWA) (2009b). Beijing water resources bulletin in 2009. Beijing, China (in Chinese)

[4] Beijing Water Authority (BWA) (2010). Beijing water resources bulletin in 2010. Beijing, China (in Chinese)

[5] Cao G, Zheng C, Scanlon B R, Liu J, Li W (2013). Use of flow modeling to assess sustainability of groundwater resources in the North China Plain. Water Resour Res, 49(1): 159–175

[6] Cheng D (2007). Hydrogeochemical Processes and Numerical Simulation of Nitrate Nitrogen and Total Hardness in Groundwater of Beijing Urban. Dissertation for PhD Degree. China University of Geosciences (Beijing) (in Chinese)

[7] Collins S, Singh R, Rivas A, Palmer A, Horne D, Manderson A, Roygard J, Matthews A (2017). Transport and potential attenuation of nitrogen in shallow groundwaters in the lower Rangitikei catchment, New Zealand. J Contam Hydrol, 206: 55–66

[8] Cui Z, Welty C, Maxwell R M (2014). Modeling nitrogen transport and transformation in aquifers using a particle tracking approach. Comput Geosci, 70: 1–14

[9] Dong S (2010). Nitrogen Numerical Simulation in Beijing Urban Groundwater and Monitored Natural Attenuation about Groundwater. Dissertation for Master Degree. China University of Geosciences (Beijing) (in Chinese)

[10] Feng S, Huo Z, Kang S, Tang Z, Wang F (2011). Groundwater simulation using a numerical model under different water resources management scenarios in an arid region of China. Environ Earth Sci, 62(5): 961–971

[11] Guo G, Li Y, Xu L, Li Z, Yang Q, Xu M (2014). Risk assessment of quaternary groundwater contamination in Beijing plain. Environ Sci (Ruse), 35(2): 562–568 (in Chinese)

[12] Hao Q, Shao J, Cui Y, Xie Z (2014). Applicability of artificial recharge of groundwater in the Yongding River alluvial fan in Beijing through numerical simulation. J Earth Sci, 25(3): 575–586

[13] He G, Liu P, Mu X, Wu Q, Liu M (2016). Identification of nitrate sources in groundwater in the Yongding River alluvial fan with isotope technology. J Hydraul Eng, 47(4): 582–588 (in Chinese)

[14] Hu H, Mao X, Barry D A, Liu C, Li P (2015). Modeling reactive transport of reclaimed water through large soil columns with different low-permeability layers. Hydrogeol J, 23(2): 351–364

[15] Jiang B, Wong C P, Lu F, Ouyang Z, Wang Y (2014). Drivers of drying on the Yongding River in Beijing. J Hydrol (Amst), 519: 69–79

[16] Lin D, Jin M, Liang X, Zhan H (2013). Estimating groundwater recharge beneath irrigated farmland using environmental tracers fluoride, chloride and sulfate. Hydrogeol J, 21(7): 1469–1480

[17] Lin J (2004). The Analysis of Pollution History for the Groundwater in Urban and Suburb Area of Beijing. Dissertation for Master degree. Changchun: Jilin University (in Chinese)

[18] Lin J, Yang Q, Yang Q, Liu F J, Chen Z, Dou W, Jiang Y, Liu Z, Zhao W, Ni T (2011). Impacts of Reclaimed Water Reuse in the Typical River Channel of the Yongding River on the Groundwater and Model Establishment. Research report, Beijing: Beijing Institute of Hydrogeology and Engineering Geology (in Chinese)

[19] Liu C, Zhang F, Zhang Y, Song S, Zhang S, Ye H, Hou H, Yang L, Zhang M (2005). Experimental and numerical study of pollution process in an aquifer in relation to a garbage dump field. Environmental Geology, 48(8): 1107–1115

[20] Liu M, Seyf-Laye A S M, Ibrahim T, Gbandi D B, Chen H (2014). Tracking sources of groundwater nitrate contamination using nitrogen and oxygen stable isotopes at Beijing area, China. Environ Earth Sci, 72(3): 707–715

[21] Liu Y (2012). Study on the Groundwater Reservoir Storage Capacity and Utilization Mode in Western Suburb of Beijing. Dissertation for Master degree. Changchun: Jilin University (in Chinese)

[22] Meng S, Fei Y, Zhang Z, Lei T, Qian Y, Li Y (2013). Research on spatial and temporal distribution of the precipitation infiltration amount over the past 50 Years in North China Plain. Advances in Earth Science, 28(8): 923–929 (in Chinese)

[23] Peng J, Li S, Qi L (2013). Study on river regulation measures of dried-up rivers of Haihe River basin, China. Water Sci Technol, 67(6): 1224–1229

[24] Qin H, Cao G, Kristensen M, Refsgaard J C, Rasmussen M O, He X, Liu J, Shu Y, Zheng C (2013). Integrated hydrological modeling of the North China Plain and implications for sustainable water management. Hydrol Earth Syst Sci, 17(10): 3759–3778

[25] Rivett M O, Buss S R, Morgan P, Smith J W N, Bemment C D (2008). Nitrate attenuation in groundwater: a review of biogeochemical controlling processes. Water Res, 42(16): 4215–4232

[26] Sheibley R W, Jackman A P, Duff J H, Triska F J (2003). Numerical modeling of coupled nitrification–denitrification in sediment perfusion cores from the hyporheic zone of the Shingobee River, MN. Adv Water Resour, 26: 977–987

[27] Soldatova E, Guseva N, Sun Z, Bychinsky V, Boeckx P, Gao B (2017). Sources and behaviour of nitrogen compounds in the shallow groundwater of agricultural areas (Poyang Lake basin, China). J Contam Hydrol, 202: 59–69

[28] Stuart M E, Lapworth D J (2016). Macronutrient status of UK groundwater: nitrogen, phosphorus and organic carbon. Sci Total Environ, 572: 1543–1560

[29] Sun F, Shao H, Kalbacher T, Wang W, Yang Z, Huang Z, Kolditz O (2011). Groundwater drawdown at Nankou site of Beijing Plain: model development and calibration. Environ Earth Sci, 64(5): 1323–1333

[30] Sun F, Shao H, Wang W, Watanabe N, Bilke L, Yang Z, Huang Z, Kolditz O (2012). Groundwater deterioration in Nankou—A suburban area of Beijing: data assessment and remediation scenarios. Environ Earth Sci, 67(6): 1573–1586

[31] Wang J, He J, Chen H (2012). Assessment of groundwater contamination risk using hazard quantification, a modified DRASTIC model and groundwater value, Beijing Plain, China. Sci Total Environ, 432: 216–226

[32] Wang L, Wang Z, Koike T, Yin H, Yang D, He S (2010). The assessment of surface water resources for the semi-arid Yongding River Basin from 1956 to 2000 and the impact of land use change. Hydrol Processes, 24(9): 1123–1132

[33] Wang S, Shao J, Song X, Zhang Y, Huo Z, Zhou X (2008). Application of MODFLOW and geographic information system to groundwater flow simulation in North China Plain, China. Environmental Geology, 55(7): 1449–1462

[34] Wang S, Zheng W, Currell M, Yang Y, Zhao H, Lv M (2017). Relationship between land-use and sources and fate of nitrate in groundwater in a typical recharge area of the North China Plain. Sci Total Environ, 609: 607–620

[35] Xu X (2006). A study of Numerical Simulation for Groundwater Flow in the Beijing Plain. Dissertation for Master Degree. China University of Geosciences (Beijing) (in Chinese)

[36] Yin S, Wu W, Liu H, Bao Z (2016). The impact of river infiltration on the chemistry of shallow groundwater in a reclaimed water irrigation area. J Contam Hydrol, 193: 1–9

[37] Yu C, Yao Y, Cao G, Zheng C (2015). A field demonstration of groundwater vulnerability assessment using transport modeling and groundwater age modeling, Beijing Plain, China. Environ Earth Sci, 73(9): 5245–5253

[38] Yu Y, Ma M, Zheng F, Liu L, Zhao N, Li X, Yang Y, Guo J (2017). Spatio-temporal variation and controlling factors of water quality in Yongding River replenished by reclaimed water in Beijing, North China. Water, 9(7): 453

[39] Zhai Y, Wang J, Huan H, Zhou J, Wei W (2013). Characterizing the groundwater renewability and evolution of the strongly exploited aquifers of the North China Plain by major ions and environmental tracers. J Radioanal Nucl Chem, 296(3): 1263–1274

[40] Zhang S, Wang C, Meng X, Hua D, Men B, Li Z (2013). Evaporation study in Beijing section of the Yongding River. Progress in Geography, 32(4): 580–586 (in Chinese)

[41] Zhou Y, Dong D, Liu J, Li W (2013). Upgrading a regional groundwater level monitoring network for Beijing Plain, China. Geoscience Frontiers, 4(1): 127–138

[42] Zhu L, Liu C, Li X, Guo G, Pan Y (2013). Precipitation infiltration change in Beijing Plain in the context of urbanization. Journal of China University of Geosciences (Earth Science edition), 38(5): 1065–1072

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature