Spatiotemporal Correlations between Water Footprint and Agricultural Inputs: A Case Study of Maize Production in Northeast China

: To effectively manage water resources in agricultural production, it is necessary to understand the spatiotemporal variation of the water footprint (WF) and the influences of agricultural inputs. Employing spatial autocorrelation analysis and a geographically weighted regression (GWR) model, we explored the spatial variations of the WF and their relationships with agricultural inputs from 1998 to 2012 in Northeast China. The results indicated that: (1) the spatial distribution of WFs for the 36 major maize production prefectures was heterogeneous in Northeast China; (2) a cluster of high WFs was found in southeast Liaoning Province, while a cluster of low WFs was found in central Jilin Province, and (3) spatial and temporal differentiation in the correlations between the WF of maize production and agricultural inputs existed according to the GWR model. These correlations increased over time. Our results suggested that localized strategies for reducing the WF should be formulated based on specific relationships between the WF and agricultural inputs.


Introduction
Freshwater has an uneven temporal and spatial distribution from the perspective of availability and quantity. A growing population and increased demand for agricultural products have led to increasing competition for freshwater resources over the past several decades [1][2][3]. As a sector of intense water usage, agriculture is the largest user of freshwater, accounting for roughly 70% of total water usage [4][5][6]. Water scarcity is a major constraint in some agricultural production activities. Water usage for producing crops is expected to increase by 0.7% per year from its estimated level of 6400 billion m 3 ·year −1 in 2000 to 9060 billion m 3 ·year −1 by 2050 as a result of the need to feed a projected global population of 9.2 billion [7,8]. As a result, it is important to evaluate water resource utilization in agricultural production to meet future freshwater supply challenges [9,10].
The water footprint (WF) concept provides a new way of evaluating water utilization in agricultural production [11]. For a given crop, the WF refers to the total amount of freshwater that is consumed and contaminated in the crop production process. It contains three components: the amount of rainwater consumed, the amount of irrigation water in the form of surface or ground water consumed, and the amount of freshwater required to assimilate the load of pollutants [12]. These components are defined as green WF, blue WF and grey WF, respectively [13,14]. The WF has led to a better understanding of the water requirement of crop production, which is attracting increased attention and has motivated a large number of studies [15][16][17][18][19]. Thus far, several studies have shed light on the spatial differences in the WFs of a variety of types of crops production in China. Sun [22]. These previous studies suggested that the WF exhibits variability in different regions at the national level. However, to the best of our knowledge, few studies have tried to explore the spatial pattern characteristics of the WF at the prefecture level. Studies regarding the WF prefecture levels are more important because they can directly relate water usage to crop production. Besides, the underlying factors influencing WFs were mostly investigated from the perspective of climate. In fact, the WF of crop production primarily depends on agricultural inputs rather than the regional climate and its variation [23]. Therefore, it is essential to investigate the relationships between the WF and agricultural inputs.
To better understand the correlations between the WF of crop production and agricultural inputs, an integrated approach is much needed. The geographically weighted regression (GWR) model is an expansion of the classical multivariate linear regression model, which considers not only the spatial non-stationary of relationships but also, in part, the spatial dependence, as it relies on local rather than global parameters [24]. Additionally, the GWR model has been applied in social economics, urban geography, meteorology, forestry and many other disciplines [25][26][27]. Few studies have been conducted on the analysis of correlations between the WF and agricultural inputs over space with the GWR method.
We selected Northeast China, which provides one-third of China's food production, as our study region. In this region, maize is one of the primary crops. The major objectives of this paper were: (1) to determine the spatial distribution of the WF of maize production at the prefecture level; (2) to identify the spatial patterns and clusters of the WF using global and local spatial autocorrelation analyses; and (3) to analyse the spatial relationship between the WF and the three agricultural inputs based on the GWR model. This study could provide valuable information for designing and implementing effective measures for reducing the WF of maize production in Northeast China.

Study Area
Northeast China is located between 118°53′-135°05′ E and 38°43′-53°33′ N, and consists of Liaoning, Jilin, and Heilongjiang Provinces ( Figure 1). The region covers an area of 7.92 × 10 5 km 2 , of which farmland occupies 2.64 × 10 5 ha, representing 16.5% of the total arable land in China [28]. The maize sowing area and yield at the prefecture level in 2012 are shown in Figure 2. The climate in the region is cold temperate and humid in the north and semi-humid in the west. The mean annual air temperature ranged from −4.2 °C to 10.9 °C, decreasing from south to north. The precipitation is generally concentrated in the summer and autumn (from May to October) [29], which coincides with the growing period of maize in Northeast China. The annual precipitation decreased from the southeast (1070 mm) to the northwest (450 mm). Maize, as the major grain crop in Northeast China, is highly seasonal with a single crop each year. The differences in the physical geography significantly affect the growth periods and crop coefficients of maize in each province [30] (Table 1).

