Mapping of Soil Total Nitrogen Content in the Middle Reaches of the Heihe River Basin in China Using Multi-Source Remote Sensing-Derived Variables

Soil total nitrogen (STN) is an important indicator of soil quality and plays a key role in global nitrogen cycling. Accurate prediction of STN content is essential for the sustainable use of soil resources. Synthetic aperture radar (SAR) provides a promising source of data for soil monitoring because of its all-weather, all-day monitoring, but it has rarely been used for STN mapping. In this study, we explored the potential of multi-temporal Sentinel-1 data to predict STN by evaluating and comparing the performance of boosted regression trees (BRTs), random forest (RF), and support vector machine (SVM) models in STN mapping in the middle reaches of the Heihe River Basin in northwestern China. Fifteen predictor variables were used to construct models, including land use/land cover, multi-source remote sensing-derived variables, and topographic and climatic variables. We evaluated the prediction accuracy of the models based on a cross-validation procedure. Results showed that tree-based models (RF and BRT) outperformed SVM. Compared to the model that only used optical data, the addition of multi-temporal Sentinel-1A data using the BRT method improved the root mean square error (RMSE) and the mean absolute error (MAE) by 17.2% and 17.4%, respectively. Furthermore, the combination of all predictor variables using the BRT model had the best predictive performance, explaining 57% of the variation in STN, with the highest R2 (0.57) value and the lowest RMSE (0.24) and MAE (0.18) values. Remote sensing variables were the most important environmental variables for STN mapping, with 59% and 50% relative importance in the RF and BRT models, respectively. Our results show the potential of using multi-temporal Sentinel-1 data to predict STN, broadening the data source for future digital soil mapping. In addition, we propose that the SVM, RF, and BRT models should be calibrated and evaluated to obtain the best results for STN content mapping in similar landscapes.


