Semicentennial Tillage Significantly Affects the Soil Evolution in Arid Regions of China

1. Introduction

Arid regions constitute approximately 20% of the Earth’s terrestrial surface and are among the ecosystems most vulnerable to climate change (Lian et al. 2021). As the largest country under threat of desertification (Guo et al. 2022), China features decertified land spanning over 2.6 million km², accounting for 27% of the whole territory (Ren et al. 2024). These areas are predominantly located in extremely arid, arid, and semiarid zones, where water resources are critically scarce (Wei et al. 2021).

The scientific reclamation of desert areas is a primary strategy to mitigate desertification and has been widely implemented in arid regions (Ma et al. 2024). While the natural weathering of soil is a prolonged process, conventional and conservation tillage practices significantly influence soil microstructure, moisture and thermal dynamics, nutrient content, and the composition of soil organic matter. These practices collectively improve soil quality in newly reclaimed cultivated land (Liu et al. 2021; Man et al. 2021). Numerous studies have demonstrated the positive impact of reclamation on desert soil quality. Firstly, reclamation reduces the risk of soil degradation and enhances resistance to wind erosion (Fallahzade et al. 2020). Secondly, it improves the physical properties of desert soil by decreasing bulk density and increasing the proportions of clay and silt particles, as well as soil organic carbon (SOC) levels. These changes enhance the soil’s ability to retain nutrients and water (Cao et al. 2021; Shang et al. 2019). Lastly, short-term agricultural activities following reclamation stimulate soil microbial activity and significantly increase nitrogen (N), phosphorus (P), and other nutrient levels (Chen et al. 2020; F.-R. Li et al. 2018). However, soil quality improvements from reclamation are not indefinite. According to Hua et al. (2008), desert soils stabilize 6–8 years after reclamation, and soil quality often begins to decline after 10 years. Furthermore, the long-term effects of reclamation, particularly for pure desert soils, remain uncertain. Most existing studies focus on semi-shrubby deserts, grasslands, and oasis forests in regions with limited agricultural histories, suggesting that at least 50 years may be required for desert soils to transform into highly productive agricultural land (Su et al. 2010). Short-term observations may fail to capture the slow and complex processes of soil evolution, obscuring the long-term trends following reclamation (Bacq-Labreuil et al. 2021; Ma et al. 2024).

Arid and semiarid regions are particularly vulnerable to desertification due to low and highly variable precipitation, high evaporation rates, nutrient-poor soils, sparse vegetation, and fragile ecosystems (Maestre et al. 2021). These regions also feature saline soils with low moisture content and extremely high electrical conductivity (Hachmi et al. 2018; Mihi et al. 2019). Such physicochemical constraints limit the transformation rate of organic matter, leading to the accumulation of untransformed materials. This phenomenon negatively affects soil fertility, impeding the recycling of essential nutrients required for plant growth (Hall et al. 2011; Sinsabaugh et al. 2009; Y. P. Wang et al. 2010). Given these complexities, research focusing on long-term agricultural chronologies is essential to comprehensively understand the impact of desert reclamation on soil properties and to discern the gradual evolution patterns of reclaimed desert soils.

Few studies have focused on subsoil in arid regions due to its limited development and low particle cohesion, which complicate the application of conventional investigative methods. In this context, digital soil mapping (DSM) offers a promising alternative for predicting the three-dimensional distribution of various soil properties. Over the past two decades, DSM techniques have gained significant traction for identifying spatial soil patterns and their associated characteristics (McBratney et al. 2003; Minasny and McBratney 2016; Poggio et al. 2021). Among these techniques, machine learning algorithms have become increasingly popular due to their capacity to model complex interactions between soil properties and environmental variables with high accuracy, particularly when applied to large soil datasets (Arrouays et al. 2020; Heuvelink and Webster 2022). Several studies have shown that machine learning techniques can accurately predict SOC at various scales and in a variety of ecosystems, including semiarid and arid regions (Wiesmeier et al. 2011; Zeraatpisheh et al. 2019). These DSM products are essential for enhancing soil management and conservation strategies, particularly given the vulnerability of semiarid regions to climate change and anthropogenic pressures (Azamat et al. 2024). Moreover, soils exhibit significant temporal and spatial variability, necessitating an investigation of vertical distributions of soil properties (Poggio and Gimona 2014). Each soil property follows a unique vertical distribution pattern, varying continuously with depth within the profile. These patterns can be modeled using depth functions, which have evolved from the freehand curves pioneered by(Bishop et al. 1999; Jenny 1941)) to more sophisticated mathematical models. Examples include polynomial functions, power functions, equal-area quadratic spline functions (EAQSFs) (Bishop et al. 1999), exponential decay functions (EDFs) (Minasny 2006), and sigmoid functions (Y. Zhang et al. 2017).

