Estimating salt content of vegetated soil at different depths with Sentinel-2 data

The accurate and timely monitoring of the soil salt content (SSC) at different depths is the prerequisite for the solution to salinization in the arid and semiarid areas. Sentinel-2 has demonstrated significant superiority in SSC inversion for its higher temporal, spatial and spectral resolution, but previous research on SSC inversion with Sentinel-2 mainly focused on the unvegetated surface soil. Based on Sentinel-2 data, this study aimed to build four machine learning models at five depths (0∼20 cm, 20∼40 cm, 40∼60 cm, 0∼40 cm, and 0∼60 cm) in the vegetated area, and evaluate the sensitivity of Sentinel-2 to SSC at different depths and the inversion capability of the models. Firstly, 117 soil samples were collected from Jiefangzha Irrigation Area (JIA) in Hetao Irrigation District (HID), Inner Mongolia, China during August, 2019. Then a set of independent variables (IVs, including 12 bands and 32 spectral indices) were obtained based on the Sentinel-2 data (released by the European Space Agency), and the full subset selection was used to select the optimal combination of IVs at five depths. Finally, four machine learning algorithms, back propagation neural network (BPNN), support vector machine (SVM), extreme learning machine (ELM) and random forest (RF), were used to build inversion models at each depth. The model performance was assessed using adjusted coefficient of determination (R2adj), root mean square error (RMSE) and mean absolute error (MAE). The results indicated that 20∼40 cm was the optimal depth for SSC inversion. All the models at this depth demonstrated a good fitting (R2adj≈ 0.6) and a good control of the inversion errors (RMSE < 0.16%, MAE < 0.12%). At the depths of 40∼60 cm and 0∼20 cm the inversion performance showed a slight and a great decrease respectively. The sensitivity of Sentinel-2 to SSC at different depths was as follows: 20∼40 cm > 40∼60 cm > 0∼40 cm > 0∼60 cm > 0∼20 cm. All four machine learning models demonstrated good inversion performance (R2adj > 0.46). RF was the best model with high fitting and inversion accuracy. Its R2adj at five depths were between 0.5 to 0.68. The SSC inversion capabilities of all the four models were as follows: RF model > ELM model > SVM model > BPNN model. This study can provide a reference for soil salinization monitoring in large vegetated area.


