Spatial prediction of soil organic matter in northern Kazakhstan based on topographic and vegetation information

Abstract

This study aimed to improve the accuracy of spatial prediction for soil organic matter, potential mineralizable carbon (PMC) and soil organic carbon (SOC), using secondary information, namely topographic and vegetation information, in northern Kazakhstan. Secondary information included elevation (ELEV), mean curvature (MEANC), compound topographic index (CTI) and slope (SLOPE) obtained from a digital elevation model, and enhanced vegetation index (VI) values obtained from a moderate resolution imaging spectroradiometer (MODIS). The prediction methods were statistical (multiple linear regression between soil organic matter and secondary information) and geostatistical algorithms (regression-kriging Model-C and simple kriging with varying local means [SKlm]). The VI, ELEV and MEANC were selected as the independent variables for predicting PMC and SOC. However, MEANC showed an opposite effect on PMC and SOC accumulation patterns. Model validity revealed that SKlm was the most appropriate method for predicting PMC and SOC spatial patterns because model validity revealed the smallest errors for this method. Maps from the kriged estimates showed that a combination of secondary information and geostatistical techniques can improve the accuracy of spatial prediction in study areas.

Key words:

• digital elevation model
• moderate resolution imaging spectroradiometer
• regression-kriging
• simple kriging
• spatial prediction

INTRODUCTION

Kazakhstan accounts for approximately 90% of the total cereal production area in Central Asia (International Maize and Wheat Improvement Center (CYMMIT) 2000), making it one of the most important regions for global food production. Northern Kazakhstan alone produces 70% of the national cereal production (Nobe 1998). This area is covered with fertile Chernozemic soils (Main Management of geodesy and cartography USSR (GUGK) 1982) and these soils are also considered to store a substantial amount of organic carbon. It is well documented that soil organic carbon (SOC) is an important indicator of soil quality and agronomic sustainability (Reeves 1997; Six et al. 2002) and functions as a huge source and sink of carbon dioxide (Lal 2004). The Chernozem soils in northern Kazakhstan are, therefore, a vital resource from both agricultural and environmental viewpoints.

Because of the political program Khrushchev's Virgin Lands Agricultural Program (1954–1960), the natural grass forb steppes of the area have disappeared over past decades and have largely been replaced by arable land, mainly for the cultivation of spring wheat (Medvedev 1987). Changes in agricultural land use can lead to changes in the amount of carbon held in soils, thus influencing the fluxes of CO₂ to and from the atmosphere (Falloon et al. 2002). Environmentally friendly land management benefits from analysis and forecast of the spatial distribution and dynamics of soil organic matter. Monitoring of the soil organic matter dynamics in northern Kazakhstan has shown that CO₂ emission under cereal cropping strongly depends on potential mineralizable carbon (PMC) content (Takata et al. 2007). To scale up this result to a statewide level, it is necessary to develop an accurate method of mapping PMC.

The spatial pattern of soil organic matter, which is controlled by environmental or state factors (Jenny 1941), is commonly studied using topographic and vegetation information. To evaluate the relationship between SOC and topographic information, quantitative topographic data have been applied to a digital elevation model (DEM), which is the basis for calculation of surface attributes, such as slope, aspect, curvature, and compound topographic index (CTI). These attributes have been used successfully to predict the spatial variability of SOC (Florinsky et al. 2002; Gessler et al. 2000; Moore et al. 1993; Muller and Pierce 2003; Terra et al. 2004). In addition to topographic information, vegetation information improves soil map accuracy for estimates of biomass variation within particular land uses (McBratney et al. 2003). This estimation method has been developed using visible and near-infrared vegetation indices, such as the Normalized Difference Vegetation Index (NDVI), which uses an Advanced Very High Resolution Radiometer (AVHRR) for natural vegetation and crops (Lobell et al. 2003; McKenzie and Ryan 1999). The Terra and Aqua satellites, launched in 1999 and 2002, respectively, are equipped with a Moderate Resolution Imaging Spectroradiometer (MODIS) that provides spatial resolution of 250 m for the red and near-infrared band and 500 m for the blue band. Huete et al. (1999) reported that the Enhanced Vegetation Index (EVI) calculated using these bands has a higher performance for predicting plant biomass than the AVHRR. Therefore, EVI appears to be a powerful tool for global environmental monitoring. The availability of secondary information, such as topographic and vegetation information, would support our understanding of the spatial dependency of soil organic matter (the primary target variable), particularly when the soil organic matter data are sparse or poorly correlated in space (Goovaerts 1997).

