Soil Salinity Mapping of Croplands in Arid Areas Based on the Soil–Land Inference Model

Abstract

Soil salinization can decrease soil productivity and is a significant factor in causing land degradation. Precision mapping of salinization in agricultural fields would improve farmland management. This study focuses on the cropland in the Manas River Basin, located in the arid region of northwest China. It explores the potential of a soil mapping method, the Soil–Land Inference Model (SoLIM), which only requires a small number of soil samples to infer soil salinization of farmlands in arid areas. The model was utilized to create spatial distribution maps of soil salinity for the years 2009 and 2017, and changes in the distribution were analyzed. The research results indicate: (1) Through the analysis of sample point data, it was observed that soil salinity in the study area tends to accumulate in the surface layer (0–30 cm) in spring and in the subsoil layer (60–90 cm) during the crop growing season, with significant spatial variability. Therefore, it is necessary to conduct detailed salinity mapping. (2) Using field measurements as validation data, the simulation results of the SoLIM were compared with spatial interpolation methods and regression models. The SoLIM showed higher inference accuracy, with R² values for the simulation results of the three soil layers all exceeding 0.5. (3) The SoLIM spatial inference showed salt accumulation in the northern part and desalination in the southern part. The findings of this study suggest that the SoLIM has the potential to effectively map soil salinization of croplands in arid areas, offering an efficient solution for monitoring soil salinity in arid oasis croplands.

1. Introduction

Soil salinization refers to the accumulation of soluble salts in the surface layer of the soil due to evaporation and transport of water [1]. Soil salinization can lead to land degradation, reduced crop yields, and deterioration of the ecological environment [2,3,4], and has become one of the major environmental problems worldwide. Soil salinization affects about 3% of the world’s land resources and more than 20% of irrigated cropland [5,6]. It is prevalent in coastal wetlands and arid–semi-arid regions [7], with China’s Xinjiang region serving as a typical arid area, where saline soils account for about 60.6% of the total saline soil area in China [8]. The accumulation of soil salinity in oasis farmlands in arid regions is closely related to crop yield [9,10]. Therefore, an accurate assessment of the spatial distribution of soil salinity can provide valuable insights for the rational allocation of soil resources and management of soil salinity.

Currently, soil salinity mapping methods mainly include spatial interpolation methods based on spatial autocorrelation theory and regression-based methods. Spatial interpolation methods have been widely used due to their simple principles and fast computation speed [11,12,13,14]. With the rapid development of remote sensing technology, the difficulty of detecting environmental factors affecting soil salinization has also decreased. Therefore, linear regression and machine learning methods have been extensively developed and applied in soil salinity mapping [15,16,17,18]. In general, the simulation accuracy of both spatial interpolation and regression methods is highly dependent on the quantity and spatial representativeness of the samples [19]. Some studies indicate that traditional methods may struggle to achieve satisfactory accuracy with a limited number of samples, which significantly limits the monitoring of soil salinity in agricultural areas [20]. In recent years, a soil inference method based on the theory of environmental similarity has emerged in soil property mapping. This method uses the geographic similarity between sampling points and prediction points for spatial prediction [21]. It demonstrates good performance with only a small number of sampling points, thus significantly reducing the cost of soil mapping.

Zhu, et al. [22] developed the Soil–Land Inference Model (SoLIM) based on the theory of geographic similarity. This model has been widely applied to infer soil types and soil properties [22]. For example, Yang, et al. [23] used the acquired knowledge of soil–environment relationships and the SoLIM to update the conventional soil-type map of Wakefield in northern Canada. Their study showed that the updated digital soil map contained more spatial details compared to traditional soil maps. Wen, et al. [24] compared the applicability of four methods, including the SoLIM, in the complex terrain of the Loess Plateau in China. The model was found to be more effective in predicting soil organic carbon. Although the SoLIM has been verified in the prediction accuracy of various soil properties, it is less involved in the prediction of soil salinity in arid zones. Therefore, it is necessary to explore its potential for inferring the soil salinity of croplands in arid regions.

