Spatial Distribution of Soil Nitrogen, Phosphorus and Potassium Stocks in Moso Bamboo Forests in Subtropical China

Abstract

Moso bamboo is famous for fast growth and biomass accumulation, as well as high annual output for timber and bamboo shoots. These high outputs require high nutrient inputs to maintain and improve stand productivity. Soil nitrogen (N), phosphorus (P), and potassium (K) are important macronutrients for plant growth and productivity. Due to high variability of soils, analysing spatial patterns of soil N, P, and K stocks is necessary for scientific nutrient management of Moso bamboo forests. In this study, soils were sampled from 138 locations across Yong’an City and ordinary kriging was applied for spatial interpolation of soil N, P, and K stocks within 0–60 cm. The nugget-to-sill ratio suggested a strong spatial dependence for soil N stock and a moderate spatial dependence for soil P and K stocks, indicating that soil N stock was mainly controlled by intrinsic factors while soil P and K stocks were controlled by both intrinsic and extrinsic factors. Different spatial patterns were observed for soil N, P, and K stocks across the study area, indicating that fertilizations with different ratios of N:P:K should be applied for different sites to maintain and improve stand productivity. The total soil N, P, and K stocks within 0–60 cm were 0.624, 0.020, and 0.583 Tg, respectively, indicating soils were important pools for N, P, and K.

1. Introduction

Soil nitrogen (N), phosphorus (P), and potassium (K) are important macronutrients which can limit or co-limit plant growth [1,2]. Human activities, such as fertilization, reclamation, and weeding, have greatly affected the biogeochemical cycling of N, P, and K, thereby altering the pattern, magnitude, and extent of nutrient limitation on land [3]. Effective, efficient, and site-specific management and estimation of soil N, P, and K have attracted great interests for scientists looking to improve nutrient input efficiency, thus increasing stand productivity and reducing environmental risks [4]. Due to various climatic conditions [5], parent materials [6], topography [7], vegetation types [8], soil texture [9], and land use [10], soils are characterized by a highly spatial and temporal variability. This has made the accurate estimation of spatial nutrient content difficult. However, rational soil management requires a deep understanding of the spatiotemporal variability of soil N, P, K, and their maps. For example, areas of particular concern, such as nutrient deficiency, can be identified. Additionally, soil N and P levels are closely related to soil organic carbon cycling [11], which may lead to dynamic effects on greenhouse gas emissions that potentially result in feedback to global climate change. Therefore, it is necessary to improve our understanding of the spatial distribution of soil N, P, and K stocks when evaluating current and potential soil productivity, identifying potential environmental protections, and assessing climate change [12].

Fast developments of geostatistical techniques, such as kriging [11], have significantly advanced the estimation of spatial variability of soil N, P, and K stocks. The most common interpolation techniques calculate the unmeasured property at a given place using measured neighbours data with a weighted mean [13]. The core of geostatistics is that experimental semivariograms are calculated to analyse the spatial autocorrelation of the edaphic variables and determine the range of spatial dependence [13]. The spatial data structure, choice of variogram models, search radius, and number of closest neighbouring points can determine the performance of Kriging approaches. Compared to other interpolation approaches (e.g., Inverse Distance Weighting and splines methods [14,15]), ordinary kriging (OK) can provide more accurate estimates of the spatial distribution of soil nutrients because: (1) OK can provide the best linear unbiased estimates and information on the distribution of the estimation error [15]; (2) OK has relatively strong statistical advantages [16]. Additionally, OK is easy to conduct with high accuracy compared to other kriging approaches, such as cokriging, universal kriging, and factorial kriging. For example, cokriging requires an additional correlated covariant, which leads to the substantial increase of field and lab work and does not increase the statistical power [15]. Thus, OK has been widely used in spatial interpolation for different soil properties. For example, Liu et al. [11] predicted the spatial patterns of soil N and P stocks for the whole Loess Plateau using OK, and total N and P stocks amounted to 0.217 and 0.205 Pg (1 Pg = 10¹⁵ g), respectively. Martín et al. [17] predicted the spatial distribution of soil carbon stocks within 0–30 cm with OK across Spain.

