Decomposition technique for contributions to groundwater heads from inside and outside of an arbitrary boundary: application to Guantao County, North China Plain

To assess the efficiency of the groundwater management of an administrative unit, we propose decomposing the groundwater head changes within an administrative unit into inside and outside contributions by using numerical models. Guantao County of Hebei Province, China, serves as an example to demonstrate the decomposition technique. The groundwater flow model of Guantao was constructed using observed heads as prescribed head boundary conditions. The model was coupled with Hydrus 1-D to calculate the groundwater recharge distribution in time reflecting the delay and damping effects of the soil column on seepage at the surface. The model was calibrated by adjusting parameters such as hydraulic conductivities, recharge infiltration ratios and specific yields. The calibrated parameters are then used in a large model with a boundary at a large distance from Guantao administrative boundary to determine the groundwater head changes due to inside drivers. The differences of the two models on the Guantao boundary serve as the specified head values on the boundary for a small-scale model, which is used to calculate the groundwater head imposed by outside drivers. To eliminate inconsistencies caused by the different types of boundary conditions of large and small models, the groundwater head changes due to inside drivers must be updated. The results indicate that the groundwater head changes in the center and south of Guantao County are influenced equally by both inside and outside contributions, while in the north outside contributions have the stronger impact. On average, 48.5 % of groundwater head changes in the whole of Guantao County is influenced by inside contributions, while 51.5 % is due to outside contributions. The sensitivity analysis shows that the groundwater head changes and their decomposition are much more sensitive to infiltration ratios than to the aquifer parameters. The parameters within Guantao have a certain influence on the net groundwater head changes, while the parameters outside of Guantao only have an influence on the decomposition.

1 Introduction

The natural shallow groundwater flow field is usually determined by infiltration from precipitation and streams, by discharge to streams or springs and by phreatic evapotranspiration. With the development of human activities, the natural groundwater levels are modified by land use change, crop planting and irrigation, groundwater pumping and more. Over the last few decades, the groundwater resources worldwide have been intensively used for household, industry and agriculture purposes (Zektser and Everett, 2004). It is estimated that of the approximately 1000 km³ of groundwater pumped worldwide per year, about one-quarter is not replenished by recharge anymore (WWAP, 2012; Wada et al., 2012). The major user of groundwater is irrigated agriculture. Especially in the North China Plain (NCP), groundwater became the major water resource for agriculture, accounting for more than 70 % of the total water supply (Wu, 2013). The continuous unsustainable abstraction of groundwater caused a dramatic groundwater-level decline: by 18 m in the piedmont region in 25 years (Chen et al., 2003), by 5–50 m in the shallow aquifer and by 10–100 m in the deep aquifer over the last 50 years (Wu, 2013). This severe groundwater-level decline causes land subsidence, drying up of the lower reaches of rivers and streams, salt-water intrusion at the Hebei coast, and municipal water shortages (Wu, 2013; Liu et al., 2001; Han et al., 2011). It also increases the cost of pumping groundwater. To support decision makers in putting forward rational water management strategies in such areas and alleviating groundwater decline or even reaching some recovery of groundwater levels, it is useful to separate impacts on groundwater head changes caused by different drivers. Hydrological models have increasingly been applied in assessing the influences on groundwater levels due to individual drivers such as precipitation, changes in cropping system, and irrigation practice, water imports and others (Shu et al., 2012; Cao et al., 2013; Li, 2013; Iwasaki et al., 2014). For a small region within a larger aquifer system it might be of interest to identify the influences of drivers within and outside of the region.

To build a groundwater model of an administrative unit is difficult as usually natural model boundaries do not coincide with political boundaries. However, we often want to assess the efficiency of groundwater management of an administrative unit. The groundwater head changes within an administrative unit are influenced not only by the groundwater head changes induced by inside drivers, but also by the groundwater head changes induced by outside drivers. Monitoring groundwater levels within the administrative unit often cannot directly and faithfully reflect the effects of water management measures implemented in that unit. In other words, the monitored groundwater head changes are caused by both inside and outside contributions. Therefore, we propose to decompose the groundwater head changes over a given time period within an administrative unit into inside and outside contributions to assess the efficiency of the unit's management of groundwater resources. To do so three groundwater models are required:

1. a model of the administrative unit with the political boundary implemented as a prescribed head boundary using observed heads,

2. a larger model containing this model, which computes only the propagation of inside changes to the outside and which may extend to physical model boundaries, and

3. a model of the administrative unit with prescribed heads obtained by subtraction from the previous two models on the political boundary, which is used to determine the groundwater head contributions controlled by outside drivers only.

Guantao County is located in the NCP and is delimited by a political administrative boundary. To test aquifer rehabilitation measures, Guantao is chosen as a pilot region for research. When applying remedial measures in Guantao, e.g., by decreasing water abstraction from the aquifer via suitable measures, the groundwater level in Guantao County could increase, decrease or show no change with time depending on the superimposed influence of the outside groundwater flow field. This impedes the assessment of the effects and the efficiency of Guantao water management strategies. The solution to fair assessment lies in decomposing the groundwater head changes in Guantao into two parts: one part is the head change due to the influence of recharge and discharge in Guantao (inside contribution) and the other is the head change contribution determined by outside drivers, i.e., the groundwater flow field without considering source and sink terms in Guantao (outside contribution).

Numerous studies were carried out in the NCP to understand the large-scale groundwater dynamics employing MODFLOW (Lu et al., 2011; Zhou et al., 2012; Qin et al., 2013). For local management purposes these models are far too crude. But groundwater simulations in smaller administrative units of the NCP have been quite rare up to now. In our study, a groundwater flow model for Guantao is constructed using MODFLOW USG under the pre- and post-processor Visual MODFLOW Flex 2014.1 (VMF) developed by Waterloo Hydrogeologic (Panday et al., 2013).

The objective of this research is to take Guantao as a case study to demonstrate how to decompose the groundwater head changes into inside and outside contributions using numerical models. First the general background and available data are introduced. Then the mathematical method of decomposing the groundwater head changes is described in detail. Following that, a numerical model for the Guantao aquifer system is developed and calibrated under both steady-state and time-varying conditions. Then, both a large-scale model and a small-scale numerical model are developed to separately determine the groundwater head changes due to inside and outside contributions. Finally, the sensitivity of the groundwater head changes and their decomposition with respect to parameters and local management measures are discussed.

2 Study area and data

2.1 Study area

The NCP is bordered by the Tai Hang Mountains in the west, the Yellow River in the south, the Bohai Sea in the east and the Yan Shan Mountains in the north, and has an area of altogether 140 000 km² (Fig. 1). Guantao County is located in the southern NCP, between longitudes 115^∘06^′ and 115^∘28^′ E and latitudes 26^∘27^′ and 36^∘47^′ N, covering an area of 456 km² (Fig. 1). There are eight townships and 277 villages under the county's administration, with a population of 340 000 (GSB, 2002–2013). From the piedmont in the west to the coast in the east, the NCP is divided into four geomorphological units, with Guantao lying in the alluvial plain zone (Wu et al., 1996). The topography of Guantao is quite flat, with an elevation of 45 m a.s.l. (above sea level) in the southwest, which is only slightly higher than the value of 36 m a.s.l. in the northeast.

Figure 1 Schematic map of the NCP (a) and administrative map of Guantao (b).

