Low sensitivity of three terrestrial biosphere models to soil texture over the South American tropics

Drought stress is an increasing threat for vegetation in tropical regions, within the context of human-induced increase of drought frequency and severity observed over South American forests. Drought stress is induced when a plant's water demand is not met with its water supply through root water uptake. The latter depends on root and soil properties, including soil texture (i.e. the soil clay and sand fractions) that determines the soil water availability and its hydraulic properties. Hence, soil clay content is responsible for a significant fraction of the spatial variability in forest structure and productivity. Soil-textural properties largely vary at the spatial resolution used by Terrestrial Biosphere Models (TBMs) and it is currently unclear how this variability affects the outputs of these models used to predict the response of vegetation ecosystems to future climate change scenarios. In this study, we assessed the sensitivity of the carbon cycle of three state-of-the-art TBMs, i.e. ORganizing Carbon and Hydrology in Dynamic EcosystEms (ORCHIDEEv2.2), Ecosystem Demography model version 2 (ED2), and Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS) to soil-textural properties at the regional level over the South American tropics using model default pedotransfer functions. For all three TBMs, the model outputs, including gross primary productivity (GPP), aboveground biomass (AGB), soil carbon content and drought stress, were shown to be mostly insensitive to soil-texture changes representative of the spatial variability in soil properties, except for a small region characterised by very low water availability in ORCHIDEEv2.2 and ED2. We argue that generic pedotransfer and simple drought stress functions, as currently implemented in TBMs, should be reconsidered to better capture the role of soil texture and its coupling to plant processes. Similarly, we suggest that better estimates of the soil-texture uncertainty resulting from soil-texture data aggregate should be considered in the future. Those steps forward are critical to properly account for future increasing drought stress conditions in tropical regions.

1 Introduction

Over the last 3 decades, the Amazon tropical forest has been facing an increase in environmental pressures, including the severity and length of drought events (Spinoni et al., 2014). This trend is projected to be exacerbated by the end of the century (Duffy et al., 2015), also resulting from the rapid deforestation rates and the regional precipitation recycling (Staal et al., 2020). Observations and manipulative field experiments have revealed a clear sensitivity of the Amazon forest to severe drought, potentially leading to large-scale increase in tree mortality and decrease in forest productivity through reduced photosynthesis (Nepstad et al., 2007; Phillips et al., 2009; Gatti et al., 2014; Doughty et al., 2015; Corlett, 2016; Feldpausch et al., 2016). Increasing tree mortality in the Amazon is thought to be induced by soil moisture deficit (low water supply) combined with low air humidity and high air temperature (high water demand), which when combined, leads to either hydraulic failure or stomatal closure that may cause carbon starvation in trees (Rowland et al., 2015). Increased vulnerability of the Amazon forest to drought stress will have large impacts on the regional and global carbon, nutrient and water cycles as well as the coupled climate system, and has been proposed as one of the factors involved in the observed decline of the Amazon carbon sink strength (Brienen et al., 2015; Maeda et al., 2015; O'Connell et al., 2018; Hubau et al., 2020).

Soil water availability is intimately related to root distribution (De Deurwaerder et al., 2021), root and soil depth (Jackson et al., 1996), as well as soil properties including texture. By modulating the retention and accessibility of water (and nutrients) to the trees (Silver et al., 2000; Laurance et al., 1999), soil texture, and especially the clay fraction, shapes forest structure, function and its spatio-temporal dynamics. At the regional level, clay and nutrient gradients were shown to explain a substantial part of the variability in forest biomass, soil carbon pools, and forest productivity across the Amazon basin (Laurance et al., 1999; Aragão et al., 2009; Jiménez et al., 2014). The intensity of the dry season and the availability of nutrients (e.g. phosphorus) affect species distribution (Condit et al., 2013; Jirka et al., 2007), while soil moisture gradients have also been shown to affect plant traits and leaf area index (LAI) (Fyllas et al., 2009; Flack-Prain et al., 2021). Furthermore, by affecting canopy conductance and hence carbon assimilation, soil moisture directly impacts the dynamics of water and carbon fluxes at the tree level (Harris et al., 2004).

In order to study the resilience of the Amazon forest to future drought and deforestation, Terrestrial Biosphere Models (TBMs) are key tools that integrate ecophysiological processes at different spatio-temporal scales and the response of ecosystems to environmental changes. In most TBMs, water availability directly affects carbon assimilation through so-called drought stress functions that modulate leaf stomatal conductance. Joetzjer et al. (2014) showed that the amplitude and timing of plant response to moisture deficit was highly sensitive to these unconstrained functions, which prevented the accurate representation of the impact of drought over the Amazon rainforest. Consequently, current TBMs are unable to simulate the spatial variability of forest productivity and biomass over the Amazon (Johnson et al., 2016). Drought stress response parameterisation and sensitivity were also shown to affect the coupling strength between the land surface and the atmospheric boundary layer and thus the performance of coupled climate models (Combe et al., 2016). To better capture the drought stress effect on vegetation, unprecedented efforts are being made in the community to improve the representation of plant hydraulics in TBMs (e.g. Christoffersen et al., 2016; Xu et al., 2016; Mencuccini et al., 2019). However, less effort has been spent on enhancing the representation of the belowground compartments, despite its key role in drought stress (Carminati and Javaux, 2020). Roots, soils, their interactions, and their effect in TBMs have not received as much attention as they should. In this study, we primarily focus on the impact of soil texture on TBMs but also call for more research into other belowground components.

In TBMs, soil moisture is determined by soil hydrology submodels which typically rely on soil-textural information and pedotransfer functions. The TBMs use a rather limited number of pedotransfer functions while soil-texture inputs often resume to a few products with different spatial resolutions (horizontally and vertically). The main products used for regional and global simulations are the FAO/UNESCO soil map of the world (Batjes, 1997), the Harmonized World Soil Database (Nachtergaele et al., 2008), and the more recent SoilGrids250m products (Hengl et al., 2017; Poggio et al., 2021). Those gridded soil maps are not independent of one another, and often aggregated at lower spatial resolution to match other model forcings (e.g. meteorological drivers) which results in average soil properties that neglect inter- and intra-grid cell heterogeneity. Yet, most TBMs do not consider the uncertainty in the influence of soil texture in vegetation activity and drought when applied at large scales.

To date, the existing evaluations of the response of TBMs to soil properties mainly focused on hydrology and water fluxes. These analyses tend to show a lack of sensitivity of TBMs to soil texture and composition. For instance, Li et al. (2012) showed that the performance of CABLE remained insensitive to parameters of soil water dynamics across three contrasting sites, even after improving critical processes related to root functioning. In line with these results, Tafasca et al. (2020) investigated the impact of soil texture on soil water fluxes and storage at different scales with ORganizing Carbon and Hydrology in Dynamic EcosystEms (ORCHIDEE), as part of the Land Surface, Snow and Soil moisture Model Intercomparison Project (LS3MIP; van den Hurk et al., 2016). They showed that, while the model exhibits realistic behaviours at the local scale, it is weakly sensitive to the choice of soil-texture maps at the global scale. The effect of soil representation on vegetation and carbon has been sporadically assessed in the literature. By applying the Ecosystem Demography model, version 2 (ED2) over the Amazon rainforest, Longo et al. (2018) showed different sensitivities of the aboveground biomass (AGB) to soil texture, depending on the rainfall regimes at two contrasting sites. While this study suggests that soil hydraulic properties mediate the Amazon ecosystem response to rainfall regimes, the authors highlighted that the current parameterisation of the model does not account for the diversity in soil types and is limited for representing certain configurations such as clay-rich soils.

Since soil is a major carbon pool and a key driver of water and nutrient availability for plants, we expect a large model sensitivity to soil properties, which should propagate into the simulated vegetation and the ecosystem biogeochemical cycles. This assumption especially applies in a tropical region like South America that frequently suffers from drought and is characterised by heavily weathered and poor soils. To test this assumption, we explored the sensitivity of the vegetation carbon dynamics to soil texture in three state-of-the-art TBMs, representative of the main classes of commonly used TBMs: Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS), ED2, and ORCHIDEE v2.2. Model sensitivity to soil texture was assessed based on the inter- and intra-grid cell variability in clay content as quantified from the SoilGrid250m database. For each simulation, we present the sensitivity of soil carbon pools, gross primary productivity (GPP), and AGB resulting from the different soil configurations, for both the conditions after the model spin-up and the historical simulation spanning the 1860–2016 period. Results of the different simulations for the three models are compared to one another and with existing observation products to assess model robustness. We finally discuss the main findings in light of implemented mechanisms and propose future development to improve the representation of the soils and drought stress in TBMs.

