Coal is predicted to continue to play an indispensable role in China’s energy structure in the future (Wang 2015). The Shen Dong Coalfield is located in the arid and semi-arid region of western China, and its coal resource reserves account for approximately one-fourth to one-third of China's total coal reserves (Wang 2017). However, this area is located in the transition zone between the Loess Plateau and the Mu Us Desert. The surface ecology, which is eroded by water and wind, is a typical ecologically fragile area (Wang et al. 2008; Zhang et al. 2012). Natural conditions of drought and water shortages make soil moisture (SM) a key factor in restricting vegetation growth and restoration of the ecological environment (Dougill et al. 1998; Bi et al. 2014). High-intensity mining causes the movement and deformation of rock formations, leading to the formation of subsidence areas and cracks, which in turn can affect groundwater, soil, and vegetation from the bottom to top (Bi et al. 2014; Wang 2017). The first result of the formation of subsidence areas is a change in SM (Zhao 2006). Changes in SM in subsidence areas lead to changes in nutrients such as nitrogen, phosphorus, potassium, and organic matter (Li et al. 2001; Cheng et al. 2014; Liu et al. 2014). These changes in turn affect the distribution and circulation of original substances in the ecosystem, thereby impacting the ecological environment. Therefore, it is of great significance to elucidate the mechanisms that affect SM in coal mining subsidence areas in order to develop effective programs for vegetation restoration and ecological reconstruction.

Some scholars believe that cracks caused by coal mining cause decreases in SM in subsidence areas (Wei et al. 2006; Wang et al. 2011; Liu et al. 2016; Guo et al. 2019). The depth, density, width, as well as distance from the cracks, all influence the changes in SM (Zou et al. 2014; Zhang et al. 2015; Ma and Yang 2019; Wu et al. 2019, 2020). Zhang et al. (2015) tracked and monitored the cracks in a fixed position in the Bu Lian Ta mine in the Shenfu-Dongsheng coalfield and found that the SM was lower in the fractured and non-fractured portions of the subsidence area than in virgin areas. Moreover, SM around ground cracks was significantly lower, with a significant downward trend over the depth range of 0–90 cm, ranging from 9.27% to 15.47%; Wu et al. (2019) studied the Fuxin subsidence area and showed that SM increased with increasing distance from cracks, but the effect was not notable at distances beyond 2 m. Other scholars believe that coal mining does not cause decreases in SM. These researchers contend that SM is supplied by atmospheric precipitation rather than groundwater, so the impact of coal mining–related disturbances on SM is very limited (Li et al. 2012; Wang et al. 2006). For example, the "Ecological restoration technology experiment and demonstration research in Shen Dong mining subsidence area" report, completed by Inner Mongolia Agricultural University and China University of Mining and Technology in 2007, showed that coal mining subsidence does not affect soil nutrients and moisture (Li et al. 2012). However, the abovementioned research only examined the effect of coal mining on SM from the perspective of the absolute value of SM, ignoring characteristics of the soil itself as a heterogeneous continuum. Even if the soil is affected by the same disturbance, the SM will differ due to differences in spatial locations; that is, SM exhibits spatial variability. At present, scholars in China and abroad focus primarily on the change in the absolute value of SM in a given subsidence area, and there has been little research on spatial variations in SM in aeration zones in subsidence areas.

The present study examined the Nalin River No. 2 mine in the eastern Mu Us Desert. A combined classical statistics and multi-dimensional geo-statistics approach was employed to analyze the variability of SM in the aeration zone in the subsidence area both temporally and spatially. The impact of coal mining subsidence on SM in the aeration zone is discussed from the perspective of spatial variability rather than that of absolute SM changes, so as to provide a scientific basis for ecological restoration efforts in mining areas.

2.1 Overview of the study area

