Comparison of Different Interpolation Methods for Prediction of Soil Salinity in Arid Irrigation Region in Northern China

: Numerous methods have been used in the spatial prediction of soil salinity. However, the most suitable method is still unknown in arid irrigation regions. In this paper, 78 locations were sampled in salt-affected land caused by irrigation in an arid area in northern China. The geostatistical characteristics of the soil pH, Sodium Adsorption Ratio (SAR), Total Salt Content (TSC), and Soil Organic Matter (SOM) of the surface (0–20 cm) and subsurface (20–40 cm) layers were analyzed. The abilities of the Inverse Distance Weighting (IDW), Ordinary Kriging (OK), and CoKriging (CK) interpolation methods were compared, and the Root Mean Square Error (RMSE) was used to justify the results of the methods. The results showed that the spatial distributions of the soil properties obtained using the different interpolation methods were similar. However, the surface layer exhibits more spatial heterogeneity than the subsurface layer. Based on the RSME, the nugget/sill value and range significantly affected which method was the most suitable. Lower nugget/sill values and lower ranges can be fitted using the IDW method, but higher nugget/sill values and higher ranges can be fitted using the OK method. These results provide a valuable reference for the prediction of soil salinity. method. This may be determined by both the nugget/sill values and value ranges of the predicted properties. Lower nugget/sill values and smaller ranges may be fitted best using the IDW method, while higher nugget/sill values and larger ranges may be fitted best using the OK method. These results provide a powerful basis for the spatial prediction of soil salinity and provide useful information for irrigation, fertilization, and other land management measures.


Introduction
Soil salinization is one of the most serious land degradation types, limiting agricultural productivity worldwide [1,2]. Higher soil salinity always results in a high osmotic pressure between the plants and soil, which leads to water loss from the plants [3]. Most crops are sensitive to soil salinity [1]. The soil salinity is mainly affected by the following factors: (1) the application of water (mainly in agriculture) containing salts; (2) the weathering of the primary and secondary minerals in the soil; (3) the decay of organic matter in the soil; and (4) instability of the underground water table [4]. Soil salinization always influences the soil pH, Sodium Adsorption Ratio (SAR), Soil Organic Matter (SOM), and other properties [5]. For example, soil salinity has been found to be negatively correlated with pH (p < 0.001); that is, lower pH soils tend to be more saline than higher pH soils [6]. In addition to pH, higher soil salinity is always accompanied by higher sodium contents [7]. Sodium ions can cause instability of soil aggregates, followed by dispersion of clay particles, which results in the clogging of soil pores. The adsorption of sodium by the soil is expressed by the SAR of the solution [8], and a higher SAR always results in a higher crop growth limitation. Soil salinity also has a significant negative impact on the SOM. This is due to the decreased plant production and/or alteration of the community composition, which changes the quality and quantity of the SOM and increases the decomposition rate [9].
Knowledge of the spatial distribution of soil salinity helps obtain information for land reclamation or restoration [10], in making effective site-specific management decisions [3], and in ecological and agricultural management [11]. In addition, it can provide guidelines for irrigation and the application of gypsum and other chemicals [10]. For soils with high salinity, more infiltration water is needed to enable salt leaching. Moreover, the application of gypsum should be increased. The soil salinity may always have a considerable nugget effect [12]. In addition to soil salinity, soil pH may also have a significant spatial dependence. However, the soil pH variation has been found to be lower than the other soil properties, and it increases slightly with soil depth [3]. The SAR is calculated from the Na + , Mg 2+ , and Ca 2+ contents, and it is considered a more appropriate property for determining soil alkalization than a single cation concentration [13]. Knowing the spatial distribution of the SAR is also important for land reclamation and agricultural development [14]. Knowing the spatial distribution of the SOM is crucial for crop growth. Studies have shown that the SOM in saline soils tends to be low and has high spatial variability.
Geostatistical analysis can be applied in decision-making regarding environmental monitoring, remediation, land management, and planning [4]. Many interpolation methods have been used for the spatial prediction of soil salinity, including the Inverse Distance Weighting (IDW), Ordinary Kriging (OK), and CoKriging (CK) methods. The IDW assumes that the soil properties of a point are the weighted average of the measured values within the neighborhood, and the weights are inversely related to the distances between the prediction location and the sampled locations [15]. Ordinary kriging is a geostatistical method that is widely used for spatial interpolation. The OK method recognizes and models spatial variables using anisotropies, declustering, and unbiasedness, and it also computes the minimized estimation variance (also known as the kriging variance) [16]. Compared to the IDW and OK methods, the CK method considers the distance between the observed variable values and the predicted variable values and the spatial correlation between one variable and another more frequently sampled variable [17]. The accuracies of these methods have been compared in different areas and used in predicting different variables. For soil pH and electrical conductivity (EC) in coastal agricultural areas in northern Iran, the OK method had a smaller error than the IDW method for both soil properties in the different soil layers [18]. For pH, EC, and SOM in crop farmland located in the Shire of Wickepin, Australia, the OK method performed best for pH in the topsoil, and the lognormal ordinary kriging method performed best for EC in the topsoil. The IDW method had the greatest accuracy for subsoil pH [19]. For EC, pH, SAR, and other soil properties in a reclaimed area of the Behera Governorate of Egypt, the IDW method had a higher prediction efficiency than the kriging method [20]. For the soil EC in a representative rural community located in the western part of Biskra, Algeria, the CK method performed better than the OK method, with RMSE values of 0.92 and 1.53, respectively [21]. As can be seen, the most suitable method for the spatial prediction of soil salinity and the related indexes is still under debate, and further study is needed in specific regions.
The Hetao Irrigation District is a typical arid irrigation region, and it is very important to the agricultural development of Inner Mongolia, a region in northern China. It is one of China's three largest irrigation districts and is an important food supply base [22]. The precipitation in this area is much lower than the crop water requirement. Irrigation is necessary for agriculture in this area. More than 75% of the basin is irrigated using Yellow River water. The average water diversion from the Yellow River is about 5.165 billion m 3 /year (from 1980 to 2000), which accounts for about 1/10 of the annual discharge of the Yellow River [23,24]. Long-term irrigation has resulted in salinization, which is becoming increasingly serious [25]. Knowing the spatial distribution of the soil salinity is crucial to precise fertilization and irrigation. The method most commonly used to study the spatial distribution of soil salinity is remote sensing [22,26]. However, the remote sense method mainly focuses on large-scale areas. Understanding the spatiotemporal characteristics of soil salinity is necessary for both large regions and the field scale. The spatial distributions of the soil salinity and crop canopy coverage are consistent at the field scale [27]. Thus, the field scale prediction of soil salinity may be more helpful for farmers when planning fertilization and irrigation. A detailed understanding of the spatiotemporal variability of soil salinity is helpful in agriculture regionalization at the regional scale. For example, where to plant salt tolerance crops and plant high yield crops without salt tolerance.
This study aims to compare the abilities of the different spatial interpolation methods to predict the spatial distributions of the soil salinity and related soil properties and to analyze their spatial distributions in the arid irrigation region, northern China.

