Global monitoring of soil multifunctionality in drylands using satellite imagery and field data

Abstract

Models derived from satellite image data are needed to monitor the status of terrestrial ecosystems across large spatial scales. However, a remote sensing-based approach to quantify soil multifunctionality at the global scale is missing despite significant research efforts on this topic. A major constraint for doing so is the availability of suitable global-scale field data to calibrate remote sensing indicators (RSI) and, to a lesser extent, the sensitivity of spectral data of available satellite sensors to soil background and atmospheric conditions. Here, we aimed to develop a soil multifunctionality model to monitor global drylands coupling ground data on 14 soil functions of 222 dryland areas from six continents to 18 RSI derived from a time series (2006–2013) Landsat dataset. Among the RSI evaluated, the chlorophyll absorption ratio index was the best predictor of soil multifunctionality in single-variable-based models (r = 0.66, P < 0.01, NMRSE = 0.17). However, a multi-variable RSI model combining the chlorophyll absorption ratio index, the global environment monitoring index and the canopy-air temperature difference improved the accuracy of quantifying soil multifunctionality (r = 0.73, P < 0.01, NMRSE = 0.15). Furthermore, the correlation between RSI and soil variables shows a wide range of accuracy with upper and lower values obtained for AMI (r = 0.889, NMRSE = 0.05) and BGL (r = 0.685, NMRSE = 0.18) respectively. Our results provide new insights on assessing soil multifunctionality using RSI that may help to monitor temporal changes in the functioning of global drylands effectively.

Introduction

Drylands, areas with a precipitation/potential evapotranspiration ratio below 0.65 (Huang et al., 2016), are essential for sustaining life on our planet, as they cover around 42% of the global land surface, produce 42% of the world's food and host 30% of the world's endangered species (Gaur & Squires, 2017). However, drylands are threatened by climate change and desertification (Burrell et al., 2020), which can induce abrupt changes in their structure and functioning. These changes have been associated with increased aridity conditions (Berdugo et al., 2020) or reduced soil fertility and multifunctionality (Berdugo et al., 2017). Soil multifunctionality is understood as the ability of soils to maintain several ecosystem functions and services simultaneously (Garland et al., 2021). Consequently, it is crucial to monitor attributes of ecosystems, such as soil multifunctionality, in order to anticipate sudden changes that may be brought about by land degradation and the effects of climate change.

Earth observation satellites are critical for monitoring temporal trends in ecosystem attributes across global drylands. Optical sensors with coarse spatial resolution, such as the National Oceanic and Atmospheric Administration Advanced Very High-Resolution Radiometer or satellite passive microwave observation, have provided valuable information on quantifying dryland biomass at the regional scale (Tian et al., 2016). However, empirically validating the data from these sensors is challenging because it requires measuring similar areas to their pixel sizes (>10 km²). Broad-scale high-temporal frequency satellite data such as Landsat or MODIS have played an essential role in monitoring dryland vegetation dynamics. They have extensive spatial coverage and frequent observations, making them useful for this purpose. Landsat has been particularly successful in monitoring dryland vegetation attributes, with reliable accuracy in retrieving fractional cover and leaf area index at the regional scale (Sonnenschein et al., 2011; Sun, 2015). One of the methods most widely used to infer vegetation attributes has been the calculation of remote sensing indicators (RSI), such as the normalized difference vegetation index (NDVI; Rouse et al., 1974). However, NDVI applicability on a global scale is limited due to the spectral influence of mixed sparse vegetation and bare soil (Huete & Jackson, 1987).

Remote sensing indicators that minimize soil background have recently received considerable attention. Indices such as the soil-adjusted vegetation index (SAVI; Huete et al., 1985), the optimized soil-adjusted vegetation index (OSAVI; Rondeaux et al., 1996), the atmospherically resistant vegetation index (ARVI; Kaufman & Tanre, 1992), the modified chlorophyll absorption in reflectance index (MCARI; Haboudane et al., 2004) or the global environment monitoring index (GEMI; Ren & Feng, 2015) are more resistant than NDVI to saturation, background reflectance conditions and atmospheric effects. For instance, ARVI has a similar dynamic range to NDVI, but on average, it is four times less sensitive to atmospheric effects than NDVI (Kaufman & Tanre, 1992). However, the sensitivity of vegetation indices has mainly been studied at local and regional scales, and no studies have evaluated their suitability across global drylands. For instance, some studies show that traditional indices like NDVI perform better than modified vegetation indices for monitoring above-ground green biomass in arid and semi-arid grasslands (Ren & Feng, 2015). In contrast, other studies show that the modified vegetation indices, such as the SAVI and L-SAVI, improved the detection of spatio-temporal changes in the vegetation in a semi-arid area (Fatiha et al., 2013). These examples underscore the current deficiency in assessing and comparing existing vegetation indices at the global scale.

