Differential Law and Influencing Factors of Groundwater Depth in the Key Agricultural and Pastoral Zones Driven by the Minimum Hydrological Response Unit

Abstract

The water cycle in the key agricultural and pastoral zones (KAPZs) is an important factor for maintaining the stability of the ecosystem. Groundwater collection and lateral seepage are indispensable parts of the water cycle, and it is difficult to monitor the groundwater situation in each area. The strength of the alternate circulation of groundwater is directly related to the utilization value and development prospects of groundwater; therefore, creating an effective method for the detection of groundwater burial depth has become an issue of increasing concern. In this paper, we attempt to create a method for the detection of groundwater burial depth that combines cokriging interpolation, spatial autocorrelation, geographically weighted regression, and other methods to construct a quantitative relationship between different land cover types and groundwater depth. By calculating the band index of the land cover type, the groundwater depthinformation of the unknown area can be obtained more accurately. Through collaborative kriging interpolation, normalized difference vegetation index (NDVI), precipitation, and hydrogeological conditions were used as covariates. The groundwater burial depth of Wengniute Banner in 2005, 2009, 2013, and 2017 was the response variable, and the groundwater burial depth in the study area was calculated. The groundwater burial depth data after the cokriging interpolation was used to transform the raster data into vector data in space using the improved hydrological response unit (HRU) model to make it more suitable for the actual groundwater confluence. Subsequently, 551 minimum response units (MHRUs) were obtained by division, and the spatial autocorrelation analysis was performed accordingly. The groundwater burial depth in the study area is spatially distinct from east to west, and the groundwater level shows a trend of being high in the west and low in the east, gradually increasing due to precipitation and rivers. The average change of groundwater depth in the time series is not significant, but it does gradually show a trend of accumulation. According to the aggregation characteristics of spatial autocorrelation analysis, a geographically weighted regression model of groundwater depth and NDVI, normalized difference drought index (NDDI), and net relatedness index (NRI) was established. The NDVI representing the forest land and the Adjusted R² of the groundwater depth is 0.67. The NRI representing the cultivated land and the Adjusted R² of the groundwater depth is 0.8675. The NDDI representing the bare land and the Adjusted R² of the groundwater depth is 0.7875. It shows that the band index representing the ground type has a good fitting effect with the groundwater burial depth.

1. Introduction

The Earth critical zone (ECZ), proposed by the National Research Council of the United States in 2001 [1], is a complex combination of geology, geochemistry, biology, hydrology, landforms, and atmospheric processes. Research on the ECZ is flourishing in a globally coordinated environment [2], which is currently a major topic of geoscience based on systematic principles. The natural geographical environment, human social behavior, and ecological service effects are unified by ECZ in this dynamic system, that is, the key belt of the earth integrates the structure-pattern-function of the system. The atmosphere–soil circle–hydrosphere–lithosphere is integrated and resolved by this critical region, in other words, ECZ is the spatial bond of these four layers [3,4]. Groundwater is the boundary between the soil layer and the rock layer in ECZ. It is not only a buffer zone for underground soil and rock, but also frequently interacts with energy on the ground [5]. As an inseparable part of the earth’s key water cycle, groundwater undertakes a variety of ecological service functions such as flood storage and replenishment, water conservation, soil conservation [6,7,8,9,10,11], and maintains the balance and stability of the underground ecosystem. As an undercurrent runoff that cannot be directly measured, groundwater is different from surface runoff, surface runoff is easy to analyze. The convergence and recharge processes of groundwater are quite different in multiple environmental modes [12,13,14]. Changes in groundwater often cause imbalances and damage to the stability of the entire ecosystem. At present, the excessive exploitation of groundwater in China has caused a dramatic drop in the groundwater level, which is a comprehensive reflection of hydrogeological elements, including the changes in groundwater recharge, runoff, excretion, as well as the environmental problems of groundwater and its severity. The shallow groundwater level has also dropped drastically, the original swamp wetlands have been drained, the original oasis has become a bare ground, and the landscape has been destroyed [15]. All of these issues have caused great damage to the buffer zone of the ECZ, which seriously affects the stability and self-healing ability of the ECZ [16,17]. The groundwater volume could quantitatively reflect its characteristics—the change of groundwater depth is the most direct manifestation of the amount of groundwater, and thus, is the most important control indicator of groundwater management [18].

Wengniute Banner is one of the severe water-deficient areas in Inner Mongolia. The water-deficient and ecologically fragile area accounts for more than 30% of the city’s total area and is a typical research area for key agricultural and pastoral zones (KAPZs) [19,20]. Due to specific natural geography and hydrogeological conditions, groundwater has always been an important source of water for the Wengniute Banner and played an important role in the development of the regional economy [19,20,21,22,23]. However, severe water shortage, coupled with long-term over-exploitation and a lack of groundwater resources has become an ecological “bottleneck” that constrains regional economic development [22]. Most researchers study the groundwater distribution mainly through two approaches. The first approach is the model mechanism for simulating the collection and convergence process of groundwater, such as the groundwater flow model [24], land surface process model [18], transient numerical model [25], global climate models [26], regional climate models [15], etc. The second approach is using a regression model to fit the groundwater data to find the law of spatial and temporal changes of groundwater [13,14]. Some researchers took a different approach and used qualitative sensing inversion methods to study groundwater characteristics qualitatively [27]; Kotowski and Satora used the Kruskal–Wallis test (K-W) for the spatial analysis of the concentrations of selected ions in groundwater. The K-W test is the non-parametric equivalent of the one-way analysis of variance. There are also many other papers presenting a similar, modified or other methodological approach [28,29,30,31]. Rahmati investigated the delineation of groundwater potential zone based on Dempster–Shafer (DS) theory of evidence and evaluate its applicability for groundwater potentiality mapping [32]. Al-Ruzouq is based on the weighted overlay analysis of the analytic hierarchy process and is supported by random forest machine learning technology to explore the potential groundwater zones in the northern United Arab Emirates [33]. Salman applies statistical and geostatistical methods to monitoring data sets to assess the pollution risk of soil and shallow groundwater [34]. However, these models and methods do not comprehensively analyze the process of collection and convergence, and the temporal and spatial variation of groundwater is affected by many factors [35]; therefore, it is urgent to find a way to effectively reflect the temporal and spatial variation of groundwater; then, the characteristics and laws of groundwater volume changes under the combined action of various driving factors could be analyzed.

In this paper, the groundwater depth is presented in the spatial range using cokriging interpolation and used as the main variables including 2005, 2009, 2013, and 2017. The covariates consist of the normalized difference vegetation index (NDVI), precipitation, and hydrogeological conditions for the corresponding year. Cokriging interpolation is based not only on the spatial geographic characteristics of groundwater depth data but also on the various factors affecting groundwater depth. The synergistic kriging interpolation results have certain limitations in the response relationship between groundwater and spatial characteristics. Subsequently, the hydrological response unit theory is used to convert the grid units in the space into vector units, so that the data are closer to reality and conforms to the state of the recharge and confluence of groundwater depth. The research goal of this paper is to explore the characteristic information of groundwater burial depth in the KAPZs and analyze its distribution law on the scale of space and evolution characteristics on the scale of time. In combination with remote sensing band information, the relationship between different types of groundwater and the buried depth of the groundwater is clarified. The minimum response units (MHRUs) that provide the groundwater depth interpolation results are spatially analyzed by the spatial autocorrelation method to study the differentiation law and characteristics of Wengniute banner on time and space. The goal of this method is to provide a realistic reference value for the rational development and utilization of groundwater resources, as well as local land planning and utilization and agricultural development planning. The MHRU is an improved hydrological response unit (HRU) with the same hydrological characteristics and is divided according to factors such as soil, slope, and vegetation in the basin. The research method based on the MHRU could play an important role in the field of rational allocation of water resources in the Earth’s key zones and water environment risk assessment. The concept of MHRU was applied to the simulation of groundwater depth and scattered point data were simulated for the entire study area. Spatial autocorrelation studies the spatial distribution characteristics of groundwater depth and explores the spatial correlation model, which can better analyze the spatial autocorrelation characteristics of groundwater depth in soil from a regional perspective.