The Nalin River No. 2 mine is located in the ecologically fragile northwestern part of China in the eastern Mu Us Desert (108° 51′ 30″ E–109° 00′ 00″ E, 37° 58′ 00″ N–38° 05′ 30″ N). The geographical location of the mine is shown in Fig. 1. The region has an arid and semi-arid continental climate, dry and rainless, windy and sandy. The average annual precipitation at the study site is 350 mm, and the average annual evaporation is 2500 mm. Additionally, the region is located in the semi-desert area of the plateau desert landform, with sparse and scattered vegetation and crescent or wavy dunes on the surface. The soil types are primarily fixed and semi-fixed aeolian sandy soils, and the vegetation consists primarily of xerophytic and semi-xerophytic sandy types. Two seasonal rivers drain the area, and the depth of the phreatic water is 19.15–25.65 m. The average depth of the main buried coal seam is approximately 602 m, and the average thickness of the seam is approximately 4.65 m. The mining method is longwall combined.

Schematic diagram of the geographic location of the study region

2.2 Collection and determination of soil samples

Based on the time of the mining completion, the study area was divided into a control area (CK), 1-year subsidence area (S1), and 2-year subsidence area (S2) by partition sampling. Due to the underground coal mining process, a checkerboard distribution method was adopted. The grid design was 75 m × 50 m, and each sampling point was in the center of the grid. The focus was on the subsidence basin areas and edge fracture development areas, such as the edge and center line of the working face. A total of 100 soil sample collection points were designated in three areas, and among these points, 20 were set in the control area, and 40 points were respectively designated in the 1-year subsidence area and 2-year settlement area, as shown in Fig. 2. Using a soil drill, 0–10 m soil samples were collected at each sampling point, with an interval of 1 m. Three composite samples were randomly collected from each layer of soil depth. After natural drying, the samples were sieved using 2-mm mesh for SM measurement. SM was measured by gravimetry, in which each fresh soil sample was dried at 105 °C for 12 h; SM was determined by dividing the difference between the fresh weight and dry weight of the soil by the dry weight (Zou et al. 2014).

Layout of sample points in the different areas

2.3 Data analysis

2.3.1 Classical statistical methods

SPSS 20.0 data analysis tool was used to conduct normal distribution K-S testing on all data (p = 0.05), and the mean, median, standard deviation (SD), coefficient of variation (CV), maximum value, and minimum value were calculated. Mean and median values reflect the centralized trend of samples, whereas the positional relationship between the mean and median values reflects the relationship between data distribution and outliers. The SD and CV reflect the variability characteristics of the data as a means of measuring the degree of data dispersion (Wang 2017; Zhang et al. 2012). It is generally recognized that a CV ≥ 100% indicates high variation, whereas 10% < CV < 100% indicates moderate variation, and CV ≤ 10% indicates low variation (Li et al. 2021).

2.3.2 Geo-statistical methods

Geo-statistical methods can be adopted to study the spatial distribution of heavy metals, nutrients, and moisture in soil (Shen 1989; Pietrzykowski and Chodak 2014; Jing et al. 2018). Semi-variance functions can be used to determine the spatial autocorrelation of variables, and a Kriging interpolation method can also be used to simulate and estimate the regional content around the sampling point.

A semi-variance function, also known as a variogram, is used to quantify the randomness and spatial structure of variables according to Eq. (1):

In this formula, γ(h) is the variogram; h is the spatial distance between two variables; N(h) is the number of sample point pairs; Z(X_i) is the observation value of the spatial position point X_i; and Z(X_i + h) is the observation value of the point h away from X_i (i = 1,2, N(h)). Commonly used semi-variograms include Gaussian, spherical and exponential models:

In these formulas, α is the range; C₀ is the nugget constant, which represents the variance caused by random error; C is the space structure value caused by systematic factors; C₀ + C is the abutment value, representing the total variance of variables; and C₀/ (C₀ + C) is the spatial structure ratio. If C₀/(C₀ + C) < 25%, it can be shown that the data exhibit a strong spatial correlation; when 25% < C₀/(C₀ + C) < 75%, it can be shown that the data exhibit a medium spatial correlation; and when C₀/(C₀ + C) > 75%, it can be shown that the data exhibit a weak spatial correlation (Cambardella et al. 1994).

