Although wastewater reuse in agriculture can ease water scarcity, this practice also alters the variation of groundwater recharge and groundwater levels. This study employed a geostatistical method to systematically investigate the spatio-temporal variations and storage fluctuations of groundwater in a wastewater irrigation area in a southeastern suburb of Beijing. Specifically, we generated an optimal geostatistical model for measuring groundwater levels. Furthermore, we proposed that universal kriging is a suitable method for examining groundwater spatial variations, whereas a raster-based model can provide high accuracy for studying groundwater fluctuations; the nugget effect value of groundwater levels increases with increasing exploitation intensity. Our results indicated that groundwater levels increased overall in the early stages of wastewater irrigation development, followed by local increases in some pockets in the middle stages of development, large-scale increases in the late stages and an increasing variation of magnitude over time. The results also showed that groundwater level declined less on farmlands than that in urban areas, suggesting that wastewater irrigation facilitates groundwater conservation by reducing groundwater exploitation and enhancing groundwater recharge. Our results are conducive to developing an effective groundwater management plan and for improving the accuracy of groundwater resource assessments.

INTRODUCTION

Groundwater is a crucial strategic resource that is particularly important in arid and semi-arid regions. Sustainable groundwater is characterized by good water quality and a stable supply; furthermore, it is relatively immune to seasonal influences (Maiti & Tiwari 2014). Investigation of spatio-temporal groundwater variations assists with effective groundwater management planning (Jan et al. 2013; Wang et al. 2014), helps to avoid excessive extraction and geological disasters (Izady et al. 2013), and effectively protects and utilizes groundwater resources.

An increased number of groundwater-monitoring sites can provide more detailed descriptions of the distribution characteristics of groundwater (Tabios & Salas 1985). However, due to limitations in funding and management, the number of groundwater-monitoring sites is often limited and the sites are typically unevenly distributed. Therefore, establishing a highly accurate interpolation method is an effective approach to evaluate the spatial variations of groundwater levels (Sreekanth et al. 2009; Izady et al. 2011; Vafakhah 2012). Spatial interpolation is a method that predicts missing data (or data at unavailable sites) based on known data from established sites, and the techniques can be divided into geostatistical interpolation and deterministic interpolation (Yao et al. 2014; Xiao et al. 2016). More specifically, geostatistical interpolation includes ordinary kriging (OK), simple kriging (SK), and universal kriging (UK), whereas deterministic kriging interpolation includes global polynomial interpolation (GPI) and local polynomial interpolation (with the most common model being inverse distance weighting (IDW)).

Kriging interpolation, also known as spatial interpolation, is widely regarded as the best interpolation model for representing regionally variable features (Ahmadi & Sedghamiz 2007; Mutua & Kuria 2012), but it has certain drawbacks such as generating less smooth contouring (Dalezios et al. 2002; Gundogdu & Guney 2007). GPI performs a trend surface analysis of regional variables to examine long-term and global trend effects and is suitable for areas characterized by regional variables with smooth surface changes (Mutua & Kuria 2012); however, GPI is vulnerable to over-influence from extreme values (Whitten 1970; Wang et al. 2014). IDW facilitates interpolation and yields high-accuracy results based on the principle of weighted values (Xie et al. 2011), but the model requires sufficient and evenly distributed samples; otherwise, interpolation accuracy may be considerably reduced (Caruso & Quarta 1998; Arif et al. 2011).

Wastewater reuse is an important measure for alleviating water scarcity in arid and semi-arid regions (Agriculture CWWS 2012; Wang et al. 2015). Using wastewater or reclaimed water in agricultural irrigation and river recharge projects not only reduces groundwater exploitation but also increases water infiltration. Therefore, areas using reclaimed water demonstrate apparently smaller marginal reductions compared with surrounding areas (Yin et al. 2012). In general, wastewater or reclaimed water comes from areas outside the irrigation area; therefore, water reuse changes the local infiltration-recharge strength and groundwater-exploitation intensity. However, there have been few reports about optimal interpolation models and variation patterns of groundwater levels in wastewater irrigation areas in semi-arid regions (Wu 2009).

