Changes of SOC Content in China’s Shendong Coal Mining Area during 1990–2020 Investigated Using Remote Sensing Techniques

: Coal mining, an important human activity, disturbs soil organic carbon (SOC) accumulation and decomposition, eventually affecting terrestrial carbon cycling and the sustainability of human society. However, changes of SOC content and their relation with influential factors in coal mining areas remained unclear. In the study, predictive models of SOC content were developed based on field sampling and Landsat images for different land-use types (grassland, forest, farmland, and bare land) of the largest coal mining area in China (i.e., Shendong). The established models were employed to estimate SOC content across the Shendong mining area during 1990–2020, followed by an investigation into the impacts of climate change and human disturbance on SOC content by a Geo-detector. Results showed that the models produced satisfactory results (R 2 > 0.69, p < 0.05), demonstrating that SOC content over a large coal mining area can be effectively assessed using remote sensing techniques. Results revealed that average SOC content in the study area rose from 5.67 gC · kg − 1 in 1990 to 9.23 gC · kg − 1 in 2010 and then declined to 5.31 gC · Kg − 1 in 2020. This could be attributed to the interaction between the disturbance of soil caused by coal mining and the improvement of eco-environment by land reclamation. Spatially, the SOC content of farmland was the highest, followed by grassland, and that of bare land was the lowest. SOC accumulation was inhibited by coal mining activities, with the effect of high-intensity mining being lower than that of moderate- and low-intensity mining activities. Land use was found to be the strongest individual influencing factor for SOC content changes, while the interaction between vegetation coverage and precipitation exerted the most significant influence on the variability of SOC content. Furthermore, the influence of mining intensity combined with precipitation was 10 times higher than that of mining intensity alone.


Introduction
Coal has been widely acknowledged as an important energy source for humankind, while coal-generated power accounts for over 60% of the total power in China [1]. However, coal mining destroys soil structure, induces changes in soil physicochemical properties, and alters the spatio-temporal distribution of soil organic carbon (SOC) [2,3]. Soil systems, the largest carbon pool in the terrestrial ecosystem, play an essential role in carbon sequestration and emission across the world [4,5]. SOC is an important component of global In this study, the variation and driving mechanisms of SOC in coal mining areas were explored by remote sensing and field sampling in Shenfu Dongshen (Shendong) mining area in China. First, prediction models were established for different land-use types (grassland, forest land, farmland, and bare land) through field sampling and Landsat data. Afterwards, the spatial distribution of the SOC content in 1990-2020 was achieved based on Landsat data and prediction models, followed by an analysis of the spatial pattern of the SOC content and its relation with influencing factors including precipitation, temperature, mining intensity, and vegetation coverage.