Developing global models of soil multifunctionality faces significant challenges, one of which is the lack of suitable field data that fits the spatial resolution of satellite imagery required to build and validate these models. To overcome this challenge, we need to move beyond remote sensing indicators (RSI) mainly related to plant cover and incorporate other indicators that can potentially analyse biophysical properties such as plant composition and functioning. For instance, Zhao et al. (2018) demonstrated a significant relationship between visible black-sky albedo and soil multifunctionality across global drylands. However, this study was limited to a selection of 61 homogeneous plots from the 224 dryland datasets compiled by Maestre et al. (2012) to avoid the mismatch between field data collected from 30 × 30 m plots and MODIS image resolution of 500 × 500 m (NASA LP DAAC, 2017). Such limitations can be overcome by combining existing global datasets of ground-collected soil data collected within 30 × 30 m plots (Delgado-Baquerizo et al., 2013; Maestre et al., 2012) with spectral data provided by Landsat 7 ETM. Landsat satellites offer the highest spatial resolution in the thermal region, with freely available imagery and the longest temporal record (roughly every 16 days) spanning the last 49 years (Wulder et al., 2022).

In comparison, new missions such as Sentinel-3 cover this region with a 1000-m spatial resolution. Landsat features spatial resolutions ranging from 15 to 80 m in the visible and infrared region and between 60 to 120 m in the thermal region, depending on the specific Landsat mission (Holden & Woodcock, 2016). Integrating medium spatial resolution image data (30 × 30 m) such as Landsat data and field-based observations could improve the assessment of soil multifunctionality worldwide in a cost-effective and accessible manner. Furthermore, the high-performance computational capacities of Google Earth Engine can be used to access and process satellite data on the cloud, providing new possibilities for analysing large volumes of data globally (Gorelick et al., 2017).

Here, we combine the use of 18 surface reflectance vegetation indices and thermal remote sensing-based indices, hereafter called remote sensing indicators (RSI), with a global assessment of 14 soil functions measured in situ in 222 dryland ecosystems on six continents (Maestre et al., 2012; Ochoa-Hueso et al., 2018). Our study has two main objectives: first, to evaluate the sensitivity of remote sensing indicators (RSI) in characterizing soil multifunctionality in dryland ecosystems worldwide, and second, to develop models to upscale ground-based observations with RSI data. By achieving these objectives, we aim to provide robust and comprehensive models that can enhance our understanding of dryland ecosystems' current status and dynamics globally.

Materials and Methods Study sites

Field data were gathered from 222 sites in 19 countries (Argentina, Australia, Brazil, Chile, China, Ecuador, Iran, Israel, Kenya, Mexico, Morocco, Peru, Spain, Tunisia, USA, Venezuela, Botswana, Burkina Faso and Ghana; Figure 1). These sites are a subset of the 236 sites used in Ochoa-Hueso et al. (2018); we had to exclude 14 sites due to the lack of cloud-free images during the inventories at these sites. The 222 sites surveyed covered all major vegetation/soil types and the wide range of environmental conditions across global drylands (UNEP-WCMC, 2007; Figure S1).

Figure 1. Distribution of the dryland sites used in this study. Dryland land areas are displayed in orange according to FAO/UNEP Land Cover Classification System (UNEP-WCMC, 2007).

At each site, field data were collected from 30 × 30 m plots between February 2006 and December 2013 using a standardized protocol described in detail in Maestre et al. (2012). Plant cover data were obtained from four 30-m-long transects using the line-intercept method (Tongway & Hindley, 2004). Soil samples were collected using a stratified procedure of five 50 × 50 cm quadrats randomly placed under the dominant perennial vegetation patch type and in open areas devoid of perennial vegetation. Five soil cores extracted at 0–7.5 cm depth from each quadrat were bulked and homogenized in the field. The composite included samples for microsites in open areas devoid of perennial vegetation and under the canopy of the dominant perennial plant species. The number of soil samples collected varied between 10 and 15 per site (depending on whether one or two dominant plant species were found at each site), accounting for more than 2600 samples. After field collection, the soil samples were taken to the laboratory, where they were sieved (2 mm mesh), air-dried for 1 month and stored for laboratory analysis. Dry soil samples were then analysed for soil functions related to the cycling and storage of carbon (Table 1), as described in Maestre et al. (2012). A plot-level estimate of all the soil functions was obtained using a weighted average of values from open and vegetated microsites weighted by their respective cover at each plot. As a soil multifunctionality index, we used the average Z-score for all soil functions estimated at the plot level (Maestre et al., 2012).

