Evaluating the interactions between surface water and groundwater in the arid mid-eastern Yanqi Basin, northwestern China

ABSTRACT

The environment of Bosten Lake in the Mid-Eastern Yanqi Basin (MEYB), an arid inland area in northwest China, has deteriorated greatly due to increasing groundwater exploitation and changes in the interactions between groundwater and surface water. This study intended to simulate the spatio-temporal variability of groundwater and surface water across the entire MEYB over the period 2000–2013. The applicable groundwater flow model and mass balance calculation method for river water were constructed to evaluate the change in groundwater recharged by and discharged to different segments of the Kaidu River. Simulation results show that the entire river seepage in the MEYB increased from 1.05 to 6.17 × 10⁸ m³/year between 2000 and 2013. The increasing river seepage, induced by increasing groundwater exploitation, plays the most important role in the water level decline in the downstream reaches of the Kaidu River and in Bosten Lake. This implies that the current utilization of groundwater resources in the MEYB is unsustainable.

KEYWORDS:

• surface water–groundwater interaction
• exploitation
• numerical modelling
• Kaidu River
• Bosten Lake

Editor D. Koutsoyiannis Associate editor T. Okruszko

1 Introduction

Groundwater and surface water are essential to agricultural, industrial, domestic and environmental uses, and thus have important effects on economic development (Glazer and Likens 2012, McCallum et al. 2012, Guay et al. 2013, Unland et al. 2013, Elsawwaf et al. 2014, Matula et al. 2014). The Mid-Eastern Yanqi Basin (MEYB) is a typical arid inland basin in Bayingolin Mongolian Autonomous Prefecture, where groundwater and surface water depletion has become a critical problem affecting major regions of Xinjiang Uyghur Autonomous Region in northwestern China (Figs 1 and 2). Due to its relatively suitable conditions, the MEYB is the main grain production region in southern Xinjiang and the main water-supply source area to Korla City (Ouyang 2008, Zhong 2008, Du 2009, Tao 2012, Wang 2013, 2014).

Figure 1. Location and topography of the Mid-Eastern Yanqi Basin. I–I′ and II–II′ are the lines of cross-sections in Figure 4.

Figure 2. (a) Elevation of the study region, and (b) distribution of diversion ditches, drainage ditches and division of the river channel in the MEYB.

Bosten Lake lies in the eastern MEYB, and is the biggest inland freshwater lake in China, covering an area of 1002.4 km². As the most important river in the MEYB, the Kaidu River flows through the MEYB into Bosten Lake, with an annual flow of 34.86 × 10⁸ m³/year. Bosten Lake is also the origin of the Kongque River (Wang et al. 2014). As the unique outflow stream of Bosten Lake, the Kongque River is the only viable source of water supply for the 0.85 million downstream population. Further south of the Kongque River, the lower reaches of the Tarim River flow through the Tarim basin and into Luobupo Lake in Xinjiang Province. Over the past half century, unbridled water utilization of the Tarim River has led to serious water scarcity and ecological damage. The most seriously adverse effect is on the lower reaches of the Tarim River, where significant groundwater level decline has occurred and the ecosystem has deteriorated. To retrieve the natural environment of the lower reaches of the Tarim River, the Ecological Water Conveyance Project (EWCP) was implemented in 1999 (Ye et al. 2009). Under the EWCP, Bosten Lake, as the emergency water source, has been diverted through the Kongque River to the lower reaches of the Tarim River since 2000. Consequently, Bosten Lake has become the main water supply source for those living in the lower reaches of the Tarim and for its ecology; thus, it is known as “the Mother Lake” of the Bayingolin Mongolian Autonomous Prefecture (Ouyang 2008, Du 2009).

With the rapid development of agriculture, water from the diversion can hardly satisfy the water consumption over the MEYB, so groundwater is exploited. However, the increased groundwater exploitation has changed the interaction between the Kaidu River and the groundwater, and has led to serious problems of water resources and the Bosten Lake environment. Therefore, it is important to study the interaction between the Kaidu River and the groundwater in the MEYB to achieve sustainable development and protect the water resources and environment.

It has been 57 years since the meteorological stations and hydrological stations in the MEYB were built. Then, many projects were launched to conduct hydrogeological reconnaissance, surface water observation, groundwater observation and hydrological budget research. Consequently, various researches were implemented to study soil water, surface water and groundwater in the MEYB. Based on the evaluation of water resources in the MEYB, Liu (2004) built a regional groundwater model to calculate the amount of water exchange between surface water and groundwater. However, this model was limited to short-term variation of interaction between surface water and groundwater, being for only one year. Soil water and salt transportation in a typical experimental area of the MEYB were also studied by Wang (2005). Furthermore, Ouyang (2008) built a water and salt balance model to evaluate the shallow groundwater and irrigation schedule in the MEYB at the field scale. Similarly, a regulation-control model for the whole Yanqi Basin was developed by Zhong (2008) to investigate the water cycle and salt transport in Bosten Lake, although the regional variation of groundwater and surface water is hard to derive from this work. Using chemical and isotopic data, Wang (2013, 2014) determined the groundwater age and its recharge conditions, and analysed the groundwater sustainability. The temporal and spatial variation of water stable isotopes in the Kaidu River was investigated by Wang et al. (2014) to identify the hydrological processes and hydrological composition of groundwater in the MEYB. In spite of extensive work and substantial progress in this area, the interaction between surface water and groundwater at the regional scale is still not clear.

Identifying the surface water–groundwater interaction in a watershed is essential for understanding the water cycle and sustainable integrated water resources utilization (Glazer and Likens 2012, McCallum et al. 2012, Guay et al. 2013, Unland et al. 2013, Elsawwaf et al. 2014, Matula et al. 2014). The Navier-Stokes-based shallow water equations for surface water flow and Darcy’s equation for groundwater flow were coupled and applied to a linear Stokes problem by Marco et al. (2002). Considering precipitation, upstream runoff and cultivated land area, a multifactor hydrological response model was developed by Wang et al. (2007) to investigate the responses of hydrological processes to land use in the middle basin of the Heihe River. The stable isotope technique was used to trace a major source of mountain runoff to the Heihe River and make important observations to identify the water cycle process (Wang et al. 2009). To investigate hydrological response to multimodel climate output, Sulis et al. (2012) conducted a process-based model of integrated surface water and groundwater. The ²²²radon technique was applied to evaluate surface water–groundwater connectivity in arid environments (Oyarzún et al. 2014). Benthic fluxes were simulated by Sawyer et al. (2013) to investigate the dynamic response of surface water–groundwater exchange in a shallow estuary. Through physico-chemical in-stream profiling, spatial groundwater–surface water interactions were delineated in small headwater catchments by Haria et al. (2013).