Guantao has the characteristics of a warm continental monsoon climate: hot and rainy in summer, and cold and dry in winter with an average annual precipitation of 532 mm, an average annual potential evaporation of 1516 mm and an annual mean temperature of 13.4 ^∘C (GEF Project, 2009). About 80 % of the annual precipitation occurs between June and September. The Weiyun River is a seasonal river, which dries up outside of the monsoon season. It flows along Guantao's eastern boundary between June and October with a flow of 0.4×10⁹ m³ a⁻¹ averaged over the period 2001–2012 (data collected by GIWP). The Weixi channel is the main irrigation channel diverting surface water from the Weiyun River to croplands. The total cultivated area is around 300 km² (data from 2013) and the main crops are winter wheat, summer corn and cotton. Winter wheat is sown in October and harvested in June of the following year. After wheat harvest, the summer corn is immediately planted on the same area and harvested in September. The cotton growing period is between April and September. The annual diversion in the channel has only been around 3.1×10⁶ m³ a⁻¹ (average between 2001 and 2012; data collected by GIWP). During the monsoon season, less irrigation is needed because the concentrated precipitation is sufficient to satisfy crop needs most of the time. The non-monsoon season coincides with the growing season of winter wheat, which means irrigation with groundwater is the only way to support winter wheat growth.

The Quaternary aquifer of Guantao consists of fine sand layers interbedded with clay or silt aquitard layers. Vertically it is divided into shallow, middle and deep layers according to different deposits. Based on previous studies and the recent hydrogeological investigation by CNACG (2015), the middle layer mainly consists of clay with an average thickness of around 120 m. It contains saline water with total dissolved solids reaching 10 000 g m⁻³ in some parts (GEF Project, 2009). The shallow layer is exclusively used for irrigation water supply, while the deep layer is mainly used for domestic and industrial water supply. Still, the deep layer has to provide some irrigation water too due to locally excessive salinity in the shallow aquifer. The groundwater-level differences between deep and shallow aquifers are huge, and the main chemical groundwater characteristics in both layers are totally different (CNACG, 2015). Hence, the middle layer is regarded as an efficient aquitard which practically does not allow a direct hydraulic connection between deep and shallow layers. The shallow aquifer is recharged by precipitation, irrigation backflow and channel and river seepage. It is discharged through pumping wells for irrigation. The deep layer receives only little recharge from upstream at the piedmont, so groundwater levels have been decreasing since pumping for domestic, industrial and irrigation uses started. Long-term intensive groundwater pumping caused many environmental problems such as groundwater depletion, land subsidence, and seawater intrusion. Stopping all abstraction from the deep layer is the only way to let groundwater levels in the deep aquifer recover. The political decision has been taken to stop almost all pumping from the deep aquifer and replace this water with water from the South–North Water Transfer scheme. By 2017 about 30 % of domestic and industrial use water was replaced, while agriculture is still pumping about 10 % of its groundwater use from the deep aquifer (local communication). In this study, the shallow aquifer is the main concern, and the modeling effort comprises only this layer.

2.2 Data

All data required to calibrate a numerical groundwater model for the time between 2003 and 2012 are available. Groundwater heads (provided by GIWP and Handan DWR) are observed at three different frequencies. There are 4 observation points with monthly measurements, 11 observation points with two data records each year and 29 observation points with four measurements each year. Daily precipitation data have been purchased from Guantao Meteorological Bureau (yearly data of precipitation between March and October shown in Fig. 2). Annual channel diversions for irrigation and monthly runoff of the Weiyun River as well as officially recorded annual pumping rates have been collected from the Water Resources Bulletin of Guantao for the period from 2003 to 2012 (Fig. 2). The pumping rate for irrigation is highly dependent on precipitation in the NCP. Pumping rates are generally less in years with higher precipitation. Hence, the data point for 2003 showing a combination of higher precipitation with a larger pumping rate is questionable. The sudden significant decrease in the reported pumping rate after 2006 is also questionable as there were hardly any changes in cropping area, crop types and irrigation methods. Therefore, the pumping rates for these years will be adjusted during manual calibration of the transient numerical model.

Figure 2 Time series of the annual pumping rates in the shallow aquifer shown together with annual precipitation (solid line: reported data; dashed line: values adjusted during calibration).

3 Methods

3.1 Equations

Groundwater heads in Guantao ($h_{\text{Guantao}}(t,x,y$)) result from the superposition of groundwater head changes due to inside drivers ($\mathrm{\Delta }{h}_{\text{Guantao}}(t,x,y$): inside drivers) and groundwater heads determined by outside drivers ($h_{\text{out}}(t,x,y$): outside drivers). The decomposition approach is written as follows:

$$\begin{array}{}\text{(1)}& h_{\text{Guantao}}\left(t,x,y\right)=\mathrm{\Delta }{h}_{\mathrm{Guantao}}\left(t,x,y\right)+h_{\mathrm{out}}\left(t,x,y\right).\end{array}$$

The governing equation for the transient groundwater flow in a phreatic aquifer can be described using a linear partial differential equation (PDE) with a time-independent transmissivity provided the head changes in a predefined period are relatively small compared to the aquifer thickness.

$$\begin{array}{}\text{(2)}& \mathit{\mu }{\displaystyle \frac{\partial h}{\partial t}}=T{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x^{\mathrm{2}}}}+T{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial y^{\mathrm{2}}}}+q^{\mathrm{in}}-q^{\mathrm{out}}\end{array}$$

h refers to the groundwater head (m), T is the transmissivity (m² s⁻¹), μ is the specific yield (–), t is the time (s), and q_in and q_out are source and sink terms, respectively (m² s⁻¹).

The groundwater heads in Guantao can be determined by the governing equation below together with a (measured) specified head boundary along the administrative boundary.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathit{\mu }{\displaystyle \frac{\partial h_{\mathrm{Guantao}}}{\partial t}}& =T{\displaystyle \frac{{\partial }^{\mathrm{2}}{h}_{\mathrm{Guantao}}}{\partial x^{\mathrm{2}}}}+T{\displaystyle \frac{{\partial }^{\mathrm{2}}{h}_{\mathrm{Guantao}}}{\partial y^{\mathrm{2}}}}\\ & +q_{\mathrm{Guantao}}^{in}-q_{\mathrm{Guantao}}^{\mathrm{out}}\end{array}\end{array}$$

$q_{\mathrm{Guantao}}^{\mathrm{in}}$ includes not only the anthropogenic sources (irrigation backflow and channel infiltration), but also the natural recharge by precipitation. River seepage is included in the prescribed head boundary as the boundary coincides with the Weiyun River.

Accordingly, the PDE to determine the groundwater head changes due to inside drivers is as follows:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathit{\mu }{\displaystyle \frac{\partial \mathrm{\Delta }{h}_{\mathrm{Guantao}}}{\partial t}}& =T{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{\Delta }{h}_{\mathrm{Guantao}}}{\partial x^{\mathrm{2}}}}+T{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{\Delta }{h}_{\mathrm{Guantao}}}{\partial y^{\mathrm{2}}}}\\ & +q_{\mathrm{Guantao}}^{\mathrm{in}}-q_{\mathrm{Guantao}}^{\mathrm{out}}.\end{array}\end{array}$$