TABLE 1 Soil variables used for the calculation of the soil functions.

Total nitrogen	TON
Available potassium	AVP
Activity of b-glucosidase	BGL
Activity of phosphatase	FOS
Organic carbon	ORC
Ammonium	AMO
Nitrate	NIT
Aminoacids	AMI
Proteins	PRO
Phenols	PHE
Aromatic compounds	ARO
Hexose content	HEX
Pentose content	PEN
Potential N transformation rate	NTR

Satellite data and processing

We used Landsat 5 TM and Landsat 7 ETM+ filtered products to obtain spectral data from 927 available images collected between 2006 and 2013. First, we selected images taken as close as possible (within 1–3 months) to the day when field surveys were conducted for the 236 sites analysed in Ochoa-Hueso et al. (2018). We then used a local filter to select the cloud-free Landsat images closest to the field surveys, reducing the dataset to 222 sites.

To correct for atmospheric gases and aerosols, which can vary in space and time and can significantly impact Landsat spectral data collected on different dates (Masek et al., 2006; Roy et al., 2014), we atmospherically corrected the Landsat imagery using the Landsat ecosystem disturbance adaptive processing system (LEDAPS, version 3.4.0; Schmidt et al., 2013). We also corrected surface reflectance to account for data from plots measured on different dates. To retrieve surface temperature, we used a methodology proposed by Jimenez-Muñoz et al. (2009) that employs a single-channel algorithm using the thermal-infrared Landsat channel (band 6). This algorithm conducts emissivity and atmospheric corrections to retrieve the surface-level temperature. The algorithm has been extensively validated in other independent studies across various land covers, showing an RMSE between 1 and 2 K, which is typical accuracy for remotely sensed land surface temperature products (Copertino, 2012; Zhang, He, et al., 2016). First, we estimated surface emissivity using a simple approach based on fractional vegetation cover and NDVI (Sobrino et al., 2008). We then calculated atmospheric functions from total atmospheric water vapour values obtained from the Copernicus Climate Change Service implemented by the European Centre for Medium-Range Weather Forecasts (Muñoz-Sabater et al., 2021). The single-channel algorithm needed these data to be applied (Hersbach et al., 2020). We applied a 30-m buffer to extract data from each study area and weighted the value of each pixel covered in this area to minimize errors in the geolocation and referencing of each pixel. An average of the pixel weights by the percentage of the area overlapping each plot was used to ensure an extensive, systematically collected sample scheme.

We obtained spectral reflectance for each TM and ETM+ reflective band and surface temperature, which we used to calculate the RSI dataset for each location. To estimate a wide range of soil and plant traits (Hernandez-Clemente et al., 2021), we evaluated a list of 18 RSI, including formulations based on the near-infrared and visible regions (NDVI, GLI, SIPI), modified vegetation indices proposed to minimize background and soil effects (SAVI, TSAVI, OSAVI, TSAVI/OSAVI, MCARI2 and GEMI), modified vegetation indices considering atmospheric corrections (ARVI, AFRI and VARI), formulations based on the short-wave infrared (SWIR) bands (S1260 and NBR2) and thermal bands (Ts-Ta) and WDI. For definitions and descriptions of acronyms, refer to Table 2.

TABLE 2 Remote sensing indicators and their formulations derived from Landsat data evaluated in this study.

