Spatial heterogeneity distribution of soil total nitrogen and total phosphorus in the Yaoxiang watershed in a hilly area of northern China based on geographic information system and geostatistics

Abstract Soil total nitrogen (STN) and total phosphorus (STP) are important indicators of soil nutrients and the important indexes of soil fertility and soil quality evaluation. Using geographic information system (GIS) and geostatistics, the spatial heterogeneity distribution of STN and STP in the Yaoxiang watershed in a hilly area of northern China was studied. The results showed that: (1) The STN and STP contents showed a declining trend with the increase in soil depth; the variation coefficients (C v) of STN and STP in the 0‐ to 10‐cm soil layer (42.25% and 14.77%, respectively) were higher than in the 10‐ to 30‐cm soil layer (28.77% and 11.60%, respectively). Moreover, the C v of STN was higher than that of STP. (2) The maximum C 0/(C 0 + C 1) of STN and STP in the soil layers was less than 25%, this indicated that a strong spatial distribution autocorrelation existed for STN and STP; and the STP showed higher intensity and more stable variation than the STN. (3) From the correlation analysis, we concluded that the topographic indexes such as elevation and slope direction all influenced the spatial distribution of STN and STP (correlation coefficients were 0.688 and 0.518, respectively). (4) The overall distribution of STN and STP in the Yaoxiang watershed decreased from the northwest to the southeast. This variation trend was similar to the watershed DEM trend and was significantly influenced by vegetation and topographic factors. These results revealed the spatial heterogeneity distribution of STN and STP, and addressed the influences of forest vegetation coverage, elevation, and other topographic factors on the spatial distribution of STN and STP at the watershed scale.


O R I G I N A L R E S E A R C H
Spatial heterogeneity distribution of soil total nitrogen and total phosphorus in the Yaoxiang watershed in a hilly area of northern China based on geographic information system and geostatistics Yu Liu 1 | Peng Gao 1 | Liyong Zhang 1 | Xiang Niu 2 | Bing Wang 2