This study displays the optimal geostatistical approach for interpolation of groundwater levels and getting better knowledge about variation patterns of groundwater levels in wastewater irrigation areas in semi-arid regions. The results can not only provide an guidance for using groundwater in the study area, but also provide an effective method for studying groundwater level variation patterns in wastewater irrigation areas globally. This study examined suburban irrigation downstream from an alluvial fan and focused on the following: (1) the most suitable interpolation model for investigating the groundwater level and the variation magnitude of groundwater in the wastewater irrigation area; (2) the temporal/spatial variations of groundwater levels and the relevant factors influencing the wastewater irrigation area over 30 years; and (3) the impacts of reclaimed water and the land-use types in the research area.

RESEARCH METHODS

General description of the investigated area

As shown in Figure 1, the research area is located in the Tongzhou District of southeastern Beijing and covers 901 km2. It has a continental monsoon climate in a warm temperate zone and a long-term average annual precipitation of 524.6 mm, with 70% of precipitation concentrated between June and September (Niu et al. 2015). The investigated area is situated on the overlapping area between the Chaobai River alluvial fan and the Yongding River alluvial fan. It has flat topography with a smooth downward tilt from the northwest to the southeast over an elevation of 1.5–40.0 m, with a gradient of 0.10%–0.96%. The quaternary burial depth gradually deepens from the northwest to the southeast, achieving a maximum depth of 600 m (Guo et al. 2014). Its aquifer is primarily composed of multiple layers of medium and fine sand that present as lamellar or lens-like structures laid out between clay layers. At a depth of 80 m, there is a continuous aquitard, the superior and inferior aspects of which are about a shallow aquifer and a deep aquifer, respectively (Wang 2014). Vertical infiltration dominates the groundwater recharge of this area.
Figure 1

Location and hydrogeological map of the study area.

Figure 1

Location and hydrogeological map of the study area.

The investigated area is a major farming region of Beijing that has agricultural, residential, urban and river acreages accounting for 58.9%, 20.2%, 15.8% and 5.1% of the total area, respectively. The research area accommodates several drainage systems, including the Chaobai River, North Canal, Liangshui River and Fenggangjian River. These drainage systems are primarily responsible for drainage from Beijing. In this region, surface-water irrigation was initially adopted, and wastewater irrigation has since been developed beginning in the 1960s, with the maximum wastewater irrigation area reaching 2.66 × 104 hm2 (Wang et al. 2015). After 1980, due to diminished surface-water runoff, the use of water storage ditches combined with irrigation was gradually implemented. Since 2004, reclaimed water irrigation has been developed using secondary effluent channeled from the Gaobeidian Wastewater Treatment Plant. In 2005, 2007 and 2008, a reclaimed water irrigation area was constructed in three stages that generated a total area of 1.89 × 104 hm2 and consumed approximately 12 × 108 m3 of reclaimed water annually (Shang et al. 2016).

This study obtained groundwater level data for 31 years from 1980 to 2010 that covered 37 groundwater level monitoring wells. Among them, 28 wells were located in the research area, including 14 in the irrigation area and 14 in the surrounding areas; the other nine monitoring wells were located outside of the research area. The groundwater level was recorded every 5 days. The land-use data were based on an aerial map with a 1 m resolution that was generated in 2005.

Research methods

  • (1) Changes in the groundwater level were calculated using the following equations: 
    formula
    1
     
    formula
    2
     
    formula
    3
    where W1988, W1990, W1998, W2000, W2008 and W2010 were the groundwater levels in 1988, 1990, 1998, 2000, 2008 and 2010, respectively, and W1990-1988, W2000-1998 and W2010-2000 were the changes in the groundwater level.
  • (2) Selection of an optimal interpolation model: different models including GPI, IDW, OK and UK in the geostatistical module of ArcGIS 10.0 were used to interpolate the groundwater level and changes in the groundwater level. The optimization process included giving formulated objective functions (Rathnayake & Tanyimboh 2014; Rathnayake 2016), choosing decision variables, inserting constrains with the constraint handling technique in the optimization approach and selecting optimal solutions(Rathnayake 2015a; Rathnayake & Tanyimboh 2015). The mean error (ME), root-mean-square error (RMSE) and correlation coefficient (R2) of the different models were calculated to select the optimal interpolation model (Guo et al. 2014; Niu et al. 2015; Rathnayake 2015b).

    The MEs and RMSEs predicted by a model can reveal the predictive accuracy of the model. The ME and RMSE can be calculated with the following equations: 
    formula
    4
     
    formula
    5
    where n represents the sample quantity, z(pi) represents the observed value of location pi, z*(pi) represents the predicted value of location pi.
    The coefficient of determination R2 between a predicted value and an observed value is also an important parameter for revealing predictive accuracy (Gundogdu & Guney 2007). The R2 was calculated using the following equation: 
    formula
    6
    where represents the average value of location pi, represents the average predicted value of location pi.
  • (3) Spatial variation characteristics: by analyzing the nugget effect and the spatial variation intensity of the groundwater, the impacts of human and structural factors on groundwater levels and variation magnitudes were assessed to analyze the spatial variation characteristics of individual water level parameters.

  • (4) Geostatistical methods for interpolation: geostatistical methods for interpolation start with the recognition that the spatial variation of any continuous attribute is often too irregular to be modelled by a simple, smooth mathematical function. The covariance or variogram of the data only depends on the distance between any two values.

