- Email: [email protected]
Landscape scale estimation of soil carbon stock using 3D modelling
Science of the Total Environment 487 (2014) 578–586
Contents lists available at ScienceDirect
Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv
Landscape scale estimation of soil carbon stock using 3D modelling F. Veronesi, R. Corstanje ⁎, T. Mayr Cranfield University, Bldg. 37, School of Applied Sciences, MK43 0AL Bedfordshire, United Kingdom
H I G H L I G H T S • • • •
Effective regional mapping of soil C must consider its behaviour over the soil profile. Regional estimates of C stock must consider the profile behaviour of soil bulk density. Depth functions combined with interpolation can successfully map stock in 3D. Spatial interpretation of performance identifies soil strata.
a r t i c l e
i n f o
Article history: Received 30 September 2012 Received in revised form 14 February 2014 Accepted 15 February 2014 Available online 11 March 2014 Keywords: Soil carbon 3 Dimensional maps Soil C landscapes Carbon stock estimates
a b s t r a c t Soil C is the largest pool of carbon in the terrestrial biosphere, and yet the processes of C accumulation, transformation and loss are poorly accounted for. This, in part, is due to the fact that soil C is not uniformly distributed through the soil depth profile and most current landscape level predictions of C do not adequately account the vertical distribution of soil C. In this study, we apply a method based on simple soil specific depth functions to map the soil C stock in three-dimensions at landscape scale. We used soil C and bulk density data from the Soil Survey for England and Wales to map an area in the West Midlands region of approximately 13,948 km2. We applied a method which describes the variation through the soil profile and interpolates this across the landscape using well established soil drivers such as relief, land cover and geology. The results indicate that this mapping method can effectively reproduce the observed variation in the soil profiles samples. The mapping results were validated using cross validation and an independent validation. The cross-validation resulted in an R2 of 36% for soil C and 44% for BULKD. These results are generally in line with previous validated studies. In addition, an independent validation was undertaken, comparing the predictions against the National Soil Inventory (NSI) dataset. The majority of the residuals of this validation are between ±5% of soil C. This indicates high level of accuracy in replicating topsoil values. In addition, the results were compared to a previous study estimating the carbon stock of the UK. We discuss the implications of our results within the context of soil C loss factors such as erosion and the impact on regional C process models. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Soil C is the largest pool of carbon in the terrestrial biosphere (Schlesinger and Bernhardt, 1997), it accounts for three times as much carbon compared to that available in the vegetation and twice that in the atmosphere (Schils et al., 2008). The capacity to predict the consequences of climate change on the soil C pool and act accordingly, depends upon our understanding of C distribution in the soil volume (Jobbagy and Jackson, 2000). These consequences are explained in the IPCC report (2007) which suggests a change in the precipitation pattern, with an increase of around 2–10% in the annual runoff. This means an increase in the frequency of exceptional weather events coupled with dry soils, giving rise to potential increase in the risk of soil erosion. For this reason an increasing percentage of soil C stocked in the very top part of the soil profile can be at risk of being lost, impacting the fertility ⁎ Corresponding author.
http://dx.doi.org/10.1016/j.scitotenv.2014.02.061 0048-9697/© 2014 Elsevier B.V. All rights reserved.
of the soils which are already highly stressed by years of intensive agriculture practices. In order to assess which areas will be most affected by these future climate pattern changes, there is a need for precisely predicting soil C stocks through the depth profile and at fine spatial resolution. However, most of the available soil C stock maps, including the UK map (Bradley et al., 2005, provide only average estimates of C available in topsoil (from 0 to 30 cm) and subsoil (from 30 to 100 cm). The spatial changes in vertical distribution of soil C is an aspect still poorly represented (Jobbagy and Jackson, 2000; Gifford and Roderick, 2003)). There have been various technical advances in estimating spatial C stocks in the soil profile. For instance, Ellert and Bettany (1995) developed the use of genetic horizons as reference levels by using an ‘equivalent soil mass’ concept. The sampling for soil C (and other elements) was carried out by depth increments corresponding to soil horizons. Zan et al. (2001) modified the approach of Ellert and Bettany (1995), replacing genetic horizons with depth segments
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
such an s 15 cm layers. Both suffer from the problem that both bulk density and soil C continuously vary down the profile rather being uniform within the genetic horizon or each 15 cm segment. Grimm et al. (2008) developed a 3D model for predicting soil organic carbon a fixed depth intervals using digital soil mapping techniques and soil units as the basis for spatialisation. Minasny et al. (2013) provided a comprehensive review of modelling organic carbon and demonstrated their various approaches with case studies. For the Edgeroi area of Australia they used ta combination of equal-area splines and a spatial prediction model to map the whole are for bulk density and organic carbon. In a typical soil profile, the first centimetres can have a very high organic matter content (Batjes, 1996), derived from plant and animal residues. With depth the percentage in volume of organic carbon decreases non-linearly with a pattern that is frequently modelled with exponential functions (Russell and Moore, 1968; Hilinski, 2001). There are also other equations proposed to describe the decrease of soil carbon with depth, but they are just a variance of the exponential model (Arrouays and Pelissier, 1994; Bernoux et al., 1998). Minasny et al. (2013) used a generalized negative exponential depth function. However, in very wet soils decay is inhibited, causing the accumulation of organic matter and the formation of peat, which can be found either continuously throughout the soil profile, or as a buried horizon. In either case, C in these soils exhibits a different behaviour by depth than the previously described typical soil. For this reason, considering only volume averages cannot be sufficient to accurately assess the potential changes driven by climate change. There are broadly three pathways through which C is lost from soil; through erosion, leaching (hydrological) and respiration (Dawson and Smith, 2007). Each of these is sensitive or depends on the vertical distribution of soil C. Soil C loss through water erosion can occur in two ways, through sheet flow over the surface of the soil or through channel erosion. Overland flow removes top soil and changes the distribution of soil C through the profile. The amount lost through erosion also directly depends on the vertical distribution of soil C. Hydrological loss of C from soil is a function of the movement of C from soil OM to the soil water, which is then lost through leaching as discharge from soil to the surface water. Dawson and Smith (2007) describe numerous studies in which
579
significant correlations between soil C and dissolved organic C in UK catchments are observed; movement and loss of C through this pathway depends on the spatial and vertical distribution of soil C (Grieve, 1990; Hinton et al., 1998). Peak C loss tends to coincide with peak discharge events, in which the 15–20 cm surface layer dominates flow characteristics (Worrall et al., 2003). This would suggest that the vertical and horizontal distribution of soil C determines both the concentration and also flow characteristics of C loss from soil. The largest amount of C is lost through organic matter mineralization. These losses are typically modelled using spatially distributed, integrated land–atmosphere process models such as JULES (Harrison et al., 2008). Improvements to the vertical distribution of soil C component in these types of models increases their response-sensitivity to changes in soil stocks and processes. Jones et al. (2005) compared the outputs from a generic Soil C model with measurements of topsoil (organic) C aggregated spatially according to both soil mapping and administration units. In a comprehensive review by Minasny et al. (2013), the outcome of various studies are tabulated and where, in general, R2 for prediction accuracy range from 0.21 to 0.98 for soil C prediction from areal imagery, the statistics obtained by cross-validation, internal or external validation. There is only one study in this list in which the C stock predictions are independently validated (R2 of 0.41). In this paper we consider a method developed previously (Veronesi et al., 2012) for in field mapping of soil compaction, and we consider this now for regional scale C stocks. The method is developed and applied within an area of the UK with a high density of soil profile observations, but validated using an independent topsoil dataset sampled on a 5 km grid. The method creates a fine vertical resolution 3D map of the soil C stock, describing its behaviour both vertically and spatially. 2. Materials and methods 2.1. Study area The study area is located in the West Midlands region, extending approximately 13,948 km2 (Fig. 1). Within the pilot area, there are 118
Fig. 1. Map of the area under study. The area extends for 13,948 km2 over the West Midlands region and on the Welsh border.
580
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
different soil associations according to the National Soil Map (NATMAP) of England and Wales with the Bromyard Association having the largest extend (8.51%). The area is covered largely by loamy soils, although the silt content gets larger towards the East and North whereas the clay content gets larger towards the South. According to the Centre for Ecology & Hydrology (CEH) Land Cover Map 2000 (Fuller et al., 2000) around one third of the area is covered by arable land, which acts as a C source (Schils et al., 2008); 50% is covered by natural, semi-natural or improved grassland and 11% is urban or suburban areas. The remaining land is classified as woodland, inland water or heath. 2.2. Dataset description 2.2.1. Training data The analytical data used for this study came from the Representative Profile data collected by the Soil Survey of England and Wales, now the National Soil Resources Institute (NSRI) to support national mapping programmes. Representative Profiles were sampled by horizons as identified in the field. The data were collected between 1969 and 1985 with more than 90% covering the period 1973 to 1983. A total of 435 profiles were available for the study areas. In addition to the temporal variation, the locations were selected to describe model profiles for individual soil series used in the national soil classification, consequently, the samples are not very well distributed. There are clusters of high density sampling and a large area in the West, which coincides with the Welsh border, with few profiles (see Fig. 1). This is very often characteristic for legacy data and highlights their limitation as they were not collected for quantitative soil mapping approaches. However, they represent the only available national dataset.
2.4.1. 1st stage At first the focus is on the selection of the depth functions that best describe the variation by depth of the two soil properties of interest (soil C and BULKD), basing the choice on a framework for testing different types of mathematical equations followed by identifying the most appropriate number of coefficients. The analysis was undertaken using DataFit from Oakland Engineering which allows to investigate a wide range mathematical equations. For the vertical distribution of soil C we considered the general category of exponential functions because they are known from the literature as being able to describe the variations by depth of soil C (Gifford and Roderick, 2003; Hilinski, 2001; Hiederer, 2009; Minasny et al., 2013). For BULKD we noticed that the general pattern was also exponential. From the analysis, two equations were selected: – A Power function for describing the vertical distribution of soil C: β
C i ¼ β 1 di 2
ð1Þ
where Ci is the value of C at the ith depth interval and di is the depth. – An inverse power function for describing the vertical distribution of soil BULKD:
BULKDi ¼
1 β
β1 di 2
ð2Þ
where BULKDi is the value of BULKD at the ith depth interval and di is the depth.
2.3. Environmental covariates In order to improve the accuracy of the soil C and BULKD predictions we implemented four grids of environmental covariates, namely: – Digital elevation model (DEM). A 50 m resolution DEM, created by the Ordnance Survey and constructed from the 1:50,000 topographic map. From this five derivatives were obtained: slope, aspect, total curvature, profile curvature and plan curvature. – Soil association map from NSRI. A 1 km grid map with information about the predominant soil association in each grid cell. – Land Cover Map 2000. Available from CEH, it has a resolution of 1 km with information about the predominant usage of the land. – Digital Geological Map (DigMap 250) from the British Geological Survey at the scale of 1:250,000.
2.3.1. Validation data The National Soil Inventory (NSI) (Loveland, 1990) was established to obtain an unbiased estimate of the distribution of the soils of England and Wales and of the chemistry of the topsoil (0–15 cm depth). It followed a non-purposive design on a 5-km grid covering all soil types, habitats and land uses. The original sampling (1978–83) was repeated on approx. 40% of sites 12–25 years later in such a way as to detect a specified change in the soil with specified confidence. Samples of air-dried soil from both samplings were archived. The NSI is unique among national soil inventories globally in having been comprehensively resampled. 2.4. Methodology The top-down mapping method (Veronesi et al., 2012) relies on soil specific depth functions to describe the variation of soil properties with depth. Its application on the soil C and BULKD dataset can be summarized in three stages:
2.4.2. 2nd stage In the second stage of the mapping process, we fitted Eqs. (1) and (2) to each available soil profile, computing the local R2 value. From the goodness of fit analysis we produced two R2 maps. 2.4.3. 3rd stage The final step of the mapping process involves applying a universal co-kriging algorithm to estimate the value of the two coefficients of the depth functions, on a regular grid of 1 km of horizontal resolution, assuming no variation within the grid cell. Universal kriging was adopted because at this scale the assumption of intrinsic stationarity cannot hold and the soil process needs to be considered as non-stationary. This interpolation method is able to work under these circumstances by modelling the trend in the variable using a multi-linear regression between the soil property and the environmental covariates and then creating the variogram from the residuals of the regression (Hengl et al., 2007). Subsequently the depth functions are solved for each interpolated point in order to create a lattice of 3D estimates with a 10 cm vertical resolution. 2.5. Validation The validation was undertaken using two methods (Minasny et al., 2013): cross-validation and an independent validation against the NSI data. The cross-validation approach relies on excluding a set percentage of soil profiles (5, 10, 20, 30%) and re-predicting their values. The observed and predicted values are then compared, calculating the R2. The second method is an independent validation accomplished by comparing the soil C percentage available in the NSI dataset, with the estimated soil C percentage. To assess the goodness of fit, the residuals between the observed values in the NSI dataset and the predicted values with the top-down mapping method were computed.
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
581
Fig. 2. General profiles of the two soil properties and depth functions fitting. Red crosses signify the mean value of the soil property in each depth layer.
Subsequently, the C stock was computed with Eq. (4) (Xu et al., 2011):
2.6. Soil C stock estimation The soil C density was calculated using the following equation, suggested by (Xu et al., 2011):
C Stock ¼ C Density ¼ C i BULKDi di
ð3Þ
where ‘C Density’ (t/ha) is the soil carbon density; Ci is the C (%) content of the soil layer i, BULKDi is the BULKD (g/cm3) value of the soil layer i and di is the depth of the soil layer i (cm). We then overlaid the resulting density estimations with the Land Cover Map 2000 from CEH (Fuller et al., 2000) and we made the following corrections, as suggested in Bradley et al. (2005) and Xu et al. (2011): – For urban areas and water bodies the density value was set to 0. – For semi urban areas the value was set to one half of the estimated value.
X C Density Area 1010
ð4Þ
where ‘C Stock’ (Tg) is the carbon stock value estimated after the corrections cited above and ‘Area’ is the area of a single prediction grid cell in m2. We calculated the C stock value for each depth interval estimated during the mapping process, namely: 0–10; 10–20; 20–30; 30–40; 40–50; 50–60; 60–70; 70–80; 80–90; and 90–100 cm. We also computed average values for topsoil (0–30 cm) and subsoil (30–100 cm). 2.7. C density variation estimation We estimated the uncertainty related to each point on the prediction grid adapting a method proposed by Goidts et al. (2009) and Schrumpf et al. (2011), for assessing the error propagation of each variable in the C
Fig. 3. R2 values for soil C and bulk density within the study area.
582
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
Table 1 Descriptive statistics of the observed and predicted datasets. Topsoil refers to the 0 to 30 cm soil layer from the surface; the subsoil refers to 0 to 100 cm soil layer. Observed dataset Min
Max
Mean
Median
St.Dev.
0.30 0.10 0.59 0.61
43.00 39.10 1.76 1.82
3.06 1.01 1.23 1.45
2.30 0.60 1.25 1.47
3.23 2.73 0.23 0.16
Predicted with top-down mapping SOC Topsoil 0.00 (%) Subsoil 0.00 BULKD Topsoil 0.27 (g cm3) Subsoil 0.74
45.75 68.08 2.82 2.49
4.36 0.45 1.07 1.35
2.95 0.24 1.14 1.35
4.24 0.89 0.26 0.11
SOC (%) BULKD (g cm3)
Topsoil Subsoil Topsoil Subsoil
Density estimation. The variance of the C Density assessment is based on the following equation:
2
VarðC DensityÞ ¼ ðC DensityÞ
! 2 2 σC σ BULKD σ C−BULKD þ þ 2 C BULKD C 2 BULKD2
ð5Þ
where Var(C Density) is the error associated with the density prediction, in t/ha; σC and σBULKD denote the standard deviation of carbon content and bulk density; and σC–BULKD is the covariance of C and BULKD. We applied this equation to two depth intervals, namely the 0–30 cm layer (topsoil) and the 0–100 cm layer (subsoil). 3. Results 3.1. Landscape level soil C maps
curves were fitted to the entire dataset, still yielded R2 values of 0.48 and 0.59 respectively. In order to numerically demonstrate the accuracy of these two equations in describing the spatial pattern of soil C and BULKD variation by depth, we cross-validated the fits to each individual profile, and obtained a profile R2 (Fig. 3). The average profile R2 for was 0.92 for soil C and 0.85 for BULKD, with the majority of profile R2 in the range of 0.90 to 0.99 for soil C and 0.80 to 0.99 for BULKD respectively. We classed profiles which were not adequately described by the depth functions (Eqs. (1) and (2)) as those that had R2 of less than 0.5 (i.e. the function, on the whole, did not describe the majority of the variation by depth). We observed no such R2 in soil C, which would suggest that in this area and at this scale the soil C profiles present a relative homogeneous pattern. For BULKD there are only two areas which present R2 values below 0.5, present in the more southern part of the study area. This is due to the fact that three profiles in that area have a linear depth pattern that the inverse power function struggles to describe. Normally, bulk density increases with depth but in these three cases it remains almost constant. These differences are present in only three samples. This is a key property of the depth function mapping method (Veronesi et al., 2012). The fact that the selected depth functions can describe only a small set of geometries means that, from the goodness of fit map, we can extract very useful information about the 3D spatial variability in the study area. In this case, we can identify an area of linear increase in BULKD values along the profile. This way this method can be used as diagnostic tool in order to detect changes in the soil profile that otherwise would be identifiable only by looking at profiles individually. We subsequently estimated the value of soil C and BULKD with a universal co-kriging algorithm. The final 3D grid of predictions has a horizontal resolution of 1 km and vertical resolution of 10 cm. The total number of points in the prediction grid was 220,928. The
The two depth functions for describing the soil profiles of C and BULKD, were selected because they can describe the majority of the profiles with an R2 above 0.6 having only two coefficients (see Eqs. (1) and (2)). These two aspects are equally important in the definition of the best function as a high number of coefficients can potentially result in a propagation of the error during the interpolation phase, decreasing the overall accuracy. The validity of the two depth functions in describing the general pattern of the two soil properties is visualised in Fig. 2, where the red crosses represent the average value at each depth layer and the function
Fig. 4. Distribution of the residuals between the soil values observed in the NSI dataset and the values predicted in this study.
Fig. 5. Spatial pattern of the residuals of the validation against the NSI data.
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
583
Fig. 6. Soil C stock maps.
computational time for the entire mapping process was around 5 min, which is a fraction of the time needed to estimate the same number of points using 3D geostatistics alone, which for this number of points can take much longer (Kerry and Hawick, 1998). To investigate the results of the mapping, the summary statistics of the original dataset were compared with the predicted datasets (Table 1). For the topsoil (0–30 cm) the predicted dataset slightly overestimates the values of both C and BULKD and underestimates the variance. The mean C value in the observed data for the topsoil is 3.06 with a standard deviation of 3.23, compared to a predicted mean of 4.36, with a standard deviation of 4.24. Regarding the BULKD values, the observed mean for the topsoil is 1.23, with a standard deviation of 0.23, while the predicted mean is 1.07, with a standard deviation of 0.26. On the other hand for the subsoil (0–100 cm), the results are less accurate. The C observed mean is 1.01, with a standard deviation of 2.73, while the estimated mean is 0.24, with a standard deviation of 0.89. For BULKD, the subsoil results are much better. The observed mean is 1.45, with a standard deviation of 0.16, while the estimated mean is 1.35, with a standard deviation of 0.11. These results indicate that the
Table 2 Carbon stock estimates divided by land use. Topsoil refers to the 0 to 30 cm soil layer from the surface; the subsoil refers to 0 to 100 cm soil layer.
Agricultural land Grassland (natural, semi-natural and improved) Woodland
Topsoil C stock (Tg)
Subsoil C stock (Tg)
20.68 ± 0.17 88.9 ± 0.38
54.42 ± 0.22 33.23 ± 0.25
8.02 ± 0.31
2.84 ± 0.16
top-down mapping method can reproduce the general distribution of the soil dataset very well, even at landscape scale. The topsoil estimates appear to be more accurate than subsoil predictions. This means that the two depth functions are very accurate in describing the topsoil spatial variability and less accurate in lower horizons. This is probably due to a decrease in the number of profiles with data at lower depths, thereby decreasing the sample support at the lower depths. 3.2. Validation The accuracy of the mapping method was then tested with two methods: cross-validation and independent validation. The cross-validation results indicate an average value of R2 for the lowest exclusion percentage of 36% for C, with a standard deviation of 0.20; and 44% for BULKD, with a standard deviation of 0.15. These results may appear relatively low, but they are in line with previous validated research of 3D mapping at the catchment or landscape scale (Malone et al., 2009; Kempen et al., 2011; Minasny et al., 2013). The problem with mapping large areas using legacy soil observations is that the dataset was not sampled with probabilistic designs, making it less suitable for assessing the accuracy of maps. In this study, the majority of the area around Birmingham is well covered by clusters of high-density data, while the remaining area was only sparsely sampled. It is therefore difficult to consider cross-validation a sound method for validating the map, as when an excluded sample is re-predicted using other observations tens of kilometres apart, it is clear the estimation may substantially differ from the observation. For this reason we also performed an independent validation, comparing the topsoil C prediction against the data in the NSI dataset. To
584
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
assess the accuracy of the prediction, we calculated the residuals between the observed and the predicted values. Fig. 4 contains the distribution of the residuals, while in Fig. 5 they are mapped. The results indicate that the residuals are generally very low, with a mean C of 0.59% and a standard deviation of 4.46%. In Fig. 4, it is possible to appreciate that the majority of the residuals are in the interval between −5% and + 5%, even though there is a small number of large percentage values, which in Fig. 5, these higher values are in areas poorly sampled in the representative profiles data. A plausible explanation for this is that there are areas of probable marshes where the topsoil C content is particularly high, but these are almost impossible to predict if those areas were not originally sampled.
relatively low level of uncertainty for the soil C predictions. In particular the majority of the topsoil predictions have an associated variance between 0 and 0.19 t/ha with peaks of 0.95 t/ha in the Welsh border area, where the sample density is very low. Moreover, the subsoil uncertainty estimations are also low, with the majority of the point uncertainties between 0 and 0.87 t/ha, and a peak on the Welsh border of 5.07 t/ha. These results corroborate the data presented in the previous sections. In the topsoil, the top-down mapping method can achieve a generally higher level of accuracy. On the other hand, from the data described in the previous sections it is impossible to give a precise indication of the accuracy of the method in the subsoil. However, this uncertainty estimation gives us a clear indication of the overall accuracy of the mapping process, which can achieve good results also in the subsoil.
3.3. Stock estimation 4. Discussion The soil C stock was calculated and corrected using the procedure described above. The results indicate that the C stock in the topsoil is 162.91 ± 3.63 × 10−5 Tg, and the C stocked in the subsoil is 66.71 ± 3.52 × 10−5 Tg for the study area respectively. The C stock estimates are presented in the maps in Fig. 6. It appears that on the higher ground on the Welsh border, the topsoil has a relatively higher C content, due the fact that the majority of this area is grassland. By contrast, the subsoil map presents a picture where the higher C content is in the arable lands around Birmingham, and this can probably be a result of leaching of C from the arable topsoil as suggested by Bellamy et al. (2005). The C stock was also calculated for land use categories and the results are shown in Table 2. The uncertainties of the C Density point estimations were computed with the approach described in the Materials and methods section. The results are presented in Fig. 7. From this figure, we can appreciate the
In general, we can divide the study area into two distinct zones: the lowlands around Birmingham, where most of the arable and urban areas are concentrated, and the higher ground in the northern and eastern parts. These two zones are characterized by almost opposite soil profile patterns, where the lowlands exhibit low C stocks in the topsoil and relatively high stocks in the subsoil, and the higher ground presents relatively high C stocks in the topsoil with a general decrease with depth. The cause for this pattern, in the lowland areas, is the land use and land management, as the predominant use is arable agriculture. Increased soil organic matter decomposition rates and top soil erosion depleted the topsoil of OM, instead leaving the subsoil untouched or enriched with C from leaching (Dawson and Smith, 2007). On the higher ground, the density of arable and urban land is much lower and this has left these areas in an almost natural state, with the
Fig. 7. Soil C stocks uncertainty maps.
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
585
Fig. 8. Soil C percentage by volume on agricultural land.
higher concentration of C in the topsoil. This is a possible risk factor particularly for the Welsh border, because the soils there have a high silt percentage and this makes them vulnerable to erosion, as mentioned in the report of the PESERA pan-European erosion model (Kirkby et al., 2004). The fact that in these areas the value of C stock is higher indicates a high risk of C losses by erosion and it would seem to suggest that erosion protection land management practices would be appropriate. Landscape level estimates of soil C stocks are usually based on soil surveys which were not designed for assess C stocks, changes and drivers (Schrumpf et al., 2011). Samples are taken at locations in the landscape associated to representative soil characteristics and at depths representing soil horizons. As others (Malone et al., 2009) the method presented in this paper is capable of taking horizon based depths, transforming this to a continuous description of the depth variation and interpolating this across the landscape. For this landscape, we compare these estimates to the currently used map of the C stock for the United Kingdom created by Bradley et al. (2005) and our results are higher for the study area. Bradley et al. (2005) calculated a topsoil stock of 91.33 Tg and subsoil stock of 48.66 Tg compared to our estimates of C stocks of 162.91 ± 3.63 × 10−5 Tg in the topsoil, and 66.71 ± 3.52 × 10 −5 Tg in the subsoil. Bradley et al. (2005) used an approach based on pedo-transfer functions, and different soil datasets: the NSI data for the topsoil, extrapolating the soil C value from 15 to 30 cm depth because the NSI data are averages of the first 15 cm, and the soil series map to extrapolate the C values for the subsoil. The method in essence produces aggregate means for soil polygons, which seems to underestimate soil C stocks.
If we consider this map as input to soil C transformation and transport models, i.e. respiratory responses such as JULES, erosion or hydrological losses, then these results would suggest that general behaviour of C and bulk density by depth is fairly consistent across this type of landscape (dominated by arable land use) and approaches such as that described by Ungaro et al. (2010), in which a mean stock is obtained by a general depth function and subsequently interpolated are adequate to describe C stocks and C concentrations at a landscape scale. We could also consider generating slices if the C transformation models need a more specified vertical distribution; and this method, as well as equal area splines as described by Malone et al. (2009), would be adequate. The advantage of the method described in this paper is that it contains an element of diagnosis; where the dominant function does not describe the depth variation adequately, stratification of the landscape can be considered, with different depth functions for different portions of the landscape (e.g. arable vs peat soils). The mapped soil C can also be considered as degradation maps, where at less than 2%, soil C volume is so low that there is a significant decline in soil structural stability and fertility (Loveland and Webb, 2003). From the mapping results, we extrapolated the predicted points on the land defined as in arable use and then we computed the percentage of C stock by volume; the resulting map is presented in Fig. 8. Fig. 8 shows a scenario where the majority of the arable land has values of C below or slightly above the 2% threshold, meaning they are currently at very high risk of loss of fertility and soil structure, i.e. a risk map of loss of soil functions due to extremely low levels of soil C and therefore do not function effectively in the terrestrial C cycle. Conversely, these are areas in which, through changes in land management practices (such as leaving straw
586
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
and stubble on the fields), significant gains can be achieved increasing soil C stocks in the future. 5. Conclusions This method of 3D mapping with depth functions, can achieve a high level of accuracy, in particular for estimating the first 30 cm of soil. This indicates that this map can be readily implemented in soil erosion model, for example, incr'easing the amount of information in the very top part of the profile. As to the subsoil, the uncertainty estimation indicates a good level of accuracy, with the majority of the soil C stock values presenting a variance between 0 and 0.87 t/ha. On the one hand, by analyzing the soil patterns and comparing them with the land use and the soil erosion model, it is possible to conclude that this area is at high risk of soil C losses. In the higher ground on the Welsh border the soils, which have a high concentration of C, are prone to erosion because of their high silt content. On the other hand, from the risk assessment it is possible to conclude that the majority of the arable land around Birmingham is in danger of losing soil functionality such as fertility, as their percentage of C is below the 2% by volume threshold. References Arrouays D, Pelissier P. Modeling carbon storage profiles in temperate forest humic loamy soils of France. Soil Sci 1994;157:185–92. Bernoux M, Arrouays D, Volkoff B, Cerri CC, Jolivet C. Bulk densities of Brazilian Amazon soils related to other soil properties. Soil Sci Soc Am J 1998;62:743–9. Batjes NH. Total carbon and nitrogen in the soils of the world. Eur J Soil Sci 1996;47: 151–63. Bellamy PH, Loveland PJ, Bradley I, Lark RM, Kirk G. Carbon losses from all soils across England and Wales 1978–2003. Nature 2005;437:245–8. Bradley RI, Milne R, Bell J, Lilly A, Jordan C, Higgins A. A soil carbon and land use database for the United Kingdom. Soil Use Manag 2005;21:363–9. Dawson JJC, Smith P. Carbon losses from soil and its consequences for land-use management. Sci Total Environ 2007;382:165–90. Ellert BH, Bettany JR. Calculation of organic matter nutrients stored in soils under contrasting management regimes. Can J Soil Sci 1995;75:529–38. Fuller RM, Smith GM, Sanderson JM, Hill RA, Thomson AG, Cox R, Brown NJ, Clarke RT, Rothery P, Gerard FF. Land Cover Map 2000: a guide to the classification system, Centre for Ecology and Hydrology, 2002. http://www.google.co.uk/url?sa=t&rct= j&q=&esrc=s&frm=1&source=web&cd=1&ved=0CC8QFjAA&url=http%3A%2F% 2Fwww.ceh.ac.uk%2Fdocuments%2Flcm2000_classification_guide.pdf&ei=wZ2UUqy_ Foi9ygOsr4GwAQ&usg=AFQjCNFtcxS7T5TDTkD7YuTma0dKrH5GLA&bvm=bv. 57155469,d.bGQ, 2000. Gifford RM, Roderick ML. Soil carbon stocks and bulk density: spatial or cumulative mass coordinates as a basis of expression. Glob Chang Biol 2003;9:1507–14. Goidts E, Van Wesemael B, Crucifix M. Magnitude and sources of uncertainties in soil organic carbon (SOC) stock assessments at various scales. Eur J Soil Sci 2009;60:723–39. Grieve IC. Seasonal, hydrological and land management factors controlling dissolved organic carbon concentrations in the Loch Fleet catchments, southwest Scotland. Hydrol Process 1990;4:231–9. Grimm R, Behrens T, Märker M, Elsenbeer H. Soil organic carbon concentrations and stocks on Barro Colorado Island — digital soil mapping using Random Forests analysis. Geoderma 2008;146:102–13.
Harrison RG, Jones CD, Hughes JK. Competing roles of rising CO2 and climate change in the contemporary European carbon balance. Biogeosciences 2008;5:1–10. Hengl T, Heuvelink G, Rossiter DG. About regression-kriging: from equations to case studies. Comp Geosci 2007;33:1301–15. Hiederer R. Distribution of organic carbon in soil profile data. EUR 23980 EN. Luxembourg: Office for Official Publications of the European Communities; 2009. Hilinski TE. Implementation of exponential depth distribution of organic carbon in the CENTURY model. Dept Soil Crop Sci, Colorado State Uni; 2001 [http://www.nrel. colostate.edu/projects/century5/reference/html/Century/exp-c-distrib.htm]. Hinton MJ, Schiff SL, English MC. Sources and flowpaths of dissolved organic carbon during storms in two forested watersheds of the Precambrian Shield. Biogeochemistry 1998;41:175–97. Intergovernmental Panel on Climate Change — IPCCClimate change: synthesis report, 2007; 2007. Jobbagy E, Jackson R. The vertical distribution of soil organic carbon and its relation to climate and vegetation. Ecol Appl 2000;10:423–36. Jones RJA, Hiederer R, Rusco E, Montanarella L. Estimating organic carbon in the soils of Europe for policy support. Eur J Soil Sci 2005;56:655–71. Kempen B, Brus DJ, Stoorvogel JJ. Three-dimensional mapping of soil organic matter content using soil type — specific depth functions. Geoderma 2011;162:107–23. Kerry KA, Hawick KE. Kriging interpolation on high-performance computers. In: Sloot P, Bubak M, Hertzberger B, editors. High-performance computing and networks, International Conference and Exhibition Amsterdam, The Netherlands, April 21–23; 1998. p. 429–38. Kirkby MJ, Jones RJA, Irvine B, Gobin A, Govers G, Cerdain O, Van Rompaey AJJ, Le Bissonnais Y, Daroussin J, King D, Montanarella L, Grimm M, Vieillefont V, Puigdefabregas J, Boer M, Kosmas C, Yassoglou N, Tsara M, Mantel S, Van Lynden GJ, Huting J. Pan-European soil erosion risk assessment: the PESERA map. Special Publication Ispra, vol. 73; 2004 [no. S.P.I.04.73]. Loveland P, Webb J. Is there a critical level of organic matter in the agricultural soils of temperate regions: a review. Soil Tillage Res 2003;70:1–18. Loveland PJ. The National Soil Inventory of England and Wales. In: Lieth H, Markert B, editors. Element concentration cadasters in ecosystems. Weinheim, Germany: VCH; 1990. p. 73–80. Malone BP, McBratney AB, Minasny B, Laslett GM. Mapping continuous depth functions of soil carbon storage and available water capacity. Geoderma 2009;154:138–52. Minasny B, McBratney AB, Malone BP, Wheeler I. Digital mapping of soil carbon. Adv Agron 2013;118:1–47. Russell JS, Moore AW. Comparison of different depth weightings in the numerical analysis of anisotropic soil profile data. Transactions of the 9th International Congress of Soil Science 1968;4:205–13. Schils R, Kuikman P, Liski J, Van Oijen M, Smith P, Webb J, Alm J, Somogyi Z, Van den Akker J, Billett M, Emmett B, Evans C, Lindner M, Palosuo T, Bellamy P, Jandl R, Hiederer A. Review of existing information on the interrelations between soil and climate change. CLIMSOIL final report; 2008. Schlesinger WH, Bernhardt ES. Biogeochemistry, an analysis of global change. San Diego, California, USA: Academic Press; 1997. Schrumpf M, Schulze ED, Kaiser K, Schumacher J. How accurately can soil organic carbon stocks and stock changes be quantified by soil inventories? Biogeosciences 2011;8: 1193–212. Ungaro F, Staffilani F, Tarocco P. Assessing and mapping topsoil organic carbon stock at regional scale: a simple kriging approach conditional on soil map delineation and land use. Land Degrad Dev 2010;21:565–81. Veronesi F, Corstanje R, Mayr T. Mapping soil compaction in 3D with depth functions. Soil Tillage Res 2012;124:111–8. Worrall F, Burt T, Adamson J. Controls on the chemistry of runoff from an upland peat catchment. Hydrol Process 2003;17:2063–83. Xu X, Liu W, Zhang C, Kiely G. Estimation of soil organic carbon stock and its spatial distribution in the Republic of Ireland. Soil Use Manag 2011;27:156–62. Zan CS, Fyles JW, Girouard P, Samson RA. Carbon sequestration in perennial bioenergy, annual corn and uncultivated systems in southern Quebec. Agric Ecosyst Environ 2001;86:135–44.
Contents lists available at ScienceDirect
Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv
Landscape scale estimation of soil carbon stock using 3D modelling F. Veronesi, R. Corstanje ⁎, T. Mayr Cranfield University, Bldg. 37, School of Applied Sciences, MK43 0AL Bedfordshire, United Kingdom
H I G H L I G H T S • • • •
Effective regional mapping of soil C must consider its behaviour over the soil profile. Regional estimates of C stock must consider the profile behaviour of soil bulk density. Depth functions combined with interpolation can successfully map stock in 3D. Spatial interpretation of performance identifies soil strata.
a r t i c l e
i n f o
Article history: Received 30 September 2012 Received in revised form 14 February 2014 Accepted 15 February 2014 Available online 11 March 2014 Keywords: Soil carbon 3 Dimensional maps Soil C landscapes Carbon stock estimates
a b s t r a c t Soil C is the largest pool of carbon in the terrestrial biosphere, and yet the processes of C accumulation, transformation and loss are poorly accounted for. This, in part, is due to the fact that soil C is not uniformly distributed through the soil depth profile and most current landscape level predictions of C do not adequately account the vertical distribution of soil C. In this study, we apply a method based on simple soil specific depth functions to map the soil C stock in three-dimensions at landscape scale. We used soil C and bulk density data from the Soil Survey for England and Wales to map an area in the West Midlands region of approximately 13,948 km2. We applied a method which describes the variation through the soil profile and interpolates this across the landscape using well established soil drivers such as relief, land cover and geology. The results indicate that this mapping method can effectively reproduce the observed variation in the soil profiles samples. The mapping results were validated using cross validation and an independent validation. The cross-validation resulted in an R2 of 36% for soil C and 44% for BULKD. These results are generally in line with previous validated studies. In addition, an independent validation was undertaken, comparing the predictions against the National Soil Inventory (NSI) dataset. The majority of the residuals of this validation are between ±5% of soil C. This indicates high level of accuracy in replicating topsoil values. In addition, the results were compared to a previous study estimating the carbon stock of the UK. We discuss the implications of our results within the context of soil C loss factors such as erosion and the impact on regional C process models. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Soil C is the largest pool of carbon in the terrestrial biosphere (Schlesinger and Bernhardt, 1997), it accounts for three times as much carbon compared to that available in the vegetation and twice that in the atmosphere (Schils et al., 2008). The capacity to predict the consequences of climate change on the soil C pool and act accordingly, depends upon our understanding of C distribution in the soil volume (Jobbagy and Jackson, 2000). These consequences are explained in the IPCC report (2007) which suggests a change in the precipitation pattern, with an increase of around 2–10% in the annual runoff. This means an increase in the frequency of exceptional weather events coupled with dry soils, giving rise to potential increase in the risk of soil erosion. For this reason an increasing percentage of soil C stocked in the very top part of the soil profile can be at risk of being lost, impacting the fertility ⁎ Corresponding author.
http://dx.doi.org/10.1016/j.scitotenv.2014.02.061 0048-9697/© 2014 Elsevier B.V. All rights reserved.
of the soils which are already highly stressed by years of intensive agriculture practices. In order to assess which areas will be most affected by these future climate pattern changes, there is a need for precisely predicting soil C stocks through the depth profile and at fine spatial resolution. However, most of the available soil C stock maps, including the UK map (Bradley et al., 2005, provide only average estimates of C available in topsoil (from 0 to 30 cm) and subsoil (from 30 to 100 cm). The spatial changes in vertical distribution of soil C is an aspect still poorly represented (Jobbagy and Jackson, 2000; Gifford and Roderick, 2003)). There have been various technical advances in estimating spatial C stocks in the soil profile. For instance, Ellert and Bettany (1995) developed the use of genetic horizons as reference levels by using an ‘equivalent soil mass’ concept. The sampling for soil C (and other elements) was carried out by depth increments corresponding to soil horizons. Zan et al. (2001) modified the approach of Ellert and Bettany (1995), replacing genetic horizons with depth segments
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
such an s 15 cm layers. Both suffer from the problem that both bulk density and soil C continuously vary down the profile rather being uniform within the genetic horizon or each 15 cm segment. Grimm et al. (2008) developed a 3D model for predicting soil organic carbon a fixed depth intervals using digital soil mapping techniques and soil units as the basis for spatialisation. Minasny et al. (2013) provided a comprehensive review of modelling organic carbon and demonstrated their various approaches with case studies. For the Edgeroi area of Australia they used ta combination of equal-area splines and a spatial prediction model to map the whole are for bulk density and organic carbon. In a typical soil profile, the first centimetres can have a very high organic matter content (Batjes, 1996), derived from plant and animal residues. With depth the percentage in volume of organic carbon decreases non-linearly with a pattern that is frequently modelled with exponential functions (Russell and Moore, 1968; Hilinski, 2001). There are also other equations proposed to describe the decrease of soil carbon with depth, but they are just a variance of the exponential model (Arrouays and Pelissier, 1994; Bernoux et al., 1998). Minasny et al. (2013) used a generalized negative exponential depth function. However, in very wet soils decay is inhibited, causing the accumulation of organic matter and the formation of peat, which can be found either continuously throughout the soil profile, or as a buried horizon. In either case, C in these soils exhibits a different behaviour by depth than the previously described typical soil. For this reason, considering only volume averages cannot be sufficient to accurately assess the potential changes driven by climate change. There are broadly three pathways through which C is lost from soil; through erosion, leaching (hydrological) and respiration (Dawson and Smith, 2007). Each of these is sensitive or depends on the vertical distribution of soil C. Soil C loss through water erosion can occur in two ways, through sheet flow over the surface of the soil or through channel erosion. Overland flow removes top soil and changes the distribution of soil C through the profile. The amount lost through erosion also directly depends on the vertical distribution of soil C. Hydrological loss of C from soil is a function of the movement of C from soil OM to the soil water, which is then lost through leaching as discharge from soil to the surface water. Dawson and Smith (2007) describe numerous studies in which
579
significant correlations between soil C and dissolved organic C in UK catchments are observed; movement and loss of C through this pathway depends on the spatial and vertical distribution of soil C (Grieve, 1990; Hinton et al., 1998). Peak C loss tends to coincide with peak discharge events, in which the 15–20 cm surface layer dominates flow characteristics (Worrall et al., 2003). This would suggest that the vertical and horizontal distribution of soil C determines both the concentration and also flow characteristics of C loss from soil. The largest amount of C is lost through organic matter mineralization. These losses are typically modelled using spatially distributed, integrated land–atmosphere process models such as JULES (Harrison et al., 2008). Improvements to the vertical distribution of soil C component in these types of models increases their response-sensitivity to changes in soil stocks and processes. Jones et al. (2005) compared the outputs from a generic Soil C model with measurements of topsoil (organic) C aggregated spatially according to both soil mapping and administration units. In a comprehensive review by Minasny et al. (2013), the outcome of various studies are tabulated and where, in general, R2 for prediction accuracy range from 0.21 to 0.98 for soil C prediction from areal imagery, the statistics obtained by cross-validation, internal or external validation. There is only one study in this list in which the C stock predictions are independently validated (R2 of 0.41). In this paper we consider a method developed previously (Veronesi et al., 2012) for in field mapping of soil compaction, and we consider this now for regional scale C stocks. The method is developed and applied within an area of the UK with a high density of soil profile observations, but validated using an independent topsoil dataset sampled on a 5 km grid. The method creates a fine vertical resolution 3D map of the soil C stock, describing its behaviour both vertically and spatially. 2. Materials and methods 2.1. Study area The study area is located in the West Midlands region, extending approximately 13,948 km2 (Fig. 1). Within the pilot area, there are 118
Fig. 1. Map of the area under study. The area extends for 13,948 km2 over the West Midlands region and on the Welsh border.
580
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
different soil associations according to the National Soil Map (NATMAP) of England and Wales with the Bromyard Association having the largest extend (8.51%). The area is covered largely by loamy soils, although the silt content gets larger towards the East and North whereas the clay content gets larger towards the South. According to the Centre for Ecology & Hydrology (CEH) Land Cover Map 2000 (Fuller et al., 2000) around one third of the area is covered by arable land, which acts as a C source (Schils et al., 2008); 50% is covered by natural, semi-natural or improved grassland and 11% is urban or suburban areas. The remaining land is classified as woodland, inland water or heath. 2.2. Dataset description 2.2.1. Training data The analytical data used for this study came from the Representative Profile data collected by the Soil Survey of England and Wales, now the National Soil Resources Institute (NSRI) to support national mapping programmes. Representative Profiles were sampled by horizons as identified in the field. The data were collected between 1969 and 1985 with more than 90% covering the period 1973 to 1983. A total of 435 profiles were available for the study areas. In addition to the temporal variation, the locations were selected to describe model profiles for individual soil series used in the national soil classification, consequently, the samples are not very well distributed. There are clusters of high density sampling and a large area in the West, which coincides with the Welsh border, with few profiles (see Fig. 1). This is very often characteristic for legacy data and highlights their limitation as they were not collected for quantitative soil mapping approaches. However, they represent the only available national dataset.
2.4.1. 1st stage At first the focus is on the selection of the depth functions that best describe the variation by depth of the two soil properties of interest (soil C and BULKD), basing the choice on a framework for testing different types of mathematical equations followed by identifying the most appropriate number of coefficients. The analysis was undertaken using DataFit from Oakland Engineering which allows to investigate a wide range mathematical equations. For the vertical distribution of soil C we considered the general category of exponential functions because they are known from the literature as being able to describe the variations by depth of soil C (Gifford and Roderick, 2003; Hilinski, 2001; Hiederer, 2009; Minasny et al., 2013). For BULKD we noticed that the general pattern was also exponential. From the analysis, two equations were selected: – A Power function for describing the vertical distribution of soil C: β
C i ¼ β 1 di 2
ð1Þ
where Ci is the value of C at the ith depth interval and di is the depth. – An inverse power function for describing the vertical distribution of soil BULKD:
BULKDi ¼
1 β
β1 di 2
ð2Þ
where BULKDi is the value of BULKD at the ith depth interval and di is the depth.
2.3. Environmental covariates In order to improve the accuracy of the soil C and BULKD predictions we implemented four grids of environmental covariates, namely: – Digital elevation model (DEM). A 50 m resolution DEM, created by the Ordnance Survey and constructed from the 1:50,000 topographic map. From this five derivatives were obtained: slope, aspect, total curvature, profile curvature and plan curvature. – Soil association map from NSRI. A 1 km grid map with information about the predominant soil association in each grid cell. – Land Cover Map 2000. Available from CEH, it has a resolution of 1 km with information about the predominant usage of the land. – Digital Geological Map (DigMap 250) from the British Geological Survey at the scale of 1:250,000.
2.3.1. Validation data The National Soil Inventory (NSI) (Loveland, 1990) was established to obtain an unbiased estimate of the distribution of the soils of England and Wales and of the chemistry of the topsoil (0–15 cm depth). It followed a non-purposive design on a 5-km grid covering all soil types, habitats and land uses. The original sampling (1978–83) was repeated on approx. 40% of sites 12–25 years later in such a way as to detect a specified change in the soil with specified confidence. Samples of air-dried soil from both samplings were archived. The NSI is unique among national soil inventories globally in having been comprehensively resampled. 2.4. Methodology The top-down mapping method (Veronesi et al., 2012) relies on soil specific depth functions to describe the variation of soil properties with depth. Its application on the soil C and BULKD dataset can be summarized in three stages:
2.4.2. 2nd stage In the second stage of the mapping process, we fitted Eqs. (1) and (2) to each available soil profile, computing the local R2 value. From the goodness of fit analysis we produced two R2 maps. 2.4.3. 3rd stage The final step of the mapping process involves applying a universal co-kriging algorithm to estimate the value of the two coefficients of the depth functions, on a regular grid of 1 km of horizontal resolution, assuming no variation within the grid cell. Universal kriging was adopted because at this scale the assumption of intrinsic stationarity cannot hold and the soil process needs to be considered as non-stationary. This interpolation method is able to work under these circumstances by modelling the trend in the variable using a multi-linear regression between the soil property and the environmental covariates and then creating the variogram from the residuals of the regression (Hengl et al., 2007). Subsequently the depth functions are solved for each interpolated point in order to create a lattice of 3D estimates with a 10 cm vertical resolution. 2.5. Validation The validation was undertaken using two methods (Minasny et al., 2013): cross-validation and an independent validation against the NSI data. The cross-validation approach relies on excluding a set percentage of soil profiles (5, 10, 20, 30%) and re-predicting their values. The observed and predicted values are then compared, calculating the R2. The second method is an independent validation accomplished by comparing the soil C percentage available in the NSI dataset, with the estimated soil C percentage. To assess the goodness of fit, the residuals between the observed values in the NSI dataset and the predicted values with the top-down mapping method were computed.
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
581
Fig. 2. General profiles of the two soil properties and depth functions fitting. Red crosses signify the mean value of the soil property in each depth layer.
Subsequently, the C stock was computed with Eq. (4) (Xu et al., 2011):
2.6. Soil C stock estimation The soil C density was calculated using the following equation, suggested by (Xu et al., 2011):
C Stock ¼ C Density ¼ C i BULKDi di
ð3Þ
where ‘C Density’ (t/ha) is the soil carbon density; Ci is the C (%) content of the soil layer i, BULKDi is the BULKD (g/cm3) value of the soil layer i and di is the depth of the soil layer i (cm). We then overlaid the resulting density estimations with the Land Cover Map 2000 from CEH (Fuller et al., 2000) and we made the following corrections, as suggested in Bradley et al. (2005) and Xu et al. (2011): – For urban areas and water bodies the density value was set to 0. – For semi urban areas the value was set to one half of the estimated value.
X C Density Area 1010
ð4Þ
where ‘C Stock’ (Tg) is the carbon stock value estimated after the corrections cited above and ‘Area’ is the area of a single prediction grid cell in m2. We calculated the C stock value for each depth interval estimated during the mapping process, namely: 0–10; 10–20; 20–30; 30–40; 40–50; 50–60; 60–70; 70–80; 80–90; and 90–100 cm. We also computed average values for topsoil (0–30 cm) and subsoil (30–100 cm). 2.7. C density variation estimation We estimated the uncertainty related to each point on the prediction grid adapting a method proposed by Goidts et al. (2009) and Schrumpf et al. (2011), for assessing the error propagation of each variable in the C
Fig. 3. R2 values for soil C and bulk density within the study area.
582
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
Table 1 Descriptive statistics of the observed and predicted datasets. Topsoil refers to the 0 to 30 cm soil layer from the surface; the subsoil refers to 0 to 100 cm soil layer. Observed dataset Min
Max
Mean
Median
St.Dev.
0.30 0.10 0.59 0.61
43.00 39.10 1.76 1.82
3.06 1.01 1.23 1.45
2.30 0.60 1.25 1.47
3.23 2.73 0.23 0.16
Predicted with top-down mapping SOC Topsoil 0.00 (%) Subsoil 0.00 BULKD Topsoil 0.27 (g cm3) Subsoil 0.74
45.75 68.08 2.82 2.49
4.36 0.45 1.07 1.35
2.95 0.24 1.14 1.35
4.24 0.89 0.26 0.11
SOC (%) BULKD (g cm3)
Topsoil Subsoil Topsoil Subsoil
Density estimation. The variance of the C Density assessment is based on the following equation:
2
VarðC DensityÞ ¼ ðC DensityÞ
! 2 2 σC σ BULKD σ C−BULKD þ þ 2 C BULKD C 2 BULKD2
ð5Þ
where Var(C Density) is the error associated with the density prediction, in t/ha; σC and σBULKD denote the standard deviation of carbon content and bulk density; and σC–BULKD is the covariance of C and BULKD. We applied this equation to two depth intervals, namely the 0–30 cm layer (topsoil) and the 0–100 cm layer (subsoil). 3. Results 3.1. Landscape level soil C maps
curves were fitted to the entire dataset, still yielded R2 values of 0.48 and 0.59 respectively. In order to numerically demonstrate the accuracy of these two equations in describing the spatial pattern of soil C and BULKD variation by depth, we cross-validated the fits to each individual profile, and obtained a profile R2 (Fig. 3). The average profile R2 for was 0.92 for soil C and 0.85 for BULKD, with the majority of profile R2 in the range of 0.90 to 0.99 for soil C and 0.80 to 0.99 for BULKD respectively. We classed profiles which were not adequately described by the depth functions (Eqs. (1) and (2)) as those that had R2 of less than 0.5 (i.e. the function, on the whole, did not describe the majority of the variation by depth). We observed no such R2 in soil C, which would suggest that in this area and at this scale the soil C profiles present a relative homogeneous pattern. For BULKD there are only two areas which present R2 values below 0.5, present in the more southern part of the study area. This is due to the fact that three profiles in that area have a linear depth pattern that the inverse power function struggles to describe. Normally, bulk density increases with depth but in these three cases it remains almost constant. These differences are present in only three samples. This is a key property of the depth function mapping method (Veronesi et al., 2012). The fact that the selected depth functions can describe only a small set of geometries means that, from the goodness of fit map, we can extract very useful information about the 3D spatial variability in the study area. In this case, we can identify an area of linear increase in BULKD values along the profile. This way this method can be used as diagnostic tool in order to detect changes in the soil profile that otherwise would be identifiable only by looking at profiles individually. We subsequently estimated the value of soil C and BULKD with a universal co-kriging algorithm. The final 3D grid of predictions has a horizontal resolution of 1 km and vertical resolution of 10 cm. The total number of points in the prediction grid was 220,928. The
The two depth functions for describing the soil profiles of C and BULKD, were selected because they can describe the majority of the profiles with an R2 above 0.6 having only two coefficients (see Eqs. (1) and (2)). These two aspects are equally important in the definition of the best function as a high number of coefficients can potentially result in a propagation of the error during the interpolation phase, decreasing the overall accuracy. The validity of the two depth functions in describing the general pattern of the two soil properties is visualised in Fig. 2, where the red crosses represent the average value at each depth layer and the function
Fig. 4. Distribution of the residuals between the soil values observed in the NSI dataset and the values predicted in this study.
Fig. 5. Spatial pattern of the residuals of the validation against the NSI data.
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
583
Fig. 6. Soil C stock maps.
computational time for the entire mapping process was around 5 min, which is a fraction of the time needed to estimate the same number of points using 3D geostatistics alone, which for this number of points can take much longer (Kerry and Hawick, 1998). To investigate the results of the mapping, the summary statistics of the original dataset were compared with the predicted datasets (Table 1). For the topsoil (0–30 cm) the predicted dataset slightly overestimates the values of both C and BULKD and underestimates the variance. The mean C value in the observed data for the topsoil is 3.06 with a standard deviation of 3.23, compared to a predicted mean of 4.36, with a standard deviation of 4.24. Regarding the BULKD values, the observed mean for the topsoil is 1.23, with a standard deviation of 0.23, while the predicted mean is 1.07, with a standard deviation of 0.26. On the other hand for the subsoil (0–100 cm), the results are less accurate. The C observed mean is 1.01, with a standard deviation of 2.73, while the estimated mean is 0.24, with a standard deviation of 0.89. For BULKD, the subsoil results are much better. The observed mean is 1.45, with a standard deviation of 0.16, while the estimated mean is 1.35, with a standard deviation of 0.11. These results indicate that the
Table 2 Carbon stock estimates divided by land use. Topsoil refers to the 0 to 30 cm soil layer from the surface; the subsoil refers to 0 to 100 cm soil layer.
Agricultural land Grassland (natural, semi-natural and improved) Woodland
Topsoil C stock (Tg)
Subsoil C stock (Tg)
20.68 ± 0.17 88.9 ± 0.38
54.42 ± 0.22 33.23 ± 0.25
8.02 ± 0.31
2.84 ± 0.16
top-down mapping method can reproduce the general distribution of the soil dataset very well, even at landscape scale. The topsoil estimates appear to be more accurate than subsoil predictions. This means that the two depth functions are very accurate in describing the topsoil spatial variability and less accurate in lower horizons. This is probably due to a decrease in the number of profiles with data at lower depths, thereby decreasing the sample support at the lower depths. 3.2. Validation The accuracy of the mapping method was then tested with two methods: cross-validation and independent validation. The cross-validation results indicate an average value of R2 for the lowest exclusion percentage of 36% for C, with a standard deviation of 0.20; and 44% for BULKD, with a standard deviation of 0.15. These results may appear relatively low, but they are in line with previous validated research of 3D mapping at the catchment or landscape scale (Malone et al., 2009; Kempen et al., 2011; Minasny et al., 2013). The problem with mapping large areas using legacy soil observations is that the dataset was not sampled with probabilistic designs, making it less suitable for assessing the accuracy of maps. In this study, the majority of the area around Birmingham is well covered by clusters of high-density data, while the remaining area was only sparsely sampled. It is therefore difficult to consider cross-validation a sound method for validating the map, as when an excluded sample is re-predicted using other observations tens of kilometres apart, it is clear the estimation may substantially differ from the observation. For this reason we also performed an independent validation, comparing the topsoil C prediction against the data in the NSI dataset. To
584
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
assess the accuracy of the prediction, we calculated the residuals between the observed and the predicted values. Fig. 4 contains the distribution of the residuals, while in Fig. 5 they are mapped. The results indicate that the residuals are generally very low, with a mean C of 0.59% and a standard deviation of 4.46%. In Fig. 4, it is possible to appreciate that the majority of the residuals are in the interval between −5% and + 5%, even though there is a small number of large percentage values, which in Fig. 5, these higher values are in areas poorly sampled in the representative profiles data. A plausible explanation for this is that there are areas of probable marshes where the topsoil C content is particularly high, but these are almost impossible to predict if those areas were not originally sampled.
relatively low level of uncertainty for the soil C predictions. In particular the majority of the topsoil predictions have an associated variance between 0 and 0.19 t/ha with peaks of 0.95 t/ha in the Welsh border area, where the sample density is very low. Moreover, the subsoil uncertainty estimations are also low, with the majority of the point uncertainties between 0 and 0.87 t/ha, and a peak on the Welsh border of 5.07 t/ha. These results corroborate the data presented in the previous sections. In the topsoil, the top-down mapping method can achieve a generally higher level of accuracy. On the other hand, from the data described in the previous sections it is impossible to give a precise indication of the accuracy of the method in the subsoil. However, this uncertainty estimation gives us a clear indication of the overall accuracy of the mapping process, which can achieve good results also in the subsoil.
3.3. Stock estimation 4. Discussion The soil C stock was calculated and corrected using the procedure described above. The results indicate that the C stock in the topsoil is 162.91 ± 3.63 × 10−5 Tg, and the C stocked in the subsoil is 66.71 ± 3.52 × 10−5 Tg for the study area respectively. The C stock estimates are presented in the maps in Fig. 6. It appears that on the higher ground on the Welsh border, the topsoil has a relatively higher C content, due the fact that the majority of this area is grassland. By contrast, the subsoil map presents a picture where the higher C content is in the arable lands around Birmingham, and this can probably be a result of leaching of C from the arable topsoil as suggested by Bellamy et al. (2005). The C stock was also calculated for land use categories and the results are shown in Table 2. The uncertainties of the C Density point estimations were computed with the approach described in the Materials and methods section. The results are presented in Fig. 7. From this figure, we can appreciate the
In general, we can divide the study area into two distinct zones: the lowlands around Birmingham, where most of the arable and urban areas are concentrated, and the higher ground in the northern and eastern parts. These two zones are characterized by almost opposite soil profile patterns, where the lowlands exhibit low C stocks in the topsoil and relatively high stocks in the subsoil, and the higher ground presents relatively high C stocks in the topsoil with a general decrease with depth. The cause for this pattern, in the lowland areas, is the land use and land management, as the predominant use is arable agriculture. Increased soil organic matter decomposition rates and top soil erosion depleted the topsoil of OM, instead leaving the subsoil untouched or enriched with C from leaching (Dawson and Smith, 2007). On the higher ground, the density of arable and urban land is much lower and this has left these areas in an almost natural state, with the
Fig. 7. Soil C stocks uncertainty maps.
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
585
Fig. 8. Soil C percentage by volume on agricultural land.
higher concentration of C in the topsoil. This is a possible risk factor particularly for the Welsh border, because the soils there have a high silt percentage and this makes them vulnerable to erosion, as mentioned in the report of the PESERA pan-European erosion model (Kirkby et al., 2004). The fact that in these areas the value of C stock is higher indicates a high risk of C losses by erosion and it would seem to suggest that erosion protection land management practices would be appropriate. Landscape level estimates of soil C stocks are usually based on soil surveys which were not designed for assess C stocks, changes and drivers (Schrumpf et al., 2011). Samples are taken at locations in the landscape associated to representative soil characteristics and at depths representing soil horizons. As others (Malone et al., 2009) the method presented in this paper is capable of taking horizon based depths, transforming this to a continuous description of the depth variation and interpolating this across the landscape. For this landscape, we compare these estimates to the currently used map of the C stock for the United Kingdom created by Bradley et al. (2005) and our results are higher for the study area. Bradley et al. (2005) calculated a topsoil stock of 91.33 Tg and subsoil stock of 48.66 Tg compared to our estimates of C stocks of 162.91 ± 3.63 × 10−5 Tg in the topsoil, and 66.71 ± 3.52 × 10 −5 Tg in the subsoil. Bradley et al. (2005) used an approach based on pedo-transfer functions, and different soil datasets: the NSI data for the topsoil, extrapolating the soil C value from 15 to 30 cm depth because the NSI data are averages of the first 15 cm, and the soil series map to extrapolate the C values for the subsoil. The method in essence produces aggregate means for soil polygons, which seems to underestimate soil C stocks.
If we consider this map as input to soil C transformation and transport models, i.e. respiratory responses such as JULES, erosion or hydrological losses, then these results would suggest that general behaviour of C and bulk density by depth is fairly consistent across this type of landscape (dominated by arable land use) and approaches such as that described by Ungaro et al. (2010), in which a mean stock is obtained by a general depth function and subsequently interpolated are adequate to describe C stocks and C concentrations at a landscape scale. We could also consider generating slices if the C transformation models need a more specified vertical distribution; and this method, as well as equal area splines as described by Malone et al. (2009), would be adequate. The advantage of the method described in this paper is that it contains an element of diagnosis; where the dominant function does not describe the depth variation adequately, stratification of the landscape can be considered, with different depth functions for different portions of the landscape (e.g. arable vs peat soils). The mapped soil C can also be considered as degradation maps, where at less than 2%, soil C volume is so low that there is a significant decline in soil structural stability and fertility (Loveland and Webb, 2003). From the mapping results, we extrapolated the predicted points on the land defined as in arable use and then we computed the percentage of C stock by volume; the resulting map is presented in Fig. 8. Fig. 8 shows a scenario where the majority of the arable land has values of C below or slightly above the 2% threshold, meaning they are currently at very high risk of loss of fertility and soil structure, i.e. a risk map of loss of soil functions due to extremely low levels of soil C and therefore do not function effectively in the terrestrial C cycle. Conversely, these are areas in which, through changes in land management practices (such as leaving straw
586
F. Veronesi et al. / Science of the Total Environment 487 (2014) 578–586
and stubble on the fields), significant gains can be achieved increasing soil C stocks in the future. 5. Conclusions This method of 3D mapping with depth functions, can achieve a high level of accuracy, in particular for estimating the first 30 cm of soil. This indicates that this map can be readily implemented in soil erosion model, for example, incr'easing the amount of information in the very top part of the profile. As to the subsoil, the uncertainty estimation indicates a good level of accuracy, with the majority of the soil C stock values presenting a variance between 0 and 0.87 t/ha. On the one hand, by analyzing the soil patterns and comparing them with the land use and the soil erosion model, it is possible to conclude that this area is at high risk of soil C losses. In the higher ground on the Welsh border the soils, which have a high concentration of C, are prone to erosion because of their high silt content. On the other hand, from the risk assessment it is possible to conclude that the majority of the arable land around Birmingham is in danger of losing soil functionality such as fertility, as their percentage of C is below the 2% by volume threshold. References Arrouays D, Pelissier P. Modeling carbon storage profiles in temperate forest humic loamy soils of France. Soil Sci 1994;157:185–92. Bernoux M, Arrouays D, Volkoff B, Cerri CC, Jolivet C. Bulk densities of Brazilian Amazon soils related to other soil properties. Soil Sci Soc Am J 1998;62:743–9. Batjes NH. Total carbon and nitrogen in the soils of the world. Eur J Soil Sci 1996;47: 151–63. Bellamy PH, Loveland PJ, Bradley I, Lark RM, Kirk G. Carbon losses from all soils across England and Wales 1978–2003. Nature 2005;437:245–8. Bradley RI, Milne R, Bell J, Lilly A, Jordan C, Higgins A. A soil carbon and land use database for the United Kingdom. Soil Use Manag 2005;21:363–9. Dawson JJC, Smith P. Carbon losses from soil and its consequences for land-use management. Sci Total Environ 2007;382:165–90. Ellert BH, Bettany JR. Calculation of organic matter nutrients stored in soils under contrasting management regimes. Can J Soil Sci 1995;75:529–38. Fuller RM, Smith GM, Sanderson JM, Hill RA, Thomson AG, Cox R, Brown NJ, Clarke RT, Rothery P, Gerard FF. Land Cover Map 2000: a guide to the classification system, Centre for Ecology and Hydrology, 2002. http://www.google.co.uk/url?sa=t&rct= j&q=&esrc=s&frm=1&source=web&cd=1&ved=0CC8QFjAA&url=http%3A%2F% 2Fwww.ceh.ac.uk%2Fdocuments%2Flcm2000_classification_guide.pdf&ei=wZ2UUqy_ Foi9ygOsr4GwAQ&usg=AFQjCNFtcxS7T5TDTkD7YuTma0dKrH5GLA&bvm=bv. 57155469,d.bGQ, 2000. Gifford RM, Roderick ML. Soil carbon stocks and bulk density: spatial or cumulative mass coordinates as a basis of expression. Glob Chang Biol 2003;9:1507–14. Goidts E, Van Wesemael B, Crucifix M. Magnitude and sources of uncertainties in soil organic carbon (SOC) stock assessments at various scales. Eur J Soil Sci 2009;60:723–39. Grieve IC. Seasonal, hydrological and land management factors controlling dissolved organic carbon concentrations in the Loch Fleet catchments, southwest Scotland. Hydrol Process 1990;4:231–9. Grimm R, Behrens T, Märker M, Elsenbeer H. Soil organic carbon concentrations and stocks on Barro Colorado Island — digital soil mapping using Random Forests analysis. Geoderma 2008;146:102–13.
Harrison RG, Jones CD, Hughes JK. Competing roles of rising CO2 and climate change in the contemporary European carbon balance. Biogeosciences 2008;5:1–10. Hengl T, Heuvelink G, Rossiter DG. About regression-kriging: from equations to case studies. Comp Geosci 2007;33:1301–15. Hiederer R. Distribution of organic carbon in soil profile data. EUR 23980 EN. Luxembourg: Office for Official Publications of the European Communities; 2009. Hilinski TE. Implementation of exponential depth distribution of organic carbon in the CENTURY model. Dept Soil Crop Sci, Colorado State Uni; 2001 [http://www.nrel. colostate.edu/projects/century5/reference/html/Century/exp-c-distrib.htm]. Hinton MJ, Schiff SL, English MC. Sources and flowpaths of dissolved organic carbon during storms in two forested watersheds of the Precambrian Shield. Biogeochemistry 1998;41:175–97. Intergovernmental Panel on Climate Change — IPCCClimate change: synthesis report, 2007; 2007. Jobbagy E, Jackson R. The vertical distribution of soil organic carbon and its relation to climate and vegetation. Ecol Appl 2000;10:423–36. Jones RJA, Hiederer R, Rusco E, Montanarella L. Estimating organic carbon in the soils of Europe for policy support. Eur J Soil Sci 2005;56:655–71. Kempen B, Brus DJ, Stoorvogel JJ. Three-dimensional mapping of soil organic matter content using soil type — specific depth functions. Geoderma 2011;162:107–23. Kerry KA, Hawick KE. Kriging interpolation on high-performance computers. In: Sloot P, Bubak M, Hertzberger B, editors. High-performance computing and networks, International Conference and Exhibition Amsterdam, The Netherlands, April 21–23; 1998. p. 429–38. Kirkby MJ, Jones RJA, Irvine B, Gobin A, Govers G, Cerdain O, Van Rompaey AJJ, Le Bissonnais Y, Daroussin J, King D, Montanarella L, Grimm M, Vieillefont V, Puigdefabregas J, Boer M, Kosmas C, Yassoglou N, Tsara M, Mantel S, Van Lynden GJ, Huting J. Pan-European soil erosion risk assessment: the PESERA map. Special Publication Ispra, vol. 73; 2004 [no. S.P.I.04.73]. Loveland P, Webb J. Is there a critical level of organic matter in the agricultural soils of temperate regions: a review. Soil Tillage Res 2003;70:1–18. Loveland PJ. The National Soil Inventory of England and Wales. In: Lieth H, Markert B, editors. Element concentration cadasters in ecosystems. Weinheim, Germany: VCH; 1990. p. 73–80. Malone BP, McBratney AB, Minasny B, Laslett GM. Mapping continuous depth functions of soil carbon storage and available water capacity. Geoderma 2009;154:138–52. Minasny B, McBratney AB, Malone BP, Wheeler I. Digital mapping of soil carbon. Adv Agron 2013;118:1–47. Russell JS, Moore AW. Comparison of different depth weightings in the numerical analysis of anisotropic soil profile data. Transactions of the 9th International Congress of Soil Science 1968;4:205–13. Schils R, Kuikman P, Liski J, Van Oijen M, Smith P, Webb J, Alm J, Somogyi Z, Van den Akker J, Billett M, Emmett B, Evans C, Lindner M, Palosuo T, Bellamy P, Jandl R, Hiederer A. Review of existing information on the interrelations between soil and climate change. CLIMSOIL final report; 2008. Schlesinger WH, Bernhardt ES. Biogeochemistry, an analysis of global change. San Diego, California, USA: Academic Press; 1997. Schrumpf M, Schulze ED, Kaiser K, Schumacher J. How accurately can soil organic carbon stocks and stock changes be quantified by soil inventories? Biogeosciences 2011;8: 1193–212. Ungaro F, Staffilani F, Tarocco P. Assessing and mapping topsoil organic carbon stock at regional scale: a simple kriging approach conditional on soil map delineation and land use. Land Degrad Dev 2010;21:565–81. Veronesi F, Corstanje R, Mayr T. Mapping soil compaction in 3D with depth functions. Soil Tillage Res 2012;124:111–8. Worrall F, Burt T, Adamson J. Controls on the chemistry of runoff from an upland peat catchment. Hydrol Process 2003;17:2063–83. Xu X, Liu W, Zhang C, Kiely G. Estimation of soil organic carbon stock and its spatial distribution in the Republic of Ireland. Soil Use Manag 2011;27:156–62. Zan CS, Fyles JW, Girouard P, Samson RA. Carbon sequestration in perennial bioenergy, annual corn and uncultivated systems in southern Quebec. Agric Ecosyst Environ 2001;86:135–44.