Remote sensing indicator	Formulation	References
GLI: green leaf index	GLI = (2·ρGreen – ρRed − ρBlue)/(2·ρGreen + ρRed + ρBlue)	(Louhaichi et al., 2001)
SIPI: structure insensitive pigment index	SIPI = (ρNIR − ρBlue)/(ρNIR + ρRed)	(Peñuelas et al., 1995)
NDVI: normalized difference vegetation index	NDVI = (ρNIR – ρRed)/(ρNIR + ρRed)	(Rouse et al., 1974)
SAVI: soil-adjusted vegetation index	SAVI = (1 + L) × ((NIR-Red)/(NIR + Red + L))	(Huete, 1988)
TSAVI: transformed soil-adjusted vegetation index	TSAVI = (a × (ρNIR – a × ρRed − b))/(ρRed + a × ρNIR − (a × b) + 0.08 × (1 + a2))	(Baret & Guyot, 1991)
OSAVI: optimized soil-adjusted vegetation index	OSAVI = (ρNIR − ρRed)/(ρNIR + ρRed + 0.16)	(Rondeaux et al., 1996)
TSAVI/OSAVI	TSAVI/OSAVI	(Baret & Guyot, 1991)
MSAVI: modified soil-adjusted vegetation index	MSAVI = 0.5 × (2 × ρNIR + 1 − (((2 × ρNIR + 1)2)0.5 – 8 × (ρNIR − ρRed)))	(Qi et al., 1994)
MCARI2: modified chlorophyll absorption ratio index 2	MCARI2 = (1.5 × (2.5 × (ρ800 − ρ670) – 1.3 × (ρ800 − R550)))/((2 × ρ800 + 1)2 − (6 × ρ800 – 5 × (ρ670)0.5) − 0.5)0.5	(Haboudane et al., 2002)
EVI: enhanced vegetation index	EVI = 2.5 × (ρNIR − ρRed)/(ρNIR + 6 × ρRed − 7.5 × ρBlue + 1)	(Huete et al., 2002)
GEMI: global environment monitoring index	GEMI = n (1–0.25 × n) − (Red − 0.125)/(1 − Red); n = 2 (NIR2 − Red2) + (1.5 × NIR + 0.5 × Red)/(NIR + Red + 0.5)	(Pinty & Verstraete, 1992)
ARVI: atmospherically resistant vegetation index	ARVI = (ρNIR − ρRed) – γ (ρRed-ρBlue)/(ρNIR + ρRed) – γ (ρRed − ρBlue)	(Kaufman & Tanre, 1992)
AFRI2100: aerosol-free vegetation index 2100	AFRI2100 = (ρNIR − 0.5 × ρ2100)/(ρNIR + 0.5 × ρ2100)	(Karnieli et al., 2001)
VARI: visible atmospherically resistant index	VARI = (ρGreen − ρRed)/(ρGreen + ρRed − ρBlue)	(Gitelson et al., 2002)
S1260: sulphur index 1260	S1260 = (ρ1260 − ρ660)/(ρ1260 + ρ660)	(Mahajan et al., 2014)
NBR2: landsat normalized burn ratio 2	NBR2 = (SWIR1 − SWIR2)/(SWIR1 + SWIR2)	(Norton et al., 2009)
Ts-Ta: surface temperature minus air temperature	Ts-Ta	(Jackson et al., 1981)
WDI: water deficit index	WDI	(Moran et al., 1994)

Modelling approach

In this study, we investigated the capacity of the RSI evaluated to predict soil multifunctionality across global drylands. The first step in data analysis entailed selecting the most significant RSI for determining soil multifunctionality, followed by an evaluation of the model's performance, as depicted in Figure 2. To ensure the interpretability of our model output, we first reduced the number of variables by using a filter-based feature selection approach (Gosiewska et al., 2021). We excluded RSI that were highly correlated with each other (r > 0.85, (Dormann et al., 2013)) and used only those with a variance inflation factor lower than 10 (Kutner et al., 2004). This resulted in an RSI selection referred to as RSI-mc. We then performed a principal component analysis (PCA) to interpret the contribution of each RSI based on the first two principal components with an importance higher than 15% (Pacheco et al., 2013). We identified the loading vectors in the biplot of the principal components explaining >60% of the variance, which were used to select the RSI variables with the highest eigenvalues per axis, resulting in an RSI selection referred to as RSI-pca. Lastly, we included a single-index selection using the RSI most correlated with soil multifunctionality, referred to as the 1-RSI model, to check the improvement achieved with the variable reduction approaches (RSI-mc and RSI-pca).

Figure 2. Data analysis workflow – Remote sensing indicators (RSI) selection and model performance evaluation for soil multifunctionality determination.