Data Sources
The climate data of 36 weather stations, one location per prefecture, were taken from the China Meteorological Data Sharing Service System [31] for the period of 1998 to 2012, including monthly average maximum and minimum temperature, relative humidity, wind speed, sunshine hours and precipitation. The agricultural data, including maize yield, maize sowing area, agricultural machinery power, effective irrigation area and consumption of chemical fertilizers, were obtained from the statistical yearbooks 1998-2012 of the three provinces.

Methods
First, the WFs of maize production were accounted for in Northeast China. Next, the spatial patterns of the WF were investigated. Finally, the relationships between the WF and agricultural inputs were explored. The overall calculation process flow chart is diagrammed in Figure 3.  Figure 3. Diagram for WF assessment and the correlations between WF and agricultural inputs.

WF Calculation Methodology
The WF of crop production was estimated following the calculation framework developed by Hoekstra et al. [14], which is the total water consumption during crop growth process, including blue WF, green WF and grey WF [32], namely: total green blue grey where WFtotal is the WF of crop production(m 3 ·ton −1 ); WFgreen is the green WF (m 3 ·ton −1 ); WFblue is the blue WF (m 3 ·ton −1 ); and WFgrey is the grey WF (m 3 ·ton −1 ). The WFgreen and WFblue are calculated as crop water use divided by crop yield. The green and blue water use is calculated by accumulation of daily evapotranspiration over the growing period. green green green where CWUgreen and CWUblue are green and blue water usage (mm); ETgreen and ETblue are green and blue water evapotranspiration (mm); Y is the crop yield (ton·ha −1 ); and the factor 10 converts water depths in millimetres into water volumes per land surface in m 3 ·ha −1 .
The WFgrey is calculated as the chemical application rate to the field per hectare (AR, kg·ha −1 ) times the leaching-run-off fraction (α) divided by the maximum acceptable concentration (Cmax, kg·m −3 ) minus the natural concentration for the pollutant considered (Cnat, kg·m −3 ) and divided by the crop yield (Y, ton·ha −1 ).
Considering the availability of local data, in this study, we quantified the grey WF related to nitrogen use only. We have further taken α as 10% following Chapagain et al. [12], Cnat as zero, and Cmax as 12 mg·L −1 following China water quality standards.

ET Calculations Based on the CROPWAT Model
The green and blue water evapotranspiration are estimated based on crop water requirement by employing the CROPWAT 8.0 model [33,34], which are calculated as: where ETc is the crop evapotranspiration and Peff is the effective precipitation over the crop growing period (mm). ETc is calculated by employing the CROPWAT model based on the FAO Penman-Monteith method [35].
where Kc is the crop coefficient; ET0 is the reference evapotranspiration, calculated as follows: where Rn is the net radiation (MJ·m −2 ·d −1 ); G the soil heat flux (MJ·m −2 ·d −1 ); T the average air temperature (°C); U2 the wind speed at 2 m height (m·s −1 ); es the saturation vapor pressure (kPa); ea the actual vapor pressure (kPa); Δ the slope of saturation vapor pressure versus air temperature curve (kPa·°C −1 ); the hygrometer constant (kPa·°C −1 ).
where P is the precipitation at a month step (mm).

Spatial Autocorrelation
Both global and local spatial autocorrelation test results are used to confirm the presence of strong spatial dependence and heterogeneity in the distribution of the WF, supporting the use of the GWR technique.
Global Moran's I is applied to detect the overall degree of spatial autocorrelation among the regions. The equation is given as: where n is the number of spatial unit; Xi and Xj are the values of WF of spatial units i and j, respectively; X is the mean of X; Wij is a matrix of spatial weight applied to the comparison between spatial units i and j. Generally, W is defined using an adjacency standard. When the units i are adjacent to units j, Wij is equal to 1, otherwise Wij is equal to 0. The value of Moran's I varies between −1 and 1. The significance of the spatial autocorrelation can be tested by Z statistics, which is given as: where E(I) and Var(I) denote the expectation value and the variance of Moran's I, respectively [37]. If Z-value is greater than 0 and statistically significant, an obvious positive autocorrelation exists in the spatial distribution [38].
The local indicator of spatial autocorrelation index (LISA) takes each spatial unit as the inspection object and analyses its properties of spatial autocorrelation [39]. The equation is given as: where Ii is the local autocorrelation index of WFand Zi and Zj are the standardization of values of WF of spatial units i and j.