The innovations of this study are: (1) The cokriging interpolation method was used to simulate the spatial distribution of groundwater depth. (2) The accumulation characteristics of groundwater are analyzed by obtaining the spatial distribution of groundwater depth. (3) The groundwater area with aggregation characteristics was analyzed by remote sensing the waveband of land type to determine the quantitative relationship between groundwater depth and waveband combination. Based on the results of this study, by obtaining groundwater information, the groundwater depth in the space can be calculated through collaborative kriging interpolation. Without groundwater information, it can be quantitatively retrieved through remote sensing bands, and groundwater depth information with a certain accuracy can also be obtained. This study has certain limitations. The accuracy of the model is higher in farming and pastoral areas with low vegetation coverage and it is also suitable for arid and semi-arid areas. The error is larger in high altitude alpine mountains and dense forests with high vegetation coverage.

2. Materials and Methods

2.1. Study Area

2.1.1. Basic Situation

Wengniute Banner is located in the central part of Chifeng City, Inner Mongolia Autonomous Region. It is connected to Aohan Banner and Songshan District of Chifeng City in the south, Keshiketeng Banner in the west, Linxi County, Balin Right Banner, and Aru Kerqin Banner in the north, and Neiman Banner in Tongliao City in the east. Wengniute Banner is located between 117°49′–120°43′ east longitude and 42°26′–43°25′ north latitude [21,22]. The west is an alpine platform, the middle is a hilly area with gentle undulations and mountains and rivers, and the east is the Horqin Sandy Land. Within the study area, there are both low-lying marshes between the sand and the meadows and plains formed by the impact of the water flowing from the intersection of the Laoha River and the Xilamu River. It forms the complex terrain known as the “Five Sands and Four Hills and One Field” [20]. The total area studied is 11,889 km². The area belongs to the temperate semi-arid continental monsoon climate with 2850–3000 h of annual sunshine duration, 5.5–6.4 °C of annual average temperature, and 300–400 mm of the average annual precipitation. The precipitation in the western mountainous area is 400 mm, the central hilly area is 360 mm, and the eastern bare ground area is 300 mm. The spatial and temporal distribution of groundwater in the area is extremely uneven, depending on atmospheric precipitation recharge, farmland irrigation, dry channel infiltration, and reservoir seepage [20,23,36,37].

2.1.2. Hydrogeological Conditions

The study area belongs to the first-level tectonic unit of the North China plate and is bounded by the Chifeng–Kaiyuan deep fault. It is divided into two second-level tectonic units, the continental margin of the northern part of North China and the North China block. The platform area is small, and the folds are not obvious. The scale, occurrence, nature, and types of fault structures are complex and changeable. According to their geological types, they are divided into fault basin boundary faults, general faults, and volcanic structure faults. The basin boundary fracture is the main boundary that divides the groundwater system or subsystem. Tensile faults and partial compressive faults in the bedrock area can store water, conduct water, or store or conduct water on one side, and loose layers or shallow overburden areas and rock formations with fractured pore water are the main water-bearing horizons.

Based on the occurrence conditions and hydraulic characteristics of groundwater, the main types of aquifers (groups) in this area include loose rock pore water, bedrock fissure water, and fracture zone vein water. The spatial distribution of water richness and water quality of each aquifer is quite different.

Loose rock pore water aquifers are mainly distributed in river valley plains and loess hilly areas. Depending on the lithology, structure, and thickness of the aquifer, their water richness varies in different regions. The first and second terraces on both sides of the Jiaolai River, Yangxu River, Mengke River, Shaolang River, and Laoha River Valley are rich in water. The water in bedrock fissures includes layered rock fissure water and massive rock fissure water. The fissure water of layered rocks is distributed in the low and middle mountains and low hilly areas in the east and middle. The vein-like water aquifers in the fault zone are mainly distributed in the low mountain and hilly area in the southern part of the area. The fault structure is relatively developed, with east–west, north–east, and north–west faults. In the fault structure, several groups of fracture zone rocks are relatively broken. The fissures are developed, the opening is good, and the filling is less.

2.2. Data Sources