Introduction
As a major component of the terrestrial nitrogen (N) pool, soil total nitrogen (STN) not only provides essential nutrients for plant growth but also affects soil function and the concentration of greenhouse gases in the atmosphere. Low  result in loss of nitrogen from the soil, causing soil fertility degradation and water pollution [1]. Soil degradation associated with soil nitrogen loss reduces soil security and greatly contributes to climate change. In total, 21% of the total annual emissions of nitric oxide (NO) are related to soil degassing [2]. On a global scale, NO emissions from soil and from fossil fuel combustion are in the same range [3,4].
In addition, as a key indicator of soil fertility and quality, STN content is closely related to agricultural productivity and food security. Therefore, reliable predictions of STN content are critical to sustaining sustainable agricultural development and understanding the regional N cycle. Up to date STN maps are important to identify spatial variation and control factors of STN, which can help maintain food security and soil security and provide a reference for managing climate change. To address environmental challenges, such as climate change and land degradation, an accurate and efficient method is needed to predict the spatial distribution of STN and improve its prediction accuracy.
Digital soil mapping provides a low-cost and efficient method for predicting the spatial distribution of soil nutrients. Most digital soil mapping methods are based on soil landscape models, which establish mathematical or statistical relationships between soil properties and related environmental variables [5,6]. Many digital soil mapping methods have been used to predict soil properties, including generalized linear model (GLM), random forest (RF), artificial neural networks (ANNs), boosted regression trees (BRTs), support vector machine (SVM), and regression kriging (RK). Although these models have achieved satisfactory results in different fields, choosing the best prediction method for a given landscape has always been a challenge. Jeong et al. [7] used RF and SVM models to conduct STN mapping studies in the Soyang lake watershed in South Korea and found that the RF model performed better while Were et al. [8] found that the latter performed better than the former in Kenya. Akpa et al. [9] reported that the RF model was superior to BRT in the prediction of soil properties in Nigeria while Yang et al. [10] found that the latter had better predictive performance in an alpine ecosystem. Because these studies obtained inconsistent results for different environments, the predictive models should be compared and evaluated for a given landscape.
STN prediction using digital soil mapping techniques requires sufficient environmental information. Commonly used environmental data are on land use/land cover (LULC), the climate, the topography, as well as remote sensing-derived variables. Remote sensing has received recent attention to improve digital soil mapping due to the extensive area it can cover and the rich spatial information it provides. Optical satellite imagery is one of the most commonly used remote sensing data sources for digital soil mapping and has contributed significantly to research on the prediction of STN. Recent studies have also noted the potential of synthetic aperture radar (SAR) for predicting soil chemical properties. For example, Bartsch et al. [11] explored the feasibility of C-band ENVISAT ASAR data in soil organic carbon (SOC) mapping and found that the SOC in the study area can be quantified by SAR images. Ceddia et al. [12] used optical and L-band ALOS PALSAR data to conduct SOC prediction studies in central Amazon, and found that the backscattering coefficient of SAR data can improve the prediction accuracy. A similar study was also reported by Ma et al. [13], who used SAR (Sentinel-1) and optical (Landsat 7) images to map soil properties in eastern China.
A wide range of information can be extracted from SAR data [14], the most common of which is the backscatter coefficient. The backscattering intensity of SAR data can be used to obtain information about surface properties and has been widely used in soil science research to obtain information about soil roughness [15], soil moisture [16], and soil salinity [17]. The application of SAR images in mapping soil properties depends on the sensitivity of backscattering intensity to changes in soil moisture and land surface conditions [18]. Yang et al. [19] found a significant correlation between SAR backscatters and various soil properties (including SOC, STN, and sand) during the growing season, and reported that multi-temporal SAR data is useful for predicting soil chemical properties because it can capture soil-vegetation relationships. Similar results have been reported by, for example, Maynard et al. [20] and Takada et al. [21]. The soil-vegetation relationships observed in these studies can explain the spatial variation of the soil and thus help remote sensing techniques to predict soil properties [22,23].
Although SAR data is a promising data source for digital soil mapping, few studies have used the backscatter coefficients of SAR data as predictor variables for STN mapping.
Arid and semi-arid areas account for about 50% of China's land surface. Due to irrational land use of these areas, as well as a decline in soil fertility and a lack of water resources, desertification in such areas is intensifying. Hence, greater attention should be given to the variation in soil properties here. Representing both arid and semi-arid regions in China, the spatial distribution of soil properties in the Heihe River Basin (HRB) has attracted the attention of researchers in recent years. However, these studies have focused primarily on SOC [24,25] and soil texture [26], with little information on STN in the area. This study used machine learning algorithms to predict the STN content in the middle reaches of the HRB in Northwest China. The specific objectives were (i) to explore whether adding SAR data improves STN prediction; (ii) to evaluate and compare the performance of RF, BRT, and SVM models; and (iii) to assess the importance of predictor variables in mapping the STN content of the study area. For these purposes, 15 predictor variables (including the climate, the topography, LULC, and remote sensing-derived variables) and 85 soil samples were obtained. The remote sensing images used in this study included optical (Landsat-8) and multi-temporal SAR (Sentinel-1A) data.

Study Area
The study area is in the middle reaches of the HRB (100.23 • E, 39.06 • N), which is located in Zhangye City, Gansu Province, in northwestern China (see Figure 1). HRB is China's second largest inland river, covering an area of about 1.3 × 105 km 2 [27], with a total population in the middle reaches of 82.55 × 10 4 in 2010 [28]. The study area exhibits a continental dry temperate climate, with a mean annual temperature (MAT) of approximately 6 to 8 • C and a mean annual precipitation (MAP) of 150 mm [29]. The elevation of the area is between 1300 and 1700 m above sea level. Major land use types are cultivated land, grassland, and barren land, and the main crops include corn and spring wheat. Precipitation in the middle reaches of the HRB is limited, with farmland accounting for about 86.85% of water consumption in the basin [30]. Gray-brown desert soil and gray desert soil are the zonal soil types in the study area [31].

Soil Data
A total of 85 soil samples (0-20 cm) were used to calibrate and validate the models (Figure 1), which were collected and provided by "Cold and Arid Regions Science Data Center at Lanzhou (CARD)". Large-scale field sampling is time-consuming and labor-intensive. In order to obtain sufficient soil data to characterize the variation of soil properties, the purposive sampling strategy of

