Large-scale spatial variability of eight soil chemical properties within paddy fields

Catena 188 (2020) 104350

Contents lists available at ScienceDirect

Catena journal homepage: www.elsevier.com/locate/catena

Large-scale spatial variability of eight soil chemical properties within paddy fields Liangxia Duana, Zhenwei Lib, Hongxia Xiea, Zhiming Lic, Liang Zhanga, Qing Zhoua,

T

⁎

a

College of Resources & Environment, Hunan Agricultural University, Changsha 410128, China Key Laboratory for Agro-Ecological Processes in Subtropical Region, Institute of Subtropical Agriculture, Chinese Academy of Sciences, Changsha 410125, China c Hunan Station of Soil and Fertilizer, Changsha 410005, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Spatial variation Soil chemical properties Geostatistics Interpolation Rice fields

Rice is one of the most important crops in the world. It provides food for 40% of the world’s population and 60% of China’s population. Soil nutrients have a significant impact on both agriculture and the environment, particularly with regards to soil fertility, soil quality, and rice production. However, minimal research has been conducted to address spatial patterns of soil nutrients at a large scale within paddy fields. This information is crucial for improving not only soil nutrient management, but also rice yields. Soil surface samples (0–20 cm of plough depth, totalling 8890) were collected from paddy fields to determine the spatial variability of eight soil chemical properties (soil organic carbon (SOC), total nitrogen (TN), available N (AN), total phosphorus (TP), available P (AP), total potassium (TK), rapidly available K (RAK), and slowly available K (SAK)). Inverse distance weighting (IDW) and ordinary kriging methods were applied to produce a continuous soil nutrient surface. Results indicated that soil chemical properties vary substantially. Paddy fields were characterized by high mean concentrations of SOC (19.7 g kg−1), TN (1.91 g kg−1), and AN (164.7 mg kg−1), implying that no additional C and N fertilizer was needed in regions with high SOC and N. SAK and AP demonstrated moderate (36.9%) and weak (96.2%) spatial dependence, respectively, while strong spatial dependence (10.1–13.5%) was observed across the remaining six soil chemical properties. Regional distribution maps of soil chemical properties were produced and ordinary kriging methods interpolated more accurately than the IDW method.

1. Introduction Rice (Oryza sativa L.) is the staple food for 40% of the world’s population, and for > 60% of the total population in China (Bi et al., 2009; Zhang et al., 2010; Li et al., 2015a). China is the largest producer and consumer of rice in the world, producing approximately 28.7% of the world’s total rice. To meet the needs of an increasing population, this production must increase by 1.21% annually (Ministry of Agriculture, 2008). However, rice production decreased from 3.37% during 1970–1990, to 0.60% during 1990–2006, which may be attributed to soil nutrient imbalances (Zhang et al., 2010; Liu et al., 2014). Knowledge of soil nutrient levels and their spatial variability is important to sustain and enhance rice production, and ensure national food security in China. Within agricultural ecosystems, soil organic carbon (SOC), nitrogen (N), phosphorus (P), and potassium (K) are the major indicators and determinants of soil quality and fertility as they are strongly linked to plant growth and productivity (Hati et al., 2008; Yuan et al., 2014;

⁎

Guan et al., 2017). Generally, SOC content is recognized as a key component of soil fertility (Xu et al., 2011). For example, SOC was one of the most important indicators of productivity within low-input farming systems (Bruun et al., 2006). Organic amendments had a positive effect on rice production as it increased SOC and soil nutrient capacity (Bi et al., 2009). Therefore, maintaining a satisfactory SOC was a major determinant of the sustainability and productivity of terrestrial ecosystems. Nitrogen and P are among the most important nutrients for crop growth and they play a vital role in maintaining rice production (Spiertz, 2009; Yuan et al., 2014). To sustain or improve rice yields, N and P chemical fertilizers are widely applied to paddy fields to increase N and P levels (Liu et al., 2014). General overuse of N and P fertilizers poses potential environmental contamination issues since excess N and P can enter into surface and ground water bodies through leaching, runoff and soil erosion (Beman et al., 2005; Chen et al., 2008; Hou et al., 2015). In contrast, K fertilizer application is becoming increasingly inadequate in China (Tan et al., 2012). A report by the National

Corresponding author. E-mail addresses: [email protected] (Z. Li), [email protected] (H. Xie), [email protected] (Q. Zhou).

https://doi.org/10.1016/j.catena.2019.104350 Received 3 July 2019; Received in revised form 29 October 2019; Accepted 2 November 2019 0341-8162/ © 2019 Elsevier B.V. All rights reserved.

Catena 188 (2020) 104350

L. Duan, et al.