To solve Eq. (4), apart from the knowledge of sources and sinks in Guantao, a proper choice of the model boundary is also important. The ideal boundary is chosen at a large distance from the Guantao County administrative boundary to make sure that the influence of sources and sinks within Guantao on this distant boundary can be assumed to be equal to zero.

The groundwater heads in Guantao caused by outside drivers only can be obtained by numerically solving the following flow equation:

$$\begin{array}{}\text{(5)}& \mathit{\mu }{\displaystyle \frac{\partial h_{\mathrm{out}}}{\partial t}}=T{\displaystyle \frac{{\partial }^{\mathrm{2}}{h}_{\mathrm{out}}}{\partial x^{\mathrm{2}}}}+T{\displaystyle \frac{{\partial }^{\mathrm{2}}{h}_{\mathrm{out}}}{\partial y^{\mathrm{2}}}}+q_{\mathrm{outside}}^{\mathrm{in}}-q_{\mathrm{outside}}^{\mathrm{out}}.\end{array}$$

$q_{\mathrm{outside}}^{\mathrm{in}}/q_{\mathrm{outside}}^{\mathrm{out}}$ are the sorce/sink terms outside of the Guantao area. The model used to solve Eq. (5) is only established for the Guantao area, where the $q_{\mathrm{outside}}^{\mathrm{in}}/q_{\mathrm{outside}}^{\mathrm{out}}$ terms are equal to zero. A specified head boundary condition is defined on the administrative boundary, which can be obtained as the difference between the boundary condition of Eq. (3) and the solution of Eq. (4) on the county boundary. Anyway, when two variables in Eq. (1) are known, the third variable is determined as well.

As mentioned above, to decompose the groundwater heads in Guantao, three numerical models are used in this research. The first is a small-scale numerical model within Guantao's administrative boundary to calculate $h_{\mathrm{Guantao}}(t,x,y)$ (Guantao flow model). It uses measured heads on its (specified head) boundary, which coincides with the county border. The second is a larger-scale numerical model extending to distant natural boundaries to calculate $\mathrm{\Delta }{h}_{\mathrm{Guantao}}(t,x,y)$ (large flow model). Its boundary conditions are zero head change induced by Guantao sources and sinks. The third one is again a small-scale numerical model within Guantao's administrative boundary to calculate $h_{\mathrm{out}}(t,x,y)$ (h_out flow model in Guantao). The specified heads on the boundary are now determined by sources and sinks outside of Guantao only. They can be obtained from the measured heads by subtracting the head changes obtained from the large flow model on the county boundary. The analysis of the decomposed groundwater heads can be carried out at any time and any location within Guantao County. For a predefined time period from t to t+Δt Eq. (1) can be restated as Eq. (6), which says that the groundwater head changes over that period ($\mathit{\delta }{h}_{\mathrm{Guantao}}(\mathrm{\Delta }t,x,y))$ can also be decomposed into head changes due to inside contributions ($\mathit{\delta }\mathrm{\Delta }{h}_{\mathrm{Guantao}}(\mathrm{\Delta }t,x,y))$ and outside contributions ($\mathit{\delta }{h}_{\mathrm{out}}(\mathrm{\Delta }t,x,y))$.

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathit{\delta }{h}_{\mathrm{Guantao}}\left(\mathrm{\Delta }t,x,y\right)& =\mathit{\delta }\mathrm{\Delta }{h}_{\mathrm{Guantao}}\left(\mathrm{\Delta }t,x,y\right)\\ & +\mathit{\delta }{h}_{\mathrm{out}}\left(\mathrm{\Delta }t,x,y\right),\end{array}\end{array}$$

in which

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}\mathit{\delta }{h}_{\mathrm{Guantao}}\left(\mathrm{\Delta }t,x,y\right)& =h_{\mathrm{Guantao}}\left(t+\mathrm{\Delta }t,x,y\right)\\ & -h_{\mathrm{Guantao}}(t,x,y),\end{array}\text{(8)}& {\displaystyle }\begin{array}{rl}\mathit{\delta }\mathrm{\Delta }{h}_{\mathrm{Guantao}}\left(\mathrm{\Delta }t,x,y\right)& =\mathrm{\Delta }{h}_{\mathrm{Guantao}}\left(t+\mathrm{\Delta }t,x,y\right)\\ & -\mathrm{\Delta }{h}_{\mathrm{Guantao}}\left(t,x,y\right),\end{array}\text{(9)}& {\displaystyle }\mathit{\delta }{h}_{\mathrm{out}}\left(\mathrm{\Delta }t,x,y\right)=h_{\mathrm{out}}\left(t+\mathrm{\Delta }t,x,y\right)-h_{\mathrm{out}}\left(t,x,y\right),\end{array}$$

for any point (x,y) within Guantao borders. As we are mainly interested in changes in heads due to both inside and outside contributions, we will in the following use the formulations of Eq. (6) through Eq. (9). Δt refers to the whole simulation period of 10 years.

3.2 Model description

3.2.1 Guantao flow model

The Guantao model boundary is the administrative boundary of Guantao, which is conceptualized as a specified head boundary condition. All available groundwater head observations close to the boundary are made use of by defining the boundary heads by interpolation. The bottom of the aquifer (lower boundary of the model) is defined by the boundary between the shallow aquifer and the aquitard. It is interpolated from the well logs of the hydrological investigation (CNACG, 2015) and is considered impervious.

Precipitation is the main recharge to groundwater in Guantao. Other natural recharge terms include river infiltration of the Weiyun River on the eastern boundary and infiltration from the Weixi channel passing through Guantao from south to north. The Weiyun River infiltration is not explicitly described in the numerical model, but is implicitly considered in the specified head boundary condition. The Weixi channel infiltration is simulated using the river package of MODFLOW. The river level is assumed to be 3 m higher than the deepest river channel bottom level. Guantao is an agricultural county. Apart from the main Weixi channel, there are still numerous smaller sub-channels connected to the Weixi channel to divert and store surface water. Infiltration from these smaller sub-channels is considered part of the areal recharge from irrigation and is simulated as a contribution to the recharge package. This is a common way to deal with the leakage to groundwater from numerous field channels in the irrigated cropland, which can also be found in other similar studies (Rejani et al., 2008; Cao et al., 2013).

Groundwater is the main source of irrigation. This is very common in the NCP (Hu et al., 2010; Zhang et al., 2010). In Guantao, due to lack of surface water, groundwater accounts for over 95 % of total applied irrigation water according to data between 2001 and 2011. Nowadays the traditional conveyance of irrigation water by canals and ditches is replaced by pipe irrigation to increase the conveyance efficiency from wells to the cropland. In the field, ridge irrigation and furrow irrigation are still the primary methods for both groundwater and surface water irrigation. Hence, the irrigation backflow is another important source of groundwater recharge. In the study, the recharge contributions by precipitation, irrigation backflow and canal seepage were added up and implemented in the recharge package. The distribution of agricultural area is obtained from satellite remote sensing images of 2014, by mapping the normalized difference vegetation index (NDVI) in winter and spring (provided by Haijing Wang, hydrosolutions Ltd.). The total agricultural area changes only slightly between 2002 and 2013 according to the statistical year books (GSB, 2002–2013). Hence, the agricultural area is assumed to be constant over the modeled time span.