Numerical models are applied to assess groundwater–surface water interactions and play an increasingly important role. Based on MODFLOW computer code, a mathematical model was used to estimate the artificial groundwater recharge quantities for a semi-arid basin (Al-Assa’d and Abdulla 2010). Similarly, groundwater evapotranspiration was estimated using MODFLOW models in a semi-arid region (Ajami et al. 2011). A groundwater flow model was developed by Yao et al. (2014) to assess the flow exchanges between different parts of the water system in the Heihe Basin, and this provided an effective means to understand the complex flow system. For underground coal mining, a flow model had an important effect on evaluating the rapid infiltration of surface water and slow release of groundwater from the mine ore body (Cravotta et al. 2014). A transient, integrated hydrological model was applied to a hard-rock catchment by Hassan et al. (2014) to study the dynamics of surface water–groundwater interactions in a complicated hard-rock system. By measuring the hydraulic head in the riparian zone near rivers and the fluid pressure in the river beds, the connectivity between rivers and aquifers was explored in a semi-arid river basin by Lamontagne et al. (2014). Furthermore, the Soil and Water Assessment Tool (SWAT) and MODFLOW were coupled to enhance the ability to understand groundwater–surface water interactions, in which the coupled SWAT-MODFLOW model was applied to simulate surface hydrology and groundwater hydrology accurately for a river system in a semi-arid region (Bailey et al. 2016).

In this study, we develop a numerical model to study surface water and groundwater interaction in the MEYB, the main irrigation area on the west bank of Bosten Lake. In contrast to the groundwater flow model on a small time scale and in limited space (cf. Liu 2004, Wang 2004, Ouyang 2008, Du 2009), we built a regional-scale, long-term model to evaluate the mass exchange between surface water and groundwater between 2000 and 2013. Investigations and field experiments were conducted to obtain the spatial distribution of the surface water system and of hydrogeological parameters. Then, a conceptual and numerical model based on MODFLOW-NWT was established and calibrated with the long series of observation data. The calibrated model was used to investigate the interactions between surface water and groundwater in the MEYB between 2000 and 2013. Meanwhile, the water balance of the whole river channel was computed using the Streamflow-Routing (SFR2) package (Niswonger and Prudic 2005) of MODFLOW-NWT in the study region to evaluate the change in groundwater recharge and discharge in various segments of the Kaidu River over a long period of time. In this study, the surface water and groundwater are treated as an entity and analysed systematically and comprehensively at regional and long-term scales. Current understanding of the MEYB is limited by spatio-temporal variability in depletion at the small scale and the short-term variation of dynamics of the surface water and groundwater. Thus, this study is essential for better understanding of the hydraulic links between the groundwater and surface water systems. Detailed terms in the water budget, the correlations among them, the interaction between surface water and groundwater in six different river segments, and sensitivity analysis of flow model were explored comprehensively to derive the main factors leading to water depletion and deterioration of Bosten Lake. An effective management tool was developed to help estimate the effects of groundwater exploitation and distribute water resources for water supply and environmental protection.

2 Site description and hydrogeology

The MEYB is the main irrigation area on the western shore of Bosten Lake. The location of the MEYB and its boundaries (41°48ʹ–41°24ʹN; 85°54ʹ–86°54ʹE) is shown in Figure 1. The basin includes Hejing, Yanqi and Bohu counties. It is located in the loose Quaternary sediments of the South Tianshan Mountains, on the middle and lower reaches of the Kaidu River. As the MEYB is situated in inner Eurasia and is remote from the ocean, the environmental and hydrological conditions are influenced by a continental climate with dry weather, dry air and low temperatures, where precipitation is quite low and evaporation is intensive.

There are significant differences in climate between the mountain and plain areas across the Yanqi Basin. There is really no summer season in the mountain area, where the lowest elevation is higher than the semi-permanent snow line and therefore perpetually covered with snow and ice. The annual precipitation in the mountain area ranges from 300 to 500 mm (Liu 2004, Du 2009). In the plain region, the seasons are demarcated by a short summer and a long and dry winter. The average annual precipitation is 71.3 mm, most of which occurs between May and September. Corresponding to this climate, however, the average pan evaporation reaches 2120.9 mm annually, most of which is concentrated between March and October. Meanwhile, the ratio of open water evaporation to precipitation is more than 30, causing water scarcity all over the MEYB (Liu 2004, Zhong 2008, Du 2009).

As the biggest inland freshwater lake in China (Fig. 1), Bosten Lake consists of three parts: big lake, small lake and swamp, with a total area of about 1002.4 km² and an amount of surface water of approximately 4.0 × 10⁹ m³. In the MEYB, diversion ditches, drainage ditches and rivers crisscross in the fields (Fig. 2(b)); the Kaidu River is the largest river, and the only perennial river that supplies Bosten Lake. There are three hydrological stations or major hydro-junctions along the Kaidu River in the MEYB: the first hydro-junction, Yanqi hydrological station, and Baolangsumu hydro-junction (Fig. 2(b)). The annual discharge from the Kaidu River to Bosten Lake is 34.8 × 10⁸ m³ (2000–2013), accounting for over 80% of the total annual recharge volume of Bosten Lake.

For the Kaidu River, monthly runoff varies tremendously. The dry season is from October to the following March, and the flood season from April to September. In April–May, the runoff amount of the Kaidu River is about 19% of the annual amount, with peak runoff in June–August (45.5%). Then, in September–November, the runoff decreases and recharge comes mainly from groundwater in the whole river basin. The driest period is between December and the following March, with runoff of 14.7% of the annual amount (Ouyang 2008, Zhong 2008, Du 2009, Qian et al. 2006, Tao 2012). The annual variation in Kaidu River runoff at different hydrological stations is given in Table 1.

Table 1. Annual variation in Kaidu River runoff at different stations (2000–2013).

	Daily runoff (10^6 m^3)				Annual runoff (10^9 m^3)