pattern was not practical as a sampling strategy. A dense soil sampling strategy of 8890 sampling sites throughout Hunan Province was selected to gain accurate soil nutrient data (Fig. 1c). The average sampling interval was approximately 1.42 km. At each sampling site, 20 soil samples were taken from the surface soil (0–20 cm, plough depth) within a roughly 200–800 m2 area in one rice field and then all samples were mixed to obtain a homogenous sample for analysis of soil chemical properties (soil organic carbon (SOC), total nitrogen (TN), available nitrogen (AN), total phosphorus (TP), available phosphorus (AP), total potassium (TK), rapidly available potassium (RAK), and slowly available potassium (SAK)), pH and cation exchange capacity (CEC) measurements. The SOC was measured using the potassium dichromate colorimetric method. TN was determined using the Kjeldahl method. AN was extracted with 1 mol L−1 KCl and measured according to the Alkali Diffusion method. TP was extracted from soil with 1 N HCl after ignition at 550 °C. AP was extracted using 0.5 mol L−1 NaHCO3 (pH = 8.5), and then the TP and AP concentrations were determined following the Mo-Sb colorimetric method. For TK determination, soil was digested in a nickel crucible with sodium hydroxide at 750◦C, and RAK, SAK were extracted with 1 mol L−1 ammonium acetate and nitric acid, respectively, and then TK, RAK, SAK were measured by using flame photometry. Soil pH was measured using the glass electrode method. Soil CEC was determined using the ammonium acetate method. Additional intact soil samples were collected using cylinders (5 cm long and 20 cm2 cross-section area) for each site to determine bulk density (BD) through the oven-drying method. Longitude, latitude, and elevation of each site were measured using GPS (5 m precision). Precipitation, accumulated temperature ≥0 °C (AT0), and accumulated temperature ≥10 °C (AT10) were derived from the China Administration of Meteorology.

Soil Survey Office (1998) showed that 40% of agricultural land in southern China was K deficient, while 70% of paddy soil was also considered K deficient (Yang et al., 2004). Potassium fertilizer does not improve K levels as quickly as N and P fertilizer application. Therefore, K deficiency is generally closely related to excess N and P application (Liu et al., 2009; Zhang et al., 2010). Knowing how to establish and maintain optimal fertilization application rates remains a challenge in paddy fields. Due to this, the identification of spatial patterns of SOC, N, P, and K is essential for farmers attempting to balance fertilizer application and increase fertilization efficiency for rice productivity. Geostatistics is an effective tool to characterize the spatial variability of soil nutrients (Webster and Oliver, 2007; Liu et al., 2009). Numerous studies have used geostatistics to analyze spatial variability and distribution of soil properties at various scales (Cambardella, 1994; Schöning et al., 2006; Huang et al., 2007; Wang et al., 2009; Liu et al., 2012a; Liu et al., 2013; Elbasiouny et al., 2014; Liu et al., 2014; He et al., 2015; Guan et al., 2017; Usowicz and Lipiec, 2017; Dai et al., 2018). Interpolation approaches such as inverse distance weighting (IDW) and ordinary kriging (OK) methods are generally used to predict spatial variability of soil properties (Liu et al., 2009; Li and Shao, 2014). For example, Liu et al. (2013) mapped the spatial distributions of soil total N and soil total P densities in surface soil layers across the Loess Plateau using the OK method. Usowicz and Lipiec (2017) used the IDW approach to generate 2D maps of soil properties and crop yields in a cultivated field. However, the IDW and OK methods usually produce considerable differences in landscape or regional estimates of soil properties (Guan et al., 2017). Therefore, it is important to compare the predictability of IDW and OK approaches for enhancing the prediction accuracy of spatial variability of soil properties. Previous studies have focused on the spatial variability of several soil nutrients measured at small scales, a single soil nutrient across large scales, or several soil nutrients at large scale but with sparse sampling. However, minimal research has been conducted at a regional scale to reliably quantify the spatial variability of SOC, N, P, and K with highdensity sampling in paddy fields, but it is important to do so as rice could sustain the growing population of China. The objectives of this study were to: (1) evaluate the SOC, N (total nitrogen and available nitrogen), P (total phosphorus and available phosphorus), and K (total potassium, rapidly available potassium and slowly available potassium) concentrations on a large scale (2.1 × 105 km2) using 8890 soil samples from paddy fields; and (2) investigate the spatial variation in soil nutrients and create a spatial distribution map.

2.3. Statistical analysis 2.3.1. Descriptive statistics Summary statistics (mean, median, minimum, maximum, interquartile range, skewness, and kurtosis) were determined to represent the general trends of the original datasets. The normal distributions of the datasets were tested using the 1-sample Kolmogorov-Smirnov test (K-S) (P < 0.05). Due to the large sample size, most of the variables were not normally distributed, and thus the interquartile range was also calculated to show the variation of the selected eight soil properties (Hu and Si, 2014). Descriptive statistical analysis was performed using SPSS 17.0 software.

2. Materials and methods 2.3.2. Geostatistical analysis To evaluate the spatial structure of the variability (i.e. heterogeneity), a semivariogram (half the expected squared difference between paired data separated by a distance, γ(h)) for each of the selected eight soil chemical properties was calculated as:

2.1. Study area The study was conducted across Hunan Province (108°47′–114°15′ E, 24°38′–30°08′ N; elevation 9–1953 m; area = 211800 km2) (Fig. 1a). Hunan Province is located in central China and features a continental, subtropical monsoon climate, with a mean annual air temperature of 16–18 °C, and an accumulated temperature (> 10 °C) of 5000–9500 °C (Chen et al., 2016). Annual precipitation ranges from 1250 to 1500 mm. The region is characterized by a basin with hills in central areas, surrounded by large mountains to the east, west, and south (Fig. 1b). The soil is derived from sandstone, shale, laterite, and purple sandstone, which are mainly acidic ferrisols and classified as Plinthudults according to US Soil Taxonomy, with a pH levels considered ‘acidic’ or ‘strongly acidic’. Hunan Province is a predominantly arable area and is one of the major rice-producing provinces in China. The planting pattern used within this region is the double-rice cropping system.

γ (h) =