When combined with machine learning algorithms, these depth functions have proven effective in modeling the three-dimensional distribution of soil properties, as demonstrated by studies on SOC (Allory et al. 2022; Fu et al. 2024). However, research on other soil properties, particularly in arid and semiarid regions, remains limited. Therefore, further investigation is urgently needed to evaluate and optimize the performance of existing modeling approaches for predicting a broader range of soil properties in these challenging environments.

In this study, a variety of soil properties were selected. Through correlation analysis, four soil properties, SOC, TP, NH₄⁺ and pH, were ultimately identified as the most critical soil properties. These properties were selected to elucidate the variation trends with depth and their spatial distribution characteristics, providing insights into the effects of semicentennial tillage on soil evolution. To determine the best approach for predicting soil properties in the study area in Xinjiang Province, Northwest China, three machine learning models—the multiple linear regression (MLR), quantile regression forest (QRF), and boosted regression tree (BRT) methods—combined with various depth functions were investigated. The objectives of this study were to (1) predicting and comparing the spatial distributions of soil properties (SOC, TP, NH₄⁺ and pH) via various DSM approaches, and (2) verify the impact of semicentennial desert reclamation on soil properties and soil evolution.

2. Materials and Methods

2.1. Description of the Study Area

The study area is situated in southern Xinjiang, China (81°17′–81°22′ E, 40°28′–40°31′ N) (Figure 1), at the confluence of the Hetian, Aksu, and Yerqiang Rivers. It lies north of the Tarim River and south of the Taklimakan Desert. The region experiences high annual sunshine exposure, with a mean duration of 2303.1 hours. Due to its arid climate, characterized by low precipitation (mean annual precipitation of 60 mm) and high evapotranspiration (1800 mm), substantial amounts of water from the Tarim River have historically been utilized to reclaim desert land and establish oasis agriculture, beginning with the founding of Aral City. Over time, the city’s reclaimed cultivated land has continued to expand, and the desert edge has continued to move towards the south. However, the effects of this reclamation on the three-dimensional distributions of the critical soil properties in this area are unknown.

2.2. Experimental Design and Soil Sampling

We selected the 12th Regiment farm, which has been cultivated in natural sandy lands for more than 50 years; additionally, six adjacent, noncultivated virgin sandy land sites served as a reference (control). Previous studies in this area revealed that the soil properties in the cropland areas before cultivation were similar to those in the adjacent uncultivated sandy areas.

In August 2019, 71 sampling sites were distributed via a regular grid layout on the farm, with intervals of 1 km from east to west and 500 m from north to south. Each sampling site was drilled to a depth of 1 m, and two soil samples from 0–10 cm, 10–20 cm, 20–40 cm, 40–60 cm, 60–80 cm and 80–100 cm were collected. The soil samples from the same level at different drilling points were mixed, and visible roots and litter debris were removed from each soil sample. Finally, the samples were air dried and sieved through a 2 mm soil sieve for physicochemical analysis.

2.3. Analysis of Soil Physical and Chemical Properties