Hydrological station	April–May	June–August	September–November	December–March	Mean	Maximum	Minimum
Dashankou	10.8	17.2	6.2	4.2	3.486	5.712	2.461
Yanqi	7.1	12.3	3.6	4.4	2.584	4.349	1.457
Baolangsumu	6.1	12.5	3.0	4.6	2.502	3.546	0.722

The MEYB is the main agricultural production region of southern Xinjiang, and irrigation water is the major source of water consumption due to the arid climate. Before the 1950s, the agricultural practice was carried out on a small scale and the hydrological cycle was kept in a nearly original or natural state. Since the 1950s, agriculture has been developing rapidly, with large tracts of wasteland being reclaimed. Due to non-uniform water distribution, the agricultural fields have to be irrigated by water diverted from the Kaidu River; as a result, the diversion system and drainage system covering the whole MEYB were built quickly (Ouyang 2008, Zhong 2008, Du 2009, Qian et al. 2006, Tao 2012). With the rapid growth of agriculture, water from diversion can hardly satisfy the water consumption over the whole MEYB, so groundwater has to be exploited.

Groundwater in the MEYB exists in Quaternary alluvial aquifers, which comprise laterally discontinuous layers of gravel, sand and clay (Figs 3 and 4). The media in the aquifers vary from coarse particles to fine particles oriented from the piedmont to Bosten Lake, while the aquifers also turn from single-layer to multiple layers. In the piedmont area, sand gravel and sand are distributed alternately, and the groundwater recharges from lateral inflow, river flow, and precipitation infiltration. An aquifer with good permeability generally has a large amount of groundwater resources and excellent water quality. In the alluvial plain, the coarse sandy gravel media are interbedded with silt and silt clay (Fig. 4), and the depth to water table is so shallow (1–3 m below ground level) that evaporation from phreatic water is quite intense and is an important part of the groundwater discharge.

Figure 3. Schematic hydrogeological cross-section and groundwater system oriented west to east through the aquifer of the entire Yanqi Basin.

Figure 4. Hydrogeological cross-sections of the MEYB along the lines I–I′ and II–II′ (see Fig. 1).

The area of irrigated farmland has been increasing since the beginning of the 21st century, and this has led to more intensive groundwater exploitation. In the MEYB, the cultivated land area has increased from 6700 to 120 000 hm² since the 1950s (Liu 2004). In those areas irrigated by groundwater, irrigation activity has caused groundwater overexploitation. Consequently, the water table level is deepened year by year. The increasing groundwater exploitation has become an important limiting factor to social development and sustainability (Ouyang 2008, Zhong 2008, Du 2009, Gleeson et al. 2010, Qian et al. 2006, Tao 2012). Water diverted from the Kaidu River for crop irrigation is also increasing year by year, and this has been more than 30% of the runoff amount, simultaneously; thus, more and more groundwater is being exploited. As a result, the interaction between the Kaidu River and the groundwater has changed, while the flow rate from the Kaidu River into Bosten Lake is decreasing continuously. The Bosten Lake covered an area of 1262.3 km² in 2002, which decreased to 800 km² in 2013, with a decline in water level of up to 4.34 m (2002–2013). Simultaneously, the salinity of Bosten Lake has gradually increased from 1.17 g/L in 2002 to more than 1.5 g/L in 2013. Now, Bosten Lake has become a brackish lake (Liu 2004, Zhong 2008). In the same period, many precious species have disappeared or come close to extinction, leading to a decline in biological diversity. Furthermore, as the origin of the Kongque River and water supply source for the lower reaches of the Tarim, the water level decline of Bosten Lake is threatening the water safety of the Tarim lower reaches and its watershed ecology. This will adversely affect the ecological health of the Tarim Basin and restrict the sustainable development of both society and the economy in southern Xinjiang (Ye et al. 2009, Pang et al. 2010, Aishan et al. 2013, Chen et al. 2013). The combination of these factors makes the interaction between the Kaidu River and groundwater in the MEYB an important research issue.

3 Methodology

3.1 Model set-up and discretization

The governing equation for the flow model in the study region can be expressed as:

(1)

where h is the hydraulic head [L]; K is hydraulic conductivity [L T⁻¹]; B is the elevation of the bottom of the aquifer [L]; q_s is the volumetric flow rate of fluid sinks/sources per unit volume of aquifer [T⁻¹]; S is the specific yield (dimensionless); and t is time [T].

The flow model is discretized spatially using a finite-difference method, as illustrated in Figure 5(a), which contains 254 columns and 120 rows with a uniform spacing of 500 m. Note that the modelled domain contains Bosten Lake, and the lake grid is set as inactive in the model while the water balance of Bosten Lake is being computed. All of the pumping wells in the study region are fully penetrating wells, penetrating multiple layers from the shallow aquifer through to the deep one so that the aquifers can be vertically grouped and considered as a single unconfined groundwater flow system. Correspondingly, the grids in model are composed of one single layer, specified as an unconfined type vertically to simulate a two-dimensional regional flow system. The Layer Property Flow (LPF) package was utilized to assign both horizontal and vertical hydraulic conductivity values to the aquifers. The horizontal distribution of hydraulic conductivity (Fig. 5(b)) is interpolated by a kriging method according to data obtained from 196 pumping tests in boreholes distributed all over the study region. In the LPF package, the vertical hydraulic conductivity values are specified as one-tenth of the values in the horizontal direction. Similarly, the specific yield ranges from 0.04 to 0.26 in the study region (Fig. 5(c)). The top and bottom elevations of the model are derived by interpolation of the data from digital elevation and the 1:200 000 hydrogeological map of the MEYB. For the transient model, the total simulation time is 14 years, discretized into 27 stress periods based on the irrigation scheme; every irrigation period in April–October is a stress period, as is November to the following March. The groundwater flow model is based on and solved by MODFLOW-NWT (Niswonger et al. 2011).

Figure 5. (a) Spatial discretization of the modelled domain and the prescribed boundary conditions; (b) distribution of hydraulic conductivity; and (c) distribution of specific yield in the model.

3.2 Boundary conditions and sink/source terms