1 2N (h)

N (h)

∑ i=1

[Z (x i + h) − Z (x i )]2

(1)

where N(h) refers to the number of pairs of points separated by h; Z (x i + h) and Z (x i ) are the sample values at the positions x i + h and x i , respectively. The empirical semivariograms were fitted by several theoretical models (e.g., exponential, Gaussian, linear, and spherical). The lowest residual sum of squares (RSS) and the greatest coefficient of determination (R2) were used to select the optimal model to derive the input parameters for spatial interpolation. Several important parameters (sill variance (C + C0), range (a), and nugget variance (C0)) were applied to explore the spatial structure of the soil chemical properties at a given scale. The theoretical variogram increased with increasing lag distance until it reached the maximum (sill, C + C0), at range (a). The ratio of nugget variance to sill variance was determined to identify the spatial dependency of the selected variables. Variables

2.2. Soil sampling and analysis Due to the large scale, heterogeneous landscape, a regular grid 2

Catena 188 (2020) 104350

L. Duan, et al.

Fig. 1. Locations of the (a) study area, (b) digital elevation model (DEM), and (c) sampling sites.

to a finite distance and estimates a variable at an unsampled location using the distance and values to nearby known data points (Phachomphon et al., 2010; Guan et al., 2017). The formula of the IDW method can be calculated as:

were regarded as ‘weak’, ‘moderate’, or ‘strongly spatial dependent’ when the calculated ratio was > 0.75, 0.25–0.75, or < 0.25, respectively (Cambardella, 1994). Only exponential (Eq. (2)) and linear (Eq. (3)) models were used and can be calculated as:

r (h) =

⎧0 h ⎨C0 + C 1 − e− a ⎩

⎧C0 ⎪ γ (h) = C0 + C ⎨ ⎪C0 + C ⎩

(

h=0

)

h>0

n

∑ z (x i ) dij−r (2)

Z ∗ (x 0) =

()

i=1

0 < h⩽a h>a

n

∑ dij−r

h=0 h a

i=1

(4)

where Z*(x0) is the predicted variables at sampling site x0; xi are the data points within a chosen neighborhood; dij is the distance between the predicted point and the data points; and r is the weight related to distance by dij. As a geostatistical method, kriging considers both the distance and the degree of variation between known data points. The interpolation procedure of the ordinary kriging started with selecting a suitable semivariogram model. Based on the calculation of semivariances and semivariogram models, ordinary kriging was applied to create prediction maps of the selected eight soil chemical properties. The prediction was calculated as the linear sum:

(3)

Geostatistical analysis was applied using GS+ software (Version 9.0). 2.3.3. Inverse distance weighting (IDW) and ordinary kriging method Inverse distance weighting (IDW) and ordinary kriging methods were used to perform spatial interpolations of the selected eight soil chemical properties. Detailed information of the underlying theories within these two approaches can be found in previous studies (Schöning et al., 2006). IDW assumes that each point affects the resulting surface 3

Catena 188 (2020) 104350

L. Duan, et al. n

∑ λi Z (xi)

Z ∗ (x 0) =

often a surplus of soil AN in croplands of China (Liu et al., 2014). Excess N exported from cropland has increased surface runoff of N which has contributed to contamination of drinking water and ecological degradation (Beman et al., 2005; Hou et al., 2015). The TP and AP both varied substantially from 0.10 to 2.5 g kg−1, and 0.7 to 315.5 mg kg−1, respectively (Table 3). The mean TP value (0.71 g kg−1) within the study region was similar to those in the cropland of the Loess Plateau (0.67 g kg−1), while the corresponding mean TN value (1.91 g kg−1) was approximately three times greater than those in the Loess Plateau (0.70 g kg−1) (Liu et al., 2013). The mean AP value (16.98 mg kg−1) within the study region was similar to the value (13.93 mg kg−1) reported by (Liu et al., 2014) in the paddy fields of southern China. The TK ranged from 0.5 to 70 g kg−1 with an average of 18.59 g kg−1 (interquartile range of 9.50 g kg−1). The mean TK within the study region was similar to that reported by Liu et al. (2012b) for the Loess Plateau (19.44 g kg−1). The RAK and SAK varied between 13 and 590 mg kg−1 (interquartile range of 62 g kg−1) and between 8 and 1671 mg kg−1 (interquartile range of 140 g kg−1), which fell within the range of reported values across China (He et al., 2015).

(5)

i=1

where λi is the coefficient or weight that can be obtained from the semivariogram. After IDW and ordinary kriging interpolation, finite discrete data points can be transformed to an infinite continuous surface with values for every possible location. The spatial distribution map of the eight soil chemical properties was produced by ArcGIS 10.2. 2.3.4. Model performance criteria Leave-one-out cross-validation was used to determine the interpolation quality of the eight soil chemical properties (Liu et al., 2012a; 2013; Guan et al., 2017). To determine the errors, one point was left out and its value is predicted based on the remaining values. Subsequently, the predicted value was compared with the observed value in the situation of the omitted point. This procedure was repeated n times until each value had been left out once. Performance of the models was evaluated using the average error (AE), mean absolute error (MAE), root mean square error (RMSE), and coefficient of determination (R2), Nash-Sutcliffe model efficiency (NSE). The following equations represent these statistical metrics:

AE =

3.2. Geostatistical analysis

n

1 n

∑ (Pi − Oi)

(6)

i=1

1 n

MAE =