The standard for selecting the semi-variance function model is that the closer to zero is the mean absolute error (MAE) and the root mean square error (RMSE) of the cross-check result, the closer the regression fitting coefficient of determination R² is to 1, and the greater will be the accuracy of the model's simulation.

The Kriging interpolation method is based on the spatial autocorrelation, using the original data of the regionalized variables and the structure of the variogram, with the adoption of linear, unbiased, and optimal interpolation methods for the unknown sampling points of the regionalized variables. The formula is as follows:

In this formula, Z(x₀) is the value of the point to be estimated; n is the number of sampling points; Z(x_i) is the value of the ith sampling point; λ_i is a group of weight coefficients; ∑λ_i = 1; and the selection of λ_i ensures that the estimation of Z(x₀) is unbiased and the estimation variance minimized.

2.3.3 Principal component analysis

Principal component analysis (PCA) is a multivariate statistical analysis method that linearly transforms multiple variables as a means of selecting important variables (Zhang et al. 2017; Chahal and Eerd 2019; Zuber et al. 2017). The mathematical model of a PCA is:

where m₁, m₂, …, m_P represent ‘P’ principal components. Using PCA, the data should be standardized, the correlation coefficient between variables should be calculated, and the eigenvalue, eigenvector, contribution rate, and cumulative contribution rate should be obtained. In this study, PCA was used to identify the major factors of coal mining that affect SM.

3.1 Statistical characteristics of SM

Descriptive statistics of SM from 0 to 10 m in the CK, S1, and S2 regions are summarized in Table 1. As shown in the table, (1) Both the mean and median of the overall SM of the CK area were higher compared with S1 and S2 (average value increased by 1.18% and 0.96%, respectively). Moreover, the mean soil moisture at each depth level was very close to the median value, indicating that the centralized distribution of the data was not dominated by outliers. (2) The trend in SM variation in the vertical direction in different regions was basically the same. SM exhibited a downward trend above 2 m and reached the lowest value in the 1–2-m soil layer (4.49%, 3.07%, and 3.02%, respectively). SM, however, increased from top to bottom with increasing soil depth below 2 m and reached a maximum value in the 9–10-m soil layer (10.19%, 8.49%, and 8.63%, respectively).

According to the K-S normal distribution test, the 0–10 m soil moisture data in the three regions exhibited a normal distribution, and the overall trend in variation coefficients was S1 > S2 > CK, with an average CV of 43.60%, 29.31%, and 25.68%, respectively. Compared with CK, the coefficients of variation of S1 and S2 were 17.92% and 3.63% higher, respectively; that is, the degree of spatial variation in SM was highest in the 1-year subsidence area.

Compared with the CV for the same soil depth in different regions, the spatial variability of CK increased at depths of 0–3 m and 4–7 m. The variability of S1 gradually increased in the range 0–5 m and 7–9 m, with the highest variability at 3–4 m and 8–9 m. S2 showed an increasing trend overall, but the variability was in between that of CK and S1, indicating that coal mining and other factors have a particular impact on SM in the 2-year subsidence area. However, as the subsidence stabilizes, the area itself may undergo certain self-repairs (Li et al. 2012) that relatively reduce the degree of variation. The above data show that the trend in variation of the subsidence and control areas was discrepant at different soil depths due to mining activities.

In summary, the trend in variation of SM in the vertical direction in each area was basically similar. However, due to the different effects of coal mining activities and other factors on the soil at different depths, there were some differences in the degree of variation. For example, the degree of variation in SM at 0–1 m and 5–10 m in the control area was significantly lower than that in the subsidence area, whereas there was no consistent trend in the 2–5 m SM variation, indicating that the surface subsidence increased the degree of variation in SM. With the subsidence entering the stable stage, the degree of variation degree declined relatively, but in the 0–1-m soil surface layer and 5–10-m deep soil layer, the degree of variation was still greater than that in the control area.