The boundary setting of the model is illustrated in Figure 5(a). The interface between the plain and mountain terrains lies on the northwestern border of the model, where a significant amount of mountain front groundwater discharges laterally into the study region; thus the northwestern border of the model is defined as a specified-flux boundary based on the Flow and Head Boundary (FBH1) package (Leake and Lilly 1997). There are seven observation wells distributed at the northwestern boundary (Fig. 5(a)), where the hydraulic gradient is 0.0021–0.0054. Also, the hydraulic conductivity values at the seven observation points range from 15.0 to 81.5 m/d. Then the specified flux can be calculated as the product of hydraulic conductivity and hydraulic gradient. On the northern and western borders, the boundaries are defined as specified heads due to the existence of some observation points. Along the northeastern border is swamp area where groundwater reaches the land surface, so specified heads of the border are based on the ground elevation. The big lake and small lake of Bosten Lake are located on the eastern and southern edges, respectively, and were defined using the LAK and GHB packages, respectively. The big lake area of Bosten Lake was modelled using the LAK package (Merritt and Konikow 2000) of MODFLOW-NWT, in which both groundwater level and lake stage respond interactively to each other until an equilibrium state is reached. It should be noted that the bottom elevations of the lake are specified in the discretization package and no upper limit of lake stage is required by LAK logically. The lake girds are assigned as inactive cells in the model and the sink/source terms (e.g. precipitation, evaporation) are input to LAK (Fig. 6(c)). The runoff entering the lake and water diversion from the lake are illustrated in Figure 6(d); the runoff of the Kaidu River into Bosten Lake is calculated using the SFR2 package in our model and the other runoffs into Bosten Lake are aggregated as a direct input to the LAK package.

Figure 6. (a) Zonation of recharge coefficient of precipitation infiltration; (b) distribution of evapotranspiration coefficient; (c) precipitation and evaporation, and inflow/outflow of Bosten Lake; and (d) inflow and diversion of main segments within the study region.

As an open border, the top boundary is implemented using the UZF package (Niswonger et al. 2006) in MODFLOW-NWT, and the areal sink/source terms, including precipitation recharge, irrigation infiltration and evaporation discharge, are imposed to the model. The saturated water content beneath the land surface is set as 0.3, and a Brooks-Corey exponent of 3.5 for the unsaturated zone is assigned to all aquifer grids. In this study area, the effective precipitation can be calculated approximately by surface pan precipitation monitored at weather stations and hydrological stations minus interception of 3 mm/month by the vegetation canopy (Du 2009). Then the actual recharge to the aquifer from precipitation is equal to the effective precipitation multiplied by the recharge coefficient of precipitation infiltration, as shown in Figure 6(a). Moreover, the return water from irrigation infiltration can be accurately calculated by multiplying the irrigation quota of 3.106–6.715 × 10⁻³ m/d by the irrigation infiltration coefficient of 0.10–0.25, as distributed in our study region (Du 2009). Therefore, the total amount of areal recharge in the irrigation area is the sum of infiltration from both precipitation and irrigation. As for the evaporation estimation, the distribution of the evapotranspiration coefficient, as shown in Figure 6(b)), is determined by land use and vegetation coverage over the whole study region and the maximum evapotranspiration (ET_max) is equal to the evaporation from the water surface multiplied by the evapotranspiration coefficient. The calculated maximum ET ranges from 7.0058 × 10⁻⁶ to 7.5981 × 10⁻³ m/d in the study region, while the ET extinction depth is assigned as 5.0 m for all grids.

Besides the areal sinks/sources, drainage from ditches and channels can be treated as linear sink/source terms. In this study, the drainage ditches are simulated by the Drain package (DRN) of MODFLOW. The depth of the drainage ditches is 1.39–4.6 m and the width is 8–12 m; the drain conductance is about 5.58–68.9 m²/d. In contrast, the Kaidu River, Wulasitai River and other diversion ditches are simulated by the SFR2 package (Niswonger and Prudic 2005) of MODFLOW-NWT. The Kaidu River is divided into five segments (Segments I–V) and the Wulasitai River (Segment VI) and each ditch are set as individual segments. The cross-section of each segment is assumed to be a wide rectangular shape, the stream depth of which is calculated by Manning’s equation with Manning’s roughness coefficient of 0.03. The width of the segments of the Kaidu River ranges from 69 to 114 m, that of the diversion ditches ranges from 3 to 19 m, while the width of the Wulasitai River is set to a constant 22 m. Furthermore, the hydraulic conductivity of the stream bed ranges from 1.39 to 49.27 m/d, and the stream-bed thickness is about 3.0 m for both the Kaidu and Wulasitai rivers and 1.0m for the diversion ditches, according to field surveys and experiments. As shown in Figure 2(b), the runoff from the upstream reaches enters Segment I of the Kaidu River at the first hydro-junction, and then the outflow at the downstream end of Segment I continues flowing into segments II–V in sequence, finally flowing into Bosten Lake. Meanwhile, the outflow of Segment VI is also entering Segment II of the Kaidu River, but Segment VI belongs to the Wulasitai River. In addition, the inflows of the diversion ditches are all diverted from segments I–IV of the Kaidu River, which is measured in the period 2000–2013 and shown in Figure 6(d)).

The bottom of the aquifer is specified as a no-flow boundary and groundwater exploitation can be clarified by the Well package, where input data are from the pumping rate at well locations throughout the study region.

3.3 Model calibration and validation

To measure the optimality of model calibration and validation, various quantitative performance criteria were used (Lemann and Casella 1998, Krause et al. 2005, Hyndman and Koehler 2006, SBTSC (State Bureau of Technical Supervision of China) 2008):

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

where O_i is the ith observation for the constituent being evaluated; S_i is the ith simulated value; e_i is the error term, which is the difference between observed and simulated values; is the mean of observed data for the constituent being evaluated; is the mean of all simulated values for the constituent being evaluated; and N is the total number of observation data. Equation (2) shows the error terms, being the difference between observed and simulated values; the relative error (RE; Equation (3)) is the ratio of the absolute value of error terms to the observed value; the mean error (ME; Equation (4)) is the arithmetic average of all error terms; the mean absolute error (MAE; Equation (5)) is another popular error measure that averages the absolute value of the error terms; the mean squared error (MSE; Equation (6)) is the arithmetic mean of the sum of the squares of the error terms; the root mean square error (RMSE; Equation (7)) is the square root of the MSE; the standard deviation of residuals (STD_res; Equation (8)) is defined as the square root of the variance, which measures the amount of variation or dispersion from the average; the coefficient of correlation (R; Equation (9)) is the ratio between the covariance and the multiplied standard deviations of the observed and simulated values (Lemann and Casella 1998, Zuo et al. 2015); and the Nash-Sutcliffe efficiency (E; Equation (10)) is the sum of the absolute difference between the simulated and observed values normalized by the variance of the observed values (Garambois et al. 2015, Zuo et al. 2015).