Geographically Weighted Regression (GWR) Model
The GWR model extends the traditional regression framework by embedding the spatial location of the data into the regression parameter [40,41]. The local estimation of the parameters with GWR is expressed as follows [42]: where yi is the dependent variable in the ith location; xik is the ith value of the kth independent variable; k represents the number of independent variables; i represents the number of samplings; (ui, vi) denotes the coordinates; β0 (ui, vi) represents the intercept value, and βk (ui, vi) is the continued function; εi represents the model residual.
To calculate the spatial distribution of the weightings, the optimal bandwidth is required. The Cross-Validation (CV) or Akaike Information Criterion (AIC) is usually used for calculation. In this paper the AIC was chosen, which was fixed by the maximum likelihood principle to determine the optimal bandwidth [43].
GWR model analysis was performed using the software ArcGIS10.1. In this study, the WF of maize production was treated as the dependent variable, and three agricultural inputs, namely, agricultural machinery power (AMP), effective irrigation area (EIA) and chemical fertilizers consumption (CFC) were chosen as the independent variables of the models.

Spatial Distribution of WF of Maize Production
From 1998 to 2012, the average annual WF of maize production was 1029 m 3 ·ton −1 in Northeast China. Green WF was the biggest part of the total WF of maize production with the proportion at 51%, followed by grey and blue WF, which accounted for 28% and 21%, respectively. This suggests that the WF of maize production in this region largely relies on green water resources, i.e., consumptive use of rainwater, which had a non-uniform seasonal distribution.
The WF of maize production exhibited significant differences on the prefecture level. As can be seen from 1998, the highest values of WF were observed in Baicheng, where it was 1608 m 3 ·ton −1 (Figure 4). This could be explained by the highest evapotranspiration. Dandong took second place with 1347 m 3 ·ton −1 . The lowest WFs were calculated for Siping with 691 m 3 ·ton −1 , and Tieling with 756 m 3 ·ton −1 , which could be explained by the higher maize yields in these two prefectures. Compared to 1998, the WFs of all prefectures showed a decreasing trend in 2012. Dandong had replaced Baicheng with the highest WF, which was 1313 m 3 ·ton −1 , followed by Baicheng with 1291 m 3 ·ton −1 . Tieling had replaced Siping as the prefecture with the lowest WF, which was 667 m 3 ·ton −1 . Siping ranked second, with a value of 672 m 3 ·ton −1 .

Spatial Pattern Characteristics of the WF of Maize Production
Both the global and local spatial autocorrelation test results confirm the presence of strong spatial dependence and heterogeneity in the distribution of the WF, supporting the use of the GWR technique.
Global Moran's I of WFs in 1998 and 2012 (Table 2) were 0.5512 and 0.5530, respectively, and all Z-values were significant (p-value = 0.01, with 999 permutations), which demonstrated a positive spatial correlation of prefecture-level WFs in the entire study region. In other words, prefecture units with higher and lower levels of WFs were spatially clustered. Meanwhile, the spatial autocorrelation of WFs in 2012 was slightly stronger than that in 1998. As for the local spatial autocorrelation, LISA statistics can be classified into four categories and mapped in a LISA cluster map ( Figure 5). For instance, prefectures that have high WFs surrounded by high WFs (HH) and those with low WFs surrounded by low WFs (LL) represent potential spatial clusters. Prefectures with high WFs surrounded by low WFs (HL) and low WFs surrounded by high WFs (LH) suggest potential spatial outliers. The cluster of high WF values (HH) was located in southeast Liaoning in 1998 (composed of the three prefectures of Dandong, Anshan and Yingkou), but it increased in 2012 (composed of the five prefectures of Dandong, Liaoyang, Anshan, Yingkou and Dalian), indicating that the WF disparity in this region decreased. The cluster of low WF values (LL) included the four prefectures of Changchun, Siping, Liaoyuan and Tieling, which were primarily located in central Jilin Province, and the number of prefectures in this cluster remained relatively stable. Other forms of local spatial autocorrelation included one HL-type prefecture (Yanbian) in both years.  Table 3 summarizes the results of the GWR between the WF and agricultural inputs. Results of the GWR models demonstrated that the R 2 and adjusted R 2 values in 1998 and 2012 were approximately equal, both representing over 63% of the total variance of the WF. The local parameters are described by six-number summaries consisting of the minimum, maximum, mean, lower quartile, median and upper quartile values, which presented large differences. In 1998, the correlation between the EIA and the WF was still positive; however, the AMP and CFC were either positive or negative. In 2012, the correlations between the AMP, EIA and the WF were positive and the CFC was negative. The relationships between the WF and agricultural inputs were spatially non-stationary (Figures 6 and 7). In 1998, most regression coefficients of the AMP were found to be negative, except for south-central Liaoning Province. The CFC had a positive correlation with the WF in most areas, other than south Liaoning Province. The EIA had a positive correlation with the WF in the entire region. Compared to 1998, the regression coefficients of all agricultural inputs increased in 2012, suggesting the correlations were becoming stronger. Overall the correlations of the AMP and CFC were changing from positive to negative. The AMP and CFC had negative correlation with the WF in the entire region, while the EIA still had a positive correlation with the WF. The absolute value of regression coefficient of the AMP, EIA and CFC all exhibited a decreasing trend from the centre to the south and north.