Most secondary information-based predictions of soil organic matter are built on statistical models, such as multiple linear regression, that describe the relationships between soil and secondary information at each point (Florinsky et al. 2002; Gessler et al. 2000; Moore et al. 1993; Muller and Pierce 2003; Terra et al. 2004). However, a major disadvantage of the multiple regression method is that the primary target variable at a particular grid node is derived only from secondary information at each point, regardless of the records at the surrounding measured data. Such an approach assumes that the predicted residual values are spatially uncorrelated. Spatial correlation of the residuals can be taken into account, however, using two types of geostatistical algorithms, namely simple kriging (Goovaerts 1997) and ordinary kriging (Odeh et al. 1994). Geostatistics is a branch of applied statistics that quantifies the spatial dependence of a measured property, and uses the spatial dependence of the measured variables to predict values of the variable at unsampled locations (Mulla and McBratney 2000).

While many studies have focused on the relationship between SOC and secondary information (Florinsky et al. 2002; Gessler et al. 2000; Moore et al. 1993; Muller and Pierce 2003; Terra et al. 2004), relatively few studies to date have characterized the spatial variability of PMC with secondary information. Moreover, spatial prediction of PMC in a low relief agricultural landscape of the Eurasian steppe is limited. Therefore, our first objective is to clarify the relationship between the distribution pattern of SOC and PMC and secondary information. Then, to advance the study of soil organic matter dynamics on a regional scale at this location, our second objective is to test the performance of the spatial predictions for cereal cropping fields on a regional scale using three methods.

MATERIALS AND METHODS

Study area

The main cereal production areas of the Republic of Kazakhstan are located in the Aqmola and Northern Kazakhstan regions. These areas are characterized by a semiarid and continental climate with high temperature fluctuations (−40°C in winter and +35°C in summer) and a short growing season (May–September). Average annual precipitation ranges from 280 to 350 mm. The study area was chosen to incorporate a central area covering both regions; 54°04′N 69°34′E and 50°43′N 70°49′E in the upper left and lower right areas of Fig. 1, respectively. According to local soil classification (Fig. 2), the chosen area is characterized by three different soils: Ordinary Chernozem (OC), Southern Chernozem (SC) and Dark Chestnut (DC), which in US Taxonomy correlate to Pachic Haplustolls, Typic Haplustolls and Entic Haplustolls, respectively. In the study area, OC is distributed in the lowest altitude zone of the northwest and the highest altitude zone of the west. Southern Chernozem is distributed in the low altitude zone of the north to northeast and the medium altitude zone of medium latitude of this area, and DC is distributed only in the southern region of the study area.

Agricultural land extraction

The boundaries of the agricultural lands were delineated by viewing interpretation using Landsat ETM⁺ and Aster images taken on 17 August 2002 and 19 June 2004, respectively. The total area of extracted arable fields was 2.5 million hectares.

Figure 1  Location of study area.

Figure 2  Study area description a. Shows the distribution of cereal cropping fields and soil sampling sites. b. Shows the elevation (m) of study area. c. Shows the soil map.

Secondary information