The study area of this research is oasis cropland in the Manas River Basin (MRB) in Xinjiang, where soil salinity is a serious problem. We used the SoLIM to simulate the spatial distribution of soil salinity. The accuracy of the predictions was evaluated using measured sample data. In addition, the results were compared with those obtained from spatial autocorrelation models (ordinary kriging and inverse distance weighting) and regression models (multiple linear regression and geographically weighted regression). Finally, based on the spatial predictions generated by the SoLIM, an analysis of the characteristics of soil salinity in different soil types within the study area was performed. The objective of this research was to explore the feasibility of the SoLIM in mapping soil salinity in regions with severe salinization, offering scientific guidance for the precise management of salinized farmland.

2. Materials and Methods

2.1. Study Area

The study area is located in the mountain-front oasis of the Manas River Basin in the arid region of northwestern China’s Xinjiang [25]. It has an elevation ranging from 300 to 600 m and a geographical extent of (85°12′2″–85°48′26″ E, 44°54′41″–44°22′1″ N). The total area covers 25.95 × 10² km². As shown in Figure 1, the study area is divided into two parts by the West Bank Canal, which runs from east to west. The northern part is the Xiayedi irrigation area, bordering the Gurbantunggut Desert and extending to the downstream of the Manas River. The southern part consists of the downstream branches of the Anjihai River, the West Main Canal, the Main Drainage Canal, and the Alkali Drainage Canal, forming the Anjihai Irrigation Area, which connects to the mountain-front oasis of the northern slope of the Tianshan Mountains. The oasis area in the basin is composed of alluvial fans, alluvial plains, and deltas [26]. It is particularly characterized by higher soil salinity in the transition zone between the alluvial fan edge and the central part of the alluvial plain. In this region, the groundwater depth is shallow, the mineralization is high, and the flow is restricted [27]. The regional mean annual temperature is about 8 °C, with an annual precipitation of 100–200 mm and a potential evaporation of 1500–2100 mm [28]. It has characteristics of a dry climate, high evaporation, and large diurnal temperature fluctuations, which belong to a temperate continental climate [29]. The dominant soil types include Gleysols, Arenosols, and Fluvisols. The main crops grown are cotton, with smaller amounts of grape, corn, and wheat cultivation [30]. Since the implementation of the Western Development Strategy and the adoption of subsurface drip irrigation technology in 2000, there has been land expansion and changes in land and water resource use. Widespread flood irrigation has led to a rise in groundwater levels [31], resulting in secondary salinization problems that significantly affect agricultural production and the ecological environment in the basin [32].

2.2. Data Source and Preprocessing

2.2.1. Soil Salinity Data

Field sampling was carried out in May 2009 and July 2017, with a total of 71 and 51 soil samples collected, respectively. Handheld GPS devices were used to determine the latitude and longitude of each sampling point, and information such as surface conditions, irrigation patterns, and vegetation cover types were recorded for each point. At each sampling point, a five-point sampling method was used to collect five soil samples within a 30 × 30 m area, with sampling depths of 0–30 cm, 30–60 cm, and 60–90 cm. After natural air drying, the soil was sieved through a 1 mm mesh, and a 1:5 soil/water ratio was used to prepare the extract solution. Soil salinity was determined by the weight method, and the soil salinity data obtained were used for model inference and validation.

In this study, we referred to the standards set by the Department of Agriculture and Rural Affairs of Xinjiang to categorize the degree of soil salinization [20]: non-saline (SSC < 3 g·kg⁻¹), slightly saline (3 < SSC < 6 g·kg⁻¹), moderately saline (6 < SSC < 10 g·kg⁻¹), strongly saline (10 < SSC < 20 g·kg⁻¹), and extremely saline (SSC > 20 g·kg⁻¹). Two-thirds of the data from each soil layer were used as the modeling set, while the remaining third was designated as the validation set. In 2009, 24 samples were randomly selected for the validation set and 47 samples for the training set. In 2017, 17 samples were selected for the validation set and 34 samples for the training set (Figure 2).

2.2.2. Environmental Variables Data