2 Material and methods

2.1 Study region

This study focuses on the South American tropical region, ranging from 90 to 30^∘ W in longitude and from 15^∘ N to 20^∘ S in latitude. The spatial resolution of the simulations was set to 1^∘. To drive the model, we used the 6-hourly CRU-NCEP v7 meteorological forcing dataset (Viovy, 2018). Climate variables include air temperature and humidity, incoming short-wave and long-wave radiation, precipitation rate, surface pressure, and winds. No land-use changes were applied to any simulation but a land-cover mask representative of the current plant functional type (PFT) distribution as derived from the ESA-CCI land-cover map (Poulter et al., 2015) is applied by default to ORCHIDEE simulations. This land-cover map corresponds to the year 2015, which is consistent with remote-sensing products used for model evaluation described below (see Sect. 2.5).

2.2 Vegetation models

In this study, we explored the sensitivity to soil texture of important model outputs (e.g. GPP, soil carbon) of three state-of-the-art TBMs with different levels of complexity. These TBMs, namely LPJ-GUESS, ED2, and ORCHIDEE v2.2, are briefly described in the next three subsections while a more detailed list of parameters, pedotransfer functions, and description of the impact of drought stress on plant productivity for each model can be found in Sect. S1 in the Supplement.

2.2.1 LPJ-GUESS

The Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS) model is a process-based dynamic vegetation model which can simulate the global vegetation distribution with its associated carbon, nitrogen, and water cycles (Smith et al., 2001, 2014; Oberpriller et al., 2022). The model has three possible modes of representing vegetation. The population mode is inherited from the LPJ model (Sitch et al., 2003), while the individual and cohort modes correspond to the vegetation representation of the GUESS model (Smith et al., 2001). For this study, the model was run in cohort mode (default vegetation representation). Cohorts represent the properties of the average individuals belonging to an age class of a given PFT. However, for herbaceous PFTs, the LPJ-GUESS model simulates only one average individual per patch. The coarsest spatial level in this model is the grid cell, for which soil texture, meteorological drivers, and nitrogen deposition should be provided. Different stands will each occupy a fraction of a given grid cell, representing different land cover and management types (natural vegetation, cropland, managed forests, etc.). Each stand contains one (population mode) or multiple (cohort and individual modes) replicate patches. The latter allows the model to account for heterogeneity in age distribution of the vegetation, due to stochastic differences in population dynamics. Within each patch, the different cohorts will grow and compete for light, water, and soil nitrogen.