Digital elevation model data at a resolution of 90 m of the study area were downloaded from the NASA website (STRM-3: http://www2.jpl.nasa.gov/srtm/). The downloaded DEM data were resampled to 250 m pixel size. The terrain attributes included elevation (ELEV), slope (SLOPE), mean curvature (MEANC) and compound topographic index (CTI). The MEANC is the average of normal section curvature. Negative MEANC values denote mean-concave landforms, while positive values denote mean-convex landforms. The CTI represents landscape areas where water accumulates (Moore et al. 1993) and is computed using algorithms and commands by Geographical Information System ILWIS 3.3 software (ITC, Enschade, the Netherlands).

In this study, we used MODIS 250 m, 16-day composite vegetation index products (MOD13Q1). The MOD13Q1 product contains two vegetation indices: the NDVI and the EVI. The NDVI is the most popular index and the most used in various remote sensing applications. This index, however, is subject to the influences of soil and air conditions. For combined correction of optical properties of soil background and the atmospheric effects, Huete et al. (1999) proposed EVI. Thus, EVI was used in this study. The MOD13Q1 data for spring (23 April to 23 May) and summer (12 July to 12 August) seasons for 2000 through 2005 were downloaded from the MODIS website (http://edcimshttp://www.cr.usgs.gov/pub/imswelcome/). We computed the 5-year mean EVI for the summer season as secondary information (VI).

To eliminate meadow fields and rangelands from the extracted agricultural lands, the phenological difference between cereal and perennial grass was used. Germination in the meadow fields and rangelands of the study region starts from early April and active growth occurs in May, while cereal seeding starts in late May with active growth occurring in July. Therefore, threshold values were set as follows:

where EVI is a mean value of the spring term (23 April−23 May). This extraction was conducted by each year. Using these thresholds, all ground-checked pasture fields (approximately 5000 ha) were extracted.

Soil sampling and soil organic matter measurement

In this study, the spatial prediction of soil organic matter is limited to cereal fields. Soil sampling (518 sites: Fig. 2) was conducted between May 2003 and August 2005, and the sampling was carried out to cover the arrangement of cereal fields. Samples of surface soil (0–15 cm) collected from six points at each site using a hand auger were mixed and used as composite soil samples. In addition, core soil samples (100 cc) were collected from the center of each site to determine soil bulk density. Latitude and longitude data of each sampling site were obtained using Global Positioning System (GPS) (Garmin eTrexVenture).

McBratney et al. (1991) suggested that when several different soils are distributed in a study area, the spatial variation of locally defined regions within the area needs to be examined separately. As three soil classes existed in the study area, the geostatistical mapping of surface soil organic matter content would likely be improved by analyzing the three local regions independently. Thus, based on the soil map of Kazakhstan, the 518 sampling sites were classified into the three soil types: 131 sites with OC, 204 sites with SC and 183 sites with DC.

As the primary target variables, total soil organic carbon (SOC) and potential minerarizable carbon (PMC) were determined. The SOC content of the 518 sites was measured using the Tyurin method, and PMC content of 434 sites (116 OC sites, 193 SC sites and 106 DC sites) was determined by incubating 15 g soil samples for 19 weeks at an air temperature of 30°C and 60% water field capacity. CO₂ was collected using an alkaline trap in a closed chamber and measured with 0.1 mol L⁻¹ HCl. Non-linear regression was used to describe PMC according to the following equation (SPSS Inc. 1998):

where C_min is the quantity of mineralized C (mg kg⁻¹ dry soil) at time t (day), C₀ is PMC (mg kg⁻¹), and k is a non-linear mineralization constant (i.e. fraction mineralized per day).

Summary statistics of the samples are given in Table 1. Mean PMC and SOC values of surface soil (0–15 cm) are highest at the OC sites and lowest at the DC sites. The coefficient of variation (CV) for PMC was higher than that for SOC.

Prediction methods

Multiple linear regression

Multiple linear regression was used to predict the dependent variable (u_α) at a certain location u (the primary targets of PMC and SOC in this study) from secondary information value X, which contains sampling plots in the pixels (pixel size is 250 m):

where y(u_α) is the dependent variable, X₁, X₂, . . , X_k are the independent variables, and e is the random error. b₀, b₁, b₂, . . . b_k are the regression coefficients, which have to be estimated from the data.

Regression-kriging Model C

The regression-kriging model C (Ahmed and Demarsily 1987; Odeh et al. 1995) combines multiple regression and ordinary kriging. In the first step, multiple regression of the primary target variables z(u_α) with secondary information was used to create the predicted value

Table 1 Descriptive statistics of potential mineralizable carbon (PMC) and soil organic carbon (SOC)

	PMC					SOC
	Mean (Mg C ha^−1)	Minimum (Mg C ha^−1)	Maximum (Mg C ha^−1)	SD (Mg C ha^−1)	CV (%)	Mean (Mg C ha^−1)	Minimum (Mg C ha^−1)	Maximum (Mg C ha^−1)	SD (Mg C ha^−1)	CV (%)
DC	2.29	1.03	4.84	0.58	25.3	30.8	17.8	41.3	4.05	13.1
SC	2.85	1.35	5.75	0.90	31.5	39.1	17.6	53.6	7.75	19.8
OC	3.17	1.22	5.85	1.05	32.9	49.6	29.3	65.5	6.34	12.8
CV, coefficient of variation; DC, Dark Chesnut; OC, Ordinary Chernozem; SC, Southern Chernozem; SD, standard deviation.

y(u_α). The residuals of the regression (Eq. 1), r(u_α), can be obtained for the sample points as follows:

In the second step, ordinary kriging was applied to the residuals r(u_α). Ordinary kriging is a univariate interpolation method based on a weighting scheme to estimate the unknown primary target variable at unsampled locations by neighboring observation. Kriging minimizes the estimator variance and ensures the unbiasedness of the estimator (Webster and Oliver 2001). As a first step in applying ordinary kriging to residuals r(u_α), the spatial dependency of each residual r(u_α) was characterized by the semi-variogram γ(h) using the following equation:

(3)

where γ(h) is the experimental semi-variogram value at distance interval h, N(h) is the number of sample value pairs within the distance interval h, and r(u_α) and r(u_α+h ) are residual values at two points separated by the distance interval h. Spherical or exponential models were fitted to the experimental semi-variogram defined by Eqs 4 and 5, respectively, using geostatistical software Gs⁺ (Gamma design software):

(4)

where a is the range, c₀ is the nugget, and c₀ + c is the sill. Nugget is the variance at zero distance; range is the lag distance between measurements at which one value of one variable does not influence neighboring values, and sill is the asymptotic plateau of the semivariogram function and is used to estimate the range. Often the sill is approximately equal to the sample variance (Journel and Huijbregts 1978). If the semivariogram function is γ(h) = c, then the slope of the semivariogram was close to zero; the soil variable is considered non-spatially correlated (pure nugget effect) and does not distribute in patches. Once the model variogram was constructed, ordinary kriging was applied to residuals r(u_α) to compute the kriging weights λ for interpolating at unsampled locations (Eq. 6) using the neighboring residuals r(u_α):

(6)

λ_α ^OK(u) is the ordinary kriging weight that minimizes the estimation variance, Var{r*_OK(u) − r(u)}, while ensuring the unbiasedness of the estimator, E{r*_OK(u) − r(u)} = 0. These weights are obtained by solving a system of linear equations that is known as the ordinary kriging system (Eq. 7):

(7)

where µ(u) is the Lagrange parameter accounting for the constraint on the weights, and n(u) is the numbers in the actual dataset. Semivariogram values are the only information required for solving the kriging system. Ordinary kriging were applied to residuals r(u_α) to interpolate at unsampled locations using the neighboring residuals r(u_α) with the weights λ_α ^OK(u).

In the third step, the spatially distributed regression result and the kriged residuals are summed to estimate the target variable Z*_Model-C:

Simple kriging with varying local means

As with the regression-kriging Model-C approach, simple kriging with varying local means (SKlm) combines spatially distributed regression results and the residuals that are interpolated by simple kriging (Goovearts 1997):

(9)

(10)

where the weights λ_α ^OK(u) are obtained by solving the simple kriging system:

(11)

where C_R(u_α − u_β) is the covariance of the actual data-to-actual data, and C_R(u_α − u) is the covariance of the actual data-to-estimated data. Unlike ordinary kriging (Eq. 7), simple kriging can be expressed in terms of only covariance because there is a lack of constraints on the kriging weights (Goovaerts 1997) and the weights do not sum to unity.

Model validity

To predict soil organic matter, we utilized the following soil class and whole datasets: OC dataset (86 sites for PMC; 111 sites for SOC); SC dataset (173 sites for PMC; 184 sites for SOC); DC dataset (96 sites for PMC; 163 sites for SOC); and the whole dataset (355 sites for PMC; 458 sites for SOC). Another group of sample points at different locations (n = 20 in each soil class) were used for testing each of the prediction methods; where we computed two indices from the observed and predicted values as used by Voltz and Webster (1990) and Odeh et al. (1994). As there are n sites belonging to the validation sample, the indices of validity computed from the observed values z(s_j) and the predicted values z*(s_j) are the mean error (ME) and the root mean square error (RMSE). The mean error is defined as:

(12)

This measures the bias of prediction and should be close to zero for unbiased methods.

The RMSE is expressed as:

(13)

This is a measure of precision of prediction and should be as small as possible for precise prediction.

RESULTS AND DISCUSSION

Relationship between soil organic matter and secondary information

The best combination of secondary information for explaining variation in primary target values was explored using stepwise multiple linear regression. Table 2 shows the equation of prediction for soil class and whole datasets. For all PMC predicting equations and the SOC soil class dataset, R² was very low, which is in agreement with the results of Florinsky et al. (2002), Hengl et al. (2004) and Lopez-Granados et al. (2005). However, all predicting equations showed a significant relationship at the 1% level, indicating that there was still a weak relationship between soil organic matter and secondary information. In the case of the SOC whole dataset, prediction performance was high, at the same level reported by Moore et al. (1993), McKenzie and Ryan (1999) and Dai and Huang (2006). Both the PMC and SOC whole datasets showed zones with higher VI values and lower ELEV values for sites with high soil organic matter accumulation (Table 2). As ELEV showed positive and negative effects for soil class datasets, ELEV for whole datasets probably represents regional variation in PMC and SOC.

The MEANC showed different influences on PMC and SOC accumulation patterns (Table 2). PMC showed high accumulation in convex zones in whole, OC and SC datasets. In contrast, SOC showed high accumulation in concave zones, as reported by other studies (Florinsky et al. 2002; Gessler et al. 2000; Moore et al. 1993). Table 3 shows the correlation matrix of PMC and secondary information, and Table 4 shows SOC and secondary information. The PMC showed a significant positive correlation with SLOPE in the whole and OC datasets, whereas SOC showed a negative correlation with SLOPE in whole and individual soil class datasets. In general, SOC is accumulated in concave zones through sediment transport by water erosion (Gregorich et al. 1998; Pennock 2003; Schimel et al. 1985;

Table 2 Multiple linear regression results of potential mineralizable carbon (PMC) and soil organic carbon (SOC)

	Dataset	N	Equation	R^2
PMC	DC	96	Y = 3.72 + 3.40VI − 0.007ELEV	0.13**
	SC	173	Y = 3.81 − 0.003ELEV**+ 93.4MEANC*	0.07**
	OC	86	Y = 2.11 + 0.004ELEV + 55.93MEANC	0.18**
	Whole	355	Y = 1.72 + 4.27VI** − 0.001ELEV** + 69.91MEANC*	0.08**
SOC	DC	163	Y = 10.78 + 43.1 VI − 420.94MEANC + 0.02ELEV	0.24**
	SC	184	Y = 37.41 + 33.0VI − 619.35MEANC − 0.02ELEV − 3.57SLOPE	0.14**
	OC	111	Y = 30.50 + 1.273CTl**	0.05**
	Whole	458	Y = 28.45 + 79.27VI** − 0.053ELEV** − 438.04MEANC*	0.44**
*P < 0.05;
**P < 0.01.
CTI, compound topographic index; DC, Dark Chesnut; ELEV, elevation; OC, Ordinary Chernozem; MEANC, mean curvature; SC, Southern Chernozem; VI, vegetation index.

Table 3 Correlation matrix of potential mineralizable carbon (PMC) and secondary information

		PMC	VI	MEANC	ELEV	CTI
DC: N = 96	SOC	0.48**
	VI	0.26**
	MEANC	0.08	−0.08
	ELEV	−0.26**	0.11	−0.07
	CTI	−0.11	−0.05	−0.64**	−0.03
	SLOPE	0.02	−0.27**	0.16	0.01	−0.40**
SC: N = 173	SOC	0.47**
	VI	−0.06
	MEANC	0.12	0.06
	ELEV	−0.25**	0.19*	0.06
	CTI	−0.13	−0.02	−0.47**	0.06
	SLOPE	−0.06	0.00	0.01	0.03	−0.33**
OC: N = 86	SOC	0.28*
	VI	0.24*
	MEANC	0.12	0.10
	ELEV	0.38**	0.46**	0.09
	CTI	−0.01	−0.20	−0.29**	−0.38**
	SLOPE	0.26*	0.25*	−0.08	0.56″	−0.45**
Whole: N = 355	SOC	0.53**
	VI	0.26*
	MEANC	0.11*	0.06
	ELEV	−0.16**	0.12*	0.05
	CTI	−0.09	−0.07	−0.47**	−0.05
	SLOPE	0.11*	0.10	0.01	0.17*	−0.37**
*P < 0.05;
**P < 0.01.
CTI, compound topographic index; DC, Dark Chesnut; ELEV, elevation; OC, Ordinary Chernozem; MEANC, mean curvature; SC, Southern Chernozem; SOC, soil organic carbon; VI, vegetation index.

Table 4 Correlation matrix of soil organic carbon (SOC) and secondary information

		SOC	VI	MEANC	ELEV	CTI
DC: N = 163	VI	0.46**
	MEANC	−0.18*	−0.07
	ELEV	0.20**	0.20**	−0.03
	CTI	0.03	0.00	−0.64**	−0.02
	SLOPE	−0.17*	−0.26**	0.19*	−0.08	−0.41**
SC: N = 184	VI	0.11
	MEANC	−0.17*	0.09
	ELEV	−0.25**	0.19**	0.07
	CTI	0.07	−0.03	−0.52**	0.03
	SLOPE	−0.25**	0.03	0.12	0.117	−0.37**
OC: N = 111	VI	0.05
	MEANC	−0.04	0.07
	ELEV	−0.15	0.44**	0.07
	CTI	0.22*	−0.18	−0.34**	−0.35**
	SLOPE	−0.21*	0.28**	−0.08	0.61**	−0.42**
Whole: N = 458	VI	0.40**
	MEANC	−0.07	0.07
	ELEV	−0.50**	0.04	0.05
	CTI	0.06	−0.08	−0.52**	−0.09
	SLOPE	−0.20**	0.13*	0.09	0.21**	−0.39**
*P < 0.05;
**P < 0.01.
CTI, compound topographic index; DC, Dark Chesnut; ELEV, elevation; OC, Ordinary Chernozem; MEANC, mean curvature; PMC, potential mineralizable carbon; SC, Southern Chernozem; VI, vegetation index.

Table 5 Geostatistical parameters for the primary target variables and the residuals

	Dataset	N	Model	Nugget	Sill	Q value	Range (km)	R^2
Primary target variables
PMC	DC	96	Exponential	0.16	0.65	0.76	198	0.72
	SC	173	Exponential	0.45	0.90	0.50	194	0.57
	OC	86	Spherical	0.57	1.07	0.47	57.2	0.41
	Whole	355	Exponential	0.46	2.60	0.82	411	0.66
SOC	DC	163	Spherical	8.34	16.9	0.51	399	0.42
	SC	184	Spherical	44.6	89.2	0.50	158	0.26
	OC	111	Nugget effect	−	−		–	−
	Whole	458	Spherical	23.0	127	0.82	203	0.85
Residuals
PMC	DC	96	Exponential	0.22	0.46	0.53	211	0.53
	SC	173	Exponential	0.59	1.18	0.50	48.4	0.73
	OC	86	Spherical	0.37	1.00	0.63	211	0.68
	Whole	355	Exponential	0.49	2.03	0.76	411	0.54
SOC	DC	163	Nugget effect	−	−		−	−
	SC	184	Spherical	43.3	86.6	0.50	311	0.24
	OC	111	Nugget effect	−	−		−	−
	Whole	458	Spherical	33.5	107	0.69	302	0.73
DC, Dark Chesnut; OC, Ordinary Chernozem; PMC, potential mineralizable carbon; SC, Southern Chernozem; SOC, soil organic carbon.

Yoo et al. 2006) and SOC, therefore, tends to show negative correlation with slope. Sediment transport by water erosion was reported to cause PMC transport to concave zones by Gregorich et al. (1998) and Jacinthe et al. (2004) and, thus, their data conflicts with our data. However, Schemel et al. (1985) and Woods and Schuman (1988) showed that wind erosion removes surface soil particles from exposed areas and deposits them in protected areas, regardless of slope position. In addition, these researchers showed that the correlation between PMC and SOC at wind erosive sites was lower than at water erosive sites. Our result also shows a lower correlation coefficient (r = 0.28) between PMC and SOC in the OC zone (Table 3), where soil might have high wind erodibility associated with a high organic matter content (Table 1) and light particles, when compared to other datasets (r = 0.48 in the DC dataset and r = 0.47 in the SC dataset).

Geostatistical analysis

Table 5 summarizes the geostatistical parameters for the raw data and residuals. The Q value (calculated as [sill-nugget]/sill) shows the degree of development of spatial dependence at the sampling scale. The Q values for the raw data of the PMC and SOC whole datasets were higher than those of the soil class datasets, as were the Q values for PMC and SOC residuals. These results indicate that the spatial dependence of the whole dataset was more developed than that of each soil class dataset. In addition, similar to Herbst et al. (2006), the Q values for residuals in the present study were found to be lower than those of the primary target variables.

Pure nugget effects were identified for the primary target SOC in the OC dataset and for the SOC residuals in the DC and OC datasets. In the case of the DC and OC datasets, the SOC residuals were spatially uncorrelated, hence, all kriging weights in Eq. 6 and 7 are zero and the Model-C and SKlm estimate is the same value as that provided by linear regression (Goovaerts 2000). The ranges of residuals tend to be longer than those of raw data, indicating that the trend (Journel and Huijbregts 1978) could be removed (Hengl et al. 2004). The same result was also reported by Baxter and Oliver (2005). As listed in Table 6, the range parameters of VI and ELEV are 118 and 390 km, respectively. In addition, the MEANC range should be shorter than 8 km because there was a nugget effect on MEANC (the least lag was set at 8 km in this study). The VI, ELEV and MEANC no longer affect the variogram of residuals in Table 5. Accordingly, the ranges of residuals should show the influence of large spatial dependence factors such as climate. In addition, a long range indicates that the soil variable is influenced by other values of this variable over greater distances than soil variables that have shorter range (Lopez-Granados et al. 2002). Accordingly, the PMC and SOC whole datasets had a longer range than the soil class datasets, indicating that the PMC and SOC whole datasets influenced neighboring values over greater distances than the soil class datasets.

Table 6 Geostatistical parameters for the topographical and vegetation indices

	N	Model	Nugget	Sill	Q value	Range (km)	R^2
VI	458	Spherical	0.001	0.004	0.750	117.9	0.871
MEANC	458	Nugget effect	−	−	−	−	−
ELEV	458	Exponential	10.00	9130	0.999	390.1	0.880
CTI	458	Nugget effect	−	−	−	−	−
CTI, compound topographic index; ELEV, elevation; MEANC, mean curvature; VI, vegetation index.

Table 7 Performance of spatial prediction methods of potential mineralizable carbon (PMC) and soil organic carbon (SOC)

		Validation points		Multiple Regression		Model−C		SKlm
		Mean (Mg C ha^−1)	SD	ME	RMSE	ME	RMSE	ME	RMSE
PMC	Soil Class*	2.82	0.85	−1.64	0.816	−0.51	0.789	−0.63	0.739
	Whole	2.82	0.85	2.30	0.741	−0.93	0.743	−1.19	0.723
SOC	Soil Class	39.4	9.15	5.31	9.853	5.66	10.34	0.81	5.718
	Whole	39.4	9.15	−0.13	7.197	0.60	5.801	0.19	5.571
*Three of the soil class datasets (Dark Chesnut, Ordinary Chernozem and Southern Chernozem) were glued to a map, and the mean error (ME) and root mean square error (RMSE) were calculated using a total of 60 validation points. SD, standard deviation.

Performance of different methods

The performance of spatial prediction methods of PMC and SOC are listed in Table 7. For PMC, ME and RMSE ranged from −1.64 to 2.30 and from 0.72 to 0.82, respectively. The most unbiased (ME-based) model was Model-C for the soil class dataset, and the most precise (RMSE-based) model was SKlm for the whole dataset. The RMSE of almost all models was lower in the whole dataset than in the soil class dataset. For SOC, ME and RMSE ranged from −0.13 to 5.67 and from 5.57 to 9.85, respectively. The most unbiased model was multiple regression for the whole dataset, and SKlm for the whole dataset showed the lowest RMSE. The RMSE of almost all models was lower in the whole dataset than in the soil class dataset. These findings suggest that a larger number of samples and larger value range improve the accuracy of the multiple linear regression model because the whole dataset had both a larger number of samples and larger value range.

When comparing the two types of regression-kriging, SKlm is more consistent than Model-C. The difference between SKlm and Model-C is the calculation constraint of weight for estimating regression residuals. In general, ordinary kriging is more robust in the presence of a trend than simple kriging (Deutsch and Journel 1992; Journel and Rossi 1989). While ordinary kriging is often used for kriging of the residuals (Hengl et al. 2004), simple kriging is more consistent when the variable has a large range without trend, because it uses all datasets to interpolate non-sampled sites without autocorrelation, which yields a zero estimated residual value (Simbahan et al. 2006). In the present study, the PMC and SOC residuals had large ranges with relative high Q values (Table 5). A longer range indicates that the soil variables, which are located over a great distance, influence the estimation of the target variable. Thus, simple kriging, which uses all datasets, was more consistent than ordinary kriging, which uses only neighboring datasets.

Figure 3 shows the predicted PMC and SOC maps using the SKlm model on whole datasets, which was found to be an adequate method with the lowest RMSE value in this study (Table 7). Also shown is the effect of different prediction methods on the PMC map. Both regression-kriging methods showed that the spatial distribution of PMC could be influenced by cereal field boundaries (Fig. 3). This finding indicates the possibility that PMC content could depend on field management, including frequency of summer fallow.

The highest SOC content was found in the northern region and the lowest content in the southern region of the study area. The PMC content, however, was highest in the mid-western region. It is noteworthy that the highest SOC content zone did not coincide with the highest PMC content zone. In other words, the zone that has the largest sink of organic matter did not coincide with the largest source of CO₂ on a state level. Conservation of soil organic matter in the highest PMC zones can require more careful land management. Takata et al. (2007) reported that in situ CO₂ emission from soil correlated with PMC content. Moreover, they identified that a summer fallow-based crop rotation system leads to soil organic matter depletion in the lowest PMC content zones (southern part of the study area).

Figure 3  Predicted potential mineralizable carbon map. a. Soil organic carbon map. b. Using simple kriging with varying local means, and the effect of different prediction methods on potential mineralizable carbon map created by multiple regression (PMC-1), Model-C (PMC-2), and SKIm (PMC-3).

For the introduction of site-specific agricultural management that reduces organic matter depletion, further studies to predict the spatial variability of actual CO₂ emission are warranted.

Conclusions

The independent variables for predicting PMC and SOC were selected as VI, ELEV and MEANC; however, MEANC showed an opposite effect on PMC and SOC accumulation patterns, suggesting that water and wind erodibility and biomass productivity might be correlated with topography. Consequently, the spatial dependence of soil organic matter was influenced by the spatial distribution of topography and vegetation. The residuals of multiple linear regressions revealed a potential for the spatial prediction of PMC and SOC. The most appropriate method for predicting the spatial patterns of PMC and SOC was SKlm using a whole dataset, as model validity revealed the smallest errors for this method. It is well known that cultivation history influences the soil organic matter content. However, most of the cereal field was reclaimed during 1954 to 1960, and nearly 50 years of monoculture of wheat, summer fallows had been practiced in crop rotations in the region. Therefore, cultivation history would be one of the most important factors for predicting the spatial pattern of soil organic matter on the other region. PMC and SOC did not exhibit similar distribution patterns, suggesting that the difference in spatial patterns of carbon sink and source should be considered in environmentally friendly agricultural land management.