Estimating surrogates, utility graphs and indicator sets for soil capacity and security assessments using legacy data

Soil assessment frameworks

The SSAF, a newly developed comprehensive framework, offers a cohesive approach to understanding soils, guiding surrogate selections for soil functions, services, and threats (Evangelista et al. 2023b). It encompasses five key dimensions: capacity, condition, capital, connectivity, and codification, enabling a holistic perspective on soil analysis. Following this concept allows us to acknowledge gaps in our understanding about what may be impacting a site’s soil security. Like the SMAF and CASH, the SSAF follows three steps: (1) selecting indicators for a management goal from the existing dataset; (2) interpreting indicators by transforming the selected indicators into a unitless value or score; and (3) where appropriate, integrating the scores into a single value. The SSAF and SMAF do not prescribe indicators to use or the type of integration, recognising that there are an increasing number of algorithms that may be used (Evangelista et al. 2023b). Although the SSAF is not prescriptive on the indicators to use, its all-inclusivity provides guidance as to what indicators to seek and how to connect different aspects of soil security together.

The fertility capability classification (FCC), created in 1975 to characterise plant growth limitations in the tropics due to its focus on the function of soil as a producer of food and biomass, now describes soil characteristics relevant to both agronomic and ecologic disciplines across the world (Sanchez 2019). This framework is different from the others discussed by its focus on quantifying soil attributes of the top 50 cm considered to be inherent, or not easily changed within a span of years to decades (Sanchez et al. 2003), and being highly relevant to understanding soil security, such as soils with high phosphorous sorption capabilities, high leaching potential, and sodic soils. A set of 31 binary indicators, presence or absence, describing 20 different soil attributes along with a textural description of topsoil and subsoil have been developed over 40 years from expert knowledge. Within the SSAF framework, these inherent indicators would be classified under the capacity dimension and distinguished by their connection to a soil function of food and biomass production.

In Europe, as part of the ENVASSO project, 27 soil indicators were identified from an initial selection of 290 to protect soils in the EU (Huber et al. 2008). Potential indicators (PI) were obtained from literature review; indicators were then linked to nine soil threats and key issues within the respective threat, and final indicators chosen by expert judgment. This framework is slightly different from the other frameworks mentioned due to its distinctiveness in choosing a more specific set of indicators comprising three indicators per threat assessed, called TOP3. These indicators were carefully selected based on expert insights, considering their alignment with EU policy objectives. Similarly, the SSAF also considers policy aspect within its codification dimension which can be evaluated across various roles.

The SMAF, introduced by Andrews et al. (2004), is a site-specific tool designed to assess and evaluate soil management practices, particularly in relation to soil quality (capacity). As described above, the SMAF introduced the three steps used in the SSAF and CASH. This flexible framework is not prescriptive; however, the authors have identified over 70 PI and developed scoring functions for over a dozen along with a protocol to create the scoring functions (Wienhold et al. 2009). The CASH framework was derived based on the SMAF (Moebius-Clune et al. 2016) and focused on overall soil health (condition) within agricultural systems through the combined use of a predefined set of 18–25 indicators (from basic to NRCS-216 level of soil health analysis) covering the physical, biological, and chemical soil properties. It is then transformed into a score following the method described in Andrews et al. (2004), with value ranging from 0 to 100, with a higher number indicating optimum or near optimum condition. The scoring functions are texture based but undergoing efforts to account for texture, suborder classes, and mean annual temperature and precipitation in the continental USA using the Soil Health Assessment Protocol and Evaluation (SHAPE; Nunes et al. 2021). The overall soil health index (SHI) is determined by taking the average of all indicators measured. Within the SSAF, these indicators which are sensitive to anthropogenic activities would be classified under the condition dimension. In the SSAF, condition indicators of a target are compared to an area with similar pedogenic origin but minimal anthropogenic impact. Pedogenon can be defined as a conceptual soil taxon created from a regionalised set of quantitative state variables representing the soil-forming factors for a given reference time (Román Dobarco et al. 2021).