RESULTS AND DISCUSSION

Normal distribution test and global trend analysis

The distribution characteristics of the groundwater level can be determined by the skewness and kurtosis of a histogram of the groundwater levels. Skewness approaching 0 and kurtosis approaching 3 indicate that the data exhibit near-normal distribution (Assaf & Saadeh 2009; Saibi et al. 2009). Our results revealed that the groundwater levels W1980, W1990, W2000 and W2010 and the variation magnitudes W1990-1980, W2000-1990 and W2010-2000 all exhibited normal distributions, as shown in Table 1. The groundwater level variation magnitudes of W1990-1980, W2000-1990 and W2010-2000 did not exhibit any dominant trends (Noshadi & Sepaskhah 2005), indicating that OK was suitable for their interpolation (Assaf & Saadeh 2009). The groundwater levels W1980, W1990, W2000 and W2010 exhibited a clear firstorder dominant trend, thus UK was suitable for predictive calculations (Gundogdu & Guney 2007).

Table 1

Statistical eigenvalues of the data

ParameterMinima/mMean/mMaxima/mStandard deviation/mVariation coefficientSkewnessKurtosisDominant trend
W1980 8.38 15.15 20.29 3.72 0.25 −0.20 2.03 1st order 
W1990 5.84 12.82 19.15 4.19 0.33 −0.11 1.95 1st order 
W2000 0.83 11.34 17.77 4.25 0.24 −0.78 3.23 1st order 
W2010 −1.37 9.95 17.51 5.29 0.53 −0.48 2.30 1st order 
W1990-1980 −5.19 −2.59 −0.003 1.29 0.50 −0.59 3.26 0th order 
W2000-1990 −5.05 −1.68 1.02 1.72 1.02 −0.46 2.55 0th order 
W2010-2000 −5.65 −1.39 4.14 2.09 1.50 0.08 3.27 0th order 
ParameterMinima/mMean/mMaxima/mStandard deviation/mVariation coefficientSkewnessKurtosisDominant trend
W1980 8.38 15.15 20.29 3.72 0.25 −0.20 2.03 1st order 
W1990 5.84 12.82 19.15 4.19 0.33 −0.11 1.95 1st order 
W2000 0.83 11.34 17.77 4.25 0.24 −0.78 3.23 1st order 
W2010 −1.37 9.95 17.51 5.29 0.53 −0.48 2.30 1st order 
W1990-1980 −5.19 −2.59 −0.003 1.29 0.50 −0.59 3.26 0th order 
W2000-1990 −5.05 −1.68 1.02 1.72 1.02 −0.46 2.55 0th order 
W2010-2000 −5.65 −1.39 4.14 2.09 1.50 0.08 3.27 0th order 

Selection of the optimal interpolation model

The suitability of an interpolation simulation can be determined by comparing the differences between the fitted values and the measured values that result from different interpolation models such as IDW, GPI and OK (Adhikary et al. 2010). As shown in Table 2, the UK model yielded the lowest MEs for the groundwater level parameters W1980, W1990, W2000 and W2010, with RMSEs close to 1. The correlation coefficient R between the predicted values and the observed values using this model was the highest (0.94) (Figure 2). Therefore, the UK model was suitable for predicting spatial changes in the groundwater level in the research area. In comparison, the MEs of the groundwater variation magnitudes W1990-1980, W2000-1990 and W2010-2000 and the RMSE values of the IDW, GPI and OK models were in the ranges 0.0–0.18 and 1.21–2.07, respectively. The correlation coefficients R of these models were relatively low (Table 2), indicating a poor fit (Figure 3).
Table 2

The errors in groundwater level derived from different interpolation models