The pH of the air-dried soil was determined by extracting 5 g of soil with 0.01 M CaCl₂ (1:5) in an end-over-end shaker for 1 h followed by centrifugation at 492 × g for 10 min. The pH of the supernatant was determined via a Thermo Orion pH meter (Thermo Orion 720A+, Beverly, MA, USA). The soil organic carbon (SOC) content was determined via the K₂Cr₂O₇-H₂SO₄ oxidation method of Walkley-Black (Nelson and Sommers 1982). The total phosphorus (TP) content was determined calorimetrically after wet digestion with H₂SO₄-HClO₄ (Parkinson and Allen 1975). Ammonium-N (NH₄⁺) was determined with an AA3 Continuous Flow Analytical System (Bran+Luebbe, Norderstedt, Germany) following extractions of fresh soil with 1 mol L⁻¹ KCl.

Notably, several soil properties, including TN, TK, and CEC, were selected. We carried out correlation analysis on all these soil properties and selected the four with the lowest correlation.

2.4. Environmental Covariates

One of the primary elements influencing the formation of soil in arid and semiarid areas is topography (Zeraatpisheh et al. 2019). A digital elevation model (DEM) with a cell size of 12.5 m x 12.5 m was thus downloaded from ALOS PALSAR in order to obtain the topographic attributes (“ASF Data Search” n.d.). Elevation (El), topographic wetness index (TWI), slope (SL), slope position (Slpposi), profile curvature (PrCu), plan curvature (PlCu) and distance to an irrigation canal (DIC) were among the terrain variables derived from the DEM.

The normalized difference vegetation index (NDVI), carbonate index (CAEX) and salinity index (SI) were among the auxiliary variables used in remote sensing. The environmental covariates were obtained using the SAGA GIS (System for Automated Geoscientific Analysis), and predictors with varying resolutions were resampled to 10 m using the bilinear interpolation approach.

2.5. Modeling Approaches

Two nonlinear machine learning models, namely, QRF, BRT, and MLR, combined with different depth functions, such as exponential functions, logarithmic functions, power functions and cubic polynomial functions, were used for the digital mapping of soil properties (SOC, TP, NH₄⁺ and pH).

BRT analysis is a parametric data mining technique that can handle both nonlinear and linear relationships (Myles et al. 2004), and it has been reported that this technique has been widely used in the field of DSM (Taghizadeh-Mehrjardi et al. 2014).

A regression forest (RF) is a classifier or regression model made up of numerous decision or regression trees, each of which is reliant on the values of a randomly selected vector that is sampled separately and has the same distribution as all the other trees in the data. (Liaw and Wiener 2002). (Meinshausen n.d.) suggested altering the RF procedure’s outputs by permitting the calculation of the targeted variables’ prediction intervals. The QRF technique takes into account the spread of the response variable from which prediction intervals are formed, whereas the RF technique merely preserves the mean of the observations that fall into each node and ignores all other information for each node and each tree.

Spatial prediction has long employed linear regression, one of the first statistical methods (Draper and Smith 1998). Multiple linear regression fits the data to a straight line (or plane in an n-dimensional space, where n is the number of independent variables) (Mbagwu and Abeh 1998). As a result, each soil property has been given a quantitative estimate using the MLR model.

Pearson correlation coefficients between the environmental covariates and soil parameters were computed in order to identify significant environmental covariates. Each soil attribute was then predicted using only relevant factors (Aksoy 2012). Based on variable relevance, the most significant variables were determined for each model. In many applied models, the quantification of variable importance is a crucial topic (Genuer et al. 2010). It helps in the interpretation of the variables and their impact on model correctness (Genuer et al. 2010) as well as the identification of the most significant environmental covariates utilized in the models (Nauman and Thompson 2014). The “caret” package in R 4.3.3 and RStudio was used for all modeling and variable significance computations.

2.6. Statistical and Validation Strategy

Leave-one-out cross-validation (LOOCV) was performed to evaluate the prediction methods. In order to ensure that every data point appears in the test set at least once, cross-validation offers a framework for generating multiple training/test splits in the dataset. This approach has the benefit of impartiality and consistent performance on smaller datasets.

The model performance was assessed using three validation criteria: the coefficient of determination (R²), the mean absolute error (MAE), and the root mean square error (RMSE).

$$RMSE=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n{\left(P_i-O_i\right)}^2}{n}}}$$

