Changes in the Impacts of Topographic Factors, Soil Texture, and Cropping Systems on Topsoil Chemical Properties in the Mountainous Areas of the Subtropical Monsoon Region from 2007 to 2017: A Case Study in Hefeng, China

Understanding the spatial pattern of soil chemical properties (SCPs) together with topological factors and soil management practices is essential for land management. This study examines the spatial changes in soil chemical properties and their impact on China’s subtropical mountainous areas. In 2007 and 2017, 290 and 200 soil samples, respectively, were collected in Hefeng County, a mountainous county in central China. We used descriptive statistics and geostatistical methods, including ANOVA, semivariance, Moran’s I, and fractal dimensions, to analyze the characteristics and spatial autocorrelation changes in soil organic matter (OM), available phosphorus (AP), available potassium (AK), and pH value from 2007 to 2017. We explored the relationship between each SCP and the relationship between SCPs with topographic parameters, soil texture, and cropping systems. The results show that the mean value of soil OM, AP, AK, and pH in Hefeng increased from 2007 to 2017. The spatial variation and spatial dependency of each SCP in 2007, excluding AP and AK in 2007, were higher than in 2017. The soil in areas with high topographic relief, profile curvature, and planform curvature had less AP, AK, and pH. Soil at higher elevation had lower OM (r = −0.197, p < 0.01; r = −0.334, p < 0.01) and AP (r = −0.043, p < 0.05; r = −0.121, p < 0.05) and higher AK (r = −0.305, p < 0.01; r = 0.408, p < 0.01) in 2007 and 2017. Soil OM and AK in 2007 were significantly (p < 0.05) correlated with soil texture (p < 0.05). In contrast, oil AP and soil pH in 2007 and all SCPs in 2017 were poorly correlated with soil texture. The cropping systems played an important role in affecting all SCPs in 2007 (p < 0.01), while they only significantly affected AK in 2017 (p < 0.05). Our findings demonstrate that both topological factors, that is, the changes in cropping management and the changes in acid rain, impact soil chemical properties. The local government should place more focus on reducing soil acid amounts, soil AP content, and soil erosion by improving water conservancy facilities.


Introduction
Soil is the foundation of life and biodiversity [1,2]. Climate change and human intervention speed up soil property change [3]. Soil degradation is a global issue that leads to crop reduction. Understanding soil properties and properly managing soil are critical for avoiding soil degradation [4,5] and achieving the United Nations' sustainable development goals [6].
In the past, soil degradation in China was serious due to intensive and unsustainable agricultural activities. Soil problems such as soil acidity in South China [7] and soil cadmium in central China [8] threaten national food security. With a population of 1.4 billion people and a rising demand for improved food quality and food structure, the Chinese The study area has a temperate continental monsoon climate. The annual precipitation (1701 mm) is distributed mainly from June to September [19]. According to the Second State Soil Survey of China (SSSSC, 1982), there are 10 main soil types in Hefeng County ( Figure 2): rice soil, tidal soil, limestone soil, purple soil, red soil, yellow soil, yellow-brown soil, meadow soil, swamp soil, and brown soil. The most widely distributed soil in this area is the mountainous yellow-brown soil, which is weakly acidic but high in soil fertility.

Sampling Design and Laboratory Analyses
Three large-scale soil sampling activities have been carried out in Hefeng County, namely, the SSSSC, the soil for formulated fertilization test (2007), and the soil series survey (2017). Because these projects were carried with different standards, we only used the The study area has a temperate continental monsoon climate. The annual precipitation (1701 mm) is distributed mainly from June to September [19]. According to the Second State Soil Survey of China (SSSSC, 1982), there are 10 main soil types in Hefeng County ( Figure 2): rice soil, tidal soil, limestone soil, purple soil, red soil, yellow soil, yellow-brown soil, meadow soil, swamp soil, and brown soil. The most widely distributed soil in this area is the mountainous yellow-brown soil, which is weakly acidic but high in soil fertility. The study area has a temperate continental monsoon climate. The annual precipitation (1701 mm) is distributed mainly from June to September [19]. According to the Second State Soil Survey of China (SSSSC, 1982), there are 10 main soil types in Hefeng County ( Figure 2): rice soil, tidal soil, limestone soil, purple soil, red soil, yellow soil, yellow-brown soil, meadow soil, swamp soil, and brown soil. The most widely distributed soil in this area is the mountainous yellow-brown soil, which is weakly acidic but high in soil fertility.

Sampling Design and Laboratory Analyses
Three large-scale soil sampling activities have been carried out in Hefeng County, namely, the SSSSC, the soil for formulated fertilization test (2007), and the soil series survey (2017). Because these projects were carried with different standards, we only used the