The depth to groundwater is around 20 m according to the available groundwater head observations, so there is no phreatic evaporation from the aquifer. The discharge is exclusively due to pumping wells extracting groundwater from the aquifer. The water abstraction by pumping is implemented in the well package and the spatial distribution is chosen according to the well locations obtained from the general investigation in 2011 (local communication). The pumping rate of the individual well was determined as a constant fraction of the estimated total annual groundwater abstraction (see Sect. 4.1 below).

The two-dimensional groundwater flow model is constructed using VMF-USG, the unstructured grid version of USGS-MODFLOW. The spatial discretization consists of 141 379 cells. The cells are refined around pumping wells. The top elevation is interpolated from SRTM data downloaded from the website at http://srtm.csi.cgiar.org (last access: 2014). The thickness of the aquifer is between 55 m and 96 m. The spatial distribution of hydraulic conductivity (K) is characterized by zonation of the model domain. Four zones are defined based on pumping tests carried out by CNACG (2015) (Fig. 3). Fourteen groundwater head observations are used to define the fixed head boundary. They divide the boundary into segments with observations at both ends of a segment. The cell values within a segment are linearly interpolated from the observations on both ends. Missing observation data in some stress periods are linearly interpolated from the two closest available data records of the time series. The adjusted annual groundwater abstraction is assigned to monthly stress periods within a year based on irrigation water application. Irrigation is mainly applied outside of the monsoon season. According to local information, winter wheat is primarily irrigated in October, March and April, summer maize in May and June and cotton in April and May, so the groundwater abstraction is set to zero in all other monthly stress periods.

Figure 3 Zonation of hydraulic conductivities and recharge infiltration ratios.

Recharge is the most difficult input to quantify in the simulation. Some researchers suggest that it is a nonlinear function of the water input at the surface (Kendy, 2003, 2004; Tan et al., 2014; Min et al., 2015). To simplify the process, the recharge infiltration ratio (R) is divided into three zones according to the soil map (CNACG, 2015) (Fig. 3). For each zone, the recharge fluxes from precipitation and irrigation are estimated as constant fractions of annual precipitation and irrigation amounts as suggested in other literature on the groundwater system in the NCP (Jia et al., 2002; Liu et al., 2008; Shao et al., 2013). These recharge infiltration ratios will be calibrated in the steady-state model.

However, for a transient calculation, the temporal distribution of recharge also has to be represented adequately to describe its influence on groundwater-level dynamics correctly. In order to better understand the time lag of the groundwater-level response to the water input at the ground surface, the groundwater movement in the unsaturated soil zone was simulated using Hydrus 1-D (Simunek et al., 2013) in this study. The soil column is modeled from the ground surface to a depth of 20 m, which is comparable to the yearly averaged depth to groundwater table in Guantao over the last 10 years. Lu et al. (2011) studied the groundwater recharge at five representative sites in the NCP (two in the piedmont plain, two in the alluvial plain and one in the coastal plain) using Hydrus. Guantao is located in the alluvial plain of the NCP with soils mainly consisting of silt, clay and silty clay. Since no field sampling, monitoring or experiments in the soil zone were carried out within the project, the soil material is assumed to be uniformly distributed in the column. The soil properties are adopted from the two representative sites in the alluvial plain cited in Lu et al. (2011). The van Genuchten parameters chosen as input into Hydrus are the average values from Lu's study (2011): θ_residual=0.08, θ_saturated=0.39, α=0.54 (1 m⁻¹); n=1.98; and k_s=0.65 (m d⁻¹). The simulated soil column is vertically discretized into 100 layers with a spacing of 0.2 m. A self-adjusting numerical time discretization was adopted in the simulation with a minimum time step of 0.0001 d and a maximum time step of 1 d. The flexible time discretization in the simulation ensured convergence of the numerical calculation.

The groundwater recharge is obtained from the calibrated steady-state model as a percentage of total precipitation plus irrigation, which enters the soil column below the root zone. (All other water applied to the soil is assumed to be consumed by evaporation and plant transpiration.) Monthly inputs are computed according to rainfall and irrigation events with this constant recharge ratio. Due to the large depth to groundwater of 20 m and more the water input is delayed and attenuated. The final temporal distribution of the flux at the groundwater table is calculated in the Hydrus simulation. The upper boundary of the column is implemented as an atmospheric boundary to the surface layer, while a free drainage condition is implemented at the lower boundary. The initial distribution of water content is obtained by running the simulation periodically for 10 years to get a relatively steady distribution.

The aim of the model calibration is to minimize the residuals between the observed and computed head values by adjusting model parameters. We start out with a steady-state model. The average behavior in the period from 2003 to 2011 was used as a quasi-steady state. This is feasible, as the average head over all observations at the end of 2002 is almost identical to that at the end of 2011. That means that on average, abstractions must have been balanced by recharge over that period. All inputs required in this case are average data sets between 2003 and 2011. The steady-state model calibration was accomplished using the automated parameter estimation code (PEST) (Doherty, 2003). The transient flow model between 2003 and 2012 is calibrated manually, as only the specific yield values have to be adjusted. The groundwater flow field obtained from the steady-state model is used as the initial condition of the transient model. The transient model time span was divided into monthly stress periods. There are in total 120 stress periods. The whole calibration process is as follows. Initially, the officially reported (but estimated) time series of the pumping rate and inferred recharge ratios are input into Hydrus 1-D to obtain the average groundwater recharge between 2003 and 2011. Hydrus 1-D is run for each township using different seepage water input values according to the collected data. Then the parameter estimation model PEST is applied to calibrate parameters (K and R). The distribution of groundwater recharge from Hydrus is input into the transient model and specific yield values are calibrated manually. The specific yield influences both the amplitudes of heads within the year and their rate of change over the longer term. When the yearly averaged calculated head value at all observation wells is far away from the corresponding observed value, the recorded pumping rate in that year is adjusted first. In this application, the biggest deviation between computed and measured heads is controlled within 1 m. The whole calibration process is implemented iteratively to converge to the best fit to observations.

After model calibration, several statistical indices were used to assess the results of the numerical model. The correlation coefficient R-square is used to measure the correlation between modeled and observed groundwater heads in the steady-state case. Root mean square error (E_RMS) is used to evaluate how closely modeled groundwater heads fit observed data time series and mean absolute error (E_MA) is used to evaluate the average magnitude of errors we can expect from modeled groundwater heads. Both E_RMS and E_MA are used to evaluate the model performance on average. They are defined as follows:

$$\begin{array}{}\text{(10)}& E_{\mathrm{RMS}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^m{\sum }_{t=\mathrm{1}}^n(Y_{i,t}-O_{i,t}{)}^{\mathrm{2}}}{m\times n}}},\end{array}$$

$$\begin{array}{}\text{(11)}& E_{\mathrm{ME}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^m{\sum }_{t=\mathrm{1}}^n\mathrm{|}{Y}_{i,t}-O_{i,t}\mathrm{|}}{m\times n}},\end{array}$$

where Y_i,t is the ith modeled groundwater head at time t; O_i,t is the corresponding observed value at the same time step; m is the number of observation wells; and n is the number of time steps.

Table 1 The automatically calibrated values of hydraulic conductivities and recharge infiltration ratios and their respective sensitivities to observed heads.