3.2 Spatial structure change in SM

Table 2 shows the best-fitting models of CK, S1, and S2 in the 0–10 m soil layers. C₀ is the variation caused by experimental errors and the microdomain structure smaller than the actual sampling scale, reflecting the size of random variation C (Cambardella et al. 1994). From Table 3, a vertical comparison of C₀ values in the subsidence area shows that the random variation in deep soil was greater than that of surface soil, which may have been caused by the difficulty associated with sampling deep soil.

In terms of C₀ + C, all soil layers in the control area were lower than those in the subsidence area, indicating that the variation in SM in the control area was small and within the range of variation. However, considering the trend in the variation of C₀ in the three areas, the C values (structural variance, representing variation from non-random causes in the subsidence area) at 0–1 m and 5–10 m were relatively larger than those in the CK, indicating that certain structural factors, including climate, coal mining subsidence, and soil texture, have a greater impact on the subsidence area.

The spatial structure ratio C₀/(C₀ + C) represents the proportion of system variation among the total variation. According to the standard reported by Cambardella et al. (1994), CK exhibited a strong spatial autocorrelation between 0–3 m and 6–8 m (C₀/(C₀ + C) < 25%), with moderate spatial autocorrelation at 3–6 m and 8–10 m (25% < C₀/(C₀ + C) < 75%). S1 had a strong spatial correlation only at 1–4 m and moderate or weak spatial correlations at other depths. S2 had a strong spatial correlation only at 1–5 m and moderate or weak spatial correlations at other depths. In particular, at a depth of 5–10 m, the C₀/(C₀ + C) values of the 1-year and 2-year subsidence areas were > 50%, and the degree of spatial variation was high. In conclusion, the SM in the 1-year and 2-year subsidence areas showed strong or moderate spatial variability at the soil surface (0–1 m) and in deep soil (5–10 m), consistent with the variation trend in the classical statistical results. The contribution of random factors to the spatial distribution of the soil surface and deep layer was shown to be relatively small, and the spatial variation was mainly caused by structural factors. The above results indicate that the overall degree of variation in the subsidence area was higher than that of the control area, and the variation caused by structural factors such as climate, mining subsidence, and soil texture has a greater impact on the subsidence area, particularly the spatial variability of surface and deep soil.

Cross-check was used to verify the accuracy of the interpolation results in the three regions, and the MAE, RMSE, and R² were calculated between the predicted and true values (Fig. 3). As indicated in Fig. 3, compared with the CK, the MAE and RMSE for S1 were 22.54% and 19.13% higher, respectively; for S2, the MAE and RMSE were 1.41% and 6.09% higher, respectively, compared with the CK but 17.24% and 10.95% lower, respectively, compared with S1. These results show that the accuracy of the Kriging interpolation for the control area was generally higher than that for the subsidence area, and the accuracy of interpolation for the subsidence area was the lowest in year 1; that is, as the spatial variability in SM increased, the accuracy of the Kriging interpolation decreased. The R² values for CK, S1, and S2 were 0.76, 0.59, and 0.71, respectively, similar to the spatial interpolation results of Zou (2014) for SM in the Shendong mining area, reflecting that the greater the spatial variability in SM, the lower the accuracy of fitting between the predicted and measured values.

Cross-check charts: a1 S1; a2 S2; a3 CK

3.3 Spatial distribution of SM

The above analyses demonstrated that soil moisture has great spatial variability at 0–1 m and 5–10 m (Fig. 4). Therefore, Kriging interpolation was used to draw a spatial distribution map of SM in the CK, S1, and S2 layers (0–1 m, 5–6 m, and 9–10 m).

Spatial variability in soil moisture (resolution ratio 2 m × 2 m)