The calibration process was divided into two steps. First, the steady state was calibrated to the known groundwater field in 1999 when groundwater exploitation was not intensive. In this step, all parameters of the groundwater flow model, such as the hydraulic conductivity distribution and the value of the specified flux in the northwestern border, are adjusted to make the calculated hydraulic heads fit well with the measurements. Since the contour map of the observed water table prior to 1999 is available, the calibration of the steady-state model focuses on visual comparison of calculated and observed contour maps (Fig. 7(a)).

Figure 7. (a) Comparison of observed and calculated hydraulic head contours in 1999; (b) comparison of observed and computed groundwater heads at all observation wells; (c) groundwater level contours in November 2013; and (d) observed vs computed heads over time at three typical observation locations, indicated in (a).

In the second step, based on the preliminary adjusted parameter distribution and boundary conditions, the initial hydraulic heads for the groundwater transient model from 2000 to 2013 were based on the simulation results from the steady-state model prior to 1999. Then the transient model was calibrated according to the groundwater levels at the observation wells distributed throughout the model region (Fig. 7(a)) between 2000 and 2010, while the observed heads in the period 2011–2012 were used in the validation process. A quantitative comparison between observed heads and calculated heads at all monitoring wells and for all observation times is illustrated in Figure 7(b)). The transient model was validated satisfactorily with the value of R² up to 0.988. It is clear that no significant bias appears in the comparison results. For all of the 46 observation wells (Fig. 7(a)) in the validation period, the values of the ME and MAE were statistically −0.127 m and 1.369 m, respectively. The calculated groundwater heads match well with the observation data, with RMSE equal to 1.94 m. The calculated and observed groundwater heads are matched at three typical observation locations, as shown in Figure 7(d)), which indicates that the model is calibrated and validated very well over a long period of time. Accordingly, the numerical model appears to be a robust simulator to delineate the groundwater dynamics in the study area.

3.4 Sensitivity analysis

The main parameters in the model include specified flux (SF), horizontal hydraulic conductivity (HK), precipitation (PRE), irrigation seepage (IS), maximum ET rate (ET_max), groundwater exploitation (GE), specific yield (Sy), and hydraulic conductivity of the stream bed (HKS). Considering the uncertainty of the model and the effect of these parameters on the runoff into Bosten Lake, sensitivity analyses were conducted to test the influence of the all the factors on the flows into Bosten Lake. Here the inflow of Bosten Lake in 2013 was chosen to estimate the sensitivity. All parameters were set to increase or decrease 5% of the value to conduct the sensitivity.

4 Results and discussion

4.1 Groundwater budget analysis

It is essential to analyse the regional groundwater budget to understand water supply and water consumption in the MEYB during the model simulation periods. From the model results, the annual average inflow and outflow in the entire MEYB from 2000 to 2013 are illustrated in Figure 8(a). The main components of the inflow terms of the groundwater budget include precipitation (PRE), irrigation seepage (IS), diversion system seepage (DS), lateral recharge (LR), water quantity from storage (WFS), and river seepage (RS); while the outflow terms include drainage (DRN), evapotranspiration (ET), the water flow into Bosten Lake (WIB) and groundwater exploitation (GE). Table 2 presents all the water budget components in the groundwater flow model. In general, the main components of inflow are the LR and RS, which make up 64.58% of the total inflow. Specifically, 38.22% of inflow comes from LR and 26.36% from RS. This means LR and RS could significantly affect the total inflow and the water budget. For the outflow, the water amount leaving the MEYB, 51.64% of the total amount is through ET and 37.26% is through GE. According to the groundwater budget, an obvious deficit appears (see Table 2) where the ET is much greater than the PRE, implying that in the MEYB the climate affects the groundwater system mainly through ET rather than PRE. Furthermore, there is an increasing trend in the water pumping quantity from storage (WFS), so storage is declining, which indicates that there is unsustainable groundwater exploitation in the MEYB. The groundwater pumping amount follows an annual cycle, which increases during the crop irrigation period and decreases at other times. It is noteworthy that after 2005 the aquifer storage decreased at a faster rate due to the increased pumping rate. Figure 8(b) shows the changes in the main components of inflow and outflow in the groundwater budget over time. It can be seen that LR, RS and GE increase continuously from 2000 to 2013. Increasing GE leads to more water from the river flowing into the aquifer (i.e. an increase in RS). Also, increasing GE results in quickening the velocity of groundwater flow normal to the lateral boundary, eventually representing the increase in LR.

Table 2. Overall water budget of the simulation period 2000–2013 (10⁸ m³/year). PRE: precipitation; IS: irrigation seepage; DS: diversion system seepage; LR: lateral recharge; WFS: water quantity from storage; RS: river seepage; DRN: drainage; ET: evapotranspiration; WIB: water flow into Bosten Lake; GE: groundwater exploitation.

	Inflow							Outflow
Time	PRE	IS	DS	LR	WFS	RS	Total inflow	DRN	ET	WIB	GE	Total outflow
2000	0.04	1.27	2.96	3.00	0.03	1.05	8.49	1.50	4.69	0.31	1.64	8.14
2001	0.07	1.23	3.30	3.02	0.21	1.36	9.20	1.57	5.70	0.20	1.79	9.25
2002	0.03	1.21	3.20	3.24	0.01	1.62	9.30	1.53	5.43	0.26	1.88	9.10
2003	0.06	1.42	3.08	2.96	0.02	1.37	8.91	1.69	4.71	0.30	1.95	8.65
2004	0.03	1.65	2.17	3.61	0.69	1.43	9.58	1.51	5.21	0.24	2.53	9.49
2005	0.01	1.42	1.48	3.44	1.81	2.17	10.03	1.40	5.89	0.18	2.96	10.43
2006	0.05	1.48	1.23	3.21	1.61	2.87	10.26	1.17	5.99	0.30	3.14	10.59
2007	0.03	1.49	1.10	4.82	0.73	2.65	11.02	1.39	5.88	0.19	3.67	11.14
2008	0.09	1.52	1.21	4.63	0.65	3.33	11.03	0.69	6.38	0.07	4.32	11.47
2009	0.09	1.53	1.20	5.94	0.54	3.54	12.43	0.52	6.58	0.17	5.39	12.66
2010	0.00	1.53	1.10	6.00	1.22	4.20	13.84	0.53	7.13	0.25	6.15	14.06
2011	0.04	1.62	1.07	6.32	1.04	5.30	15.29	0.58	7.11	0.18	7.67	15.54
2012	0.07	1.71	0.99	6.28	1.52	6.23	16.80	0.54	7.23	0.21	8.76	16.74
2013	0.01	1.92	1.03	6.34	1.40	6.17	16.86	0.56	6.68	0.15	9.20	16.59