Bamboo is an important forest type in Southern China, representing an area of 6.16 million ha, more than 70% of which are Moso bamboo (Phyllostachys heterocycla (Carr.) Mitford cv. Pubescens) forests. Moso bamboo forests are famous for their rapid growth and fast biomass accumulation. However, annual timber harvest and bamboo shoots remove a large amount of nutrients, which means that Moso bamboo forests require more nutrient input compared to other forest types in order to maintain high stand productivity. Consequently, intensive management with fertilization, reclamation, and regular understory removal is becoming increasingly popular in Southern China, especially in the main bamboo producing provinces, such as Zhejiang and Fujian Provinces [18,19,20]. These activities have significantly changed the nutrient levels by affecting the microbial processes, soil structures, and chemical compositions [21,22]. Although many studies focus on the spatial patterns of soil N, P, and K stocks in different ecosystems [23,24], site specific maps of soil N, P, and K stocks for scientific management of Moso bamboo forests are still lacking in our study area. In contrast, a better understanding of the spatial distribution of soil N, P, and K stocks is strongly required for land management and for maintaining and improving stand productivity of Moso bamboo forests. Therefore, the objectives of this study were to: (1) predict the spatial variability of soil N, P, and K stocks of Moso bamboo forests in Yong’an City, China and (2) estimate total soil N, P, and K stocks within 0–60 cm.

2. Materials and Methods

2.1. Study Area

The study was conducted in the Moso bamboo forests (Figure 1, polygons) across Yong’an City, Fujian Province, Southern China (117°56′–117°47′ E, 25°33′–26°12′ N). It has a subtropical southeast monsoon climate with a mean annual temperature of 19.3 °C, while the lowest temperature is −11°C and the highest temperature is 40 °C. The average annual precipitation is about 1600 mm [25,26]. The elevation of the study area ranged from 580 m to 1605 m above sea level [25]. The accumulated temperature of ≥10 °C is 4520–5800 °C, lasting for 225–250 days, and the relative humidity is about 80% [25]. The main soil type is Oxisol. Yong’an City has a forest cover of 82%, including an area of 5.85 × 10⁴ ha of Moso bamboo forests [26]. Moso bamboo forests are mainly distributed below 800 m, most of which are pure stands and sometimes mixed with Keteleeria cyclolepis Flous., Cunninghamia lanceolata (Lamb.) Hook., Myrica rubra (Lour.) S. et Zucc., Choerospondias axillaris (Roxb.) Burtt et Hill., Liriodendron chinense (Roxb.) Burtt et Hill., and Schima Superba Gardn. et Champ., etc. To improve the stand productivity and increase the income, fertilization treatment has been widely applied in most of the Moso bamboo forests. The fertilization treatment was conducted across the whole study area using the same protocol established by the local Forest Bureau. Similar fertilization types—organic fertilizers with N, P, and K addition, were offered every year by the Forest Bureau because the residents can get an economic subsidy if they follow the protocol and use the recommended fertilizers.

2.2. Soil Sampling

The representative soils were sampled in the sub-compartments of forest resource management sites in Fujian province, China. These sample sites are established by the local Forest Bureau for soil mapping according to the protocol of forest management inventory [27] (Figure 1). In the targeted sub-compartments, a total of 138 sites were selected and a cluster of three circular plots with a size of 33.3 m² were established for each site. Due to the unfavourable access of sampling plots, most of the samples were collected in the East and Southern portions of the study area, which covered 78% of the total Moso bamboo forests. Thus, it is acceptable with these plots to study the spatial distribution of soil N, P, and K stock across the whole study area using the unbiased interpolation approach of OK (see below).

In each plot centre, soils were sampled down to 60 cm by the consideration of three layers: 0–20 cm, 20–40 cm, and 40–60 cm. In the field, soil samples in the sample layer were mixed for one cluster of the three circular plots. In the laboratory, soil samples were air-dried at room temperature and prepared for sieving through 2-mm and 0.15-mm for N, P, and K content analysis. Identifiable plant residues and root materials were removed during sieving. Because the majority of bamboo roots were distributed within 0–40 cm [28], soil samples down to 60 cm were able to meet the research purpose to study spatial patterns of soil N, P, and K stocks. Soil bulk density was determined by a cutting ring approach for each soil layer [29]. During the implementation of fertilization, stones and rocks were removed in Moso bamboo forests. As a result, few stones and rocks were observed in core samples, thus, we did not correct for gravel content. Additionally, site information about elevation, coordinate, soil depth, soil type, soil organic matter content, and bamboo diameters were recorded and determined.