Soil Data
A total of 85 soil samples (0-20 cm) were used to calibrate and validate the models (Figure 1), which were collected and provided by "Cold and Arid Regions Science Data Center at Lanzhou (CARD)". Large-scale field sampling is time-consuming and labor-intensive. In order to obtain sufficient soil data to characterize the variation of soil properties, the purposive sampling strategy of Zhu et al. [32] was adopted. Soil sample sites were selected based on major soil formation factors, including topographical conditions, the climate, and LULC. Soil sampling was carried out from 2011 to 2014, and a 1-m deep soil pit was dug at each location. Coordinate information for each sample point was collected using a handheld global positioning system receiver, and environmental characteristics (including land use, vegetation, and elevation, etc.) were recorded. The samples were air dried, ground, and sieved at 2 mm, after which the STN content was determined using the Kjeldahl method.

Environmental Variables
The environmental variables used in this study for STN mapping included LULC, remote sensing-derived variables, and topographic and climate variables. These environmental data were collected from various sources and used to generate a total of 15 predictor variables. ArcGIS 10.2 was used to convert all environment variables to raster formats with the same 30-m resolution. These predictor variables and observed STN data were then imported into a geographic information system in a common coordinate reference system for future STN mapping. The pixel values of the predictor variables corresponding to each soil sample point based on these raster layers were calculated to build the models.

Remote Sensing Variables
The remote sensing images used in this study included SAR and optical data. The SAR image used in this study was Sentinel-1A that was downloaded from the ESA (European Space Agency), whereas the optical data was Landsat-8 OLI that was downloaded from the USGS (US Geological Survey). Designed by the ESA, Sentinel-1 is composed of two satellite constellations, including Sentinel-1A and Sentinel-1B [33]. The Sentinel-1A carries a C-band SAR imaging instrument that provides four imaging modes. In this study, 4 Sentinel-1A images (single-look complex (SLC) products) in the IW (interferometric wide swath) mode were obtained. More detailed information of the SAR data used in this study is provided in Table 1. We preprocessed these SAR data with SARscape 5.2, including multi-look, coregistration, speckle filtering (a 5 × 5 window Lee filter [34]), geocoding, and radiometric calibration [14]. These images were geocoded using the ASTER GDEM and their digital numbers (DN) were converted to decibel (dB) scale backscatter coefficients. Landsat-8 OLI images with cloud cover <10% from July to September 2015 were obtained; to represent vegetation intensity and type [6,10], the normalized difference vegetation index (NDVI) was calculated using bands 4 and 5. The relief correction of all images was performed based on the polynomial geometric correction method. Because Landsat-8 OLI's red band 4 (B4)-(0.64-0.67 µm), near-infrared band 5 (B5)-(0.85-0.88 µm), and shortwave infrared band 6 (B6)-(1.57-1.65 µm) represent vegetation growth, coverage and biomass [35], respectively, the three bands of Landsat-8 OLI were extracted as predictor variables to map STN. A total of eight environmental covariates were extracted from remote sensing images, four of which were from optical images and the remainder from SAR data.