The zonal values of hydraulic conductivities and recharge infiltration ratios obtained from PEST in the calibration of the steady-state model are presented in Table 1. The values of hydraulic conductivities indicate a decreasing trend from south to north. The relatively less permeable zone between K1 and K3 is consistent with the local hydrogeological knowledge. The distribution of the recharge infiltration ratios shows that the values tend to decrease from east to west. Tan et al. (2014) show that in the NCP the average recharge ratio from tracer experiments is 0.16, ranging from 0.11 to 0.21, which is quite close to the results from another tracer test done in the NCP by Hu et al. (2016). The average value from numerical simulations is 0.14, ranging from 0.08 to 0.2. Kendy (2003) points out that the recharge ratio obtained from soil water balances at 16 sites in the NCP ranges from 0.09 to 0.24. Kendy (2004) concluded that the recharge ratio in the NCP varied between 0.2 and 0.6 in numerical model simulations. Min et al. (2015) showed that the recharge ratio varied between 0.07 and 0.53 in a 30-year Hydrus simulation. Lu (2011) concluded that the recharge ratio is 0.18 in the alluvial plain according to Hydrus simulations. The former studies reported a wide range of values for the recharge infiltration ratio in the NCP, with an average value around 0.23. This further proves that the calibrated recharge ratios in our study are reasonable and acceptable. During the auto-calibration process, R3 tends to be extremely small, with a value even less than 1 %. To avoid this unrealistic number, R3 is assumed to be fixed in the calibration at a value of 0.09, which is the average of the minimum recharge ratios obtained from the previous studies in the NCP.

The sensitivities of different parameters to observed heads are shown in Table 1. The parameter (composite) sensitivities are automatically computed by PEST and saved in the file *.sen. They are calculated according to the weighted Jacobian matrix and the number of observations (Doherty, 2003). The table indicates that compared to hydraulic conductivities, the recharge infiltration ratios have a much larger influence on groundwater heads. The sensitivity of hydraulic conductivity is higher in the south than in the north. The most sensitive recharge infiltration ratio is found in the east. The correlations between different parameters are presented in Table 2. Correlation coefficients with a modulus above 0.5 indicate that there is a non-negligible correlation among the respective parameters. The highest (anti-)correlation coefficient between two zonal hydraulic conductivities is −0.58 (between K2 and K3). This shows that their relative values cannot be uniquely determined by calibration. The recharge ratios R1 and R2 are also correlated, with a correlation coefficient of −0.57, which means that the sum of the recharge in zones 1 and 2 is more important than the individual values. The R-square of the comparison between simulated groundwater heads and observations is around 0.83, which implies that 83 % of the variance of observations could be explained by the numerical model. Although the fit to the observed data is very good considering the data situation, it has to be noted that the result is still uncertain. An increased pumping rate with a simultaneously increased recharge ratio will lead to a similar result. For a unique result, accurate pumping rates are required.

The recharge fluxes at different soil depths for one of the eight townships are shown in Fig. 4. They capture the important characteristics of the soil water flux at different depths. The response time becomes longer with increasing depth to groundwater, while the peak flux becomes smaller with a more even distribution over the year. Similar results can be observed in other areas of the NCP (Li et al., 2017; Min et al., 2017). Although the water input on the soil surface is quite different from year to year, few peaks are observed when the soil depth is larger than 10 m, due to damping.

Table 2 Correlation coefficients among different parameters.

Figure 4 Calculated monthly recharge flux at different soil depths.

The peaks vanish more quickly with increasing soil depth in dry years than in wet years. Despite the averaging, the bottom flux distribution over a year is different from year to year. Even at 20 m depth the biggest difference between the highest and lowest fluxes is still around 30 % of the lowest flux. The groundwater recharge is a complex nonlinear process which mixes all the irrigation and precipitation events occurring at different times. Through the analysis of results from Hydrus, it is obvious that the best way to accurately describe the groundwater recharge in the groundwater model is to couple the model of the unsaturated zone with the groundwater model. The recharge fluxes for other townships demonstrate a similar distribution.

Figure 5 Calculated and observed monthly time series of groundwater heads at four available long-term observation wells.

The specific yield has the same zonation as the hydraulic conductivity. Four zones of specific yield, denoted by S_y1, S_y2, S_y3 and S_y4, are calibrated manually to reflect the correct head amplitudes. The specific yields decrease from south to north with the respective values of 0.05 (S_y1), 0.03 (S_y2), 0.02 (S_y3) and 0.01 (S_y3). The manually calibrated values are to some extent subjective. The adjusted pumping rates are shown in Fig. 2. They are consistent with the common experience that less groundwater is pumped in wet years, while the water demand increases in dry years.

Figure 6 Discretization map of the large model.

Figure 7 Parameter distributions map for the NCP (modified based on figures in the literature; Liu et al., 2011).

The goodness of fit between modeled and observed groundwater heads is presented in Fig. 5. Some computed heads underestimate observations and others overestimate observations. There is no general bias. The seasonal pattern of groundwater-level dynamics is reproduced. The groundwater level starts to decrease from March due to irrigation and reaches its lowest value in June or July. From June or July on, groundwater levels start to recover and increase in response to termination of pumping and the delayed recharge by precipitation and irrigation. The large deviations between observations and computation during the first year are caused by inconsistent initial groundwater heads, which are “forgotten” after about 1 year. Remember that the steady-state solution taken as the initial head distribution corresponds more to the average situation over the 9-year period than to the initial distribution in 2003. The E_RMS and E_AM between observations and simulated groundwater heads are 1.62 and 1.18 m, respectively. The relatively small values of E_RMS and E_AM imply that there is no big variance in the individual errors between observations and calculations. These errors are still acceptable when compared to the groundwater-level variations in Guantao. Although the fit is very good due to the adjustment of pumping rates, it has to be kept in mind that there are many ways to do this adjustment, and final results are uncertain. It is clear that the calibrated model used here is not unique due to the resulting correlation among parameters and pumping rates. But this calibrated model is a completely consistent model considering the iterative calibration process. Due to the constraints on parameters, the coefficients of the covariance matrix stay relatively small (the largest matrix element in Table 2 is 0.58 in modulus). Therefore, the calibrated model can be regarded as a representative realization of reality and used for further analysis. All aspects of uncertainty will be treated in a real-time model version in which uncertainty is taken into account by ensemble techniques (paper in preparation).

Figure 8 Groundwater contour map of the steady state (in m a.s.l.).

3.2.2 Large flow model

Sources and sinks in the large flow model are conceptualized similarly to the Guantao model. The infiltration from precipitation and infiltration from applied irrigation water were calculated with Hydrus 1-D and combined in the groundwater recharge package. The Weixi channel infiltration is simulated as a defined flux boundary instead of using the river package. The infiltration flux is extracted from the Guantao model and then assigned to corresponding cells in the large model. The groundwater discharge from pumping wells is simulated in the well package. An adequate zero-flux boundary of the large flow model is defined by testing the sensitivity of the heads at the boundary location to sources and sinks in Guantao. The boundary is sufficiently far away from Guantao if heads induced by Guantao sources and sinks are close to zero. As a result, the following boundary was chosen: in the west the model is bounded by the mountains, while it is bounded by the Yellow River in the east. These boundaries coincide with the natural boundaries of the NCP. The northeast and southwest are chosen approximately parallel to groundwater contour lines (Fig. 6a), which are sufficiently far away to exclude any influence of Guantao sources and sinks. As only the effect of Guantao drivers is of interest in the large model, all sources and sinks outside of Guantao are set equal to zero.