Soil N content (g kg⁻¹) was analysed using the Kjeldahl digestion procedure (5 mL concentrated H₂SO₄, 1 g of soil passed through a 0.15-mm sieve, heated for 90 min, titrated with 0.2 mol L⁻¹ standardized Na₂B₄O₇). Soil total phosphorus concentration (g kg⁻¹) was determined by alkaline digestion (NaOH 0.2 g, <0.25 mm soil 0.25 g, 300 °C for 15 min and then 750 °C for 15 min) followed by molybdate colorimetric measurement. Soil K content (g kg⁻¹) was determined by sodium hydroxide flame photometer (NaOH 0.2 g, <0.25 mm soil 0.25 g, 300 °C for 15 min and then 750 °C for 15 min, 80 °C water 10 mL). All the analyses followed the national standard protocol of China [29].

Thus, soil N, P, and K stocks were calculated as [30,31,32]:

$$SON\left(\mathrm{or}\mathrm{P}\mathrm{or}\mathrm{K}\right)\mathrm{stock}=SON\left(\mathrm{P},\mathrm{K}\right)\times BD\times D/10$$

(1)

where SON (or P or K) are the soil N, P, and K contents (g kg⁻¹); BD is bulk density (g cm⁻³); D is the depth of soil layer (cm). Total soil N, P, and K stocks within 0–60 cm were summed for all soil layers.

2.3. Statistical and Geostatistical Analyses

A standard statistical analysis (mean, standard deviation, first quartile, third quartile, etc.) was conducted in R to illustrate trends of the soil N, P, and K stocks [33]. The coefficient of variation (CV) was used to describe the degree of general variation. To meet the assumption of normality for geostatistical analysis, the raw data was log-transformed and transformed back by weighted mean in GS + 10.0 [34]. OK was performed in GS + 10.0 and the distribution maps of soil N, P, and K stocks were produced in ArcGIS 10.2 [35].

2.4. OK

In the last several decades, many spatial interpolation approaches have been developed, such as kriging; OK normally performed better and is easier to apply compared to other approaches, thereby it has been widely used in the spatial interpolation of soil properties [11,17,23]. In the OK process, semivariograms are used to describe the spatial autocorrelation and provide parameters for optimal spatial interpolation based on the theory of regionalized variables [36]. Well-known theoretical models (e.g., exponential, Gaussian, and spherical models) are commonly used to fit the experimental semivariograms [37]. The experimental semivariograms are expressed as a function of the distance between the sampled points and calculate the integrity of spatial continuity in one or multiple directions using the following expression [38]:

$$\gamma \left(h\right)=\frac{1}{2N\left(h\right)}{\displaystyle \sum }_{i=1}^{N\left(h\right)}{[z\left(x_i\right)-z\left(x_{i+h}\right)]}^2$$

(2)

where z (x_i), and z (x_{i + h}) are values of z at locations x_i and x_{i +h}, respectively; h is the lag, and N(h) is the number of pairs of sample points separated by h. In this study, spherical, exponential, linear and Gaussian models were used to attempt to describe the semivariograms of soil N, P, and K stocks within 0–60 cm, then the best models determined by the smallest residuals and the largest determination coefficients were selected for spatial interpolation.

In the semevariograms, three major parameters were derived—the nugget (C₀), the sill (C + C₀), and the range (A₀)—in order to identify the spatial structure of soil N, P, and K stocks at a given scale. The sill (C + C₀) represents the total variation and the ratio of nugget and sill is considered as a criterion to classify the spatial dependence [17]. The range (A₀) represents the separation distance when semivariogram is stabilized, beyond which the measured data are spatially independent [39]. More details about semivariograms and the kriging approach can be found in Goovaerts [37]. In our study, OK was applied for the spatial patterns of soil N, P, and K stocks. The most likely value which could be expected in a particular grid cell when using m neighbouring observations was defined as [17]:

$$R\left(x\right)={\displaystyle \sum }_{j=1}^m{\mathsf{\delta }}_{j}z\left(x_i\right)$$

(3)

where ${\delta }_j$ means the optimal weight under the condition of ${\displaystyle \sum }{\mathsf{\delta }}_j=1$, m = 16 in this study.