Figure 8. (a) Annual mean groundwater budget for MEYB (10⁸ m³/year); and (b) variation of the main components of inflow and outflow in the groundwater budget.

Correlations calculated using Equation (9) between all components of the groundwater budget are listed in Table 3. There are several highly positive correlations (R > 0.9) between RS–LR, DRN–LR, GE–LR and RS–GE. These highly positive correlations imply that the RS is mainly related to LR and GE, while LR is slightly increasing due to groundwater exploitation (GE). Due to the increased groundwater pumping rate in the study region, the hydraulic connection between river and aquifer system has been changed so that more and more surface water leaks into the groundwater system across the interface along the river bank. Therefore, it can be concluded that the increasing trend of RS absolutely occurs due to the intensive groundwater exploitation in the MEYB.

Table 3. Coefficients of correlation (R) between all components in the groundwater budget. See Table 2 for explanation of abbreviations. Bold values indicate highly positive correlations (R > 0.9).

Item	PRE	IS	DS	LR	WFS	RS	DRN	ET	WIB	GE
PRE	1
IS	–0.140	1
DS	0.046	–0.735	1
LR	–0.032	0.747	–0.777	1
WFS	–0.288	0.606	–0.784	0.454	1
RS	–0.043	0.803	–0.776	0.914	0.629	1
DRN	–0.128	–0.675	0.793	–0.919	–0.507	–0.894	1
ET	0.038	0.600	–0.806	0.889	0.642	0.895	–0.916	1
WIB	–0.286	−0.411	0.522	−0.508	–0.265	–0.453	0.529	–0.539	1
GE	–0.078	0.835	–0.755	0.941	0.589	0.990	–0.891	0.873	–0.455	1

Considering the water scarcity in the arid region, the sustainable exploitation and protection of groundwater in the MEYB is especially important. Figure 7(c) shows the spatial distribution of hydraulic head in 2013 based on the calibrated groundwater flow model. Compared with head contours in Figure 7(a)), some cones of depression occurred in the eastern area of the MEYB due to increasing groundwater exploitation. The groundwater within 2 km of the river bank is obviously recharged by surface water in the lower reaches of the Kaidu River.

The previous analysis and discussion reveal that, in this study area, the dynamic flow system has varied significantly since 1999 in response to groundwater exploitation and interactions between surface water and groundwater. Due to the dry climate in the MEYB, the amount of recharge from PR is very limited and even negligible. The ET is the dominant outflow, although it only increased slightly from 1999 to 2013. As shown in Table 2, the modelling results clearly indicate that the current utilization of the groundwater resources in the MEYB is unsustainable, as intensive pumping makes the RS increase dramatically, leading to a decrease in the runoff from the Kaidu River into Bosten Lake. Due to intensive pumping activities without effective supervision in the whole study region, the water resource shortage will become more serious in the future, causing harm to the environment of Bosten Lake. The salinity of Bosten Lake will keep increasing, having a more serious effect on the water supply source of the Tarim lower reaches, which will produce adverse impacts on the ecological environment in the Tarim basin and make the development of society and the economy unsustainable.

The water budget analysis results indicate that more attention should be paid to sustainable water management in the MEYB. To achieve sustainable development in the MEYB, it is necessary to consider many factors related to the water budget. The ET is the most important consumption of groundwater, so efficient reduction of ET might be the most feasible approach for dealing with water scarcity in the MEYB. Meanwhile, in this arid area of the MEYB, surface water and groundwater exchange frequently, so the groundwater is easily influenced by water use and hydrological conditions. It is shown that pumping water takes up about 37.26% of the total outflow from the groundwater system and has become the main reason for the groundwater resource decrease and river seepage increase in the MEYB, so pumping rates and well locations are the two most important factors that need to be optimized in the entire MEYB.

4.2 Interactions between surface water and groundwater

Surface water interacts with groundwater based on their hydraulic connections. To evaluate river seepage in detail, the whole river channel of the Kaidu and Wulasitai rivers in the study region was divided into six different segments (see Fig. 2(b)), as simulated by the SFR2 package (Niswonger and Prudic 2005) of MODFLOW-NWT. The individual net seepages along the six river segments during the simulation period 2000–2013 are shown in Table 4 and Figure 9. During the entire simulation period, the river loses its surface water along Segments I–V, while it mainly gains its recharge from the aquifer along Segment VI between 2000 and 2009. In general, the total average leakage of the Kaidu River is 3.12 × 10⁸ m³/year, most of which is through Segment I, taking up to 53.19% of the total leakage. Figure 9 further illustrates the intensively increased river leakage mainly in Segment I. Compared with other segments (see Table 4), Segment I plays the most important role in the interaction between surface water and groundwater in the study region. However, the tendency of the seepage in Segment I to increase indicates that the water exchange is unsustainable. As may be seen in Figures 8(b) and 9, groundwater exploitation has the same variation trend as river seepage in Segment I, which implies that groundwater exploitation leads to the river seepage increase in Segment I.

Table 4. Interaction between different river segments and groundwater system. A positive value represents the aquifer gaining water from the river; and a negative value represents the aquifer losing water to the river.

	Exchange water quantity (10^8 m^3/year)