The SQAPP was formed to account for the impact that agricultural land management has had on soil properties and functions to be accessed via a mobile application. Six soil quality indicators were ultimately chosen based on work by Spiegel et al. (2015) and Bünemann et al. (2018): soil organic matter (SOM) content, pH, aggregate stability, water-holding capacity, and number of earthworms. The indicators are meant to show gradual soil quality and fertility changes in a span of more than 5 years, be soil and site specific, related to potential changes in soil functions and threats, and be easily interpretable by farm and land managers (Bai et al. 2018). Indicators in the SQAPP are presented as response ratios (RR) where the treatments were paired; for example, crop rotation versus monoculture. Rather than providing an integrating score, the SQAPP presents the results in radar charts, consistent with the idea that each indicator needs to be looked at individually. With the use of an application, this framework is focused on improving the connectivity different land managers have with their soil; connectivity is another dimension recognised under the SSAF and under which SQAPP data can be used to inform researchers and legislators.

The uniqueness of the SSAF lies in its ability to encompass multifaceted aspects of soil security. Unlike its counterpart’s framework that focuses solely on a certain dimension, the SSAF adopts a holistic approach. In the next section, methods implemented in the SSAF framework will be described.

The aims of this paper are to:

Study area and dataset

The study area is located in the Hunter Wine Country Private Irrigation District (HWCPID), New South Wales, Australia, covering an area of approximately 220 km² (Fig. 1). The climate in this region is temperate, with annual rainfall of 750 mm and seasonal temperature ranging from 4 to 30°C. The soil samples were collected in October 2019 across from the native forests belonging to three dominant soil orders (Chromosols, Kurosols, and Calcarosols) according to the Australian Soil Classification system (ASC). Following the soil map developed by Huang et al. (2018), the soil types within the study area can be mainly categorised into: Brown Chromosol, Brown and Red Kurosol, and Red Calcarosol.

Fig. 1.

Location of the sampling points for case studies within the Hunter Wine Country Private Irrigation District (HWCPID) shown in the red border on a May 2023 Landsat-9 image courtesy of the U.S. Geological Survey. The inset shows the location of the HWCPID within Australia. The Australian map shows territory and state divisions.

The dataset used corresponds to a subset of the data collected for a previous study by Xue et al. (2023). Eighteen sampling sites were selected, in which three soil cores (10 cm diameter) were collected as replicates. Within each soil core, the soil samples were collected across three depths (5–15 cm, 45–55 cm, and 90–100 cm) from the soil surface to 1 m depth.

A series of soil analyses were performed. For each soil sample, half was air-dried and ground to <2 mm prior to analysis through a commercial laboratory for 16 physico-chemical properties following methods described in Rayment and Lyons (2011). Analysis measured include (with method codes from Rayment and Lyons 2011): pH (1:5 CaCl₂) – 4B2; organic carbon (OC; Walkley and Black 1934) – 6A1; total nitrogen (N) (combustion) – 7A5; electrical conductivity (EC; 1:5 water) – 3A1; exchangeable aluminium (exch. Al; 1 M KCl) – 15G1; exchangeable sodium (exch. Na; 1 M NH₄CH₃CO₂) – 15D3; cation exchange capacity (CEC) – 15J1; and total (acid digest) phosphorus (P), copper (Cu), calcium (Ca), magnesium (Mg), manganese (Mn), nickel (Ni), potassium (K), sulfur (S), and zinc (Zn) – 17B1. Any elemental values that fall below the detection limit were recorded as zero (absent). The other half of the soil sample was subsampled and stored in the freezer (~ −20°C) before soil DNA extraction. Soil DNA was extracted from 10 g soil following the protocol provided in the supplementary (Supplementary A) in Metagen Lab®. DNA concentration was quantified using the Quantifluor dsDNA system (Promega, WI, USA). Primer sets of ArBa515F/Arch806R were used for polymerase chain reaction amplifications. Amplicon pyrosequencing was then sequenced on the Illumina MiSeq platform (PE 300). The sequencing data was filtered and merged using USEARCH (v.7.0) (Edgar 2010) and the merged sequences were clustered into operational taxonomic units (OTUs) by UPARSE (Edgar 2013). Soil microbial alpha diversity was calculated by Shannon index using the rarefied OTU table with vegan package (v2.6.2; Dixon 2003). For detailed information about the methods, we refer the reader to the original study (Xue et al. 2023). Two examples from this dataset are presented to exemplify the utilisation of the SSAF.