2.5. Model Validation

To evaluate the prediction accuracy of soil N, P, and K stocks, leave-one-out cross-validation was used [11,17,40]. In the validation process, one datum was omitted and this value was estimated by the remaining data. Thereafter, the estimated value was compared with the real value of the omitted point [17]. This process was repeated for all the observations. Four commonly used indices (i.e., absolute mean error (AME), mean error (ME), root mean square error (RMSE), and determination coefficient (R²)) were used to compare the interpolation accuracy. These indices were calculated as follows [11,15]:

$$\mathrm{AME}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left|\left(P_i-M_i\right)\right|$$

(4)

$$\mathrm{ME}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left(P_i-M_i\right)$$

(5)

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(P_i-M_i\right)}^2}$$

(6)

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(P_i-M_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(M_i-\overline{M}\right)}^2}$$

(7)

where P_i, M_i, and $\overline{M}$ are predicted values, measured values, and the mean values of the measured data, respectively.

3. Results and Discussion

3.1. Descriptive Statistics

The summary of statistics of soil N, P, and K stocks within 0–60 cm is shown in

Table 1. Summary of statistics of soil N, P, and K stocks within 0–60 cm (Mg ha⁻¹).

Nutrient	Mean	Minimum	Maximum	Median	SD	CV (%)	1st Qu	3rd Qu	Skewness	Kurtosis	p of S-W
N Stock	11.57	5.30	20.20	10.98	3.23	27.92	8.81	13.66	0.42	−0.33	0.053
P Stock	0.38	0.18	0.68	0.36	0.11	29.10	0.30	0.45	0.60	−0.14	0.010
K Stock	10.80	5.73	18.72	10.81	3.15	29.17	8.11	13.03	0.32	−0.72	0.014

SD = standard deviation; CV = coefficient of variance; 1st Qu = 25% quartile; 3rd = 75% quartile; S-W test = Shapiro-Wilk test.

. N stock in Moso bamboo forest increased from a minimum of 5.30 Mg ha⁻¹ to a maximum of 20.20 Mg ha⁻¹ with an average of 11.57 Mg ha⁻¹. These values fell within the range of N stocks in different forest types across China (1.0–49.9 Mg ha⁻¹) [41]. However, the mean N stock was even higher than the average N stock of China’s forests within one meter (8.4 Mg ha⁻¹) [41] and the Moso bamboo in Jiangxi Province (2.82 Mg ha⁻¹ within 60 cm) [42]. These differences were mainly attributed to the fertilization application in Moso bamboo forests in the study area, whereas no stand management occurs in most of China’s forests.

Soil P stock ranged from 0.18 to 0.68 Mg ha⁻¹ and the mean value was 0.38 Mg ha⁻¹. Although this range lied within the reported values of main Chinese soil types, the mean value was much lower than the average across China (6 Mg ha⁻¹ within 60 cm) [43]. Similarly, average soil K was 10.80 Mg ha⁻¹ and the lowest soil K stock was 5.73 Mg ha⁻¹ and highest soil K stock was 18.72 Mg ha⁻¹, which were lower than soil K stocks in subtropical forests in China [44]. These low values of soil P and K stocks may be related to the high output of P and K due to the high output of timber and bamboo shoots every year. In Moso bamboo forests, P is the most limited nutrient for stand productivity [4]. Although fertilization was applied in Moso bamboo forests, P and K input could not compensate for the P and K output from the stand. This result further indicated that more P and K fertilization should be applied in Moso bamboo forest in order to maintain and improve stand productivity.

A CV of 10% indicates a low variability and 10%–90% indicates a moderate variability, and CV >90% indicates higher variability [45]. In this study, the CVs of soil N, P, and K stocks were 28%–29%, indicating a moderate variability of soil N, P, and K stocks within 0–60 cm. p values of Shapiro-Wilk test were 0.053 for N stock, 0.010 for P stock, and 0.014 for K stock, indicating a normal distribution of N stock and a non-normal distribution of P and K stocks at a significance level of 0.05. However, in order to make the comparison of spatial interpolation of soil N, P, and K stocks consistent, a natural log-transformation was conducted for the three stocks in order to meet the assumption of normal distribution.

