The study area is the Shendong coal mining area located at the border between Shaanxi Province and Inner Mongolia Autonomous Region, China (1000-1500 m a.s.l; Fig. 1). The total area is approximately 3.214× 10³ km². The western and southern parts are encompassed by the Hobq Desert and the Mu Us Sandy Land, respectively (Bian et al., 2009a). Average annual precipitation is 360 mm and average annual evaporation is approximately 2300 mm. It should be noted that about 70% of the annual total precipitation occurs in August and September. The primary soil type is sand and the area is geographically characterized by an aeolian landform with sparse vegetation (Lei et al., 2010, 2014).

Fig. 1

	Figure Option ViewDownloadNew Window
	Fig. 1 Location of the study area in context of China (a) and in context of Shaanxi Province and Inner Mongolia Autonomous Region (b), and the borders of Shendong coal mining area and Daliuta Coal Mine (c)

Shendong coal mining area is affiliated to the Shendong Coal Company. There are a total of 15 coal mines distributed in the area (such as Bulianta Coal Mine, Daliuta Coal Mine, Huojitu Coal Mine, etc.), with the total production capacity of over 200× 10⁶ t/a. It should be mentioned that Daliuta Coal Mine, located in Shenmu County of Shaanxi Province, is the first underground coal mine built in the Shendong coal mining area (in 1986). It has a production capacity of 20× 10⁶ t/a with a designed service life of 119 years.

The Landsat series of satellites provides the longest continuous record of satellite-based observations (Qin et al., 2004; Sobrino et al., 2008; Chander et al., 2009; Jimé nez-Muñ oz et al., 2009). As such, Landsat is an invaluable resource for monitoring environmental changes (Fuller et al., 1994; Vogelman et al., 2001; Goward et al., 2006; Wulder et al., 2008). In this study, seven Landsat 5 (or 8) TM/OLI images spanning from 2001 to 2015 were selected (Table 1). By interpreting the satellite images, we classified land-cover types in the study area into four categories: bare land, bare rock/impervious surface, vegetated land, and water body. We then calculated the vegetation coverage of each land-cover type (Table 1).

Table 1

Table 1

Table 1 Characteristics of the selected seven Landsat 5 (or 8) images Scene Landsat scene ID Acquisition date Path/Row Vegetation coverage (%) Bare land and bare rock/impervious surface Vegetated land Water body 1 LT51270332001111BJC00 21 Apr, 2001 127/33 86.04 11.68 2.28 2 LT51270332003101BJC00 11 Apr, 2003 127/33 91.55 7.69 0.76 3 LT51270332005122BJC00 2 May, 2005 127/33 82.72 16.47 0.81 4 LT51270332007096BJC00 6 Apr, 2007 127/33 83.13 16.32 0.55 5 LT51270332009309BJC00 5 Nov, 2009 127/33 92.59 6.93 0.54 6 LC81270332014307LGN00 3 Nov, 2014 127/33 83.29 16.03 0.67 7 LC81270332015070LGN00 11 Mar, 2015 127/33 88.81 10.26 0.94

Table 1 Characteristics of the selected seven Landsat 5 (or 8) images

To validate the estimations of SMC using TM/OLI images, we conducted in situ measurements of SMC in the Daliuta Coal Mine on 26 September, 2013. A portable WET-2 sensor (Delta-T Devices, UK) was used for measuring the volumetric SMC. The positions of the sampling sites were recorded with a GPS receiver.

