Spatially varying parameters improve carbon cycle modeling in the Amazon rainforest with ORCHIDEE r8849

Uncertainty in the dynamics of the Amazon rainforest poses a critical challenge for accurately modeling the global carbon cycle. Current dynamic global vegetation models (DGVMs), which use one or two plant functional types for tropical rainforests, fail to capture observed biomass and mortality gradients in this region, raising concerns about their ability to predict forest responses to global change drivers. Here we assess the importance of spatially varying parameters to resolve ecosystem spatial heterogeneity in the ORCHIDEE (ORganizing Carbon and Hydrology in Dynamic EcosystEms) DGVM. Using satellite observations of tree aboveground biomass (AGB), gross primary productivity (GPP), and biomass mortality rates, we optimized two key parameters: the alpha self-thinning (α), which controls tree mortality induced by light competition, and the nitrogen use efficiency of photosynthesis (η), which regulates GPP. The model incorporating spatially optimized α and η parameters successfully reproduces the spatial variability of AGB (R² = 0.82), GPP (R² = 0.79), and biomass mortality rates (R² = 0.73) when compared to remote sensing observations in intact Amazon rainforests, whereas the model using spatially constant parameters has R² values lower than 0.04 for all observations. Furthermore, the relationships between the optimized parameters and ecosystem traits, as well as climate variables, were evaluated using random forest regression. We found that wood density emerges as the most important determinant of α, which is in line with existing theory, while water deficit conditions significantly impact η. This study presents an efficient and accurate approach to enhancing the simulation of Amazonian carbon pools and fluxes in DGVMs by assimilating existing observational data, offering valuable insights for future model development and parameterization.

1 Introduction

More than half of the global rainforests are located in the Amazon basin, storing approximately 140 PgC in their living biomass (Pan et al., 2011; Vancutsem et al., 2021). They constitute an important carbon sink in intact forests, but deforestation and forest degradation can turn these areas into a carbon source (Brienen et al., 2015; Gatti et al., 2021). The Amazon rainforest is characterized by high levels of biodiversity, rainfall variability, and soil diversity (Malhi et al., 2004; ter Steege et al., 2006; Flores et al., 2010; Quesada et al., 2010; Castanho et al., 2013). These factors are only partly resolved in dynamic global vegetation models (DGVMs), which are widely used to predict vegetation responses to global environmental changes (Prentice et al., 2007). Therefore, current DGVMs can hardly capture continental-scale carbon fluxes and stock spatial gradients (Johnson et al., 2016), leading to even higher uncertainty in projecting future carbon dynamics in the Amazon rainforests.

The factors driving the spatial variability of carbon fluxes and stocks in the Amazon rainforests are complex and remain largely unresolved (Johnson et al., 2016; Muller-Landau et al., 2021). Field forest plots reveal a clear aboveground biomass (AGB) gradient: higher AGB in the northeast over the Guiana Shield, where nitrogen-fixing, large-seeded legume trees with lower growth rates are dominant, and lower AGB in the drier southwest, where non-legume trees thrive (ter Steege et al., 2006; Malhi et al., 2006; Mitchard et al., 2014). The biomass gradient results from the balance between the spatial variation in woody productivity and carbon losses, mainly through tree mortality. Field observations suggest a strong association between the spatial variation of tree mortality and the gradients of AGB (Malhi et al., 2015; Johnson et al., 2016). In western and southern regions, where natural disturbances are more frequent, dominant species tend to have higher growth rates, lower wood density, and higher mortality rates, supporting the “grow fast, die young” life history strategy, and maintaining stable biomass (Baker et al., 2004; Malhi et al., 2006; Keeling and Phillips, 2007; Esquivel-Muelbert et al., 2020; Brienen et al., 2020). Woody net primary productivity (NPP) is influenced by both climate conditions and soil fertility, particularly total soil phosphorus content, which is higher in the western region. The mortality rate could be driven by the intensity of forest disturbances, particularly windthrow and drought. In the western region, windstorms are the primary cause of tree death, while in the east-central and southern regions, drought is the dominant driver (Esquivel-Muelbert et al., 2020). Soil physical properties including soil depth, soil structure, topography, and anoxic conditions also influence tree mortality (Malhi et al., 2004; Quesada et al., 2012; Sousa et al., 2022; Feng et al., 2023).

Current DGVMs cannot capture the spatial heterogeneity in productivity, biomass, and mortality rates in the Amazon. While most models include detailed photosynthetic modules, they often ignore nutrient constraints on plant growth (e.g., phosphorus; Reed et al., 2015), and they simplify varied life history strategies to a few plant functional types with static parameters (e.g., mortality rates; Johnson et al., 2016). Although an increasing number of models have incorporated phosphorus cycles (Goll et al., 2012; Wang et al., 2010), simplified tree demography (Moorcroft et al., 2001; Rödig et al., 2017; Koven et al., 2020; Naudts et al., 2015), and drought-induced mortality (Yao et al., 2022), their capability to capture these spatial heterogeneities remains limited. More importantly, most models generally use one or two plant functional types (PFTs) with uniform parameters across the entire Amazon forest. Building upon these setups, in DGVMs, variations in biomass carbon pools across different locations are primarily driven by differences in climate forcing and soil properties. However, biomass variability in the Amazon rainforests cannot be fully explained by climate factors alone, and soil properties are inadequately represented in DGVMs due to the coarse resolution and simplified processes (Quesada et al., 2012; Galbraith et al., 2013; Saatchi et al., 2015; Joetzjer et al., 2022). Given the availability of various data for parameter estimations, applying a single set of PFT-specific parameters uniformly across the entire basin is no longer justifiable (Butler et al., 2017).