3.2. Geostatistical Analysis of Soil N, P, and K Stocks

The semevariograms are presented in Figure 2 and their parameters are shown in

Table 2. Models and their parameters fitted for the semivariograms of soil N, P, and K stocks within 0–60 cm.

Nutrient	Models	Nugget (C0)	Sill (C0 + C)	Nugget/Sill	Range (A0, m)	Determination Coefficient	Residuals
N stock	Exponential	0.0111	0.0912	12.17	6100	0.642	0.00458
P stock	Exponential	0.0570	0.1340	42.54	30,570	0.741	0.00110
K stock	Spherical	0.0343	0.0986	34.79	25,700	0.643	0.00469

. The exponential model performed best for soil N and P stocks while the spherical model was best for soil K stock based on the smallest residuals and the highest determination coefficient. These models produced high determination coefficients, ranging from 0.64 to 0.74.

Nugget values present undetectable experimental errors, field variation within the minimum sampling space, and inherent variability [23]. In this study, nugget values were lowest for soil N stock (0.0111), and highest for soil P stock (0.0570). These positive values suggested a positive nugget effect, a sampling error, or random and inherent variability of soil N, P, and K stocks [23]. Sill values represent total spatial variation [11]. The sill values ranged from 0.0912 for N stock to 0.1340 for P stock.

The ratio of nugget-to-sill represents a spatial dependence. If the ratio is lower than 25%, it indicates a strong spatial dependence. If the ratio is higher than 75%, it indicates a weak spatial dependence, while a ratio between 25% and 75% indicates a moderate spatial dependence [46]. A strong spatial dependence is attributed to soil intrinsic properties, such as soil parent material, soil texture, topography, and vegetation [23,47]. A weak spatial dependence indicates that the spatial variability is mainly regulated by extrinsic variations, such as soil fertilization and cultivation practices [23,46]. Therefore, a moderate spatial dependence is controlled by both intrinsic and extrinsic factors. In this study, the nugget-to-sill ratio was 12% for soil N stock, indicating that soil N stock in the study area was mainly controlled by intrinsic factors and that external factors exerted little effects on soil N stock. This may be related to N fertilization application in Moso bamboo forests in that the N input from fertilization can meet the requirement of N output from timber and bamboo shoot harvest. This conclusion was also consistent with the above result (Table 1) that N stock in this study area was higher than other forest types. The nugget-to-sill ratios were 43% for soil P stock and 35% for soil K stock, indicating a moderate spatial dependence that was controlled by both intrinsic and extrinsic factors. This was evidenced by local stand management and complex topography. Intensive managements, such as fertilization and weeding, were widely conducted in Moso bamboo forests in order to improve timber and bamboo shoot output in the study area [19,25]. Meanwhile, the study area was characterized by a high variability of elevations varying from 580 to 1605 m above sea level [25], which could directly affect the soil processes and nutrient contents [11]. In the study area, intensive fertilization treatment was conducted every year, which was expected to be the dominant factor that led to the change of the spatial variability of soil N, P, and K stocks. However, our results (nugget-to-sill ratio) indicated that intrinsic factors played a more important role in controlling spatial patterns of soil N, P, and K stocks.

The ranges (A₀) indicate different influence zones of environmental factors at different scales [23]. Within the range, soil properties are not spatially independent, while beyond the range, soil properties were spatially independent [39]. In this study, the smallest range was found for soil N stock (6100 m), and increased to 25,700 m for soil K stock and 30,570 m for soil P stock. The results indicated that soil N stock was more heterogeneous compared to soil P and K stocks, which was associated with soil processes because soil N stock was mainly controlled by intrinsic factors (see above). Generally, the spatial range was larger than the sampling intervals (Figure 1), which suggested that the sampling strategy in this study was appropriate for studying spatial patterns of soil N, P, and K stocks. However, only a small number of sampling plots were located in the northwest of the study area due to a low cover of bamboo forests. In contrast, the majority of the sampling sites were located in the southern and eastern part of the study area, where most of the bamboo forests were distributed, thus the sampling strategy could meet our research purpose and could predict the spatial distributions of soil N, P, and K stocks accurately.

3.3. Cross-Validation of OK