$$MAE={\displaystyle \frac{{\sum }_{i=1}^n\left|P_i-O_i\right|}{n}}$$

$$R^2={\left({\displaystyle \frac{{\sum }_{i=1}^n\left(O_i-O_{avg}\right)\times \left(P_i-P_{avg}\right)}{\sqrt{{\sum }_{i=1}^n{\left(O_i-O_{avg}\right)}^2}\times {\left(P_i-P_{avg}\right)}^2}}\right)}^2$$

where P_i, O_i, O_avg, P_avg, n, and p stand for the model’s total number of explanatory variables, number of data points, average of the observed and predicted soil property values at the ith point, and predicted and observed values, respectively.

3. Results and Discussion

3.1. Correlation Analysis and Descriptive Statistics

We measured a variety of properties and selected four soil properties, SOC, pH, NH₄⁺ and TP, which are related to other soil properties, but their correlations with each other are low (Figure 2). The soil properties not selected, such as TN and SOC (r = 0.93), CEC and SOC (r = 0.67), and TK and CEC (r =0.54), were highly correlated with the four selected soil properties. In this way, we ensured the maximum representativeness and heterogeneity of the selected attributes.

The summary statistics of the soil properties in the study area are presented in Table 1, while those of the CK are presented in Table 2. Scatter diagrams of SOC, pH, NH₄⁺ and TP are shown in Figure 4. Notably, we removed the ten maximum values and ten minimum values from each layer for the sake of aesthetic considerations and to ensure that these depth functions were distinguished to the greatest extent. In the study area, the mean values and standard deviations of the SOC content decreased with depth, and the mean values were 4.49, 4.19, 3.61, 2.69, 2.31 and 2.19 g/kg, respectively. The coefficients of variation (COVs) of the SOC content at the 0–20 cm depth were less than 35%, indicating moderate variation, and those at the 20–100 cm depth ranged from 35% to 43%, indicating high variation. The mean values of pH increased with increasing soil depth and were 8.83, 8.86, 8.91, 8.99, 9.02 and 9.07. The COVs of the pH values at the six depths ranged from 2.54% to 3.11%, indicating low variation. The mean values of NH₄⁺ increased with depth from 0–50 cm and then decreased with depth from 50–100 cm, with values of 3.06, 2.80, 2.43, 2.07, 2.33 and 2.48 mg/kg, respectively. The COVs of the NH₄⁺ content at the 0–20 cm depth were less than 35%, indicating moderate variation, and those at the 20–100 cm depth ranged from 35% to 62%, indicating high variation. The mean values and standard deviations of the TP content decreased with depth, and the mean values were 1.04, 1.05, 0.90, 0.78, 0.71 and 0.70 g/kg, respectively. The COVs of the TP content at depths of 0–20 cm were less than 35%, indicating moderate variation, and those at depths of 20–100 cm were less than 25%, indicating little variation.

In the CK treatment, there was no significant trend in the values of SOC, pH, NH₄⁺ and TP with increasing depth. We plotted the changes in the values of the four soil properties at different depths after semicentennial tillage, as shown in Figure 3, which revealed the vertical changes in the soil properties following semicentennial tillage in desert areas. The observations revealed that the accumulation of SOC is higher than that of the CK at 0–60 cm but lower in the 60-100 cm layer. Additionally, the pH of the reclaimed soils was lower than that of the CK soils, while the NH₄⁺ and TP contents were greater than those of the CK soils in both the topsoil and subsoil. These findings indicate that reclamation has led to increases in N and P contents. Although reclamation promoted an increase in SOC in the subsoil, it also resulted in a decline in the SOC accumulation rate within this layer; notably, this rate decreased with increasing depth.

Figure 3. Histograms of the SOC, pH, NH₄⁺ and TP values in the experimental group (TG) minus the control group (CK).

Figure 3. Histograms of the SOC, pH, NH₄⁺ and TP values in the experimental group (TG) minus the control group (CK).

Figure 4. Fitted vertical distributions of SOC, pH, NH₄⁺ and TP. The red circles indicate the mean values of SOC, pH, NH₄⁺ and TP in each layer. All the presented formulas are inverse functions.