INTRODUCTION
Soil salinization has been an important factor leading to crop yield reduction and land degradation in arid and semiarid areas (Harti et al., 2016). Efficient and accurate monitoring of soil salt content (SSC) on a large scale is the key to tackle this problem. Among the varied monitoring methods, satellite remote sensing has become increasingly prevailing.
So far, different satellites and methods have been applied to soil salinization monitoring. Lobell et al. (2010) first used MODIS data for regional-scale soil salinity assessment and reduced the effect of temporally dynamic factors using the mean of the enhanced vegetation index (EVI) for a 7-year period. With consideration of the effects of precipitation, crop type, and soil texture, Scudiero, Skaggs & Corwin (2014) assessed the SSC based on the average of multi-year Landsat 7 data, and obtained reliable results. Wu et al. (2014b) mapped soil salinity mainly with Landsat ETM+ and MODIS multi-year data, and achieved reliable salinity prediction results in vegetated and non-vegetated areas, respectively. IKONOS data were used to analyze the Pearson correlation coefficient between broadband indices and soil salinity, the results indicated that the correlation depended on the environmental conditions (soil, vegetation cover and density), and vegetation indices performed better in densely vegetated areas (Allbed, Kumar & Aldakheel, 2014). Landsat 8 data were used to construct 12 VI-SI (vegetation indices-salinity indices) feature spaces based on the information of bare soil and vegetation. Results showed that MSAVI-SI 1 (modified soil adjust vegetation index-salinity index) can greatly improve the dynamic and periodical monitoring of soil salinity . These studies on the relationship between multiple satellite data and soil salinity have provided a good basis for regional SSC assessment. However, each of the above satellites has demonstrated such defects as low spatial resolution or small spectral range. Sentinel-2 has shown certain advantages because it simultaneously has high temporal and spatial resolution, which enable more detailed and higher-frequency monitoring for practical applications. Additionally, Sentinel-2 can obtain the red-edge region of vegetation spectrum, which can provide more effective data for vegetation growth monitoring.
Scholars have conducted some research on SSC inversion with Sentinel-2 data. Wang et al. (2020) estimated soil salinity using the machine learning model, Cubist. By comparing the two SSC distribution maps (at the depth of 0∼20 cm), they found that Sentinel-2 outperformed Landsat 8 in accuracy. Davis, Wang & Dow (2019) and Gorji et al. (2020) also discovered that Sentinel-2 had great potential for SSC inversion. Taghadosi, Hasanlou & Eftekhari (2019) established two models (multiple linear regression and support vector regression) using Sentinel-2 images, which had good performance in SSC inversion in the unvegetated areas. Wang et al. (2019) created multiple spectral indices based on Sentinel-2 data and developed an RF-PLSR model to estimate SSC. The above studies on SSC inversion with Sentinel-2 data were mostly concentrated in the surface soil. Ramos et al. (2020) evaluated soil salinity at the depth of 0∼1.5m via multiple stepwise regression based on multi-year Sentinel-2 data, and obtained relatively high prediction accuracy (the coefficient of determination ranged from 0.63 to 0.91). This study lays a groundwork for soil salinity estimation at root depth based on Sentinel-2 data. However, the evaluation of SSC at different root depths in the vegetated soil remains to be investigated.
In addition, the machine learning algorithms have been widely used for SSC inversion and water resources management . Chen et al. (2015) studied the accuracy of multiple linear regression (MLR), back propagation neural network (BPNN), and support vector machine (SVM) in soil salinity estimation using hyperspectral data, and found that BPNN and SVM were more accurate than MLR. When studying the hybrid particle swarm optimization with extreme learning machine (ELM) for daily reference evapotranspiration (ET 0 ) prediction from limited climatic data, Zhu et al. (2020) explored the ability of artificial neural networks (ANN), random forests (RF) and other empirical algorithms in estimating daily ET 0 . The results indicated that the machine learning models outperformed the corresponding empirical algorithms. Machine learning algorithms are available for SSC inversion, and yet the accuracy of each algorithm using Sentinel-2 data to estimate SSC needs more in-depth comparison.
SSC evaluation at different root depths in the vegetated soil and the algorithm accuracy in SSC estimation via Sentinel-2 data both demands further research. Therefore, this study used the Sentinel-2 images of the Jiefangzha Irrigation Area (JIA) in the vegetated area to construct the set of independent variables (IVs, including 12 bands and 32 spectral indices). Next, the optimal combinations of IVs at five depths (0∼20 cm, 20∼40 cm, 40∼60 cm, 0∼40 cm, and 0∼60 cm) were obtained using full subset selection. Finally, four machine learning algorithms, back propagation neural network (BPNN), support vector machine (SVM), extreme learning machine (ELM) and random forest (RF), were used to construct inversion models (models for SSC estimation via satellite data) and evaluate the sensitivity of Sentinel-2 to SSC at different depths and the inversion capability of the models.

Study area
The study area, JIA, is located in the northwest of Hetao Irrigation District (HID), Inner Mongolia, China, between 106 • 34 ∼107 • 34 E, and 40 • 26 ∼41 • 13 N, which is the same as that of Qiu et al. (2019). JIA, an oblique triangular area about 86 km long and 81 km wide, is the second largest irrigation area in HID. With an altitude of about 1,030 m to 1,046 m, this plateau is high in the southwest, low in the northeast, and relatively flat on the whole. It is located in an arid and semiarid area with a temperate continental climate. The annual mean temperature, annual precipitation and evaporation are 4∼6 mm, 66.3∼200 mm and 1,920∼3,450 mm, respectively. Moreover, the annual precipitation distribution is uneven (the precipitation in summer accounts for about 70% of the year). JIA is 2156.7 km 2 in total, about 57.5% of which is irrigatable. The crop planting structure is complex, and the crops mainly include corn, wheat, sunflower and so on. The local climate and hydrological conditions determine the dependence of the crop growth in this area mainly on the irrigation from the Yellow River. The annual water diversion in JIA is about 1.2 billion m 3 , and the land is mainly irrigated in summer and autumn. The water diversion in October accounts for about 30% of the annual amount. Figure 1A shows the specific geographical location of the study area.