Soil hydrology is represented by a multi-layer bucket model, where water can percolate between the different soil layers and drains at the bottom (Gerten et al., 2004). Soil depth is hard-coded to 1.5 m and subdivided into 15 layers of 10 cm thickness each. Soil moisture in the top two layers (20 cm) is available for surface evaporation. Yet, only 2 large layers are defined for percolation: excess water from the top layers (down to 50 cm) percolates into the bottom layers (remaining 100 cm), where it is distributed between the 10 layers depending on their water capacity. Soil hydraulic properties are derived from pedotransfer functions that require sand and clay contents for each grid cell (Cosby et al., 1984; Prentice et al., 1992; Haxeltine and Prentice, 1996). These are assumed to remain constant over the complete soil column. Soil water content is given as a fraction (0–1) of the available water capacity, which is in turn defined as the difference between the volumetric water content at field capacity and wilting point. Plant drought stress is expressed by the ratio between water supply and atmospheric water demand. If water supply is smaller than water demand, the PFT will be drought-stressed and canopy conductance will be reduced. Water supply is calculated as the product of a PFT-specific daily maximum transpiration rate (e_max), daily maximum root water uptake, and a factor which represents the leaf-phenological status as a fraction of the potential leaf cover. Daily maximum root water uptake is given as a function of the fractional root distribution and plant-available water content, summed over all soil layers. In the standard LPJ-GUESS parameterisation, this function is simply the product of both factors, further scaled by the total foliar projective cover in order to account for spatial overlap between cohorts.

2.2.2 ED2

The Ecosystem Demography model, version 2 (ED2) is a cohort-based vegetation model that simulates the energy, water, and carbon cycles of terrestrial ecosystems while accounting for their horizontal and vertical heterogeneities (Medvigy et al., 2009). The model was designed to be compatible with multiple configurations: it can be run as a stand-alone TBM over a single location, over a regional grid, or coupled with an atmospheric model distributed regionally (Knox et al., 2015). The coarsest hierarchical level of ED2 is the polygon within which time-varying meteorological forcing above the canopy is assumed uniform. Each polygon is subdivided into one or multiple sites with the aim of representing landscape-scale variations in abiotic properties like soil texture. Within the simulated sites, the horizontal heterogeneities in the ecosystem are simulated through a set of patches that represent the aggregation of all areas with a similar disturbance history. Finally, in each patch, the plant community population is tracked as a collection of plant cohorts, defined by their functional type and size. The ED2 has a typical time step of 10 min for the energy and water fluxes but can simulate succession and demography over larger (i.e. century) timescales.

In ED2, plant water availability is determined through a physically based soil-hydrology submodel, which encompasses heat, enthalpy, and water fluxes between different soil layers and the potentially existing temporary surface water. Water flux between soil layers is based on Darcy's law (Darcy, 1856; Bonan, 2008), surface runoff of water is simulated using a simple extinction function while subsurface drainage depends on the bottom boundary condition (e.g. free drainage, zero-flow, saturated water table). In ED2, soil depth, the number of soil layers, and layer thickness can be prescribed by the user but in the tropics, the soil is typically discretised into 16 layers along a 8 m depth soil profile with increasing layer thickness from top to bottom (Longo et al., 2019a). Most of the soil-hydraulic properties in ED2 are derived from the LEAF-3 model (Walko et al., 2000) and follow the parameterisation by Cosby et al. (1984) which is based on the soil volumetric fraction of sand and clay. Soil water retention and hydraulic conductivity curves are respectively based on Clapp and Hornberger (1978) and Brooks and Corey (1964), corrected for partially or completely frozen soil water. Drought stress negatively impacts plant productivity through a non-linear, soil-dependent wilting function, based on the ratio of water demand (plant transpiration) and supply (root water uptake). The latter is proportional to the soil water field capacity minus the soil water at the wilting point, integrated from the deepest soil layer accessible by plant roots to the soil surface. Rooting depth is related to plant height through an allometric relationship and root biomass is distributed over the soil layers according to their relative thickness.

2.2.3 ORCHIDEE v2.2

The process-based gridded vegetation model ORCHIDEE (ORganizing Carbon and Hydrology in Dynamic EcosystEms) is designed to simulate the fluxes of matter and energy, as well as the vegetation dynamics at the regional level (Krinner et al., 2005). ORCHIDEE v2.2 is the land component of the IPSL (Institut Pierre Simon Laplace) climate model developed for the Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al., 2016; Peylin et al., 2022). For a given vegetation, soil type, and climatic conditions, the model simulates physiological processes of an average ecosystem on a half-hourly time step, based on a combination of a dozen PFTs representing the major biomes on Earth.