Climate Variables
The MAP and MAT data for the study area from 1961 to 2010 were obtained from CARD (http://card.westgis.ac.cn/). We calculated the average MAP and MAT for 50 years (1961-2010) as predictor variables. These data were interpolated (high accuracy surface modeling (HASM) method [36]) from the measured values of 34 meteorological stations, 21 of which are conventional meteorological stations in the HRB and its surrounding areas, whereas the remainder are national reference stations around the Heihe River. The HASM method has been reported to improve the interpolation of climate variables compared to other conventional techniques [37].

Land Use/Land Cover Data
The LULC data were generated using Landsat TM and ETM remote sensing data from 2011 combined with field surveys and validation. The data were provided by CARD, classified into the following types of LULC: Farmland, grassland, barren, wetlands, forests, and villages.

Topographic Variables
A total of four topographic variables were calculated from the 30-m resolution ASTER GDEM (including elevation, slope, aspect, and topographic wetness index (TWI)) using SAGA GIS and ArcGIS 10.2. SAGA TWI was reported to predict more realistic soil wetness than the traditional TWI [10]. In this study, TWI was obtained using SAGA GIS and the remaining topographic variables were generated by ArcGIS 10.2.

Prediction Models
We applied three machine learning algorithms (i.e., RF, BRT, and SVM models) to predict STN. These models have a strong ability to model complex nonlinear relationships between soil properties and environmental variables. We optimized the parameters in these models through the grid search approach using the 'caret' package in the R software.
First introduced by Vapnik et al. [38], the SVM model originated from statistical learning theory and was applied to classification or regression. Based on the kernel function, the SVM model projects the input data onto the new hyperspace, optimally dividing all the data into different classes [39]. The four commonly used kernel functions are as follows: Linear, sigmoid, polynomial, and radial basis function (RBF). The selection of kernel functions can affect the prediction accuracy of the SVM model. The RBF kernel performs best at capturing nonparametric features [40] and has been successfully applied to soil mapping studies [5]. In this study, RBF was selected as the kernel model of the SVM model.
As a nonparametric method, the RF model is an ensemble machine learning algorithm used to perform classification and regression. During model training, the RF model generates a large number of random trees that are combined into a single prediction [41]. A unique bootstrap sample from the original training data was used to build each tree in the forest [42]. The samples in the RF method are divided into "in-bag" samples and "out-of-bag (OOB)" samples, and the latter are used to estimate errors [43].
Developed by Elith et al. [44], the BRT model is a machine learning technique that combines the advantages of boosting and regression trees. The regression tree is a model based on a decision tree that analyzes the variation of the response variables of a set of predictor variables [45]. As a machine learning technique similar to model averaging [46], boosting is a forward and stage-wise procedure that uses an iterative approach to develop the final model, gradually adding the tree to the model [47]. The BRT model relies on a stochastic gradient boosting procedure to improve model prediction accuracy and prevent data overfitting [48].

Statistical Analyses
This study used SPSS 24.0 software to perform a descriptive statistical analysis of STN and environmental variables to summarize these data. The following packages were used in the R software for STN mapping: 'gbm'-package (BRT), 'randomForest'-package (RF) and 'kernlab'-package (SVM), and 'raster'-package and 'maptools'-package.

Model Validation
To explore the effectiveness and potential contribution of the backscatter coefficient of multi-temporal Sentinel-1 images to STN mapping, different combinations of environmental variables were constructed ( Table 2). We then compared the predictive performance of different combinations based on RF, BRT, and SVM methods. To compare and evaluate the prediction accuracy of these models, a 10-fold cross-validation technique was used to calculate the following three commonly used validation indices: The root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R 2 ). These validation indices were defined as shown in Equations (1)-(3): where n represents the number of sample points; and P i and O i represent the estimated and observed STN content at site i, respectively.

Descriptive Statistics
The statistical results of the STN content and the values of the predictor variables corresponding to the samples are shown in Table 3. The STN content in the topsoil ranged from 0.06 to 1.68 g kg −1 (the mean and median were 0.87 and 0.96 g kg −1 , respectively), with a standard deviation (SD) of 0.34 g kg −1 . The SD value of the STN content was lower than the mean value, indicating a moderate variation in its distribution [49]. The STN content (with a skewness coefficient of 0.33) had a slightly skewed distribution.

Model Performance
We constructed five STN content models based on different combinations of environmental variables: Model A represents environmental variables without SAR data while Model E was a combination of all environmental variables; Models B and C represent optical and SAR remote sensing variables, respectively; and Model D included variables derived from Landsat-8 and Sentinel-1 images. Table 4 shows the model performance statistics for the BRT, RF, and SVM techniques using these STN models. The BRT model obtained the highest accuracy to predict STN, with the highest R 2 (0.57) value and the lowest RMSE (0.24) and MAE (0.18) values. The predictive performance of the BRT and RF models was significantly better than the SVM model. However, the prediction accuracy of the BRT model was closely followed by the RF model (R 2 = 0.55, RMSE = 0.24, and MAE = 0.19). Overall, for the RF, BRT, and SVM methods, Model C (combination of backscatter coefficients from multi-temporal SAR data) had a higher prediction accuracy than Model B (single-day Landsat-8-derived variables). This result indicates that SAR data not only provides important information for STN mapping but can also be as effective as optical images. Moreover, the predictive power of multi-temporal Sentinel-1A data was not inferior to the single-day Landsat-8 OLI image. Therefore, SAR data has broad application prospects in digital soil mapping, especially in those areas that are susceptible to cloud. For the three C models, the RF technique obtained the highest R 2 (0.44) and the lowest RMSE (0.27) and MAE (0.22). The R 2 value indicated that Model C, which was constructed using multi-temporal Sentinel-1A images, explained approximately 44% of the STN variation.
Using three different machine learning techniques, STN prediction accuracy was improved when multi-temporal Sentinel-1A images were combined with Landsat-8 images. The addition of SAR data using the BRT method resulted in 48.6%, 17.2%, and 17.4% improvements for R 2 , RMSE, and MAE, respectively. Improvements in the prediction accuracy of the STN were also observed when multi-temporal SAR data were added to Model A. For example, the addition of multi-temporal Sentinel-1A images using the BRT and RF models increased R 2 from 0.52 to 0.57 and 0.52 to 0.55, respectively. The combination of all environmental variables (Model E) obtained the highest prediction accuracy compared to the other four models (Models A, B, C, and D), with R 2 values of 0.57 and 0.55 for the BRT and RF modeling techniques, respectively. The R 2 values indicated that the RF and BRT modeling techniques can explain 55% and 57% of the STN variation, respectively.