This study referenced relevant research on soil salinity mapping and selected environmental variables for salinity inference in conjunction with the characteristics of the study area (

Table 1. Environmental variables data information.

Environment Variables	Temporal Resolution	Spatial Resolution	Dataset
LST	8 d	1 km	MOD11A2
Soil Moisture	1 d	1 km	A 1 km daily soil moisture dataset over China based on situ measurement (2000–2020)
Slope		30 m	ASTER-GDEM
TWI		30 m
NDVI	16 d	30 m	Landsat 5 TM, Landsat8 OLI
NIR	16 d	30 m
Soil texture		250 m	SoilGrids250m 2.0
Soil Classification		250 m

In the table, LST represents surface temperature; TWI is terrain wetness index; NDVI is normalized vegetation index; NIR is near-infrared surface reflectance.

). Figure 3 shows the spatial distribution of environmental variables in 2017 (Figure S1 shows environmental variables in 2009). Temperature influences the transport processes of soil water and salt [33], while terrain factors are the main influencing factors on soil formation and development [34]. As soil development is the result of the synergistic effects of environmental variables, the selection of terrain variables with causal relationships was based on the consideration of both hydrological conditions, which are closely related to soil salinity, and the distribution of surface vegetation, which is affected by salinity [35]. Therefore, four categories of environmental variables were considered in this case study: climate, terrain, vegetation, and soil texture. We resampled the spatial resolution of all environmental variables to 30 m.

Climate data: Temperature affects the process of soil water movement and therefore the distribution of soil salinity, so the land surface temperature (LST) was chosen to represent the effect of temperature on soil salinity [36]. On the other hand, the study area receives limited precipitation and experiences high evapotranspiration. Soil moisture, as an important environmental factor, is influenced by long-term subsurface drip irrigation, which leads to leaching of soil salinity with water percolation [37]. Many studies have used soil moisture data to assess land degradation and ecological risks associated with soil salinization [38,39]. Therefore, in this study, we adopted soil moisture (SM) data from Li, et al. [40] as a climatic variable.

Terrain data: Terrain can be calculated directly with a digital elevation model (DEM). According to the characteristics of gentle terrain in this study area and the study of similar areas, the factors such as altitude and aspect with weak changes and small differences in the study area are eliminated, and the slope and terrain humidity index (TWI) are selected as the inference environment variables [18,41,42]. The TWI is calculated as follows:

$$\mathrm{T}\mathrm{W}\mathrm{I}=\ln \left(\frac{\mathsf{\alpha }}{\tan \mathsf{\beta }}\right)$$

(1)

where α represents the accumulated upslope area for a given catchment and β represents the slope. As the study area is predominantly characterized by large flat plains, the multiple flow direction- function of gradient (MFD-fg) [43] was used in this study, which is more suitable for calculating the drainage area on uphill slopes in flat terrain.

Land surface cover data: We selected the normalized difference vegetation Index (NDVI) as an environmental variable to characterize surface vegetation cover [44,45,46]. Salt-affected soils with crust formation have high reflectance in the near-infrared (NIR) band and have been widely used in remote sensing studies of soil salinization [47,48,49,50]. Therefore, we chose it as an environmental variable for inferring soil salinity.

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{E}\mathrm{D}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{E}\mathrm{D}}$$

(2)

where NIR and RED represent the near-infrared and red bands of Landsat imagery, respectively.

Soil texture data: Soil texture controls the movement of water and salts in soil through soil pores [51,52,53]. Therefore, in this study, we selected soil texture and soil type data from the SoilGrids 2.0 product as environmental variables.

2.3. Methods

2.3.1. Mapping Methods for Salinization

To better demonstrate the application potential of the SoLIM in soil salinity mapping, we selected four models from traditional soil mapping methods for comparison. These include ordinary kriging and inverse distance weighting, belonging to spatial autocorrelation models, and multiple linear regression and geographically weighted regression, belonging to regression models. Spatial autocorrelation models have the advantage of simple computation but require a high spatial distribution of sampling data. Regression models generally perform well but have higher demands on the number of sampling points. Below are the basic principles and structures of the five models used in this study.

1.

Ordinary Kriging (OK)