Year	Segment I	Segment II	Segment III	Segment IV	Segment V	Segment VI	Total exchange quantity (10^8 m^3/year)
2000	0.70	0.10	0.13	0.10	0.21	−0.18	1.05
2001	0.79	0.12	0.18	0.14	0.23	−0.11	1.36
2002	0.71	0.15	0.25	0.13	0.37	0.01	1.62
2003	0.64	−0.03	0.27	0.19	0.33	−0.03	1.37
2004	0.86	−0.06	−0.01	0.27	0.38	−0.01	1.43
2005	1.01	0.18	0.19	0.22	0.50	0.06	2.17
2006	1.40	0.37	0.13	0.35	0.68	−0.04	2.87
2007	1.05	0.65	0.23	0.31	0.40	0.00	2.65
2008	1.31	0.88	0.16	0.22	0.80	−0.05	3.33
2009	1.68	0.79	0.20	0.25	0.66	−0.04	3.54
2010	1.73	1.26	0.21	0.29	0.70	0.01	4.20
2011	3.31	0.73	0.26	0.36	0.65	−0.01	5.30
2012	4.27	0.72	0.33	0.21	0.67	0.03	6.23
2013	3.77	0.89	0.30	0.33	0.89	−0.01	6.17

Figure 9. Variation in leakage for different river segments over time (2000–2013). Positive values represent the aquifer gaining water from the river and negative values represent the aquifer losing water to the river.

The water levels in the first hydro-junction, Yanqi hydrological station and Bosten Lake during the observation period are shown in Figure 10. The river level in the first hydro-junction remains nearly unchanged with the same degree of yearly fluctuation, while that in the Yanqi hydrological station decreases continuously. According to the interactions between different segments and groundwater systems shown in Table 4 and Figure 9, it can be inferred that increasing seepage in Segments I and II led to the decrease of the river level in Yanqi station. Consequently, increasing seepage in Segments I–V makes the water level of Bosten Lake decline.

Figure 10. Changing water levels observed at the first hydro-junction (top), and in Yanqi hydrological station and Bosten Lake (bottom).

The SFR2 package uses the recharge and discharge of all water items for a river segment and computes flow across the riverbed based on stream depth calculated at the midpoint of a reach. All recharge and discharge items and runoff in the Yanqi hydrological station and Baolangsumu hydro-junction are calculated using the SFR2 package. The flow rates into Bosten Lake can be obtained directly using the SFR2 package. The hydrographs for Yanqi hydrological station and Baolangsumu hydro-junction on the Kaidu River are illustrated in Figure 11, with the measured runoff also plotted for comparison purposes. The performance criteria for evaluating the difference between the computed results and the measurements are listed in Table 5. For example, compared with actual observed runoff in Yanqi hydrological station (Fig. 11(a)), average relative error (RE) between the observed and calculated runoff is 15.92%, less than the upper limit of 20% required by the standard for hydrological information and hydrological forecasting (GB/T 22,482–2008 2008), which indicates this is an acceptable result. The calculation errors for Baolangsumu hydro-junction are much smaller than those for Yanqi hydrological station: the qualification rate (QR) of the calculated runoff at Yanqi is 0.667, while that for Baolangsumu hydro-junction is up to 0.833 (Table 5). This result shows that calculation results for runoff at Yanqi hydrological station and Baolangsumu hydro-junction are both accurate.

Table 5. Error analysis of the calculated runoff. RE: relative error; ME: mean error; MAE: mean absolute error; RMSE: root mean square error; E: Nash-Sutcliffe efficiency; R: correlation coefficient; QR: qualification rate.

Station	RE (%)	ME (10^8 m^3/year)	MAE (10^8m^3/year)	RMSE (10^8m^3/year)	E (-)	R (-)	QR (-)
Yanqi	15.92	3.01	3.58	3.87	0.71	0. 942	0.667
Baolangsumu	10.27	1.19	2.17	2.73	0.83	0.974	0.833

Figure 11. Comparison of observed and computed runoff of Kaidu River at (a) Yanqi hydrological station and (b) Baolangsumu hydro-junction.

The high-level agreement between the calculated and measured runoff over time for both Yanqi hydrological station and Baolangsumu hydro-junction shows the feasibility and applicability of the methodology used in this study to evaluate the interactions between surface water and groundwater. Based on the analysis of river seepage (Table 4) and water balance from the SFR2 package, the flow rates of the Kaidu River (Segments I–V) to Bosten Lake are illustrated in Figure 12. Thus a declining trend is the common feature of the runoff at Yanqi hydrological station and Baolangsumu hydro-junction (Fig. 13) and the flow rates into Bosten Lake (Fig. 12). Specifically, after 2010, the flow rates into Bosten Lake were decreasing with a faster declining trend. Again, taking the groundwater budget (see Table 2 and Fig. 8) and the river seepage analysis (Table 4 and Fig. 9) into consideration, it is obvious that intensively increasing groundwater exploitation and river seepage has caused the steady decrease of flow rates into Bosten Lake. This indicates that groundwater exploitation plays the most important role in the water level decline of Bosten Lake, implying that the current utilization of the groundwater resources in the MEYB is unsustainable.

Figure 12. Flow rates from Kaidu River to Bosten Lake in the period 2000–2013.

Figure 13. Sensitivity analyses of inflow to Bosten Lake with respect to different model input parameters.

4.3 Uncertainty analysis and some countermeasures

The results of sensitivity analyses are listed in Table 6 and illustrated in Figure 12. From the results one can tell that the hydraulic conductivity of the stream bed is horizontal, which means that the aquifer parameters have the most obvious effects on the inflow of Bosten Lake. Besides, human activities in the study region also influence the flow rates into Bosten Lake significantly. Thus the inflow into Bosten Lake is very sensitive to groundwater exploitation and irrigation infiltration, indicating that increasing these activities will affect the river seepage and play a key role in reducing the inflow into Bosten Lake, which has decreased from 38.68 × 10⁸ m³ in 2000 to 15.37 × 10⁸ m³ in 2013 (Fig. 13).

Table 6. Results of the parameter sensitivity analysis. SF: specified flux; HK: horizontal hydraulic conductivity; PRE: precipitation; IS: irrigation seepage; ET_max: maximum ET rate; GE: groundwater exploitation; Sy: specific yield; HKS: hydraulic conductivity of stream bed. The ranking of sensitivity is relative: 1 indicates the highest sensitivity and 8 the lowest sensitivity.

	+0.05% variation of parameter		–0.05% variation of parameter