Sampling Design and Laboratory Analyses
Three large-scale soil sampling activities have been carried out in Hefeng County, namely, the SSSSC, the soil for formulated fertilization test (2007), and the soil series survey (2017). Because these projects were carried with different standards, we only used the official report of the SSSSC in Hefeng and the soil chemical data with a location in 2007 and 2017. Moreover, the common indicators in the database of chemical properties of the soil sample sites for these two years were soil organic matter content, soil available phosphorus content, soil available potassium content, and soil pH. Therefore, we only studied the spatial distribution of these soil chemical properties in 2007 and 2017 and took the SSSSC report in Hefeng as a reference.
Soil sample locations were selected based on spatial homogeneity, soil type, and crop type. Sampling locations were determined using a Garmin GPS Etrex10 receiver. Samples were taken after harvest in November 2007 (290 samples) and November 2017 (200 samples) when no fertilizer or pesticides were applied in the cropland, which would ensure that environmental noise and systematic errors in the soil samples were eliminated as much as possible ( Figure 2). Interviews with landowners were conducted at each sample site, and information was collected on their agriculture management practices such as crop rotation, yield, and irrigation.
At each location, a soil drill was used to collect five topsoil samples (0-20 cm below the surface) according to the "S" method. Five topsoil samples were mixed to obtain composite soil samples. Each soil sample (>1 kg) was placed in an aluminum box, which was then sealed. Soil organic matter (OM) was tested by dichromate wet combustion. Available phosphorus (AP) was determined by extracting samples with 0.5 mol L −1 sodium bicarbonate (NaHCO 3 ) and subsequent colorimetric analysis [21]. Available potassium (AK) was extracted with 1 mol/L NH 4 Ac, and then measured by an atomic absorption spectrometer [22]. Soil pH was determined by a pH meter (Sartorius Basic pH meter PB-10, Gottingen, Germany) with a soil/water ratio of 1:2.5 [23]. These same methods were used for both 2007 and 2017.

Descriptive Statistics and Difference Tests
The soil sample collection and agrochemical analysis are affected by uncertain factors, resulting in gross errors in soil's measurement. Gross errors are often manifested as high or low values that clearly deviate from the data, thereby affecting the overall distribution characteristics and statistical analysis of the sample. We used the 3σ rule for outlier detection [24]: where X i,j is the value of sample i and soil property j; E X j is the average value of the whole sample of the soil property j; and σ is the standard deviation of the whole sample of the soil property j. If X i,j is lower than E X j − 3σ or higher than E X j + 3σ, then it is deleted when analyzing the data. The sample numbers after outlier detection were 271, 263, 265, and 271 in 2007, and 192, 191, 183, and 184 in 2017, for OM, AP, AK, and pH, respectively. We used SPSS 25.0 for soil sample statistics, including mean, maximum, minimum, standard deviation, coefficient of variation, bias, kurtosis, and significance tests (Kolmogorov-Smirnov test, p < 0.05). Before spatial interpolation, we needed to assure that the soil variables were continuous and normally distributed. If the variables are not normally distributed, the semivariance function will be distorted, which reduces the precision of the interpolation results. Before spatial interpolation, we used a histogram, a normal Q-Q plot, and Kolmogorov-Smirnov tests to check whether the data were normally distributed. When the test results did not conform to the normal distribution, we performed either a logarithmic transformation, square root transformation, or Box-Cox transformation on the data. When the data were close to the normal distribution, we performed spatial interpolation.
Soil texture is influenced by topography, soil texture, and human activities. Five topographic factors were selected in this study, namely, elevation, slope, topographic relief, profile curvature, and planform curvature, to analyze the correlation between topography and soil texture by Pearson's correlation [25] in SPSS 25. We used a 100 × 100 m moving window in the focal statistics toolkit in ArcGIS 10.7 to generate the raster for maximum and minimum elevation within a 100 × 100 m moving window. The topographic relief of cell i was the maximum elevation in cell i minus the minimum elevation value in cell i.
The profile curvature affects the acceleration and deceleration of flow and; therefore, influences erosion and deposition. The planform curvature influences convergence and divergence of flow. A negative value indicates that flow will be decelerated on the surface, while a positive profile indicates that the flow will be accelerated. A value of zero indicates that the surface is linear. The planform curvature is perpendicular to the direction of the maximum slope. The planform curvature relates to the convergence and divergence of flow across a surface. A positive value indicates that the surface is laterally convex at that cell. A negative plan indicates that the surface is laterally concave at that cell. A value of zero indicates that the surface is linear [26][27][28].
The topographic factors were extracted from the digital elevation model (Shuttle Radar Topography Mission, 30 m [29]) by the Spatial Analysis Extension in ArcGIS 10.7; four main soil textures of Hefeng County and eight cropping systems were used to analyze their relationship with soil chemical properties by one-way ANOVA using SPSS 25 [30,31]. The selected soil texture includes light loam, medium loam, sandy loam, and sandy soil. The cropping system is indicated in Figure 2.