The predicted values of OK were plotted against the measured values to evaluate the interpolation performance (Figure 3). The linear model and 1:1 dashed line intersected for soil N, P, and K stocks. Before the intersection, the linear model overestimated soil N, P, and K stocks, and vice versa. This conclusion supported the findings of many previous studies because of the nature of the algorithms of OK, which aimed to achieve unbiased predictions of mean values [11,15]. AME, ME, and RMSE were calculated and presented in

Table 3. Cross-validation indices for ordinary kriging (OK).

Nutrient	AME	ME	RMSE
N stock	2.1506	−0.2115	2.5557
P stock	0.0626	−0.0022	0.0797
K stock	2.0250	−0.1634	2.3942

AME = absolute mean error; ME = mean error; RMSE = root mean square error.

. The closer that the AME, ME, and RMSE values are to zero, the better the model performed. ME of soil N, P, and K stocks ranged from −0.0022 to −0.2115, which were close to zero, indicating that OK produced relatively unbiased values for spatial interpolation. The ME values were negative, illustrating that OK generally underestimated soil N, P, and K stocks.

Determination coefficients ranged from 0.37 to 0.47 (Figure 3), which may appear relatively low, but they were similar to many previous studies regarding the spatial interpolation of soil properties [47,48]. This problem is associated with the dataset, which was not sampled with a probabilistic design [48]. For example, the majority of sampling plots were allocated in the southern and eastern parts of Yong’an City in the current study (Figure 1). Therefore, a dataset with a probabilistic design is recommended for spatial interpolation of soil properties. Another possible explanation might be the strong local variation of soil N, P, and K stocks caused by the variability in environmental conditions and fertilization practices. On the other hand, the low correlation between the predicted and measured data indicates that a better methodology, such as one using Artificial Neural Network and Random Forest, should be developed to improve the accuracy of spatial interpolation of soil N, P, and K stocks, including consideration of both intrinsic and extrinsic factors.

3.4. Spatial Prediction of Soil N, ,P and K Stocks

The spatial patterns of soil N, P, and K stocks predicted by OK are shown in Figure 4. Predicted soil N, P, and K stocks ranged from 7.64 Mg ha⁻¹ to 16.02 Mg ha⁻¹, from 0.19 Mg ha⁻¹ to 0.63 Mg ha⁻¹, and from 6.85 Mg ha⁻¹ to 17.13 Mg ha⁻¹, respectively. The lowest soil N, P, and K stocks were found in the middle or southern part of study area, where the city centre was located with a high population, indicating that human activities led to a significant decrease of soil N, P, and K stocks. The highest soil N stock was observed in the east of the study area, the highest soil P stock was found in the southwest and north, and the highest soil K stock was in the north. However, plant growth was constrained by the most limited nutrient. These different spatial distribution patterns suggested that fertilizers with different N:P:K ratios should be applied to maintain and improve the stand productivity. For instance, a relatively lower ratio of P but higher K ratio fertilizer should be applied in the south of the study area (Figure 4b,c), while a relatively lower ratio of N fertilization should be implemented in the east (Figure 4a). Further studies are strongly encouraged to study the N:P:K ratio for the whole study area in order to manage Moso bamboo forests scientifically. The total soil N, P, and K stocks of the whole study area were 0.624, 0.020, and 0.583 Tg (1 Tg = 1 × 10¹² g), respectively, indicating soils were important pools of soil N, P, and K.

4. Conclusions

In this study, the current status and spatial interpolation of soil N, P, and K stocks within 0–60 cm were analysed using the measured data from 138 locations in Moso bamboo forests in Yong’an City, China. OK was applied for spatial interpolation of soil N, P, and K stocks across the whole study area. Exponential and spherical models performed well in describing the spatial distribution of soil N, P, and K stocks with determination coefficient from 0.64 to 0.74, and cross-validation demonstrated that OK performed well for the spatial interpolation of soil N, P, and K stocks. Soil N stocks showed a strong spatial dependence, indicating soil N stocks were mainly controlled by intrinsic factors. Soil P and K stocks showed a moderate spatial dependence, suggesting that soil P and K stocks were controlled by both intrinsic (e.g., soil parent material, soil texture, topography) and extrinsic (e.g., soil fertilization and cultivation practices) factors. Soil N, P, and K stocks showed different spatial patterns across the whole study area, indicating that fertilizers with different ratios of N:P:K should be applied for different sites.