Comparing the interpolation results for CK, S1, and S2 at 0–1 m, 5–6 m and 9–10 m, we found no significant spatial variability at 0–1 m, and there was one low-value area and two high-value areas of SM. The low-value area was located in the 1-year subsidence area, whereas the high-value area was located in the control area and 2-year subsidence area. This was primarily because the control area was not affected by the ground fissures caused by coal mining, whereas the 2-year subsidence area was gradually stabilizing, such that the SM was being restored. At 5–6 m, the high-value area was located primarily in the control area and 1-year subsidence area, but the SM distribution in the subsidence area had poor spatial continuity, indicating that the subsidence caused by coal mining has a particular impact on the SM distribution. At 9–10 m, the distribution of high-value and low-value areas of SM was similar to that of the 0–1-m soil layer. The high-value area was located primarily in the northwest part of the control area and the 2-year subsidence area. However, the low-value area and the high-value area in the 2-year subsidence area had poor spatial continuity, leading to high spatial variability. On the whole, the distribution of SM was such that the high-value area was located in the control area, followed by the 2-year subsidence area and 1-year subsidence area, demonstrating that surface subsidence has a greater impact on SM. Areas with obvious spatial variability in SM were primarily in the subsidence area, especially in the deep soil, because the control area had not been mined. The surface subsidence and ground fissures caused by disturbances associated with coal mining changed the soil structure and water transport channels, such that the SM decreased and the spatial variability increased.

Comparison of the results of SM interpolation of CK, S1, and S2 at 0–1 m, 5–6 m, and 9–10 m vertically revealed that the SM in the three regions decreased initially and then increased, consistent with the water change trend in the results of the classic statistical analysis. In the CK area, weak spatial variability was observed at each depth, and the vertical distribution of SM exhibited clear regularity, showing a decreasing trend from northwest to southeast. In S1 and S2, there was less SM in each depth layer than in the CK area, there was no clear regularity in the vertical distribution, and the degree of variation in SM at each depth layer was also greater than that in the CK area. The reason for this result is that coal mining–associated collapse causes the soil layer in the subsidence area to be inverted and reorganized vertically, such that it causes changes in soil texture, bulk density, porosity, and other physical properties of the soil, resulting in a significant decrease in the subsidence area in terms of SM compared with the control area. Moreover, a strong variability in the spatial distribution was detected.

By comparing the SM distribution in three regions vertically and horizontally, we found that the results of Kriging interpolation were generally consistent with the results of classical statistics and model fitting. The SM in the control area was higher than that in the subsidence area, whereas the degree of variation in the subsidence area was higher than that in the control area. We speculate that the primary reason for this difference is that the surface subsidence and ground fractures caused by mining-related disturbances altered the soil structure and soil physical properties. These changes in turn affected the water migration channels, leading to a decrease in SM and an increase in spatial variability.

4.1 Analysis of influence factors

As indicated in Sects. 3.1 and 3.2 above, the spatial variability in SM in the 1-year subsidence area was relatively high. Therefore, a principal component analysis of SM at 0–10 m in the 1-year subsidence area was conducted to identify the primary influencing factors. The results are shown in Tables 3 and 4. In this study, three principal component factors (Y1, Y2, and Y3) were obtained, and their corresponding variance contribution rates were 41.88%, 19.39%, and 14.84%, respectively. The total cumulative variance rate was 76.1%.

From an examination of the matrix after rotation, the following conclusions can be drawn regarding the main factors influencing the vertical distribution of SM in the 1-year subsidence area.