Study Area
Shendong mining area (38 • (Figure 1a), is the largest underground coal mining base in China. It has abundant coal reserves, with proved reserves of 2.3 × 10 9 t. The mining area is located in the transitional zone between the Mu Us Desert and the Loess Plateau. It spans horizontally across 3355 km 2 and vertically 945-1419 m above sea level. Gentle topography is commonly found in the western part of the study site, while hilly areas with fragmented and complex landforms dominate the east [61]. The study site is characterized by a typical arid desert plateau climate, with an average annual temperature and precipitation of 6.6 • C and 362 mm, respectively. The rainfall during July-September accounts for more than 70% of the annual precipitation, which is dominated by shortduration intensive rainstorms [62].

Data Acquisition and Processing
In the study, we firstly developed the SOC content prediction models based on field sampling and remote sensing images ( Figure 2). The established models were then used to estimate the SOC content for the whole of Shendong coal mining area in 1990,1995,2000,2005,2010,2015, and 2020. The changes in SOC content and their relations with different influential factors were also examined. Soil types in the study site mainly include loess, aeolian sandy soil, chestnut soil, and dark loessial soil. They are characterized by loose structure, poor erosion resistance, low fertility, and low pH. The study area is also a transitional zone between grasslands and forest-steppes, with sparse native vegetation covering less than 30% of the land surface [63].
Desert steppe, deciduous broad-leaved shrubs, and psammophytes are major vegetation types. The coal seam in the Shendong mining area has a simple structure, which is mainly exploited by the underground mining method. Since the end of the 20th century, the annual coal output of the Shendong area has gradually risen. In 2005, it became the largest coal production base in China by delivering 100 million t. In 2015, the output reached 7.55 × 10 8 t, accounting for 20.13% of the coal production in China [63].

Data Acquisition and Processing
In the study, we firstly developed the SOC content prediction models based on field sampling and remote sensing images ( Figure 2). The established models were then used to estimate the SOC content for the whole of Shendong coal mining area in 1990,1995,2000,2005,2010,2015, and 2020. The changes in SOC content and their relations with different influential factors were also examined. In October 2020, SOC field sampling was carried out. According to traffic, land-use, and by considering the uniformity of geographical distribution, a total of 300 field sampling were collected (Figure 1c), including 75 points in farmland, 54 points in forest land,

Data Collection and Preprocessing
(1) SOC content data In October 2020, SOC field sampling was carried out. According to traffic, landuse, and by considering the uniformity of geographical distribution, a total of 300 field sampling were collected (Figure 1c), including 75 points in farmland, 54 points in forest land, 119 points in grassland, and 52 points in bare land. In order to ensure the accuracy of SOC content prediction in the study area, particularly the edge of the study area, the soil was sampled in some areas outside the study area. The location of each sampling point was determined by the global positioning system (GPS), while sampling time and soil surface condition were also recorded. To ensure soil sampling precision, composite samples (about 1 kg) were collected in ziplock bags, consisting of five soil samples in the surface layer (0-20 cm) within a 30 m × 30 m rectangle area using the diagonal sampling method. After drying in the shade, grinding, and sieving, the SOC content was measured by the potassium dichromate oxidation-external heating method [64].
(2) Remote sensing images In order to facilitate the development of SOC perdition models, remote sensing images with wide spectral range, enough bands and high resolution should be preferred. Meanwhile, a 30-year study period means that it was difficult to cover the whole period with the remote sensing data from a single sensor. Given the above, Landsat5 TM and Landsat8 OLI remote sensing images were selected in our study with the former being used for 1990-2010 and the latter being used for afterwards. The images were obtained from the official web-site of the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/ (accessed on 15 February 2021)). Owing to the fact that (1) the SOC content prediction model should be constructed according to the land use type and (2) the maximum value of Normalized Difference Vegetation Index (NDVI) in each year is required, a total of 21 images from July to September from 1990 to 2020 were selected (orbit No. 127-033, cloud cover <2%). Images were preprocessed using radiometric calibration, atmospheric correction, orthorectification, and cropping using ENVI 5.3.