∑ |Pi − Oi|

(7)

i=1

1 n

RMSE =

n

∑ (Pi − Oi)2

(8)

i=1

n

R2 =

Theoretically, the above classical analysis of the selected eight soil chemical properties assumed a randomized sampling design and independency of soil nutrients without spatial structure (Li et al., 2015b; Dai et al., 2018). Therefore, spatial dependency was considered using geostatistical analysis. Fig. 2 presents the geostatistical model of best-fit from each variogram, and Table 2 exhibited the best-fit of model parameters and spatial structural indices. A linear semivariogram model (R2 = 0.43) was the best fit for AP of the 8890 sampling sites, while the other seven soil chemical properties were fitted with an exponential model (R2 ranging from 0.53 to 0.84) (Table 2; Fig. 2). It is important to note that the R2 values are below 0.7 in most case (62.5%, Table 2). Due to the large sample size (8890), various transformations of the raw data did not improve the degree of normality. Hence, the raw data was used for the semivariogram analysis. Moreover, the variation of the selected soil chemical properties was large (Table 1). Therefore, the large sample size and large variability may exert substantial influence on the prediction results of the semivariogram analysis. For all the selected soil chemical properties, a strong spatial correlation was generally found, while SAK and AP demonstrated a moderate (0.37) and weak (0.96) spatial dependence (Table 2). The strong spatial dependence of the other six soil chemical properties may have been dominated by intrinsic characteristics such as soil type, soil forming processes, and geomorphologic properties (Wang et al., 2009; Guan et al., 2017; Liu et al., 2018). This contrasted with a weak spatial dependence of AP, which may have been controlled by extrinsic variations such as climatic conditions, soil fertilization, cultivation practices, and land use change (Kilic et al., 2004). The intrinsic and extrinsic variations were both responsible for the moderate spatial dependence of SAK. Range values (the distance beyond which the variable was spatially independent), varied from 17.1 km (TP) to 314.5 km (AP). The ranges for the selected eight soil chemical properties were generally greater than our sampling interval (Fig. 1c; Table 2), indicating that our sampling system was sufficient to detect the spatial heterogeneity of soil chemical properties within our study area.

n

2

[∑i = 1 (Pi − Pm)(Oi − Om)] n ∑i = 1

(Pi − Pm)2 (Oi − Om)2

n

(9)

n

∑ (Oi − Om)2 − ∑ (Pi − Oi)2 NSE =

i=1

i=1 n

∑ (Oi − Om)2 i=1

(10)

where Pi and Oi are the predicted and observed values, respectively, at location i, and Pm and Om are the mean of the predicted and observed values. respectively. 3. Results and discussion 3.1. Descriptive statistics Descriptive statistics of the eight selected soil chemical properties and their potential influencing factors are exhibited in Table 1. A wide range of values were observed for each of the selected soil chemical properties. The SOC of the 8890 soil samples varied considerably from 2.61 to 35.67 g kg−1 with a mean value of 19.72 g kg−1 (interquartile range of 7.95 g kg−1) (Table 1). The variation of the SOC was within the range (2.2–45.8 g kg−1) in the paddy fields of southern China (Liu et al., 2014). However, the mean SOC was higher when compared to data obtained from other regions such as the Loess Plateau, Yangtze River Delta Region, east China, and north China (Pan et al., 2003; Huang et al., 2007; Liu et al., 2012a) implying a high capacity to support rice productivity within the study region (Pan et al., 2009). The TN and AN of the paddy fields across the study area were 0.23–3.82 g kg−1 and 10–562 mg kg−1 (Table 1), which were within the scope of the TN values across China (0.15–4.18 g kg−1) (Yang et al., 2007), and similar to the scope within the paddy fields of southern China (0.29–4.20 g kg−1) (Liu et al., 2014). It is important to note that the AN values within this study were relatively high. Previous studies have also indicated that soil AN is not a limiting factor and there is

3.3. Comparison of interpolation performance Predicted and observed values were plotted to evaluate the interpolation performance (Figs. 3 and 4). The scatter point distribution patterns were different, indicating that the IDW and ordinary kriging methods may predict different values for the same sampling site. Generally, interpolation of both methods exhibited a centralization effect 4

Catena 188 (2020) 104350

L. Duan, et al.

Table 1 Summary statistics for eight soil chemical properties and their potential influencing factors (n = 8890). Variable SOC TN AN TP AP TK RAK SAK Longitude Latitude Elevation AT0 AT10 Precipitation Bulk density CEC pH

Unit −1

g kg g kg−1 mg kg−1 g kg−1 mg kg−1 g kg−1 mg kg−1 mg kg−1 ° ° m o C o C mm g cm−3 cmol kg−1 –

Mean

Median

Min.

Max.

Interquartile range

Skewness

Kurtosis

19.72 1.91 164.65 0.71 16.98 18.59 96.25 199.96 111.85 27.57 259.10 6258 5382 1408 1.17 27.07 6.06

19.55 1.89 162.00 0.71 12.30 17.70 83.00 160.00 112.1 27.56 190.00 6210 5353 1393 1.16 26.60 6.00

2.61 0.23 10.00 0.10 0.70 0.50 13.00 8.00 108.8 24.55 25.00 5609 4644 1200 0.80 2.10 4.10

35.67 3.82 562.50 2.50 315.50 70.00 590.00 1671.00 114.26 29.98 1308 6983 5841 1700 1.53 53.40 8.40