Relative Importance of Environmental Data
The relative importance of each predictor variable calculated using Model E in the RF and BRT modeling techniques is shown in Figure 2. In this study, the relative importance of predictor variables was normalized to 100% to facilitate comparisons among predictor variables. The relative importance of the same environmental variables in the RF and BRT models was slightly different. The five most important predictor variables for the RF model were LULC, BC_3, NDVI, Band_6, and Band_4 while in the BRT model, they were LULC, BC_3, BC_2, NDVI, and elevation. In both models, LULC and BC_3 had the same relative importance rankings, namely, first and second, respectively. These findings indicate that these environmental variables are the dominant predictor variables of STN mapping in this study area. In the RF model, Landsat-8-derived variables had the highest relative importance (with a relative importance of 31%), followed by Sentinel-1-derived variables (28%), LULC (15%), topographic (14%), and climate (12%) variables. The relative importance of Sentinel-1-derived variables, LULC, Landsat-8-derived variables, and topographic and climate variables in the BRT model were 31%, 26%, 19%, 19%, and 5%, respectively. Moreover, the remote sensing variables in the RF and BRT models had a relative importance of 59% and 50%, respectively.

Spatial Prediction of STN Content
The spatial distribution maps and corresponding descriptive statistics of STN content predicted by the RF, BRT, and SVM methods using Model E are shown in Figure 3 and Table 5, respectively. Figure 4 shows STN maps predicted using Model D constructed by using only remote sensing variables. The mean (±SD) values of the STN content predicted by the RF, BRT, and SVM methods were 0.81 (±0.23), 0.83 (±0.30), and 0.83 (±0.19) g kg −1 , respectively. The average and SD values of the STN content predicted by the three methods were lower than those of the observed STN content data. This result indicates that the predicted STN content is less variable than the measured STN. In addition, this result is consistent with the findings of Wang et al. [49] and Adhikari et al. [50]. Since the performance of tree-based models (RF and BRT) was significantly better than SVM, we used RF and BRT methods to further explore the differences between STN content predicted by Model E and Model D. Figure 5a,b show the difference between STN content predicted using Model E and Model D (Model E-Model D). The mean (±SD) of these differences in STN content predicted by the RF and BRT methods was −0.03 (±0.07) and −0.01 (±0.13), respectively, indicating that there was only a small difference between the Remote Sens. 2019, 11, 2934 9 of 18 STN content predicted by Model E and Model D in most of the study areas. There was a relatively high difference between the STN content predicted by Model E and Model D in barren land and vegetation coverage areas (see Figure 5c). important predictor variables for the RF model were LULC, BC_3, NDVI, Band_6, and Band_4 while in the BRT model, they were LULC, BC_3, BC_2, NDVI, and elevation. In both models, LULC and BC_3 had the same relative importance rankings, namely, first and second, respectively. These findings indicate that these environmental variables are the dominant predictor variables of STN mapping in this study area. In the RF model, Landsat-8-derived variables had the highest relative importance (with a relative importance of 31%), followed by Sentinel-1-derived variables (28%), LULC (15%), topographic (14%), and climate (12%) variables. The relative importance of Sentinel-1derived variables, LULC, Landsat-8-derived variables, and topographic and climate variables in the BRT model were 31%, 26%, 19%, 19%, and 5%, respectively. Moreover, the remote sensing variables in the RF and BRT models had a relative importance of 59% and 50%, respectively.

