Whole Soil Bulk Density

Overview

An attempt was made to update digital soil mapping of whole soil bulk density for Australia. This was an update of first attempt by Viscarra Rossel et al. (2014). Based on model evaluations using a dataset not included in any modelling, the updated version (2nd Version) represents a demonstrable improvement on the 1st version.

Since the first version, more measured site data has been made available and retrievable via the Australian SoilDataFederator. In 2014 there were 3776 sites with measured whole soil bulk density. For the new update, 6116 sites had measured data. Because of usually strong empirical relationships between bulk density, soil texture and soil carbon, the use of pedotransfer functions (to predict bulk density from soil texture and soil carbon) was performed with the intention of increasing data density and spatial coverage of data that would ultimately improve digital soil mapping prediction skill. This added a further 15735 sites after building a spatial pedotransfer function using a dataset of 12308 cases (3939 sites with bulk density, soil carbon and soil texture data).

The basic steps of the work entailed.

Use soil data federator to get pertinent soils observation data
Develop spatial pedotransfer function prediction whole soil bulk density using soil carbon and texture data.
Compile measured and inferred whole soil bulk density data (86306 cases), then setting aside a dataset of 7500 cases for external model evaluation.
Predictive models using random forest algorithm with 78806 data cases fitted. To account for uncertainties in pedotransfer function inferred data, Monte Carlo simulations were performed from the pedotransfer function model. Simulation was repeated 100 times.
Predictive model uncertainties quantified using UNEEC approach (Uncertainty Estimation based on local errors and Clustering).
Quantification of model extension limits derived using hybrid method involving multivariate convex hull analysis and count of observations.
Digital soil maps with quantified uncertainties (5th and 95th prediction interval limits) and assessment of model extrapolation risk were produced at 90m resolution for the following depths: 0-5cm, 5-15cm, 15-30cm, 30-60cm, 60-100cm, 100-200cm.

The Data

All soil data (whole soil bulk density and soil carbon and texture) were extracted from the soil data federator. The data was subjected through serval checks and processes to ensure data integrity and any spurious cases detected were removed from any further analysis. Raw data were then further processed or harmonized in terms of consistent depth support by fitting mass preserving splines then returning data for the following depth intervals: 0-5cm, 5-15cm, 15-30cm, 30-60cm, 60-100cm, 100-200cm.

A suite of predictive covariates were compiled for the use in national digital soil mapping and related tasks (Searle et al. 2022). For the modeling tasks in the present work, rather than sourcing from more than 100 thematic maps, principal component layers derived for each of the key soil forming factors (climate, relief, parent material, and organism) were used. Using the PCA layers, balanced up the number of contributing layers from each soil forming factor.

Figure 1. Distribution of site data with some measure of whole soil bulk density. Measures were wither performed in field or lab, or inferred from pedotransfer function.

Deriving pedotransfer function of whole soil bulk density

Using the available 12308 cases for model building, the Cubist machine learning algorithm learning relationships between the target variable (whole soil bulk density) and covariates which included measures soil carbon and clay. Models were also fitted with (spatial context) and without (aspatial context) environmental covariate information. Models fitted with spatial contextual information would be thought to capture local variations in data, making the pedotransfer function more extensible.

A key aspect of deriving a pedotransfer function is the quantification of uncertainties. These were done using bootstrapping which involves repeated data re-sampling (with replacement) and model fitting. 100 repeated model systems were derived, and the combined model systems are evaluated based on the out-of-bag data that are not used in each model system. Both the structural and systematic errors of the bootstrap model were used to derive for each data case a probability density function which is summarized in terms of a mean and variance. Goodness of fit model e valuations based on out-of-bag cases are shown in the following table and plots (Figure 2).

Figure 2. Goodness of fit model evaluations and modelled vs. observed plots for both pedotransfer function types to predict whole soil bulk density using soil carbon and clay content and with or without environmental covariates.

The models were extended to all cases where only soil carbon and texture data were available. For each data case, the mean and variance were used as inputs for a random number generator to simulate individual outcomes [1000]. 100 outcomes were saved to serve as data inputs for spatial modelling.

Predictive spatial soil modelling and uncertainty quantification

Modelling was underpinned by random forest model. While random forest model already has built into its algorithm the bootstrap resampling of available data, it is not completely appropriate in cases where input data have measurable uncertainty, although some weighting mechanisms may be explored if deemed appropriate. Instead, to account for data ‘measurement errors’ due to much of the available data being inferred via pedotransfer function, a further bootstrap routine was implemented on top of the built-in random forest one. This was to ensure appropriate handling of the prediction variations brought about by the pedotransfer model simulations. 100 bootstrap modelled were fitted. After each model fitting, the model was extended to the out-of-bag data where a model residual was calculated and saved. After the 100 realisations, and average model residual (derived from out-of-bag estimates) was calculated. The model residuals were to be used as input into the modelling uncertainty quantification.

Rather than fitting models for each depth interval, depth layer was used as a predictive covariate, to streamline model building and to allow soil vertical co-variation to be indirectly included in the modelling process. In total, of the 78806 cases available for model building, the following proportions comes from each of the depth intervals: 0-5cm (27%), 5-15cm (27%), 15-30cm (19%), 30-60cm (12%), 60-100cm (10%), 100-200cm (5%). The distribution of data over the different depths highlights the disparity of observation from the soil surface and near surface relative to deeper layers. Model evaluation was performed for each of the depth intervals.