Several studies have attempted to simulate the spatial variations of biomass and gross primary productivity (GPP) in the Amazon forest by linking different processes or forcing the model with observed parameters. In models that apply a constant tree mortality rate, higher NPP is typically associated with higher biomass (Keeling and Phillips, 2007). Delbart et al. (2010) demonstrated that linking higher tree mortality rates to higher NPP in a DGVM can improve the simulation of spatial biomass variations, but may fail to capture current biomass sinks in pristine tropical forests (Brienen et al., 2015). Castanho et al. (2013) forced another DGVM with an interpolated map of observed woody biomass residence times and correlated maximum carboxylation rate (V_cmax) with a soil total phosphorus map. While this study successfully reproduced site-level woody NPP and AGB, the use of interpolated woody biomass residence times and soil phosphorus data, based on limited observations, remains questionable (Saatchi et al., 2015; He et al., 2023).

The increasing availability of remote sensing data provides wall-to-wall observations of key forest properties such as AGB, GPP, tree height, leaf area index (LAI), and normalized difference vegetation index (NDVI), which all remain underutilized for parameterizing and evaluating models. Rödig et al. (2017) developed an individual-based model, using a remote sensing canopy height map as a constraint, and linked mortality rates to soil clay fraction and precipitation to simulate the spatial gradient of biomass. Ma et al. (2024) incorporated more observations from satellites, including GPP, LAI, AGB, and forest age, and applied a stepwise calibration at the global scale. Despite these advancements, few studies have simultaneously employed multiple remote sensing datasets to optimize their models in the Amazon forest.

In this study, we aim to improve the modeled spatial gradient of biomass, productivity, and mortality across the intact Amazon rainforest in the ORCHIDEE (ORganizing Carbon and Hydrology in Dynamic EcosystEms) model. This model simulates forest demography based on tree cohorts and incorporates tree mortality mechanisms, including self-thinning and drought-induced mortality (Yao et al., 2022). We assessed whether spatial variation in model parameters can improve the ability of the ORCHIDEE model to simultaneously reproduce multiple observed patterns, including AGB, GPP, and biomass mortality rates derived from remote sensing data. Finally, we discuss mechanistic links between environmental factors and spatial gradients in model parameters by exploring the results of an explainable machine-learning method, in order to guide future model developments.

2 Materials and methods

In this study, we describe the optimization of the model parameters for the Amazon forest (plant functional type: tropical broadleaf evergreen tree) using remote sensing data. We first set a baseline model with optimal, but spatially constant, parameters at the Amazon basin scale. Next, we allowed parameters to vary in space with a spatial resolution of 1° × 1°. Finally, we used a random forest model to derive relationships between the potential influencing factors and the optimal spatially varying parameters. This approach enabled us to explain the emerging spatial variations of the parameters across the Amazon basin and prescribe their spatial variability in the model (Fig. 1).

Figure 1 Flowchart for the parameter optimization. η: nitrogen use efficiency; α: alpha self-thinning; LHS: Latin hypercube sampling; GPP: gross primary productivity.

2.1 Study area

This study focuses on the Amazon basin (79–43° W, 18° S–10° N), predominantly covered by old-growth, intact tropical evergreen rainforest (Flores et al., 2024). We excluded regions based on the following criteria: (1) areas with mean annual precipitation below 1500 mm and mean annual temperature below 18 °C based on CRUJRA v2.4 from 1901–1920 (Climatic Research Unit and Japanese Reanalysis data; Harris et al., 2020) to keep the focus on the warm and humid tropics; (2) grid cells (1° × 1°) with less than 10 % forest cover, based on forest pixels that are labeled as “undisturbed tropical moist forest” in TMF (tropical moist forest; Vancutsem et al., 2021) to exclude the pixels that were likely to be regrowing or degraded.

2.2 The ORCHIDEE model

ORCHIDEE is a global terrestrial ecosystem model that simulates energy, water, and carbon cycles (Krinner et al., 2005; Naudts et al., 2015; Vuichard et al., 2019). In this study, we used ORCHIDEE r8240 and further incorporated soil–tree–leaves hydraulics and drought-induced tree mortality (Yao et al., 2022). We also adjusted the NPP allocation parameters and forced the model with a spatially explicit wood density map (See Text S1, Table S1 in the Supplement; Chave et al., 2010; Yang et al., 2024a; Chen et al., 2020; Yang et al., 2021). These developments collectively constitute the new ORCHIDEE version r8849.

The input data for the ORCHIDEE model include soil texture and meteorological forcing (air temperature, precipitation, shortwave downward radiation, longwave downward radiation, specific humidity, wind speed, and atmospheric pressure). We used CRUJRA v2.4 with a 6 h temporal resolution as the source of meteorological forcing (Harris et al., 2020). To reduce computational time, we resampled the original CRUJRA data from a resolution of 0.5° × 0.5° to 1° × 1°. Consequently, all results are presented at 1° × 1° resolution. For soil texture data we used the Harmonized World Soils Database version 2.0 (HWSD2; Fig. S1 in the Supplement; FAO and IIASA, 2023), a comprehensive global soil database that includes various national datasets at a resolution of 1 km. In the ORCHIDEE model, soil texture determines the fractions of clay, silt, and sand, with other soil properties derived from a lookup table (Table S2) that controls soil hydraulic conductivity and soil water holding capacity. Soil depth is prescribed as 2 m in the model. Given that more than 99 % forest in our study area is tropical broad-leaved evergreen trees (Harper et al., 2023), we forced the model with only the tropical broad-leaved evergreen trees PFT in each pixel. Here, we provide a summary of the model processes. All the changed parameters are summarized in Table S1.

2.2.1 GPP

In the ORCHIDEE model, GPP is calculated based on the leaf scale model of Farquhar et al. (1980), further simplified by Yin and Struik (2009). The key parameters driving the photosynthesis process are V_cmax (µmol CO₂ m⁻² s⁻¹) and J_max (the maximum electron transport rate, µmol CO₂ m⁻² s⁻¹), with J_max typically assumed to be proportional to V_cmax. In the ORCHIDEE model, V_cmax is dynamically calculated as follows:

$$\begin{array}{}\text{(1)}& V_{\mathrm{cmax}}={\displaystyle \frac{\mathit{\eta }\times M_{\mathrm{n},\mathrm{l}}^{\mathrm{active}}}{{\sum }_{i=\mathrm{1}}^n\left(\right(\mathrm{1}-\mathrm{0.7}\times (\mathrm{1}-{\mathrm{Light}}_i))\times {\mathrm{LAI}}_i)}},\end{array}$$

where $M_{\mathrm{n},\mathrm{l}}^{\mathrm{active}}$ (gN m⁻²) is the active nitrogen content in the leaves, Light_i is the fraction of cumulative light transmitted to canopy level i, and LAI_i (m² m⁻²) is the leaf area index at level i. The ORCHIDEE model utilizes a dynamic three-dimensional canopy structure model to simulate gap probability (Haverd et al., 2012). The crown is divided into n levels (n=10). At each level, the model calculates the leaf volume based on tree crown morphology and density. LAI_i is then calculated by leaf mass, specific leaf area (SLA = 0.0153 m² gC⁻¹), and leaf volume for this layer. The canopy gap is calculated based on LAI_i, tree height, tree density, and solar altitude angle. Light_i decreases progressively from 100 % at the top level to the gap fraction at the bottom level. Given the limited observations available to constrain the nitrogen cycle in the Amazon forest and prevailing limitation by geogenic soil nutrients like P or K, we instead prescribe leaf nitrogen concentration for the whole Amazon forest (Sitch et al., 2003). In addition, we deactivate nitrogen controls on ecosystem processes and associated soil–plant feedbacks (Text S2).

2.2.2 Biomass mortality rates

We defined biomass mortality rates as the ratio of biomass loss to the total biomass. The ORCHIDEE model incorporates 20 circumference tree size classes within each forest PFT (Naudts et al., 2015), thus facilitating the representation of annual recruitment of young individual trees, and tree mortality driven by light competition (i.e., self-thinning). Self-thinning-induced tree mortality is calculated as

$$\begin{array}{}\text{(2)}& N_{max}={\left({\displaystyle \frac{D_{\mathrm{g}}}{\mathit{\alpha }}}\right)}^{\frac{\mathrm{1}}{\mathit{\beta }}},\end{array}$$

where N_max is the maximum number of trees per hectare for a given quadratic mean tree diameter (D_g, cm), and α (cm) and β are the self-thinning parameters. The model calculates the relative density index (RDI) as the ratio of the number of trees to N_max. When the RDI exceeds an upper threshold (RDI_upper), the excess trees are removed through mortality evenly in each circumference class, until the RDI reaches a lower threshold (RDI_lower). Both RDI_upper and RDI_lower are cubic functions of D_g, with user-defined parameters. This self-thinning Eq. (2) indicates that the growth of some trees occurs at the expense of others, which die due to competition for limited resources such as water, light, or nutrients (Yu et al., 2024). The power-law function assumed in the model has been widely observed in the field plots (Enquist and Niklas, 2001; Muller-Landau et al., 2006). Higher values of α and β suggest greater environmental carrying capacity and lower biomass mortality rates. In the model, the default value for β is −0.73, which is close to the average value reported by Yu et al. (2024) for tropical forests. However, the default value for α has not been evaluated and, to our knowledge, no empirical estimates for α have been provided in previous studies. We thus estimate α based on field forest inventory data and find the average value is 1941 cm (Brienen et al., 2015).

The model also incorporates a drought-induced tree mortality module (Yao et al., 2022) to resolve periodic increases in mortality associated with droughts. When drought-induced mortality exceeds the mortality from the self-thinning process, the model assigns the total mortality rate to the drought-induced value; otherwise, it defaults to the self-thinning mortality rate. The drought mortality parameters were calibrated previously for the drought experiments in Caxiuanã National Forest (Yao et al., 2022). Although drought-induced mortality can reduce biomass periodically in addition to self-thinning, AGB and biomass mortality rates are predominantly controlled by α rather than drought-induced mortality at the basin scale during our study period (Text S1.3; Fig. S2). Considering that self-thinning is the dominant and regularly occurring mortality process, we treat self-thinning mortality as background mortality driven by resource competition among trees in this study.

2.2.3 Aboveground biomass

The ORCHIDEE model has nine biomass pools: aboveground sapwood (Sap_ab), belowground sapwood, aboveground heartwood (Heart_ab), belowground heartwood, leaves, roots, fruits, labile carbohydrate (Carbo_lab), and reserved carbohydrate (Carbo_res). The modeled aboveground biomass is calculated as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{AGB}& ={\mathrm{Sap}}_{\mathrm{ab}}+{\mathrm{Heart}}_{\mathrm{ab}}+\mathrm{Leaf}+\mathrm{Fruit}\\ & +{\displaystyle \frac{{\mathrm{Carbo}}_{\mathrm{lab}}+{\mathrm{Carbo}}_{\mathrm{res}}}{\mathrm{2}}}.\end{array}\end{array}$$

2.3 Observational data

We used GPP observations from the GOSIF dataset for data assimilation at the grid scale (Li and Xiao, 2019). GOSIF products estimated GPP using sun-induced chlorophyll fluorescence (SIF) from GOSAT, GOME-2, and OCO-2 at a 0.05° resolution from 2001 to 2022. Given the high uncertainty in GPP estimates, particularly in tropical forests where GPP observations are scarce, we also evaluated other GPP products derived from flux tower upscaling (FLUXCOM_RS and FLUXCOM_RS_METEO) and process-based models (BESSv2 and FORMIND) as part of our sensitivity tests (Jung et al., 2020; Li et al., 2023; Rödig et al., 2019). Across these datasets, the average GPP values for the Amazon basin range from 24.0 to 33.4 MgC ha⁻¹ yr⁻¹. Although GOSIF provides the highest estimation, it appears to be more consistent with observed GPP values (Fig. S3, Marthews et al., 2012; Malhi et al., 2015). We compared GPP products with model outputs from 2001 to 2015 for GOSIF, FLUXCOM_RS, FLUXCOM_RS_METEO, and BESSv2. For FORMIND, the GPP represents average values from 2003 to 2006.