Figure 4. Fitted vertical distributions of SOC, pH, NH₄⁺ and TP. The red circles indicate the mean values of SOC, pH, NH₄⁺ and TP in each layer. All the presented formulas are inverse functions.

In arid and semiarid regions, nutrient accumulation occurs primarily through artificial means. The unique natural characteristics of these areas, particularly the ratio of evaporation to precipitation, pose challenges for the accumulation of essential nutrients such as C, N, and P, which are necessary for the growth of plants (Gamalero et al. 2020).

The contents of SOC are determined primarily by the balance between exogenous crop carbon inputs and the decomposition of existing SOC (F. Zhang et al. 2022). In arid and semiarid regions, differences in temperature between the topsoil and subsoil affect microbial activity, which subsequently influences SOC mineralization at different depths (X. Li et al. 2020; C. Wang et al. 2021; Yu et al. 2022). After reclamation, topsoil has a relatively high clay content that facilitates the formation of mineral-associated organic carbon compounds (Song et al. 2024). The minimal variation observed in SOC at different depths in the CK can be attributed to the high and uniform sand content.

The reduced levels of SOC in the 60–100 cm layer relative to those in the CK soil can be explained by the evaporation rate exceeding the precipitation rate, which limited the downwards migration of soluble organic carbon. The study area is located on the northern edge of the Taklimakan Desert and is characterized by scarce precipitation and rapid evaporation. This environment has led to a gradual accumulation of alkaline ions and soluble salts (Peng et al. 2024). The low vegetation cover results in minimal root activity, which adversely affects the decomposition of organic matter. Cotton, a tap-rooted crop, possesses a main root that can extend to depths of 1–2 m; however, its numerous lateral roots are primarily concentrated within the 0.2–0.6 m range (Kang et al. 2021; Z. Zhang et al. 2020). These lateral roots play crucial roles in promoting SOC decomposition and reducing soil alkalinity, thereby contributing to differences in pH between the 20–60 cm and 60–100 cm soil depths.

NH₄⁺ and TP are both essential nutrients; however, natural processes in desert soils are insufficient for their effective accumulation. In addition to precipitation and evaporation, other soil activities exert minimal influence on the distribution of these nutrients (Shao et al. 2022), which accounts for the limited variation in NH₄⁺ and TP across the CK, with both nutrients showing relatively low concentrations. A high rate of evaporation compared to precipitation typically leads to elevated concentrations of soluble ions in the topsoil relative to the subsoil. However, the lack of significant depth variation in NH₄⁺ and TP in the CK soil may be attributed to moderate rainfall occurring two days prior to sampling. After reclamation, although the surface clay content increased, the soils predominantly retained sandy textures in the profile.

Following irrigation and fertilization, NH₄⁺ was evenly distributed across various soil depths due to water movement. Conversely, phosphorus remains concentrated at or near the surface layer owing to its strong adsorption properties and low solubility, which tend to bind tightly with soil particles. Additionally, under alkaline conditions, phosphate ions readily form insoluble calcium phosphates with calcium ions (Geng et al. 2022; Luo et al. 2023; Yahaya et al. 2023).

3.2. Optimal Depth Functions Selection

Under normal circumstances, there is little difference in the vertical distribution of desert soil properties, whereas watering and fertilization under cultivation produce a series of effects. Consequently, semicentennial tillage is expected to have a longer and more lasting effect on the vertical distribution of desert soil properties. To demonstrate these changes, we used depth functions to simulate the variation trend of desert soil properties with depth after semicentennial tillage. SOC, pH, NH₄⁺ and TP within 1 m were fitted by different depth functions, including the exponential function, power function, logarithmic function and cubic polynomial function. The results revealed that different function types had significantly different fitting effects for each soil property (Table 3 and Figure 4). The binomial exponential function had the best fitting effect for SOC (R²=0.90) and TP (R²=0.77), and the cubic equation had the best fitting effect for pH (R²=0.92) and NH₄⁺ (R²=0.89). These results showed that semicentennial tillage strongly changed the vertical distribution of critical soil properties and transformed the uniform distribution of desert soil properties from topsoil to subsoil, which may affect soil evolution.