There are 40 monitoring sites in Wengniute Banner monitoring well (see the Appendix A

Table A1. Groundwater depth monitoring well data.

Number	2005	2009	2013	2017	Average	Max	Min	Elevation	Longitude	Latitude
1	2.078	2.887	4.589	2.998	3.138	4.589	2.078	642	118.9792	42.91196
2	2.865	2.999	4.865	2.445	3.2935	4.865	2.445	635	118.9726	42.9413
3	3.125	3.485	5.004	2.754	3.592	5.004	2.754	636	118.9824	42.96041
4	2.587	2.658	4.689	2.885	3.20475	4.689	2.587	635	119.0721	42.96381
5	3.004	2.778	5.114	3.005	3.47525	5.114	2.778	620	119.0986	42.933
6	2.745	3.658	5.554	3.251	3.802	5.554	2.745	641	119.0953	42.91452
7	3.124	2.635	3.554	2.041	2.8385	3.554	2.041	663	118.9579	42.86041
8	3.21	2.953	3.675	2.554	3.098	3.675	2.554	684	118.9959	42.74533
9	3.458	3.758	4.441	3.335	3.748	4.441	3.335	668	119.0193	42.72186
10	2.425	2.554	4.587	2.561	3.03175	4.587	2.425	760	118.7696	42.83954
11	3.425	2.865	4.521	3.741	3.638	4.521	2.865	618	119.281	42.88964
12	3.147	3.963	4.895	3.612	3.90425	4.895	3.147	617	119.2517	42.81263
13	2.945	2.547	4.521	2.714	3.18175	4.521	2.547	657	118.9185	42.64127
14	2.741	2.587	3.987	2.885	3.05	3.987	2.587	665	118.8963	42.60969
15	2.947	2.494	4.101	2.774	3.079	4.101	2.494	782	118.7143	42.59284
16	3.001	2.864	3.225	2.414	2.876	3.225	2.414	998	118.3486	42.79786
17	3.879	4.341	5.789	4.256	4.56625	5.789	3.879	1385	118.1936	42.78704
18	3.457	4.145	5.424	4.741	4.44175	5.424	3.457	1528	118.1171	42.77843
19	2.589	2.014	1.471	1.552	1.9065	2.589	1.471	639	118.9356	43.07449
20	2.886	2.226	1.487	2.654	2.31325	2.886	1.487	573	119.0221	43.12981
21	2.587	2.784	1.875	2.419	2.41625	2.784	1.875	601	119.0019	43.10315
22	3.784	2.994	2.789	3.197	3.191	3.784	2.789	631	118.6987	43.18637
23	3.147	3.797	2.741	2.778	3.11575	3.797	2.741	635	118.6464	43.17077
24	4.854	4.112	7.001	5.841	5.452	7.001	4.112	534	119.3318	43.11926
25	3.647	3.695	4.365	5.014	4.18025	5.014	3.647	472	119.4026	43.23153
26	3.124	3.521	3.145	3.058	3.212	3.521	3.058	488	119.3483	43.23171
27	2.457	2.471	2.774	3.21	2.728	3.21	2.457	428	119.6696	43.24446
28	2.741	2.114	2.014	2.874	2.43575	2.874	2.014	408	119.7832	43.28839
29	3.854	3.482	3.417	4.978	3.93275	4.978	3.417	403	119.8577	43.24789
30	3.004	3.854	4.147	3.21	3.55375	4.147	3.004	385	119.9504	43.30288
31	3.147	3.664	3.564	3.684	3.51475	3.684	3.147	380	120.0072	43.32779
32	4.21	3.245	4.568	3.447	3.8675	4.568	3.245	361	120.1048	43.29964
33	3.254	3.789	4.68	3.658	3.84525	4.68	3.254	345	120.2033	43.34744
34	2.471	2.541	3.221	1.995	2.557	3.221	1.995	340	120.3062	43.36716
35	3.244	3.015	3.965	2.994	3.3045	3.965	2.994	321	120.3937	43.33634
36	3.785	4.174	4.621	3.024	3.901	4.621	3.024	318	120.4803	43.32451
37	3.417	3.784	3.689	3.254	3.536	3.784	3.254	299	120.5838	43.35076
38	2.11	2.056	2.254	1.598	2.0045	2.254	1.598	291	120.6672	43.38781
39	3.658	3.624	4.102	3.641	3.75625	4.102	3.624	310	120.5412	43.3055
40	4.995	5.841	5.669	5.604	5.52725	5.841	4.995	315	120.4593	43.21318
41	4.251	5.364	5.107	5.108	4.9575	5.364	4.251	312	120.5051	43.24379
42	4.025	4.851	4.952	5.11	4.7345	5.11	4.025	310	120.5436	43.25948
43	5.841	5.106	6.521	5.669	5.78425	6.521	5.106	352	120.2958	43.10887
44	3.669	3.554	3.964	3.051	3.5595	3.964	3.051	338	120.2719	43.00264
45	4.551	3.945	4.957	4.116	4.39225	4.957	3.945	415	119.9013	43.02574
46	5.471	4.865	6.665	6.332	5.83325	6.665	4.865	722	119.4228	42.94605
47	3.654	3.147	4.102	3.169	3.518	4.102	3.147	439	119.6713	42.81367
48	3.547	3.697	4.657	3.654	3.88875	4.657	3.547	452	119.5606	42.75841
49	4.19	4.614	4.152	4.251	4.30175	4.614	4.152	494	119.5032	42.67225
50	3.621	3.487	3.981	3.547	3.659	3.981	3.487	437	119.6149	42.69795
51	3.254	3.597	4.514	3.642	3.75175	4.514	3.254	693	119.0223	42.74746
52	3.147	3.417	3.654	3.067	3.32125	3.654	3.067	721	118.8507	42.81537

), as shown in Figure 1. The monitoring wells have been observed since 2005, and the observation data until December 2017 have relatively complete information. When monitoring groundwater data, follow the groundwater monitoring well construction specification (DZ/T 0270-2014) and the “Technical Specification for Long-term Groundwater Dynamic Observation (MT/T 633-1996)”. Regarding the accuracy of water level monitoring, when measuring the static water level, the maximum error of two consecutive measurements should not be greater than ±1 cm/10 m. Every time the water level is monitored, whether the observation well has ever been pumped is recorded, and whether it has been affected by the pumping of nearby wells. The Landsat TM images in 2005 and 2009 and Landsat OLI images in 2013 and 2017 (http://www.gscloud.cn/) were selected to perform band synthesis, image enhancement, and geometric correction using Environment for Visualizing Images (ENVI) 5.1 software. ENVI, as a complete remote sensing image processing platform, is the flagship product of Exelis Visual Information Solutions in the United States. Then, through the establishment of interpretation impact and field investigation, the land-use status and NDVI in 2005, 2009, 2013, and 2017 were obtained. The DEM (digital elevation model) was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn/). Data from the precipitation site near Wengniute Banner were obtained; data for Bahrain Zuoqi, Linxi, Wengniute Banner, Baoguotu, and Chifeng can be obtained from the National Meteorological Information Center (http://www.nmic.cn/). The soil lithology, soil particle size, topography, and other data in hydrogeological conditions were obtained from the Chengdu Institute of Mountains, Chinese Academy of Sciences.

2.3. Pretreatment

The groundwater depth data for the summer of 2005, 2009, 2013, and 2017 were selected as the basic data for spatial interpolation (

Table 1. The 2005, 2009, 2013, and 2017 groundwater information.

Years	Average	Max	Min
2005	3.391	5.841	2.078
2009	3.435	5.841	2.014
2013	4.179	7.001	1.471
2017	3.449	6.332	1.552

).

The seasonal changes in the depth of groundwater at various points are shown in Figure 2. The seasonal changes are more obvious. Due to more rain in spring and summer, the depth of groundwater is generally higher than in autumn and winter. The error of seasonal variation of groundwater depth at most points is within ±2 m, but from the data in the figure, the difference between the general minimum and maximum is too large, and the average value can accurately represent the characteristics of groundwater most of the time. Due to the limitations of the experiment in this study, it is impossible to use all the data for the year. The maximum and minimum values are not too different between months; therefore, the average value of the groundwater depth in a year is selected to represent the characteristics of the groundwater depth.

2.4. Data Processing

In this paper, groundwater depth data were used as the host variable, and NDVI, precipitation, and hydrogeological conditions were used as covariates. The groundwater depth data within the study area were obtained using cokriging interpolation. At the same time, MHRU was divided according to DEM and land cover. Using MHRU as the carrier, through spatial autocorrelation analysis of three types of land with obvious spatial aggregation characteristics: forest land, cultivated land, and bare land. Then, the corresponding remote sensing inversion band index NDVI, net relatedness index (NRI), and normalized difference drought index (NDDI) were selected for the geographically weighted regression (see Figure 3). This model is mainly suitable for bedrock and loess hilly areas. The accuracy is higher in areas where the aquifer types are loose rock pore water, bedrock fissure water, and fracture zone vein water. Fine-grained lithological combinations and arid and semi-arid areas with poor diving are more suitable.

The originality of this paper is to construct a spatial simulation method of groundwater depth and to present groundwater information in space through interpolation analysis. Compared with traditional kriging interpolation, cokriging interpolation combines hydrogeological conditions, precipitation and vegetation coverage, and other influencing factors, highlighting anisotropy in the interpolation model, making the simulation results more consistent with the actual groundwater cycle characteristics, and at the same time, also solving the unrealistic limitation of the traditional interpolation model due to single variable interpolation. At the same time, in subsequent research, a regression model between groundwater depth and remote sensing bands of different land types was constructed through geographically weighted regression. Considering a situation where researchers lack groundwater data, groundwater depth information can be quantitatively simulated; however, because the construction of the groundwater spatial interpolation model uses local hydrogeological information, this model is more suitable for arid and semi-arid farming-pastoral interlaced areas, especially in areas with finer lithology and poor diving.

2.5. Methods

2.5.1. Cokriging Model

As shown in Figure 1, the groundwater depth is distributed in the form of point data. To analyze the temporal and spatial differentiation of groundwater depth (annual average groundwater depth) quantitatively and qualitatively, it is necessary to spatially present the groundwater depth data. Compared with other interpolation methods, cokriging is a strong spatial correlation between host variables interpolated within a region and multiple attribute variables [25] and adds multiple regionalized covariate variables based on the main variable. The correlation is constructed by fitting the model, which is helpful to estimate the distribution of dependent variables [38,39,40].

Cokriging interpolation simulates the distribution of groundwater depth, which uses the vegetation characteristics, precipitation, and hydrogeological conditions to influence the static characteristics of groundwater depth as covariates. The distribution of topography and the coverage of surface vegetation, which can reflect the soil’s water-fixing capacity, have a crucial impact on the recharge range, mode, and degree of groundwater depth. The groundwater depth is generally stable in places with high vegetation coverage and has not changed much over the years, but it has affected the time variation characteristics of groundwater depth.

Three major factors affect groundwater:

• The terrain influences the convergence of groundwater, which is the most significant variable affecting the spread of groundwater.
• Hydrogeological conditions affect the depth of groundwater.
• Precipitation is the main source of groundwater runoff and hydrogeological conditions are an important carrier for groundwater diversion.

These three jointly determine the spatial variation characteristics of groundwater depth; therefore, this paper selects NDVI, precipitation, and hydrogeological conditions as covariates, and performs cooperative kriging interpolation on groundwater depth.

In ArcGIS10.4 (ESRI, Redlands, CA, USA, 2016), ordinary kriging interpolation, simple kriging interpolation, pan kriging interpolation, and extract kriging interpolation were selected for interpolation comparison. Ordinary kriging interpolation is mostly used for the unbiased estimation of a single variable. Pan kriging interpolation needs to know the overall trend of the entire interpolation space. Extract kriging interpolation is a nonlinear interpolation model. The most effective interpolation method can be chosen through parameter comparison.

The impact factor data corresponding to the covariate K of the collaborative regionalization is $Z_k$(x), (k = 1, 2, …, x), where $k_0$ is the main variable groundwater depth change rate. The mathematical expectation of the rate $Z_{k_0}\left(x_i\right)$ (i = 1, 2, 3, …, n) of change of the buried depth at each groundwater survey site exists and is a fixed value $m_{k_0}$.

$Z_{Vk0}$ is the average value of the underground burial depth change rate $k_0$.

$$Z_{Vk0}=\frac{1}{V_{k_0}}\underset{V_{k_0}}{\overset{}{{\displaystyle \int }}}{Z}_{k_0}(x)dx$$

(1)

$Z_{\alpha k}$ is the average value on the carrier $V_{\alpha k}$:

$$Z_{\alpha k}=\frac{1}{V_{\alpha k}}\underset{V_{k0}}{\overset{}{{\displaystyle \int }}}{Z}_k(x)dx$$

(2)

The estimated rate $Z_{V_{k_0}}^{*}$ of groundwater depth change over the entire region is a linear combination of the effective values of the k coordinated regional variations.

$$Z_{V_{k_0}}^{\text{*}}={\displaystyle \sum }_{k=1}^K{\displaystyle \sum }_{\alpha k}^{nk}{\mu }_{\alpha k}{Z}_{\alpha k}$$

(3)

In the formula, $Z_{V_{k_0}}^{*}$ is the estimated value of the change rate of the buried depth of the main variable; ${\mu }_{\alpha k}$ is the cokriging interpolation factor weight coefficient; $Z_{\alpha k}$ is a synergistic region variable.

In this paper, the covariates are DEM, hydrogeological conditions, and NDVI, respectively using $u_i(i=1,2,\dots n)$, $v_j(j=1,2,\dots m)$, $p_s(q=1,2,\dots r)$ to represent:

$$Z_{V_{k_0}}^{*}={\displaystyle \sum }_{i=1}^n{a}_i{u}_i+{\displaystyle \sum }_{j=1}^m{b}_j{v}_j+{\displaystyle \sum }_{q=1}^r{c}_s{p}_s$$

(4)

In the formula, $a_i$, $b_j$, $c_s$ are the precipitation, hydrogeological conditions, and the weight coefficient of the NDVI, respectively [41,42,43,44].

Interpolation analysis can present the buried depth of groundwater in space and simulate the buried depth of groundwater at a specific location; however, kriging interpolation has low accuracy when there are fewer original points and uneven distribution; therefore, the use of cokriging interpolation combines multiple factors that affect the depth of groundwater to improve the model to improve the accuracy of groundwater depth prediction.

2.5.2. Minimum Hydrological Response Unit

The result of groundwater depth is spatially presented as raster data, which is interpolated according to the cokriging interpolation; however, the division unit of groundwater depth is not similar to the regular shape of the grid in the actual situation due to the spatial difference of multiple factors and the irregularity of land use types; therefore, in this paper, we borrow the concept of MHRU to transform raster data into vector data, so that the subsequent spatial analysis would be more suitable for the actual situation in the study area.

The minimum hydrological response unit (MHRU) is based on the hydrological response unit (HRU), which is the basic unit of simulation in the SWAT (Soil and Water Assessment Tool) model [45,46]. The minimum response unit is defined by land use, soil type, and slope based on HRU. Combined with the moderate index method, the adaptive index threshold system is constructed, and the sub-basin can be better divided by hydrological analysis based on DEM. Based on the generation of the sub-basin, the hydrological response unit is determined and relevant research is carried out at this scale [47,48,49].

The basin calculation can be achieved using hydrological analysis in ArcGIS. The D8 single flow algorithm principle defines the eight directions around each grid unit of the DEM [50,51]. Assuming that the rain falls on a certain grid in the terrain, the water flow of the modified grid will flow to the grid with the lowest of the eight grid tops. If the maximum descending direction of multiple cell grids stays the same, then the adjacent cell range should be expanded until the steepest descent direction is found [52,53,54]; therefore, it is feasible to use the hydrological response unit to express the variation of groundwater depth.

The HRU is delineated using the hydrological analysis tool in ArcGIS. The watershed delineated by DEM can only represent one feature of the terrain. There are a variety of land types in the same watershed, and the groundwater content of the storage is also different; therefore, the land cover data are used to divide the watershed twice, and small patches of the boundary are combined to obtain the final MHRU. The land cover and watershed division of each MHRU, which can evaluate the groundwater depth and analyze the spatial and temporal differentiation of groundwater depth, are different and have unique characteristics.

First, we assign a value to the divided MHRU. Then, through the partition module in the analysis tool, the average groundwater depth of each hydrological response unit in the space can be statistically analyzed, and then space depth can be analyzed through the attribute connection.

2.5.3. Spatial Autocorrelation

It is necessary to conduct a qualitative and quantitative study on the spatial and temporal differentiation characteristics of groundwater by coupling the well-defined MHRU and groundwater depth interpolation results [55,56]. The groundwater depth value is presented through the MHRU, and the results are subjected to spatial autocorrelation analysis. Spatial autocorrelation refers to the potential interdependence of observational data between variables in the same distribution [57]. Tobler once pointed out that “the first law of geography: anything is related to something else, but things in the neighborhood are more relevant than things in the distance.” Statistically, correlation analysis can be used to detect whether there is a correlation between the changes of the two phenomena—it is the same attribute variable of different observation objects, which is called autocorrelation [58,59,60,61]. The parameters describing spatial autocorrelation include global autocorrelation and local autocorrelation. Local autocorrelation was used in this study to express the recharge and convergence process of groundwater.

Local spatial autocorrelation is a correlation index used to measure the extent to which a local spatial unit affects the spatial autocorrelation of the overall study [62]. Local spatial autocorrelation represents the degree of correlation of the specific hydrological response unit with the groundwater depth of the adjacent hydrological response unit and is calculated as follows [63,64]:

$$\mathrm{I}=\frac{n\left(x_i-\overline{x}\right){\sum }_{j=1}^m{W}_{ij}\left(x_j-\overline{x}\right)}{{\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}$$

(5)

where $m$ is the number of hydrological response units adjacent to the hydrological response unit $i$.

Based on the groundwater depth data for several years, the spatial autocorrelation analysis can be performed according to the area divided by the minimum hydrological response unit. In this way, the trend of groundwater depth in the past two decades is obtained and presented in the form of spatial plots, and the clustering degree and related situation of groundwater depth in the actual situation can be analyzed.

2.5.4. Geographically Weighted Regression

Based on the previous period, the spatial and temporal differentiation of the groundwater depth was qualitatively analyzed using spatial autocorrelation analysis. In the process of the confluence of groundwater depth, the land cover type is an important influencing factor. Through geographically weighted regression, the response mechanism of recharge and convergence processes of different land types and groundwater depth was quantitatively analyzed. The specific remote sensing band index was selected to represent the land type that has a great influence on the groundwater depth and is represented by the MHRU. A geographically weighted regression model of the index corresponding to different land types and groundwater depth was constructed. Combined with spatial characteristics, quantitative analysis of the influence of surface characteristics of different soils on groundwater depth [65,66,67] was then carried out.

Geographically weighted regression is based on the ordinary linear regression model, which embeds the spatial position of the data into the regression equation. The formula is as follows:

$$y_i={\beta }_0\left(u_i,v_i\right)+{\displaystyle \sum }_{k=1}^P{\beta }_k\left(u_i,v_i\right)x_{ik}+{\epsilon }_i$$

(6)

where $\left(u_i,v_i\right)$ is the sampling point $i$ coordinate, ${\beta }_k\left(u_i,v_i\right)$ is the kth regression parameter at the sampling point $i$, which is a function of the geographical position, and the weight function is obtained in the estimation process.

The Gaussian kernel function to determine the bandwidth is as follows [68,69,70]:

$${\omega }_{ij}=exp\left(-{\left(\frac{d_{ij}}{b}\right)}^2\right)$$

(7)

where b is a non-negative attenuation parameter describing the functional relationship between weight and distance, called bandwidth. The larger the bandwidth, the slower the weight is attenuated with distance, and the faster the weight is attenuated.

3. Results

3.1. Cokriging Interpolation of Groundwater Depth

Based on the geostatistical guidance tool in ArcGIS, the groundwater depth was analyzed by collaborative kriging interpolation, and then different covariance functions for function fitting accuracy analysis were selected. The comparison of 11 functional models in the ring model, spherical model, three-ball model, K-Bessel model, J-Bessel model, and stability model is shown in

Table 2. Interpolation accuracy evaluation of different semivariogram models.

	Mean Error	Root Mean Square Error	Mean Standard Error	Standardized Average Error	Nugget	Range	Sill
Stable model	0.417	0.425	0.236	1.065	0.02428	2167	7.324
Ring model	0.598	1.352	0.524	1.654	0.03412	3121	6.895
Spherical model	0.889	0.578	0.634	2.006	0.1839	5083	7.119
Three-ball model	0.647	0.694	1.032	1.417	0.09647	4898	8.886
Five-ball model	1.057	1.152	0.559	1.362	0.1299	4152	8.912
Index model	1.115	0.893	0.487	0.997	0.1154	6893	7.912
Gaussian model	0.784	0.741	0.964	1.074	0.07341	3211	7.112
Rational quadratic model	0.965	1.52	0.835	1.254	0.09985	3421	6.789
J-Bessel model	1.321	0.665	1.23	1.367	0.1231	4976	6.661
K-Bessel model	−0.963	0.783	1.008	1.634	0.09763	2341	7.731
Hole model	0.742	1.005	0.641	2.104	0.06712	2283	7.901

. In the accuracy verification of the covariance model, the stable model has the best effect of the 11 models.

To verify the effect of cokriging interpolation in different experimental well densities in the study area, 10 sets of training samples were selected for comparison and verification, of which three samples were located in the eastern region, three samples were located in the central region, and four samples were located in the western region (Figure 4). The samples in the eastern region show that the groundwater depth is shallower, and the model accuracy is higher. The samples in the central region show that the groundwater is buried deep, and the model accuracy is average. The samples in the western region show that the depth of groundwater is deeper, and the accuracy of the model is higher (

Table 3. Model partition accuracy evaluation.

	Number	2005			2009			2005			2005
		Actual	Ordinary	Cokriging	Actual	Ordinary	Cokriging	Actual	Ordinary	Cokriging	Actual	Ordinary	Cokriging
East	1	2.578	3.443	2.909	2.893	3.321	2.423	2.886	3.887	3.212	2.798	3.321	2.894
	2	3.112	3.332	3.354	3.339	3.997	3.412	3.321	3.339	3.091	2.995	3.339	3.521
	5	3.098	2.778	2.898	2.997	2.897	3.119	3.712	3.012	3.445	3.213	3.681	3.102
Central	3	3.764	3.213	3.675	3.908	3.117	3.871	3.965	3.121	3.665	3.589	3.191	3.612
	4	2.997	3.664	3.398	3.776	2.987	3.879	3.121	3.985	3.996	3.302	3.441	2.891
	7	4.007	3.553	3.678	4.112	4.321	3.888	3.873	3.117	3.211	3.995	3.337	4.203
West	6	3.512	3.379	3.288	3.612	2.678	3.381	3.779	3.123	3.339	3.478	2.895	3.546
	8	3.641	3.354	3.562	3.978	3.693	3.337	3.312	3.785	2.908	3.357	3.021	3.555
	9	3.162	2.886	3.112	4.119	3.278	3.903	2.997	2.775	3.019	3.412	3.004	3.336
	10	2.998	3.339	3.095	2.895	3.871	2.669	3.109	3.861	3.797	3.112	3.658	3.389

).

3.2. Spatial Distribution of Groundwater Depth

Compared with a single groundwater depth as an interpolation factor, the cokriging interpolation makes the result smoother in space, and has no abrupt point-like aggregation characteristics, which can better reflect the spatial distribution of groundwater (Figure 5 and Figure 6). The results are more obvious in spatial differentiation. Among them, the difference in groundwater depth in the study area is obvious, and the depth of groundwater is greatly affected by land-use type. In the four years, the groundwater depth on the east side is relatively large, mainly distributed in the farmland centered on Shengli Village in the lower reaches of the Xilamulun River and the dune wasteland near Dule Benji. The central and eastern parts of Wengniute Banner are mainly plain dunes, covered with vegetation such as wild grasses. Years of drought and lack of water have resulted in lower groundwater depths within this range. Although large areas of farmland are adjacent to the Xilamulun River and there is a certain amount of crop cover, the water demand is also huge due to intensive crop planting. Residents generally dig wells around the farmland to irrigate the land nearby. The amount of water causes the groundwater depth in large-scale farmland to be generally deeper. The farmland concentration area in the lower reaches of the Laoha River is similar. The area with better groundwater depth is mainly distributed in the west. According to the results of four years’ interpolation, the vegetation coverage in the western mountainous areas is relatively high, and there are many tall and dense shrubs. The diversity of plants with improved, developed root systems and good soil environment are all conducive to the maintenance of soil moisture. A large amount of water could be fixed in the soil by the rainfall; therefore, the groundwater depth is generally shallow in the eastern part of the study area. The severely water-deficient areas in Wengniute Banner are mainly low-mountain hills, loess hills, and loess terraces, and their geology, geomorphology, and hydrogeological conditions are extremely complex. By analyzing the geological and hydrogeological environment of groundwater in the study area, the water quantity is abundant and the water quality is good, which can meet the needs of human and animal life, industry, agriculture, and animal husbandry in the area; however, water shortages are more serious in low hills, loess hills, and loess plateau areas. Factors such as precipitation, topography, lithology, and groundwater circulation are the main reasons for the serious water shortage in this area. Stratum lithology and neotectonic movement are the decisive factors affecting the enrichment of groundwater.

Comparing the groundwater depth interpolation results in 2005, 2009, 2013, and 2017, the variation law of groundwater depth in time was analyzed with four years as the boundary. In 2005, the groundwater depth was between 2.7 and 4.1 m, and the overall distribution was relatively uniform. The groundwater depth was mainly distributed in the dune area of the northeastern part of the study area, and the groundwater depth in most areas exceeded 3 m. It indicates that groundwater diversion and convergence are relatively smooth in the past year, and there are few obstacles in the transmission process. In 2009, the groundwater depth was between 2.7 and 4.3 m, and its situation was not significantly different. The area with large groundwater depth was concentrated in the northeast and southern farmland. It indicates that the groundwater depth this year was greatly affected by location coupled with low precipitation. The groundwater depth in the river area and the farmland area is relatively close, and the farmland needs excessive extraction of groundwater for agricultural irrigation. In 2013, the groundwater depth was between 2.9 and 4.5 m, and the difference was small. The groundwater diversion and convergence process in the soil was unobstructed, but the overall groundwater depth was low, and the distribution trend was similar in 2009. The precipitation is relatively average but not sufficient, and surface runoff plays a major regulatory role. In 2017, the groundwater depth was between 2.5 and 4.4 m, which was relatively evenly distributed. Based on the results of the above four years, it can be concluded that the average value of groundwater depth will not fluctuate too much in a short period, but local areas such as northeast farmland and central and eastern dunes are still changing significantly. The groundwater depth is greatly affected by factors such as precipitation in the current year.

3.3. Spatial Autocorrelation Analysis

According to the hydrological analysis based on the DEM data in ArcGIS, 6614 HRUs are obtained, as shown in Figure 7. The traditional HRU has a large number of small plaques at the boundary, and the small plaques of less than 75 hectares are eliminated by the elimination tool, and 83 HRUs are obtained. The groundwater characteristics of HRU less than 75 hectares do not significantly differ in land type. The HRU is only based on the topographical division of the study area. The spatial division of the groundwater depth is also affected by the type of land use. Cultivated land, forest land, and construction land have different groundwater diversion and confluence and impacts. Different land types have characteristics such as soil types, vegetation coverage, soil water content, and water use; therefore, the hydrological response unit based on the topography is divided according to the land-use situation. After the division, the small shift is eliminated, and finally, 551 MHRUs are obtained. The MHRU represents different land-use types and groundwater confluence directions, which is of great significance for spatially presenting the recharge and convergence of groundwater depth. The response unit is divided according to the actual situation. The transmission of groundwater depth is not bounded by a fixed grid cell size, and the MHRU is a more realistic response unit constructed according to the actual situation. A single MHRU has consistent groundwater depth and spread with the same type of land use, approximate soil type, vegetation cover, terrain slope characteristics, and groundwater flow direction; therefore, the analysis of the spatial differentiation law of groundwater depth based on the MHRU is highly feasible, scientific, and implementable.

The distribution of the MHRUs in the study area is shown in Figure 8. Different types of MHRUs have certain distribution characteristics. Among them, 50 min response units of forest land are distributed in the west and southwest of Wengniute Banner, and there are national highways as boundaries. The west of the national highway is mainly forest land, and the elevation is generally above 1000 m. The distribution of MHRUs in forestland is relatively concentrated in one large area. In the forest, the aggregation characteristics of the MHRU are prominent, and the groundwater depth recharge and convergence characteristics are consistent. The distribution of water bodies in the study area is also concentrated. The MHRU is divided into the main distribution of water in the Xilamulun River Basin in the north and the Hongshan Reservoir in the lower reaches of the Laoha River. The MHRUs of water are divided into 26 sections with small areas. The MHRUs of cultivated land are divided into 158 sections with small areas, scattered distribution, and unspecific characteristics. One MHRU for construction land was obtained and is located in the city center of Wengniute Banner. There are 230 MHRUs for grassland, which are the most common of the six land types, and have the smallest distribution unit with the widest distribution. The only characteristic is the distribution around the bare ground area. There are 86 MHRUs in the bare ground, which are mainly distributed in the middle and east of the study area. Among them, five are larger and more concentrated. The law of groundwater diversion in such minimum response units should be most pronounced.

The results of collaborative kriging interpolation are assigned to the MHRU, and then the local autocorrelation under the spatial autocorrelation in the spatial statistical tool is used to analyze the distribution and evolution of groundwater depth in four years. The spatial autocorrelation results are shown in Figure 9.

The spatial autocorrelation characteristics presented on the MHRU by four-year groundwater interpolation can be obtained. The groundwater burial depth of the 135 MHRUs in the northeast of the study area showed high accumulation characteristics from 2005 to 2017. The characteristics of high aggregation are spatially expressed as high-value areas within the range—specifically a combination of significant aggregation and a degree of homogeneity among the values in the range. Most of the high-gathering areas are arable land, and the gathering characteristics are obvious, indicating that within a large area of farmland, groundwater connectivity is better. The standardized planting crops are more consistent in water demand and water fixing capacity, combined with relatively flat terrain, making this part of the contiguous area appear as a high-value gathering area. The dune plains in the middle and east of the study area also exhibited large-scale high-value accumulation, which showed non-large-value accumulation characteristics in individual years, indicating that the spatial distribution of groundwater depth in the plain dune areas is very irregular and not stable enough. There are some areas in the plain dune area where the groundwater depth is relatively high, and there are obstacles to groundwater interaction. The flow and recharge of groundwater are affected to a certain extent, and the aggregation characteristics are not as obvious as those of cultivated land. There are 114 MHRUs in the western part of the study area, all of which are low-value aggregation areas in four years. Among them, the MHRUs are divided by forest land. The terrain in the west is too large, in the same area, the number of MHRUs in low-value aggregation areas is less than that in the east. From the distribution of the low-value gathering area, it can be found that it is mainly distributed in the range of forest land, rather than the MHRU with a high groundwater depth such as water body. It shows that the buried depth of groundwater can show large-scale accumulation in the mountainous area in the west, and the well-developed tree roots in the forest have a strong water conservation function, which promotes the accumulation and confluence of groundwater. There is no obvious accumulation characteristic in the middle of the study area, no regular characteristic for the distribution of groundwater burial depth high- and low-value areas, and no accumulation phenomenon. Bare ground and grassland are mixed in the range and the groundwater confluence process in these two types of areas is hindered greatly.

The correlation between high-value aggregation areas and low-value aggregation areas is greater than 99%. The correlation between the second-highest aggregation area and the second-lowest aggregation area is greater than 95%. The correlation between the high-value dispersion area and the low-value dispersion area is greater than 90%. Comparing the spatial differentiation characteristics of groundwater depth in four years, it is clear that there is little difference in the change of the low-value accumulation area in the MHRU of the woodland in the western mountainous area; however, in the short term, the low-value gathering areas showed obvious aggregation characteristics, and are mainly distributed in areas where forest land is concentrated. Some sub-low-value accumulation areas and low-value dispersion areas gradually decreased or disappeared. The expansion of the forest land area caused the accumulation characteristics in the west to increase year by year. In the eastern part of the study area, the overall distribution and changes from 2005 to 2017 are similar to those in the west, except for 2013, when the changes in the east are weaker, and the characteristics of cultivated land aggregation are relatively stable. This shows that in recent years, the overall groundwater depth has displayed aggregation characteristics, and the groundwater confluence process has shown an orderly state.

3.4. Geographically Weighted Regression

The groundwater depth exhibits low-value aggregation characteristics in forest land, whereas the eastern cultivated and bare land areas exhibit high-value aggregation characteristics. Based on this, a geographically weighted regression analysis of the three land types and groundwater depth was performed, and NDVI, NRI, and NDDI were used to represent forest land, cultivated land, and bare land, respectively. NDVI, as a spectral band index, mainly reflects the characteristics of vegetation cover on the surface. For mountain forest land with high coverage, it is feasible to use NDVI as the index of the forest land. For the cultivated land in the study area, repeated tests on NRI, RVI (Ratio Vegetation Index), SAVI (Soil-Adjusted Vegetation Index), and GNDVI (Green Normalized Difference Vegetation Index) and other band indexes that can reflect the growth of crops could be implemented. NRI is the band that reflects the index of the nitrogen element of the crop. NRI can be used as the band index to reflect the cultivated land due to nitrogen fertilizer application. The central part of the study area is bare land and there are significantly large areas of desertification. NDDI effectively refers to the central bare land.

The results of the geographically weighted regression model of cultivated land in the four years of 2005, 2009, 2013, and 2017 are shown in both Figure 10 and

Table 4. Net relatedness index (NRI) geographic weighted regression fitting parameters.

Year	AICc	R2	R2 Adjusted
2005	−86.14	0.9119	0.8756
2009	−43.06	0.9185	0.8890
2013	−484.0	0.9704	0.9586
2017	−61.83	0.8199	0.7543

, which shows the standardized residuals of the geographically weighted regression model of cultivated land and NRI. Patches with absolute residuals greater than 2.5 have a poor fitting effect. The closer the absolute value comes to 0, the better the fitting effect of the geographically weighted regression model. According to the graph and table analysis, the result of the geographic weighted regression of cultivated land is an improvement when compared to other methods. In four years, the coefficient R is greater than 0.8 and the R² is greater than 0.75, which shows that the geographically weighted regression model has a strong correlation with the fitting of cultivated land NRI and groundwater depth. To further analyze the response mechanism, the standardized residual error is lower in the cultivated land in the middle of the study area, adjacent to the forest land in the west and the bare land in the middle and east, and the fitting effect of geographically weighted regression is also better than the other methods. The cultivated land in the southern and northeast corners of the study area is close to the boundary of the study area, resulting in the deviation of the fitting effect. The independent variables of geographically weighted regression of cultivated land in four years are between −4 and 10. The fitting coefficient of cultivated land and groundwater depth is mainly negative in the central and western parts of the study area. The higher the NRI, the smaller the groundwater depth. The cultivated land in the area is close to the river and the woodland in the western mountain area. The developed root system of the woodland allows for water conservation. The NRI coefficient of the cultivated land is mainly positive in the northeast of the study area. Although the cultivated land in this area is also close to the Xilamulun River, the large area of cultivated land requires more water, and the wasteland on the west side will also divert groundwater runoff. The water conservation function is not high enough to maintain a better groundwater environment in the case an increase in the range of water demand.

Figure 11 and

Table 5. Normalized difference vegetation index (NDVI) geographic weighted regression fitting parameters.

Year	AICc	R2	R2 Adjusted
2005	−48.69	0.8325	0.7889
2009	−30.98	0.6041	0.4918
2013	−99.59	0.7958	0.7387
2017	−41.39	0.7422	0.6698

show the geographically weighted regression results of NDVI of forest land and groundwater depth. During 2005, 2009, 2013, and 2017, the NDVI of the forest land and the groundwater depth R² were all above 0.6, and the adjusted R² was above 0.5. NDVI and the groundwater depth have a strong correlation. The forest land is distributed in the western part of the study area, in which the residuals of some patches in the northeast, southeast corner, and west edge of the forest land are large, whereas the residuals in most other areas are within a normal function fitting interval for four years. The independent variables of the geographic weighted regression of forest land in four years are between −5 and 0.5, and the NDVI of vegetation coverage index is negatively correlated with groundwater depth. The areas with larger NDVI coefficients are mainly distributed in the smaller patches near the borders in the north and northwest of the woodland. It shows that in these areas, the impact of the tall trees is more obvious in the mountain area on the buried depth of groundwater. Part of the cultivated land in the southern woodland of the mountainous area affects the conservation of groundwater, and part of the groundwater converges to the cultivated land to supplement the water absorbed by the crops.

Figure 12 and

Table 6. Normalized difference drought index (NDDI) geographic weighted regression fitting parameters.

Year	AICc	R2	R2 Adjusted
2005	−143.7	0.8153	0.7766
2009	−44.38	0.8318	0.7823
2013	−282.2	0.9398	0.9170
2017	−54.83	0.7267	0.6805

are the geographically weighted regression results of bare ground NDDI and groundwater depth. In 2005, 2009, 2013, and 2017, the bare ground NDDI and groundwater depth R² were all above 0.72, and the adjusted R² was above 0.68. NDDI and groundwater burial depth have a strong correlation. In accordance with the analysis above, the fitting result of the geographically weighted regression model of the bare ground area is very unstable, and the inter-annual variability is large. The standardized parameters of the middle and eastern bare ground in 2005 and 2017 are relatively large. In 2009 and 2013, the normalized residuals were relatively stable in this area. The fitting of the groundwater depth of the bare ground to NDDI is greatly affected by other factors. The independent variables of bare ground geographic weighted regression in the four years range from −2 to 6. The NDDI of the land drought index is positively correlated with the groundwater depth. The fitting coefficient of NDDI is greater than two in most bare land areas, and the impact of land drought on groundwater depth increases with the increase of NDDI.

4. Discussions

The KAPZs is an important component of the ECZ. Exploring its water cycle characteristics is the basis for optimizing the use of water resources. In this study, we used a variety of theories and methods such as spatial and statistical analysis under geography to construct a groundwater response unit for typical areas in KAPZs. The response unit could be used to study the trend and law of the collection and divergence of underground runoff in space, and analyze its distribution law in two dimensions of space and time. Furthermore, the intrinsic driving force of land cover could be determined through the construction of a geographic regression model combining spatial information.

To study the characteristics of groundwater burial depth in the KAPZs, the data of groundwater burial depth needs to be displayed in space first. Compared with ordinary kriging interpolation results, the use of cokriging interpolation methods combined with precipitation, vegetation coverage, and hydrogeological conditions interpolation results can effectively avoid errors caused by the spatial distribution of observation points (Figure 5 and Figure 6), to make the spatial distribution more regular, changes more gentle and more in line with the actual differentiation of groundwater depth. Seyed Hamid Ahmadi et al. used the Dalabu Plain as the study area and compared the accuracy of kriging interpolation and synergistic kriging interpolation in groundwater depth simulation in arid and semi-arid areas, indicating that kriging is more accurate [71]. The performance of the collaborative kriging interpolation method is superior to other geostatistical methods. Adhikary applied ordinary synergy kriging to the calculation of groundwater results in two watersheds in Australia [6]. In the age structure of groundwater, modern groundwater is the uppermost and most active part of groundwater, mainly distributed within 3 m from the surface. Inverse distance interpolation has poor interpolation results for regional studies, and ordinary collaborative kriging is determined as the best interpolation method for spatial distribution.

In this paper, the ordinary interpolation results are extremely irregular in space, the spatial and temporal distribution differences are too large and scattered, and there are no obvious rules to follow. The results of the groundwater depth of the coordinated kriging interpolation have certain interannual and spatial changes over the years, but there are rules to follow. The overall distribution characteristics show a relatively stable trend, although it fluctuates slightly. The groundwater depth is between two and five, and the overall distribution shows a trend of high in the east and low in the west; therefore, it is more effective to use the kriging interpolation to analyze the results of groundwater depth affected by various factors. Compared with previous studies on groundwater, we focused on the research of spatial analysis and vector patch division. After studying 605 catchment areas, Zhang et al. concluded that the forest ratio and soil type have a significant effect on runoff and groundwater characteristics [72]. In this paper, under different land types and soil environments, the study area was divided into many response units of different sizes according to the topographic fluctuation and the direction of groundwater confluence. Through the integration and removal of fine patches, a total of 551 HRUs including forest land, cultivated land, grassland, water land, bare land, and construction land were delineated (Figure 7). The HRU can fully reflect the characteristics and laws of land distribution in the KAPZs, and provide a data carrier for the following correlation analysis and regression analysis.

According to the MHRU, the spatial autocorrelation analysis of the groundwater depth in the agricultural–pastoral ecotone was carried out. The spatial autocorrelation can effectively reflect the spatially high and low concentration characteristics of groundwater depth. The results (Figure 8) show that the groundwater burial depth shows an obvious step-like distribution, and the east shows a high-value aggregation. The wide distribution of cultivated land and bare land leads to a significant drop in groundwater level compared to the west; therefore, in the spatial autocorrelation analysis, forest land, bare land, and cultivated land have a strong correlation with groundwater depth. Yan demonstrated the correlation between cultivated land and groundwater level through linear regression and the groundwater level in areas where the proportion of cultivated land increased [73]. Scanlon et al. used the unsaturated zone model to simulate the high correlation between the bare ground and groundwater level [74]. From the qualitative research of groundwater depth to quantitative research, the driving force of groundwater depth change is to represent NDVI, NDDI, and NRI of three types of land and establish ground-weighted regression with groundwater depth. Based on the fitting coefficient R² greater than 0.8, the forest land and groundwater burial depth showed a negative correlation trend, and the bare land and groundwater burial depth showed a positive correlation trend; therefore, the study of groundwater burial depth can predict large-scale groundwater burial depths and changes based on geographically weighted regression models through NDVI and NDDI. Cultivated land in mountainous and plain areas presents distinctly different distribution rules. Due to the proximity of forest land and rivers within the range of cultivated land in mountainous areas, there is a negative correlation with groundwater depth. The cultivated land in the eastern part of the study area has a positive correlation with groundwater burial depth due to its proximity to bare land. Using the geographically weighted regression model, the quantitative relationship between ground-type band index and groundwater depth is perfectly established, and the band index can be used to predict the convergence and recharge of groundwater depth in unknown areas. Shahabeddin also studied the relationship between groundwater depth and land cover in Iran’s Hammirza Plain based on geographically weighted regression, indicating that land types have a huge impact on groundwater depth, but the study did not deal with different land cover types separately [75]. This article highlights the geographically weighted regression analysis of the forest land, bare land, and cultivated land with the greatest impact on groundwater depth; however, this study is limited to the groundwater depth of the KAPZs. There may be certain errors in high-altitude areas, densely-vegetated areas, and hilly mountain areas. Because of the uniqueness of KAPZs, we found it to be an interesting and relevant research area.

The groundwater burial depth in the KAPZs is an indicator that is hard to detect in the underground soil layer, and plays an important role as a buffer carrier and maintains the water cycle in ECZ [45,46,76,77]. The spatial–temporal differentiation of groundwater depth in KAPZs has obvious characteristics. The types of land cover with obvious hierarchical distribution in the key agricultural and pastoral belts have a decisive influence on the groundwater depth. Based on this, the research on the driving analysis of the ground cover type on groundwater burial depth in KAPZs was studied. This study fits perfectly with the confluence and recharge process of groundwater burial depth. This process implies dynamic changes of the structure-pattern-function of multi-layer intercommunication in the principle of the system and provides a certain theoretical basis and material guarantee for the study of the circle structure of the KAPZs, the water cycle, and the effective use of groundwater resources in KAPZs.

In our research, we focused on the spatial simulation of water depth. Since traditional groundwater depth monitoring uses monitoring wells to collect data, and groundwater is of great significance to human survival and life, it is necessary to obtain groundwater depth information at specific locations in space in today’s society. The groundwater system is unmonitored; therefore, we needed to study a model that can simulate the depth of groundwater in a certain area under specific conditions. In this study, two prediction models were set up to detect the presence or absence of local groundwater depth data. Based on the groundwater depth data of some points in the area, the cokriging interpolation model can be used to simulate the groundwater depth within the spatial range of the area. If there is no groundwater depth point data, the combination of remote sensing image bands can also simulate the distribution of groundwater to a certain extent. Effectively estimating the buried depth of groundwater in a specific area is of great significance for explaining the water cycle of the earth’s sphere and human development and utilization.

5. Conclusions

In this paper, we used the cokriging interpolation to calculate the groundwater depth data in the study area. At the same time, MHRU was divided by DEM and land coverage. Following this, spatial autocorrelation analysis was carried out using MHRU as the carrier, and the corresponding remote sensing inversion index was selected for geographically weighted regression. The main variables were groundwater depth in 2005, 2009, 2013, and 2017, and data that could spatially affect groundwater depth, such as NDVI, precipitation, and hydrogeological conditions, were used as covariates. The groundwater depth data are presented in space through collaborative kriging interpolation. In the four-year interpolation results, the groundwater depth data displayed little change on average, and the groundwater depth varies between 2 and 5 m.

The groundwater burial depth data after coordinated kriging interpolation was used to convert raster data into vector data spatially using the hydrological response unit model. Through hydrological analysis, the HRUs were demarcated and re-divided according to the land use type of the study area, which was the MHRU, totalling 551. Spatial autocorrelation analysis was carried out based on the assignment of the MHRU. In the four-year data, there were 135 high-value clusters and 114 low-value clusters. The overall groundwater burial depth shows aggregation characteristics, and the groundwater confluence process shows a gradually orderly state.

The overall analysis of the spatial–temporal differentiation characteristics of groundwater depth in the study area was then carried out. The groundwater has an obvious spatial differentiation pattern between the east and west, mainly from the land cover of the mountainous forest belt in the west and the farmland and dune belt in the east, as well as the difference in topography and landform. The groundwater table shows a trend of high in the west and low in the east. The average change of groundwater depth in the time series is not large but shows a gradual trend of accumulation, indicating the obstacles that hinder the recharge of groundwater are gradually being eliminated.

According to the accumulation characteristics through spatial autocorrelation analysis, the land cover type is associated with the accumulation characteristics of groundwater depth, and the characteristics of forest land, cultivated land, and bare land are prominent; therefore, a geographically weighted regression model of groundwater depth and NDVI, NDDI, and NRI was established and the spatial location information was combined. The NDVI representing the forest land and the adjusted R² of the groundwater depth is 0.67. NRI representing the cultivated land and the adjusted R² of the groundwater depth is 0.8675. NDDI representing the bare land and the adjusted R² of the groundwater depth is 0.7875. The band index representing the land type has a good fitting effect with the groundwater burial depth, and the groundwater resources that cannot be directly monitored could be simulated and predicted by diversified land type characteristics; therefore, the development and utilization of groundwater resources in KAPZs should be combined with the above characteristics, rational planning, key development, and effective use. The effective use of groundwater can better maintain the stability of the water circulation system and provide better ecological service functions, which is of great significance for optimizing the ecological pattern of ECZ.