We use the biomass map from Yu et al. (2023) to benchmark AGB in the ORCHIDEE model in 2020 (hereafter, Yu-Biomass). The Yu-Biomass map is based on a large number of lidar-biomass models, field plot data, and L-band radar data, with an original resolution of 100 m. Several other biomass maps are available for this region (Santoro and Cartus, 2024; Avitabile et al., 2016; Baccini et al., 2012). These maps generally incorporate field observation data and multiple remote sensing datasets, including radar backscatter, lidar, and optical indices. Among these maps, Yu-Biomass map shows the highest correlation with independent airborne lidar estimates in the Amazon forest (Longo et al., 2016; Yu et al., 2023). We also validated these datasets using field observations (Mitchard et al., 2014) and found that the Yu-Biomass map captures the biomass gradient more accurately than the others (Fig. S4).

We optimized the modeled biomass mortality rates using biomass mortality data derived from Planet data (Dalagnol et al., 2023, 2021). This dataset estimates forest gaps based on very-high-resolution (about 4–5 m) remote sensing and airborne lidar scanning data processed through deep learning methods. The identified gaps are then converted to biomass mortality estimates. The dataset provides annual biomass mortality data from 2016 to 2019 aggregated at a 1 km resolution. Given that this study focuses on background mortality, to minimize the potential influence of the 2015/16 El Niño drought, we excluded the year 2016 and used the average from 2017 to 2019. Biomass mortality rates were derived by dividing biomass mortality by the total aboveground biomass from Yu-Biomass. We also compared the remote sensing based biomass mortality rates to stem mortality rates from 189 long-term RAINFOR forest plots, with data primarily collected between 1981 and 2010 (Esquivel-Muelbert et al., 2020). This comparison yielded an R² of 0.39 and an RMSE of 0.53 % (Fig. S5).

2.4 Simulation protocol

We began with a spin-up simulation to bring the vegetation biomass carbon pool to equilibrium, i.e., a linear trend in vegetation carbon over 40 years of <0.03 MgC ha⁻¹ yr⁻¹ for 99 % of pixels. We recycled the climate forcing from CRUJRA from 1901–1920 for 600 years with a constant CO₂ concentration of 296.8 ppm (i.e., the CO₂ concentration in 1901). The equilibration of the soil biogeochemistry was calculated by an analytical solution (Vuichard et al., 2019). Starting from this equilibrium state, we proceeded with the transient simulation using climate forcing and CO₂ concentration data for the period 1901–2020. We did not account for land cover change, as our analysis is limited to undisturbed forests. Model parameters remain unchanged during the spin-up and transient simulations. Each time parameters are modified, we re-run the spin-up from the beginning to maintain consistency. This simulation protocol is widely used in previous model optimization studies (Ma et al., 2024; Peylin et al., 2016).

2.5 Optimization processes

Based on the description of the ORCHIDEE processes in Sect. 2.2, we conducted a sensitivity analysis of various model parameters on AGB, GPP, and biomass mortality rates using a one-at-a-time approach. This method is widely used for sensitivity analysis, where each parameter is varied individually while keeping the others unchanged (Liu et al., 2015). We define VP_i (variation percentage,  %) as an indicator for the sensitivity of model output variables i (i.e., AGB, GPP, and biomass mortality rates). Following a previous study (Liu et al., 2015), we varied each parameter by ±10 %, and VP is calculated as follows:

$$\begin{array}{}\text{(4)}& \mathrm{VP}={\displaystyle \frac{\left|{\mathrm{run}}_{+\mathrm{10}\mathit{\%}}-{\mathrm{run}}_{-\mathrm{10}\mathit{\%}}\right|}{{\mathrm{run}}_{\mathrm{ref}}}}\times \mathrm{100},\end{array}$$

where run_+10 % and run_−10 % are the model results with +10 % and −10 % variations in a given parameter while fixing other parameters, and run_ref is the reference model result with default parameters. We selected parameters related to self-thinning, photosynthesis, carbon allocation, and turnover rate in the model. AGB and the mortality rate are the most sensitive to β, α, and η, while GPP is most sensitive to η, SLA, and k_lsmin (Table 1).

Table 1 Variation percentage (VP) of AGB, GPP, and mortality rate for different parameters.

To optimize the model parameters for matching the spatial pattern of GPP, AGB, and biomass mortality rates, we constructed a loss function to optimize two key parameters: α and η. We chose the two parameters for the following reasons. (1) Although β is more sensitive to AGB and mortality rate than α, α has a wider range of values compared to parameter β (Brienen et al., 2015; Yu et al., 2024). Regardless of which parameter is optimized, the observations can be matched. To avoid the issue of equifinality when optimizing both parameters, we chose to adjust α while keeping β fixed at its default value, considering that α is less constrained than β in previous studies (Yu et al., 2024). (2) η has the highest sensitivity to GPP in the model and is found to be highly variable within the moist tropical biome (Kattge et al., 2009; Ellsworth et al., 2022). (3) We did not try to optimize more parameters because the number of optimized parameters should not exceed the number of observations (GPP, AGB, and biomass mortality rates) to avoid the equifinality issue.