The saved model residuals for each of the 78806 data cases were combined with the covariate data, then clustering was performed using fuzzy k-means clustering. To determine the optimal number of clusters for each depth interval, three criteria were calculated for each cluster number configuration in size steps from 5 to 1000 (with a step size of 2) with a fuzzy exponent value of 1.1. These criteria were:

The fidelity of the 1:1 relationship between confidence interval and coverage probability.
The prediction interval range at the 90th confidence interval.
The fidelity of the 1:1 relationship between the 90th and 99th confidence levels and coverage probabilities.

The optimal class number size for each depth was: 0-5cm (103), 5-15cm (103), 15-30cm (65), 30-60cm (679), 60-100cm (721), 100-200cm (749). With the optimal class numbers derived, the model error distributions for each class were calculated, then with particular interest in recording the 5th and 95th percentiles in order to calculate the associated prediction interval maps using the following equations:

Lower prediction interval limit.

Upper prediction interval limit.

where and correspond to the weighted lower and upper limits for the ith observation. and are the lower and upper limits for each cluster j, and is the membership grade of the ith observation to cluster j.

Figure 3. Table at top of image collection shows the model evaluations from external data for both the updated (SLGA Version 2) and original (SLGA Version 1) whole soil bulk density models for each of the specified depths. The middle xy-plots of predicted vs. observed is shown for the 0-5cm depth interval. A similar visual pattern is observed for the other depths and thus not shown in this document. Similarly, the prediction interval coverage probability plot is for the 0-5cm depth interval showing a well constructed quantification of model uncertainties and also similar for the other depth intervals.

Figure 4. National digital soil mapping of whole soil bulk density (0-5cm) show both the mean prediction and upper and lower prediction limits. The maps on the left column are for SLGA Version 1, while to ones on the right are for SLGA Version 2. These maps show similar general spatial pattern in the mean prediction at the national scale. The range of the uncertainties do appear smaller in general for SLGA Version 2, improving the use and reliability of the mapping.

Quantification of model extension limits

Understanding the extent of model extrapolation was done using a combination of geometrical and distance-based calculation using the available model data and associated environmental covariates. Geometrical method was based on multidimensional convex hull analysis as described on van den Hoogen et al. (2021). This is a relatively simple check to determine whether new cases based on their given covariate data value is in or out of the multidimensional space. If outside the space, the assumption is that model extrapolation to those cases will be unreliable. The distance-based calculation is a follow up metric assessing the number of observed data that are similar in terms of covariate data similarity. Once some distance threshold is met the metric is simply the number of observed data within the defined threshold that are similar to a given case which in mapping terms is a pixel (Malone et al. 2019). Interpretation of the mapping of model extension limits should look less at the number but rather the relative values across the mapping domain as this serves to highlight general regions where model reliability is ok or not.

The 2 maps below provide an indication of model extrapolation risk for the 0-5cm and 100-200cm depth intervals for a contrasting effect. Clearly for the 0-5cm depth interval, a large proportion of the continent and in key agricultural areas, existing observational data coverage appears to provide useful information when models are extrapolated to areas where data are limiting. Exceptions exist in northern Australia, western Tasmania, and the complex terrain along the eastern and western margins of the continent. Some areas in central Australia are not well captured by the existing set of available data. By contrast, for the 100-200cm depth interval, models are largely limited in their extrapolation, due largely to the sparsity of data at this depth. The disparity of observation between top and deeper than 30cm or 40cm is common in most soil databases (especially soil bulk density), so depictions shown in the maps below would be expected. Methodological, the work of quantification of model extrapolation risk is still incomplete and needs further attention as currently the method does not take into consideration correlation of soil properties with depth i.e., one could confidently predict whole soil bulk density at depth given some measurements made higher up in a soil profile. Because of this, assessments of model extrapolation might be of a more positive outlook than currently depicted. However, the fact there is massive disparity between top and subsoils in terms of observational data is a fundamental issue that still needs to be addressed irrespective of improvements to the method framework for quantifying model extrapolation risk.

Figure 5. Assessment of model data coverage, and by association, model extrapolation risk for whole soil bulk density for the 0-5cm and 100-200cm depth intervals. With increasing soil depth, the risk of extrapolation errors increases due to an increasing sparsity of observed data.

References

Malone, B.P., Minansy, B., Brungard, C., 2019. Some methods to improve the utility of conditioned Latin hypercube sampling. PeerJ 7, e6451.

Searle, Ross; Malone, Brendan; Wilford, John; Austin, Jenet; Ware, Chris; Webb, Mathew; Roman Dobarco, Mercedes; Van Niel, Tom (2022): TERN Digital Soil Mapping Raster Covariate Stacks. v2. CSIRO. Data Collection. https://doi.org/10.25919/jr32-yq58

Viscarra Rossel, Raphael; Chen, Charlie; Grundy, Mike; Searle, Ross; Clifford, David; Odgers, Nathan; Holmes, Karen; Griffin, Ted; Liddicoat, Craig; Kidd, Darren (2014): Soil and Landscape Grid National Soil Attribute Maps - Bulk Density - Whole Earth (3" resolution) - Release 1. v5. CSIRO. Data Collection. https://doi.org/10.4225/08/546EE212B0048

van den Hoogen, J., Robmann, N., Routh, D., Lauber, T., van Tiel, N., Danylo, O., Crowther, T.W., 2021. A geospatial mapping pipeline for ecologists. bioRxiv, 2021.2007.2007.451145.