7.95 0.77 61.00 0.43 12.40 9.50 62.00 140.00 1.96 1.96 151.00 368.00 339.00 144.00 0.14 26.00 1.50

0.09 0.11 0.71 0.55 4.49 2.04 2.08 2.88 −0.44 −0.14 1.46 0.16 −0.14 0.49 0.63 0.07 0.29

−0.04 0.06 3.08 0.77 39.05 8.30 7.81 13.24 −0.87 −0.91 1.68 −0.18 0.11 −0.16 0.56 −1.23 −0.93

SOC, soil organic carbon; TN, total nitrogen; AN, available nitrogen; TP, total phosphorus; AP, available phosphorus; TK, total potassium; RAK, rapidly available potassium; SAK, slowly available potassium. Bold type denotes the selected eight soil chemical properties. Table 2 Modelled semivariogram parameters for the eight soil chemical properties. Variable

Models

Nugget (C0)

Sill(C0 + C)

Proportion [C0/(C + C0)]

Range (km)

R2

SOC TN AN TP AP TK RAK SAK

Exponential Exponential Exponential Exponential Linear Exponential Exponential Exponential

4.00 0.036 310 0.012 291.10 8.40 432 8380

37.30 0.355 2929 0.093 302.53 65.88 3189 22,740

0.107 0.101 0.106 0.126 0.962 0.128 0.135 0.369

23.10 21.00 25.2 17.1 314.47 27.00 26.10 70.50

0.527 0.527 0.838 0.578 0.426 0.735 0.621 0.844

SOC, soil organic carbon; TN, total nitrogen; AN, available nitrogen; TP, total phosphorus; AP, available phosphorus; TK, total potassium; RAK, rapidly available potassium; SAK, slowly available potassium. Table 3 Cross-validation indices for inverse distance weighted (IDW) and ordinary kriging methods. Longitude

Methods

SOC

IDW Ordinary IDW Ordinary IDW Ordinary IDW Ordinary IDW Ordinary IDW Ordinary IDW Ordinary IDW Ordinary

TN AN TP AP TK RAK SAK

kriging kriging kriging kriging kriging kriging kriging kriging

AE

MAE

RMSE

R2

NSE

P

−0.046 0.013 −0.005 0.001 −0.861 0.023 −0.001 0.001 −0.027 0.079 −0.026 −0.021 −0.028 0.045 1.564 0.235

4.369 4.310 0.436 0.427 35.191 34.340 0.255 0.237 10.487 10.197 5.662 5.387 36.297 35.292 73.568 72.278

5.692 5.472 0.565 0.540 49.600 47.022 0.329 0.297 18.499 16.982 8.201 7.617 54.673 51.264 119.282 115.333

0.215 0.219 0.191 0.198 0.231 0.256 0.031 0.062 0.041 0.044 0.101 0.131 0.153 0.175 0.415 0.431

0.148 0.213 0.115 0.191 0.167 0.251 −0.188 0.032 −0.173 0.012 −0.027 0.114 0.048 0.163 0.391 0.431

< 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01

SOC, soil organic carbon; TN, total nitrogen; AN, available nitrogen; TP, total phosphorus; AP, available phosphorus; TK, total potassium; RAK, rapidly available potassium; SAK, slowly available potassium. AE, average error; MAE, mean absolute error; RMSE, root mean square error; R2, coefficient of determination; NSE, NashSutcliffe model efficiency. IDW, inverse distance weighting method.

To compare the performance of the spatial interpolation, the AE, MAE, RMSE, R2 and NSE were all determined for the IDW and ordinary kriging methods (Table 3). In most cases, the AE was close to 0, implying that the IDW and ordinary kriging methods were all relatively unbiased in interpolating the selected soil chemical properties (Table 3). Interestingly, the lower MAE and RMSE values were observed in the ordinary kriging method, while lower R2 and NSE values were attained by the IDW method. These results indicated that the ordinary

where the large values were underestimated and low values were overestimated (Figs. 3 and 4). This may have been due to the nature of the algorithms that were used to achieve unbiased predictions of the mean values (Liu et al., 2013). However, the points were more widely spread when using the IDW method than that of the ordinary kriging method. This difference may be attributed to the commonly observed smoothing effect inherent in the ordinary kriging estimation method (Wang et al., 2009). 5

Catena 188 (2020) 104350

L. Duan, et al.

Semivariance

45

0.45

(a) SOC

30

0.30

15

0.15

0.00

0 0

Semivariance

4500

120000

240000

0

360000

0.12

(c) AN

3000

0.08

1500

0.04

0 0

Semivariance

450

120000

240000

120000

240000

360000

(d) TP

0.00

360000

0 75

(e) AP

300

50

150

25

120000

240000

360000

120000

240000

360000

240000

360000

(f) TK

0

0 0

120000

4500

Semivariance

(b) TN

240000

0

360000

30000

(g) RAK

3000

20000

1500

10000

0

(h) SAK

0

0

120000

240000

360000

0

Separation distance (m)

120000

Separation distance (m)

Fig. 2. Semivariograms (dots) and fitted model (line) of (a) soil organic carbon (SOC), (b) total nitrogen (TN), (c) available nitrogen (AN), (d) total phosphorus (TP), (e) available phosphorus (AP), (f) total potassium (TK), (g) rapidly available potassium (RAK), and (h) slowly available potassium (SAK).

6

Catena 188 (2020) 104350

L. Duan, et al.