Spatial Prediction of STN Content
The spatial distribution maps and corresponding descriptive statistics of STN content predicted by the RF, BRT, and SVM methods using Model E are shown in Figure 3 and Table 5, respectively. Figure 4 shows STN maps predicted using Model D constructed by using only remote sensing variables. The mean (±SD) values of the STN content predicted by the RF, BRT, and SVM methods were 0.81 (±0.23), 0.83 (±0.30), and 0.83 (±0.19) g kg −1 , respectively. The average and SD values of the STN content predicted by the three methods were lower than those of the observed STN content data. This result indicates that the predicted STN content is less variable than the measured STN. In addition, this result is consistent with the findings of Wang et al. [49] and Adhikari et al. [50]. Since the performance of tree-based models (RF and BRT) was significantly better than SVM, we used RF and BRT methods to further explore the differences between STN content predicted by Model E and Model D. Figures 5a,b show the difference between STN content predicted using Model E and Model D (Model E-Model D). The mean (±SD) of these differences in STN content predicted by the RF and BRT methods was−0.03 (±0.07) and -0.01 (±0.13), respectively, indicating that there was only a small difference between the STN content predicted by Model E and Model D in most of the study areas.
There was a relatively high difference between the STN content predicted by Model E and Model D in barren land and vegetation coverage areas (see Figure 5c).

Model Performance
Our findings showed that tree-based models (RF and BRT) can achieve better STN prediction accuracy than SVM. This was consistent with the results of Wang et al.'s [5] soil property prediction study in Australia, which reported that the prediction accuracy of RF and BRT models was higher than that of SVM. Ottoy et al. [51] compared the accuracy of BRT, SVM, and ANN models to predict soil properties in Belgium, and found that the tree-based model (BRT) performed best. Tziachris et al. [52] used different models to map soil properties in the Kastoria area and found that the RF model improved the prediction accuracy compared to ordinary kriging (OK). However, the results of Ding et al. [53] showed that the SVM model was superior to RF in SOC prediction. Gomes et al. [41] found that the RF model was superior to Cubist in SOC mapping in Brazil while Sorenson et al. [54] reported that Cubist performed better than RF in SOC prediction. Although these studies have different comparison results, other digital soil mapping studies in the HRB have also found that the prediction performance of RF models was better than SVM [55,56]. Based on these results, no single machine learning technique is most suitable for all landscapes. Therefore, it is necessary to evaluate and compare the prediction capabilities of different models under different landscape and environmental input variables.