| INTRODUCTION
As the basic building blocks of all known forms of life on Earth, soil total nitrogen (STN) and total phosphorus (STP) are significant evaluation indexes of soil fertility in the terrestrial ecosystem (Gao, Yang, & Liu, 2015;Zhang, Cha, & Shen, 2011). Fully understanding the characteristics of their spatial heterogeneity distribution at the watershed scale is the basis of reasonable evaluation to soil nutrients. In addition to providing nutrient sources for vegetation growth, they also affect the formation of the soil structure, soil biodiversity, and soil physical stability of the resistance to erosion, and so on (Kay & Rainer, 2008).
Soil nitrogen and phosphorus are important parts of the geobiochemical circulation. The unreasonable distribution of soil nitrogen and phosphorus is an important factor causing eutrophication in permanent wetlands such as rivers and lakes (Li et al., 2008). The enrichment nutrients in water body such as N and P, lead to some algae proliferate abnormally. As the results of transparency of water and dissolved oxygen decreasing, a large number of aquatic plants died. With the problems of nonpoint source pollution and eutrophication becoming increasingly serious, research on the spatial heterogeneity distribution features of STN and STP at the watershed scale has became a focus in the current environmental and agroforestry field (Dobos, Micheli, & Montanarella, 2007). In recent years, the employment of geographic information system (GIS) and geostatistics to study the spatial distribution characteristics and to reveal the spatial variation laws and its influencing factors has garnered more attention.
Researchers have conducted many studies about the spatial distribution characteristics of soil nutrients, and the results showed that the spatial distribution of STN and STP exhibited random or structured spatial variation characteristic because of the significantly different soil physical, chemical, and biological processes in different directions (Wang, 1999;Wang, Zhang, Yu, & Zhang, 2006). The research scope involved different scales. Initially, geostatistics was applied to study the spatial heterogeneity distribution and characteristics of soil nutrients by researches. The spatial correlation of soil nutrients in Hawaii was analyzed by geostatistics, and the results showed that the spatial distance of soil phosphorus, potassium, calcium, and magnesium content was between 32 and 42 km (Yost, Uehara, & Fox, 1982). In 1997, geostatistical method was used to study the soil spatial changes in the mid-west Taiwan (Chien, Lee, Guo, & Houng, 1997). Soil zinc map of the USA was analyzed using geostatistics, GISs, which concluded that the correlation distance of total zinc was 480 km (White, Welch, & Norvell, 1997). One hundred and six soil samples were collected and analyzed for soil organic matter, STN, available N, available P and available K, and the spatial variation in these soil nutrients was quantitatively analyzed by geostatistical methods, and the contour maps of each nutrient were drawn (Zhang, 1998). The spatial pattern characteristics of soil carbon and nitrogen in the turkey lakes watershed were studied employing the regression analysis (Creed, Trick, Band, & Morrison, 2002). In 2005, many studies based on the geostatistics were conducted to study the soil spatial variations in the Senegal Valley and found clear spatial variations in soil nutrients (Haefelel & Wopereis, 2005). Gradually, the GIS was used to predict the site indices in the studies of the spatial distribution of different soil nutrients. The spatial heterogeneity distribution of soil nutrients was researched employing GIS and geostatistics, and the results indicated that different soil nutrients had obvious different spatial variation patterns and different spatial correlations (Cambardella et al., 1994). In 2006, the digital elevation model (DEM) was employed to evaluate the spatial heterogeneity distribution of soil carbon and nitrogen in a Japanese cedar forest, and the researchers found that using DEM greatly facilitated the research on soil carbon and nitrogen (Kotaro, 2006). The GIS based on geostatistical method was applied to analyze the spatial variance characteristics of STP in Chaohu lake watershed and the contribution factors. The results showed that the semivariogram model of STP followed exponential model and the STP had a moderate spatial autocorrelation (Zhou, Gao, Sun, Zhao, & Zhang, 2007). In the study of soil nutrients of forest in Badaling, the results showed that different vegetation types corresponded to different spatial distributions of the soil nutrients (Yu, Zhang, & Zhu, 2009). The spatial distribution of soil organic carbon (SOC) in the Yaoxiang small watershed diminished gradually from the eastern and northwestern parts to central and southern parts, which was basically consistent with the trend of DEM in the Yaoxiang small watershed (Zhang, Gao, Wang, Liu, & Li, 2015). In brief, scholars have focused on the spatial distri- The study area belongs to a typical mountainous area in northern China, with serious soil and water erosion. Since the 1980s, with the implementation of the Return Farmland to Forests Project, many afforestation vegetation types have been developed, and many researches have been written on ecological afforestation technologies and soil hydrological benefits (Wan, Ji, & Gu, 2014;Yang, Hu, & Luo, 2011).
However, there are generally few reports on the soil improvement benefits of the Return Farmland to Forests Project, and it is not beneficial to assess the ecological effects of the Return Farmland to Forests Project.
In this paper, we selected the Yaoxiang watershed in southern Shandong Province as a research site. Using GIS and geostatistics, our goal was to achieve the following objectives: (1) to analyze the influences of slope, slope direction, and elevation on STN and STP, and establish the regression equation between elevation and the content of STN and STP;(2) to discuss the influences of forest vegetation coverage, elevation, and other topographic factors on the spatial heterogeneity distribution; and (3) to reveal the spatial heterogeneity distribution of STN and STP at the watershed scale using Kriging interpolation. We hope this can provide the theoretical basis for the construction of the Return Farmland to Forests Project and the evaluation of the forest ecological service function in the mountainous lands of northern China.

| Study area condition
The Yaoxiang watershed is located in Tai

| Sampling and processing
This work was based on the "Observation Methodology for Long-term Forest Ecosystem Research," which was the Forestry Standards of the People's Republic of China (LY/T 1952-2011) (Wang, Lu, Li, & Niu, 2011). The sampling plots were designed based on the topographic features and vegetation types (Table 1) of the Yaoxiang watershed, according to a DEM grid graph generated by a 1:1,000 topographic map (10 × 10 m) and combined with land-use data interpreted from remote sensing images in the Yaoxiang watershed. In November 2014, a total of 77 sampling plots were selected in the study area (Fig. 2), and a total of 154 soil samples were collected from the depth of 0-10 and 10-30 cm.
Meanwhile, GPS was used to record the longitude, latitude, and elevation values. After collecting the samples, they were taken to the laboratory by the processes of air drying, fine grinding, and sieving (2.000-, 0.250-, and 0.149-mm soil sieve). STN was measured using the semimicro-Kjeldahl determination (Li, Chen, & Zhao, 2006;Xie et al., 2012), and STP was measured using fused sodium hydroxide with the molybdenum stibium antireagent color method (Nelson & Sommers, 1982).    The formula of the semivariance function was calculated with equation (1) (Jiang, Liang, & Li, 2005). The semivariance function has five theoretical models and four indicators to describe the spatial correlation of variables, including the exponential model, gaussian model, spherical model, pure nugget effect model, and linear model. C 0 is the nugget, C 1 is the partial sill,