Fig. 3. Cross-validation results of inverse distance weighting (IDW) interpolation of (a) soil organic carbon (SOC), (b) total nitrogen (TN), (c) available nitrogen (AN), (d) total phosphorus (TP), (e) available phosphorus (AP), (f) total potassium (TK), (g) rapidly available potassium (RAK), and (h) slowly available potassium (SAK).

7

Catena 188 (2020) 104350

L. Duan, et al.

Fig. 4. Cross-validation results of ordinary kriging (OK) interpolation of (a) soil organic carbon (SOC), (b) total nitrogen (TN), (c) available nitrogen (AN), (d) total phosphorus (TP), (e) available phosphorus (AP), (f) total potassium (TK), (g) rapidly available potassium (RAK), and (h) slowly available potassium (SAK).

8

Catena 188 (2020) 104350

L. Duan, et al.

Fig. 5. Spatial distribution maps of (a) soil organic carbon (SOC), (b) total nitrogen (TN), (c) available nitrogen (AN), (d) total phosphorus (TP), (e) available phosphorus (AP), (f) total potassium (TK), (g) rapidly available potassium (RAK), and (h) slowly available potassium (SAK) using the inverse distance weighting (IDW) method.

Fig. 6. Spatial distribution maps of (a) soil organic carbon (SOC), (b) total nitrogen (TN), (c) available nitrogen (AN), (d) total phosphorus (TP), (e) available phosphorus (AP), (f) total potassium (TK), (g) rapidly available potassium (RAK), and (h) slowly available potassium (SAK) using the ordinary kriging method.

kriging method performed better than the IDW method in interpolating soil chemical properties. This is consistent with previous studies (Liu et al., 2012a; Guan et al., 2017). Although the prediction values were significant for all of the eight soil chemical properties using the ordinary kriging method, the R2 (0.04–0.43) and NSE (0.01–0.43) values were relatively low (Table 3). Even though the sampling sites were large, by not using a probabilistic design, predictability of the spatial interpolation was generally reduced (Veronesi et al., 2014). Furthermore, some soil chemical properties such as AP had strong spatial variability which exerted substantial influence on the predicted results. The low R2 was also observed in previous studies where there was a large dataset or the variable had strong local variation (Liu et al., 2012a; Liu et al., 2013; Fu et al., 2014; Guan et al., 2017). These low R2 values indicate that a better methodology, such as state-space model, a

least square support vector machine, and partial least squares regression, should be developed to enhance the predictability of soil properties. Using the semivariogram parameters, prediction maps were produced for the eight soil chemical properties (Figs. 5 and 6). The SOC and TN varied widely across Hunan Province and showed a similar spatial pattern. This was characterized by large values in the central and south, and low values in the north. Available nitrogen also demonstrated substantial spatial variation across the study area and relatively high concentrations were observed in the southwest (Figs. 5 and 6). High levels of SOC, TN, and AN were found across almost the entire study area (Liu et al., 2014). High concentrations were observed during the non-planting season in some regions, which may have substantial consequences for environmental quality and human well-being 9

Catena 188 (2020) 104350

L. Duan, et al.

(Ersahin, 2001; Vitousek et al., 2009). Therefore, nitrogen fertilization should be reduced or correct fertilizer application management (e.g. time of nitrogen fertilizer application should coincide with rapid nitrogen uptake of crops) should be adopted (Beman et al., 2005; Leip et al., 2011). Unlike the SOC or TN, the TP and AP were relatively low across the entire study area and high values were only found sporadically (Figs. 5 and 6). Uneven phosphorus fertilizer application may be responsible for this spatial pattern. During the rice growing season, soil temperatures were high which can increase plant available phosphorus levels. Therefore, phosphorus fertilizer should not be applied in the regions with high TP or AP concentrations (Giardina et al., 2000; Lawrence and Schlesinger, 2001; Bruun et al., 2006). In contrast, to enhance or maintain paddy yields, phosphorus fertilizer should be continually applied in low TP or AP regions due to weak mobility of phosphorus in soils (Liu et al., 2014; Fischer et al., 2018). Total potassium concentrations were high in the south, while high RAK values were observed in the northwest and east. The SAK was generally low with scattered high levels (Figs. 5 and 6) implying widespread potassium depletion. An intensive multiple cropping system, coupled with a low potassium fertilizer application rate, led to a severe potassium deficiency (Liu et al., 2009; Zhang et al., 2010; Tan et al., 2012). Yang et al. (2004) also indicated that approximately 70% of paddy soil is considered potassium deficient in China. Some remedial measures such as rice straw incorporation can enhance potassium levels, thus increasing rice grain yield, particularly in regions of low potassium concentrations (Yuan et al., 2014). Knowledge of soil nutrient levels is required to make informed decisions regarding fertilizer application to support crop requirements.