Model Performance
Our findings showed that tree-based models (RF and BRT) can achieve better STN prediction accuracy than SVM. This was consistent with the results of Wang et al.'s [5] soil property prediction study in Australia, which reported that the prediction accuracy of RF and BRT models was higher than that of SVM. Ottoy et al. [51] compared the accuracy of BRT, SVM, and ANN models to predict soil properties in Belgium, and found that the tree-based model (BRT) performed best. Tziachris et al. [52] used different models to map soil properties in the Kastoria area and found that the RF model improved the prediction accuracy compared to ordinary kriging (OK). However, the results of Ding et al. [53] showed that the SVM model was superior to RF in SOC prediction. Gomes et al. [41] found that the RF model was superior to Cubist in SOC mapping in Brazil while Sorenson et al. [54] reported that Cubist performed better than RF in SOC prediction. Although these studies have different comparison results, other digital soil mapping studies in the HRB have also found that the prediction performance of RF models was better than SVM [55,56]. Based on these results, no single machine learning technique is most suitable for all landscapes. Therefore, it is necessary to evaluate and compare the prediction capabilities of different models under different landscape and environmental input variables.
The prediction accuracy (Table 4) of Model C constructed by multi-temporal Sentinel-1A data indicates that the model could explain 44% of the STN variation. Compared with a model constructed by RapidEye (optical imaging) to perform an STN prediction study in India (R 2 = 0.41) [57], the prediction performance of the model generated by multi-temporal SAR data in this study was not inferior. Tian et al. [58] reported a study in Tibet that found that NDVI explained 6% to 29% of the soil nutrient stoichiometry.
In our study, the inclusion of multi-temporal Sentinel-1A data improved the prediction accuracy of the STN. This is expected because the addition of more useful information improves the prediction accuracy of the model. Therefore, our results demonstrate the effectiveness of multi-temporal SAR data on STN mapping in an agro-ecosystem and has the potential to improve prediction accuracy. Some previous studies have also demonstrated the effectiveness of adding other useful information to digital soil mapping. For example, Wang et al. [5] evaluated the effect of the addition of seasonal fractional cover data on SOC mapping in eastern Australia, achieving an improvement in the RMSE of 2.8% to 5.9%. The studies carried out by Wang et al. [49] and Yang et al. [59] respectively evaluated the effects of cultivation history and crop rotation information on predicting soil properties and demonstrated the effectiveness of such information on soil mapping. The most accurate model (Model E) of this study used RF and BRT modeling techniques to explain 55% and 57% of the STN variation, respectively (Table 4). Although the environmental conditions, sampling strategies, and verification methods of this study were different from those of the previous STN prediction studies, the prediction performance of the RF and BRT methods in this study was not inferior compared to these studies. For example, Forkuor et al. [60] also developed an RF model for soil property prediction in Burkina Faso, but it only explained less than 40% of the STN variation. Xu et al. [61] used the RF model for STN mapping in southern India to explain less than 50% of STN variation. In an STN prediction study in Fuyang, Zhejiang Province, China, He et al. [62] found that the RF and BRT models explained approximately 50% of the STN variability. In contrast, Wang et al. [63] used the RF model to predict STN content in northeast China, obtaining R 2 values that were higher than this study that explained 69% of the STN variation. These different predictive performances may be due to differences in the type and quality of the ancillary data, the study area, and the number of field observations.

Importance of Predictor Variables
Our results not only indicate the importance of remote sensing data to predict STN but also emphasize the need to incorporate SAR images. Radar can provide spectral information beyond vegetation cover and soil surface [12], and backscatter signals from SAR images are used to retrieve target properties, such as forest above-ground biomass, soil texture, soil moisture, and salinity. Many studies have shown that information from SAR data, such as the backscatter coefficient, can detect vegetation [64] and soil moisture [65]. There is a significant positive correlation between STN and vegetation in the topsoil [66], and the backscattering coefficient is an important indicator representing vegetation density and biomass. Yang et al. [23] reported that Sentinel-1 images successfully predicted soil properties because of their ability to capture characteristics of short-term vegetation changes. Anne et al. [22] found that soil properties were related to vegetation canopy detected by remote sensing images because the soil was strongly affected by the vegetation. Although Sentinel-1 images have rarely been used as predictor variables for mapping STN, previous studies using these images to monitor vegetation demonstrated the ability of SAR data to capture vegetation information, which can be further applied to map STN because of the relationships in the soil-vegetation system. For example, Guo et al. [67] found that the combination of Sentinel-1 and Sentinel-2 obtained the best above-ground biomass monitoring results while Navarro et al. [68] used Sentinel-1 to obtain satisfactory soil moisture inversion results on the Loess Plateau of China.
A variety of optical satellite images have been applied to soil mapping, such as Landsat, Sentinel-2, and Pleiades-1A. The optical remote sensing variables commonly used for mapping STN are spectral reflectance and derived vegetation indices [69]. Spectral reflectance, which is closely related to vegetation density and biomass, and LULC are important environmental variables that affect STN content. Soil and vegetation have interactive effects, such as spatial and temporal changes in soil nutrients due to the accumulation and decomposition of vegetation biomass. Xu et al. [61] found that spectral reflectance is one of the main factors affecting STN variation. Wang et al. [63] reported that optical remote sensing-derived variables were the most important environmental variables affecting STN variation compared to topographic and climate variables. However, this finding differed from our results, mainly due to the addition of multi-temporal Sentinel-1A images to map STN in this study. LULC patterns are one of the most direct and important factors affecting soil nutrient changes. Many studies have shown that different LULC patterns have a greater impact on soil nutrient content [70][71][72]. Consistent with the findings of Martin et al. [43], it also reported that the LULC data obtained the highest relative importance. In addition, the results of Su et al. [73] and Genxu et al. [74] showed that soil nutrient content in the HRB has strong spatial variation in different LULC types.
As one of the five soil formation factors, terrain indirectly affects the soil by causing redistribution of matter and energy. Therefore, topographic variables are closely related to the spatial variation of soil properties and are often used as a key predictor for digital soil mapping [75,76]. Among all topographic variables, elevation achieved the highest relative importance. Elevation can affect the microclimate at local scales and indirectly affect microbial activity, affecting the decomposition and transformation of STN [77]. A significant correlation between elevation and MAT and MAP was observed in a study of soil mapping in the HRB [24], which further demonstrates that elevation affects regional climate variables. In previous digital soil mapping studies, elevation was also reported to be the most effective topographical variable [55,78]. In addition, slope, aspect, and TWI were also found to be important environmental variables affecting STN distribution in previous studies [79,80].
As the most commonly used predictors to map STN, temperature and rainfall are important climate variables that affect the distribution of STN on regional and continental scales [81,82]. By affecting soil water content, climate variables have been reported to affect not only microbial activity but also plant-microbial interactions that influence N availability [83]. In addition, precipitation indirectly alters N cycling by affecting plant N uptake and plant productivity. This has been supported in other studies that found changes in precipitation caused a shift in the plant community structure [84], which in turn alters the N cycling of the ecosystem [85]. In a study in central New Mexico, Cregger et al. [86] reported that precipitation changes have both direct and indirect impacts on the N cycling of this semi-arid forest.