The hourly air temperature data were obtained from the China Meteorological Forcing Dataset (http://westdc.westgis.ac.cn/data/) and the adjacent meteorological stations around the study area. It should be noted that the hourly air temperature data from 2013 to 2015 were collected from the Yulin Meteorological Station and the Dongsheng Meteorological Station, because China Meteorological Forcing Dataset did not extend beyond 2012. Since the time interval of hourly air temperature data from the China Meteorological Forcing Dataset was 3 hours, we used the cubic spline interpolation (Borak and Jasinski, 2009) to scale the time interval down from 3 hours to 1 hour. Yearly precipitation and air temperature data provided by the Shendong Coal Company were obtained from local observation in Ulan Moron Town.

All visible light wave bands of the Landsat data were firstly resampled to a resolution of 30 m. As the atmospheric scattering effects should be removed prior to analysis to extract land cover information efficiently from Landsat data, we used the method of FLAASH (Fast Line-of-sight Atmospheric Analysis of Hypercubes) atmospheric correction to correct the visible light wave bands. FLAASH contains the MODTRAN 4+ radiation transmission model that has a high precision (Yuan et al., 2009). The Mono-window algorithm (Qin et al., 2001) was used for land surface temperature inversion. The influence of the atmosphere was considered in the process of calculation, therefore, atmospheric correction was unnecessary for the thermal infrared bands (Qin et al., 2003). Geometric calibration of Landsat images was then performed for the bands on scene 2 to scene 10 using scene 1 as a benchmark to ensure that the overall error would be controlled within 0.5 pixels. Then, the vector boundary of the study area was overlaid on all images to mask the study area.

For the thermal inertia model, it is critical to obtain the soil temperature difference. Only a few satellite sensors can provide the diurnal soil temperature difference, such as the MODIS sensor and NOAA/AVHRR sensor (Verstraeten et al., 2006). However, these sensors have a moderate spatial resolution (1 km). The Landsat TM/OLI images have higher spatial resolution (30 m) than the MODIS and NOAA/AVHRR images, but they cannot provide the day-night soil temperature difference (Lei et al., 2014). To solve this problem, we developed an improved thermal inertia model by integrating the land surface temperature derived from the TM/OLI images and a simple sine model to obtain the soil temperature difference.

2.3.1 Thermal inertia model

The thermal inertia (TI) can be calculated using the thermal inertia model developed by Cai et al. (2007) that is based on the MODIS images. The solution of TI was used in the first-order approximation of the Fourier series. It should be noted that following the real thermal inertia model of Xue and Cracknell (1995a), the thermal inertia model of Cai et al. (2007) was modified to be suitable in situations where either the satellite overpass time coincides with the local maximum and minimum temperature time or not.

The thermal inertia (TI) was calculated as follows (Cai et al., 2007):

$TI=\frac{(1-A){{S}_{0}}{{C}_{t}}{{A}_{1}}[\cos (\omega {{t}_{2}}-{{\delta }_{1}})-\cos (\omega {{t}_{1}}-{{\delta }_{1}})]}{\Delta T\sqrt{\omega }\sqrt{1+\frac{1}{b}+\frac{1}{{{b}^{2}}}}}, $ (1)

$b=\frac{\tan (\omega {{t}_{\max }})}{1-\tan (\omega {{t}_{\max }})}, $ (2)

${{\delta }_{1}}=\frac{b}{1+b}, $ (3)

where, A is the surface albedo; S₀, the solar constant (S₀=1367 W/m²); C_t, the atmospheric transmittance in the visible spectrum, typically equal to 0.75 (Xue and Cracknell, 1995b); A₁, the coefficient of the Fourier series, which can be derived from solar declination and local latitude; ω (rad/s), the angular velocity of Earth’ s rotation; t₁ and t₂(s), the time of the MODIS satellite overpassing at day and night, respectively; δ ₁, the phase difference in the first-order approximation of the Fourier series, which can be expressed as a function of b (Eq. 3); ∆ T(K), the land surface temperature difference during time t₁ and t₂; b, a process variable, which can be calculated from Equation 2; and t_max(s), the time of maximum temperature during the day.

2.3.2 Surface albedo

Before the determination of surface albedo, spectral reflectance in short wave radiation should be obtained (Cai et al., 2007). For TM/OLI images, the spectral reflectance in each band can be acquired using the FLAASH model. The broadband surface albedo (A) is the weighting function for converting spectral reflectance and can be calculated from Equation 4 (Liang, 2001):

$A=0.356{{\rho }_{ch.1}}+0.130{{\rho }_{ch.3}}+0.373{{\rho }_{ch.4}}+0.085{{\rho }_{ch.5}}+0.072{{\rho }_{ch.7}}-0.0018, $ (4)

For TM images, ρ _ch.1, ρ _ch.3, ρ _ch.4, ρ _ch.5 and ρ _ch.7 refer to the surface reflectance from bands 1, 3, 4, 5 and 7, respectively; and for OLI images, ρ _ch.1, ρ _ch.3, ρ _ch.4, ρ _ch.5 and ρ _ch.7 refer to the surface reflectance from bands 2, 4, 5, 6 and 8, respectively.

2.3.3 Diurnal land surface temperature difference (∆ T)

In the study of Cai et al. (2007), ∆ T is the diurnal land surface temperature difference during the MODIS satellite overpass time t₁ and t₂ at day and night. However, TM/OLI images were used in this study, we thus introduced a new method to obtain ∆ T and t_max.

A simple sine model that describes the diurnal land surface temperature difference was established in this study. It has been widely reported that the diurnal cycles of land surface temperature and air temperature can be approximately described as sine curves (Xue and Cracknell, 1995a; Campbell and Norman, 1998; Cai, 2006; Stisen et al., 2007; Jiang et al., 2010; Ma et al., 2013). Therefore, we assumed that both diurnal air temperature and land surface temperature are the sine functions of time t.

${{T}_{\text{air}}}={{M}_{1}}\sin ({{a}_{1}}t+{{b}_{1}})+{{c}_{1}}, $ (5)

${{T}_{ls}}={{M}_{2}}\sin ({{a}_{2}}t+{{b}_{2}})+{{c}_{2}}, $ (6)

where, T_air (K) is the air temperature at time t; T_ls(K) is the land surface temperature at time t; M₁ and M₂ (K) are the amplitude of the sinusoidal components; a₁, a₂, b₁ and b₂ are the parameters to calculate the phase; and c₁ and c₂ (K) are the offsets.

According to above assumption, we obtained ∆ Tand t_maxthrough the following three steps.

First, the air temperature sine model (Eq. 5) was fitted by inputting air temperature data obtained from the local meteorological stations and the China Meteorological Forcing Dataset, and then the parameters a₁, b₁, andc₁were obtained (Fig. 2).

Fig. 2

	Figure Option ViewDownloadNew Window
	Fig. 2 Observed hourly data of air temperature and fi tted air temperature curve after sine function (11 November, 2015)

Second, the land surface temperature sine model (Eq. 6) was derived from the fitted air temperature sine model. The prerequisite for estimating diurnal land surface temperature is that the land surface temperature and air temperature are not out of phase. The prerequisite is fulfilled from the morning to the midday. It was reported that the phase shift only occurs approximately 2.5 h after solar noon (Rosenberg et al., 1974). In addition, the diurnal minimum air temperature and land surface temperature are rather close (Stisen et al., 2007), implying that the offset of the air temperature sine model (Eq. 5) can be approximately substituted for the land surface temperature. Therefore, the parameters a₂, b₂, and c₂ are approximately and respectively equal to the parameters a₁, b₁, and c₁, which have been restricted to the time period of 00:00-14:30 (Beijing time). Thus, Equation 6 can be rewritten as Equation 7:

${{T}_{ls}}={{M}_{2}}\sin ({{a}_{1}}t+{{b}_{1}})+{{c}_{1}}.$ (7)

Since the T_ls at the time of TM/OLI image acquisition can be obtained from the land surface temperature inversion, the amplitude (M₂) of land surface temperature sine model was calculated from Equation 7.

Third, the ∆ Tand t_max were obtained from Equation 7. The value of t_max can be directly obtained from Equation 7. When the t₁ and t₂ are known, the value of ∆ Tcan be calculated from Equation 7 because ∆ T is equivalent to T_ls at t₁ minus T_ls att₂. It should be noted again that the model developed by Cai et al. (2007) is suitable in situations where either the satellite overpass time coincides with the local maximum and minimum temperature time or not. Thus, we consistently set t₁ at 11:00 (the general overpass time of Landsat) and t₂ at 02:00 (the general time of minimum land surface temperature) for all the calculation data to reduce errors.

2.4.1 SMC estimation

The relationship between thermal inertia and SMC in the M model is as follows (Ma, 1997):

$TI={{\left. {{\left\{ 2.1{{d}_{s}} \right.}^{(1.2-0.02{{d}_{P}}{{\omega }_{w, s}})}}\times \exp [-0.007{{({{d}_{P}}{{\omega }_{w, s}}-20)}^{2}}]+{{d}_{s}}^{(0.8+0.02{{d}_{P}}{{\omega }_{w, s}})} \right\}}^{\frac{1}{2}}}\times (0.2+\frac{{{d}_{P}}{{\omega }_{w, s}}}{100}){{d}_{s}}\times \sqrt{\frac{1}{1000}}, $ (8)

where, d_s (g/cm³) is the soil particle density; d_p is the relative soil particle density, which is equal to the average soil density divided by the water density, i.e., it is a dimensionless physical constant (2.65); and ω _{w, s} (%) is the weight percentage of SMC. Then, a lookup table can be established, which is composed of thermal inertia, soil particle density and SMC. Finally, the SMC can be estimated from thermal inertia by means of the lookup table. It should be particularly noted that the improved thermal inertia model is only suitable for bare soils or sparsely vegetated areas.

2.4.2 SMC validation

SMC was in suit measured by a portable WET-2 sensor at the No. 52303 Mining Face in the Daliuta Coal Mine. To minimize the impacts of land topography and geomorphology, we selected ten sampling sites with uniform slope (about 0.30-0.35) and with similar soil types and vegetation coverages. To avoid the possible disturbances of mining subsidence, we selected 5 sampling sites in the subsidence basin (No. 1-5 sampling sites in Table 2) and 5 sampling sites in the un-subsidence area (No. 6-10 sampling sites in Table 2). The distance between two adjacent sampling locations was about 30 m. SMC was measured at 5-cm depth with three replicates in each sampling site.

Table 2

Table 2

Table 2 Comparisons between measured soil moisture content (SMC) data and estimated SMC data from Landsat TM/OLI images Sampling site No. Location Soil particle density (g/cm3) Measured SMC (%) Estimated SMC (%) Difference (%) Latitude Longitude 1 39° 17′ 51′ ′ N 110° 18′ 59′ ′ E 2.70 12.51 12.01 0.50 2 39° 17′ 51′ ′ N 110° 19′ 02′ ′ E 2.69 5.12 10.09 -4.97 3 39° 17′ 50′ ′ N 110° 19′ 05′ ′ E 2.68 10.50 12.19 -1.69 4 39° 17′ 49′ ′ N 110° 19′ 06′ ′ E 2.68 5.55 10.42 -4.87 5 39° 17′ 49′ ′ N 110° 19′ 08′ ′ E 2.70 20.09 15.63 4.46 6 39° 17′ 50′ ′ N 110° 18′ 59′ ′ E 2.70 11.21 12.10 -0.89 7 39° 17′ 49′ ′ N 110° 19′ 01′ ′ E 2.68 10.50 12.40 -1.90 8 39° 17′ 49′ ′ N 110° 19′ 02′ ′ E 2.69 16.21 14.09 2.12 9 39° 17′ 48′ ′ N 110° 19′ 05′ ′ E 2.69 16.74 13.45 3.29 10 39° 17′ 48′ ′ N 110° 19′ 06′ ′ E 2.69 12.80 12.32 0.48

Table 2 Comparisons between measured soil moisture content (SMC) data and estimated SMC data from Landsat TM/OLI images

We first established the statistical relationship between the modelled data using the improved thermal inertia model and the in situ SMC data. Again, it should be noted that the in situ SMC measurements were conducted at approximately the same time as the image was taken (i.e., on 29 September, 2013) for validation, i.e., on 26 September, 2013. The estimated data from the RS-derived thermal inertia were gravimetric SMC values, which are different from the in situvolumetric SMC values. Thus, the estimated SMC data should be inverted into the volumetric values by multiplication of soil particle density (see Table 2).

Regression analysis was used for the validation (data not shown). Generally speaking, the estimated SMC agreed well with thein situ SMC (with R² value of 0.9398), indicating that the improved thermal inertia model performed rather well. Furthermore, the average difference betweenin situ SMC and estimated SMC from the TM/OLI images was 2.52% (Table 2), being lower than the result from the MODIS images of Cai et al. (2007).

In this study, all statistical analyses were performed using Microsoft Excel 2010 (Microsoft Inc., Redmond, Washington, USA). The simple linear regression method and the Pearson correlation coefficient were both used to describe the relationship between SMC and influencing factors (e.g., precipitation and air temperature). Difference was considered significant at P< 0.05 level. Furthermore, all figures were plotted using Origin 9.1 software.

3.1.1 Spatial variations of SMC at Testing Ground I

To detect whether underground mining is the dominant factor impacting the SMC at Testing Ground I, we compared the difference of mean SMC between mined area and un-mined area from 2001 to 2015. As shown in Table 3, mean SMC values in the mined area were slightly lower than those in the un-mined area during the period from 2001 to 2015.

Table 3

Table 3

Table 3 Comparison of mean SMC between mined area and un-mined area at Testing Ground I Date Mean SMC in the mined area (%) Mean SMC in the un-mined area (%) Difference (%) 21 Apr, 2001 7.6795 7.7840 -0.1045 11 Apr, 2003 7.9582 8.0628 -0.1046 2 May, 2005 8.0331 8.1605 -0.1274 6 Apr, 2007 9.8082 9.9383 -0.1301 5 Nov, 2009 10.9331 11.0605 -0.1274 3 Nov, 2014 12.6614 12.8005 -0.1391 11 Mar, 2015 12.7624 12.8991 -0.1367 Note: Testing Ground I refers to the Shendong coal mining area testing ground.

Table 3 Comparison of mean SMC between mined area and un-mined area at Testing Ground I

3.1.2 Spatial variations of SMC at Testing Ground II

In this study, we selected the No. 22617 Mining Face in the Daliuta Coal Mine to delineate the spatial variations of SMC at Testing Ground II. Furthermore, spatial distributions of SMC in the pre-mining period and post-mining period were analyzed. Because No. 22617 Mining Face was active in 2011, SMC values in 2009 were selected as the pre-mining data and SMC values in 2014 were selected as the post-mining data. Consequently, we calculated the mining disturbance zones of No. 22617 Mining Face based on the theory of mining subsidence (Guo and Chai, 2013). As shown in Figure 4, the mining disturbance zones at No. 22617 Mining Face consisted of the neutral zone, compression zone, internal tensile zone, and external tensile zone. To analyze the effect of these mining disturbance zones, we set the mining face as the boundary using ArcGIS software. Finally, 3 buffer zones were established outward and 4 buffer zones were established inward at an interval of 20 m. Each buffer zone approximately corresponded to a specific mining disturbance zone (Fig. 4).

Fig. 4

	Figure Option ViewDownloadNew Window
	Fig. 4 Buffer zones of No. 22617 Mining Face in the Daliuta Coal Mine (a) and an overall view of subsidence at No. 22617 Mining Face in the Daliuta Coal Mine (b). Numbers 1-7 in the lower figure represent the 1st-7th layer of buffer zones in the upper figure, respectively.

Average SMC values in each layer of buffer zone at No. 22617 Mining Face in the Daliuta Coal Mine in the pre-mining period (taken 2009 as the representative) and post-mining period (taken 2014 as the representative) are shown in Figure 5. Average SMC values showed an irregular trend in the pre-mining period and an increasing trend in the post-mining period along the layers of buffer zones. Due to the underground mining activities, the SMC clearly decreased from the un-mined zone (i.e., the 7^th layer of buffer zone) to the neutral zone (i.e., the 1^st layer of buffer zone). These results indicated that SMC was significantly impacted by underground mining activities, not only inside the mined area but also outside the mined area (near the boundary of the mined area in the external tensile zone).

Fig. 5

	Figure Option ViewDownloadNew Window
	Fig. 5 Average SMC at each layer of buffer zone at No. 22617 Mining Face in the Daliuta Coal Mine on 5 November, 2009 and 3 November, 2014. It should be noted that 2009 was the representative year of the pre-mining period, while 2014 was the representative year of the post-mining period.

Inter-annual variations of area-averaged SMC in the mined area and un-mined area at two testing grounds (Testing Ground I and Testing Ground II) from 2001 to 2015 are shown in Figure 6a. Values of area-averaged SMC at Testing Ground I and Testing Ground II were nearly identical. During the entire study period (2001-2015), a pronounced increasing trend of area-averaged SMC was observed. The difference of area-averaged SMC between mined area and un-mined area at the two testing grounds could be used to examine the impacts of underground mining activities on SMC. It can be seen from Figure 6b that with increased underground mining activities, the absolute value of the difference in area-averaged SMC between mined area and un-mined area gradually increased during the entire period. The results revealed that underground mining activities have significant impacts on SMC, especially at a mining face (see section 3.1.2).

Fig. 6

	Figure Option ViewDownloadNew Window
	Fig. 6 (a) Inter-annual variations of area-averaged SMC in the mined area and un-mined area at Testing Ground I and Testing Ground II during the period of 2001-2015 and (b) difference of SMC between mined area and un-mined area at Testing Ground I and Testing Ground II during the period of 2001-2015

To estimate long-term trend variation of SMC, we calculated the linear trend for every pixel of SMC using the linear regression model (Wei, 2007). The simple linear regression model was built with the parameters of time and SMC. The linear trend of SMC from 2001 to 2015 was calculated using the generalized linear equation (Wei, 2007):

$SMC=slope\times year+b, $ (9)

$Slope=\frac{\sum\nolimits_{i=1}^{n}{SM{{C}_{i}}{{T}_{i}}}-\left( \sum\nolimits_{i=1}^{n}{SM{{C}_{i}}} \right)\times \frac{\left( \sum\nolimits_{i=1}^{n}{{{T}_{i}}} \right)}{n}}{\sum\nolimits_{i=1}^{n}{T_{i}^{2}}-\frac{{{\left( \sum\nolimits_{i=1}^{n}{{{T}_{i}}} \right)}^{2}}}{n}}, $ (10)

where, slope (%/a) is the tendency value of SMC during the period of 2001-2015; b (%) is the intercept; n is the number of years (equal to 8 here); T_i is the i^th year beginning from 2001; and SMC_i (%) is the soil moisture content in the i^th year. The value of slope represents the trend of SMC during the period of 2001-2015. The positive slope values represent increases of SMC (wetting trend) while the negative slope values represent decreases of SMC (drying trend). The spatial patterns of trends in SMC during the period of 2001-2015 are shown in Figures 7a-1 and a-2. There was a significant spatial heterogeneity of SMC throughout the study area during the period of 2001-2015. Overall, areas with signifi cant decreasing trend of SMC were distributed primarily in the surrounding region of the water body and also in the grasslands of the southeastern and northeastern regions.

Fig. 7

	Figure Option ViewDownloadNew Window
	Fig. 7 Spatial patterns of trend in SMC from 2001 to 2015 (a-1 and a-2) and percentages of area with wetting trend and drying trend (b) for Testing Ground I and Testing Ground II. The positive slope values represent increases of SMC (wetting trend) while the negative slope values represent decreases of SMC (drying trend).

To further detect the impacts of underground mining activities on SMC during the period of 2001-2015, we calculated the proportional areas of positive and negative slope pixels (Fig. 7b). For Testing Ground I, the wetting trend was dominant both in the mined area and un-mined area, and the percentage of area with wetting trend was higher in the un-mined area than in the mined area. However, the percentage of area with drying trend was 3.09% higher in the mined area than in the un-mined area. For Testing Ground II, the percentage of area with drying trend in the mined area was approximately 1.69-fold higher than that in the un-mined area.