We used the Latin hypercube sampling (LHS) method to select the initial parameter values, a technique widely used in model parameterization (McKay et al., 2000; Yan et al., 2023). LHS can subdivide the full parameter space evenly to generate representative parameter values. This procedure ensures an efficient sampling of the parameter distribution with a suitable number of parameter samples, thus avoiding impracticable computational efforts due to numerous model simulations (McKay et al., 2000). The recommended number of sampled points is typically 10 times the number of dimensions (Loeppky et al., 2009). Accordingly, we used 20 samples as the initial parameter values, generated using LHS. The default η value is 14.08 µmolCO₂ gN⁻¹ s⁻¹ based on the biome average of trait data for tropical forests on highly weathered soils (Kattge et al., 2009), and the α value is 1941 cm (Brienen et al., 2015). In this study, the initial uniform prior distribution ranges for α and η were set to 1000–2800 cm and 8–18 µmolCO₂ gN⁻¹ s⁻¹, respectively, to ensure that all optimized values fell within these boundaries (Fig. S7).

After the initial ORCHIDEE run based on the 20 sets of parameter values, we interpolated the results for GPP, AGB, and biomass mortality rates over the whole domain using quadratic splines (Fig. S8), considering that these variables change monotonically and continuously with α and η. The interpolation achieved R² values of 0.99, 0.98, and 0.93 for AGB, GPP, and biomass mortality rates, respectively, based on leave-one-out validation, where 19 simulations were used for interpolation and the remaining one for validation.

Next, for each 1° × 1° grid cell, we constructed a quadratic loss function to identify the best estimates of α and η that minimize the loss (Groenendijk et al., 2021), defined as:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathrm{Loss}}_{\mathrm{agb}}(\mathit{\alpha },\mathit{\eta })={\left({\displaystyle \frac{{\mathrm{AGB}}_{\mathrm{sim}}\left(\mathit{\alpha },\mathit{\eta }\right)-{\mathrm{AGB}}_{\mathrm{obs}}}{\stackrel{\mathrm{‾}}{{\mathrm{AGB}}_{\mathrm{obs}}}}}\right)}^{\mathrm{2}},\text{(6)}& {\displaystyle }{\mathrm{Loss}}_{\mathrm{gpp}}(\mathit{\alpha },\mathit{\eta })={\left({\displaystyle \frac{{\mathrm{GPP}}_{\mathrm{sim}}\left(\mathit{\alpha },\mathit{\eta }\right)-{\mathrm{GPP}}_{\mathrm{obs}}}{\stackrel{\mathrm{‾}}{{\mathrm{GPP}}_{\mathrm{obs}}}}}\right)}^{\mathrm{2}},\text{(7)}& {\displaystyle }{\mathrm{Loss}}_{\mathrm{mor}}(\mathit{\alpha },\mathit{\eta })={\left({\displaystyle \frac{{\mathrm{Mor}}_{\mathrm{sim}}\left(\mathit{\alpha },\mathit{\eta }\right)-{\mathrm{Mor}}_{\mathrm{obs}}}{\stackrel{\mathrm{‾}}{{\mathrm{Mor}}_{\mathrm{obs}}}}}\right)}^{\mathrm{2}},\text{(8)}& {\displaystyle }\begin{array}{rl}\mathrm{Loss}\left(\mathit{\alpha },\mathit{\eta }\right)& ={\displaystyle \frac{{\mathrm{Loss}}_{\mathrm{agb}}\left(\mathit{\alpha },\mathit{\eta }\right)}{{\mathrm{CV}}_{\mathrm{AGB}}}}+{\displaystyle \frac{{\mathrm{Loss}}_{\mathrm{gpp}}\left(\mathit{\alpha },\mathit{\eta }\right)}{{\mathrm{CV}}_{\mathrm{GPP}}}}\\ & +{\displaystyle \frac{{\mathrm{Loss}}_{\mathrm{mor}}\left(\mathit{\alpha },\mathit{\eta }\right)}{{\mathrm{CV}}_{\mathrm{mor}}}},\end{array}\end{array}$$

where AGB_sim(α,η), GPP_sim(α,η), and Mor_sim(α,η) are the AGB, GPP, and biomass mortality rates simulated by ORCHIDEE model after interpolation; AGB_obs, GPP_obs, and Mor_obs are the AGB, GPP, and biomass mortality rates from the observations; $\stackrel{\mathrm{‾}}{{\mathrm{AGB}}_{\mathrm{obs}}}$, $\stackrel{\mathrm{‾}}{{\mathrm{GPP}}_{\mathrm{obs}}}$, and $\stackrel{\mathrm{‾}}{{\mathrm{Mor}}_{\mathrm{obs}}}$ are the spatial mean AGB, GPP, and biomass mortality rates for the whole basin; and CV_AGB, CV_GPP, and CV_Mor are the coefficients of variation for AGB, GPP, and biomass mortality rates from the observations. We used coefficients of variation as the weight to balance different losses from AGB, GPP, and biomass mortality rates at the regional scale. We selected the α and η values that minimized Loss(α,η) for each grid cell and generated spatial maps of both parameters across the entire Amazon region. Then, we re-ran the ORCHIDEE model forced by the α and η spatial maps obtained from the optimization. For comparison, we also performed the same optimization using regional mean AGB, GPP, and biomass mortality rates from both the observations and the simulations to derive spatially constant α and η values across the entire Amazon region as a baseline (Fig. 1).

2.6 Analyze the influencing variables