Ordinary Kriging (OK) interpolation is an interpolation method based on spatial autocorrelation. It determines the weight through the semi-variable function, constructs the covariance matrix, accurately describes the spatial correlation between the known positions, and estimates the unknown position through the semi-variable function and the covariance matrix to obtain the optimal estimate [54]. It is suitable for spatial data analysis and has been used in soil salinity mapping for a long time [11,55,56]. The semivariogram is used as a basic tool to test the spatial distribution structure of soil properties. Based on the theory of regionalized variables and intrinsic assumptions [57], the semivariogram is expressed as follows:

$$\mathsf{\gamma }\left(\mathrm{h}\right)=\frac{\mathrm{I}}{2\mathrm{N}\left(\mathrm{h}\right)}\sum _{\mathrm{i}=1}^{\mathrm{N}\left(\mathrm{h}\right)}{\left[\mathrm{z}\left(\mathrm{x}\mathrm{i}\right)-\mathrm{z}(\mathrm{x}\mathrm{i}+\mathrm{h})\right]}^2$$

(3)

where $\mathsf{\gamma }\left(\mathrm{h}\right)$ is the semi-variance, h is the lag distance, z is the soil attribute parameter; $\mathrm{N}\left(\mathrm{h}\right)$ is the number of sample pairs within the distance interval h; $\mathrm{z}(\mathrm{x}\mathrm{i}+\mathrm{h})$ and $\mathrm{z}\left(\mathrm{x}\mathrm{i}\right)$ are the sample values of two points separated by distance h.

2.

Inverse Distance Weighting (IDW)

Inverse distance weighting (IDW) is an interpolation method based on the weight distribution of distance. The weight is assigned by calculating the distance between the unknown point and the known sample. The estimated value of the unknown position is obtained by weighted average of the known sample values. The efficiency is better, and it is widely used in soil salinity mapping [58,59]. The interpolation function is as follows:

$$\mathrm{Z}\left(\mathrm{x}\right)=\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}{\mathrm{Z}}_{\mathrm{i}}}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}}$$

(4)

$${\mathrm{w}}_{\mathrm{i}}={\mathrm{d}}_{\mathrm{i}}^{-\mathrm{u}}$$

(5)

where $\mathrm{Z}\left(\mathrm{x}\right)$ is the predicted value of the interpolation point, ${\mathrm{Z}}_{\mathrm{i}}$ is the salinity value at the known point, n is the total number of known points used for interpolation, ${\mathrm{d}}_{\mathrm{i}}$ is the distance between point i and the predicted point, and ${\mathrm{w}}_{\mathrm{i}}$ is the weight assigned to point i [60]. The closer the distance is, the greater the weight of the point is [61].

3.

Multiple Linear Regression (MLR)

The multiple linear regression (MLR) model is characterized by its simplicity, ease of computation, and interpretability, making it one of the most commonly used models in digital soil mapping [62,63,64]. MLR fits a linear equation with quantitative coefficients using two or more covariates to predict the outcome of the target variable [65]. It assumes a linear relationship between the target variable (Y) and the covariates (x) [66]. The structure of the MLR model can be expressed as:

$$\mathrm{Y}=\mathrm{a}+\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{b}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{i}}+\mathsf{\epsilon }={\mathrm{b}}_1{\mathrm{x}}_1+{\mathrm{b}}_2{\mathrm{x}}_2+\cdots +{\mathrm{b}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{i}}+\mathsf{\epsilon }$$

(6)

where Y is the target variable, x_i represents the covariates, a is the intercept, b_i are the regression coefficients, and ε denotes the regression residuals.

4.

Geographically Weighted Regression (GWR)

Geographically weighted regression (GWR) is a local linear regression model that emphasizes the influence of spatial location on the model for sample points. It fits an independent linear regression model to each sample point [67]. The basic principle is to use a kernel function to partition the dataset into subsets and to calculate the regression equation using a decay function [68]. The number of points for each regression equation is determined by the bandwidth of the kernel function, and the optimal bandwidth is determined by comparing the values of the Akaike information criterion (AIC) [69]. In this study, we initially employed partial least squares regression to examine variables with multicollinearity. Additionally, we used global Moran’s I to assess whether there was spatial clustering of residuals at a global level in the OLS model. The structure of the GWR model can be expressed as follows:

$${\mathrm{Y}}_{\mathrm{j}}={\mathsf{\beta }}_0\left({\mathrm{u}}_{\mathrm{j}},{\mathrm{v}}_{\mathrm{j}}\right)+\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\beta }}_{\mathrm{i}}\left({\mathrm{u}}_{\mathrm{j}},{\mathrm{v}}_{\mathrm{j}}\right){\mathrm{x}}_{\mathrm{i}\mathrm{j}}+{\mathsf{\epsilon }}_{\mathrm{j}}$$

(7)

where $\left({\mathrm{u}}_{\mathrm{j}},{\mathrm{v}}_{\mathrm{j}}\right)$ are the coordinates of the j-th sample. ${\mathrm{Y}}_{\mathrm{j}}$ is the value of the dependent variable for the j-th sample, ${\mathrm{x}}_{\mathrm{i}\mathrm{j}}$ is the observed value of the independent variable set for the j-th sample, ${\mathsf{\beta }}_0$ is the regression intercept at the location $\left({\mathrm{u}}_{\mathrm{j}},{\mathrm{v}}_{\mathrm{j}}\right)$, ${\mathsf{\beta }}_{\mathrm{i}}$ is the regression coefficient at the location $\left({\mathrm{u}}_{\mathrm{j}},{\mathrm{v}}_{\mathrm{j}}\right)$, and ${\mathsf{\epsilon }}_{\mathrm{j}}$ is the random error.

5.

Soil–Land Inference Model (SoLIM)

The Soil–Land Inference Model (SoLIM) is a model used for soil mapping, which is built upon the integration of GIS, fuzzy logic, and expert knowledge [70]. Its theoretical foundation is based on the principle that “the more similar the geographic configuration of two points (regions), the more similar the attributes (values) of the target variable at these two points (regions)”. It applies this principle to infer soil properties at unsampled locations, achieving good practical results [71]. The process of soil attribute mapping based on the SoLIM involves three steps. In the first step, the similarity ${\mathrm{S}}_{\mathrm{i},\mathrm{j},\mathrm{E}}^{\mathrm{k}}$ between the soil environment of the unknown sample point (i,j) and known soil sample point k is determined based on the fuzzy logic theory [72]. In the second step, ${\mathrm{S}}_{\mathrm{i},\mathrm{j},\mathrm{E}}^{\mathrm{k}}$ is calculated as the similarity between the soil salinity at the known sample point k and the soil salinity at the position to be speculated (i,j). Finally, the SSC at the position to be speculated (i,j) is calculated by Equation (8) [73]. The equation is as follows:

$${\mathrm{S}}_{\mathrm{i},\mathrm{j}}=\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}·{\mathrm{S}}_{\mathrm{k}}$$

(8)

$${\mathrm{w}}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}=\frac{{\mathrm{S}}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{S}}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}}$$

(9)

$${\mathrm{S}}_{\mathrm{i}\mathrm{j}}^{\mathrm{k}}=\underset{\mathrm{v}=1}{\overset{\mathrm{m}}{{\displaystyle \mathrm{P}}}} \left({\mathrm{E}}^{\mathrm{v}}\left({\mathrm{e}}_{\mathrm{i}\mathrm{j}}^{\mathrm{v}},{\mathrm{e}}_{\mathrm{k}}^{\mathrm{v}}\right)\right)$$

(10)

where ${\mathrm{S}}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}$ measures the similarity between the soil at pixel (i,j) and the representative known soil sample point k, ${\mathrm{e}}_{\mathrm{i}\mathrm{j}}^{\mathrm{v}}$ is the value of the v-th environmental variable at location (i,j) and ${\mathrm{e}}_{\mathrm{k}}^{\mathrm{v}}$ is the value of the v-th environmental variable associated with sample k, ${\mathrm{E}}^{\mathrm{v}}$ is the function that evaluates the similarity on the individual environmental variable level, The P-function is used to integrate the variable-level similarities into the location level. n is the total number of typical sample points within the region, ${\mathrm{S}}_{\mathrm{i},\mathrm{j}}$ represents the estimated soil salinity value at the unknown pixel (i,j) within the region, ${\mathrm{S}}_{\mathrm{k}}$ represents the soil salinity value at the known typical sample point k, ${\mathrm{w}}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}$ represents the weight assigned to k by the unsampled point [74].