To solve the large model, its parameter distributions are required. The absolute values of parameters outside of Guantao in the large model are not decisive in the study; only the parameter values close to the Guantao boundary are important. The hydraulic conductivity and specific yield within Guantao have the same values as the calibrated parameters from the Guantao flow model, while the model parameters outside of Guantao are taken from the literature (Liu et al., 2011). The parameter distributions used in the large model are shown in Fig. 7. The sensitivity of the results with respect to the large model's parameters will be explored later. The large groundwater flow model is also constructed using VMF-USG. The modeled area is discretized into 144 829 cells horizontally (Fig. 6). The transmissivity is calculated by multiplying the hydraulic conductivities by the aquifer thickness obtained from the Guantao model in each time step, which means that the transmissivity is constant in each time step but is updated from time step to time step to take into account the changing saturated aquifer thickness. The aquifer thickness outside of Guantao is assumed equal to the average aquifer thickness of the Guantao area. The initial condition is set to zero, which is based on the assumption that the inside drivers of Guantao are in a quasi-steady state. This is consistent with the Guantao flow model, which uses the results from the steady-state model as the initial condition in the transient model.

3.2.3 h_out flow model

The model for outside contributions is a small-scale model with the Guantao administrative boundary as the model boundary. The model structure and parameters are defined in the same way as in the Guantao flow model. Since this model is used to determine the groundwater head within Guantao caused by outside drivers, the source and sink terms inside are set to zero. The transmissivity used in the model is defined in the same way as in the large model. The specified head values on the Guantao border are obtained by subtracting the head changes determined by the large flow model from the measured heads on the boundary. Since the initial heads for the large model are zero everywhere, the initial heads for the h_out flow model are equal to the initial heads of the Guantao flow model.

4 Results and discussion

4.1 Groundwater flow characteristics in Guantao

The groundwater contour map calculated from the steady-state model is presented in Fig. 8. Two groundwater depression cones have formed in Guantao, which have a dominant influence on the groundwater flow direction, with head gradients from the outside towards the centers of the depression cones. The large-scale regional groundwater flow is directed across the NCP from southwest to northeast. Due to over-pumping, the groundwater level has been declining continuously since the 1960s (Cao et al., 2013). The intensive local pumping caused the development of numerous depression cones in the NCP. Analysis of past data shows that the depression cone in Guantao formed in the early 2000s (Zheng et al., 2010; Shao et al., 2013; Cao et al., 2013). The groundwater flow characteristics are consistent with the conclusions from former studies. Groundwater-level gradients are steeper in the east than in the west, which is caused by the higher observed head values on the eastern boundary, especially in its southeastern and northeastern sections. The groundwater recharge from precipitation and irrigation amounts to 141 000 m³ d⁻¹, which is over 87 % of the total input into Guantao. The net inflow from the specified head boundary accounts for only 12 % of the total input. The Weiyun River is located on the eastern boundary of Guantao County. Since the Guantao boundary is simulated as a specified head boundary, the infiltration from the Weiyun River is contained in the net inflow from the head boundary. The infiltration from the Weixi channel constitutes about 1 % of the total inflow. The main sink is due to pumping wells extracting groundwater from the aquifer. It amounts to over 96 % of the total outflow. The rest of about 4 % is due to boundary flow directed to the outside. The imbalance between inflow and outflow through the specified head boundary amounts to about 8 % of the total inflow.

Figure 9 Annual water budget for the period 2003–2012.

The groundwater balance of the transient model is presented in Fig. 9. The boundary inflow and main channel infiltration are not shown in the figure because they occupy only a small part of the total groundwater recharge. The infiltration from precipitation and irrigation fluctuates gently with time due to the delayed and damped groundwater recharge from precipitation and irrigation backflow through the long soil column as discussed in Sect. 3.2.1. The pumping rate fluctuates strongly in the first 6 years, while after that its variability decreases. This is due to the temporal behavior of precipitation, which directly (and without delay) influences the groundwater withdrawals for irrigation in summer (Fig. 2). The aquifer storage increases when groundwater recharge exceeds groundwater withdrawals. This is the case in 2003, where groundwater storage recovers by 19.6×10⁶ m³ a⁻¹. In the following years, it fluctuates negatively and positively and reaches a smaller value of recovery of 2.7×10⁶ m³ a⁻¹ at the end of the simulation period. Between 2003 and 2011 the groundwater storage is depleted at a rate of 0.18×10⁶ m³ a⁻¹ only, which confirms that using the data set of 2003–2011 for the steady-state model was justified. The dashed line in Fig. 9 represents the annually averaged groundwater-level dynamics in Guantao County. It is obvious that the groundwater head increases with decreasing groundwater withdrawals.

Figure 10 The spatial distribution of the groundwater head and its decomposition in Guantao at the final time step. (a) refers to the spatial distribution of head changes due to inside drivers calculated by the large model; (b) is a zoom into the large model showing only Guantao County (Δh_Guantao); (c) is the groundwater head map of Guantao obtained from the Guantao model (h_Guantao); (d) is the groundwater head distribution determined by outside drivers and obtained from the h_out flow model (h_out); (e) shows the inconsistency of Δh_Guantao in the large model; (h) refers to the distribution of the updated Δh_Guantao; (g) shows the error of the updated Δh_Guantao.

4.2 Groundwater-level changes due to inside and outside drivers

After running the three corresponding numerical models, the spatial distribution of the groundwater level and its decomposition at the final time step are presented in Fig. 10. The large light green area in Fig. 10a with Δh_Guantao values around zero indicates that the large model boundary is sufficiently far away from Guantao to make the chosen boundary conditions irrelevant. Δh_Guantao gradually changes to non-zero values around and in the Guantao area, where inside drivers dominate. Within Guantao the values of Δh_Guantao change from −5 m in the south to 4 m in the north (Fig. 10b), which means that the driving forces inside Guantao did not contribute to groundwater head changes (Δh_Guantao) in the central part of the county. The groundwater depression cone in the south is attributed to the locally intensive groundwater exploitation. The groundwater head contour map in Fig. 10c shows that groundwater in Guantao flows from northeast and southeast to the groundwater depression cone, which is consistent with the result obtained from the steady-state model. The groundwater head gradients are steep in the east, while they are flat in the west. The groundwater flow field determined by outside drivers is more uniform than changes superimposed by inside drivers (Fig. 10d). Its flow is from east and south to northwest. The regional groundwater flow directions in Guantao are governed by both inside and outside drivers, while the groundwater head depression cones formed in Guantao are caused by inside drivers.