Previous empirical studies have found that GPP/NPP and mortality rates are influenced by various factors including climate conditions, wood density, water table depth, and soil properties such as soil texture and soil fertility (Malhi et al., 2004; Quesada et al., 2012; Vicca et al., 2012; Castanho et al., 2013; Campioli et al., 2015; Sousa et al., 2022). We thus selected eight variables to explain the spatial distribution of the optimal α and η values (Fig. S9, Table S3). Mean annual temperature (MAT), mean annual precipitation (MAP), and downward shortwave radiation (SWdown) are obtained from CRUJRA v2.4 (Harris et al., 2020). Maximum cumulative water deficit (MCWD) is a proxy for climatic water supply and indicates the intensity of drought. MCWD is calculated as

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{CWD}}_t=\left\{\begin{array}{l}{P}_t-E_t,\\ \mathrm{if}{P}_t<\mathrm{100}\mathrm{mm}\mathrm{and}t=\mathrm{1}\\ {\mathrm{CWD}}_{t-\mathrm{1}}+P_t-E_t,\\ \mathrm{if}{P}_t<\mathrm{100}\mathrm{mm}\mathrm{and}t>\mathrm{1}\\ \mathrm{0},\mathrm{if}{P}_t\ge \mathrm{100}\mathrm{mm}\end{array}\right.\text{(10)}& {\displaystyle }\mathrm{MCWD}=min\left({\mathrm{CWD}}_t\right)\end{array}$$

with t being the month $\mathrm{1},\mathrm{\dots },\mathrm{12}$, P_t the precipitation for month t calculated from CRUJRA v2.4, and E_t the monthly evapotranspiration, which is fixed at 100 mm (Aragão et al., 2007). MCWD was calculated based on the hydrological year, which begins in May (t=1) and ends in April (t=12) of the following year (Chen et al., 2024). Wood density was sourced from Yang et al. (2024a). Clay fraction from HWSD2 serves as a proxy for soil texture (FAO and IIASA, 2023). Soil total phosphorus content, a key indicator of soil fertility, was derived from Darela-Filho et al. (2024). Water table depth, representing both anoxic conditions and additional water availability beyond precipitation, was obtained from Fan et al. (2013). To avoid multicollinearity, precipitation was excluded from the analysis due to its high correlation with MCWD (Fig. S10).

We used random forest regression to analyze the potential influences of the eight variables on α and η. Random forest regression is a widely used machine learning algorithm that builds multiple decision trees to improve prediction accuracy (Breiman, 2001). We randomly selected 75 % of the data to train the model, leaving the remaining 25 % for validation. The random forest model was optimized with the following parameters for both α and η: ntree = 100, max_depth = 3. We constructed 50 models using different random seeds and averaged their outputs to derive the final results.

SHAP (SHapley Additive exPlanations) values were used to evaluate the contribution of each explanatory variable in the random forest regression (Lundberg et al., 2020). Based on Shapley values in game theory, SHAP values represent the marginal contribution of each variable, considering all possible combinations of feature values (Shapley, 1953). A SHAP value, whether positive or negative, indicates a corresponding positive or negative effect on the model's output. This method offers an interpretable insight into how individual features impact model predictions. We built the random forest model using the Python package “sklearn” and computed SHAP values with the Python package “shap.”

3 Results

3.1 Simulations with optimal spatially constant parameters

The optimal spatially constant parameter values for α and η are 1898 cm and 12.14 µmolCO₂ gN⁻¹ s⁻¹, respectively, to match the AGB, GPP, and biomass mortality rates from observations. Basin-scale averaged values of AGB, GPP, and biomass mortality rates closely aligned with observations (126 ± 7 MgC ha⁻¹, 33 ± 3 MgC ha⁻¹ yr⁻¹, and 2.2 ± 0.2 % vs. 127 ± 29 MgC ha⁻¹, 33 ± 2 MgC ha⁻¹ yr⁻¹, and 2.0 ± 0.6 %, respectively). However, simulations failed to capture spatial heterogeneity, such as the higher biomass in the Guiana Shield or the higher biomass mortality rates in the southeastern Amazon (Fig. 2). The R² values between the simulated AGB, GPP, and observations data were less than 0.04, and the correlation between simulated biomass mortality rates and observations was even negative (Fig. 3a–c). Climate forcing, wood density, and soil texture had a limited impact on the modeled spatial variation of AGB, resulting in flat gradients of forest AGB and biomass mortality rates. In the ORCHIDEE model, higher MCWD is generally associated with higher GPP (Fig. 3b), but in reality GPP does not vary significantly across the entire Amazon region (Fig. 2b). On the contrary, MCWD has little influence on modeled tree mortality, despite its observed association with higher mortality rates (Fig. 3c). These discrepancies suggest that spatially varying parameters are necessary to better match the modeled AGB, GPP, and mortality rates with observational data.

Figure 2 Spatial patterns of AGB, GPP, and biomass mortality rates from observations (a–c), from the model with spatially constant parameters (d–f), and from the model with spatially varying parameters (g–i).

Figure 3 Comparison of model-simulated AGB, GPP, and mortality rates with observations. Results from model simulation with spatially constant parameters (a–c) and with spatially varying parameters (d–f). The dashed red line is the 1 : 1 line. The black solid line is the best fit between modeled results and observations. The color of the dots represents the value of MCWD to show its correlation with AGB, GPP, and biomass mortality rates.

3.2 Simulation with optimal spatially varying parameters

The model with optimal spatially varying parameters successfully captured the spatial patterns of AGB, GPP, and mortality rates, including the higher biomass in the Guiana Shield, the lower GPP in the northwestern Amazon, and the higher mortality in the southern Amazon (Fig. 2). The R² values between the observed data and model-simulated AGB, GPP, and mortality rates improved substantially, increasing from near zero to 0.82, 0.79, and 0.73, respectively (Fig. 3). The spatially varying α values were higher in the northeast and lower in the southeast, while η values were higher in both the northeast and southern regions (Fig. S11).

Optimizing both α and η is necessary. When only α is optimized, the model-simulated AGB and mortality rates improved a lot compared with observations, but GPP shows a weak correlation with observed GPP (R² = 0.01, Fig. S12a–c). This is mainly because GPP is not sensitive to α (Table 1). Similarly, the VP of η is larger for GPP and biomass mortality rates than AGB (Table 1). Optimizing only η results in no correlation between model-simulated AGB and observed AGB (R² = 0.0, Fig. S12d–f), although simulated GPP and mortality rates had a stronger correlation with the observations. The spatial pattern of α shows a strong positive correlation with Yu-Biomass and a negative correlation with observed mortality rates. In contrast, η is positively correlated with GOSIF GPP and mortality rate data (Fig. S13).

3.3 Spatial correlation between optimal parameter maps and influencing factors

We further applied explainable machine learning (i.e., random forest regression and SHAP) to find the influence of various factors on the spatial pattern of optimal parameters. The seven variables explained 47 % of the variation in α and 48 % in η on testing data. The Shapley results show non-linear relationships between the optimal parameters and influencing factors, including climatic, plant, and edaphic properties. For example, when MCWD is below −300 mm, α increases sharply with rising MCWD, but once MCWD exceeds −300 mm, its impact on α becomes minimal. Similarly, the negative impact of MCWD on η becomes positive when MCWD falls below −100 mm (Fig. 4b). Clay fraction has a greater impact on η when it exceeds 40 %, particularly in soils classified as clay or silty clay, which accounts for 35 % of all grids (Figs. 4e, S1; Table S2). Wood density shows a strong positive correlation with α when it exceeds 0.61 gC cm⁻³ (Fig. 4d), while the effect of water table depth on α follows a dome-shaped pattern, where both deep and shallow water table depths reduce α (Fig. 4g). Wood density shows the highest feature importance for α, based on mean absolute SHAP values (Fig. 5a). Additionally, MCWD emerges as the most important feature for η (Fig. 5b).

Figure 4 SHAP values of α and η for seven explanatory variables: (a) temperature; (b) maximum cumulative water deficit (MCWD); (c) downward shortwave radiation (SWdown); (d) wood density; (e) clay fraction; (f) total available phosphorus; (g) water table depth. The shaded area indicates the standard deviation of SHAP values. The gray bars represent the distributions of explanatory variables.

Figure 5 Variable importance characterized by the mean SHAP values of seven explanatory variables affecting α (a) and η (b).

4 Discussion

4.1 Implications from the inconsistency between model results and observations

By incorporating spatially varying parameters, the simulation shows much better agreement with the observations. However, some grid cells still exhibit underestimation or overestimation of AGB, GPP, or biomass mortality rates (Fig. S14), reflecting trade-offs made during the optimization of the loss function. In the ORCHIDEE model, GPP, biomass, and mortality rates are tightly coupled. The amount of dead woody biomass is determined by both biomass and mortality rates, which should equal the fraction of NPP allocated to wood under an equilibrium state. While carbon use efficiency (CUE = NPP/GPP: 33.7 % ± 1.2 %) and the woody NPP fraction (f_NPP: 30.8 ± 1.1 %) do not vary a lot in the model, a grid with low biomass and low mortality rate should exhibit a correspondingly low GPP (Fig. S15). When observations are inconsistent with each other (e.g., low biomass, low biomass mortality rate, but high GPP from observations), optimizing α and η alone is insufficient to accurately simulate AGB, GPP, and biomass mortality rates simultaneously.

There are two potential reasons for the remaining model–observation inconsistency. First, the model may overlook important processes. For instance, the model currently only accounts for woody biomass loss through tree mortality. However, Zuleta et al. (2023) found that 42 % of woody biomass loss occurs due to damage to living trees, such as branch falls, trunk breakage, and wood decay. In the ORCHIDEE model, aboveground woody biomass is divided into two pools: sapwood and heartwood. To improve the accuracy of biomass mortality estimates, a separate biomass pool for branches should be incorporated in future model developments to better account for biomass losses from branches. Besides, while CUE and f_NPP fall within observed ranges, the variability of these factors may be underestimated. Some studies suggest that CUE ranges from 23 % to 45 % (Doughty et al., 2018), whereas our model estimates it at 32 % to 42 %. Similarly, the spatial standard deviation of observed f_NPP (2.9 %) exceeds that of our modeled results (1.1 %; Yang et al., 2021).

Second, the observational data may contain inherent biases. Mortality rates derived from field observations may be overestimated because field plots are often located near roads and rivers, where disturbances are more frequent (Saatchi et al., 2015). Additionally, satellite AGB may be underestimated in dense forests due to sensor saturation (Yu et al., 2023). GPP data carries the greatest uncertainty, largely due to the limited observations in the Amazon basin. We compared five GPP datasets across the Amazon basin (see Sect. 2.3). While the GOSIF and FLUXCOM_RS datasets show higher GPP in the southwestern Amazon and lower GPP in the central region, the spatial patterns differ across other regions and datasets (Fig. S16). The highest correlation coefficient between the GPP datasets is 0.46, observed between GOSIF and FLUXCOM_RS, while correlations among the other datasets are generally weaker (Fig. S3b).

We conducted sensitivity tests using additional GPP and biomass datasets. Four other GPP datasets were used to optimize the parameters α and η. After optimization, both AGB and GPP showed strong agreement between the model and observations, demonstrating the robustness of our approach (Fig. S17). The linear correlations between the eight influencing factors and the optimal parameters are generally consistent with previous results (Figs. S18, S19), with only a few exceptions. For example, when using GPP data from FLUXCOM_RS+METEO, MCWD and wood density are positively correlated with η, which is different when using GOSIF and other GPP datasets.

4.2 Implications for the ORCHIDEE model from the random forest regression

The uncertainties in simulated carbon stocks and fluxes arise from model structure, parameterization, and forcing data (Bonan and Doney, 2018). During the optimization of α and η, we assumed that all uncertainties stem from parameterization. By analyzing the correlation between the optimal parameters and influencing factors, our results showed great consistency with existing knowledge and highlighted important processes that are either missing or poorly represented in the ORCHIDEE model when simulating the Amazon forest.

For instance, a threshold of −300 mm in MCWD was identified in the partial dependence analysis between α and MCWD. When MCWD drops below −300 mm, further decreases lead to lower α, indicating higher mortality (Fig. 4b). This result aligns with previous findings showing a negative correlation between AGB mortality rates and MCWD only when MCWD is particularly low (Sousa et al., 2022). Although our model incorporates drought-induced tree mortality, the higher tree mortality under lower MCWD conditions is likely underestimated. The dome-shaped relationship between water table depth and α is also in line with earlier studies (Fig. 4g). Shallow water table depths are associated with low biomass, as excessive water can limit oxygen flow to roots, restricting plant growth (Costa et al., 2023). In waterlogged conditions, weakened root systems and reduced soil cohesion increase the risk of tree falls (Ferry et al., 2010). Conversely, deeper water table depths elevate the risk of drought-induced mortality, thereby increasing mortality rates (Chen et al., 2024). However, water table depth is not explicitly represented in our model. We also found that higher wood density is associated with higher α, which is in line with previous findings that forests with higher wood density tend to have lower mortality rates and higher biomass (Fig. 4d; Keeling and Phillips, 2007; Brienen et al., 2020). However, our model does not account for the feedback between wood density and tree mortality rates.

Furthermore, we observed that η increases when MCWD falls below −100 mm. In the model, reduced water availability limits transpiration, leading to a decrease in GPP (Yao et al., 2022). However, the observed GPP exhibits a weak negative correlation with MCWD (Fig. S20). This suggests that grid cells with more negative MCWD require a higher η to simulate higher GPP. Given the unclear mechanisms driving the correlation between MCWD and GPP in the Amazon forest, it remains uncertain whether η should increase as MCWD decreases. We also found that η increases when the clay fraction exceeds 40 %. Although precipitation is distributed evenly across soil types with varying clay fractions, the runoff fraction increases sharply when the clay fraction surpasses 40 % (Fig. S21), due to the very low hydraulic conductivity of clay soils (Table S2). Besides, while available phosphorus is thought to influence woody production (Quesada et al., 2012), total phosphorus does not seem to be a critical factor in explaining the variation in optimal α and η (Fig. 5). In fact, the correlation between total available phosphorus and GPP is quite weak (R² = 0.05), and the correlations with AGB or mortality rates are even weaker (Fig. S22). This may be due to uncertainties in the GPP data (Figs. S3, S16).

Our random forest model results show that the seven selected variables explained 47 % of the variation in α and 48 % in η on testing data. The limited explanatory power of the random forest model may be attributed to the uncertainties in the current environmental driver datasets. When we re-ran the model using parameters predicted by the random forest regression, it still captured the spatial gradients of carbon pools and fluxes, but with reduced accuracy compared to the optimized results. The modeled AGB, GPP, and mortality rates with spatially varying parameters from the random forest explain 47 %, 25 %, and 22 % of the spatial variation in the observations, respectively (Fig. S23). In the future, by including more relevant and accurate influencing factors, we may be able to enhance the explanatory power of the random forest model and improve ORCHIDEE by using the random forest algorithm to predict parameter variation across space, rather than relying on prescribed parameter maps. Ultimately, linking environmental drivers directly to model parameters in ORCHIDEE could enable a more mechanistic parameterization.

4.3 Strengths of the optimization method

In this study, we significantly improved the simulations of spatial gradients and mean values of AGB, GPP, and biomass mortality rates in the Amazon forest by optimizing two key parameters related to photosynthesis and self-thinning. Our method accounted for parameter variations within a PFT (tropical broad-leaved evergreen trees), an important factor that has largely been overlooked in other models (Butler et al., 2017). Besides, unlike previous studies that optimized numerous parameters simultaneously, our approach provided an efficient parameterization framework by focusing on the sensitive parameters for each pixel. In the future we can incorporate additional observations, such as LAI and tree height, as their accuracy improves. This method is also applicable to other DGVMs. Except for optimizing the parameters, this approach also helps identify the missing or underrepresented processes in the model. We also tried to link the optimal parameters with other plant traits or environmental factors. By doing this, we may eventually improve the model by incorporating relationships between different traits (e.g., linking wood density to α).

4.4 Limitations of this work

Our parameter optimization framework does not aim to optimize temporal variation or long-term trends (see Text S3). Besides, the input datasets used for model calibration carry significant uncertainties. More importantly, even with accurate observation data, our optimized parameter could be incorrect, as the mismatch between model results and observations can arise from the model structure and forcing data (Bonan and Doney, 2018).

In addition to the missing processes mentioned earlier, this version of the ORCHIDEE model also lacks several key processes, such as phosphorus–carbon interactions and species-based functional group diversity (Reed et al., 2015; Sun et al., 2021; Rüger et al., 2020). We ignored the spatial variability of leaf nitrogen content by fixing it uniformly across the Amazon forest. However, V_cmax depends on both η and the leaf $\mathrm{C}/\mathrm{N}$ ratio, the latter also exhibiting significant variability (Ellsworth et al., 2022).

In the context of climate change, the Amazon forest is increasingly exposed to disturbances such as droughts, fires, and windthrow events (Chen et al., 2024; Flores et al., 2024; Feng et al., 2023). Although the model accounts for drought-induced mortality, it does not include other disturbances like fire, windthrow, and logging in the Amazon basin. In this study, we focused exclusively on intact forests, as defined by TMF data (Vancutsem et al., 2021). Incorporating natural and anthropogenic disturbances into this ORCHIDEE version is underway (Naudts et al., 2015; Yue et al., 2014; Marie et al., 2024; Chen et al., 2018) and might turn out to be essential in future work to fully capture the complexity of ecosystem dynamics across the entire Amazon rainforest.

5 Conclusion

Our study developed an efficient parameter optimization framework that integrates observation-based data for AGB, GPP, and mortality rates to optimize the spatial variation of parameters in the ORCHIDEE model. By incorporating spatially varying parameters, the model effectively simulates the spatial variation of AGB, GPP, and mortality rates. The spatial patterns of the optimal parameters are generally reasonable, revealing previously underrepresented processes in the model. Given the complexity of the Amazon forest and the practical limits on model complexity, applying spatially varying parameters can help accurately represent its spatial variability.