Site Description
The study area is located in the eastern part of the Hetao Irrigation District, northern China ( Figure 1). It has a typical arid climate, with low precipitation and very high evapotranspiration. The annual precipitation is about 130-350 mm, and the annual temperature is about 5.6-7.4 °C [22]. The soil type is solonchak, according to World Reference Base for Soil Resources. The soil texture is clay loam with sand, silt, and clay contents of 24.11%, 43.60%, and 32.29%. The underground water table is high due to the frequent irrigation conducted every year. The average amount of irrigation water derived from the Yellow River every year is about 30-40 cm. The total dissolved salt in the Yellow River is about 0.5%. Thus, a great deal of salt has been added to the irrigated fields every year through irrigation. However, the drainage system is low compared to the irrigation system, i.e., only 10%. Thus, high irrigation and low drainage is a serious problem in the study area. The study area is next to Ulansuhai Lake, a drainage lake in the Hetao Irrigation District. Therefore, the water table in the study area is low and only varies by 0.5-1.5 m throughout the year. This has resulted in significant salinization in the area.
. Figure 1. Locations of the study area and sampling points. The number beside the point is the sampling ID

Soil Sampling and Laboratory Measurements
From the study field site (size in ha), soil bag samples from the topsoil (0-20 cm) and the subsoil (20-40 cm) at 78 sampling points were taken (Figure 1), which adds up to 156 soil bag samples. Then, they were taken to the laboratory and were air-dried and ground in preparation for analysis. The soil properties, including the pH, Na + , Ca 2+ , Mg 2+ , total salt content (TSC), SOM, sand content, silt content, and clay content, were measured. The pH was measured using a pH meter (Sartorius PB10). The Na + content was measured using a flame photometer, and the Ca 2+ and Mg 2+ contents were measured via ethylene diamine tetraacetic complexometric acid titration [28]. The SOM was measured using the oxidation of the potassium dichromate external heating method [29]. The SAR was calculated using the Na + , Ca 2+ , and Mg 2+ contents and the following formula, in which the unit of the ion concentrations is mmol L −1 [13]. Soils with SARs of greater than 12, exchangeable sodium percentages (ESPs) of greater than 15, and EC values of greater than 5 pose a very high hazard to plant growth [7]. (1)