2.3.2. Precision Validation

In order to verify the accuracy of the model, we selected three indicators to evaluate the accuracy of the inference results: the coefficient of determination (R²), which indicates the correlation between the observed and simulated data, with higher R² values indicating higher simulation accuracy, and the root mean square error (RMSE) and the mean absolute error (MAE), which represent the errors between the observed and predicted values and indicate the stability of the model, with smaller values indicating higher accuracy and stability. The formulas used to calculate these three indicators are as follows:

$${\mathrm{R}}^2=1-\frac{{\sum }_{\mathrm{i}}{\left(\widehat{{\mathrm{y}}_{\mathrm{i}}}-{\mathrm{y}}_{\mathrm{i}}\right)}^2}{{{\sum }_{\mathrm{i}}\left(\overline{{\mathrm{y}}_{\mathrm{i}}}-{\mathrm{y}}_{\mathrm{i}}\right)}^2}$$

(11)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{\mathrm{n}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\widehat{{\mathrm{y}}_{\mathrm{i}}}\right)}^2}$$

(12)

$$\mathrm{M}\mathrm{A}\mathrm{E}=\frac{1}{\mathrm{n}}\sum _{\mathrm{i}=1}^{\mathrm{n}}\left|\left({\mathrm{y}}_{\mathrm{i}}-\widehat{{\mathrm{y}}_{\mathrm{i}}}\right)\right|$$

(13)

where y_i represents the measured soil salinity value, ȳ represents the mean of the measured soil salinity values, ŷ represents the predicted soil salinity value, n represents the number of sample points.

3. Results

3.1. Statistical Analysis of Soil Salinity at Different Depths

To compare soil salinity at different depths, we performed statistical analysis on the sampled data (

Table 2. Statistical characteristics of soil salinity content (g·kg⁻¹) at different depths.

Year	Depth (cm)	Max	Min	Mean	SD	CV	n
2009	0–30	15.25	0.95	4.40	3.18	0.72	71
	30–60	12.5	1.2	4.08	2.95	0.72	71
	60–90	12.05	1.17	3.83	2.40	0.63	71
	0–90	35.3	3.37	12.31	7.43	0.60	71
2017	0–30	9.67	0.68	2.33	1.91	0.82	51
	30–60	19.33	0.57	4.87	4.85	0.99	51
	60–90	20.15	0.68	5.33	5.17	0.97	51
	0–90	41.38	2.47	12.44	10.54	0.85	51

SD represents standard deviation; CV represents coefficient of variation; n represents sample size.

). In the spring of 2009, the mean soil salinity content (SSC) in the surface soil (0–30 cm) was higher than in the middle and lower layers, ranging from 0.95 to 15.25 g·kg⁻¹, with an average of 4.40 g·kg⁻¹. In the summer of 2017, the mean SSC in the lower soil layer (60–90 cm) was the highest, ranging from 0.68 to 20.15 g·kg⁻¹, with an average of 5.33 g·kg⁻¹. Soil salinity increased with soil depth. In May 2009, salt accumulation was mainly in the surface soil layer, while in July 2019, salt accumulation was observed in the lower soil layer.

The coefficient of variation (CV) reflects the degree of variability in the samples. In 2009, the CV values for soil salinity in different layers of the study area ranged from 0.60 to 0.72, with significant variability observed in the surface and middle layers (0–60 cm). In 2017, the CV values for soil salinity in each layer ranged from 0.82 to 0.99, with the highest variability observed in the middle layer (30–60 cm).

3.2. Accuracy Evaluation of the Five Models