To evaluate the suitability of our model, we compared three different approaches: two artificial intelligence methods, an evolutionary algorithm model (EAM) and a random forest model (RF), and simple linear regression (LR). The comparison between these three methods was made to analyse the suitability of each approach, which varies with the variability and dispersion of the data (Franklin & Miller, 2010). The EAM model was based on a genetic algorithm used to generate high-quality solutions to optimize model accuracy in computer science, known as evolutionary algorithms (Vikhar, 2016). We computed the models with the Eureqa software package v1.24 (Datarobot Inc., Boston, USA), combining the arithmetic, trigonometric and exponential building blocks for the best model accuracy. Eureqa uses evolutionarily search to determine the best predictive models, simplifying the final calculated model. The RF model was built by using the Caret library (Kuhn et al., 2020) within the R environment (R Core Team, 2013) and the package ‘caret’. The adjustment parameter mtry (Randomly Selected Predictors) was established initially by iterating over the whole range of values. Then, a pre-processing transformation was applied by centring and scaling the training data. Comparatively, we also tested a simple method based on a linear fitting (Freedman, 2009) between RSI and soil multifunctionality. Finally, the models were trained considering a resampling method of five k-folds and three repetitions. The model accuracy was evaluated with a cross-validation bootstrap procedure (Austin & Tu, 2004). For doing so, data were randomly split into K = 500 sets, selecting 80% of our dataset to generate each predictive model, and the remaining 20% was set aside for validation purposes. The average and standard deviation from this cross-validation bootstrap procedure was used for validation. We calculated the R-squared (R²) and the normalized root-mean-square error (NRMSE) by contrasting predicted versus observed values.

Results Soil multifunctionality RSI determination

The feature reduction simplified the number of variables included in the modelling process. For example, the first reduction, RSI-mc, resulted in a list of 10 RSI: MCARI2, NBR2, MSAVI, GLI, S1260, AFRI22, TSAVI_OSAVI, GEMI, Ts-Ta and WDI. In the PCA biplot, we observed that these RSI were grouped into three clusters of loading vectors, each associated with climate (Figure 3A) and vegetation type (Figure 3B) and enclosed by concentration ellipses.

Figure 3. Standardized principal components (PC1 vs. PC2) plot of the 10 remote sensing indicators (RSI) less correlated among them. Biplot vectors are RSI loadings, whereas the position of the 222 sites is shown within each climate (arid, semi-arid and dry-subhumid, A) and vegetation type (grasslands, shrublands, open forests with shrubs and savannahs, B).

The modified vegetation indices MCARI2, MSAVI, GLI, TSAVI_OSAVI, S1260, NBR2 and AFRI22 were negatively correlated with the same principal component (PC1), while GEMI was positively related to PC1. On an orthogonal axis, the third group of vectors, Ts-ta and WDI, were the best contributors to PC2. In the PCA biplot, the contribution of WDI and Ts-Ta was quite similar, with the eigenvalue being slightly higher for Ts-Ta. However, it should be noted that WDI was significantly more correlated than Ts-Ta with the soil functions evaluated (Figure 4).

Figure 4. Pearson correlation coefficients between soil multifunctionality (M) and individual soil functions (Table 1) and the different vegetation indices used (Table 2); P > 0.05 value are shown with an ‘x’ symbol (n = 222).

The ellipses in the Standardized Principal Components (PC1 vs. PC2) plot serve as visual representations of the distribution and variability of data points associated with the same climate (Figure 3A) and vegetation type (Figure 3B). In these plots (Figure 3A and B), large ellipses centred within the graph represent semi-arid climates, grasslands and open forests. These ellipses demonstrate a high degree of variability in the data for each cluster and a consistent representation across the 10 RSI-mc selected. Conversely, certain ellipses are associated explicitly with particular RSI. For example, GEMI and water stress indicators (Ts-Ta and WDI) contribute more to sites in arid areas. At the same time, GEMI and soil-adjusted, atmospherically resistant RSI (MCARI2, MSAVI, GLI, TSAVI_OSAVI, S1260, NBR2 and AFRI22) are more prominent in dry-subhumid areas (Figure 3A and Figure S1). In the standardized principal components plot of vegetation types (Figure 3B), GEMI and water stress indicators mainly represent savannahs. Shrublands form a cluster to the left, characterized by soil-adjusted, atmospherically resistant RSI and water stress indicators (Figure 3B).

We selected one RSI per group with the highest eigenvalues from the three main groups of eigenvectors in the standardized principal components (PC1 vs. PC2) plot (Figure 3). As a result, RSI-pca selection reduced the predictors to MCARI2, GEMI and Ts-Ta. Finally, we compared the RSI reductions, RSI-mc and RSI-pca, to the 1-RSI with the highest correlation for estimating soil multifunctionality. According to Figure 4, MCARI2 strongly correlates with soil multifunctionality (R = 0.54, R² = 0.28, P < 0.01).