SSAF procedure

To evaluate the soil security assessment, a combination of a single or multiple roles (function, service, or threat) along with a single or multiple dimensions per role (capacity, condition, capital, connectivity, and codification) needs to be identified, as outlined in the SSAF (Evangelista et al. 2023b).

Once the combination of role(s) and dimension(s) have been identified, the SSAF procedure shown in Fig. 2 can be undertaken.

The subsequent sections will elaborate on each step of the SSAF procedure.

Fig. 2.

Conceptual framework for proposed soil security assessment framework (SSAF) method to create minimum dataset using (IVa) ordination and regression (IVb) method. U* is the utility of the transformed surrogate, PI is potential indicators, U SPI* is the utility of F-test screened potential indicators (SPI), and U ^* is the estimated utility. The obelisk (†) shows that the choice of regression will impact the outputs. In the case of neural networks, U ^* is the main output.

Potential indicator utilisation and weighting

Two of the most commonly used methods of weighting PI were explored: ordination and regression.

Ordination by principal component analysis

Principal component analysis is used as the ordination technique to assign weights to the PI, following the methods described in Andrews and Carroll (2001), Andrews et al. (2002a), and Andrews et al. (2002b) with slight modifications (Fig. 2 – section IVa). First, the transformed USPI* was used instead of using the SPI directly, and also normalised the final weights to add up to one. Details of this method are described in Fig. 3.

Fig. 3.

Conceptual framework using PCA for dimension reduction with steps proposed in this paper indicated by an obelisk (†). PC is the principal component, USPI* is the utility of F-test screened potential indicators, W_n = % variation explained by nth selected PC/∑ % variation explained by all selected PCs.

In short, the principal components (PC) with large eigenvalues (≥1) were selected (Andrews and Carroll 2001; Andrews et al. 2002a, 2002b). However, if less than three PC had eigenvalues of ≥1, another PC that explained ≥5% variations within the dataset was also included (Andrews et al. 2002a). Within a particular PC, any SPI that had absolute values of loadings within the range of 10% of the highest factor loading were kept. The following steps depend on the number of identified properties:

1. If there is only one property identified, the property is included in the selected indicators.

2. If there are two properties identified, a correlation analysis is performed. If the properties are uncorrelated (|r| < 0.6), both properties are included in the selected indicators. (Andrews et al. 2002a). However, if both properties are correlated, the property that has the higher loading is included.

3. If there are more than two properties identified, a correlation analysis is also performed. Properties that are uncorrelated (|r| < 0.6) are included in the selected indicators. From the correlated properties, if only two are left, the one with the higher loading is included. If more than two are left, only the property that has the highest absolute values of sum of correlation coefficient is included.

Each selected indicator has a weight corresponding to the PC from which it was selected. The weights for each PC can be obtained by dividing the percent of variation explained by the PC by the sum of percent variation explained of all selected PC (eigenvalues ≥ 1). Note that an indicator can be selected in more than one PC. This happens when there is a non-linear relationship with the rest of the properties and a single PC is not capable of capturing all the variation. We could not find examples of this in previous studies. Additionally, unlike Andrews et al. (2002b), the final weights were normalised to sum to one, instead of having an unbounded score, greater than one.

The final estimated utility of a role, (U^role*) then corresponds to the weighted sum of all the selected indicators (Eqn 1).

(1)U^role*=∑i=1nWi× U ^i*