TypeParameterInterpolation modelMERMSEEquation of best fitaR
Groundwater level W1980 IDW −0.29 1.59 y = 0.6587x +4.8775 0.93 
GPI −0.02 1.85 y = 0.7856x +3.2237 0.86 
UK 0.25 1.56 y = 0.8320x +2.4652 0.94 
W1990 IDW −0.10 1.62 y = 0.7073x +3.6474 0.94 
GPI 0.01 2.26 y = 0.7581x +3.1083 0.84 
UK −0.00 1.35 y = 0.8605x +1.7845 0.94 
W2000 IDW 0.12 2.04 y = 0.6465x +4.1306 0.89 
GPI 0.05 2.59 y = 0.6621x +3.8831 0.79 
UK 0.08 1.72 y = 0.8337x +1.9609 0.91 
W2010 IDW 0.23 2.77 y = 0.6230x +3.9895 0.87 
GPI 0.07 3.45 y = 0.6212x +3.8505 0.76 
UK 0.16 2.32 y = 0.8196x +1.9648 0.90 
Groundwater level marginal change W1990-1980 IDW 0.10 1.22 y = 0.1279x −2.1570 0.30 
GPI 0.02 1.42 y = −0.0674x −2.7375 0.19 
OK 0.03 1.31 y = 0.0076x −2.5355 0.03 
W2000-1990 IDW 0.18 1.33 y = 0.3520x −0.9038 0.62 
GPI 0.07 1.52 y = 0.2914x −1.1193 0.46 
OK 0.07 1.21 y = 0.4445x −0.8624 0.70 
W2010-2000 IDW 0.15 2.07 y = 0.0893x −1.1143 0.21 
GPI 0.03 2.10 y = 0.0500x −1.2915 0.13 
OK 0.13 2.06 y = 0.0806x −1.1485 0.21 
TypeParameterInterpolation modelMERMSEEquation of best fitaR
Groundwater level W1980 IDW −0.29 1.59 y = 0.6587x +4.8775 0.93 
GPI −0.02 1.85 y = 0.7856x +3.2237 0.86 
UK 0.25 1.56 y = 0.8320x +2.4652 0.94 
W1990 IDW −0.10 1.62 y = 0.7073x +3.6474 0.94 
GPI 0.01 2.26 y = 0.7581x +3.1083 0.84 
UK −0.00 1.35 y = 0.8605x +1.7845 0.94 
W2000 IDW 0.12 2.04 y = 0.6465x +4.1306 0.89 
GPI 0.05 2.59 y = 0.6621x +3.8831 0.79 
UK 0.08 1.72 y = 0.8337x +1.9609 0.91 
W2010 IDW 0.23 2.77 y = 0.6230x +3.9895 0.87 
GPI 0.07 3.45 y = 0.6212x +3.8505 0.76 
UK 0.16 2.32 y = 0.8196x +1.9648 0.90 
Groundwater level marginal change W1990-1980 IDW 0.10 1.22 y = 0.1279x −2.1570 0.30 
GPI 0.02 1.42 y = −0.0674x −2.7375 0.19 
OK 0.03 1.31 y = 0.0076x −2.5355 0.03 
W2000-1990 IDW 0.18 1.33 y = 0.3520x −0.9038 0.62 
GPI 0.07 1.52 y = 0.2914x −1.1193 0.46 
OK 0.07 1.21 y = 0.4445x −0.8624 0.70 
W2010-2000 IDW 0.15 2.07 y = 0.0893x −1.1143 0.21 
GPI 0.03 2.10 y = 0.0500x −1.2915 0.13 
OK 0.13 2.06 y = 0.0806x −1.1485 0.21 

ax represents the observed value and y stands for the predicted value.

Figure 2

Correlation between the predicted and calculated values for groundwater level W1980 derived from the (a) IDW, (b) GPI and (c) OK models.

Figure 2

Correlation between the predicted and calculated values for groundwater level W1980 derived from the (a) IDW, (b) GPI and (c) OK models.

Figure 3

Correlation between the predicted and calculated values of the groundwater level variation magnitude W1990-1980 derived from the (a) IDW, (b) GPI and (c) OK models.

Figure 3

Correlation between the predicted and calculated values of the groundwater level variation magnitude W1990-1980 derived from the (a) IDW, (b) GPI and (c) OK models.