Drought stress effect on vegetation is simulated through a physically based soil hydrology scheme and saturation-based Richards equation (Richards, 1931). The soil is discretised into 11 layers along a 2 m depth profile with increasing layer thickness from the top to the bottom (de Rosnay et al., 2002). Infiltration is processed before soil moisture redistribution while unsaturated values of hydraulic conductivity and diffusivity follow the models of Mualem (1976) and van Genuchten (1980). Soil parameters are set constant for each dominant USDA soil-texture class (Carsel and Parrish, 1988) provided as input. In ORCHIDEE, the soil texture is uniform over the soil column and only the saturated hydraulic conductivity decreases exponentially with depth to account for soil compaction and bioturbation (d'Orgeval et al., 2008). To compute infiltration and surface runoff, the model also accounts for horizontal variations in soil hydraulic conductivity (Vereecken et al., 2019). Soil evaporation and transpiration depend on soil moisture and properties, and transpiration is limited by stomatal resistance which increases when soil moisture drops from field capacity to wilting point. For each PFT, the root density decreases exponentially with depth up to 2 m, thus influencing the drought stress factor on transpiration. The drought stress factor will impact stomatal and mesophyll conductance at the leaf level and hence carbon assimilation. Finally, assimilated carbon is dynamically allocated to several vegetation pools, including the leaves. This will directly influence the seasonality in leaf area index (LAI), which has a feedback on the partitioning between soil evaporation and transpiration, and thus the resulting soil moisture.

2.3 Soil scenarios and simulation protocol

We performed three regional simulations with each model, using different soil-texture maps. These three soil maps represent the soil texture corresponding to the average (Mean clay), minimum (Min. clay), and maximum (Max. clay) topsoil clay content of each $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell from SoilGrids250m (Poggio et al., 2021), see Fig. 1. Soil texture and hence hydraulic properties were assumed to be vertically uniform for each model and simulation. For each model, both the soil depth and the number of soil layers were set up according to the most default model configurations for the tropics (Table 1).

Figure 1 Topsoil (0–5 cm) clay fraction spatial distribution as defined by the latest version of SoilGrids with a magnified example of a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell (a). The Min., Mean, and Max. clay scenarios are those soil types that are characterised respectively by the minimum, average, and maximum clay fraction content in each grid cell, excluding those without soil-textural information (for instance the rivers as illustrated in the magnified map). In panel (a), the density plot reveals the clay fraction distribution at the regional and local (magnified grid cell) levels. Subplot (b) shows both the intra-grid cell (i.e. mean of the standard deviation (SD) of the clay content) and the inter-grid cell (the standard deviation (SD) of the mean clay content) variability as a function of the spatial resolution. Subplot (c) is the resulting soil-texture distribution for each scenario. The arrows in subplot (c) show the average change of soil texture moving from the Min. clay to Mean clay scenario (black) and observed (red) changes in soil texture after deforestation (Eleftheriadis et al., 2018).

Table 1 Summary of the representation of soil, roots, and soil drought stress for each model in their default version.

For each scenario, long-term spin-up with pre-industrial atmospheric CO₂ concentration from 1860 (287.14 ppm) were performed to each grid cell starting from near-bare ground conditions. This was achieved by recycling the initial 10 years of CRU-NCEP until soil and vegetation carbon pools reached an equilibrium. The spin-up runs were continued with historical simulations from 1860 to 2016 for each model using the full CRU-NCEP forcing dataset and varying atmospheric CO₂ concentration according to Friedlingstein et al. (2020).

We did not change any model parameters (with the exception of the soil textural information) compared to the model default parameterisation for the tropics. More detail on each model's parameter sets can be found in the aforementioned references describing the three models. In ED2, we simulated four competing PFTs (grass, early-, mid-, and late-successional tropical trees) similarly to Longo et al. (2019a) for Amazon regional runs. To facilitate inter-model comparisons, we chose to run multiple ED2 simulations for each polygon rather than simulating multiple sites per polygon. For the LPJ-GUESS model, we activated all PFTs as for global simulations (Sitch et al., 2003; Ahlström et al., 2012), but due to bioclimatic limits only tropical broadleaf evergreen and raingreen trees as well as C₄ grasses emerged from the simulation with nonzero biomass.

2.4 Evaluation datasets

We used three different datasets to assess the model robustness and performance under the three soil scenarios. Firstly, we compared the spatial distribution of aboveground biomass (AGB) as simulated by the models with the integrated biomass map of Avitabile et al. (2016) over the study region. This biomass map is one of the reference products used by the International Model Benchmarking system (ILAMB; Collier et al., 2018) to evaluate TBMs, e.g. in the global carbon budget exercise (Friedlingstein et al., 2020). Secondly, we contrasted the model outputs of the average ecosystem gross primary production (GPP) with a moderate-resolution dataset of vegetation GPP derived from MODIS satellite data (MOD17A2), see Running et al. (2015). Finally, we differentiated the soil organic carbon stocks produced by each vegetation model/scenario with those derived from local observations upscaled to the globe in SoilGrids. We used those datasets for qualitative comparison only and not with the objective to improve model accuracy and/or precision.

2.5 Analyses

To assess the relative importance of the intra-grid cell variability, we aggregated the topsoil clay fraction in SoilGrids from its finest (250 m) to a much coarser (5^∘) resolution with the R package “raster”. For each resolution, we compared the intra-grid cell variability with the inter-grid cell variability. We defined the intra-grid cell variability as the average of the clay fraction standard deviation (SD) within each grid cell, and the inter-grid cell variability as grid-level SD of the average clay fraction of each grid cell.

All the results from the vegetation model simulations presented below are either (i) averages of the last 10 years of either the spin-up or the historical period (2006–2016) or (ii) the averages of the very last year of the historical period (2016). We particularly focused on the inter-model and inter-scenario comparison of GPP at the ecosystem and the PFT levels, as well as the resulting ecosystem-level AGB and soil carbon. For each simulation of each model, we also computed the normalised soil drought stress index (SDI) from the model outputs whose definition is model-specific (see Sect. S1) but always normalised between 0 (full stress) and 1 (no stress). We related SDI to the ecosystem GPP through quantile regression analyses using the R package “quantreg”. We used a quantile regression given the nature of the vegetation productivity response to SDI: for SDI close to 1, the GPP variability is high (other resources can limit GPP) while for SDI to 0, the GPP is necessarily low. To evaluate the model performance, we averaged the model outputs (GPP, AGB, soil carbon) for each grid cell and vegetation model/scenario over the period of observation of the remote-sensing products (see previous section). We compared the resulting maps through a correlation analysis to compare their spatial distribution, and compared their density distribution through standard metrics (mean, root mean squared error).

The soil textures across the 3 scenarios were also classified according to the 12 major soil-texture categories defined by the United States Department of Agriculture (Soil Survey Manual, 2002), using the “soiltexture” R package (Moeys, 2018), which allowed us to quantify soil class frequencies for each soil scenario and to define transition matrix when switching from 1 soil-textural map to another. All analyses and plots were performed in R version 3.6.3.

3 Results

3.1 Intra- and inter-grid cell variability in topsoil clay content

With the native spatial resolution of the SoilGrids product (250 m), we observed a wide distribution in clay content over the South American tropics ranging from nearly 0 % to 74 %, with a median around 28 % and a standard deviation of 7 % (Fig. 1a). Such extreme clay content values can also be found within grid cells when using the spatial resolution typically applied in TBMs: in the magnified 1^∘ grid cell (∼111 km at the Equator) of Fig. 1a, clay fraction varied between 0 % and 58 % with a median and a standard deviation of 23 % and 7 %, respectively. When SoilGrids was aggregated from its finest to coarser spatial resolutions, we observed a rapidly increasing intra-grid cell variability in soil texture: over the whole region, the average intra-grid cell variability (i.e. the mean standard deviation) in clay content strongly raised by 4 % from 250 m to 1^∘, and kept increasing up to 6 % at 5^∘. At 1^∘ resolution and coarser, the variability within and between grid cells reached similar orders of magnitude (Fig. 1b).

A clear shift from sandy/silty soils toward more clayey soils can be observed when moving from the Min. Clay to the Max. Clay scenario (Figs. 1c and S1A in the Supplement). The mean clay fraction reached 17 %, 28 %, and 34 % for the Min. clay, Mean clay, and the Max. clay scenario, respectively (Fig. 1c). The resulting changes in sand fraction density distributions were less marked, except for the Min. clay scenario (Fig. S1B).

3.2 Models performance in default configuration

In their most default configuration, all three models showed poor performances in capturing the spatial variability in AGB (Figs. 2 and S2 in the Supplement), GPP (Fig. 3), and soil carbon content (Figs. S3 and S4 in the Supplement) as estimated from independent products, regardless of the soil scenarios.

Figure 2 Aboveground biomass (AGB) spatial distribution, as generated by Avitabile et al. (2016) (a) or as predicted at the end of the historical period (average over the 2006–2016 period) by the three terrestrial biosphere models used in this study for the Mean clay scenario (b). The upper-right corners in each plot represent the AGB density distributions over the simulated region for all three scenarios (coloured lines) and the observations (black). Note that the land cover was prescribed in the ORCHIDEE model, while it was an emergent property of the ED2 and LPJ-GUESS models.

Figure 3 Ecosystem GPP as a function of the soil drought stress index (SDI) as predicted at the end of the historical period (2016) by each terrestrial biosphere model used in this study. The SDI values range between 0 and 1, with no stress represented by SDI=1, full stress conditions represented by SDI=0. The boxplots represent the distributions of the stress index for each scenario and the coloured lines are the 95 % quantile regression per scenario (with the same colour legend). Each dot is a grid cell (1^∘ resolution).

The reference AGB map from Avitabile et al. (2016) shows a bimodal distribution in biomass over the South American tropics as a reflection of the distribution in forest (12.5 ± 2.7 kg_C m⁻²) versus non-forest biomes (2.4 ± 2.3 kg_C m⁻², see Fig. 2). Both the ORCHIDEE v2.2 and ED2 models not only reproduced this bimodal distribution (non-forest peak at 3.1 ± 2.3 and 0.7 ± 1.5 kg_C m⁻²; forest peak at 12.3 ± 1.8 and 17.5 ± 2.4 kg_C m⁻², respectively), but also overestimated the overall AGB on average (8.0 and 11.4 kg_C m⁻² for both models, respectively while data average is 6.7 kg_C m⁻²). On the contrary, LPJ-GUESS simulated a unimodal biomass distribution with an overestimated average biomass of 12.0 ± 5.5 kg_C m⁻². The spatial correlation of AGB between the map of Avitabile et al. (2016) and the models varied between 0.35 (LPJ-GUESS) and 0.82 (ORCHIDEE), with an intermediate performance for ED2 (0.67).

When compared to remote-sensing estimates over the 2006–2016 period, ORCHIDEE v2.2 and ED2 overestimated the GPP with simulated average values of 2.4 ± 1.0 and 3.3 ± 1.5 ${\mathrm{kg}}_{\mathrm{C}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, respectively, compared to reference values of 2.2 ± 0.3 ${\mathrm{kg}}_{\mathrm{C}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. Only LPJ-GUESS simulated similar average values of 2.2 ± 0.6 ${\mathrm{kg}}_{\mathrm{C}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$.

Finally, all three models overestimated soil carbon content compared to the information from the SoilGrid database with LPJ-GUESS, ED2, and ORCHIDEE v2.2, respectively, simulating mean values of 14.8 ± 6.3, 10.6 ± 8.2, and 7.6 ± 2.9 kg_C m⁻² for a reference value of 4.3 ± 1.7 kg_C m⁻². Also the simulated spatial distributions of the soil carbon content were drastically different from the reference one (Figs. S3 and S4), regardless of the soil scenario. The spatial correlations of soil carbon between SoilGrids250m and the models were systematically low: 0.10 for ORCHIDEE, 0.25 for ED2, and 0.26 for LPJ-GUESS.

The relative better performance of the ORCHIDEE model in capturing the spatial variability in vegetation and average soil carbon stocks can be partly explained by the use of a land-cover map to constrain vegetation type distribution compared to ED2 and LPJ-GUESS for which the PFT distribution is an emergent property of the models.

3.3 Model sensitivity to clay content variability

Large differences in all investigated model outputs existed between models for the same soil scenario. However, the performance of each model was almost independent of the soil scenario for all investigated products (AGB, GPP, soil carbon). All three models exhibited a strong correlation between the soil drought stress index (SDI) and the overall ecosystem productivity, as illustrated by the quantile regression analysis (Fig. 3). Across the three scenarios, we observed that increasing clay content slightly increased drought stress (i.e. decreased SDI) by 2.6 %, 0.7 %, and 1.5 % (change of the drought stress index from the Min. clay to the Max. clay scenario) for ORCHIDEE, ED2, and LPJ-GUESS, respectively (horizontal boxplots in Fig. 3). This increase in simulated drought stress was accompanied by a decrease in productivity for all three models, respectively by 2.7 %, 1.9 % and 3.2 % (vertical boxplots in Fig. 3).

Nonetheless, we observed substantial changes in PFT-level GPP simulated for the three scenarios for some grid cells (Fig. 4), especially for the ORCHIDEE model, which indicates some shifts in the simulated PFT composition. This situation occurred in about 3 % of the grid cells for ED2, 6 % of the grid cells for LPJ-GUESS, and 7 % for ORCHIDEE when switching from the Mean clay scenario to the Min. clay or the Max. clay scenario (Fig. S5 in the Supplement). Yet, these PFT-level shifts in GPP compensate for each other when aggregated at the ecosystem level, resulting in similar total GPP and spatial distributions that remain almost unaffected by shifts in soil composition (Fig. 3). Not only did the simulated SDI and GPP not change substantially, but we also observed very limited shifts in soil carbon content (Figs. S3 and S4) and AGB (Figs. 2 and S2) in response to changes in soil clay content (Fig. 2). Between the Min. clay and Max. clay scenarios, we observed a 3.0 %, 0.7 %, and 4.2 % increase in the average simulated AGB, and a −11.9 %, 10.1 %, and 7.6 % change in soil carbon content, as simulated by ORCHIDEE, ED2, and LPJ-GUESS. All the aforementioned observations also apply to state conditions resulting from the spin-up phase (as exemplified for the AGB spatial distribution at the end of the spin-up for all three models and all three scenarios, see Fig. S6 in the Supplement).

Figure 4 Scenario (Max. clay: a–c and Min. clay: d–f) versus reference (Mean clay scenario) GPP for each TBM used in this study. Each dot is the PFT-level GPP over a specific grid cell (1^∘ resolution) at the end of the historical period (2016).

We observed some substantial impacts of the scenario on the ecosystem GPP (up to a 100 % change of ecosystem GPP) for some of the soil-textural class transitions (which represent the frequency of soil class changes when moving from one soil-textural map to another), but those transitions were rather rare events and hence limited to a small area of the simulated region (Fig. 5). The most frequent transitions were within the same soil class (the diagonal of the soil transition matrix of Fig. 5a): those represented 32 % of all transitions between the Mean clay and the Min. clay scenarios, and 43 % of all transitions between the Mean clay and the Max. clay scenarios and were almost unaffected by the soil clay content (relative change of ecosystem GPP between −0.6 % and 2.1 % for all models and scenarios). For both ED2 and ORCHIDEEv2.2, the most important changes occurred over (very) low water-availability regions (Mean Annual Precipitation or MAP<1000 mm for ED2, MAP<2000 mm for ORCHIDEE) while the sensitivity to soil texture was independent of the water availability in LPJ-GUESS (Fig. 6). Yet, only a small fraction of the water-limited area was concerned by changes of simulated state variables: 20 % of the grid cells with MAP<1000 showed a relative change of the ecosystem GPP larger than 10 % in both ED2 and ORCHIDEEv2.2. Moreover, positive and negative shifts balanced one another and hence had a very limited impact on the regional ecosystem productivity (Figs. 3 and 6). This was partly due to a lack of sensitivity of the pedotransfer functions to the explored range of variations of soil texture (Fig. S7 in the Supplement), leading to narrower soil hydraulic parameter distributions compared to previously generated at this spatial scale (Montzka et al., 2017).

Figure 5 Soil transition matrix representing the frequencies of moving from one soil class to another when changing the soil-textural map (a) and the relative change of ecosystem gross primary production (GPP) for each category of transition and terrestrial biosphere model (b) as predicted by each TBM at the end of the historical period. In (a), the colour intensity represents the frequency of each transition. The grey cells are transitions that did not occur in the simulated scenarios.

Figure 6 Relative change of the annual ecosystem GPP at the end of the historical simulation with the mean annual precipitation (MAP) for both scenarios (shapes) and all three TBMs considered in this study, across the entire simulated region (each point is a grid cell). The MAP is the annual average over the last 10 years (2006–2016) of the CRU-NCEP dataset forcing. R² of all linear models were lower than 0.01.

Finally, we note that the LPJ-GUESS model crashed for some specific soil textures. Those soil textures occur naturally in the field and were relatively frequent in our simulations, especially in the Min. clay scenario (5.4 % of the grid cells). The problem occurred for silty soils with low fractions of both sand (<12 %) and clay (25 %). The default pedotransfer functions applied to those specific soil textures led to volumetric water content at field capacity larger than the water content at saturation (see Fig. S8 in the Supplement), which caused the model to crash during model initialisation for those particular grid cells.

4 Discussion

The South American tropics frequently suffer from (mega)droughts such as the 2015–2016 El Niño event, which severely impacted ecosystems on the continent. These droughts result from the increasing environmental pressures, including climate change and deforestation (Staal et al., 2020). In a recent study, Yang et al. (2022) linked patterns of forest biomass changes with drought severity and duration as well as soil clay content, which indicates that both water demand and supply (or both climate and soil/roots) influence forest functioning. In contrast to that study, we found here that the aboveground biomass simulated by three state-of-the-art TBMs was mostly insensitive to soil texture, except for some limited areas with low water availability in ED2 and ORCHIDEEv2.2.

In TBMs, plants and soils are coupled by a drought stress function which depends on soil moisture. We suggest that the overall lack of sensitivity that we observed originates from the combination of two main limitations in the current implementation of the hydrology submodel in TBMs: (1) the shifts in soil texture resulting from the spatial variability in soil clay content does not translate into realistic shifts in soil hydraulic properties (Fig. S7) and (2) the implemented drought stress functions do not properly capture the effect of changes in soil hydraulic properties on vegetation. Both limitations should be rapidly tackled in order to improve TBM performance (Fisher and Koven, 2020) and are briefly discussed below.

Current TBMs use a limited number of generic, widespread pedotransfer functions, which can be class-based or continuous. However, most of these functions were developed and calibrated decades ago (1984 and 1988, see Table 1) with fewer and less geographically spread calibration data than what are available today. On top of the limited size of the training data (especially for the tropical regions), the main drawback of these pedotransfer functions resides in their inability to capture the variability and non-linearity of many parameters for given soil-textural classes. For example, in their review, Van Looy et al. (2017) highlighted large differences in saturated hydraulic conductivity within each soil class derived from different data sources and location. As a result, by using generic and global functions, soil parameters in TBMs are substantially different from region-specific observations (Kishné et al., 2017; Van Looy et al., 2017). Such generic, global functions could therefore lead to inaccurate characterisation of the soil properties, as illustrated by the LPJ-GUESS crashes with realistic soil compositions (Fig. S8). Since no generic functions are able to properly capture soil properties at the global scale (Patil and Singh, 2016), intermediate solutions should be implemented in TBMs for a better representation and scaling of soil properties. For example, region-specific pedotransfer functions, regional calibration, ensemble simulations using multiple pedotransfer models, or the combination of regional pedotransfer functions could be used to estimate the uncertainties that soil properties are responsible for (Hodnett and Tomasella, 2002; Barros and de Jong van Lier, 2014; Medeiros et al., 2014).

Soils have a direct, strong role in the response of plants to drought (Carminati and Javaux, 2020). All three vegetation models used in our study apply simple drought stress functions that depend on the available water and the water demand (Table 1 and Sect. S1). We observed that shifting soil properties from low to high clay content barely affected the simulated soil drought stress despite substantial changes in soil texture and classification (Figs. 1, 3 and S1). Even if the sensitivity to soil texture might slightly increase with drought stress (Fig. 6) and hence under future climate change scenarios, it confirms that generic drought stress functions are not suitable to capture the impact of changes in plant water availability on plant processes, as suggested from previous studies (Uribe et al., 2021; Joetzjer et al., 2014; Combe et al., 2016) and argues for a better representation of root–soil coupling in TBMs. To reproduce the annual pattern of net ecosystem exchange of carbon over the Amazon with the Simple Biosphere Model (SiB3), Baker et al. (2008) demonstrated the importance of combining multiple mechanisms, not only related to soil water distribution but also on root dynamic schemes. Indeed, the Amazon forest was shown to have high GPP during the dry season (Green et al., 2020; Negrón Juárez et al., 2007; Saleska et al., 2003), and the role of water uptake with deep roots is currently not properly simulated in TBMs, e.g. because of the shallow soils (and hence root systems) that are simulated (Table 1), (Verbeeck et al., 2011; Nepstad et al., 1994). Recent developments in the TBM community have focused on improved plant hydraulics but to a lesser extent to the root–soil interactions (Xu et al., 2016). Recent studies have demonstrated the need for a better representation of root water uptake in drying soils to simulate plant response to drought stress and its impact on biomass. For instance, the new dynamic root scheme (Joetzjer et al., 2022) coupled to explicit plant hydraulic processes in ORCHIDEE managed to reproduce observed water and carbon dynamics at the Caxiuanã throughfall exclusion field experiment in eastern Amazonia (Yao et al., 2021). Although they better capture biomass and flux dynamics at the site level, the new implementations of plant hydraulics is empirical and complex, and leads to an increased number of model parameters and hence to a larger required size of field observational data to calibrate it. Such data are not readily available for a large number of sites or a larger area and/or might be difficult to measure, especially in complex ecosystems like tropical forests. The process of model complexification might also result in over-fitted simulations at the site level, mainly focusing on climate factors (e.g. drought) while overlooking the unconstrained soil and root properties. Better estimate of water demands based on eco-evolutive optimality theories (Prentice et al., 2014), root biomass (Franklin et al., 2012), and soil–root interactions (Lu et al., 2020; Vanderborght et al., 2021) could help bridge the gap between complex and over-parameterised models on the one hand and simple unrealistic model-specific functions on the other.

Soil texture and clay content have a direct, strong impact on the distribution and mineralisation of carbon and nutrients in soils (e.g. Hassink, 1992; Telles et al., 2003; Plante et al., 2006; Zinn et al., 2007). In our simulations, the three TBMs could not reproduce SoilGrid soil carbon distribution and showed very low sensitivity to changes in clay content (Figs. S3 and S4) despite long-term spin-up during which we expected large differences between equilibrium states induced by different soil composition. This further highlights the poor representation of soil processes in TBMs and their coupling to vegetation dynamics. As for the relationship between soil hydraulic properties and texture, we argue that current development focusing on soil processes in TBMs (e.g. nutrient mineralisation, soil organic carbon) should systematically assess model sensitivity to soil properties and texture parameters.

By selecting the dominant soil-texture class, a significant fraction of the soil spatial heterogeneity is omitted in TBMs running at coarse spatial resolution. This effect has already been documented in Tafasca et al. (2020), who suggested that spatial aggregation statistically enhances medium textures, leading to excessive evapotranspiration and insufficient total runoff. Accounting for subgrid variability in soil texture and moisture through systematic sensitivity analysis, or directly representing this effect in TBMs with models (Qu et al., 2015) could alleviate these uncertainties and improve model performance. Intra-grid cell variability in soil texture might have large impacts on simulating vegetation dynamics, especially in demographic models for which plant competition and access to resources drive ecosystem composition and dynamics (i.e. growth/mortality) (Rowland et al., 2015; Johnson et al., 2016). In addition to the aggregation bias, we also expect substantial biases in simulated ecosystem properties resulting from intrinsic uncertainties from soil products. In SoilGrids, maps of soil properties are generated using machine-learning methods that account for direct soil observations and environmental variables describing vegetation, climate, topography, geology, and hydrology. However, the number of soil observations available over the Amazon tropics is very low (10 or fewer soil-textural observations for grid cell of 70 000 km² for many grid cells in the studied area), potentially leading to high uncertainties in regional soil properties at fine resolution (see Fig. S6 from Poggio et al., 2021).

Realistically reproducing the South American tropics in silico suppose that TBMs accurately represent the most essential processes driving ecosystem functioning. In the context of increasing drought intensity and severity, a better representation of the impacts of soils on plant productivity and status is urgently needed. More generally, the belowground compartment should receive more attention from the vegetation modelling community. It has been long known that root and soil depth critically influence drought tolerance (Nepstad et al., 1994; Fan et al., 2017), yet this knowledge has not been integrated in most TBMs (Verbeeck et al., 2011). In this study, we focused on soil texture only, but further research should include other soil properties (e.g. soil structure see Fatichi et al., 2020) and root traits. Other critical processes are impacting this region and require as much attention. The Amazon basin has gone through intense deforestation activities for more than 3 decades. Among the consequences of such environmental pressure, deforestation increases drought (Staal et al., 2020) and modifies soil properties (Veldkamp et al., 2020). Interestingly, the impact of deforestation on soil texture is of the same order as the one investigated in this study (Fig. 1c). These processes and their interactions should be accounted for in the next generation of TBMs.

5 Conclusion

The TBMs are keystones of global carbon and water budget assessments. Past developments strongly focused on representing plant processes and their response to climate. Despite their importance and recent efforts from the TBM community, belowground processes have remained overlooked. Here, we showed that the carbon-related processes are mostly insensitive to soil texture over the South American tropics for all three investigated TBMs (LPJ-GUESS, ORCHIDEE, and ED2). These results suggest a poor representation of the soil–vegetation coupling in TBMs, mainly because of inadequate pedotransfer functions and soil drought stress definitions. To date, the use of generic pedotransfer and drought stress functions is common in the TBMs which are used for carbon and water budget assessments, as well as future projections, which leads to large errors in the model predictions. Appropriately representing soil spatial heterogeneity (through better estimates of the impact of soil-texture uncertainties) and soil–plant coupling, such as the non-linearity of soil–root resistance, is a major challenge that needs to be urgently addressed in TBMs. This will lead to a better representation of the effect of drought stress on vegetation and to a reduction of the carbon budget uncertainties, which is particularly needed in complex and heterogeneous ecosystems such as tropical forests.