Where U^role* is the estimated utility of soil for the given role, W_i is the weight of indicator i, and U^* is the utility of indicator i to support a particular certain role.

In the ordination case, obtaining best-fit curves for each of the indicators corresponds to generating a series of univariate regressions to predict the utility. Each indicator can then generate its own estimate of the utility and all the estimates are combined using a weighted average where the weights are obtained by PCA. This approach helps to rationalise the effect of each individual indicator but disregards interactions between indicators.

Regression

To overcome the limitations of the ordination method, the use of multivariate regression approach using a neural networks model is explored. In this context, neural networks present two desirable properties: (a) they can fit non-linear relationships between indicators and the utility, and (b) they can also consider interactions between indicators (Fig. 2 – section IVb).

Here, two types of regression methods were explored: (1) with the complete set of PI and (2) only using indicators selected by the F-test screening. Since the number of observations was limited, a small neural network architecture was utilised with two hidden-layers and tanh activation functions. The number of neurons was set depending on the number of inputs, with 15 and five neurons for the two hidden layers in the case of using the 16 PI and four or five, and two for the hidden layers in the case of only using selected indicators. For the output of the model, a sigmoid activation function, which limits the output to the 0–1 range, was used. The models were trained for 200 epochs, using a learning rate of 0.0005 and a batch size of 10 samples. The hyper-parameters were found using a grid search using two randomly selected sites as validation set. The model was also trained with early stopping in place to avoid overfitting.

In the case of regression, it is not always possible to obtain explicit weights for each indicator as in the case of ordination (PCA). If that is a priority, multiple linear regression can be utilised to obtain weights for each indicator. In the case of neural networks, the final estimated utility U^role* is predicted directly using the internal weights of the neurons which are adjusted during the training process.

Model evaluation

The performance of the ordination and regression model was evaluated using coefficient of determination (R²), and root mean square error (RMSE) based on the agreement of the estimated utility from the surrogates and those from PI and SPI. R² represents the proportion of total variance in the target variable explained by the model. It assesses how well the utility is approximated by model predictions. RMSE utilises the square root of the sum of the squares of the residuals to evaluate the performance of the model, where it is adjusted by the number of observations. It describes how close the estimated utility from the model is to the observed utility from the surrogates.

Implementation

The data analyses were conducted in Python (v3.6.9; Python Software Foundation 2021) using packages from matplotlib (v3.6.2; Hunter 2007), statsmodels (v0.13.5; Seabold and Perktold 2010), scikit-learn (v1.2.0; Pedregosa et al. 2011), and Tensorflow (v2.4.1; Abadi et al. 2016).

Case study 1. Function: a habitat for, and of, biodiversity

Fitting a utility graph for the surrogate

The utility graph of Shannon index shows increasing utility as the Shannon index increases, as a diverse microbial community is desired. It is also assumed that this relationship is sigmoidal, reflecting the idea that synergistic interactions between soil fauna improves the utility up to a certain point. Overall, the sigmoidal shape (Eqn 2) of our utility graph covers the range in our data and has the inflection point at the mid-point (Fig. 4).

(2)Ubiodiversity*=11+e−H′+3.5

Fig. 4.

Utility (U*) graph for Shannon index (H′) as a surrogate for soil as a habitat of biodiversity function.

Fitting utility graphs and screening

After selection of the surrogate property, the remaining datasets are the PI. Upon screening with the F-test, only six models fit the population of their respective property well enough to be passed on to the next step: Mg, EC, exch. Na, OC, S, and pH (Fig. 5). The following variables were thus not part of the subsequent PCA: total N, P, Cu, Ca, Mn, Ni, K, and Zn, as well as exch. Al and CEC.

Fig. 5.

Screened potential indicators to estimate the utility (U*) of soil as a habitat of biodiversity function.

Weighting using ordination (PCA)

Out of the six SPI, Table 1 shows the PCA indicators selected for the utility of biodiversity, in decreasing order of importance: EC, OC, pH, and Mg. For this area, Na and S were highly correlated with the selected indicators, hence redundant to assess biodiversity.