Figure 3 indicates that direct employment of an interpolation model yielded poor accuracy in predicting the groundwater level variation magnitude. In Figures 4 and 5, data at the monitoring sites in the irrigation area were subjected to different interpolation models to evaluate the correlation between the predicted values and the measured values at the other nine monitoring sites. Specifically, in Figure 4, the UK model was used to directly interpolate the groundwater level variation magnitude, yielding a correlation coefficient R between the predicted and measured values of only 0.084. In Figure 5, after obtaining 2 years of UK interpolation-based grid maps of groundwater levels, the differences in the groundwater levels in each grid map included the variation magnitude between the 2 years, which produced a spatial distribution map of the groundwater level variation magnitude. The correlation coefficient R in this approach was relatively high at 0.51. Therefore, when examining the pattern of spatial fluctuation in groundwater levels, it is recommended that the subtraction approach be adopted using 2 years of groundwater level data.
Figure 4

Correlations between the predicted and measured values for groundwater level variation magnitudes derived from the UK model.

Figure 4

Correlations between the predicted and measured values for groundwater level variation magnitudes derived from the UK model.

Figure 5

Correlations between the predicted and calculated values using a grid-subtraction approach for displaying groundwater level variation magnitudes derived from the UK model.

Figure 5

Correlations between the predicted and calculated values using a grid-subtraction approach for displaying groundwater level variation magnitudes derived from the UK model.

Groundwater dynamic variation characteristics

As shown in Figure 6, the UK model was used to interpolate the spatial distribution of groundwater levels. UK can be seen as a point interpolation, and is defined as a Kriging method with a local trend, which is a continuous varying trend surface on top of which the variation to be interpolated is superimposed. The geostatistical module in Arcgis 10.0 was used to create the spatial results. In 1980, the groundwater level of the investigated area was in the range 8.4–20.2 m, with an average of 15.1 m. The northern and northwestern areas had relatively higher levels, whereas the southeastern areas had shallow burial depths and the lowest water level. In 2010, the groundwater levels were in the range −1.4–17.5 m, with an average of 9.95 m. Globally, the groundwater levels gradually decreased from north to south and from west to east. Moreover, the groundwater levels displayed a diminishing trend in the northeastern urban development area as well as in southern and southeastern agricultural areas, reflecting a significant increase in groundwater exploitation.
Figure 6

Spatial distribution map of groundwater levels derived from the UK model.

Figure 6

Spatial distribution map of groundwater levels derived from the UK model.

Analysis of the causes of groundwater dynamic variation

The nugget effect of the regional groundwater level semi-variogram can be used to determine the degree of impact from anthropogenic or structural factors on groundwater spatial variations. The nugget effect value is in the range 0–1. If the value is below 0.25, the groundwater level is primarily influenced by anthropogenic factors; if it is in the range 0.25–0.75, the level is affected by both structural and anthropogenic factors; if it is above 0.75, the level is primarily impacted by structural factors.

The nugget effect values of groundwater levels W1980, W1990 and W2000 were all 0, indicating that the spatial distributions of the groundwater levels were almost completely influenced by structural factors during these years. The nugget effect value of W2010 was 0.054, indicating that the spatial distribution of the groundwater level was primarily affected by structural elements, with some impact from anthropogenic factors (Carreira et al. 2010). The nugget effect values of the groundwater level variation magnitudes W1990-1980, W2000-1990 and W2010-2000 were 0.074, 0.135 and 0.457, respectively, following a gradually increasing trend. This result indicated that the groundwater level variation magnitude was increasingly influenced by anthropogenic factors. Between 1980 and 1999, the variation magnitude was only 2.4 m, indicating that extraction and recharge were basically balanced. However, since 1999, the research area has entered a period of relatively insufficient precipitation, receiving only 75% of the long-term average annual precipitation. The lack of precipitation in conjunction with increased exploitation of groundwater has led to an imbalance between exploitation and recharge. As a consequence, the groundwater level decreased from 12.4 m in 1999 to 9.9 m in 2010. The main causes of the increase in the nugget effect values included the spatial variation in exploitation intensity and the different land uses associated with reclaimed water irrigation (Table 3).

Table 3

Geostatistical parameters for various factors affecting the groundwater level