Semivariance Modeling
Semivariance modeling is widely used in geostatistical analysis. In soil science, the soil properties are always considered to be regionalized variables. The semivariance (γ(h)) of a considered soil property can be calculated as follows [30]: where N(h) is the pair of points with lag distance h; Z(i) is the value of the considered property at point xi, and Z(xi + h) is the value of the property at a point a distance h away from point xi. Soil is a continuous spatial variable, and the semivariogram of a soil property is a continuous function. However, the semivariance is always calculated from a series of soil samples that contain discontinuities. Thus, the semivariance of these points can be simulated using curves, and the equation used for the curve fitting is called the theoretical model of the semivariogram. The commonly used models are the spherical model, the exponential model, and the Gaussian model: where γ(h) is the semivariance; h is the lag distance; C0 is the nugget, and C is the sill. C0/(C0 + C) can reflect the degree of spatial dependence. The spatial dependence is divided into strong (C0/(C0 + C) < 25%), moderate (25% ≤ C0/(C0 + C ≤ 75%), and weak (C0/(C0 + C) > 75%) dependence [31]. The simulation is done in professional geostatistical software (GS + 9.0). The software simulated the data with all three models and compared the R 2 of them. The best model is chosen with the highest R 2 value.

Inverse Distance Weighting Interpolation
The IDW method is a widely used, straightforward, noncomputational interpolation method. The predicted value (Pj) for location j is calculated based on the distance weight [15].
Where Oi is the observed value at location i; and λi is the weight, which can be calculated as In this equation, dij is the distance between i and j; and α is a specified value, which is always larger than 1 and is always taken to be 2 [32]. As can be seen, points closer to the target location have greater weights.

Ordinary Kriging Interpolation
Ordinary kriging is a widely used interpolation method, which can be used to obtain an optimal, linear, and unbiased prediction for a nonsampled point based on the values of neighboring sample points. The OK method assigns weights to these points based on their distances from the point being estimated and the spatial variability structure, and the formula is as follows [21,33]: In this equation, Z(x0) is the predicted value; and wi is the weight for ith observation value Z(i). To obtain an optimal and unbiased prediction, the sum of the weights must be equal to 1, and the estimation or kriging variance is minimized.

Cokriging Interpolation
In the IDW and OK predictions, the only parameter used to determine the weight is the distance. Variables that have spatial relationships with two or more variables can be predicted using the cokriging method: In this equation, Zu(b) is the predicted variable u; and V is the number of related variables. N is the number of pairs of points. wiI is the weight of the Ith variable at point i, and ZI(xi) is the observed matrix value for the Ith variable at point i. To obtain an optimal and unbiased prediction, the sum of the weights of the predicted properties must be 1, and that of the other properties must be 0 [33]: In addition, the estimation or the cokriging variance should be minimized. In this study, soil texture (sand content, silt content, or clay content) was selected as the covariance in the cokriging prediction. This is mainly because the soil texture forms the basis of the soil structure and strongly affects the leaching of salt [14].

Validation of the Different Interpolation Methods
A total of 78 points were sampled in the study area. Thirty of them were selected as validation data, and the other 48 points were used as interpolation data. The IDW, ordinary kriging, and cokriging methods were used to predict the TSC, pH, LogSAR, and SOM values for both the surface and subsurface layers in ArcMap 10.6. Then, the prediction maps were converted to raster style, and the predicted values of the validation points were extracted using the "Extract value to points" tool in ArcMap. Thereafter, the Root Mean Squared Error (RMSE) of the predicted value for each soil property was calculated as where m is the total number of validation points, yi is the observed value, and yj is the predicted value. For a given property, the lower the RMSE, the more accurate the method.
The regression analysis of the predicted and observed values were also conducted in SPSS 22.0. A higher coefficient of determination (R 2 ) means more variance was interpreted using the predicted method.

Statistical Analysis and Data Distribution of the Studied Soil Properties
The mean TSC values are 16.4 g/kg and 13.42 g/kg for the surface and subsurface soil layers, respectively ( Table 1). The TSC of the surface soil layer is significantly higher than that of the subsurface soil layer. However, the coefficients of variation (CVs) for the soil layers are similar. The TSC is significantly higher in the surface layer than in the subsurface layer. The SAR and SOM are also significantly higher in the surface layer than in the subsurface layer. The CVs are also similar for the two layers, both of which indicate moderate variation. Unlike these properties, pH has both similar mean values and CVs for the two soil layers. Moreover, the CVs indicate weak variation. The surface layer's sand, silt, and clay contents were 22.99%, 45.58%, and 31.43%. However, they were 25.27%, 41.02%, and 33.70% for the subsurface layer. Although the sand, silt, and clay contents all exhibit moderate variation, the variations in the silt contents of both layers are close to being weak. The data for the soil properties exhibit normal distributions, except for SAR for both layers and the sand content of the subsurface layer, which obey normal distributions after log transformation (Table 2). CV-coefficients of variation, SD-standard deviation, SAR-sodium adsorption ratio, TSC-total salt concentration, SOMsoil organic matter. The letter after the mean value means the difference is significant at 0.05 level.

Geostatistical Analysis of the Soil Salinity and Related Soil Properties
The considered soil properties exhibit a clear nugget effect, except for the pH of the surface layer, SAR of the surface layer, and the SOM in both layers (Figure 2). The exponential model best fits the TSC for the surface layer and the subsurface layer's spherical model. The nugget/sill values for the surface and subsurface soil layers are 11.76% and 36.12%, respectively (Table 3). This suggests a strong spatial dependence for the surface layer but a moderate spatial dependence for the subsurface layer. The range for the surface layer (170.4 m) is much smaller than that for the subsurface layer (291.7 m). R 2 of the subsurface layer (0.81) is higher than the surface layer (0.62). For the other soil properties, the values for the surface layer are best fitted by the spherical model and those of the subsurface layer by the exponential model.
Moreover, the nugget/sill values for the pH and SAR for the surface layer are much higher (moderate spatial dependence) than those for the subsurface layer (strong spatial dependence). However, for the SOM, both layers exhibit a strong spatial dependence. The ranges of the pH, SAR, and SOM values for the surface layer are much higher than those for the subsurface layer, indicating a longer spatial correlation for the subsurface layer.   LogSAR means the SAR has been transferred to a normal distribution using a log transformation, C0-nugget value, and C-sill value.

Correlations between Soil Salinity and Other Soil Properties
In the surface layer, the TSC was significantly negatively correlated with the pH and the clay content, but it was very significantly positively correlated with the silt content at the p < 0.01 level ( Table 4). The TSC was significantly correlated with LogSAR at the p < 0.05 level. No significant correlations were found between the TSC and the SOM and sand content. The subsurface layer exhibits very similar correlations between the TSC and other soil properties. However, the correlation coefficients were higher than those for the surface layer, especially for the pH and LogSAR. The other soil properties exhibited significant pair correlations, especially with soil texture. The silt and clay contents were significantly correlated with the TSC and pH in both the surface and subsurface layers. The SOM was significantly correlated with the sand content in the surface layer and correlated with the sand content and clay content in the subsurface layer. Although the soil texture was not significantly correlated with LogSAR, the correlation between the clay content and LogSAR for the surface layer was significant at the p < 0.1 level. Therefore, the soil texture, which exhibits significant correlations with the considered soil properties, was input as a covariable in the CoKriging prediction. * means the correlation is significant at the 0.5 level, and ** means the correlation is significant at the 0.01 level.

Spatial Distributions of The Soil Salinity and Other Soil Properties Predicted Using the Different Methods
The general trends of the spatial distributions of the soil properties obtained using the different methods were similar. However, the IDW method showed a clearer patch than the other methods (Figures 3-6). The cokriging method was the most homogeneous compared to the other two interpolation methods. For the TSC of the surface layer, all methods produced higher value zones in the north and south and a lower value zone in the middle of the study area ( Figure 3). The higher and lower value zones are clearer for the subsurface layer than for the surface layer. The higher value zone is mainly located in the central-north and southwest areas, but the lower value zone is located in the southeastern part of the study area. The pH values exhibited a different spatial distribution than the TSC, mainly decreasing from southeast to northwest (Figure 4).
Moreover, the two soil layers exhibit similar spatial distributions. The spatial distribution of the SAR firstly increases and then decreases from west to east. The highest value zones appear in the middle part for the surface layer and in the southeast part for the subsurface layer ( Figure 5). The SOM showed a similar trend with the pH value. However, SOM exhibits different spatial distributions in the two layers ( Figure 6). The higher value zone is mainly located in the northeast for the surface layer, but in the subsurface layer, the higher value zone is mainly located in the southwest.

Comparison of the Different Prediction Methods
For the TSC of the surface layer, the cokriging method performed the best, but the kriging method performed the best (Table 5). For the pH and SAR of the surface layer and the SOM in both layers, the IDW method performed the best. For the pH and SAR of the subsurface layer, the kriging and cokriging performed the best, respectively. The scatter plot shows that the TSC for both layers, the pH for the subsurface layer, and the SOM for the surface layer exhibited a more decentralized trend (Figure 7), indicating the larger uncertainty of the prediction.

Soils Salinity of the Different Soil Layers and Its Correlation with the Other Properties
The soil salinity is usually higher for the surface soil layer than for the subsurface layer, especially for secondary salinization (or man-made soil salinization). Our study area's TSC of the surface layer is 16.40 g/kg, significantly higher than that of the subsurface layer (13.42 g/kg) ( Table 1). The same phenomenon was observed in a study in a reclaimed area of the Behera Governorate of Egypt. The concentrations of the main ions were higher in the surface layer than in the subsurface layer [20]. Still, the studied soil depths were different than those used in our study. In the Igdir Plain (Turkey), researchers have also found a decrease in EC with depth in the soil profile [3]. Another study found that the EC was higher in the surface layer than in the deeper layers, and the deeper layers did not exhibit significant differences [34]. In Xinjiang Province, northern China, the EC values in the surface soil layer (0-10 cm) were found to be significantly higher than at the 10-20 cm and 20-30 cm soil layers [35]. The general trend indicated salt aggregation in the surface soil. This was mainly caused by shallow underground saline water, which moved up through the soil due to capillary action and evaporation [10]. When the saline water reached the soil surface, the water evaporated into the air, and the salt was left behind in the soil. The soil salinity exhibited a moderate variation in our study. The soil EC was also found to exhibit a moderate variation (CV of EC was 47.4%) in the southeastern part of Algeria [21] and in the Manas River watershed in northern China (CV of EC was 33.1-63.0%) [10], but their CV values were higher than those in our study (Table 1). Another study reported a strong variation (CVs of 1.31-1.57) in the different soil layers [35]. The CVs of the different EC datasets may also differ. For example, in the coastal agricultural areas of northern Iran, the CV of the 0-50 cm layer was 56.7%, indicating moderate variation, but it was 114.7% for the 50-100 cm layer (strong variation). In our study, the variation in the two layers was the same. This may be caused by the higher sampling thickness (20 cm). Therefore, a thinner sample should be conducted in future studies.
Soil salinity can exhibit a positive [21] or negative [36] correlation with pH, but the correlation is generally not significant. In our study, a significant negative correlation was found between TSC and pH in both soil layers. This was mainly caused by the fact that the main cation in the soil in our study area is sodium. When the TSC decreases, the sodium ion concentration mainly decreased. This caused the proportion of divalent cations (Ca 2+ and Mg 2+ ) to increase. In addition, the proportion of bicarbonate in the anions increased, while that of sulfate decreased. This caused an increase in the soil pH [37]. The soil salinity is positively correlated with the SAR in our study area (Table 3), and similar results have been found in other studies [6,7,10]. This indicated that soil salinity might increase soil alkalinity. A previous study found that saline soils are characterized by relatively high SAR values (7.31-52.10), which are significantly higher than those of non-saline soils (2.66-9.53) [5]. The SAR values in our study were much higher than the values in their study, indicating more severe soil alkalinity. The soil texture had a significant influence on the soil salinity, and it may predominantly control the soil salinity distribution. The silt content had a significant positive correlation with soil salinity, but the sand content had a significant negative correlation with soil salinity [38]. Another study reported a significant positive correlation between soil salt content and sand content, but it was negatively correlated with the clay content [14]. We found that the TSC is significantly positively correlated with the silt content, but it is significantly negatively correlated with the clay content (Table 4). This may be because the capillary force is higher in silt, which may increase the accumulation of salt on the surface. However, the capillary force is lower in clay soils, so salt accumulation on the surface may be low.

Geostatistical Characteristic of Soil Salinity in the Arid Irrigation Region in Northern China
The semivariance of the salinity of the surface soil is fitted best by the exponential model, but the subsurface soil is fitted best by the spherical model (Table 3). A previous study reported different results and showed that the spherical model was best for the surface soil layer and the exponential model was best for the subsurface soil layer [3]. Studies have also reported that all of the soil profile layers are fitted best by the spherical model [4,18]. These different results may be caused by different sampling times and sampling strategies. The nugget effect increases with soil depth, although both of the layers in our study exhibit a moderate spatial dependence. One study reported that the surface layer has a strong spatial dependence, but the bottom layer has a moderate spatial dependence [18]. Other results indicate that, in general, the nugget/sill value decreases with soil depth [39]. The spatial dependence of soil salinity is different in different regions and at different times, affecting the spatial distribution predicted using the kriging method. Except for the TSC, all of the other soil properties showed a lower nugget effect for the surface soil layer than for the subsurface soil layer. In addition, the nugget/silt value was much lower for the surface layer. This indicates that the spatial dependence of the soil properties was higher in the surface soil layer.

The Best Method for Predicting Soil Salinity in the Arid Irrigation Distinct
Many researchers have evaluated the performances of the different interpolation methods. The results almost always differ. In a newly reclaimed area in the southern part of the Behera Governorate of Egypt, the IDW method outperformed the other methods, including the ordinary kriging, regression kriging, and cokriging methods [20]. In the agricultural coastal areas of northern Iran, the OK method had the minimum error compared with the IDW and conditional simulations methods [18]. In the southern part of the Arkansas River Basin in Colorado, USA, the performances of the different geostatistical models were as follows: ordinary kriging > regression kriging > cokriging [40]. In other words, the applicability of the spatial prediction methods to a specific soil in a specific region are still unknown. In our study area, different methods were found to be optimal for the different soil properties and soil layers (Table 5). We found that the optimal methods are closely related to the relationship between the nugget/sill value and the range of the considered soil properties values (Figure 8). The soil properties with lower nugget/sill values and smaller ranges are best predicted using the IDW method. The soil properties with intermediate nugget/sill values and value ranges are best predicted using the CK method. The soil properties with higher nugget/sill values and larger ranges are best predicted using the OK method. This may be caused by the different calculation methods used by the different prediction methods. The kriging method is a function of the statistical stationarity, while the IDW is a function of the distance between the sampling points [32]. When both the nugget/sill value and the range are lower, the property has a strong spatial dependence in a short distance. This condition may be more in accord with the assumptions of the IDW method. Figure 8. The relationship between range and nugget/sill value for different methods. Four points are in the IDW circle, including pH for the surface layer, SAR for the subsurface layer, and SOM for both layers. One point is in the CK circle, including TSC for the surface layer. The other three points are in the OK circle, including TSC for the subsurface layer, pH for the subsurface layer, and SAR for the subsurface layer.

Conclusions
The average TSCs were 16.40 g/kg and 13.42 g/kg for the surface (0-20 cm) and subsurface (20-40 cm) soil layers in a typical piece of salt-affected land in the arid irrigation region, northern China. In addition, the SAR was also extremely high, i.e., 55.64 for the surface soil and 49.79 for the subsurface soil layers. This suggests that both soil salinization and alkalization are serious problems in this area. Compared to the other soil properties, the TSC and SAR have higher CV values, indicating a higher variance in the soil salinity and alkalinity. The nugget/sill values of the TSC are 11.76% and 36.12% for the surface and subsurface soil layers, respectively. This indicates that the soil salinity has a moderate spatial dependence. The TSC is significantly correlated with the soil texture (mainly the silt and clay contents), but the LogSAR is not correlated with the soil texture. This indicates that the soil texture may more easily influence the soil salinity. The general trends of the prediction maps of the considered soil properties obtained using the three different methods (IDW, OK, and CK) are similar. However, the OK and CK methods produced more homogeneous results than the IDW method. The RMSEs of the different methods showed that the TSC was predicted best using the cokriging (surface) and kriging (subsurface) methods. The pH and SAR were predicted best using the IDW (surface) and kriging (subsurface) methods, and the SOM was predicted best using the IDW (both layers) method. This may be determined by both the nugget/sill values and value ranges of the predicted properties. Lower nugget/sill values and smaller ranges may be fitted best using the IDW method, while higher nugget/sill values and larger ranges may be fitted best using the OK method. These results provide a powerful basis for the spatial prediction of soil salinity and provide useful information for irrigation, fertilization, and other land management measures.