Global models of soil multifunctionality

The scaling-up of soil multifunctionality across a wide range of climates and vegetation types with EAM, RF and LM models showed R² and NMRSE values ranging from 0.17 to 0.55 and from 0.25 to 0.15, respectively, which depended on the type of model and number of RSI considered (Figure 5A). The highest accuracy for soil multifunctionality was obtained with models computed with EAM using the RSI-pca selection, which improved NRMSE by 25% from RF and 27% from LM models respectively (Figure 5B). The EAM analysis reduced the NMRSE in soil multifunctionality estimations for the three variable selection methods followed (RSI-pca, RSI-mc and 1-RSI). The EAM-driven analysis utilizing MCARI2 resulted in a 22% reduction in NMRSE compared to the linear analysis derived from RSI-pca. Furthermore, it was observed that employing models based on RSI-mc was unnecessary, as the RSI-pca produced the most reliable outcomes when processed through EAM (Figure 5A).

Figure 5. (A) Accuracy results in the prediction of soil multifunctionality according to the coefficient of determination (R2) and the normalized root mean square error (NMRSE). Results are shown for the linear model (LM), random forest (RF) and evolutionary algorithm model (EAM) modelling approaches used, and for the three variables reduction sets used: RSI-mc based on MCARI2, NBR2, MSAVI, GLI, S1260, AFRI22, TSAVI_OSAVI, GEMI, Ts-Ta and WDI; RSI-pca based on MCARI2, GEMI and Ts-Ta; and 1-RSI based on MCARI2−. (B) Observed versus predicted soil multifunctionality for the RSI-pca selection with the EAM model using the remote sensing indicators MCARI2, GEMI and Ts-Ta and EAM analysis (n = 222). The dashed line represents the 1∶1 line. See Table 2 for the acronyms of the indices used.

The accuracy of the soil multifunctionality model based on RSI-pca using EAM analyses (r = 0.733, NMRSE = 0.15) also shows consistency across soil functions, with upper and lower values obtained for AMI (R = 0.889, NMRSE = 0.05) and BGL (R = 0.685, NMRSE = 0.18) respectively (Figure 6).

Figure 6. Comparison of the accuracy of EAM-based predictions of soil multifunctionality (M) and soil functions (TON, AVP, BGL, FOS ORC, AMO, NIT, AMI, PRO, PHE, ARO, HEX, PEN and NTR described in Table 1). The results are shown in terms of the correlation coefficient (R) represented with blue columns and the normalized root mean square error (NMRSE) represented with red error bars. The predictions were made using three remote sensing indicators (MCARI2, GEMI and Ts-Ta) selected through principal component analysis (RSI-pca).

Discussion

We developed and validated a model to estimate soil multifunctionality across global drylands using a comprehensive global field survey and satellite imagery. Our results highlight the reliability of RSI, such as MCARI2, NBR2, MSAVI, GLI, S1260, AFRI22, TSAVI_OSAVI or GEMI, to model dryland soil multifunctionality. These RSI were developed to reduce the influence of soil background and atmospheric effects on the regions with low-density vegetation cover. In contrast, other simpler RSI formulations, such as the widely used NDVI, showed about 50% lower coefficient of determination values than MCARI2 (see Figure S2) to model soil multifunctionality. These results may be related to a higher sensitivity of NDVI to soil brightness effects or to the presence of senesced vegetation and standing litter (Baret et al., 1993). However, NDVI could still provide complementary information to modelling soil multifunctionality with MCARI2 as a proxy for primary production (Prince, 1991; Tucker et al., 1983) or ecosystem structure and functioning (Gaitán et al., 2013) over large scales. In contrast to studies mainly based on the relationship between fractional vegetation cover and NDVI (Song et al., 2017), here we developed a model to monitor soil multifunctionality, a key feature of dryland ecosystems (Maestre et al., 2016).