ParametersVariation coefficientNuggetPartial sillSillNugget effect
W1980 0.25 0.00 3.28 3.28 0.00 
W1990 0.33 0.00 5.40 5.40 0.00 
W2000 0.24 0.00 7.13 7.13 0.00 
W2010 0.53 0.70 12.34 13.04 0.054 
W1990-1980 0.5 0.13 1.60 1.73 0.074 
W2000-1990 1.02 0.52 3.35 3.88 0.135 
W2010-2000 1.5 2.08 2.47 4.55 0.457 
ParametersVariation coefficientNuggetPartial sillSillNugget effect
W1980 0.25 0.00 3.28 3.28 0.00 
W1990 0.33 0.00 5.40 5.40 0.00 
W2000 0.24 0.00 7.13 7.13 0.00 
W2010 0.53 0.70 12.34 13.04 0.054 
W1990-1980 0.5 0.13 1.60 1.73 0.074 
W2000-1990 1.02 0.52 3.35 3.88 0.135 
W2010-2000 1.5 2.08 2.47 4.55 0.457 

Urbanization and land-use variations directly influence the spatial distribution patterns of groundwater (Jat et al. 2008; Choi et al. 2012; Iwasaki et al. 2013; Moukana et al. 2013). Due to considerable differences in groundwater-recharge coefficients, different land-use types cause significant groundwater spatial variation (Salama et al. 1999; Zhang & Schilling 2006). As indicated in Table 4, between 2004 and 2010, the groundwater level in the research area decreased by 0.92 m, decreasing by 1.25 m, 0.81 m and 0.59 m in the urban, farming, and river and pond areas, respectively. These results indicated that the urban area had the greatest degree of exploitation intensity that contributed to the highest levels of water reduction. This level of exploitation intensity was followed by farming areas, whereas the river and pond areas had the lowest water level reductions with relatively higher levels of infiltration recharge. The reduction in groundwater was primarily influenced by urban and farming areas.

Table 4

Changes to groundwater levels under different land uses

No.Land typeW2010-2004
Urban −1.25 
Farmland −0.81 
Rivers and ponds −0.59 
Research area −0.92 
No.Land typeW2010-2004
Urban −1.25 
Farmland −0.81 
Rivers and ponds −0.59 
Research area −0.92 

As shown in Figure 7(a), a comparison between 2006 and 2004 revealed that the variation magnitude of the groundwater levels fluctuated in the range −6.9–1.0 m, following an overall decreasing trend. The average reduction was 1.1 m, and most areas had reductions under 2 m. Overall, the western area and some southwestern pockets displayed a significant reduction in groundwater levels, including 1.42 m in the irrigation area and 0.99 m in the surrounding area. A significant reduction in the irrigation area indicated extraordinary groundwater exploitation. Although the first stage of the irrigation project was accomplished in 2005, reclaimed water was not adequately supplied during the test run, which indicated little effect on the groundwater level fluctuations in the irrigation area. After the first and second stages of the project were completed in 2005 and 2007, the systems were operating normally, which correlated with apparent changes in the groundwater level fluctuations inside and outside of the irrigation area. As shown in Figure 7(b), a comparison of 2006–2008 revealed that the groundwater level variation magnitude was in the range −2.1–1.9 m, with an average decrease of 0.1 m. In other words, the area exhibited reduced fluctuations. The overall water level of the reclaimed water irrigation area displayed an increasing trend with a 0.35 m increase in the irrigation area and a 0.02 m decrease in the surrounding area. These results revealed that construction of the reclaimed water irrigation area clearly reduced groundwater extraction. In contrast, in the southeastern farming area, where reclaimed water irrigation was undeveloped, a continuous decline in the groundwater level occurred with a maximal variation magnitude of 2.1 m. Figure 7(c) illustrates that after the third stage of the project in 2008, the 2008–2010 comparison revealed an increasing trend in the groundwater level in the investigated area that continued to progress. Specifically, the groundwater level variation magnitude was in the range −6.7–2.1 m, with an overall increasing trend. The average increase was 0.23 m, with a 0.35 m increase in the irrigation area and a 0.50 m reduction in the surrounding area. Together, our data indicated that the construction of the reclaimed water irrigation area was highly effective in reducing groundwater exploitation and promoting infiltration recharge.
Figure 7

Inter-annual variability calculation models of groundwater spatial data over different time periods.

Figure 7

Inter-annual variability calculation models of groundwater spatial data over different time periods.

CONCLUSIONS