Parameter	Sensitivity (10^8 m^3/year)	Rank	Sensitivity (10^8 m^3/year)	Rank
SF	0.000256	8	−0.008910	8
HK	0.014338	3	−0.014353	3
PRE	0.001310	7	−0.009530	7
IS	0.014527	2	−0.011850	6
ETmax	−0.010106	5	0.012080	5
GE	−0.013210	4	0.016560	1
Sy	−0.004705	6	0.013208	4
HKS	−0.015360	1	0.015164	2

Analysis of the water balance for the whole channel of the Kaidu River (Table 4 and Fig. 9) indicates that the annual water leakage mostly occurs in Segment I of the river and this has is increased from 0.70 × 10⁸ m³ in 2000 to 3.77 × 10⁸ m³ in 2013. Therefore, it is necessary to reduce the river seepage of Segment I to protect the water supply and environment of Bosten Lake and the Tarim lower reaches. As the main cause for the river seepage increase, groundwater exploitation reached 9.20 × 10⁸ m³/year in 2013 (Table 2), which indicates a rise of 460% in 14 years. The pumping rate of groundwater exploitation and the well locations in the field are the two most important factors that need to be optimized to reduce the river seepage for the entire MEYB. In addition, in allusion to the large amount of evapotranspiration, with a high value of 6.68 × 10⁸ m³/year in 2013 (Table 2), some measures to protect and improve the regional environment, such as adopting water-saving technologies, changing the crop rotation type, adjusting irrigation and surface water diversion schemes, would be beneficial for relieving the pressure of water shortage in the MEYB and increasing the inflows into Bosten Lake.

5 Conclusions

The Mid-Eastern Yanqi Basin (MEYB) is a typical arid inland basin in Northwestern China, an important location in the Silk Road Economic Belt, where the interaction between surface water and groundwater is greatly changed by groundwater overexploitation. Groundwater and surface water measurement data and a developed flow model based on MODFLOW-NWT have provided valuable insight into the interaction between surface water and groundwater in the basin. The developed model, at the regional spatial scale and long temporal scale, is a useful method for the exploration of sustainable groundwater–surface water resources in the MEYB. According to the results of the calibrated groundwater flow model, groundwater levels have declined significantly in some areas and formed hydraulic depression cones between 2000 and 2013. After 2005 the aquifer storage decreased at a faster rate due to increased groundwater pumping, implying unsustainable groundwater overexploitation in some local areas. Strong positive correlations (>0.9) were obtained between: river seepage and lateral recharge (RS–LR), drainage and lateral recharge (DRN–LR), groundwater exploitation and lateral recharge (GE–LR) and river seepage and groundwater exploitation (RS–GE). Besides, the variations in groundwater exploitation and river seepage have the same trends, which suggests that the river seepage increases to compensate for the deficit in water resources caused by groundwater exploitation.

To evaluate river seepage in different parts of the MEYB, both the Kaidu and Wulasitai rivers in the study region were divided into six segments (Fig. 2(b)) and the water balance was calculated directly using the SFR2 package. The calculated runoff amounts in Yanqi hydrological station and Baolangsumu hydro-junction were well qualified against the actual runoff measurements. As shown in Table 4 and Figure 9, the entire river seepage increased from 1.05 × 10⁸ m³/year in 2000 to 6.17 × 10⁸ m³/year in 2013, with the increased river leakage occurring mainly in Segment I. Furthermore, the calculated runoff rates in Yanqi hydrological station and Baolangsumu hydro-junction are also well matched with the measured values. Declining trends may be seen for runoff at Yanqi hydrological station (Fig. 11(a)) and Baolangsumu hydro-junction (Fig. 11(b)), and the runoff into Bosten Lake (Fig. 12). Based on the water budget in the whole region, it is apparent that groundwater overexploitation and river seepage have led to the steady decreasing trend for the runoff into Bosten Lake, indicating groundwater exploitation is the key factor for the water level decline of Bosten Lake.

Both water budget analyses of the groundwater flow model and the water balance for different Kaidu River segments suggest that the current utilization of the groundwater resources in the MEYB is unsustainable. For the groundwater system, evapotranspiration is an important water consumer, so it is direct and feasible to find approaches to reduce ET in the MEYB to alleviate the problem of water scarcity. With over 53% of the average total seepage of the Kaidu River (3.12 × 10⁸ m³/year) being from Segment I, a reduction in the river seepage at this segment is very important and will be a practical way to protect the water supply and environment of Bosten Lake and the Tarim lower reaches. Among the measures that can be taken to prevent groundwater exploitation in the MEYB, the pumping rate and locations of well fields are the two most important factors to be optimized in the entire MEYB to reduce the river seepage. Sensitivity analysis also indicates that human activities such as groundwater exploitation and agricultural irrigation greatly affect the runoff of the Kaidu River into Bosten Lake, which suggests the major factor for water-level depletion is groundwater and surface water exploitation. To make water consumption in agriculture more efficient, it is necessary to adopt water-saving technologies, change the crop rotation type, adjust irrigation and implement surface water diversion schemes to improve efficiency of water use. Combining all four strategies may be the most effective method to solve the water resources shortage and environmental problems in the MEYB. These management strategies are essential for solving problems and improving sustainable development in the MEYB. Simultaneously, they are also valuable and effective approaches for developing sustainable management strategies for water resources in other hotspots with depletion of the groundwater−surface water resource. The regional-scale and long time-scale model is very important for studying the interactions between groundwater and surface water systems.

Acknowledgements

The authors appreciate the support from the Nanjing Hydraulic Research Institute and Xinjiang Tarim River Basin Management Bureau, Xinjiang Uyghur Autonomous Region of China. The numerical calculations in this paper were implemented on the IBM Blade cluster system in the High Performance Computing Center of Nanjing University. The authors have benefited from communications with Dr Xiaoqing Shi and Dr Xing Min in establishing numerical modelling. The authors are also profoundly grateful to the Co-editor Professor Demetris Koutsoyiannis, Associate Editor Dr Tomasz Okruszko, and the anonymous reviewer whose invaluable suggestions and constructive comments helped to improve the manuscript significantly.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The authors gratefully acknowledge the financial support of the National Key Research and Development Plan of China (2016YFC0402800), the National Natural Science Foundation of China (U1503282, 41772254 and 41372235), and the China Postdoctoral Science Foundation (2017M622905).