The study area is located in Shanxi Province (110°14′-114°33′ E, 34°35′-40°45′ N), Central China (Figure 1A). This region is in the semi-arid continental monsoon climate zone, which is divided into three climate zones, namely, cool-temperate, mid-temperate, and warm-temperate zones based on the latitudinal extent. The mean annual temperature varies from 8°C in the north to 12°C in the south, and the accumulated temperature (≥10°C) ranges from 1,850 to 5,480°C. The accumulated precipitation varies from 410 to 690 mm. The major soil types in Shanxi Province are cinnamon soil, fluvo-aquic soil, castanozem soil, and castano-cinnamon soil according to Chinese Soil Taxonomy (Figure 1B), with the pH level as “alkaline.” The elevation ranges from 245 m to 3,058 m, and more than 70% of the total area is covered with mountains (Figure 1C). High-quality farmland is mainly located in the six basins of Shanxi Province, which are the Datong basin, Xinding basin, Taiyuan basin, Changzhi basin, Linfen basin, and Yuncheng basin (Figure 1C).

A total of 302,045 points of soil samples were sampled during 2006–2008 (Figure 1D), and the temporal change during the period of sampling was ignored in the study. At each location, a composite sample of 0–20 cm soil layer from three sampling points was collected within 1 m, and an undisturbed surface soil was collected using a metallic core cylinder of 100 cm³ volume. Soil total nitrogen (SN) was measured using the Kjeldahl method, standing for the SN content of around 2007. Soil bulk density (SBD) was measured using the oven-dry method. Additionally, a total of 6,128 points of topsoil samples (0–20 cm) sampled during 2017 were collected (Figure 1E). The SN and SBD were also measured using the Kjeldahl method and the oven-dry method, respectively. The SND (kg m⁻²) of 2007 and 2017 were calculated as follows (Tian et al., 2006): (1)SND= L×SN×SBD×(1-F/100)100, where L is soil thickness (cm), SN is soil total nitrogen concentration (g kg⁻¹), SBD is soil bulk density (g cm⁻³), and F is gravel (>2 mm) percentage. In this study, F was taken as 0 because large gravel rarely appeared in homogeneous loess (Liu et al., 2011). Subsequently, the temporal changes of SND during 2007–2017 (SND_T, kg m⁻² 10a⁻¹) could be calculated by the difference between SND of 2017 and SND of 2007. Finally, the spatial distribution of SND_T with 30 m resolution in Shanxi Province was calculated using the ordinary kriging without anisotropy in ArcGIS.

In the study, three soil sampling transects, which were mainly under the Datong basin located in the cool-temperate zone, the Taiyuan basin located in the mid-temperate zone, and the Linfen basin located in the warm-temperate zone, were established to 200 km long and with 1 km intervals. The SND_T content along the three transects was extracted from the spatial distribution of SND_T.