(3) Land use and vegetation coverage data
Land use data with a 30-m spatial resolution for 1990,1995,2000,2005,2010,2015, and 2020 were obtained from the Data Center for Resources and Environmental Sciences, the Chinese Academy of Sciences (https://www.resdc.cn/ (accessed on 30 June 2021)). The data was derived based on Landsat imagery using the human visual interpretation method. The data have been acknowledged as the most accurate land use data production in China, and have been widely applied in the national land resource survey, hydrology and ecological research [65]. They were classified into six first-level classes according to the actual conditions in the study site.
The maximum NDVI (NDVI max ) was derived based on the Landsat images to represent the spatial pattern of vegetation coverage. The NDVI max of the study area from 1990 to 2020 was calculated for 5-year intervals using Equation (1) based on the Landsat imagery from July to September of each year.
where, NIR is the reflection value in the near-infrared band, and R is the reflection value in the red band.
• Mining intensity data for the Shendong mining area The Shendong coal mining area is dominated by the monoclinal structure with less fault development. The coal mining is dominated by an underground large-scale, highintensity mechanized manner. Coal mining intensity can be determined by underground plane mining intensity and vertical mining intensity. In other words, the mining intensity mainly depends on the length and thickness of the coal seam being taken [66]. According to Fan et al. [66], the mining intensity in the Shendong mining area was classified into five levels: unmined, low-intensity mining, moderate-intensity mining, high-intensity mining, and extremely high-intensity mining (Figure 1b). The high-intensity mining area (28.9%) and unmined area (27.03%) accounted for the largest proportion, followed by extremely high-intensity mining, moderate-intensity mining, and low-intensity mining areas. The extremely high-intensity mining area was located between the Ulan Mulun River and the Kuye River, where 10-million-ton mines such as Daliuta and Huojitu were distributed. The high-intensity mining area was characterized by relatively continuous working faces, large coal mining height (1.3-4.5 m), and high coal output per unit area. In the moderate-and low-intensity mining areas, mechanized mining was not fully achieved. As a result, there were a small number of integrated coal mines with low capacity and relatively tattered working faces.

(4) Meteorological and topographic data
We used the mean annual temperature and annual precipitation as covariates. Data on the temperature, precipitation, latitude, longitude, and altitude of the 24 meteorological stations around the Shendong mining area in 1990,1995,2000,2005,2010,2015, and 2020 were obtained from the China Meteorological Data Network (http://data.cma.cn/ (accessed on 30 June 2021)) ( Figure 1a). Mean annual temperature and annual precipitation were interpolated to continuous surfaces using kriging [67]. Cross-validation analysis was performed to verify the interpolation results, and the accuracy met the criteria of the optimal model [68] (for example, the correlation coefficient between measured precipitation and interpolation results in 2020 was 0.722).
In addition, Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) data with a 30-m resolution were obtained from the Geospatial Data Cloud (https://www.gscloud.cn/ (accessed on 10 October 2020)), and a digital elevation model (DEM) of the study site was obtained through data coordinate transformation and clipping.

Establishment of Models of SOC Content
(1) Band selection and analysis of SOC spectral properties Based on the land-use status and field survey data, the atmospherically corrected images were subjected to analysis for SOC spectral curve to grasp the SOC spectral properties over the research site. In order to expand the difference between spectral properties, improve the correlation between SOC and reflectance data, while eliminating some noise that could not be eliminated in the preprocessing, we combined or mathematically transform the spectral data, including reciprocal (1/B i ), logarithmic (lg B i ), the logarithm of reciprocal (lg 1 Bi ), ratio (B i /B j ), difference (B i -B j )), subtraction-addition ( Bi−Bi+1 Bj+Bj+1 ), and addition-subtraction ( Bi+Bi+1 Bj−Bj+1 ), where B i , B j represent the reflectivity of bands i and j (i, j = 2, 3, 4, 5, 6, 7. i = j). Then, a correlation analysis was conducted between the SOC content of the soil samples and the spectral reflectivity or its mathematical transformation at each sampling point. The characteristic bands that were selected during modeling were extracted from the Landsat8 OLI image at each sampling point. Since both Landsat5 TM and Landsat8 OLI images needed to be used, it was necessary to select the common bands of the two images for modeling. The wavelength ranges of Band1-Band5 and Band7 bands in the Landsat5 TM image are similar to that of the Band2-Band7 bands in the Landsat8 OLI image, respectively. Therefore, only the characteristic bands (Band2-Band7) of Landsat8 OLI remote sensing images and their mathematical transformations were considered for modeling. When the constructed model was applied in periods with the Landsat5 TM image, the selected bands of Landsat 8 images were replaced by the corresponding bands of the Landsat5 TM image. The variables with higher correlation coefficient (R, Equation (2)) were selected as model covariates.
where x i and x are the measured and mean value of the spectral reflectivity (mathematical transformation of reflectivity), y i and y are the measured and mean value of the SOC content.
(2) SOC content prediction models For each land-use type, seventy percent of the sampling points were employed for training the model, while 30% were used for validation. The highly relevant spectral reflectivity (mathematical transformation form of reflectivity) were selected for modeling. The SOC remote sensing prediction models were fitted using multiple linear regression equations. Finally, the model accuracy was assessed using the validation set, with the coefficient of determination R 2 (Equation (3)) and the root mean square error (RMSE) (Equation (4)) being used as performance metrics.
where y i is the measured value, y is the mean value,ŷ i is the predicted value, and n is the number of samples.

Changes of SOC Content and Its Relation with Impacting Factors (1) Changes in SOC content
We used the SOC prediction model of each land use type to spatially map the SOC content in Shendong mining area. Therefore, vectors of mining areas and the SOC content maps were used to calculate the mean SOC content in mining areas with different mining intensities. According to the classification standard of SOC in the second general survey of soil in China [69], the SOC content was classified into 6 levels: I > 23.20 gC·kg −1 , II = 17.40-23.20 gC·kg −1 , III = 11.60-17.40 gC·kg −1 , iv = 5.80-11.60 gC·kg −1 , v = 3.48-5.80 gC·kg −1 , and vi = 0.58-3.48 gC·kg −1 .
Here, the changes in SOC content were studied using the transfer matrix method, which can describe the loss direction and source of SOC of each level within any two adjacent time periods in detail. In addition, the SOC content for the two time periods was mapped, resulting in a SOC content classification transformation map. Afterwards, the change process of SOC content was analyzed. The land-use transfer matrix used is given by the Equation (5): where M ij is the area where the i-th level of SOC content of the previous year was transformed into the j-th level of SOC content of the following year, and n is the SOC content classification number in the study area.
(2) Impacting factors of SOC content changes Geo-detectors are powerful tools for detecting and analyzing spatial heterogeneity of specific variables [70]. The factor detector in a geo-detector is able to detect the explanatory power of the independent variable on the spatial heterogeneity of the dependent variable. The explanatory power is measured by the q value (Equation (6)), which ranges between [0,1]. The higher the value, the stronger the explanatory power of the independent variable on the dependent variable. The interaction detector can identify the interaction between different risk factors, i.e., whether the interaction between factors X 1 and X 2 will increase or decrease the explanatory power on the dependent variable, or whether these factors have an independent influence on Y.
where h = 1, . . . , L is the stratum of the variable Y or factor X, i.e., classification or partition; N h and N are the number of units in the h layer and the whole area; σ 2 h and σ 2 are the variance of the value in the h layer and the whole area, and SSW and SST are the sum of squares within and the sum of squares total, respectively.
The driving factors of the spatial heterogeneity of SOC content were explored using the geo-detector. Based on natural and human factors, six SOC influential factors were selected, including temperature, precipitation, elevation, NDVI, land use, and mining intensity. Since the geo-detector requires the input variables to be categorical, the factors were classified using equal-interval to divide the elevation, annual precipitation, and annual average temperature into nine categories. Moreover, NDVI was divided into eight levels (<0.3, 0.3-0.4, 0.4-0.5, 0.5-0.6, 0.6-0.7, 0.7-0.8, 0.8-0.9 and 0.9-1), and land-use type and mining intensity were classified according to their original categories.

SOC of Soil Samples
The SOC content of the soils at sampling points ranged between 0.66 gC·kg −1 and 31.99 gC·kg −1 , with an average of 6.03 gC·kg −1 and a standard deviation of 5.01 gC·kg −1 . The mean SOC content was at its highest in the farmland, followed by grassland, while it reached the lowest in the bare land class (Table 1).

Established SOC Remote Sensing Prediction Models
Correlation analysis showed that 1/( were highly and significantly correlated with the SOC content (R > 0.5; p < 0.05). The SOC prediction models were established through linear regressions using the highly-correlated bands ( Table 2).   A validation of the prediction models demonstrated that R 2 values for grassland, forest, and farmland were all significant and could reach above 0.73. Furthermore, the RMSEs were lower than 2.74 gC·kg −1 , indicating that the models predict sufficiently accurate results. Both RMSE (1.71) and R 2 (0.69) of the prediction model for bare land were slightly lower than that of other models (Figure 3).    The spatial pattern of SOC content was roughly the same for each study period. That is, it was the lowest in the southeast and gradually rose from southeast to northwest (Figure 5). In 1990, the SOC content in the mining area was lower, mainly at level V. Compared to 1990, an increased SOC content was observed in 1995, and the area with increased content was located primarily in the south of the mining area. In 2000, the SOC content rose significantly, especially in the north. It continuously increased in 2005 and 2010 but did not increase significantly in the northwest, where there were dense regions of mining pits. The spatial pattern of SOC content was roughly the same for each study period. That is, it was the lowest in the southeast and gradually rose from southeast to northwest ( Figure 5). In 1990, the SOC content in the mining area was lower, mainly at level V. Compared to 1990, an increased SOC content was observed in 1995, and the area with increased content was located primarily in the south of the mining area. In 2000, the SOC content rose significantly, especially in the north. It continuously increased in 2005 and 2010 but did not increase significantly in the northwest, where there were dense regions of mining pits. In 2010, the SOC content in the mining area reached its historical peak, after which the local SOC content in the mining area declined in 2015. In 2020, the SOC content in the whole area significantly decreased compared with that in 2015, mainly in the dense region of pits. Compared with that in 1990, the increase in SOC content was not apparent in the northwest, but rather in the south. Overall, the local SOC content in the Shendong mining area was gradually increasing but did not increase in the dense region of pits. The results of the transfer matrix analysis (Figure 6, Tables 3 and 4) revealed that the SOC content in the study area was initially low. In 1990, the SOC content was mainly at level VI, accounting for 67.6% of the site. In 2010 and 2020, SOC content was mainly at level IV, accounting for 45.28% and 33.84% of the area, respectively. More specifically, from 1990 to 2010, the overall trend in SOC content variation was from low to high. In 2010, the area with a SOC content of level VI accounted for 89.95% of the site, while the area with a level IV SOC content increased to 1097.04 km 2 . From 1990 to 2010, the area of SOC content change was 2652.4 km 2 , accounting for 88.15% of the total soil area. The area with the increased SOC content was 2092.1 km 2 , accounting for 70.8% of the total soil area. However, the SOC content continuously declined due to the increase in construction land and the massive exploitation of coal resources from 2010 to 2020. The area where the SOC content decreased was 78.09% of the total changed area and 50.9% of the total soil area. This decrease included a 214.69 km 2 area where the SOC content dropped to 0 gC·kg −1 , accounting for 14.27% of the area with decreased SOC content. The results of the transfer matrix analysis (Figure 6, Tables 3 and 4) revealed that the SOC content in the study area was initially low. In 1990, the SOC content was mainly at level VI, accounting for 67.6% of the site. In 2010 and 2020, SOC content was mainly at level IV, accounting for 45.28% and 33.84% of the area, respectively. More specifically, from 1990 to 2010, the overall trend in SOC content variation was from low to high. In 2010, the area with a SOC content of level VI accounted for 89.95% of the site, while the area with a level IV SOC content increased to 1097.04 km 2 . From 1990 to 2010, the area of SOC content change was 2652.4 km 2 , accounting for 88.15% of the total soil area. The area with the increased SOC content was 2092.1 km 2 , accounting for 70.8% of the total soil area. However, the SOC content continuously declined due to the increase in construction land and the massive exploitation of coal resources from 2010 to 2020. The area where the SOC content decreased was 78.09% of the total changed area and 50.9% of the total soil area. This decrease included a 214.69 km 2 area where the SOC content dropped to 0 gC·kg −1 , accounting for 14.27% of the area with decreased SOC content.

SOC Distribution under Different Mining Intensities
The mean SOC content was the highest in the high-intensity mining area. This was followed by the areas of extremely high-intensity mining, unmined, low-intensity mining, and moderate-intensity mining. In addition, the mean SOC content was on the rise in 1990-2010 in each mining area, except in moderate-intensity and extremely high-intensity mining areas. It declined in 2000-2005 in the moderate-intensity mining area and in 2005-2010 in the extremely high-intensity mining area. In 1990-2010, the mean SOC content in the mining area continuously rose, with significant differences in the SOC content among mining areas that were observed in 2005. The mean SOC content in the high-intensity and extremely high-intensity mining areas was significantly higher than in the mining area and other areas. Meanwhile, in other areas it was lower than the average in the mining area.
Furthermore, the mean SOC content in the moderate-intensity mining area displayed a decreasing trend, and reached a peak in the extremely high-intensity mining area. In 2010, the mean SOC content in the extremely high-intensity mining area declined (Figure 7). However, it was still at the average level of the mining area. The mean SOC content was still significantly higher in the high-intensity mining area than in the mining area and other areas. In 2010-2020, the mean SOC content in each mining area declined. The mean SOC content in unmined areas slowly increased from 1990-2005, but it remained unchanged in 2005-2015, whereas the mean SOC content in other areas experienced significant changes. The mean SOC content was the highest in the high-intensity mining area. This was followed by the areas of extremely high-intensity mining, unmined, low-intensity mining, and moderate-intensity mining. In addition, the mean SOC content was on the rise in 1990-2010 in each mining area, except in moderate-intensity and extremely high-intensity mining areas. It declined in 2000-2005 in the moderate-intensity mining area and in 2005-2010 in the extremely high-intensity mining area. In 1990-2010, the mean SOC content in the mining area continuously rose, with significant differences in the SOC content among mining areas that were observed in 2005. The mean SOC content in the high-intensity and extremely high-intensity mining areas was significantly higher than in the mining area and other areas. Meanwhile, in other areas it was lower than the average in the mining area.
Furthermore, the mean SOC content in the moderate-intensity mining area displayed a decreasing trend, and reached a peak in the extremely high-intensity mining area. In 2010, the mean SOC content in the extremely high-intensity mining area declined ( Figure  7). However, it was still at the average level of the mining area. The mean SOC content was still significantly higher in the high-intensity mining area than in the mining area and other areas. In 2010-2020, the mean SOC content in each mining area declined. The mean SOC content in unmined areas slowly increased from 1990-2005, but it remained unchanged in 2005-2015, whereas the mean SOC content in other areas experienced significant changes.

Influential Factors of SOC Changes
The geo-detector revealed that the explanatory power of land-use on the SOC content was the strongest (q = 0.165), followed by NDVI (q = 0.149), precipitation (q = 0.138), temperature (q = 0.115), elevation (q = 0.111), and mining intensity (q = 0.041), where all influential factors passed the significance test of 99% CI.
The statistical results of interaction detector (Table 5) showed that the interaction between the influential factors was stronger than that of a single factor. The interaction of NDVI with precipitation and elevation ranked first, with its interaction with precipitation being the mostly obvious (q = 0.555). The interaction of mining intensity with other factors

Influential Factors of SOC Changes
The geo-detector revealed that the explanatory power of land-use on the SOC content was the strongest (q = 0.165), followed by NDVI (q = 0.149), precipitation (q = 0.138), temperature (q = 0.115), elevation (q = 0.111), and mining intensity (q = 0.041), where all influential factors passed the significance test of 99% CI.
The statistical results of interaction detector (Table 5) showed that the interaction between the influential factors was stronger than that of a single factor. The interaction of NDVI with precipitation and elevation ranked first, with its interaction with precipitation being the mostly obvious (q = 0.555). The interaction of mining intensity with other factors was significantly improved, with the q value being 8 times that of a single factor. In particular, the q value of the mining intensity-precipitation interaction reached 0.408, which was lower than that of the NDVI-precipitation interaction and NDVI-elevation interaction, but 10 times that of a single factor.

SOC Prediction Models of Coal Mining Area
So far, remote sensing techniques have been widely applied in the dynamic monitoring of physicochemical properties of regional soils [71][72][73], but there have been few studies on soil properties in coal mining areas, specifically at a large scale and within complex landscapes. Due to long-term mining of different intensities in coal mining areas, varying degrees of land surface subsidence have occurred, and land use has significantly changed through time. It resulted in the complex topography and land use distribution. It was expected that the remote sensing spectral information would be affected to a certain extent [17,60]. Different from the SOC prediction model that was solely constructed using image information, we incorporated some additional data during the development of our prediction models. Compared to other SOC remote sensing prediction models, the established models accounted for the impact of land-use types on SOC sensitive bands, and thereby improved the model precision. Several studies have confirmed that the use of ancillary data can improve the accuracy of soil prediction models [74,75].
The number and heterogeneity of the sample points are amongst the key factors that affect the performance of prediction models [76]. Here, we collected a large number of soil samples across sampling sites that covered different landscapes properties such as topography, altitude, soil types, and land use types. Although the constructed model takes the interference of land use types on spectral information into account, it does not fully consider the influence of soil types and soil properties on spectral information. This might be a reason for the low R 2 of the prediction model across bare land type. The soil parent material and land-use type have a greater impact on the visible-near infrared spectrum of the soil [77]. Besides, we found that pulverized coal is also an important factor that affected the accuracy of SOC prediction in coal mining areas. The development and utilization of coal resources generates a large amount of dust, which makes the coal mine area seriously polluted by pulverized coal. As a result, the spectral reflectance of pulverized coal has a certain spectral response corresponding to the 4 th and 5 th bands of the TM imagery, which makes the original spectral curve subject to a certain degree of interference [78]. Therefore, soil types and soil attributes should be used as auxiliary data for remote sensing prediction modeling in future research. Besides, the interference of pulverized coal on spectral information should be fully considered.

Spatiotemporal Change Properties of SOC
The mean SOC content in the Shendong mining area slowly rose in 1990-2010, while it rapidly declined in 2010-2020, showing an overall decreasing trend at a rate of 0.018 gC·kg −1 /a. In 1990, the study site had not yet been exploited, and the natural ecosystem was less disturbed by humans. Still, aeolian sandy soil was dominant due to the poor soil, thus the SOC content was not high. However, with coal mining and land reclamation, the SOC content continuously increased during 1990-2010, which was closely related to the mining area's policy that gave equal importance to coal mining and the ecological management of the environment [63]. One may note that the area's vegetation cover was still extremely sparse though it had not been exploited in 1990. During the mining period, vegetation restoration was adopted, which improved the vegetation coverage rate from 3-11% to more than 60% in the initial stage of development. The increase in vegetation is especially important for the accumulation of SOC.
The SOC content was spatially consistent in each period. That is, it gradually rose from southeast (with the lowest content) to northwest. A possible reason is that the land is flat in the northwest and locally covered with grassland and forest land, which can inhibit soil erosion and increase carbon sink to a certain extent. In contrast, the gully area is located in the southeast, which is more prone to soil erosion than the northwest in the mining area. This resulted in low SOC content [63].
The above analysis on the changes in SOC content under different mining intensities demonstrated that the different mining intensities impacted the changes in the SOC content across the study area. The SOC content was lower in low-and moderate-intensity mining areas, mainly because of mostly occurring small and medium-sized coal mines or small collieries with low productivity, great long-term destruction, and no advanced mining technology. Although technological advancements in coal integration and industrial structure in recent years have, to a certain extent, curbed the development of high-energy-consuming enterprises, the recovery from the local ecological impact caused by many years of industrial coal activity is a long-term process. In contrast, the SOC content in high-and extremely high-intensity mining areas was higher than in other mining areas. The main reason is that the vegetation in these areas has not been regularly damaged due to mining-induced geological disasters, such as collapse and groundwater subsidence. Mining with a superlarge working face has also been confirmed in many studies to cause less damage to the environment, which reduces the range of uneven ground subsidence and stretching. It also causes spontaneous restoration of the surface vegetation to a large extent [62,63]. The mean SOC content was also higher in these areas. The catchment area is formed due to surface subsidence, and SOC will migrate to this area with soil erosion, which leads to higher SOC content [79,80]. This study did not consider the impact of topographical changes caused by coal mining on SOC content. We only classified the mining intensity in our large-scale study. However, the impact of coal mining on the redistribution of SOC remains unclear, thus further research on topographical changes is required in future studies.

Influential Factors of SOC Changes
SOC content is mainly affected by climate, topography, vegetation, land use and soil erosion [81]. In the coal mining area, SOC content is also affected by coal mining activities. Here, various influential factors were analyzed using geo-detectors.
The land-use type possessed the highest explanatory power for the spatial heterogeneity of SOC content, since there is a difference in carbon sequestration capacity of different land-use types. Furthermore, water supply and soil water retention vary due to the differences in evapotranspiration and sediment retention [25,82,83]. In the Shendong mining area, the cultivated land showed the highest average SOC content, followed by grassland, while the bare land showed the lowest value. This was different from previous studies [27,81,84]. We argue this to be related to people fertilizing the originally poor soil in order to increase the yield of crops. The soil itself in Shendong mining area is infertile, the native vegetation is sparse, and the coverage is less than 30% [63]. Therefore, SOC is difficult to reach high values under natural conditions. Although the cultivated land has been greatly disturbed, the SOC content is much higher than other land use types, which was apparently due to fertilization. The reason why the average SOC content of grassland was higher than that of forest may be due to the occurence of immature forests in the study area. Studies have shown that the soil organic carbon content of immature forests(10 years old) can be for example 32.25% lower than that of 30-year-old forest land [81].
Further, carbon sink capacity is relatively high in regions with high NDVI. These regions that facilitate the increase in SOC content through resisting rainfall erosion and other physical forms of erosion [85]. The surface subsidence caused by mining disturbs the vegetation root system and reduces the carbon sequestration capacity, which intensifies soil erosion and SOC mineralization [86][87][88]. In addition, precipitation and temperature have an important impact on the photosynthesis of plants, thereby affecting the accumulation of SOC. In the meantime, water and temperature have a huge impact on the activities of soil microorganisms, which control the rate of transformation and decomposition of plant residues that enter the soil by microorganisms, and thereby change both input and output of soil carbon [89]. Therefore, the explanatory power of temperature and precipitation for the spatial heterogeneity of SOC content was lower than that of NDVI. Furthermore, the difference in altitude affects the annual precipitation and temperature. In the meantime, high-altitude areas receive less human disturbance and are easier to accumulate SOC. In our study the effect of altitude on SOC was less than that of annual precipitation, which was in line with previous studies [90,91].
Mining intensity as an individual influencing factor had no apparent effect on the change of SOC content in our study site. Therefore, the focus should be on the optimization of land-use structures that aim to improve vegetation coverage. Further recommendations include that the direction and speed of vegetation succession should be artificially controlled. The forest structure in such areas should also be optimized by combining grasses as primary vegetation with shrubs as secondary vegetation type. In addition, scientific methods such as microbial reclamation should be used to reclaim the land subsidence in mining areas, while the greening of bare land should be enhanced.
The interaction between various factors that influence SOC content was higher than that of a single factor. The interaction between mining intensity and precipitation displayed the most significant effect on SOC content, with an explanatory power that was ×10 of a single factor. As mentioned above, mining intensity had no pronounced effect on the change in SOC content, which was also consistent with the result from the geo-detector. Moreover, the explanatory power of mining intensity on the change of SOC content was relatively low. It was found that the ground deformation caused by various mining intensities was different, which may be the main influential factor on SOC content.
Coal mining disturbs vegetation and soil, which may lead to a decrease in soil water retention capacity. This then intensifies erosion when precipitation increases. The SOC will then migrate [79,80], leading to a relatively high SOC content in the bottom of the collapsed basin caused by coal mining. This can also explain the high explanatory power of the interaction between mining intensity and precipitation, as observed in this study.

Implications
Our study has both scientific and practical implications. Firstly, changes in SOC and their relation with influential factors are globally important, given its close connection with global climate change. The coal mining area shows a typical example for areas with intensive human disturbance, which facilitated a comparison of the impacts of climate change and human activities. We found that the interaction between vegetation and precipitation dominated the SOC changes in the Shendong mining area, although intensive coal mining activities exerted negative effects on SOC accumulation. This implied that vegetation restoration is useful for the control of SOC loss induced by coal mining activities. This is particularly important given that the Loess Plateau has been projected to experience a warmer and wetter climate in the future [92]. Our findings are beneficial for an optimization of coal mining activities in the context of global climate change as well as land use planning to achieve a low SOC loss in areas with intensive human disturbances. Secondly, the model development in our study provides a useful reference for the establishment of SOC prediction models in other regions, particularly for areas with intensive human activities. The established models also provide a powerful tool for an effective assessment of SOC over a large area, which is then beneficial for an investigation of driving mechanisms of SOC changes, and thus the its conservation.

Conclusions
In order to explore changes of SOC content in Shendong mining area, this study established prediction models through an analysis of field measured SOC content and spectral information of Landsat images for different land use types including grassland, forest, farmland and bare land. The established models were validated to be reasonable with the R2 being higher than 0.69 and RMSE being no more than 2.74.
The SOC content of the Shendong mining area was found to increase from 5.67 gC·kg −1 in 1990 to 9.23 gC·kg −1 in 2010, and reduced to 7.92 gC·kg −1 in 2015 and 5.31 gC·kg −1 in 2020, which showed a declining trend at a rate of 0.018 gC·kg −1 /a. Spatially, the SOC content reached a minimum in the southeast of the Shendong mining area and generally rose gradually from southeast to northwest.
Coal mining affects the spatial distribution of SOC in the Shendong mining area to a certain extent, the high-intensity mining area had the highest SOC content, followed by the extremely high-intensity mining area, unmined area, low-intensity mining area, and moderate-intensity mining area. In the high-intensity and extremely high-intensity mining areas, the inhibition of SOC accumulation reached the lowest level. In contrast, the inhibition peaked in the low-and moderate-intensity mining areas.
Although the coal mining intensity has an influence on the SOC content, SOC content changes were dominated by other factors. The geo-detector analysis demonstrated that the most influential factor on the SOC content was land use, followed by NDVI, precipitation, temperature, elevation, and mining intensity. The interaction between the various influential factors on SOC content was higher than that of a single factor, while the interaction between NDVI and precipitation had the strongest influence on SOC content. The influence of the interaction between mining intensity and other factors increased. The interaction of coal mining and precipitation has a greater impact on SOC content than coal mining has on SOC content. Our study provided an enhanced understanding of SOC content changes in mining areas and the results are beneficial for the conservation of terrestrial carbon storage, particularly for areas with intensive human activities. Acknowledgments: Special thanks go to Fan for the mining intensity data of the Shendong mining area. United States Geological Survey, Data Centre for Resources and Environmental Sciences of Chinese Academy of Sciences, China Meteorological Data Network and Geospatial Data Cloud are acknowledged for providing remote sensing images, land use datasets, meteorological data and DEM data.

Conflicts of Interest:
The authors declare no conflict of interest.