This study paid particular attention to the statistical features and recharge/depletion patterns of groundwater in a reclaimed water irrigation area. By conducting normal distribution tests and global trend analyses of individual groundwater parameters, appropriate geostatistical interpolation models for groundwater levels were selected; we proposed that UK interpolation should be adopted for analyzing the spatial variation of groundwater levels and that a raster-based model should be chosen for evaluating groundwater level variation magnitude with a high accuracy. The nugget effect value of the groundwater level in the investigated area increased from 0 in 1980 to 0.054 in 2010, which was correlated with the nugget effect value of groundwater level variation magnitude, which rose from 0.074 for W1990-1980 to 0.457 for W2010-2000. These results suggest that groundwater level variation transitioned from being mostly dictated by structural factors to being primarily governed by anthropogenic factors. Reduced precipitation recharge, reclaimed water irrigation and intense exploitation in urban areas were the major factors contributing to the elevated nugget effect of groundwater level variation. The construction of the reclaimed water irrigation area led to decreased groundwater exploitation and increased local infiltration recharge. As such, from 2008 to 2010, increased groundwater levels in the irrigation area were 0.69 m more than in the surrounding area, indicating that the constructed project promoted groundwater conservation. In summary, our results recommend the development of an effective plan for managing groundwater and improving the accuracy of groundwater resource assessments.

ACKNOWLEDGEMENTS

This study was supported by the Fundamental Research Funds for the Central Universities (Nos 2652016020 and 2652016022) and the National Natural Science Fund (51339007) and the National Science & Technology project (2011BAD25B00 and 2012BAD08B02).

REFERENCES

