## Materials and Methods
### 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  
Figure 3: Processing flow chart. 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,\ldots,x)\), where \(k_{0}\) is the main variable groundwater depth change rate. The mathematical expectation of the rate \(Z_{k_{0}}(x_{i})\) (\(i=1,2,3,\ldots,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}}}\int\limits_{V_{k_{0}}}Z_{k_{0}}(x)dx \tag{1}\]  
\(Z_{ak}\) is the average value on the carrier \(V_{ak}\):  
\[Z_{ak}=\frac{1}{V_{ak}}\int\limits_{V_{k_{0}}}Z_{k}(x)dx \tag{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}}}^{*}=\sum\limits_{k=1}^{K}\sum\limits_{ak}^{nk}\mu_{ak}Z_{ak} \tag{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_{ak}\) is the cokriging interpolation factor weight coefficient; \(Z_{ak}\) is a synergistic region variable.