Discussion
In this study, we mapped the distribution of the WF of maize production at the prefecture level in Northeast China. A combination of the climatic condition, growing period and fertilizer consumption lead to the diversity of the WFs among 36 prefectures. The distribution of the WF in Northeast China was clustered, indicated by the significantly global Moran's I values. The LISA method detected the regions with high or low concentrations of WF. The spatial clustered regions with high WF of maize production were distributed in southeast Liaoning Province. The soils in these regions are poor and essentially unsuitable for crop production. The spatial clustered regions with low WF of maize production were distributed in central Jilin Province, which is located on the Songnen Plain and considered the Golden-Maize-Belt in China. It can be attributed to such factors as favourable climate condition, soil quality, proper irrigation facilities and relatively higher maize yield. These findings can help policy-makers prioritize maize production area and optimize maize production distribution.
The factors impacting the WF are miscellaneous. To explore the WF and its influence mechanisms, scholars began to study the relationships between WF and individual factors, such as climatic factors and agricultural input [23,44,45]. The methodology employed in this study only included correlation analysis. However, in reality, the relationships between the WF and influential factors are not always consistent in different areas. In particular, we employed the GWR technique to examine the spatial and temporal differences in the correlation between three agricultural inputs and the WF of maize production, including the AMP, EIA and CFC. The regression coefficients of the GWR also displayed a high variability in space, which could better explain the spatial differences of the correlation between the agricultural inputs and the WF. The spatial variations in the relationships are affected by the level of agricultural development.
In addition, the GWR analysis can help decision-makers clearly see the results of past relationships between the WF and agricultural inputs on a local scale and make appropriate decisions for future planning to reduce the WF of maize production. According to the results of the GWR model, the correlation between the AMP and WF was mainly positive in the study area. Attention should be paid to improving the mechanical technologies. The AMP in south-central Liaoning shifted from a positive to a negative correlation with the WF, so the amount of machinery should be adjusted appropriately. The correlation between the EIA and WF was positive in the study area. Thus, irrigation water should be reduced in this region. Water-saving irrigation techniques, such as drip or sprinkler irrigation, should be employed appropriately to inhibit the growing WF. The CFC had a negative correlation with the WF, which indicated that properly managed fertilizer use contributed to reducing the WF. Meanwhile, the CFC shifted from a positive to a negative correlation with the WF in south Liaoning. In these areas, increasing fertilizer usage or adjusting the fertilizer proportions and increasing the usage of organic fertilizers should be considered.
To the best of our knowledge, this is the first study focusing on the spatial pattern of the WF and the relationships between the WF and agricultural inputs. Nevertheless, there are still limitations in this study. There are some uncertainties in the quantification of the WF of maize production. On the one hand, the agricultural data were primarily derived from statistical data, and there are discrepancies between these values and reality. On the other hand, the grey WF was estimated based on a simplified approach, which left out local factors that influence the precise leaching rates. In addition, our spatial analysis was based on the prefecture scale, but the distribution of the WF was also not homogeneous at the prefecture level. Therefore, this study cannot identify the spatial heterogeneity within each unit.

Conclusions
The spatial distribution and cluster areas of the WFs of maize production in Northeast China were investigated and the relationships between the WF and agricultural inputs were explored. The results demonstrated that the spatial differences of the WFs of maize production were obvious in Northeast China. Two cluster regions were detected. The prefectures with higher WFs were formed in southeast Liaoning Province while the prefectures with lower WFs were formed in central Jilin Province. From the perspective of the GWR results, the WF and agricultural inputs relationships were significantly spatially non-stationary. Finally, the combination of the GWR along with the GIS approaches constitutes a powerful tool for better displaying spatiotemporally varying relationships between the WF and agricultural inputs. Extension of this study could focus on adding other driving forces, such as population, plantation structure and policy in the spatiotemporal analysis of the WF in the studied areas.