Figure 4a–c illustrate the simulation accuracy of five models for inferring SSC in the three soil layers during the spring of 2009. The SoLIM outperforms the two other types of models in simulating SSC, with R² values for the SoLIM in the three layers of 0.70, 0.65, and 0.59, respectively. In addition, the RMSE and MAE for the SoLIM are smaller than those for the other four models (

Table 3. Simulation error for the five models.

Year	Depth (cm)	Indicators	SoLIM	OK	IDW	MLR	GWR
2009	0–30	RMSE	1.95	2.74	2.34	3.16	3.04
		MAE	1.61	2.26	1.93	2.51	2.36
	30–60	RMSE	2.05	2.56	2.35	2.78	2.74
		MAE	1.62	2.04	1.84	2.22	2.08
	60–90	RMSE	2.10	2.41	2.47	2.32	2.33
		MAE	1.51	1.97	1.99	1.75	1.78
2017	0–30	RMSE	0.73	1.80	1.51	1.37	1.28
		MAE	0.55	1.01	0.98	1.12	0.91
	30–60	RMSE	3.00	4.00	3.57	4.50	2.85
		MAE	2.17	2.98	2.68	4.60	2.46
	60–90	RMSE	4.32	5.88	6.10	4.77	5.16
		MAE	3.00	4.78	4.58	3.84	3.86

RMSE is root mean square error, MAE is mean absolute error, and the unit is g·kg⁻¹.

). When comparing regression models and spatial autocorrelation models, it was observed that the performance of regression models is superior to spatial autocorrelation models. However, both types of methods underestimate the predicted values to some extent and failed to discriminate areas of severe salinization.

Figure 4d–f show the comparison between the simulated data from the five models and the measured data in July 2017. The SoLIM shows the highest simulation accuracy for the surface soil salinity (SSC), with an R² of 0.83 and RMSE and MAE of 0.73 g·kg⁻¹ and 0.55 g·kg⁻¹, respectively. In the lower soil layer, the simulation accuracy is slightly lower, with an R² of 0.64. The simulation values of the regression model have an R² of around 0.5, and the spatial autocorrelation model performs the weakest, showing a significant underestimation of the SSC predictions. Overall, the SoLIM shows high simulation accuracy for SSC at different time periods, with superior performance in the surface layer compared to the middle and lower layers. The other two methods show larger simulation errors for high SSC values.

3.3. Spatial Analysis of Soil Salinity Inference with SoLIM

Since the SoLIM has higher inference accuracy, we used the simulation results of the SoLIM to analyze the spatial distribution characteristics of soil salinity. Figure 5a–c illustrate the spatial distribution of SSC in the MRB region in May 2009. High salinity values are mainly concentrated in the central part of the study area, in the 0–60 cm soil layer with salt accumulation, where SSC reaches 6–10 g·kg⁻¹. In the northern oasis–desert transition zone with lower elevation, SSC ranges from 3 to 6 g·kg⁻¹, while in the southern region with higher elevation and less salinization influence, SSC is below 3 g·kg⁻¹. Figure 5d–f show the spatial distribution of SSC in July 2017. It is clear that salinization is less severe in the 0–30 cm soil layer, with SSC below 3 g·kg⁻¹ in most areas. Salt accumulation is pronounced in the 30–90 cm soil layer, especially in the downstream region of the Manas River in the north, where SSC in the 30–60 cm soil layer reaches 20 g·kg⁻¹. Extensive salt accumulation is observed in the central cultivated area, where the SSC reaches 15 g·kg⁻¹, and the 60–90 cm soil layer shows a distribution pattern with higher values in the north and lower values in the south.

3.4. Distribution Characteristics of Salinity for Different Soil Types

We analyzed the variations in salt content under different soil types during two time periods. The proportions of the three main soil types, Gleysols, Arenosols, and Fluvisols, in the study area were 71%, 14%, and 13%, respectively (Figure S2). As shown in Figure 6a, the soil type with the greatest change in salt content is Gleysols. The area of slightly saline regions increased by 264.85 km², while the areas of moderately saline and strongly saline regions decreased by 123.01 km² and 150.93 km², respectively. Slightly saline areas in the other two soil types also increased. Overall, compared to 2009, salt accumulation in 2017 showed a decrease in areas with moderately saline or highly saline characteristics and an increase in areas with non-saline and slightly saline characteristics. The distribution of the difference in soil salt content between 2009 and 2017 in the study area is shown in Figure 6d, where the areas with decreased and increased salt content account for 55.13% and 44.87%, respectively. The areas with decreased salt content are widely distributed in the southern cropland and central grassland regions of the study area, while the most salt-affected areas are concentrated in the northern part of the study area, especially in the lower-altitude cropland downstream of the Manas River (Figure 6b).