The groundwater head in Guantao and its decomposition obtained from numerical models are supposed to satisfy Eq. (1). The difference between the Δh_Guantao calculated by subtracting Fig. 10d from Fig. 10c according to Eq. (1) and the Δh_Guantao obtained from the large flow model is shown in Fig. 10e. Without an error it should be zero everywhere. Apparently, the error is quite large in the over-pumped area. It varies from −1.7 m in the south to 0.8 m in the north. This error is caused by the difference in the types of boundary conditions used in the three numerical models. For the large model, there are no specified heads on the Guantao border. However, if one superimposes two solutions, the boundary has to be superimposed consistently as well. To prove this, another small-scale numerical model within Guantao's administrative boundary is used to update Δh_Guantao within Guantao County. This Δh_Guantao flow model within Guantao has the same model structure, parameters and sources and sinks as the Guantao model. The only difference is that the values for Δh_Guantao along the specified head boundary are obtained from the large model. After running the Δh_Guantao flow model within Guantao with this boundary condition, the new distribution of Δh_Guantao is shown in Fig. 10f. The distribution of the updated Δh_Guantao is the same as shown in Fig. 10b, but the absolute values are different, varying from −7 m in the south to 5 m in the north. This shows that the specified head condition on the Guantao border limits the influence of inside drivers. The difference between the Δh_Guantao calculated according to Eq. (1) and the Δh_Guantao obtained from the small-scale model is shown in Fig. 10e. The error is uniformly very small, with a maximum of 0.03 m (Fig. 10g). The error distribution in other time steps can also be calculated in the same way. The largest error over the whole time period is around 0.04 m, which is acceptable given the limited accuracy of the numerical solution. This means the correct results are obtained if we calculate Δh_Guantao according to Eq. (1). This method is chosen in all the following discussions.

The groundwater head in the Guantao flow model can be decomposed at any grid point and at any time. The time series of groundwater heads and their decomposition at three chosen locations are shown in Fig. 11 (two are located within Guantao and one is on the boundary). The black solid line refers to h_Guantao obtained from the Guantao model, the dashed line represents h_out obtained from the h_out flow model in Guantao and the dotted line is Δh_Guantao, which is calculated using Eq. (1). Panels (a) and (c) show that h_out changes more smoothly in time than h_Guantao and Δh_Guantao. The dynamics of h_Guantao follow the same pattern as Δh_Guantao, which means that the temporal groundwater head fluctuations within Guantao are caused by inside drivers. It is obvious that the h_out in the south increases quickly at the beginning of the simulation period and later levels off at some value. The fluctuations of h_out on the boundary are comparable to those within Guantao, because the specified boundary heads, containing both the influence from inside and outside drivers, do not fluctuate much (panel b).

Figure 11 Time series of groundwater heads in Guantao at three selected locations and their decomposition.

4.3 Inside and outside contributions

The contour color map of groundwater head changes between the last and first time steps and its decomposition into inside and outside contributions are shown in Fig. 12. Figure 12a shows that the groundwater heads decrease in the south and increase in the north, with δh_Guantao varying from −2 to 1.5 m. A similar trend is also observed in the distribution of δΔh_Guantao (Fig. 12b). The value of δh_Guantao varies from −6 m in the south to 4 m in the north, while the distribution of δh_out shows the reverse trend of increasing towards the south and decreasing towards the north, with values varying from 5 to −3 m (Fig. 12c). Therefore the groundwater head changes in Guantao are determined by both inside contributions (δΔh_Guantao) and outside contributions (δh_out). The inside contribution shows a more pronounced difference between south and north, which is largely compensated by the outside contribution.

Figure 12 Spatial distribution of changes in groundwater head over the modeled time span and its decomposition.

To explore the groundwater head change and its decomposition in detail, three specific locations are chosen (Fig. 11). The groundwater head change in the north of Guantao County increases by 0.9 m over the modeled time span, resulting from a decrease in δΔh_Guantao by 0.2 m and an increase in δh_out by 1.1 m (Fig. 11a). This indicates that the outside contribution (δh_out) dominates the groundwater head changes (δh_Guantao). δh_Guantao in Fig. 11b only increases by 0.02 m, with the value of δΔh_Guantao implying an increase of 0.36 m compensated by δh_out with a decrease of 0.34 m. In this case, both inside and outside contributions play an important role in determining the final groundwater head changes in the central part of Guantao County. In the south, δh_Guantao is the result of a decrease of 5.5 m from δΔh_Guantao and an increase of 4.7 m from δh_out (Fig. 11c). In brief, the groundwater head change over time in the whole of Guantao County is around 0.12 m (δh_Guantao), of which −1.94 m is caused by inside contributions (δΔh_Guantao) and +2.06 m is due to outside contributions (δh_out). This amounts to an aggregated change in groundwater head over time of 4 m, of which 48.5 % is contributed by inside drivers while 51.5 % is contributed by outside drivers. This further proves that considering the inside contribution only is not sufficient to assess the efficiency of groundwater management in Guantao. In such areas, decomposing groundwater head changes into inside and outside contributions is particularly important. In order to do justice to the work of the Guantao Water Department, the head changes due to the inside contribution should be the correct criterion in the assessment. Its trend will show how the internal gap between pumping and recharge develops distributed in space. Only this quantity can be influenced by measures taken in Guantao.

4.4 Sensitivity of head change and their decomposition to model-calibrated parameters

The sensitivity of the model to parameter changes is analyzed to test its robustness. The parameters, hydraulic conductivities and specific yields are perturbed by an increase of 50 % in the individual runs, while the recharge infiltration ratios are perturbed by an increase of 20 %, respectively. (Parameters include four hydraulic conductivity values K1–K4 and four specific yield values S_y1–S_y4 for the Guantao model, and six more hydraulic conductivity values K5–K10 as well as six more specific yield values S_y5–S_y10 for the large model.) When hydraulic conductivities and specific yields are perturbed, the groundwater heads on the boundary in the Guantao flow model are not changed, while the Δh_Guantao on the Guantao administrative boundary obtained from the large model will be influenced. Thus the h_out values on the boundary in the h_out flow model in Guantao are also changed according to Eq. (1). As stated in Sect. 3.2.1, the recharge infiltration ratios are used to calculate the groundwater recharge with Hydrus 1-D. Unlike the hydraulic conductivities and specific yields, the perturbed recharge infiltration ratios change sources within Guantao, which will have an impact on the boundary heads. To make the boundary heads for the Guantao model consistent with the perturbed inside sources, the following procedure is implemented: first Δh_Guantao is computed with the large model using the groundwater recharge from Hydrus 1-D and the perturbed recharge infiltration ratios. The new groundwater heads on the Guantao boundary are then calculated based on Eq. (1) considering the h_out obtained in Sect. 4.2. Then the Guantao groundwater model runs using the consistent values for boundary heads and sources to calculate the groundwater heads in Guantao. In the next step, the h_out flow model in Guantao is run using the consistent values for parameters to calculate the groundwater heads in Guantao determined by outside drivers. Then Δh_Guantao is updated based on Eq. (1). Finally, the groundwater head change and its decomposition according to Eqs. (6) to (9) can be determined.

The model sensitivity to parameters is expressed by the normalized composite sensitivity (S_NC) which is defined in Eq. (12).

$$\begin{array}{}\text{(12)}& S_{\mathrm{NC},i}={\displaystyle \frac{\sqrt{\sum _j^n(\left(\frac{Y_j^i-Y_j^{i,\mathrm{0}}}{P_{\mathrm{per},i}-P_i}\right)\times P_i{)}^{\mathrm{2}}}}{n}}\end{array}$$

$Y_j^i$ is any of the modeled groundwater head changes (referring to δh_Guantao, δΔh_Guantao, and δh_out in the study) in observation well j over the modeled time span using the perturbed model parameters P_per,i (m); $Y_j^{i,\mathrm{0}}$ is the modeled groundwater head change in observation well j using the unperturbed model parameters P_i; n is the total number of observations.