Geostatistical Analysis
(1) Semivariance We used GS + 9.0 (Gamma Design Software, LLC, Plainwell, MI, USA) [32] to calculate each indicator's semivariance and used ArcGIS 10.7 (Esri's GIS mapping software, Esri, Redlands, CA, USA) to map the soil chemical properties. The semivariance function is used in geostatistics to quantitatively describe the spatial variation structure. The equation of the semivariance function is as follows: where γ(h) represents the semivariance with the lag distance (h), N(h) is the sample pair's number within the lag distance h, z(x i ) is the value of sample i, and z(x i + h) is the value of the sample at distance h to sample i. We used the best-fit semivariogram models for each soil property, such as a stable model, circular model, Gaussian model, exponential model, or spherical model. Values were calculated for each possible pair of observation points, and the mean values of the semivariance were displayed for chosen distance intervals to produce the experimental semivariogram [33]. The sill (C0 + C) is the constant value when the semivariance function γ(h) increases with the lag distance h and reaches a relatively stable constant from the nonzero value. The nugget is the value when the lag distance h = 0, γ (0) = C0. The largest variation in the system and the lag distance α when the semivariance function γ(h) reaches the base station value are called the base station value and range, respectively. C0 represents the spatial heterogeneity of the random part, and C represents the spatial heterogeneity of system variation; the nugget-to-sill ratio (C0/C0 + C) was used to define different spatial class dependencies for the soil properties. C0/C0 + C was used to define different classes of spatial dependence for the soil properties. C0/C0 + C < 0.25, 25-0.75, and >0.75 were classified as having strong, moderate, and weak spatial dependence, respectively [11]. The stronger the spatial dependence in one soil property, the lower the C0/C0 + C value and the less random factor interruption (such as changes in cropping management and fertilization) experienced by the soil property [11].
(2) Spatial Autocorrelation Analysis is a method of examining the autocorrelation between a regional variable and its neighboring variables. It determines whether there is spatial autocorrelation by detecting the dependence of a position variation on its neighboring positions [34]. We used global Moran's I index and the Z scores to detect the spatial correlation between the variables. The equation of Moran's I index is as follows: where I is global Moran's I index, n is the number of samples, X i and X j are the samples' property value at locations i and j, X is the mean of samples, and W i,j is the weight based on the distance between samples i and j. The value of Moran's I varies from −1 (negative autocorrelated) to 1 (positive autocorrelated). When the value of I equals 0, there is no autocorrelation between sample i and sample j. Z scores (standard scores) can test whether the spatial autocorrelation is significant [35]. The formula is as follows: where x is the value of sample, µ is the mean value of the population, and σ is the standard deviation of the population. The Z score can determine the degree of spatial autocorrelation, and the larger the absolute value, the greater the spatial correlation; on the contrary, the smaller the absolute value, the smaller the spatial correlation. When the Z score is close to 0, it indicates that the space of the regional variable is randomly distributed. (

3) Fractal Dimensions
We used fractal dimensions (FD) to describe the spatial variability of soils. Fractal dimension is a measure of self-similar complexity [36]. The higher the FD value, the more complicated the landscape. It is used to describe the spatial correlation characteristics of soil chemical properties. Burrough introduced this fractal theory into soil science in 1983 and identified that different soil properties have different fractal characteristics [37]. The soil fractal dimension (D) is calculated as follows: where FD means fractal dimension, FD ∈ [1,2], and H is the slope of the linear regression of log γ(h) ∝ log h in the scale h,H ∈ (0,1). FD = 2 indicates that there is no spatial autocorrelation in the distribution of the soil property, while FD = 1 indicates that there is a strong spatial auto-correlation in the distribution of the soil property [38].
In this study, we combined semivariance, spatial autocorrelation, and the fractal dimension to find the spatial variability of soil properties from different perspectives. There are differences in the calculation methods of semicovariance function and spatial autocorrelation. By calculating variance, the semivariance function quantifies the degree of spatial autocorrelation and the scope of the spatial variability scale of regional variables, and it provides the basis for kriging interpolation. The disadvantage of the semivariance function, however, is that it cannot provide statistical tests for significant and positive or negative spatial correlation [33]. This problem can be solved by the spatial autocorrelation index, which calculates covariance and makes up for the lack of significance tests in the semivariance model. Since the fractal dimension is a dimensionless composite indicator [38], it can be used to reveal spatial heterogeneity and provide empirical evidence for spatial autocorrelation and variance functions. Each of these three analytical tools has its own advantages and disadvantages, and they are used in combination to complement each other. The three tools are combined to comprehensively analyze the spatial structural characteristics and variation patterns of soil chemical properties.