A total of 179 meteorological stations in and around Shanxi Province from 2007 to 2017 were collected from Climatic Data Center, National Meteorological Information Center. The average temperature (AAT), accumulated precipitation (AAP), averaged sunshine during (SSD), and relative humidity (RH) during each year were obtained, and averaged annual values of the four factors during 2007–2017 were calculated. The spatial distributions of AAT, AAP, SSD, and RH in Shanxi Province were predicted by the ANUSPLIN Meteorological Interpolation Model with 30 m resolution, and the four climatic factors at the three sampling transects were extracted. Altitude was obtained from the digital elevation model (DEM) with 30 m resolution, and slope was calculated from the altitude. The topographic factors, including altitude and slope along the three transects, were also extracted. Net primary productivity (NPP) during 2007–2017 from a moderate resolution imaging spectrometer (MODIS) was downloaded from the website (https://modis.gsfc.nasa.gov/), and the averaged NPP during 2007–2017 was calculated to represent the vegetation conditions.

The land cover types in both 2007 and 2017 were obtained from a moderate resolution imaging spectrometer (MODIS), and the changed and unchanged regions of land-cover types could be acquired. The land cover types and soil types under the three sampling transects were obtained to analyze the SND_T content under different land types and soil types. The traditional statistics of SND_T and its partial correlation coefficients with environmental factors were calculated.

The PWC between SND_T and environmental variables was performed, and the parameter of percentage area of significant PWC (PAS_PWC) was calculated to evaluate the strength of their relations. According to previous studies (Si and Zeleke, 2005; Si, 2008; Zhu et al., 2016), their relationships were recognized as statistically significant, if PAS_PWC is >5%. The influencing factors were decomposed into multiple-scale components by CWT, and the correlation coefficients between SND_T and the environmental factors or their decomposed-scale components were calculated.

Along each transect, the odd points from north to south were selected as the calibration datasets, and all even points were selected as the validation dataset. First, the calibration data of SND_T and environmental factors at the sampling scale were used to develop the model of SMLR to estimate SND_T for validation. Second, the models of partial least square regression (PLSR) were developed from the environmental factors at the sampling and multiple scales, and leave-one-out cross-validation was used to determine the optimal number of latent variables (ONLV) during the calibration procedure. During the validation procedure, the PLSR model was used to predict SND_T for the validation dataset with the ONLV value determined in the calibration. In addition, the predicting models established by each transect from one basin were also validated by the transect from the other basins.

The performance of the SMLR and PLSR models for the calibration and validation data were evaluated using the parameters of R², root mean square error (RMSE), and the ratio of performance to deviation (RPD). R² was calculated as follows: (9)R2=1-∑i=1n(Y^i-Yi)2∑i=1n(Yi¯-Yi), where n is the number of samples and p is the number of the variables in the regression equation. (10)RMSE=1n∑i=1n(Y^i-Yi)2 (11)RPD=SDRMSE, where Y_i, Y^i, and Yi¯ are the measured, predicted, and averaged values of SND_T, respectively. SD is the standard deviation of the measured SND_T. Based on the previous studies (Rossel et al., 2006; Kassim et al., 2021), the performance of the predicting models could be classified into six classes: very poor model, which is not recommended (RPD < 1.0), poor model, in which only high and low values can be distinguishable (RPD = 1.0–1.4), fair model, which may be used for prediction (RPD = 1.4–1.8), good model, which can be used for quantitative predictions (RPD = 1.8–2.0), very good model (RPD = 2.0–2.5), and excellent model (RPD > 2.5).

Descriptive statistics of SND_T under different soil types along three transects are given in Table 1. From north to south, the mean SND_T values were positive and were 0.04, 0.09, and 0.04 kg m⁻² 10a⁻¹, respectively, for cool-temperate, mid-temperate, and warm-temperate transects. Thus, the increment of SND during the period of 10 years along the mid-temperate transect was the greatest among the three transects. The variability of SND_T along cool-temperate and mid-temperate transects were classified as strong because the coefficients of variation (CVs) values were >90%; the variability along the warm-temperate transect was classified as moderate because CV values were between 10 and 90%. Meanwhile, CVs for the three transects were ranked as cool-temperate > mid-temperate > warm-temperate transect. Therefore, the spatial variation of SND_T decreased greatly from cool-temperate to warm-temperate basins.

For fluvo-aquic soil, the mean SND_T values were 0.01, 0.03, and 0.05 kg m⁻² 10a⁻¹ for cool-temperate, mid-temperate, and warm-temperate transects, respectively, which gradually increased from north to south. For castano-cinnamon soil, the mean SND_T values in both cool-temperate and mid-temperate transects were 0.07 kg m⁻² 10a⁻¹. For cinnamon soil, the mean SND_T values were 0.13 and 0.04 kg m⁻² 10a⁻¹ for mid-temperate and warm-temperate transects, respectively. For cinnamon soil and litho soil, the increment of SND during the period of 10 years was the greatest along the mid-temperate transect.

Descriptive statistics of SND_T under different land types along the transects are shown in Table 2. The land type was dominant in grassland along the cool-temperate transect and in cropland along the mid-temperate and warm-temperate transects. For cropland, the mean SND_T values were 0.05, 0.09, and 0.05 kg m⁻² 10a⁻¹ for cool-temperate, mid-temperate, and warm-temperate transects, respectively, and their CV values were around 80%. The results indicated that the spatial variation of SND_T was identical under cropland, even in different climate zones. For grassland and the changed land types, the CVs of SND_T varied greatly under different climate zones.

The spatial distribution of SND_T and their local variance are presented in Figure 2. The spatial distribution of SND_T under different climate zones varied greatly, and the significant variation in the scale-location domain was locally distributed. The areas of significant variation were 36, 31, and 21% for the three basins, respectively. Therefore, the significant variations of SND_T were climate-zone dependent.

The correlation coefficients between SND_T and the influencing factors at the sampling or decomposed scales are given in Figure 3. If the correlation coefficients were >0.18, their relationships were significant at the 99% confidence level. The correlation types between SND_T and environmental factors were different at different scales, and the influencing intensities of environmental variables on SND_T were also different. For cool-temperate transect, the significant correlation was mainly located at scales of 13–40 km. For mid-temperate transect, the significant relationships was mainly located at scales of 18–68 km. However, for warm-temperate transect, the significant correlations were shown at different scales.

The PWCs between SND_T and environmental factors (slope for cool-temperate transect and SSD for mid-temperate and warm-temperate transects were selected; others are not shown) are shown in Figure 4, which indicates the partial correlations between scale components of SND_T and environmental factors at the corresponding scale. Along cool-temperate transect, the scale- and location-specific relationships of SND_T with environmental factors were greater at scales around 60 km, and their significant relations were local at scales <32 km. Along the mid-temperate transect, the scale- and location-specific relationships of SND_T with environmental factors demonstrated that their relations were generally greater than those along the other two transects. Along warm-temperate transect, the correlation of SSD with SND_T was dominant among these influencing factors.

At the original scale, the partial correlation coefficients between SND_T and influencing factors indicated that AAT, AAP, SSD, and NNP were significantly correlated with SND_T along the cool-temperate transect; AAP, SSD, and RH were significantly correlated with SND_T along the mid-temperate transect; and AAP, SSD, RH, and altitude were significantly correlated with SND_T along the warm-temperate transect (Table 3). Therefore, AAP and SSD significantly controlled the spatial distribution of SND_T at the sampling scale of 1 km at different climatic basins.

The evaluation parameters of PWC between SND_T and influencing factors are also given in Table 3. Along the cool-temperate transect, their significant relationships were mainly located at scales of 50–60 km, and the significant areas were more stable at different locations except for NPP. Along the mid-temperate transect, their significant relationships were mainly located at scales of 30–68 km, while the significant areas were local at some locations and scales, which were the greatest among the three transects. Along the warm-temperate transect, the relationships of SND_T with environmental factors were not prominent except for SSD. In brief, the partial correlations between SND_T and influencing factors at the sampling and decomposed scales were different, which indicated that the information contained in SND_T at the sampling and decomposed scales should be complementary.

The SND_T values at the sampling scale were predicted by PLSR, and the evaluating parameters at different numbers of latent variables (NLV) are presented in Figure 5. During the calibration procedure, the RMSE gradually decreased and the R² gradually increased when the NLV ranged from 0 to 50. This indicated that the ONLV in the calibration procedure should be 50. However, the RMSE gradually decreased and then increased, and R² gradually increased and then decreased with the NLV increasing in the validation procedure. This demonstrated that the predicting performance for the validation data was not always improved with the NLV increasing, and the ONLV in the validation procedure was not 50. In the study, the ONLV for the PLSR model was determined by the leave-one-out cross-validation in the validation procedure, which is presented in Figure 5. Although the NLV with the best predicting performance in the validation model was not always in line with the ONLV determined in the calibration procedure, both NLVs were close to each other. Thus, the ONLV determined in the calibration procedure was suitable for predicting the purpose of the environmental factors at the decomposed scales.

The scatterplots observed by predicted SND_T are shown in Figure 6, and evaluating parameters of the predicting models established by both methods of SMLR and PLSR are presented in Table 4. By SMLR from environmental factors, although the included factors were different along different transects, both the factors of AAP and SSD were included for SND_T prediction in the basins. The calibration model could be classified as a fair model because RPD values ranged from 1.4 to 1.8. However, during the validation procedure, only the model for the mid-temperate transect could be considered a fair model, and the other two models were only poor (RPD ranged from 1.0 to 1.4).

The method of PLSR from the environmental factors at the sampling and decomposed scales could significantly improve the predicting accuracy of SND_T, especially for the cool-temperate transect, which had the worst predicting performance by SMLR, and the best performance by PLSR among the three transects. Along the cool-temperate and mid-temperate transects, the predicting performance for the calibration and validation data was excellent, because RPD was >2.5. However, the predicting accuracy along the warm-temperate transect was good, because the RPD values in the calibration and validation procedures ranged from 2.0 to 2.5. The good performance in the cool-temperate and mid-temperate transects could also be proved by their significant PWC between environmental factors and SND_T at different scales (Figure 4). The study also indicated that the predicting models built in each basin were not reliable in the other basins because of the different effects of environmental factors on SND_T in different climate zones (data not shown).

In the spatiotemporal domain, the surface SND varied greatly due to the environmental conditions and human activities. In the study, the indicator of SND_T was established to assess the spatial dependence of temporal changes in SND. Three typical basins located in different climate zones were selected as the study area because of the high-quality arable land dominantly distributed. The mean value of SND_T was positive in all three basins and had the greatest increase in the mid-temperate basin. The result was in line with a previous study (Yang et al., 2015), which proved that the temporal changes of SN presented an increasing trend with time. Previous studies indicated that the rate of denitrification (conversion to N gas) generally increased with the increasing temperature from cool temperature to warm temperature due to the biological activity (Veraart et al., 2011). Meanwhile, the SN input was great in the mid-temperate and warm-temperate basins from fertilizer application or crop straw return because cropland dominated in the two basins. Thus, the high N input and low rate of denitrification in the mid-temperate basin led to a high level of SND_T.

The SND_T demonstrated different statistical characteristics under different soil types, because the predominant N forms were different, resulting in different N cycle processes, including fixation, mineralization, ammonification, nitrification, denitrification, plant assimilation, and so on (Powlson, 1993). The mean SND_T gradually increased from north to south under fluvo-aquic soil, which is developed from river alluvium and widely distributed in these basins. The result might be due to the reason that organic matter decomposes releasing N more quickly in warm humid climates and slower in cool dry climates (Parton et al., 2007). However, the mean SND_T was similar under castano cinnamon soils that resided in different climatic basins. Additionally, a previous study indicated that SNS was influenced by the vegetation type (Marty et al., 2017). In the study, the variability of SND_T was similar under cropland of different climatic basins (around 80%), and also similar under grassland of different climatic basins (around 165%). However, its variability under the changed land types varied greatly resulted from the different land transformations by human. In each basin, the mean SND_T was greater under cropland than that under grassland, and the extract of more SND_T might come from N input from tillage management. For example, the extract of 0.01 kg m⁻² 10a⁻¹ in the cool-temperate basin (0.05 kg m⁻² 10a⁻¹ in cropland compared with 0.04 kg m⁻² 10a⁻¹ in grassland) and the extract of 0.05 in cool-temperate basin (0.09 kg m⁻² 10a⁻¹ in cropland compared with 0.04 kg m⁻² 10a⁻¹ in grassland) might come from tillage practices, such as fertilizer or crop straw return.

In the scale-location domain, the significant variations of SND_T were local and different among different transects. This indicated that the dominant pedogenic processes were different under different climatic basins. The area of significant variation was ranked as cool-temperate > mid-temperate > warm-temperate transects, which is identical to the order of CV values (Table 1). The greatest variation of SND_T in the cool-temperate basin (160%) might be attributed to the various rates of denitrification resulting from temperature differences across the basin and the low human impact resulting from dominant grassland. Therefore, the spatial variation of SND temporal changes decreased greatly with decreasing latitude. However, many errors ignored in the study should be considered in the future, such as the discrepancy in the number of samples between 2007 and 2017, the longer period of sampling time around 2007, and the like.

Some studies hypothesized that each scale of soil properties was controlled by one or more environmental variables at the original sampling scale. For example, the controls of SOC variation at different scales were unraveled by the environmental variables at the original scale by Zhou et al. (2021). However, SND_T may be controlled by the environmental factors at the other scales that differ from the sampling scale. Some studies identified the controlling factors by separating variations in soil properties as well as in environmental factors at similar scales (Hu and Si, 2013; Zhu et al., 2020a). However, the bivariate relationships might be misleading when both variables are dependent on other variables (Hu and Si, 2021).

In the study, partial correlation at the sampling scale and PWC at the decomposed multiscale were used to unravel the controls of environmental factors on SND_T because of the interdependence between predictor and other variables. The partial correlations between environmental factors and SND_T at the sampling scale differed from their partial relations at the decomposed scales. For example, altitude and slope had an insignificant effect on SND_T at the sampling scale along the cool-temperate transect but had a significant effect at the scales > 1 km. The results indicated that the controls of environmental variables on SND_T cannot be captured only at the sampling scales, which was in line with the characteristics of other soil indicators (Si, 2008; Hu and Si, 2013; She et al., 2016; Liu et al., 2018). In the study, the selected environmental factors had all significant influences on SND_T in the three basins except the factor of NPP in the warm-temperate transect (PAS_PWC = 3.73). The insignificant relationships between NPP and SND_T in the warm-temperate basin might be attributed to the double-cropping that resulted in intensive human activities.

The significant relationship between SND_T and temperature might result from the fact that temperature was related to N mineralization rate, nitrification, and denitrification (Guntiñas et al., 2012). The study indicated that the correlation type between SND_T and precipitation was unstable at different scales and locations, which might be due to the reason that precipitation could add N inputs from the atmospheric N and may also remove N by leaching and surface runoff of soluble forms (Jarvis et al., 1996). The duration of sunshine was related to crop photosynthesis, which resulted in vegetation N absorption from the soil, and the rate of straw decomposition, which resulted in soil N input (Tian et al., 2019). Relative humidity, which might be linked with gaseous transfer such as N gas and N oxides during nitrification and denitrification processes (Weitz et al., 2001; Abdelaziz and Kim, 2020), could also influence soil N content. The topographic factors of altitude and slope could lead to N loss and deposition, resulting from the surface runoff of soluble N (He et al., 2016; Madrigal Reyes et al., 2019). Vegetation conditions could promote N storage or N absorption or dry deposition from vegetation. In the study, the results indicated that the types of relationships between environmental variables and N greatly differed at different locations or scales.

To predict SND_T at the sampling scale, the environmental factors were decomposed into multiscale components for covariates. The study revealed that some multiscale components of factors were significantly correlated with SND_T at the sampling scale. Thus, more covariates, which might be beneficial for SND_T prediction, could be generated by the CWT. Because the number of covariates is great and the variables might be cross-correlated, PLSR might be suitable for SND_T prediction (Abdi and Williams, 2013). The performance of PLSR is associated with the ONLV, which is usually determined by the leave-one-out cross-validation with calibration data (Hu et al., 2015). Whether ONLV obtained from the calibration procedure is suitable for prediction in each basin or in other basins should be explored. In the study, the results demonstrated that the ONLV determined in the calibration procedure may not be the best but it is very close to the best one, and the prediction model of PLSR could be validated for the dependent data in each basin rather than in other basins of different climatic basins.

The study revealed that the combined method of a multiscale analysis and PLSR could enhance the prediction accuracy of the SMLR method from environmental factors. In the cool-temperate basin, the predicting performance for the calibration and validation samples was excellent, which might be attributed to the reason that the significant PWC between SND_T and environmental variables was stable at scales of around 50–60 km. In the mid-temperate basin, although the prediction of SND_T was excellent, there was a large discrepancy in prediction accuracy between the calibration and validation procedures. This might be due to the local significant PWC between SND_T and environmental variables. In the warm-temperate basin, the prediction accuracy in the calibration and validation procedures was only good, which could be proved by their poor relationships with PWC. This resulted from intensive human activities during double cropping.

Based on the temporal changes of SND, the highlight of the study was to establish the predictive models of its spatial distribution. If the relationships between SND_T and environmental factors were stable as time goes by, the temporal changes of SND across the spatial areas or temporal changes of SNS could be predicted based on the results of the study, and also the content of SND and SNS could be obtained based on their background. However, the correlations between SND_T and climatic factors or vegetation would change with climatic changes. Therefore, the climatic changes should be considered for predicting temporal changes of SND with time in the future.