The models based on the combination of MCARI2, GEMI and Ts-Ta (RSI-pca) improved the accuracy in estimating soil multifunctionality compared to models using a single predictor. While Ts-Ta alone cannot serve as a reliable indicator of soil multifunctionality, it can complement other RSI, such as MCARI2 and GEMI, to enhance the accuracy of global models. When combined with MCARI2 and GEMI, Ts-Ta improves the accuracy of soil multifunctionality models, as demonstrated by a 12% decrease in NMRSE, and enhances the prediction of specific soil functions by 8%–18% (Table S1). This outcome may be attributed to the combination of RSI-pca selections (MCARI2, GEMI and Ts-Ta), which more accurately represent the variability observed in drylands distributed across diverse climates and vegetation types worldwide. In addition to the strong correlations found between MCARI2, GEMI and Ts-Ta and soil multifunctionality, further analyses using structural equation modelling (Table S2 and Figure S3) demonstrate that these RSI provide the most reliable models for estimating soil multifunctionality in drylands. These findings align with prior research, indicating that MCARI2 and GEMI indices (Haboudane et al., 2004; Pinty & Verstraete, 1992) exhibit lower sensitivity compared to other vegetation indices such as NDVI in detecting fractional cover variations ranging from 2% to 83% across analysed dryland locations (Maestre et al., 2012).

The RSI-pca model improves estimates of soil multifunctionality and individual soil functions. While the correlations of the 1-RSI model are significantly related to most of the soil functions (TON, BGL, ORC, PRO, PHE, ARO, HEX and NTR), the RSI-pca model can generate models with errors of only 5–18% for all variables. This can be explained by the indirect relationship that many soil components have with different types of variables, such as variations in biomass, soil moisture and primary production (Liu et al., 2023). The combined use of a model that absorbs this variability can reflect the specific variations of these compounds, as demonstrated in this work for AVP, NIT, AMI and PEN.

Our study emphasizes the importance of avoiding models based solely on best-fitting indices (Hornero et al., 2021). The PCA reduction method improves the results' interpretability by evaluating the RSI loading vectors used to assess soil multifunctionality and functions per climate and vegetation type. Among the selected RSI, MCARI2 and GEMI are used as a proxy for fractional cover (Haboudane et al., 2002; Pinty & Verstraete, 1992), where MCARI2 reduces the RSI's sensitivity to soil and background effects, and GEMI minimizes the impact of undesirable atmospheric perturbations. Additionally, Ts-Ta provides an indicator of the water stress condition of the vegetation linked to stomatal conductance and transpiration (Morillas et al., 2013), contributing to representing semi-arid dryland sites.

Our study provides compelling evidence that EAM methods are a reliable tool for accurately upscaling ground-based observations of soil multifunctionality on a global scale. The EAM models developed in this study showed significant improvement in NRMSE values by 37% and 33%, respectively, compared to RF and LR models for quantifying soil multifunctionality (Table S1). Furthermore, the accuracy obtained for predicting soil multifunctionality using the 1-RSI (r = 0.66, P < 0.01) and RSI-pca (r = 0.73, P < 0.01) models with Landsat data and EAM models represents a significant improvement compared to results from previous studies. For instance, Zhao et al. (2018) reported a correlation between soil multifunctionality and MODIS land surface albedo of only r = −0.314. These findings align with recent efforts to apply deep learning approaches to quantify soil organic carbon composition at the national level, as reported by Odebiri et al. (2022). These results demonstrate the potential of EAM models for providing reliable estimates of soil multifunctionality and support their application for global-scale monitoring and management of soil resources in drylands.

Biocrusts are essential components of drylands globally, significantly regulating their structure and functioning (Bowker et al., 2013; Maestre et al., 2011, 2013). Biocrusts fix substantial amounts of atmospheric CO₂ (over 2.6 Pg of C per year; Elbert et al., 2012) and impact the temporal dynamics of soil CO₂ efflux and net CO₂ uptake. Additionally, biocrusts influence soil enzyme activity (Miralles et al., 2012), nitrification (Castillo-Monroy & Maestre, 2011) and runoff-infiltration rates (Zaady et al., 2013), all of which contribute to soil multifunctionality. Remote sensing provides a valuable and reliable method for mapping biocrusts. Nevertheless, due sto the spectral resemblance between predominant dryland surface elements and biocrusts, it is necessary to utilize mapping indices based on hyperspectral data to identify areas dominated by biocrusts at the ecosystem level accurately (Rodríguez-Caballero et al., 2017). This limitation hinders the ability of most satellite imagery products, such as Sentinel, Landsat or MODIS, to effectively detect biocrusts (Rozenstein & Adamowski, 2017). Because of this, we could not consider biocrusts explicitly in our analyses. However, they have been shown to influence the soil functions we evaluated in drylands significantly (Bowker et al., 2011), and, as such, they could have also influenced our results. Nevertheless, we do not expect biocrusts to invalidate our results for two main reasons: (i) we measured soil functions at 0–7.5 cm depth, and biocrusts affect soil functions largely at the 0–2 cm depth (Pointing & Belnap, 2012), and (ii) the positive impacts of perennial vegetation on soil functions such as those studied here extend beyond plant canopies to influence adjacent open areas devoid of perennial vegetation (Maestre et al., 2009).