rice cropping systems in subtropical China. Agric. Ecosyst. Environ. 129, 534–541. Bruun, T.B., Mertz, O., Elberling, B., 2006. Linking yields of upland rice in shifting cultivation to fallow length and soil properties. Agric. Ecosyst. Environ. 113, 139–149. Cambardella, C.A., 1994. Field-scale variability of soil properties in central Iowa soils. Soil Sci. Soc. Am. J. 58, 1501–1511. Chen, M., Chen, J., Sun, F., 2008. Agricultural phosphorus flow and its environmental impacts in China. Sci. Total Environ. 405, 140–152. Chen, L., Wang, S., Wang, Q., 2016. Ecosystem carbon stocks in a forest chronosequence in Hunan Province, South China. Plant Soil 409, 217–228. Dai, W., Zhao, K.L., Fu, W.J., Jiang, P.K., Li, Y.F., Zhang, C.S., Gielen, G., Gong, X., Li, Y.H., Wang, H.L., Wu, J.S., 2018. Spatial variation of organic carbon density in topsoils of a typical subtropical forest, southeastern China. Catena 167, 181–189. Ersahin, S., 2001. Assessment of spatial variability in nitrate leaching to reduce nitrogen fertilizers impact on water quality. Agr. Water Manage. 48, 179–189. Elbasiouny, H., Abowaly, M., Abu Alkheir, A., Gad, A., 2014. Spatial variation of soil carbon and nitrogen pools by using ordinary kriging method in an area of north Nile delta. Egypt. Catena 113, 70–78. Fischer, P., Pöthig, R., Gücker, B., Venohr, M., 2018. Phosphorus saturation and superficial fertilizer application as key parameters to assess the risk of diffuse phosphorus losses from agricultural soils in Brazil. Sci. Total Environ. 630, 1515–1527. Fu, W., Jiang, P., Zhao, K., Zhou, G., Li, Y., Wu, J., Du, H., 2014. The carbon storage in moso bamboo plantation and its spatial variation in Anji County of southeastern China. J. Soils Sed. 14, 320–329. Giardina, C.P., Sanford, R.L., Døckersmith, I.C., Jaramillo, V.J., 2000. The effects of slash burning on ecosystem nutrients during the land preparation phase of shifting cultivation. Plant Soil 220, 247–260. Guan, F., Xia, M., Tang, X., Fan, S., 2017. Spatial variability of soil nitrogen, phosphorus and potassium contents in Moso bamboo forests in Yong'an City, China. Catena 150, 161–172. Hati, K.M., Swarup, A., Mishra, B., Manna, M., Wanjari, R., Mandal, K., Misra, A., 2008. Impact of long-term application of fertilizer, manure and lime under intensive cropping on physical properties and organic carbon content of an Alfisol. Geoderma 148, 173–179. He, P., Yang, L., Xu, X., Zhao, S., Fang, C., Li, S., Tu, S., Jin, J., Johnston, A.M., 2015. Temporal and spatial variation of soil available potassium in China (1990–2012). Field Crops Res. 173, 49–56. Hou, X., Feng, Z., Leip, A., Fu, B., Hui, Y., Yan, C., Gao, S., Shang, Z., Ma, L., 2015. Spatial patterns of nitrogen runoff from Chinese paddy fields. Agric. Ecosyst. Environ. 231, 246–254. Hu, W., Si, B., 2014. Can soil water measurements at a certain depth be used to estimate mean soil water content of a soil profile at a point or at a hillslope scale? J. Hydrol. 516, 67–75. Huang, B., Sun, W., Zhao, Y., Zhu, J., Yang, R., Zou, Z., Ding, F., Su, J., 2007. Temporal and spatial variability of soil organic matter and total nitrogen in an agricultural ecosystem as affected by farming practices. Geoderma 139, 336–345. Kilic, K., Ozgoz, E., Akbas, F., 2004. Assessment of spatial variability in penetration resistance as related to some soil physical properties of two fluvents in Turkey. Soil Tillage Res. 76, 1–11. Lawrence, D., Schlesinger, W.H., 2001. Changes in soil phosphorus during 200 years of shifting cultivation in Indonesia. Ecology 82, 2769–2780. Leip, A., Britz, W., Weiss, F., Vries, W.D., 2011. Farm, land, and soil nitrogen budgets for agriculture in Europe calculated with CAPRI. Environ. Pollut. 159, 3243–3253. Li, D., Shao, M.A., 2014. Soil organic carbon and influencing factors in different landscapes in an arid region of northwestern China. Catena 116, 95–104. Li, T., Hasegawa, T., Yin, X., Zhu, Y., Boote, K., Adam, M., Bregaglio, S., Buis, S., Confalonieri, R., Fumoto, T., 2015a. Uncertainties in predicting rice yield by current crop models under a wide range of climatic conditions. Glob. Change Biol. 21, 1328–1341. Li, Z.W., Zhang, G.H., Geng, R., Wang, H., 2015b. Spatial heterogeneity of soil detachment capacity by overland flow at a hillslope with ephemeral gullies on the Loess Plateau. Geomorphology 248, 264–272. Liu, X., Zhang, W., Zhang, M., Ficklin, D.L., Wang, F., 2009. Spatio-temporal variations of soil nutrients influenced by an altered land tenure system in China. Geoderma 152, 23–34. Liu, Z.P., Shao, M.A., Wang, Y.Q., 2013. Spatial patterns of soil total nitrogen and soil total phosphorus across the entire Loess Plateau region of China. Geoderma 197–198, 67–78. Liu, Z., Zhou, W., Shen, J., He, P., Lei, Q., Liang, G., 2014. A simple assessment on spatial variability of rice yield and selected soil chemical properties of paddy fields in South China. Geoderma 235–236, 39–47. Liu, Z.P., Shao, M.A., Wang, Y.Q., 2012a. Large-scale spatial variability and distribution of soil organic carbon across the entire Loess Plateau, China. Soil Res. 50, 114–124. Liu, Z.P., Shao, M.A., Wang, Y.Q., 2012b. Spatial simulation of soil total potassium in regional scale for Loess Plateau region. Trans. CSAE 28, 132–140 (in Chinese). Liu, Z., Liu, Y., Baig, M., 2018. Biophysical effect of conversion from croplands to grasslands in water-limited temperate regions of China. Sci. Total Environ. 648, 315–324. Ministry of Agriculture, 2008. China agricultural yearbook. China Agriculture Press, Beijing (in Chinese). National Soil Survey Office, 1998. Soils of China. China Agricultural Press, Beijing, pp. 1016 (in Chinese). Pan, G., Li, L., Wu, L., Zhang, X., 2003. Storage and sequestration potential of topsoil organic carbon in China's paddy soils. Glob. Change Biol. 10, 79–92. Pan, G., Smith, P., Pan, W., 2009. The role of soil organic matter in maintaining the productivity and yield stability of cereals in China. Agric. Ecosyst. Environ. 129, 344–348.