REFERENCES
Adhikary
P. P.
Chandrasekharan
H.
Chakraborty
D.
Kamble
K.
2010
Assessment of groundwater pollution in West Delhi, India using geostatistical approach
.
Environmental Monitoring & Assessment
167
(
1–4
),
599
615
.
Agriculture CWWS
2012
Coping with Water Scarcity: An Action Framework for Agriculture and Food Security
.
FAO Water Reports 38, Food and Agriculture Organization of the United Nations, Rome, Italy.
Ahmadi
S. H.
Sedghamiz
A.
2007
Geostatistical analysis of spatial and temporal variations of groundwater level
.
Environmental Monitoring & Assessment
129
(
1–3
),
277
294
.
Arif
M.
Hussain
J.
Hussain
I.
Kumar
S.
Bhati
G.
2011
Effect of GIS interpolation techniques on the accuracy of the spatial representation of groundwater monitoring data in Gaza Strip
.
Journal of Environmental Science & Technology
4
(
6
),
99
102
.
Carreira
P. M.
Marques
J. M.
Pina
A.
Gomes
A. M.
Fernandes
P. A. G.
Santos
F. M.
2010
Groundwater assessment at Santiago Island (Cabo Verde): a multidisciplinary approach to a recurring source of water supply
.
Water Resources Management
24
(
6
),
1139
1159
.
Caruso
C.
Quarta
F.
1998
Interpolation methods comparison
.
Computers & Mathematics with Applications
35
(
12
),
109
126
.
Dalezios
N. R.
Loukas
A.
Bampzelis
D.
2002
Spatial variability of reference evapotranspiration in Greece
.
Physics & Chemistry of the Earth Parts A/B/C/
27
(
23–24
),
1031
1038
.
Gundogdu
K. S.
Guney
I.
2007
Spatial analyses of groundwater levels using universal kriging
.
Journal of Earth System Science
116
(
1
),
49
55
.
Guo
G.
Ju
Y.
Zhai
H.
Xu
L.
Shen
Y.
Ji
Y.
2014
Assessment of groundwater quality of different aquifers in Tongzhou area in Beijing plain and its chemical characteristics analysis
.
Environmental Science
35
(
6
),
2114
2119
.
Izady
A.
Davary
K.
Alizadeh
A.
Ghahraman
B.
Sadeghi
M.
Moghaddamnia
A.
2011
Application of ‘panel-data’ modeling to predict groundwater levels in the Neishaboor Plain, Iran
.
Hydrogeology Journal
20
(
3
),
435
447
.
Izady
A.
Davary
K.
Alizadeh
A.
Nia
A. M.
Ziaei
A. N.
Hasheminia
S. M.
2013
Application of NN-ARX model to predict groundwater levels in the Neishaboor plain, Iran
.
Water Resources Management
27
(
14
),
4773
4794
.
Jan
C. D.
Chen
T. H.
Huang
H. M.
2013
Analysis of rainfall-induced quick groundwater-level response by using a Kernel function
.
Paddy & Water Environment
11
(
1–4
),
135
144
.
Mutua
F.
Kuria
D.
2012
A comparison of spatial rainfall estimation techniques: a case study of Nyando river basin Kenya
.
Journal of Computer Science
3
(
2
),
69
92
.
Niu
Y.
Yin
S.
Liu
H.
Wu
W.
Li
B.
2015
Use of geostatistics to determine the spatial variation of groundwater quality: a case study in Beijing's reclaimed water irrigation area
.
Polish Journal of Environmental Studies
24
(
2
),
611
618
.
Noshadi
M.
Sepaskhah
A. R.
2005
Application of geostatistics for potential evapotranspiration estimation
.
Iranian Journal of Science & Technology Transaction B Engineering
29
(
3
),
343
355
.
Rathnayake
U.
2015a
Enhanced water quality modelling for optimal control of drainage systems under SWMM constraint handling approach
.
Asian Journal of Water Environment & Pollution
12
,
81
85
.
Rathnayake
U.
Tanyimboh
T.
2014
Multi-objective optimization of combined sewer systems using SWMM 5.0
. In:
7th International Young Water Professionals Conference
, pp.
175
187
.
Rathnayake
U. S.
Tanyimboh
T. T.
2015
Evolutionary multi-objective optimal control of combined sewer overflows
.
Water Resources Management
29
(
8
),
2715
2731
.
Salama
R.
Hatton
T.
Dawes
W.
1999
Predicting land use impacts on regional scale groundwater recharge and discharge
.
Journal of Environmental Quality
28
(
2
),
446
460
.
Shang
F.
Ren
S.
Yang
P.
Li
C.
Xue
Y.
Huang
L.
2016
Modeling the risk of the salt for polluting groundwater irrigation with recycled water and ground water using HYDRUS-1 D
.
Water, Air, & Soil Pollution
227
(
6
),
article 189
.
Sreekanth
P. D.
Geethanjali
N.
Sreedevi
P. D.
Ahmed
S.
Kumar
N. R.
Jayanthi
P. D. K.
2009
Forecasting groundwater level using artificial neural networks
.
Current Science
96
(
7
),
933
939
.
Tabios
G. Q.
Salas
J. D.
1985
A comparative analysis of techniques for spatial interpolation of precipitation
.
JAWRA: Journal of the American Water Resources Association
21
(
3
),
365
380
.
Wang
S.
2014
Determination of 12 Isomers of Nonylphenol in Reclaimed water by Gas Chromatography–Mass Spectrometry
.
Master thesis
,
China University of Geosciences (Beijing)
,
Beijing
,
China
.
Wang
S.
Huang
G. H.
Lin
Q. G.
Li
Z.
Zhang
H.
Fan
Y. R.
2014
Comparison of interpolation methods for estimating spatial distribution of precipitation in Ontario, Canada
.
International Journal of Climatology
34
(
14
),
3745
3751
.
Wang
S.
Wu
W.
Liu
F.
Yin
S.
Bao
Z.
Liu
H.
2015
Spatial distribution and migration of nonylphenol in groundwater following long-term wastewater irrigation
.
Journal of Contaminant Hydrology
177–178
,
85
92
.
Whitten
E. H. T.
1970
Orthogonal polynomial trend surfaces for irregularly spaced data
.
Journal of the International Association for Mathematical Geology
2
(
2
),
141
152
.
Wu
W.
2009
Research on Groundwater Vulnerability Experiment of Reclaimed Wastewater District and Irrigation Allocation
.
China University of Geosciences (Beijing)
,
Beijing
,
China
.
Yao
L. Q.
Huo
Z. L.
Feng
S. Y.
Mao
X. M.
Kang
S. Z.
Chen
J.
Xu
J. J.
Steenhuis
T. S.
2014
Evaluation of spatial interpolation methods for groundwater level in an arid inland oasis, northwest China
.
Environmental Earth Sciences
71
(
4
),
1911
1924
.
Yin
S.
Wu
W.
Liu
H.
Zhang
X.
2012
Spatial variability of groundwater nitrate-nitrogen and cause analysis of its pollution for irrigation area with reclaimed water
.
Nongye Gongcheng Xuebao/Transactions of the Chinese Society of Agricultural Engineering
28
(
18
),
200
207
.