Spatial Prediction of STN Content
The STN content maps predicted by the three methods had similar distribution patterns and exhibited strong spatial variability. Agro-ecosystems along rivers that are largely affected by humans had higher STN content. Agriculture in this study area relies on irrigation, which is mainly from groundwater or rivers (Heihe River). Song et al. [24] conducted SOC prediction in the HRB and found that the mid-stream farmland ecosystem had a relatively high SOC content. In addition, the agro-ecosystem in the southeast had a slightly higher STN than the northwest. This result can be explained by the relatively high rainfall and lower-than-average temperatures in the southeastern part of the study area. In contrast, other areas of the study had lower STN levels, especially in desert areas without vegetation cover. Previous studies have confirmed the spatial relationship between STN and vegetation [87,88]. For example, Wang et al. [63] used the RF method to predict STN content in northeastern China, and found that areas with dense vegetation cover had higher STN. Similar findings were also reported in the STN mapping studies by Zhang et al. [89] and Wang et al. [90].

Conclusions
In this study, we predicted the spatial distribution of STN content in the middle reaches of the HRB in northwestern China by comparing the BRT, RF, and SVM models. We were able to come to the following main conclusions from this study: (1) Compared to the SVM model, tree-based models (RF and BRT) performed better at predicting STN content; (2) the combination of the multi-temporal Sentinel-1A and Landsat-8 OLI data using the BRT model compared to the application of optical data alone improved the RMSE and MAE by 17.2% and 17.4%, respectively. The inclusion of multi-temporal Sentinel-1A data in the RF and SVM models achieved similar improvements in the prediction of STN. The predictive power of multi-temporal SAR images was as strong as the optical remote sensing variables and was identified as one of the most important predictor variables for the best prediction of STN content in this study area; (3) the combination of all environmental variables achieved the highest prediction accuracy, with the BRT model having the highest R 2 (0.57) value and the lowest RMSE (0.24) and MAE (0.18) values; (4) the most important environmental variables affecting the spatial distribution of STN were the predictor variables extracted from remotely sensed images (including SAR and optical data). These predictor variables in the RF and BRT models explained 59% and 50% of the STN variation, respectively; and (5) the STN distribution maps obtained using the three machine learning techniques had similar distribution patterns, and the agro-ecosystem along the river had a higher STN content. Based on the results obtained in this study, we recommend using multi-temporal Sentinel-1 images to map STN content, especially in those areas that are susceptible to cloud. Although the use of variables derived from remote sensing improved the prediction performance, this study did not achieve a high prediction accuracy. Future research may explore the feasibility of incorporating new environmental variables to improve the prediction accuracy.