This study demonstrates the potential of Landsat images and EAM-based models to assess soil multifunctionality over large areas, but several limitations must be acknowledged. First, the temporal resolution of the sensor (one or two images per month) limits the estimations to monthly or yearly intervals, and advanced filters cannot be applied to select images with similar weather conditions within the same month. Second, the spectral resolution of the images, with spectral bands of ~30 nm on average in the VIS–NIR region, restricts the quantification of biocrusts, as discussed above, and of critical biophysical variables that evaluate the status of dryland ecosystems, such as the chlorophyll fluorescence or pigment contents of vegetation (Smith et al., 2018; Zhang, Xiao, et al., 2016). Third, the spatial resolution of the images, with pixels of 30 × 30 m, cannot capture the fine-scale spatial heterogeneity that characterizes dryland ecosystems (Smith et al., 2019), as well as that of biocrusts (Maestre & Cortina, 2002). However, new satellite missions will overcome some of these limitations. For instance, Sentinel 2 provides 13 spectral bands and a spatial resolution ranging from 10 m to 60 m, the NASA mission EMIT provides hyperspectral data from 400 to 2500 nm with a daily temporal resolution and a spatial resolution of 5 m, and the enhanced spectral resolution of the upcoming Landsat next missions. In addition, the future satellite mission FLEX will provide a single platform of a fluorescence-dedicated imager at an unprecedented spatial resolution of 300 m (Meng et al., 2022). The implementation of these new missions will enhance the ability to seamlessly integrate field data, such as those used in this study, with high-resolution indicators of photosynthetic activity and soil properties, such as texture, organic carbon and moisture. This will improve the accuracy of global models for soil multifunctionality.

Conclusions

The combined use of a unique global field dataset including 14 soil functions and a wide range of RSI calculated from Landsat has enabled us to develop a predictive model for soil multifunctionality in drylands based on three RSI: MCARI2, a soil-atmosphere resistant VI; GEMI, an atmospherically resistant VI; and Ts-Ta, a proxy of water stress conditions. Our findings demonstrate that RSI, such as MCARI2, performs better than NDVI. These findings imply that NDVI is more sensitive to the variability of global dryland covers, a crucial factor in developing comprehensive models for soil multifunctionality in sparsely vegetated regions. To the best of our knowledge, our study is the first to use and demonstrate that thermal-based indicators such as Ts-Ta, which are related to the evapotranspiration rate and water deficit, can improve global models of soil multifunctionality in combination with other RSI. Future research to improve our understanding of dryland dynamics should include EAM methods for accurately upscaling ground-based observations. The soil multifunctionality models developed in this study open the possibility of accurately mapping regional- to global-scale essential soil processes at spatiotemporal resolutions relevant to land managers across drylands worldwide.

Acknowledgements

Field data were obtained with the support of the European Research Council (ERC) grant agreement 242658 (BIOCOM). Hernández-Clemente R was supported by the Ramón y Cajal program (RYC2020-029187-I) and the State Plan for Scientific and Subprogram for Knowledge Generation (PID2021-124058OA-I00) from the Spanish Ministry of Science and Innovation (RYC2020-029187-I). Maestre FT acknowledges support from Generalitat Valenciana (CIDEGENT/2018/041) and the Spanish Ministry of Science and Innovation (EUR2022-134048).

Author Contributions

Maestre FT and Hernández-Clemente R conceived the study. Hernández-Clemente R designed the theoretical and computational framework, defined the modelling approaches, and took the lead in writing and preparing the manuscript. Hornero A contributed to the analytical framework, image processing and analytical calculations, Gonzalez-Dugo contributed to the thermal data analysis, Berdugo M and Quero JL contributed to the theoretical framework and discussion of the manuscript, Jiménez JC contributed to the thermal image pre-processing, Maestre FT contributed to the theoretical framework, provided the experimental data and had a major role in writing the manuscript.