1. (1)

   The vertical distribution of 0–1-m SM is mainly affected by factor Y3. The factor load is − 0.80, which may be indicative of surface factors such as vertical cracks and subsidence, evaporation, and vegetation cover caused by coal mining subsidence. Through the comparative analysis of 0–60 cm SM in crack and non-crack areas in the Shigetai coal mine of the Shenfu Dongsheng coalfield, Guo et al. (2019) found that the low-value area of SM was located in a crack area, and the high-value area was located in the non-crack area covered by vegetation. By analyzing the SM change and associated influencing factors in the 0–60 cm coal mining subsidence area of a mine in Huaibei City, Anhui Province, China, Xu et al. (2018) reported that changes in precipitation and evaporation conditions exert a tremendous impact on the SM distribution in the aeration zone. The findings of the above studies demonstrate to a certain extent that cracks and subsidence caused by coal mining, vegetation cover, and evaporation affect the SM at 0–1 m.

2. (2)

   The vertical distribution of SM at 5–10 m was mainly affected by factor Y1, which may be due to the change in SM caused by rocks and soil deformation and the hydraulic connection between the aeration and saturated zones. The respective factor loads were 0.69, 0.82, 0.61, 0.80, and 0.82. In the process of studying the spatial variability of SM in the Daliuta coal mining subsidence area, Zhao et al. (2010) found that due to coal mining subsidence, the soil horizon in the subsidence area was inverted and reorganized vertically, resulting in a change in SM, which significantly reduced the SM in the subsidence area compared with that in the non-subsidence area. Zhao (2006) also reported that after a mine collapse in the area, the geotechnical structure of the vadose zone was damaged, which in turn altered the distribution and migration mechanism of water in the vadose zone accordingly. Therefore, the predicted SM influencing factors at 5–10 m in this study were consistent with the conclusions of Zhao et al.

3. (3)

   Factors Y1 and Y3 had minimal effect on the vertical distribution of SM in the range of 1–5 m, as the vertical distribution of SM in this area was primarily affected by factor Y2, which might be the water holding capacity and infiltration capacity of the soil itself, with factor loads of 0.92, 0.91, 0.92, and 0.76, respectively. Zhang et al. (2021) systematically examined the impact of coal mining subsidence on SM changes in arid and semi-arid areas and reported that high-intensity and large-scale underground mining in arid and semi-arid areas causes surface subsidence, destroys the original landform, and changes various soil physical properties such as soil particle size composition, porosity, compactness, permeability, and bulk density, thus affecting SM holding capacity. Zang et al. (2012) and Bi et al. (2014) and others have reported that the infiltration rate in coal mining subsidence areas is greater than that in non-subsidence areas, that the SM holding capacity in subsidence areas is decreased, and that the vertical infiltration depth of water is significantly increased. The above research shows that subsidence causes changes in soil texture, pore structure, and other properties, resulting in changes in SM holding capacity and infiltration capacity, thereby affecting SM.

In summary, cracks formed by coal mining activities and the resulting changes in soil porosity, texture, and hydraulic pressure lead to the formation of preferential channels for water transport that ultimately affect the distribution of water and its circulation in migration channels.

4.2 Effect of coal mining on SM

Underground coal mining leads to the formation of underground goafs and causes surface collapse. The tensile action associated with the process of collapse leads to the development of a large number of surface fissures that cause changes in the original texture and structure of the soil and the preferential flow characteristics. Therefore, the hydraulic conductivity of the soil as well as the evaporation area and intensity of soil water are increased. Concurrently, both SM and the field water holding capacity, which is the water holding capacity of the entire aerated zone, are reduced.

Preferential flow is the rapid and non-equilibrium seepage flow of water in the soil, which is a sign of the transformation from homogeneous to heterogeneous soil water movement (Meirvenne et al. 2003). Fissures are a common phenomenon in soil and exert complex effects on the generation and process of preferential flow. When there are no cracks, rainfall first supplements moisture loss in the upper soil and then infiltrates downward, primarily by uniform flow. When cracks appear on the surface, the rainfall infiltration mode changes from "piston" infiltration in non-collapse areas to a "shortcut" type of infiltration. The surface water infiltrates directly along large cracks and can even connect with the groundwater, changing the original soil water circulation path. After the preferential flow is generated, the fissures increase the surface area of water infiltration, which leads to an increase in the transport speed of water and solutes (Ma 2007).