(4) Spatial Interpolation
Spatial interpolation is the main way to identify the changes in the spatial distribution pattern of soil nutrient content. This research uses kriging interpolation [39,40], the most commonly used in soil science [41], to explore the spatial distribution of soil chemical properties in Hefeng. The method is based on semivariogram theory and structural analysis to predict unsampled points. We used GS+9.0 to fit the optimal half-variance theoretical model combined with ordinary kriging interpolation in ArcGIS 10.7 for spatial interpolation. We verified the interpolation accuracy of this study, as shown in Table 1.
The results indicate that the interpolation accuracy is high. The root mean square error (RMSE) of each property was higher than one in the study period, which means that the interpolated value is positively related to the real value. The standard error of mean (SEM) of every property is close to 0. The value of the root mean square error (RMSE) of each indicator is similar to the value of the mean error (MAE) [42,43], and the standard root mean square error (RMSS) of different soil properties is close to one. The values of RMSE, SEM, RMSS, and MAE all prove that the interpolation has reasonable accuracy ( Table 1). We classified each soil property in the maps according to the standard of the SSSSC [44,45] ( Table 2), which is the official standard of soil property classification in China. It has been often referred to by other researchers [46][47][48][49][50] and is the geochemical survey standard used in China [51].  Note: OM-organic matter, AP-available phosphorus, AK-available potassium The unit of OM is g/kg; the unit of AP and AK is mg/kg.

Descriptive Statistics
In the K-S (Kolmogorov-Smirnov) test, only when the p-value is higher than 0.05 are the data normally distributed. We transformed all original soil data with a natural logarithm transformation, because in the K-S test, the p-value of all of the chemical soil properties was lower than 0.01.
The mean value of soil OM, AP and pH in Hefeng increased from 33.64 g/kg, 39.56 mg/kg, and 5.36 in 2007 to 35.49 g/kg, 58.51 mg/kg, and 5.99 in 2017, respectively ( Table 3). The mean value of soil AK decreased from 147.82 mg/kg in 2007 to 118.02 mg/kg in 2017. Moreover, the maximum value of all indexes except pH decreased. The minimum value of all indexes increased except AK. The standard deviation (SD) and coefficient of variation (CV) value of all of the soil properties declined in the period of 2007-2017, which means that the spatial variation of each soil property decreased. The CV values of all indexes in both years are higher than 10%, which means that each index is moderately variable [52]. The AP in 2007 had the highest variability among all of the indicators, with a CV value of 76.46%. In 2007, soil pH is poorly correlated with other soil chemical properties. Nevertheless, soil OM is significantly positively correlated with soil AP (r = 0.192, p < 0.01) and soil AK (r = 0.245, p < 0.01). Soil AP is significantly positively correlated with soil AK (r = 0.241, p < 0.01). In 2017, AP and AK no longer correlate well with OM, and AP is poorly correlated with AK. The correlations between other indices are significant but generally weak (p < 0.05) ( Table 4).

Pearson's Correlation between Topographic Factors and Soil Chemical Properties in Hefeng
According to Pearson's correlation analyses, in 2007 and 2017, OM (r = −0.197, p < 0.01) and AP (r = −0.043, p < 0.01) had a negative relationship with elevation. However, the other topographic factors had limited effects on OM and, with both negatively related to the slope (r = −0.159, p < 0.05) and planform curvature of land (r = −0.219, p < 0.05) in 2007 and AP negatively related to the slope in 2017 (r = −0.154, p < 0.01). AK positively was related to elevation in 2007 (r = 0.305, p < 0.01) and 2017 (r = 0.408, p < 0.01) and negatively related to topographic relief (r = −0.173, p < 0.0) and profile curvature (r = −0.175, p < 0.05) in 2017. The pH was negatively related to planform curvature of land in 2007 (r = −0.227, p < 0.01) and 2017 (r = −0.128, p < 0.01). In conclusion, the soil at a higher altitude was lower in OM and AP but high in AK. The soil located on steeper land had a lower AP and lower AK, while soil on land with high topographic relief, profile curvature, and planform curvature had less AP, AK, and pH (Table 5).

Soil Chemical Properties in Different Soil Texture
The soil chemical properties in the various soil textures were different (Table 6). Soil OM and AK in 2007 are significantly (p < 0.05) correlated with soil texture (p < 0.05). On the other hand, soil AP and soil pH in 2007 and all soil chemical properties in 2017 are poorly correlated with soil texture. The OM content was highest in medium loam and lowest in sandy soil. The soil AP contents were highest in medium loam and lowest in sandy soil. The correlation between different soil textures and soil fertility content decreased in the period of 2007-2017 (Table 6).