Soil sample collection and analysis
The Hetao irrigation district administration gave field permit approval to us (No. 2017YFC0403302). Among the main crops in JIA, sunflower, corn and zucchini are usually harvested around September 20. Therefore, August 25∼30, 2019, when the crops were flourishing and the vegetation information was abundant, was chosen as the study period. We collected 117 soil samples (83 covered with sunflowers, 10 cabbages, 10 corns, 2 wheats, 2 vegetables, and 10 bare soil) (Fig. 1B) when the types of underlying surface, salinization degree, and evenness of point distribution were taken into consideration. According to the root depth of the crops and related research results (Qiao, 2005;Zhang et al., 2019b), we selected 0∼20 cm 20∼40 cm and 40∼60 cm as the sampling depth. We adopted the five-point sampling method so that each of the 117 soil samples was the mixed soil from one center and four corner points in a cell (0.5 m × 0.5 m). The GPS data and environment information of the center point were recorded during the sampling.

Measurement of SSC
The soil samples were first dried, ground, and then screened with a 2.0-mm sieve to remove the small stones and wood pieces. The processed samples were mixed with water to make soil solution (the ratio of soil to water is 1:5). Eight hours after the solution was prepared, its electrical conductivity (EC 1:5 , ds/cm) was measured by conductivity meter (DDS-307A; Shanghai Youke Instrument Branch, Shanghai, China), and then the SSC (%) was calculated by the empirical formula (1) , which was developed by Huang et al. (2018) in their research in HID.

Statistical characteristics of SSC
We selected five depths (0∼20 cm, 20∼40 cm, 40∼60 cm, 0∼40 cm and 0∼60 cm) in this study, and the SSC of 0∼40 cm was the mean value of 0∼20 cm and 20∼40 cm and that of 0∼60 cm was the mean value of 0∼20 cm, 20∼40 cm and 40∼60 cm. The sample points were sequenced according to the SSC, and then one of every three samples were selected as the validation dataset so that the ranges of calibration and validation datasets were consistent and evenly distributed. The statistical characteristics of SSC are shown in Table 1.