| Spatial heterogeneity analysis
C 0 + C 1 is the sill, and C 0 /(C 0 + C 1 ) is the ratio of the nugget to the sill.
Along with the range, they are the parameters that characterize the spatial structure of a soil property.
With the increase in the distance interval h, the semivariance increases from a nonzero value to a constant value. The value of C 0 / (C 0 + C 1 ) indicates the degree of variation in the soil. The higher the value of C 0 /(C 0 + C 1 ) is, the higher degree of variation that is caused by the random part; on the contrary, the variation is mainly caused by the structural part (Peter & Andreas, 2015;Stacey, Lark, Whitmore, & Milne, 2006). If the value of C 0 /(C 0 + C 1 ) is smaller than 25%, the variable has intense spatial correlation, between 25% and 75%, the variable shows medium spatial correlation; if the value is larger than 75%, the spatial correlation is weak (Li et al., 2008;Mackey et al., 2000;Zhao, Liu, Sui, Zhang, & Meng, 2006). The range indicates the spatial dependent distance; when the distance of the samples is close to it, the variable has spatial correlation; on the contrary, the spatial correlation does not exist between the samples if the distance exceeds the range (Liu, Li, Yang, & Xun, 2012;Moore, Gessler, Nielsen, & Peterson, 1993).
The soil fractal dimension (D): D is used to calculate the complexity of the variable factors (Li, Lei, & Wang, 2000). There is linear relation- where D is the soil fractal dimension indicating the degree of structural heterogeneity between soil samples. The larger the value of D is, the smaller the degree of heterogeneity between the soil samples and the greater the degree of homogeneity. Otherwise, the degree of heterogeneity between samples is great and the homogeneity is low.

| Conventional statistical analysis of STN and STP in the Yaoxiang watershed
According to the statistical results of STN and STP in the Yaoxiang watershed (Table 2), Table 2 showed that the STN and STP contents in the 0-to 10-cm soil layer were greater than those in the 10-to 30-cm soil layer. The soil variation coefficient (C v ) of STN and STP (42.25% and 14.77%, respectively) in the 0-to 10-cm soil layer was greater than those in the 10-to 30-cm soil layer (28.77% and 11.60%, respectively). The analysis of variance indicated that the STN content was significantly different between the 0-to 10cm and 10-to 30-cm soil layers (p < .01), and the STP content was different between the 0-to 10-cm and 10-to 30-cm soil layers (p < .05).

| Semivariance analysis of STN and STP in the Yaoxiang watershed
The choice of semivariance model is the key to analyze the spatial heterogeneity distribution of STN and STP. By comparing different parameters, the spherical model and exponential model were selected as the perfect models of STN and STP, respectively. From the results of the semivariance model analysis (Table 3, Fig. 3), it could be concluded that the maximum C 0 /(C 0 + C 1 ) of STN and STP in the 0-to 10-cm and 10-to 30-cm soil layers was 23.2%. Moreover, the values of C 0 /(C 0 + C 1 ) of STN and STP in the 0-to 10-cm soil layer were all higher than those in the 10-to 30-cm soil layer, and the value of C 0 /(C 0 + C 1 ) of STN was higher than that of STP. These results indicated that there was strong spatial distribution autocorrelation in STN and STP, and less spatial correlation existed in the surface soil (0-10 cm) than in the 10-to 30-cm soil layer, and STP showed more intense spatial correlation than STN. So it was consistent with the conventional statistical analysis of STN and STP above ( Table 2). In addition, the step length of the semivariance T A B L E 2 Conventional statistical analysis of soil total nitrogen (STN) and total phosphorus (STP) in Yaoxiang watershed Coefficient of variation (C v ) is the ratio of standard deviation to average, capital letters indicate the significant differences between the 0-to 10-cm and 10-to 30-cm soil layer (p < .01); lower case letters mean significance between the 0-to 10-cm and 10-to 30-cm soil layer (p < .05). The same below. analysis in our study was 50 m, and the range (159-239 m) was much greater than the step length. Therefore, the results obtained from the ordinary kriging interpolation were relatively accurate.