Soil Chemical Properties in the Different Cropping Systems
With Figure 3, we were able infer that farmers plant the crop according to the soil chemical properties, and that farming activities also influence the variation of soil properties. The changes in soil chemical properties vary from crop to crop. The cropping systems were well correlated with soil chemical properties in 2007 (p < 0.01). On the other hand, we only observed a significant correlation between AK and cropping systems in 2017 (p < 0.05). Based on cropping system types, the highest mean OM values were recorded for the soil where vegetables (43.79 g/kg) and rice (37.

Semivariogram and Spatial Autocorrelation
The most appropriate model was selected based on the maximum coefficient of determination (R2). R2 in this study ranges from 0.679 to 0.875, and the whole residual sum of squares (RSS) is lower than 0.001, which indicates that the semivariance models of each soil property are fitted (Table 8). With the exclusions of log (OM) in 2007, which was best fitted with the spherical model, and log (pH) in 2007 and 2017, which was best fitted with the Gaussian model, the soil chemical properties were best fitted with the exponential

Semivariogram and Spatial Autocorrelation
The most appropriate model was selected based on the maximum coefficient of determination (R2). R2 in this study ranges from 0.679 to 0.875, and the whole residual sum of squares (RSS) is lower than 0.001, which indicates that the semivariance models of each soil property are fitted (Table 8). With the exclusions of log (OM) in 2007, which was best fitted with the spherical model, and log (pH) in 2007 and 2017, which was best fitted with the Gaussian model, the soil chemical properties were best fitted with the exponential model in 2007 and 2017. The stronger the spatial dependence in one soil property, the lower the C0/C0+C value and the less random factor interruption (such as changes in cropping management and fertilization) experienced by the soil property. The spatial dependencies for all soil chemical properties were at the middle level, except AP (2007), pH (2007), and AK (2017). The value of C0/C0+C of these soil properties ranged from 0.350 to 0.500, which means that they were affected by both stable factors, such as topography and soil parent materials, and random factors, such as farming activities and climate. The nugget/sill values of AP (2007), pH (2007), and AK (2017) were lower than 0.17, which means that these factors were high in spatial dependency and were most likely affected by the natural environment. The semivariogram range value is the maximum distance within which autocorrelation or spatial dependence exists [53]. indicates that the spatial autocorrelation of OM and AK increased while that of AK and pH decreased. The Z value of all of the properties was higher than 1, which means that the soil chemical properties in Hefeng were spatially autocorrelated. In conclusion, the results of semivariograms, Moran's I index, and FD are consistent (Table 8, Figure 4).

Spatial Distribution of Soil Properties
In this study, based on the classification of soil chemical properties in the SSSSC (Table 3) In 2007, approximately 0.30%, 0.51%, 5.89%, and 23.63% of the land were at the lowest, low, lower-medium or upper-medium level of soil AP, respectively. Moreover, 42.63% and 27.04% of the land were at the high and highest soil AP levels, respectively. In 2017, the proportion of land at the high level of soil AP increased by 26.15% of the total area in Hefeng, and the land area at the other level decreased. The CV for soil AP in Hefeng County revealed high-to-moderate variability (34.83% in 2007, 25.44% in 2017; Table 3).
The spatial distribution patterns of soil organic matter in 2007 and 2017 were inconsistent ( Figure 5), indicating that soil organic matter content was more likely influenced by the natural environment and human activities, such as fertilizer application and agricultural management ( Table 8) Table  4).
In 2007 (Figure 6a), the soil AP in Hefeng County was lower in the southwest, southeast and north, including Taiping (TP), Tielu (TL), and northern Wu Yang (WY), followed by Rongmei (RM). The soil AP in Hefeng County was higher in Yanzi (YZ), Wuli (WL), Xiaping (XP), Zhongyin (ZY), and Zouma (ZM). In 2017 (Figure 6b), the soil AP was lower in northern Dongyang, southern Rongmei (RM), central Zouma (ZM), and eastern Taiping (TP). This part of the soil should be supplemented with phosphorus fertilizer; otherwise, the AP may become a limiting factor for crop growth. The soil with lower AP or higher AP in 2017 was associated with sandy loam or medium loam, respectively (Table 6).   Table 4).
In 2007 (Figure 6a), the soil AP in Hefeng County was lower in the southwest, southeast and north, including Taiping (TP), Tielu (TL), and northern Wu Yang (WY), followed by Rongmei (RM). The soil AP in Hefeng County was higher in Yanzi (YZ), Wuli (WL), Xiaping (XP), Zhongyin (ZY), and Zouma (ZM). In 2017 (Figure 6b), the soil AP was lower in northern Dongyang, southern Rongmei (RM), central Zouma (ZM), and eastern Taiping (TP). This part of the soil should be supplemented with phosphorus fertilizer; otherwise, the AP may become a limiting factor for crop growth. The soil with lower AP or higher AP in 2017 was associated with sandy loam or medium loam, respectively (Table 6).
The distribution of soil pH levels in Hefeng County was generally acidic, with pH levels at a neutral level in the southeast, north, and southwest in 2007 (Figure 8a). In 2017, 51.40% of the soil pH was at a neutral level, mainly in Zhongyin (ZY), Xiaping (XP), Rongmei (RM), and Zouma (ZM); 34.46% of the soil was weakly acidic. Soil below pH 5.5 was located in the northeastern part of Zhongyin (ZY), and the land of soil pH below 5.5 decreased from 50.37% of the total area in 2007 to 12.79% of the total area in 2017 (Figure 8b).