3.3. Model Performance

Semicentennial tillage tremendously changed the soil properties in the study area, not only in the soil profiles but also in different regions, and these changes might have different rates with increasing depth. Two nonlinear machine learning models, namely, QRF and BRT, and a linear model, namely, MLR, were used to predict the distributions of the four soil properties. The performances of the models were evaluated via leave-one-out cross-validation (Table 4). On the basis of the validation criteria, the nonlinear models (QRF and BRT) were more accurate than MLR. The findings demonstrated that the QRF and BRT techniques performed the best in predicting every property among the models under study.

The overall trends of the SOC distribution at the six depths were significantly different, with a spatial pattern of low values in the southeast and west and high values in the middle and north at depths of 0–40 cm. At depths of 40–100 cm, the spatial distribution pattern is characterized by high values in the southeast and west and low values in the middle and north, which was completely opposite to the predicted results for the surface layer (Figure 5).

The overall variation trends of the spatial distributions of pH at the six depths were obviously different, and the pH was generally low at depths of 0–20 cm, roughly following the distribution characteristics of farmland. At depths of 20–100 cm, the spatial distribution pattern is characterized by low values in the southeast and west and high values in the middle and north (Figure 6).

There was no significant difference in the overall trend in the NH₄⁺ spatial distribution at the six depths (Figure 7). The overall change trend of the spatial distribution of TP at the six depths was consistent, with low values in the southeast and west and high values in the centre and north. The distribution characteristics at 0–40 cm were roughly consistent with those of farmland (Figure 8).

3.4. Confidence Interval of Prediction

The QRF was used to map the predicted uncertainties. The QRF can be used to estimate the value of a target variable for any required quantile. The upper (95%) and lower (5%) quantile maps can then be used to compute prediction intervals (e.g., 90%). The 5% quantile value indicates that 5% of all possible predictions are lower than this value. This means that under the same conditions, one can predict with 95% confidence that the result will be higher than this value. The 95% quantile value means that 95% of all possible predictions are lower than this value. This means that under the same conditions, one can predict with 5% confidence that the result will be higher than this value. In this study, we used 5% and 95% to represent the confidence interval of prediction (Supplemental Figs. 1-8).

4. Conclusions

Through correlation analysis, this study selected the four soil properties (SOC, pH, NH4+, and TP) with the lowest correlations with each other and high correlations with the properties that were not selected to represent the critical soil properties in the study area. Depth functions were used to assess the variation trends with depth, and nonlinear and linear techniques were used to predict the spatial distributions of SOC, pH, NH4+, and TP. The results showed that the binomial exponential function has the best fitting effect for SOC and TP, and the cubic equation has the best fitting effect for pH and NH₄⁺. The QRF model has the best performance in predicting the SOC, pH, NH₄⁺, and TP contents.

According to the predictions, all critical soil properties changed with depth along the profile, which revealed that semicentennial tillage had a long and lasting effect on the vertical distribution of critical soil properties. In arid and semiarid regions, different temperatures between topsoil and subsoil affect microbial activities, which subsequently influence SOC mineralization at different depths. The clay content, evaporation rate, low vegetation cover and cotton lateral roots played crucial roles in promoting SOC decomposition and reducing soil alkalinity, thereby contributing to variations in the pH between the 20–60 cm and 60–100 cm soil depths. Additionally, under alkaline conditions, phosphate ions readily form insoluble calcium phosphates with calcium ions. Irrigation and fertilization strongly affect the movement of soil water, which results in an even distribution of NH₄⁺. All the changes in soil properties after semicentennial tillage confirmed the occurrence of rapid soil evolution if deserts are cultivated. Thus, desert reclamation was scientifically confirmed to be suitable for agricultural use, which will ease the food production crisis, protect the ecology, and promote soil evolution.

Upon reflection, we conclude that the study area could be expanded, and more chronological observations should be carried out in the expanded area to assess the effects of desert reclamation on soil properties. All of these factors will lead to more precise conclusions and provide more specific guidance for local agricultural development.