Table 1.Standardised weights for the utility (U*) of screened potential indicators identified using the principal component analysis for the soil function: a habitat for, and of, biodiversity.

Property	$U_{\text{EC}}^{*}$	$U_{\text{OC}}^{*}$	$U_{\text{pH}}^{*}$	$U_{\text{Mg}}^{*}$	$U_{\text{Na}}^{*}$	$U_{\text{S}}^{*}$
Weight	0.53	0.22	0.17	0.09	–	–

Overall, the PCA led to overestimation of low utility scores (Fig. 6d) of biodiversity with R² = 0.35 and RMSE = 0.23. In the case of biodiversity, the final index can be obtained by the weighted sum of the indicators (Eqn 3):

(3)U^biodiversity*=0.53×UEC*+0.22×UOC*+0.17×UpH*+0.09×UMg*

Fig. 6.

Comparison of utility graphs (U*; blue circle) and predicted utility (U^ *; orange triangle) for soil as a habitat of biodiversity using (a) ordination by principal component analysis of screened potential indicators (SPI), (b) regression of all potential indicators (PI), and (c) regression of SPI. Model agreement between utility graphs and predicted utility from (d) principal component analysis of SPI, (e) regression of all PI, and (f) regression of SPI.

Case study 2. Function: a store and regulator of nutrients

Fitting a utility graph for the surrogate

Most elements are plant available within a certain range around the pH of 7 while the same elements may become unavailable towards both extremes of the pH scale or become toxic at those pH extremes. As such, a maximal utility graph best describes the utility of pH. A Gaussian model was fitted with the parameters of μ = 6.5 and σ = 1. This produced a similar utility graph to the scoring function devised by Andrews et al. (2002b) (Fig. 7).

Fig. 7.

Utility (U*) graph for pH as a surrogate for soil as a store and regulator of nutrients.

Fitting utility graphs and screening

Upon screening with the F-test, only five models fit the population of their respective property well enough to be passed on to the next step: total N, exch. Al, EC, total S, and CEC (Fig. 8). The following variables were thus not part of the proceeding PCA: total OC, P, Cu, Ca, Mn, Ni, K, and Zn, as well as exch. Na.

Fig. 8.

Screened potential indicators to estimate the utility (U*) of soil as a store and regulator of nutrients.

Weighting using ordination (PCA)

Out of the five SPI, Table 2 shows the PCA indicators selected for the utility of nutrient storage, in decreasing order of importance: S, Al, N, and CEC. For this area, EC is not an important factor for biodiversity.

Table 2.Standardised weights for the utility (U*) of screened potential indicators identified using the principal component analysis for the soil function: for store and regulator of nutrients.

Property	$U_{\text{S}}^{*}$	$U_{\text{Al}}^{*}$	$U_{\text{N}}^{*}$	$U_{\text{CEC}}^{*}$	$U_{\text{EC}}^{*}$
Weight	0.52	0.25	0.15	0.07	–

Overall, the PCA led to overestimation of the utility of nutrient storage and regulation at pH below 5.5, underestimation of the utility between 5.5 and 7.5, and overestimation of utility above 7 with R² of 0.39 and RMSE of 0.22 (Fig. 9d).

Fig. 9.

Comparison of utility graphs (U*; blue circle) and predicted utility (U^ *; orange triangle) for soil as a store and regulator of nutrients using (a) principal component analysis of screened potential indicators (SPI), (b) regression of all potential indicators (PI), and (c) regression of SPI. Model agreement between utility graphs and predicted utility from (d) principal component analysis of SPI, (e) regression of all PI, (f) and regression of SPI.

Similar to the previous case study, the final index (Eqn 4) for nutrient cycling and storage can be obtained by the weighted sum of the indicators:

(4)U^nutrient*=0.52×US*+0.25×UAl*+0.15×UN*+0.07×UCEC*