The presence of fissures also increases the evaporation area and intensity of soil water, resulting in changes in SM. The study on the characteristics of soil priority flow in the Shen Dong subsidence area shows that the soil preferential flow in the coal mining subsidence area primarily involves soil macropore flow. The soil water infiltration process in the studied area exhibits non-equilibrium characteristics. The steady rate of infiltration and outflow of a soil layer under preferential flow is greater than that of a soil layer exhibiting uniform flow (Su 2018). Therefore, the change in soil surface water can be attributed to subsidence caused by coal mining collapse as well as the soil priority flow caused by cracks. In addition, deep SM is also indirectly affected by coal mining subsidence and cracks, resulting in preferential flow. Su (2018) pointed out that as soil depth increases, the spatial differentiation and variation in the preferential flow of the soil in a coal mining subsidence area continue to increase, and the spatial morphology of flow in areas where preferential flow occurs becomes more complex. After an underground coal seam is mined out or the ore body aquifer is discharged, the stress equilibrium of the rock mass around the goaf changes. These changes inevitably cause deformation, breakage, and movement of the rock layer (Zhu 2017). This results in a change in soil structure and texture due to increasing porosity and decreasing bulk density, which ultimately affects the water holding capacity and infiltration capacity of the soil and consequently the movement of soil water. Changes in soil texture and pore structure are related to the physical, mechanical, and hydraulic properties of the soil, as well as to SM in the tensile deformation stage (Cheng 2016). Thus, the changes in soil texture and pore structure caused by rock and soil deformation cause changes in the preferential flow of soil fissures, the soil water transport channel, and the spatial variability of SM.

Based on the above data, the mechanism by which coal mining affects SM can be summarized as follows (Fig. 5): the subsidence and cracks caused by underground coal mining alter the surface soil texture to a certain extent, causing preferential flow that in turn affects the evaporation, infiltration, and migration of soil water. However, due to the limited range of influence of subsidence and cracks, even though the average SM in the subsidence area decreases, there is no consistent change in the overall factors governing spatial distribution, resulting in an increase in spatial variability. In addition, changes in the soil pore structure and texture caused by rock and soil deformation resulting from the effects of coal mining alter soil hydraulic properties, significantly disturbing the water in the deep aeration zone of the soil and increasing the spatial variability of SM.

Schematic diagram of spatial variation in soil moisture in the coal mining subsidence area

In this study, the spatial variability of SM at 0–10 m in subsidence areas of differing ages was analyzed in comparison with a control area using classical statistics and multi-dimensional geo-statistics methods. The following conclusions were drawn:

1. (1)

   The SM in the subsidence area and control area initially increased, then decreased, and then increased at 0–10 m. Compared with the control area, the average SM in the subsidence areas of different ages decreased by 1.18% and 0.96%, and the coefficient of variation increased by 17.92% and 3.63%, respectively. That is, the SM in the control area was significantly higher than that in the subsidence area, but the variability in the subsidence area increased.

2. (2)

   There was no difference in the degree of variation in SM at different levels in the control area, and there was weak spatial variability. The degree of variation in SM in each layer of the 1-year subsidence area was quite different, with strong spatial variability. The spatial variability in SM in the 2-year subsidence area was between that of the control area and 1-year subsidence area.

3. (3)

   The accuracy of SM interpolation in the control area was greater than that in the subsidence area, especially in the 1-year subsidence area, with high variability. Compared with the cross-validation results for the 1-year and 2-year subsidence areas, the RMSE relative to the control area declined by 19.13% and 6.09%, respectively, whereas the R² increased by 0.17 and 0.05, respectively.

4. (4)

   Changes in soil texture and pore structure caused by surface subsidence and surface cracks resulting from coal mining alter the soil water infiltration from uniform flow to preferential flow, resulting in an increase in spatial variability of SM.