Sentinel-2 image data
Sentinel-2 consists of two satellites (Sentinel-2A and 2B), and both provide multi-scale (10m, 20m and 60m) remote sensing images via MultiSpectral Instrument (MSI). It has a 5-day revisiting period when the two satellites are combined. It has 13 bands (440 nm∼2200 nm), including visible light, near infrared and short-wave infrared bands. Three of the bands can obtain the red-edge region of vegetation spectrum (670 nm∼760 nm), which can provide more effective data for vegetation growth monitoring. The remote sensing image data for this study were taken from Sentinel-2A, and its parameters of each band are shown in Table 2.
The satellite images in this study were obtained from the website (https://scihub. copernicus.eu/) of the European Space Agency (ESA), on August 28, 2019, which was basically synchronous with the sampling time and there was no cloud in the study area. The Sentinel-2 data used in this study were Level-2A products (Bottom-Of-Atmosphere reflectance images), which were produced by the plug-in Sen2cor (Level-2A product can also be obtained directly from ESA since December 2018). Then all the bands were resampled to 10m via the S2 Resampling Processor in the software SNAP (the bicubic method was used for resampling). As band B10 was not available when the images were processed to Level-2A, the other 12 bands were used in this study.

Selection of spectral index
In this study, we selected 32 widely used spectral indices, including salinity index, vegetation index and drought index. The indices and the relevant formulae are shown in Table 3.

Soil line fitting
As has been shown in the Nir-Red scatterplot of several studies, when the horizontal and vertical coordinates are the red and Nir bands, respectively, a series of corresponding points of the digital number values of the red and Nir infrared wavelengths of the bare soil Table 3 Spectral indices. M is the slope of soil line, I is the intercept of soil line, A is the PVI maximum point. The fitting result of soil line is shown in Fig. 2.

Spectral index formula Reference
Salinity Index (SI)

Scudiero, Skaggs & Corwin (2015)
Intensity index 1 (Int1) Fourati et al. (2015) Normalized Vegetation Index (NDVI) Normalized Difference Vegetation Index red-edge 2 (NDVIre2) Khaled (2017) approximate to fit into a straight line, which is called soil line (Wu et al., 2014a). Three spectral indices (PVI, PDI and VAPDI) used in this study were based on the concept of soil line. Eight hundred pure bare pixels (NDVI <0.1) were identified by visual interpretation

Full subset selection
According to the least squares method, the full subset selection is a method to select the optimal combinations by traversing all possible combinations of the full IV set (the IV set includes 12 bands and 32 indices). The calculation of the full subset selection takes the following steps: When the number of IVs is P, P models can be built according to K, the number of IVs input into the model (1 ≤ K ≤ P, K is integer), and there are C K P combinations of IVs for each model (Zhang et al., 2019b). Therefore, based on the calibration dataset, the optimal combination of the IVs for each model was selected according to the maximum of the adjusted coefficient of determination (R 2 adj ). Afterwards, based on the R 2 adj , root mean square error (RMSE), mean absolute error (MAE), Akaike information criterion (AIC) and Bayesian information criterion (BIC), the optimal model from P models at each depth was selected on the ground of the validation dataset. Considering the computational magnitude problem of full subset selection, the value of K was taken from 2 to 6.
Among the five criteria, R 2 adj can improve the accuracy of the comparison between the models with different numbers of IVs and samples. As the number of IVs in the model increases, R 2 adj will not necessarily increase (Srivastava, Srivastava & Ullah, 1995), which mitigates the difference among the coefficient of determination caused by the number of IVs. RMSE and MAE are indicators to evaluate the model inversion error; AIC and BIC measure the goodness of model fit. Smaller values of AIC and BIC mean the model can explain the dependent variable (DV) with fewer IVs (Atkinson et al., 2012). The equations are shown in Eqs. (2)-(6) whereŷ i , y i and y are the predicted, measured, and the average of measured values of the model, respectively; n is the number of samples; k is the number of free parameters in the model; RSS is the squared sum of the residuals between the measured and predicted data; σ 2 is the error variance.

Construction of machine learning models
Four machine learning algorithms, BPNN, SVM, ELM and RF were selected for SSC estimation. Figure 3 is the flowchart of the proposed methodology of SSC estimation in this study.

BPNN Model
BPNN algorithm, proposed by Rumelhart, Hinton & Williams (1986), has a strong nonlinear mapping capability and can adjust the internal parameters of the system according to the error between the output and actual value via the error back propagation algorithm.
Topologically, the BPNN model consists of three layers: the input, hidden, and output layers (Fig. 4) (Wang et al., 2018). After extensive pre-testing, the BPNN model in this study used the optimal combination of IVs as the input layer, SSC as the output layer, and the number of hidden layers was set as 2. The transfer functions of the input and output layers were linear, and the hidden layers were tangent-S. The target error and network learning rate were 0. 65×10 −3 and 0.05, respectively. In order to eliminate the effect of different dimensions on data analysis, the input layer and output layer data were normalized (so were the other three models). MATLAB was used to build the BPNN model (so were the other three models). Details of BPNN can be found in Xiao et al. (2020) and Chen et al. (2015).

SVM Model
SVM is a machine learning algorithm based on the principle of structural risk minimization. It focuses on transforming the input data into a high-dimensional feature space using nonlinear transformations for classification and regression (Chen et al., 2015). This algorithm enjoys such advantages as avoiding discrete values, mitigating overlearning and reducing computation. This study adopted the widely used radial basis kernel (RBF) (Zhang et al., 2019b) as the kernel function of SVM. The penalty parameter (C) and the nuclear parameter (g ) of the RBF have a great effect on the model stability, so a grid-searching technique was used to find the best parameters of C and g (at 0∼20 cm 20∼40 cm 40∼60 cm ∼40 cm and 0∼60 cm, the C was 724, 1024, 1024, 1024 and 1024, respectively, the g was 0.0313, 0.0028, 0.0039, 0.0156 and 0.01, respectively).

ELM Model
ELM is a machine learning algorithm proposed by Huang, Zhu & Siew (2006). It is the same as BPNN in structure, a traditional three-layer neural network (Fig. 4). Its main difference is that the execution does not require adjustments to the input weights of the network or the biases of the hidden elements. Therefore, ELM can reduce the influence of such subjective factors as choice of parameters, and speed up computation while ensure accuracy (Ahila, Sadasivam & Manimala, 2015;Prasad et al., 2019). After extensive pre-testing, this study adopted sigmoid as the activation function, and the number of hidden nodes was set as 6. Detailed principles of ELM can be found in Zhu et al. (2020).

RF Model
RF is a machine learning algorithm proposed by Breiman (2001). Based on multiple decision tree theory, this algorithm can be used for classification and regression. RF uses the bootstrap method to extract training sets from the input data, and randomly generates variables to build decision tree models. Thus, the decision made by a random forest model is based on the ensemble of decisions made by numerous decision trees (Fig. 5) (Zhou et al., 2020;Du et al., 2015). After extensive pre-testing, the minimum number of observations per tree leaf (minleaf) and the number of decision trees (ntree) were set as 5 and 10, respectively (the two parameters were determined by the out-of-bag errors and training set cross-validation) (Zhang et al., 2019a). Model accuracy evaluation R 2 adj , RMSE, and MAE were used to evaluate the inversion of calibration and validation model. Among the three criteria, R 2 adj can avoid the errors caused by the different number of IVs in each model. The closer its value is to 1, the better fit the model has. The smaller the RMSE and MAE are, the smaller the deviation between the predicted and measured values are. The equations are shown in Eqs. (2)-(4).

Analysis of correlation between SSC and IVs
We selected 12 bands of Sentinel-2 and 32 spectral indices to form the IV set for SSC inversion. Based on the calibration dataset, the correlation between the IVs and SSC was analyzed, as is shown in Fig. 6.
The significance level between the IVs and SSC was tested according to the significance testing table of correlation coefficient. When the degree of freedom was 78 and the absolute value of the correlation coefficient (R) was greater than 0.221, the significance level reached 0.05, and when R was greater than 0.288, the significance level reached 0.01. At all five depths (0∼20 cm, 20∼40 cm, 40∼60 cm, 0∼40 cm and 0∼60 cm), the significance level between eight IVs (B6, SI2, S1, S2, S6, BI, SMMI, and PDI) and SSC failed to reach 0.05, that between six IVs (B7, B8, B8A, Int2, NSDSI3, and VSDI) and SSC reached the 0.05, and that between the remaining 30 IVs and SSC reached 0.01.

Optimal combination of IVs based on full subset selection
In order to identify the optimal combinations of IVs, we took the SSC as the DV, and the least squares as the method for data fitting on the basis of the calibration dataset. Then we obtained the optimal combination when the number of IVs was 2 to 6 at the five depths (0∼20 cm, 20∼40 cm, 40∼60 cm, 0∼40 cm and 0∼60 cm) and calculated the R 2 adj , RMSE, MAE, AIC and BIC based on the validation dataset (Table 4). R 2 adj demonstrated a tendency of rising and then falling when the number of IVs increased at most depths, and this tendency was especially obvious at the depths of 20∼40 cm, 40∼60 cm and 0∼60 cm. This indicated that the model fitting became better when the number of IVs was on the rise, but too many IVs would complicate the model, leading to some overfitting. RMSE and MAE displayed the opposite tendency of R 2 adj : they first fell and then rose. This indicated that a local optimal solution might be generated and the validation model error might be increased when there were too many IVs. The tendency of AIC and BIC was almost the same, when the number of IVs increased, there was a first-falling-and-then-rising tendency at 20∼40 cm, 40∼60 cm and 0∼60 cm, and a tendency of increase at 0∼20 cm and 0∼40 cm.
At 0∼20 cm, when the numbers of IVs were 2 (IV-2) and 6 (IV-6), the R 2 adj was the highest (both were 0.44), and those of the RMSE and MAE were relatively close. However, the AIC of IV-2 was half of that of IV-6, suggesting that IV-2 was able to explain the DV with fewer IVs. Therefore, B1 and NSDSI3 were the optimal combination of IVs at 0∼20 cm. At 20∼40 cm, the R 2 adj of the IV-3 was 0.58, which was significantly higher than that of the other combinations. The remaining criteria of IV-3 were the lowest, indicating a small error and a high goodness of fit. Therefore, S5, SI-T and NDDI were the optimal combination of IVs at 20∼40 cm. The tendency at 40∼60 cm was the same as that at 20∼40 cm. IV-3 had the highest R 2 adj , relatively lower errors and higher goodness of fit at the same time. Therefore, S6, VSDI and NDDI were the optimal combination of IVs at 40∼60 cm. At 0∼40 cm, the R 2 adj of IV-2 to IV-4 were all above 0.45 but the AIC and BIC of Notes. *Significance level of 0.05. ** Significance level of 0.01. Bold text represents the optimal combination out of the optimal combinations of IVs at each depth.
IV-2 were the lowest among the three. IV-2 needed the smallest number of IVs to achieve a better fit. Therefore, B1 and NSDSI3 were the optimal combination of IVs at 0∼40 cm. At 0∼60 cm, the R 2 adj of IV-3 to IV-6 were all above 0.44, and the AIC and BIC of IV-3 were the lowest (3.6 and −133, respectively). The RMSE and MAE were relatively low. Therefore, B1, PVI and NSDSI3 were the optimal combination of IVs at 0∼60 cm.

Analysis of BPNN model
The optimal combination (after the full subset selection) of IVs at each depth was input into the BPNN model for training to obtain the SSC inversion model. The accuracy of calibration and validation models are shown in Table 5.
The model performance was optimal at 20∼40 cm, with the R 2 adj of around 0.6 for both calibration and validation models, indicating a good fitting and generalization performance. Its RMSE (0.15%) and MAE (0.11%) were relatively low, indicating their good control of inversion errors. The model had the worst performance at 0∼20 cm, with the lowest R 2 adj (below 0.48) of the calibration and verification model, and the RMSE (both above 0.17%) and MAE (0.13%) were the highest at the five depths. At 40∼60 cm, the model performance was second only to 20∼40 cm, with a relatively good fitting and low inversion error. The performance of the model at 0∼40 cm and 0∼60 cm was similar. The R 2 adj and MAE were around 0.51 and 0.12%, respectively, but the RMSE at 0∼40 cm was greater than that at 0∼60 cm, suggesting that some samples had some relatively big errors.
Overall, the BPNN model worked best at 20∼40 cm, followed by 40∼60 cm, and worked worst at 0∼20 cm. The other two depths had better similar results.

Analysis of SVM model
The optimal combination (after the full subset selection) of IVs at each depth was input into the SVM model for training to obtain the SSC inversion model. The accuracy of calibration and validation models are shown in Table 6.
The model performance was optimal at 40∼60 cm, with R 2 adj of 0.58 and 0.55 for the calibration and validation models, respectively, and its RMSE and MAE were the lowest. The performance of the model at 20∼40 cm and 0∼40 cm was similar, second only to 40∼60 cm. Comparatively, 20∼40 cm was slightly better in fitting of validation model. At 0-20 cm, the model performance was still the worst. Its RMSE (above 0.16%) and MAE (above 0.12%) indicated a bad control of inversion error. But its R 2 adj was around 0.5, the model at this depth still had a good fitting. At 0∼60 cm, the R 2 adj of the calibration and verification models were 0.59 and 0.52, respectively, showing a slight overfitting.
Overall, the SVM model worked best at 40∼60 cm, followed by 20∼40 cm and 0∼40 cm. The worst performance of the model was found at 0∼20 cm.

Analysis of ELM model
The optimal combination (after the full subset selection) of IVs at each depth was input into the ELM model for training to obtain the SSC inversion model. The accuracy of calibration and validation models are shown in Table 7.
The model performance was optimal at 20∼40 cm, with the highest R 2 adj (above 0.6) and lowest inversion error. The model was still the worst at 0∼20 cm, mainly because the inversion error was the largest. But it still had a good fitting (R 2 adj ≈ 0.5). The performance of the model at 0∼60 cm was slightly better than that at 0∼20 cm. The model performance was satisfactory and similar at 40∼60 cm and 0∼40 cm.
In general, the ELM model worked best at 20∼40 cm, but worst at 0∼20 cm and 0∼60 cm (0∼60 cm was slightly better). This model had relatively good performance at the other two depths.

Analysis of RF model
The optimal combination (after the full subset selection) of IVs at each depth was input into the RF model for training to obtain the SSC inversion model. The accuracy of calibration and validation models are shown in Table 8.
The RF model performed well at all five depths, with R 2 adj all above 0.5. The model performance was still optimal at 20∼40 cm. The R 2 adj of the calibration and validation models were 0.68 and 0.63, respectively, which were significantly better than at the other depths. Its RMSE (below 0.14%) and MAE (below 0.11%) indicated that the inversion error was well controlled. At 0∼40 cm, the model performance was second only to that at 20∼40 cm. The result was worst at 0∼60 cm. Its RMSE and MAE were around 0.15% and 0.12%, respectively, indicating a relatively high inversion error. However, its R 2 adj also reached 0.53. The model performance at 0∼20 cm was slightly better than that at 0∼60 cm because the former's R 2 adj of the calibration and validation models were higher (0.03 and 0.01 higher, respectively) than that of the latter. There was some overfitting at 40∼60 cm.
In general, the RF model had a good performance at all depths, working best at 20∼40 cm. The performance was relatively poor at 0∼20 cm and 0∼60 cm (0∼20 cm was slightly better). The scatterplot of measured and predicted SSC of the four models at the best depth are shown in Fig. 7.

Evaluation of inversion depths
The inversion performance of each model at different depths has been discussed in detail (the Section of Model calibrations and validations). In this section, the sensitivity of Sentinel-2 to SSC at different depths was evaluated by analyzing the combined performance of all models at each depth (Fig. 8). At 20∼40 cm, the R 2 adj was significantly higher than that of other depths (Figs. 8A-8E), and each model was able to achieve a good fitting and generalization performance. At this depth, the RMSE and MAE of the calibration models were the lowest, and those of the validation models were also relatively low (Figs. 8F-8O). It adj of the models at 0∼20 cm, 20∼40 cm, 40∼60 cm, 0∼40 cm and 0∼60 cm, respectively. (F-J) RMSE of the models at 0∼20 cm, 20∼40 cm, 40∼60 cm, 0∼40 cm and 0∼60 cm, respectively. (K-O) MAE of the models at 0∼20 cm, 20∼40 cm, 40∼60 cm, 0∼40 cm and 0∼60 cm, respectively.
Full-size DOI: 10.7717/peerj.10585/ fig-8 indicated that all models were able to control the inversion errors at 20∼40 cm. Therefore, 20∼40 cm was the optimal depth for the Sentinel-2 data for SSC inversion. The R 2 adj at 0∼20 cm and 0∼60 cm was relatively low overall, and the inversion error at 0∼20 cm was the highest among all depths (Fig. 8). The overall fitting at 40∼60 cm and 0∼40 cm was satisfactory, and the inversion error at 40∼60 cm was better controlled (Fig. 8). Therefore, Sentinel-2 had the lowest sensitivity at 0∼20 cm. At 0∼60 cm, it was slightly better than 0∼20 cm. The inversion performance was good at 40∼60 cm and 0∼40 cm, and that at 40∼60 cm was relatively better.

Evaluation of inversion models
In this section, the SSC inversion ability of each model was evaluated by analyzing the combined performance at all depths of each model (Fig. 9). The overall R 2 adj of the RF model was higher than that of the other models (Figs. 9A-9D). Although some overfitting occurred at 40∼60 cm, the model still had a good fitting and generalization performance. The RF model had the lowest RMSE and MAE at all depths (Figs. 9E-9L), and the inversion error was especially well controlled at 0∼20 cm where the inversion errors of other models were all high. Therefore, the RF model was the optimal model for SSC inversion. The BPNN model fitted relatively poorly at most depths (Figs. 9A-9D), and the MAE of BPNN model was high at all depths (Figs. 9I-9L). Therefore, the BPNN model had a relatively slightly poor SSC inversion ability. The SVM and ELM models had better fitting and generalization performance, and their inversion errors were relatively low (Fig. 9). The ELM model was obviously better for the inversion at 20∼40 cm than that of SVM model. Therefore, both of them have satisfactory SSC inversion performance though ELM model was slightly better.
When R 2 adj , RMSE and MAE were all taken into consideration, the SSC inversion capability of all models was as follows: RF model > ELM model > SVM model > BPNN model.

SSC distribution of JIA
The optimal model (RF) was used to estimate the SSC distribution at 5 depths (0∼20 cm 20∼40 cm 40∼60 cm ∼40 cm and 0∼60 cm) of JIA (Fig. 10). The study area was interspersed with salinized soil in different degree. It was dominated by non-saline (SSC < 0.2%) and slightly saline soil (0.2% < SSC < 0.5%). The severely saline (0.5% < SSC < 1%) and saline soil (SSC > 1%) only account for a small portion and mainly distributed in the northwest of the area. The salinization in the south of the area was lower than that in the north, which may be related to the irrigation method of JIA (the water was drained from the south to the north, and salt was accumulated in the north). There was more slightly and severely saline soil at 0∼20 cm than at 20-40 cm and 40-60 cm, and more non-saline soil at 20∼40 cm than at other depths. Overall, the estimated SSC distribution of JIA in this study was consistent with the actual measured information (Huang et al., 2018).

DISCUSSION
In this study, the full subset selection was used to select the optimal combination of IVs (included 12 bands and 32 spectral indices) at five depths, which mitigated the subjectivity of IV selection. In addition, R 2 adj was used to evaluate the fitting so as to mitigate the difference among the coefficient of determination caused by the number of IVs. Therefore, the reliability of the comparison of inversion performance at different depths was improved. This study showed that Sentinel-2 was most sensitive to SSC at 20∼40 cm, followed by 40∼60 cm, and the sensitivity at other depths from high to low was 0∼40 cm, 0∼60 cm and 0∼20 cm. A similar result was obtained by Zhang et al. (2019b) when studying the sensitivity to SSC at different depths based on GF-1. It has been found that SSC around crop roots can affect the crop growth by producing osmotic stress (Chen et al., 2003). When studying the water absorption model of crop root system in salinization soil in HID, Qiao (2005) found that the main water absorption layer of sunflowers (accounting for 70% of the samples in our study) was at 0∼50 cm, and the peak of maximum water absorption was at 20∼40 cm. When the sunflower was in bloom (mid to late August), the surface soil moisture content could not meet the root demand, and the peak of maximum water absorption shifted to 35 cm. During this period, the water absorption rate at 20∼40 cm and 40∼60 cm was 2 and 1.5 times of that at 0∼20 cm, respectively, and the water increment mainly came from the deep soil layer. Therefore, the SSC at 20∼40 cm had the strongest effect on crop growth (crop growth can be reflected indirectly via remote sensing data), followed by 40∼60 cm, which was basically consistent with the results in our study. Four machine learning algorithms (BPNN, SVM, ELM, and RF) were used for SSC inversion, and the RF model was found to perform well at all depths and was the optimal