4. Conclusions The current status and regional spatial variation of eight soil chemical properties (SOC, TN, AN, TP, AP, TK, RAK, and SAK) within paddy fields of Hunan Province were investigated. Compared with all of China, levels of SOC, TN, and AN were relatively high. AP exhibited strong variability, while the other seven soil chemical properties exhibited moderate variability. A weak and moderate spatial dependence was observed for AP and SAK, respectively. The remaining six soil chemical properties demonstrated strong spatial dependence. Generally, the ordinary kriging method performed better than the IDW method for spatial interpolation of eight soil chemical properties. However, a relatively low R2 and NSE were observed between the predicted and measured data. The large sampling strategy, coupled with a strong local variation of eight soil chemical properties with variability in environmental conditions, may have contributed to the low R2 and NSE values. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgments This study was supported by the Natural Science Foundation of Hunan Province (2019JJ50251), the National Natural Science Foundation of China (41807076), and the Farmland Quality Protection Special Project of Ministry of Agriculture (No: 10172130352501916) for the financial support. References Beman, J., Arrigo, Michael, K.R., Matson, P.A., 2005. Agricultural runoff fuels large phytoplankton blooms in vulnerable areas of the ocean. Nature 434, 211–214. Bi, L.D., Zhang, B., Liu, G.R., Li, Z.Z., Liu, Y.R., Ye, C., Yu, X.C., Tao, L., Zhang, J.G., Yin, J.M., 2009. Long-term effects of organic amendments on the rice yields for double

10

Catena 188 (2020) 104350

L. Duan, et al.

China. Geoderma 150, 141–149. Webster, R., Oliver, M.A., 2007. Geostatistics for Environmental Scientists. John Wiley & Sons. Xu, M., Lou, Y., Sun, X., Wang, W., Baniyamuddin, M., Zhao, K., 2011. Soil organic carbon active fractions as early indicators for total carbon change under straw incorporation. Biol. Fertil. Soils 47, 745. Yang, X., Liu, J., Wang, W., Ye, Z., Luo, A., 2004. Potassium internal use efficiency relative to growth vigor, potassium distribution, and carbohydrate allocation in rice genotypes. J. Plant Nutr. 27, 837–852. Yang, Y.H., Wen-Hong, M.A., Mohammat, A., Fang, J.Y., 2007. Storage, patterns and controls of soil nitrogen in China. Pedosphere 17, 776–785. Yuan, L., Zhang, Z., Cao, X., Zhu, S., Zhang, X., Wu, L., 2014. Responses of rice production, milled rice quality and soil properties to various nitrogen inputs and rice straw incorporation under continuous plastic film mulching cultivation. Field Crops Res. 155, 164–171. Zhang, H., Xu, M., Shi, X., Li, Z., Huang, Q., Wang, X., 2010. Rice yield, potassium uptake and apparent balance under long-term fertilization in rice-based cropping systems in southern China. Nutr. Cycl. Agroecosyst. 88, 341–349.

Phachomphon, K., Dlamini, P., Chaplot, V., 2010. Estimating carbon stocks at a regional level using soil information and easily accessible auxiliary variables. Geoderma 155, 372–380. Schöning, I., Kai, U.T., Kögel-Knabner, I., 2006. Small scale spatial variability of organic carbon stocks in litter and solum of a forested Luvisol. Geoderma 136, 631–642. Spiertz, J., 2009. Nitrogen, sustainable agriculture and food security: a review. Sustain. Agric. 635–651 Springer. Tan, D., Jin, J., Jiang, L., Huang, S., Liu, Z., 2012. Potassium assessment of grain producing soils in North China. Agric. Ecosyst. Environ. 148, 65–71. Usowicz, B., Lipiec, J., 2017. Spatial variability of soil properties and cereal yield in a cultivated field on sandy soil. Soil Tillage Res. 174, 241–250. Veronesi, F., Corstanje, R., Mayr, T., 2014. Landscape scale estimation of soil carbon stock using 3D modelling. Sci. Total Environ. 487, 578–586. Vitousek, P.M., Naylor, R., Crews, T., David, M.B., Drinkwater, L.E., Holland, E., Johnes, P.J., Katzenberger, J., Martinelli, L.A., Matson, P.A., 2009. Nutrient imbalances in agricultural development. Science 324, 1519–1520. Wang, Y., Zhang, X., Huang, C., 2009. Spatial variability of soil total nitrogen and soil total phosphorus under different land uses in a small watershed on the Loess Plateau,

11