The normalized sensitivities of the Guantao head change (δh_Guantao), the inside contribution (δΔh_Guantao) and the outside contribution (δh_out) to parameters and recharge infiltration ratios are presented in Fig. 13. The higher values of sensitivity to recharge infiltration ratios indicate that recharge has the strongest influence on the groundwater head changes and their inside contribution. The hydraulic conductivities and specific yields (K5 to K10 and S_y5 to S_y10) in the large model exert the same influence on the inside and outside contributions since all the parameters are located outside of the Guantao area. The sensitivities of the models to the thickness shown in Fig. 13 are obtained by increasing the aquifer thickness outside of Guantao by a factor of 2. The aquifer thickness outside of Guantao influences the values of transmissivity in the large model; therefore, the thickness influences both inside and outside contributions and has no influence on the Guantao head changes. The highest sensitivity to specific yields within Guantao (S_y1 to S_y4) was observed in the inside contribution, followed by the outside contribution and the Guantao head changes. The smallest sensitivity to hydraulic conductivities within Guantao was observed for the outside contribution. K1 and K2 exert a larger influence on the Guantao head changes (δh_Guantao) than on the inside contribution (δΔh_Guantao).

Figure 13 The normalized sensitivity of head change and contributions to parameters and recharge ratios.

4.5 Inside and outside contributions affected by the pumping rate

Similar simulations are also carried out regarding the Guantao pumping rates. In a sensitivity test the pumping rate is increased by 20 %. The procedures to guarantee consistency are similar to those used for changing the recharge infiltration ratios (see Sect. 4.4). The spatial distributions of groundwater head changes and their decomposition with the perturbed pumping rate are shown in Fig. 14. The differences in groundwater head changes range from −5 m in the southwest to −4 m in the northeast (Fig. 14a). δΔh_Guantao is around −11 m in the southwest and increases to a maximum of 1 m in the northeast. In contrast, the outside contribution δh_out is negative in the northeast and increases to 6 m in the southwest. After increasing the pumping rate by 20 %, the groundwater head change over time averaged over the whole of Guantao County is −3.32 m (δh_Guantao), of which −5.37 m is caused by inside contributions (δΔh_Guantao), while +2.05 m is due to outside contributions (δh_out). Compared to the results from Sect. 4.3, the increased pumping rate increased the inside contribution from 48.5 % to 72.4 % when compared to the aggregated groundwater head change.

Two specific locations (one in the south and the other in the north) are chosen to analyze the influences from inside and outside contributions on the groundwater head changes. The δΔh_Guantao in the north decreases by 5.1 m, while δh_out only increases by 1.1 m. The two components determine that δh_Guantao decreases by 4 m over the modeled time span. This indicates that inside contributions (δΔh_Guantao) play a much more important role in determining the groundwater head changes (δh_Guantao). In the south the groundwater head changes show a decrease of 3.9 m, of which + 4.7 m is from the outside contribution and −8.6 m from the inside contribution, which means that in this region δh_Guantao is influenced by both inside and outside contributions.

Figure 14 Spatial distribution of groundwater head change over the modeled time span and its decomposition: influence of decreased pumping rate.

4.6 Inside and outside contributions affected by specified head boundary conditions

To explore the impact of the boundary conditions on the groundwater head change and its decomposition, the Guantao model is run with specified heads on the boundary increased by 2 m uniformly. The differences between the models with and without changing the values of the boundary heads are shown in Fig. 15. Positive values indicate an increasing trend in the component compared to the results without increased specified head values. Since there are no changes in sources and sinks in Guantao, the difference in δΔh_Guantao is around zero everywhere (Fig. 15b). The differences in δh_Guantao and δh_out are increased by 2 m almost everywhere except on the boundary, where the differences are zero (Fig. 15a and c). It is clear that changing specified heads on the boundary will only have a strong impact on the outside contribution. It is also important to keep in mind that the Weiyun River infiltration is implicitly considered in the specified head boundary condition, which is supposed to influence the inside contribution.

Figure 15 Differences between the model runs with and without increased specified heads on the boundary.

5 Conclusions

This study demonstrates how to decompose the groundwater head changes within a political boundary into inside and outside contributions using numerical models. The groundwater head in Guantao was calculated using the Guantao flow model. The groundwater head changes on the Guantao boundary caused by inside drivers were computed in a large model with a boundary far from the Guantao administrative boundary. The difference between values on the Guantao boundary from the two models is assigned as a specified head boundary condition for the h_out flow model in Guantao to calculate the groundwater head contribution in Guantao determined by outside drivers. There are no sources and sinks in this model. To eliminate inconsistencies caused by the different types of boundary conditions in the small and big models, the Δh_Guantao within Guantao can either be calculated according to Eq. (1) or recomputed by using a small-scale groundwater flow model with the Guantao administrative boundary as the prescribed head boundary with boundary values for Δh_Guantao obtained from the large flow model.

The steady-state model of Guantao was used to calibrate hydraulic conductivities and recharge infiltration ratios using data averaged over the 9-year time interval from 2003 to 2011. All hydraulic conductivities and recharge ratios are sensitive to head observations, with recharge infiltration ratios having higher sensitivities than hydraulic conductivities. The identified recharge infiltration ratios are correlated with each other, with correlation coefficients around −0.6. The calibrated hydraulic conductivities and recharge infiltration ratios are kept fixed in the transient model, and four zonal-specific yield values are also calibrated. Hydrus 1-D is used to generate the groundwater recharge at the groundwater table since the soil column leads to considerable delay and damping of the seepage at the surface. The calibrated groundwater heads fit observations quite well, with an E_RMS of less than 2 m.

Based on numerical models discussed in the study, the groundwater head in Guantao at any time and at any point can be decomposed into the groundwater head changes determined by inside drivers and the groundwater head determined by outside drivers. The groundwater head change over the whole modeled time span can correspondingly be decomposed. The groundwater head in Guantao decreased by −2 to 1.5 m. The inside contribution induced a decrease by −6 to 4 m, while the outside contribution caused a variation between a decrease of 3 m and an increase of 5 m. The groundwater head change and its decomposition in three specific locations show that the groundwater head change in the north is impacted mainly by the outside contribution, while in the center and south both inside and outside contributions determine the development of groundwater head changes. On average, 48.5 % of groundwater head change over time is contributed by inside drivers, while 51.5 % is contributed by outside drivers. We showed that a separation of inside and outside contributions to groundwater head changes is feasible. Only by separation of contributions can one prove whether the groundwater control efforts in Guantao are successful or not. It could well be that Guantao's management is functioning but that due to the outside contribution, the observed groundwater head changes are still negative.

The sensitivity of groundwater head change and its decomposition to various model parameters was analyzed by perturbing parameters by 50 % and 20 %. The results show that model outputs are sensitive to hydraulic conductivities, specific yields and recharge infiltration ratios within Guantao itself. Parameters outside of Guantao have an influence on the large model and the h_out flow model only. The recharge infiltration ratios have a stronger influence on the groundwater head change (δh_Guantao) and the inside contribution (δΔh_Guantao) than hydraulic conductivities and specific yields. The same is true for the pumping rates. After increasing the pumping rate by 20 % everywhere, the inside contribution is increased from 48.5 % to 72.4 % in terms of the aggregated groundwater head change. Guantao shallow aquifer is currently not or only very slightly over-pumped. It is feasible to get the aquifer into equilibrium at present groundwater levels.