Soil layer(cm) Soil nutrient Mean (g/kg) Minimum (g/kg) Maximum (g/kg) Standard deviation C v (%)
Moreover, the determination coefficients (R 2 = .874-.912) of the semivariance model analysis of STN and STP were large, which demonstrated that the fitting results of the semivariance theory model were very ideal.

| Correlation analysis between topographic indexes and STN and STP contents in the Yaoxiang watershed
From the Table 4 and Fig. 4 Using multiple linear regression method, topographic indexes were used to explain and predict soil nutrient spatial variation. The regression equations are as follows: where Y STN is the STN content (g/kg) and Y STP is the STP content (g/ kg); X 1 is the elevation, X 2 is the sine slope, and X 3 is the cosine slope.
From the regression equation, we could know that the topographic indexes, such as elevation, cosine slope, and sine slope, were all related with the spatial distribution of STN and STP (fitting decision (3) Y STN = 2.238X 1 + 2.055X 2 + 1.968X 3 + 7.851, r 2 = .688 * * (4) Y STP = 1.183X 1 + 1.002X 2 + 0.577X 3 + 4.354, r 2 = .518 * * T A B L E 3 Semivariogram analysis of soil total nitrogen (STN) and total phosphorus (STP) in Yaoxiang watershed by ordinary kriging method C 0 is the nugget, C 1 is the partial sill, C 0 + C 1 is the sill, and C 0 /(C 0 + C 1 ) is the ratio of the nugget to the sill. R 2 is the determination coefficients.

| Soil fractal dimension (D) of STN and STP in the Yaoxiang watershed
From the full-range view, Table 5    The values of D indicate the spatial heterogeneity distribution degree in the different section. S-N, south-north section; NE-SW, northeastsouthwest section; E-W, east-west section; NW-SE, northwest-southeast section.

| Spatial distribution of STN and STP in the Yaoxiang watershed
the Yaoxiang watershed. Moreover, STN had high correlation with STP, and the correlation coefficient (r 2 ) was .711. The regression analysis result reflected a significant linear relationship between STN and STP, and the regression equation is: where Y STN is the STN content (g/kg) and X STP is the STP content (g/kg).

| Variation characteristics of STN and STP in different soil layers in the Yaoxiang watershed
There have been many studies on the distribution of soil nutrients in different soil depths (Hengl, Heuvelink, & Stein, 2004;Kay & Rainer, 2008;Kong, Zhang, Kou, Liu, & He, 2014). Studies about the SOC, STN, and STP distribution of typical vegetation riparian zones in the upper reaches of the Hun River have been conducted by Kong et al. (2014), and they showed that the SOC, STN, and STP contents decreased with the increase in soil layer depth. For the spatial variation in soil, the spatial heterogeneity distribution of SOC and other nutrients in the forestlands of the mountainous area of Changbai Mountain was studied, and the results showed that the SOC content varied with the soil depth and that the spatial variability of soil on the surface was greater than that in deep soil (Liu, 2013). In our study, we found that the STN and STP contents in the Yaoxiang watershed varied with the soil layer, and the STN and STP contents in the 0-to 10-cm soil layer were greater than those in the 10-to 30-cm soil layer. This result can be potentially explained by the large amount of forest litter in the Yaoxiang watershed. Litter is an important material base of the circulation of matter and the main supply of nitrogen and phosphorus in the ecosystem. The 0-to 10-cm soil layer covered more litter that falling from the forest vegetation, which were decomposed in the soil and releasing more nitrogen and phosphorus nutrients. They had beneficial effects on the soil and increased the source of nutrients in the soil, which leads to the conclusion that STN and STP contents in the 0-to 10-cm soil layer were greater than those in the 10-to 30-cm and had a significant effect. Moreover, the C v of STN was greater than that of STP. In general, when C v < 10%, there is weak variability; when C v is between 10% and 100%, there is moderate variability; and when C v > 100%, there is strong variability (Li et al., 2008). The maximum C 0 /(C 0 + C 1 ) of STN and STP in the soil layers was 23.2%; therefore, STN and STP had moderate variability in the Yaoxiang watershed. Our study also found that the C v of STN and STP in the 0-to 10-cm soil layer was greater than that in the 10-to 30-cm soil layer. This could be potentially explained by the effects of environmental factors such as topographic and human factors, which greatly damaged the surface soil structure (0-10 cm). However, the 10-to 30-cm soil layer experienced less interference by artificial factors; therefore, the value of C v of soil in the 10-to 30-cm soil layer was smaller.

| Spatial distribution analysis and the factors influencing STN and STP in the Yaoxiang watershed
Soil nutrients, such as nitrogen and phosphorus, are comprehensively affected by natural and human factors such as topography, (5) Y STN = 3.569 × X STP − 0.140, r 2 = .711 * * F I G U R E 5 Distribution map of soil total nitrogen (A, B) and soil total phosphorus contents (C, D) of different soil layers in Yaoxiang watershed by ordinary kriging method vegetation, climate, land use and management measures. Various factors jointly determine the spatial distribution pattern and the formation of the soil, and they also determine the transformation direction and rate of change of soil nutrients . In the research on the nutrient characteristics of forest soil at different elevations in the Fanjingshan Nature Reserve, researchers found that the nutrients showed clear vertical distribution characteristics with the increase in altitude (Yan, Zhang, Shi, Lin, & He, 2015). Through the research about the spatial variability of soil nitrogen and the affecting factors in the hilly area of the Mid-Sichuan Basin, it was found that the spatial distribution of STN showed similar trend with STP; both presented a low trend from north to south as the topography changed (Luo et al., 2015). In our research, the spatial distribution tendency of STN and STP decreased from the northwest to the southeast, and it was consistent with the DEM map (Fig. 2) in the Yaoxiang watershed. Therefore, the spatial distribution of STN and STP had close correlation with the elevation in the Yaoxiang watershed.
Moreover, because the maximum C 0 /(C 0 + C 1 ) of STN and STP in the soil layers is less than 25%, this indicates that a strong spatial distribution autocorrelation existed in STN and STP in the Yaoxiang watershed. Structural factors such as topography and vegetation were the main factors influencing the spatial heterogeneity distribution of STN and STP, whereas the influence of random factors was smaller. The values of C 0 /(C 0 + C 1 ) of STN and STP in the 0-to 10-cm soil layer were all higher than those in the 10-to 30-cm soil layer, indicating that there was less spatial correlation in the surface soil (0-10 cm) than in the 10-to 30-cm soil layer. From Table 3, it can also be observed that the value of C 0 /(C 0 + C 1 ) of STN was higher than that of STP. Therefore, STP showed more intense spatial correlation than STN. The spatial heterogeneity distribution of STP was mainly affected by structural factors, and its change was relatively stable (Peter & Andreas, 2015;Stacey et al., 2006 Fig. 1; Table 1). Therefore, northwest had high vegetation coverage, suitable temperature, and developed humidity plant roots. Litter decomposition in the northwest was beneficial and conducive for the accumulation of soil nitrogen and phosphorus (Liang, Yuan, & Lin, 2006;Liu, Fu, & Wu, 2003). Meanwhile, the southeast part of the watershed was mainly covered with C. mollissima and R. pseudoacacia and was disturbed by human activities in this area, resulting in much lower STN and STP contents than in the northwest part.
Thus, we should develop comprehensive measures to reduce human interference and destruction to the forest vegetation and strengthen its preservation and restoration; it is very important for the storage and accumulation of nitrogen and phosphorus in the mountain watershed.

| CONCLUSIONS
1. The STN and STP contents in different soil layers of the Yaoxiang watershed were significantly different, and they were all of medium variation. In particular, the variation in STN was greater than that of STP. Restoration.

CONFLICT OF INTEREST
The authors have no conflict of interests to declare.

FUNDING INFORMATION
The major water conservancy scientific research and technology promotion projects of Shandong province (2014-2016),