Discussion
Previous research shows that land-use change alters the ecosystem, resulting in iation in soil properties. Many factors influence the condition of soil properties, inclu human activities and natural elements. Natural conditions, such as soil parent mate topography, climate, and soil type, determine the main characteristics of soil prope Human activities, such as arable land management, including fertilizer application, cropping, and irrigation conditions, influence soil spatial-temporal varia [10,11,17,54,55]. Moreover, the methods of soil sample collection and processing als fect the precision of soil chemical property analysis by uncontrollable mistakes [5,56]

Possible Reasons for Soil Chemical Property Changes
Organic matter content is one of the core indicators that sustain soil life, and it essential determinant of soil fertility levels. According to the Hefeng soil chronicles before 1970, farmers mixed river and lake mud or livestock manure with the soil to prove soil fertility. After reform and development, chemical fertilizers began to be u and by 1990, soil fertility had declined due to the heavy use of chemical fertilizers. finding shows that the average soil OM in Hefeng was previously at a high level ( Ta  3 and [58]. In their study, the spatial heterogeneity of soil OM in western bei, where Hefeng is located, was high. There are multiple possible reasons for the cha in soil OM from 2007 to 2017. The first likely reason is the decrease in farmland-use in sity. According to the yearbook of Enshi [59,60], the multiple crop index (the acreag sowing divided by the acreage of arable land) dropped from 2.83 in 2007 to 2.39 in 2 which means that the frequency of agricultural activities has also dropped. Soil OM be influenced by fallow frequency [61,62]. The second possible reason is that, in re

Discussion
Previous research shows that land-use change alters the ecosystem, resulting in variation in soil properties. Many factors influence the condition of soil properties, including human activities and natural elements. Natural conditions, such as soil parent material, topography, climate, and soil type, determine the main characteristics of soil properties. Human activities, such as arable land management, including fertilizer application, field cropping, and irrigation conditions, influence soil spatial-temporal variation [10,11,17,54,55]. Moreover, the methods of soil sample collection and processing also affect the precision of soil chemical property analysis by uncontrollable mistakes [5,56].

Possible Reasons for Soil Chemical Property Changes
Organic matter content is one of the core indicators that sustain soil life, and it is an essential determinant of soil fertility levels. According to the Hefeng soil chronicles [57], before 1970, farmers mixed river and lake mud or livestock manure with the soil to improve soil fertility. After reform and development, chemical fertilizers began to be used, and by 1990, soil fertility had declined due to the heavy use of chemical fertilizers. Our finding shows that the average soil OM in Hefeng was previously at a high level (Tables 3 and 5 Soil rich in OM (mean value of OM higher than 34.56 g/kg) may explain the type of terrain in Hefeng County, where elevation is high and evaporation is low. With high humidity in Hefeng, soil microbial activity is hindered. This result is similar to that of Yu Fang et al. (2019) [58]. In their study, the spatial heterogeneity of soil OM in western Hubei, where Hefeng is located, was high. There are multiple possible reasons for the changes in soil OM from 2007 to 2017. The first likely reason is the decrease in farmland-use intensity. According to the yearbook of Enshi [59,60], the multiple crop index (the acreage of sowing divided by the acreage of arable land) dropped from 2.83 in 2007 to 2.39 in 2017, which means that the frequency of agricultural activities has also dropped. Soil OM may be influenced by fallow frequency [61,62]. The second possible reason is that, in recent years, the Hefeng County agricultural department has implemented a soil organic matter improvement subsidy program to improve soil fertility, and it strongly supports organic fertilizers [63]. These actions may have changed the soil OM content. The third reason may be the variation in cropping systems. According to the yearbook of Enshi [59,60], the acreage of each crop has changed, in which the grain area and sowing area for vegetables have risen from 239.1 to 243.1 km 2 and from 75.6 to 85.1 km 2 , respectively. The variation in sowing acreage of different crops was affected by the agricultural market and national policies [64,65]. The soil in crop systems, including maize, rice, tobacco leaves, and vegetables, has a relatively high OM content (Table 7, Figure 3), which is attributed to the vegetation residual and cropping management [66,67]. The frequency and intensity of activities such as fertilization, seeding, and loosening of farmland soil differ by crop, thereby causing the differences in soil OM in farmland. The increased soil OM may also be related to the increased acid rain in western Hubei, according to Wu et al.'s research (2016) [55].
Average soil AP value increased from 39.56 mg/kg in 2007 to 58.51 mg/kg in 2017. This trend corresponds with the soil AP in paddy fields in southern China [11]. A total of 86.72% of the land in Hefeng had AP content higher than 40 mg/kg in 2017. The Zouma (ZM) phosphate mine in Hefeng has a reserve of 1.178 billion tons, making it one of the four largest phosphate mines in China [68]. The wastewater, exhaust gas, and sludge from phosphate mines have a negative impact on the surrounding environment, especially in mining areas with steep slopes and topography, which are prone to landslides and severe soil erosion during mining, and this could consequently pollute the surrounding land [69]. According to the yearbook of Enshi [59,60], AP fertilizer application decreased from 6233 tons in 2007 to 5488 tons in 2017, but the pollution is still severe.
Soil in Hefeng was at the upper-middle level of AK with a mean value at 147.82 mg/kg in 2007, and it dropped to 118.02 mg/kg in 2017. We learned from a survey of farmers in Hefeng that although little potash is contained in the straw returned to the field, farmers have long neglected the use of potash, leading to a reduction in soil AK. The low application rate of K fertilizer is also reported in previous research [11,70]. According to these studies, intensive farming was also one of the main causes of soil potassium depletion [71].
Soil pH affects soil microbial activity, plant growth, and nutrient transfer [72]. However, in our study, soil pH in Hefeng is poorly correlated with AK, AP, and OM ( Table 4). The influence of coal mines and the pollution of exhaust vapor from factories (such as sulfur dioxide from sulfur and phosphate fertilizer plants) cause low soil acidity in nearby land. For example, near the sulfur plant in Erdeng village, the soil acidity is lower than 4.5, and no plants grow there except for millennium dwarfs and tea trees [57]. Soil acidity affects the effectiveness of nitrogen, phosphorus, and potassium nutrients. The AP in the soil is highest when the pH value ranges from 6.0 to 7.5, and AK is highest when the pH value between 6 and 10. When the soil is too acidic, plant roots rot [72].
Soil acidification is the main problem in Hefeng County, which is different from soil pH in the grain production area of China's subtropical plain, where soil is alkaline [11]. However, soil pH changes in the study area and the soil in central China's plain area were similar. The pH value of these two regions increased during the period of 2007-2017 (Table 3) [11]. The average soil pH in Hefeng increased by 0.63 or 11.75% (from 5.36 to 5.99) ( Table 3), and the soil changed from acidic to weakly acidic. The reason for this phenomenon is because the pH of rainfall in Hubei province has increased. According to Wang et al.'s (2016) research [73], acid rain (the pH value of the rain less than 5.6) frequency in the southwestern part of Hubei province dropped by more than 50% during the period of 2007-2014. The acreages of ecological tea plantations in Hefeng County are similar to those of arable land, and tea plantations require acidic soil. Our field survey [71] showed that in recent years, farmers have been neutralizing the land where soil acidification is severe. Local farmers use lime to moderate the soil's acidic quality and ease the problem of soil acidification, thereby allowing the soil to become loose and breathable, and providing it with an increased ability to retain water, warmth, and fertilizer.
We explored the interaction among soil chemical properties, topographic parameters, and soil texture in Hefeng County using SPSS 25. The finding shows ( Table 4) that in 2007, soil pH did not exhibit a correlation with OM, AP, or AK. Soil OM was significantly positively correlated with AP (r = 0.192, p < 0.01) and AK (r = 0.245, p < 0.01). Soil AP was significantly positively correlated with AK (r = 0.241, p < 0.01). In 2017, AP and AK no longer correlated well with OM. AP was poorly correlated with AK. The correlations between other indices were significant but generally weak (p < 0.05, R2 < 0.6). The soil OM (r = −0.197, p < 0.01) and soil AP (r = −0.043, p < 0.01) decreased at higher altitudes but exhibited higher AK. The soil located in steeper land had lower AP and AK. The soil in areas with high topographic relief, profile curvature, and planform curvature had less AP and AK, and lower pH ( Table 5).
The OM and soil AK in 2007 were significantly correlated (p < 0.05) with soil texture (p < 0.05), while soil AP and soil pH in 2007 and all soil chemical properties in 2017 were poorly correlated with soil texture ( Table 6). The OM content was highest in medium loam and was lowest in sandy soil. The soil AP content was highest in medium loam and lowest in sandy soil. The correlation between soil textures and soil chemical properties decreased from in the period of 2007-2017. The cropping systems played an important role in affecting soil chemical properties in 2007 (p < 0.01), while they only significantly affected AK in 2017 (p < 0.05).

Suggestions for Local Farming Management
Since soil chemical properties not only interact with each other, but are also influenced by soil texture, soil cultivation activities, cropping systems, and topographical factors, this research leads to the following suggestions:

Moderate Use of Fertilizer
Fertilizing according to the soil condition-it is advisable to apply semi-rotten fertilizers on sandy soils, and quick-acting fertilizers should be applied several times in small amounts to prevent fertilizer shedding.
Maintaining and improving soil organic matter [62]-the application of organic fertilizers or the application of organic fertilizers with chemical fertilizers can maintain the soil acid-base balance and slow soil acidification, because organic fertilizers contain alkaline substances. However, as livestock manure may contain pollutants such as heavy metals and antibiotics [74,75], care should be taken to avoid environmental and health risks when choosing organic fertilizers and to avoid the addition of harmful substances into acidic soils.
Neutralizing acidified soil-the common method to improve acidic soil is to apply alkaline substances such as lime [76] or straw-return to the field [77,78], but the cost of this is high [79]. From the long-term perspective, large applications of lime can lead to soil slumping and nutrient imbalance, because lime only provides the nutrient calcium, and large amounts of calcium can decrease soil effectiveness of magnesium (Mg), potassium (K), and phosphorus (P) [80].
Despite the slow decrease in the amount of phosphorus fertilizer applied in recent years (from 13,000 in 2007 to 9900 t in 2017, according to the yearbook of Enshi [59,60]), it is still applied in generally large quantities. The excessive use of phosphorus fertilizer leads to increased costs for farmers and soil environmental pollution problems; thus, it is recommended to use phosphorus fertilizer sparingly to reduce the accumulation of soil AP. Moreover, it is also necessary to control soil contamination caused by factories [81].
At present, there is a problem of soil acidification in Hefeng County, and it should be noted that the tea industry in Hefeng County ranked number one among the tea industries in Hubei Province. The soil of tea plantations is acidic; thus, soil acidification in this study area should be considered with the local situation when improving the pH of the soil.

Crop Selection According to Soil Suitability and Topographic Factors
While acidic (or alkaline) fertilizer can neutralize the properties of soil, it has a high cost [82]. In acidic soil, farmers can plant tea, oil tea, sweet potatoes, yams, buckwheat, and citrus fruits [83]. In land with different topographic relief properties, food crops can be planted on flat fields. In middle hills, farmers can plant tea, oil tea, and other specialty products. In general, famers can plant green manures and winter crops, increase potassium fertilizer, and coordinate the nutrient ratio in order to improve soil fertility.

Construction of Water Conservancy Facilities
Improving water irrigation facilities can reduce the changes in drought and flooding, and also curb soil erosion on sloping land. The shallower the cultivation layer, the thinner the soil layer. For places that are less steep, terraces can be built.
Although the amount of applied potassium fertilizer is high (an increase of from 1753 tons in 2007 to 4505 tons in 2017), soil erosion in mountainous areas may lead to potassium loss [84] (Table 5). Soil potassium remains a limiting factor for crop growth in Hefeng.
Local agricultural departments should help famers to choose proper amounts of fertilizer and construct more facilities for farmland water conservation. The government should regulate mining facilities to prevent them from negatively affecting soil.

Conclusions
This research investigated the spatial changes in soil OM, AP, AK, and pH content in Hefeng in 2007 and 2017 by using the geostatistical methods incorporated in Gs + 9.0 and ArcGIS 10.7.
The results show that the mean value of soil OM, AP, AK, and pH in Hefeng increased from 33.64 g/kg, 39.56 mg/kg, 147.82 mg/kg, and 5.36 in 2007 to 35.49 g/kg, 58.51 mg/kg, 118.02 mg/kg, and 5.99 in 2017, respectively. The spatial variation and dependency of each soil property were at the middle level, which means that they were affected by both stable factors, such as topography and soil parent materials, and random factors, such as farming activities and climate.
In conclusion, the findings of this study in combination with those of others suggest that the topography and soil texture, acid rain, and farmers' fertilization practices in Hefeng are all related to the chemical properties of soil and its temporal and spatial variation.
Because of the different soil survey standards in Hefeng, only a few items of variable soil chemical data can be represented for a long-time scale. Our research is limited to analyzing the impact of soil properties from the climatic prospective due to lack of temperature and precipitation data at different elevations in Hefeng County. Moreover, we did not analyze the spatial variance of the correlation of topological factors, soil types, soil texture, and soil chemical properties. Further work could be conducted by quantifying the relationship between changes in soil chemical properties with climatic factors and social economic factors. However, the existing data are still essential for understanding soil chemical property changes in this subtropical mountainous region. With stricter procedures of soil sampling, testing, and additional regular soil surveys and soil monitoring, our future studies in Hefeng will be more in-depth, precise, and complete.