Abstract: This study aimed at quantifying the spatial variability of Soil Organic Carbon (SOC), estimating SOC at unsampled locations and comparing the spatial variability of SOC between young and mature oil palm stands. Two study sites were chosen to represent two different palm age groups, i.e., 5 Years after Planting (YAP) and 17 YAP. A systematic sampling design was employed for soil sampling at the 0-20 cm depth based on a cluster of four palms that comprised three operational zones: Weeded Circle (WC), Frond Heap (FH) and Harvesting Path (HP). A total of 60 sampling clusters were obtained for each site. Soil samples were analyzed for SOC by dry combustion method. All measurement points were geo-referenced by differential Global Positioning System (dGPS). The SOC data were first explored using descriptive statistics, normality check, outlier detection and data transformation, followed by variography and interpolation. Spatial variability of SOC was mapped based on measured and kriged values. Results showed that all operational zones exhibited a definable spatial structure, which were described by either spherical or exponential models. All operational zones exhibited strong spatial dependence. Operational zones of 5-year old palms exhibited a shorter effective range than those of 17 year old palms. Additionally, SOC heterogeneity was evident among operational zones at both sites, where FH registered the highest SOC, followed by WC and HP. SOC concentration at 17 year old palms was found to be more stable than that from 5 year old palms. This study suggests spatial variability assessment appears to be a feasible technique to quantify the variability of SOC in oil palm cultivation.
INTRODUCTION
Studies on soil Carbon (C) have generated much interest because understanding the dynamics of soil C is an important part of sustainable agriculture strategies to increase soil fertility and accelerate C sequestration (Liu et al., 2006). In terrestrial ecosystem, soil plays an essential role in global C balance because of its ability to offset greenhouse gas emission through C sequestration. The C stored in soil is estimated to be four times greater than the total available in living vegetation. (Lal, 2004). Numerous studies have been done on measuring the quantity and the spatial distribution of Soil Organic Carbon (SOC) in order to refine agricultural management strategies that promote sustainable land use (Law et al., 2009; Jian-Bing et al., 2006; Liu et al., 2006; Franzluebbers, 2005; Follett, 2001; Chevallier et al., 2000). Nevertheless, SOC studies specifically aimed at oil palm are few and inadequately understood.
Large-scale planting of oil palm and its large biomass production has a high potential for soil C sequestration within its economic life span, which typically ranges from 25 to 30 years (Henson, 1999). In Malaysia, the total oil palm average in 2008 was 4.48 million hectares (Mohd-Basri, 2009) and the total oil palm biomass collected in 2005 was 55.73 million tonnes (Shuit al., 2009). Henson (2009) stated that the mean in situ standing plantation C stock over 30 years was 35.4 t ha-1 based on the national average biomass curve. The literature indicates that most studies have focused on the estimation of total C accumulated in oil palm standing biomass (Henson, 2009; Germer and Sauerborn, 2008; Henson, 2004; Henson, 1999; Khalid et al., 1999a, b) while information about the C stored in soil that act as C sink or source are limited. Henson (2009) suggested that no major changes of SOC occurred in oil palm plantations developed on mineral soils as compared to peat soil. Soil C responds much more slowly to the effects of land management as compared to many other soil chemical and physical properties. Significant differences in SOC concentration in various agricultural ecosystem can only be detected after at least 3 years or more typically 5 to 10 years (Kravchenko et al., 2006; Smith, 2004). However, changes in SOC over different oil palm ages are not clearly understood. It is essential to study such changes in order to understand the SOC sequestration processes within the life cycle of oil palm.
Apart from the above, significant soil variation can occur from land use and management practices, particularly in oil palm plantation that consists of several operational zones, such as Weeded Circle (WC), Frond Heap (FH) and Harvesting Path (HP). A recent study showed that WC accounted for the highest SOC content, followed by FH and HP in a 27 year-old oil palm stand (Law et al., 2009). Spatial variability of SOC within the WC, FH and HP was also reported. Spatial variability assessment provides a valuable baseline against which subsequent and future measurements can be evaluated (Liu et al., 2006). It enables faster and efficient detection of changes in SOC over a large planting area (Kravchenko et al., 2006). As such, classical and geo-spatial statistical techniques have been widely employed in order to study the spatial distribution patterns of SOC (Jian-Bing et al., 2006; Zhang and McGrath, 2004; McGrath and Zhang, 2003; Chevallier et al., 2000). Geo-spatial statistics are based on the concept of regionalized variables (Webster and Oliver, 2007). This concept provides advanced tools to quantify the spatial features of soil parameters and to carry out spatial interpolation (Liu et al., 2006). Hence, this study was aimed at quantifying the spatial variability of SOC, estimating SOC at unsampled locations using interpolation technique and comparing the spatial variability of SOC between young and mature oil palm stands.
MATERIALS AND METHODS
Site Attributes
This study was carried out in a commercial oil palm plantation located at
Port Dickson, Negeri Sembilan, Peninsular Malaysia. This study was conducted
between December 2007 and April 2008. The study site was chosen to represent
two palm age groups, i.e., 5 Years after Planting (YAP) and 17 YAP. The site
of 5 YAP is geographically located at 02° 38’ North and 101° 48’
East, with a nearly flat topography while the site of 17 YAP is 02° 37’
North and 101° 49’ East, with a gently sloping to undulating topography
(0-7% slope).
| Fig. 1: | Sampling distance within a cluster |
The total planted area of 5 YAP and 17 YAP was 42 ha and 60 ha, respectively. The soil of both sites was dominated by Tai Tak series (Typic Kandiudults). The palm densities for 5 YAP and 17 YAP were 148 palms and 144 palms ha-1, respectively, with palms spaced in a triangular pattern at a distance of 9.1x9.1x9.1 m. The annual rainfall ranges between 1700 and 2100 mm.
Sampling Protocol
A systematic design was employed for soil sampling at the 0-20 cm depth
based on a cluster of 4 palms (a, b, c and d) that included three operational
zones: Weeded Circle (WC), Frond Heap (FH) and Harvesting Path (HP), as shown
in Fig. 1. A total of 60 sampling clusters were established
for each site. The sampling area for 5 YAP and 17 YAP was 5.7 ha and 7.2 ha,
respectively. Total sampling area differed between both sites due to tree death
and existence of rocks in the sampling area, particularly at the 17 YAP site.
Thus, the sampling area for 17 YAP was larger than that of 5 YAP in order to
cater for 60 sampling clusters. Within each cluster, soil samples were obtained
from three operational zones using a core auger. Sampling points at the WC were
0.6 m apart from the palm base (a), while sampling points at the FH and HP were
both 4.0 m apart from palm base (b, d). All measurement points were geo-referenced
using a differential Global Positioning System (dGPS).
Soil Analysis
Soil samples were air-dried and ground to pass through a sieve with a 1
mm mesh. Soil samples were analyzed for pH in 1:2 soil:water suspension (Tan,
2005), total Nitrogen (N) using the Kjeldahl method (Bremmer
and Mulvaney, 1982), available Phosphorus (P) using the Bray II method (Olsen
and Sommers, 1982), exchangeable Potassium (K), Calcium (Ca) and Magnesium
(Mg) using 1N NH4OAc followed by Atomic Absorption Spectrometry and
Cation Exchange Capacity (CEC) using the leaching method (Thomas,
1982). Soil bulk density was determined using a core sampler. Pipet method
was used to determine soil texture (Gee and Bauder, 1986).
Soil C concentrations were determined using a LECO CR-412 carbon analyzer with
dry combustion method (Nelson and Sommers, 1982). Since
CaCO3 contribution was minimal, the results were interpreted as SOC.
The soil chemical and physical characteristics of the 5 YAP and 17 YAP sites
are shown in Table 1.
| Table 1: | Soil characteristics of 5 YAP and 17 YAP sites |
| WC: Weeded circle, FH: Frond Heap, HP: Harvesting Path, 1:Expressed in the following units: % (for N), μg g-1 (for P), cmol(+) kg-1 (for K, Ca, Mg and CEC) and g cm-3 (for bulk density), 2:Bulk density | |
Classical Statistical Analysis
SOC data were first analyzed using Exploratory Data Analysis (EDA), which
included descriptive statistics, normality check and outlier detection (Balasundram
et al., 2008). Descriptive statistics and normality check were computed
using Statistix Version 8.1 (Analytical Software). Normality check was performed
using the Shapiro-Wilk Test. Where necessary, non-normal data were transformed
using the appropriate function. Outlier testing was carried out using the Extreme
Studentized Deviate (ESD) method, which is also known as Grubb’ Test.
Geo-Spatial Statistical Analysis
Following the EDA, SOC data were spatially analyzed using variography and
interpolation techniques. Variography uses semivariogram to quantify and model
the spatial variance of data. The semivariogram determines the decrease in variance
between samples collected at decreasing separation distances from one another.
It can be estimated using the following formulae (Isaaks and
Srivastava, 1989):
where, h is the separation distance between location xi or xi+h, zi or zi+h are the measured values for the regionalized variable at location xi or xi+h and n(h) is the number of pairs at any separation distance h.
The semivariogram is generally fitted with an authorized model, such as spherical, exponential or Gaussian model (Webster and Oliver, 2007). These models provide information about spatial structure and spatial attributes such as nugget, sill and effective range. Nugget is defined as a measure of the amount of variance due to errors in sampling, measurement and other unexplained sources of variance. Sill is defined as the total vertical scale of the variogram that the semivariance becomes constant when the distance between sample locations increases. The distance at which the variogram reaches the sill is called the range of spatial dependence, that is, the distance at which the samples become spatially independent and uncorrelated with one another. This indicates that sampled points are no longer spatially correlated at separation distances greater than the range. These attributes are the main input parameters for interpolation.
Interpolation is used to predict values in areas that have not been sampled and is often carried out using kriging. Kriging uses the modeled variance to estimate the values between samples (Balasundram et al., 2008). Isaaks and Srivastava (1989) stated that kriging is unbiased because it tries to nullify the mean residual error and it is considered the optimal method of spatial prediction as it minimizes the variance of the error. It is obtained from the following expression:
where, z(x0) is the value to be estimated at the location of x0, z(xi) is the measured data at location neighboring the interpolation point xi, ëi is a weight factor which depends on the semivariogram model and n is the number of neighboring measured data points used for interpolation (Balasundram et al., 2007). Variography and kriging were computed using GS+ Version 7.0 (Gamma Design Software). Point kriging method was used to estimate SOC at unsampled location. Measured and kriged values were mapped using Surfer Version 8.06 (Golden Software).
Cross Validation
Cross-validation is used to assess accuracy of the interpolated values through
the comparison between estimated and actual (measured) values using information
available in the sample data set. It allows the selection of the best model
as well as compares methods of estimation and justify the use of kriging as
an estimation method (Webster and Oliver, 2007). Kriged
values were cross-validated (Balasundram et al.,
2008; Webster and Oliver, 2007; Balasundram
et al., 2006; Isaaks and Srivastava, 1989)
using the following procedure: i) Firstly, the interpolated Mean Error (ME)
should be close to zero. The ME is calculated as follows:
where, n is the number of sample points,
Secondly, the Mean Squared Error (MSE) should be less than the sample variance. The MSE is given by:
Thirdly, the ratio of the theoretical and calculated variance, called the Standardized Mean Squared Error (SMSE), should be approximately close to one. The SMSE is given by:
where, σ2is the theoretical variance.
Validation of Measured and Predicted Values
Generally, validation of measured and predicted values were done by identifying
geographical positions from a kriged map in a random manner and using those
positions to determine the sampling points at the study area. The sampled values,
also known as measured values, were cross-validated with the predicted values
obtained at the same positions from the kriged map. However, this procedure
cannot be performed in an oil palm plantation because it comprises different
operational zones with different land management practices: WC, FH and HP. All
the three operational zones have different SOC levels and should be analyzed
separately (Law et al., 2009). Therefore, only
one operational zone was presented in a kriged map. Identification of geographical
positions from the kriged map in a random manner will cause overestimation or
underestimation of SOC concentration as the selected position might not be the
actual position of the operational zone. Moreover, based on the basic principles
of kriging, a value equal to the measured value will always be generated by
kriging when interpolating at a location where a measurement exists (Mulla
and McBratney, 1999). Accordingly, an inverse validation procedure was carried
out in this study in order to validate the measured and predicted values of
SOC in oil palm.
At first, a total of 20 soil samples were collected from each operational zone at 5 YAP and 17 YAP within the study area in a random manner and were analyzed for SOC concentration. All measurement points were geo-referenced using a dGPS. Positions of those measurement points were identified and the predicted SOC values were obtained from the kriged map. Both of the measured and predicted SOC values were cross-validated. However, the exact position of measurement points could not be obtained from the kriged map because of the errors in sampling, accuracy of dGPS and accuracy of interpolation method. Nearest neighbour sampling technique was used to locate the positions that were closest to the measurement point from the kriged map. A total of 7 positions were selected for each measurement point and the predicted values were obtained for outlier detection using Grubbs’ test. The mean of these predicted values were computed and were compared with the measured value. Unpaired t-test was performed to identify the significant difference or vice versa of the measured and predicted values using Statistical Analysis System (SAS Institute, 2002). Differences were considered significant at p<0.05 and results were reported as Means±Standard Deviation.
RESULTS AND DISCUSSION
Variation among Operational Zones and Comparison between Age Groups
SOC data from the three operational zones of both age groups were examined
using Analysis of Variance (ANOVA). Mean and standard error values are given
in Fig. 2. The results showed significant differences between
operational zones for both age groups. For 5 YAP, the highest SOC content was
found in FH, followed by WC and HP, while for 17 YAP, both WC and FH were not
significantly different in SOC content and HP accounted for the lowest SOC content.
Comparatively, the study by Law et al. (2009)
showed that the highest SOC content was found in WC, followed by FH and HP in
a 27 year old palm stand.Higher SOC content was found in the FH zone as compared
to other operational zones, indicating that accumulated oil palm frond at this
area plays an essential role in soil C sequestration. Decomposition of oil palm
frond results in increasing SOC levels. Henson (2004)
noted that C sequestration in FH of oil palm was higher than other ground vegetation.
The number of fronds per palm in FH zones has a linear relationship with palm
age, where the frond number increases with palm age (Henson,
1994). The SOC content at WC increased with increase in palm stand where
it accounted for lower SOC content than FH at 5 YAP, but similar amount of SOC
with FH at 17 YAP. This indicates that root growth and turnover, as a result
of palm growth, increases SOC content (Rasse et al.,
2005). The HP of both age groups had the lowest SOC content and 17 YAP accounted
for slightly lower amount of SOC than that of 5 YAP.
| Fig. 2: | Soil Organic Carbon (SOC) content (%) at the three operational zones (WC, FH and HP) for (a) 5 YAP and (b) 17 YAP at 0-20 cm depth |
| Table 2: | Descriptive statistics of SOC (%) at the three operational zones for 5 YAP and 17 YAP |
| *: WC: Weeded Circle, FH: Frond Heap, HP: Harvesting Path, 1: numbers less than 60 indicate that non-spatial outliers were removed from the data set. Non-spatial outliers were detected using the Extreme Studentized Deviate (ESD) method, 2: Significant if the absolute value of skewness or kurtosis is > 2 times its standard error. The standard error of skewness = (6/n)0.5 while the standard error of kurtosis = (24/n)0.5, 3: Estimated using the Shapiro-Wilk test. If the test stastistic W is significant (p<0.05), the distribution is not normal (for FH at 5 YAP, FH and HP at 17 YAP, the W was computed after data transformation using natural log), ns: Not significant (p>0.05) | |
It could be mainly attributed to higher bulk density in HP at 17 YAP than that of 5 YAP, as shown in Table 1. SOC demonstrates an inverse relationship with surface bulk density (Mzuku et al., 2005). Results suggest that land management practices had a significant effect on the SOC variation. This study also demonstrated that there is a change of SOC in oil palm across different palm age and operational zones.
Classical Statistics
Descriptive statistics of SOC (%) at the three operational zones for 5 YAP
and 17 YAP are given in Table 2. SOC mean at the three operational
zones for 5 YAP and 17 YAP ranged from 1.13 to 1.76% and 1.06 to 1.86%, respectively.
The Coefficient of Variation (CV), which is defined as the ratio of standard
deviation to mean, from all three operational zones of 5 YAP and 17 YAP ranged
from 18 to 26% and 26 to 37%, respectively. Both sites indicated low to moderate
variability in data set. Based on the Shairo-Wilk statistic, data from all three
operational zones of 5 YAP and 17 YAP were normally distributed. For FH at 5
YAP, FH and HP at 17 YAP, the normality test was computed after data transformation
using natural log because they initially showed a non-normal distribution. The
coefficients of skewness and kurtosis describe the shape of the sample distribution.
All three operational zones for both sites showed a positive coefficient of
skewness, indicating that asymmetry in the sample distribution with the higher
values tailing to the right. The coefficient of kurtosis from all three operational
zones for both sites, with the exception of HP of 5 YAP, was positive, indicating
that the distribution is peaked. The sample distribution of HP of 5 YAP was
relatively flat with a negative kurtosis.
Spatial Structure and Attributes
The semivariograms of SOC for WC, FH and HP at 5 YAP and 17 YAP are shown
in Fig. 3a-f and the spatial attributes of these fitted semivariograms
are summarized in Table 3.
| Fig. 3: | Semivariograms of SOC for (a) WC, (b) FH and (c) HP at 5 YAP and (d) WC, (e) FH and (f) HP at 17 YAP |
| Table 3: | Spatial structure and attributes of SOC at three operational zones for 5 YAP and 17 YAP |
| *: WC: Weeded circle, FH: Frond heap, HP: Harvesting path, 1:S-Spherical, E-Exponential | |
For the WC, FH and HP at 5 YAP, the semivariograms of SOC were constructed based on an active lag of 187 m and lag class interval of 19, 18.7 and 19.5 m, respectively. The semivariograms of SOC for WC, FH and HP at 17 YAP were fitted with an active lag of 216 m and lag class interval of 18.5, 21.6 and 18 m, respectively.
All three operational zones from both 5 YAP and 17 YAP exhibited a definable spatial structure and were described by either spherical or exponential models. Spatial dependence was classified using nugget to sill ratio (Cambardella et al., 1994) whereby:
Based on the above definition, the nugget to sill of SOC from WC, FH and HP for 5 YAP were 0.01, 0.16 and 0.09, respectively. This indicates that the explainable proportion of the total variation of SOC from WC, FH and HP were 99, 84 and 91%, respectively, while the remaining variations can be attributed to random sources. Both WC and FH registered an Effective Range (ER) of 45.8 and 43.5 m, respectively, while HP had an ER of 66.9 m. For 17 YAP, the nugget to sill of SOC from WC, FH and HP were 0.12, 0.11 and 0.02, respectively. This indicates that the explainable proportion of the total variation of SOC from WC, FH and HP were 88, 89 and 98%, respectively, while the remaining variations can be attributed to random sources. FH exhibited the longest ER of 86.1 m, followed by WC and HP with an ER of 69.0 and 50.9 m, respectively. Comparatively, operational zones at 17 YAP exhibited a longer ER than those at 5 YAP. This indicates that SOC content at 17 YAP is more stable than that from 5 YAP.
All the operational zones exhibited strong dependence and had similar results as reported in the previous study (Law et al., 2009). It has been noted that soil C content exhibited stronger spatial dependence with relatively large spatial correlation range as compared to other soil properties (Kravchenko et al., 2006; Wang et al., 2002; Cambardella et al., 1994). ER has a great implication on sampling design. Sampling design should use separation distances that are shorter than ER because sampling points will not be subjected to spatial correlation at separation distances greater than ER. It has been recommended that spacing between sampling points should be 0.25 to 0.5 of the ER in order to capture the spatial pattern of a given property (Mulla and McBratney, 1999).
Interpolation and Spatial Variability
Interpolation was performed on SOC data that exhibited strong spatial dependence.
The natural log transformed data sets of FH at 5 YAP, FH and HP at 17 YAP were
used for interpolation and then the predicted values were back transformed to
produce the spatial distribution map of SOC concentration using the following
exponential function:
SOC = e1n SOC |
The distribution and pattern of both measured and interpolated values for each operational zone of 5 YAP and 17 YAP are illustrated as a contour map in Fig. 4.
SOC levels (%) for each operational zone were classified based on mean and Standard Deviation (SD) values (Law et al., 2009). The classes generated were as follows:
From the contour map, it is clear that SOC for all operational zones at 5 YAP and 17 YAP were variable. At 5 YAP, SOC for WC was spatially clustered with 46% of the study area showing an average value of 1.68% of SOC, 5% with values above the average and 49% with values below the average. In contrast, 41% of the study area for FH were clustered at the average value of 1.99% SOC, 4% with values above the average and 55% with values below the average. For HP, 41% of the study area were clustered at average value of 1.25% SOC, 6% with values above the average and 53% with values below the average.
At 17 YAP, SOC for WC was spatially clustered with 55% of the study area showing an average value of 2.04% of SOC, 9% with values above the average and 36% with values below the average. For FH, 37% of the study area were clustered at the average value of 2.11% SOC, 54% with values below the average, while the remaining area was above the average. At HP, 35% of the study area were clustered at the average value of 1.26% SOC, 5% with values above the average and 56% with values below the average.
Cross-validation statistics showed that all operational areas at 5 YAP and 17 YAP, with the exception of FH at 5 YAP satisfied the interpolation accuracy criteria (Table 4). MSE, that included the mean or bias and the spread of the error distribution, should be less than the sample variance (Webster and Oliver, 2007; Wackernagel, 2003). MSE of FH at 5 YAP was found to be greater than the sample variance. It could be attributed to high variation of SOC in FH that was still in the early stage of establishment. Validation of measured and predicted values of SOC at WC, FH and HP of 5 YAP and 17 YAP are given in Table 5. Results of the WC, FH and HP for both age groups were statistically insignificant. This verifies that kriging is an appropriate interpolation technique for predicting SOC levels at unsampled location.
| Table 4: | Cross-validation statistics of kriged values for SOC (%) at the three operational zones for 5 YAP and 17 YAP |
| *: WC: Weeded circle, FH: Frond heap, HP: Harvesting path, 1: Mean error , 2:Mean squared error, 3: Standardized mean squared error | |
| Table 5: | Mean±Standard Deviation (SD) and unpaired t-test for measured and predicted values of SOC (%) for WC, FH and HP at 5 YAP and 17 YAP |
| *: WC: Weeded circle, FH: Frond heap, HP: Harvesting path, ns: Not significant (p>0.05) | |
| Fig. 4: | Spatial distribution of SOC (based on measured and kriged values) for (a) WC, (b) FH and (c) HP at 5 YAP, and (d) WC, (e) FH and (f) HP at 17 YAP |
CONCLUSIONS
This study compared the spatial distribution pattern of SOC between young and mature palms, 5 Years After Planting (YAP) and 17 YAP, respectively. Spatial variability of SOC at the three operational zones: Weeded Circle (WC), Frond Heap (FH) and Harvesting Path (HP) in 5- and 17-year old oil palm stands were manifested in this study. All three operational zones at 5 YAP and 17 YAP exhibited a definable spatial structure with strong spatial dependence, which were described by either spherical or exponential models. Operational zones of 5 YAP exhibited a shorter effective range than operational zones of 17 YAP. This indicates that the distance between sampling locations for 5 YAP are closer as compared to 17 YAP. Contour maps of SOC for the three operational areas, which were constructed based on measured and kriged values, clearly showed spatial clustering of SOC values. All operational areas satisfied the interpolation accuracy criteria with the exception of FH at 5 YAP.
SOC concentration at 17 year old palms was found to be more stable than that from 5 year old palms. The results also showed that there is not much increment of SOC between young and mature palms. Differences between measured and predicted SOC values of WC, FH and HP for both age groups were statistically insignificant. As a conclusion, this study verifies that spatial variability assessment appears to be a feasible technique to quantify the variability of SOC in oil palm cultivation.
ACKNOWLEDGMENTS
The authors are grateful to Sime Darby Plantations Sdn. Bhd. for logistical support and technical collaboration. We also thank Mr. Junaidi Jaafar, Mr. Alias Tahar, Mr. Mohamad Zaidi Dan and Mr. Law Chu Chien for assisting with field work and Ms. Sarimah Hashim for help with laboratory analysis.