4. Discussion

4.1. The Relationship between Soil Salinization in Arid Zones and Environmental Factors

The spatial distribution of soil salinity is influenced by both structural factors, such as climate, parent material, terrain, and soil type, as internal drivers [75,76], and stochastic factors, such as fertilizer application, management practices, and cropping systems [77], as external influences. The interaction of these factors determines the overall spatial pattern of soil salinity [78]. Structural factors can enhance the spatial correlation of soil salinity, while stochastic factors tend to weaken its spatial correlation towards homogenization [79]. Arid zones have a higher rate of evapotranspiration, and changes in salinity are strongly influenced by both climatic conditions and human activities [80]. During the spring season, the accelerated rate of evaporation due to rising temperatures causes salts to accumulate in the surface soil as they move with the water. However, during the crop growing season, the salt in the surface soil migrates downward with the irrigation water, resulting in salt accumulation in the subsoil. The simulation results of this study effectively capture this characteristic during both seasons. In recent years, local agricultural management authorities have taken various measures to deal with soil salinity, such as constructing drainage ditches and modifying irrigation systems. The overall soil salinity in the study area decreased in 2017 compared to 2009, indicating the effectiveness of these mitigation measures. Therefore, in the future, fine mapping of soil salinity should be incorporated into agricultural management practices to help control soil salinization.

4.2. Comparison of Different Methods for Mapping Soil Salinity

The OK and IDW are two common soil mapping methods. Kulmatov, et al. [81] used the IDW method to produce soil salinity maps of the Sirdarya province, and their study showed that the IDW method has great potential for mapping soil salinity in irrigated areas. However, these two methods only focus on the spatial location of existing soil attribute samples and their own data [82], without considering the influence of environmental variables. Therefore, traditional spatial interpolation methods are more suitable for mapping soil attributes with a large number of evenly sampled points and spatial continuity [83]. Yuan, et al. [84] compared the applicability of linear regression and GWR in mapping soil salinity in arid areas and found that the GWR model, which considers surface conditions, performed better. Li, et al. [85] compared the performance of MLR, GWR, and random forest (RF) models in estimating soil salinity in a semi-arid oasis. The results showed that RF outperformed GWR, which in turn outperformed MLR. However, the simulation accuracy of machine learning and regression methods heavily relies on the number of samples available. In situations with a limited number of samples, these models tend to perform poorly [19,86]. On the other hand, the SoLIM is capable of acquiring soil–environment relationship knowledge and inferring soil properties at the watershed or even larger scales [24,87,88], while ensuring a certain level of inference accuracy with a limited number of samples [89,90,91]. Our study demonstrates that the SoLIM exhibits good applicability in predicting soil salinity in arid regions. Therefore, we used the simulated results of the SoLIM to analyze the spatial distribution of soil salinity in the study area. However, there still exists some deviation between the model-simulated spatial distribution and the observed soil salinity, which may introduce a certain level of uncertainty to our results.

5. Conclusions

This study used the Soil–Land Inference Model (SoLIM) to infer soil salinity content (SSC) of oasis croplands in the Manas River Basin, a region in arid northwestern China. The spatial distribution characteristics of SSC were analyzed for two periods: May 2009 and July 2017. Based on the inference results, the SoLIM showed higher accuracy in each soil layer compared to spatial autocorrelation models and regression models, indicating that the SoLIM can achieve acceptable accuracy with a limited number of sampling points. In the spring of 2009, soil salinity was concentrated in the 0–30 cm layer, while in July 2017, salinity accumulated in the 60–90 cm layer. The degree of soil salinity in the study area gradually increased with decreasing elevation. In terms of total salt accumulation, soil salinity was lower in 2017 than in 2009. Overall, the SoLIM demonstrated its ability to capture the soil environmental characteristics of oasis areas in arid regions based on sparse samples, demonstrating its good potential for inferring the soil salinity of croplands in arid regions.