Leaf area index is an important parameter that consistently monitors vegetation structure dynamics over large spatial and temporal scales. This study uses three satellite-derived long-term global LAI products: GIMMS LAI3g, GLASS LAI, and GLOBMAP LAI. We re-sampled all three LAI datasets from their native resolution into 0.5° × 0.5° special resolution using the nearest neighbor algorithm and interpolated to the hourly temporal resolution to force the model. All three datasets have been validated and widely used to monitor terrestrial vegetation dynamics (Dardel et al., 2014; Piao et al., 2015; Zhu et al., 2016, 2017; Jiang et al., 2017).

We developed a theoretical framework for decomposing the terrestrial carbon, nitrogen, and phosphorus stock into some traceable components based on the biogeochemical principles of the terrestrial carbon cycle. For example, the terrestrial carbon cycle can be generally described as the following processes, which includes carbon enters the ecosystem via plant photosynthesis; photosynthetic carbon is then allocated among plant pools; part of this carbon is consumed by respiration, and the remainder is further transferred to litter and soil carbon pools; lastly, the carbon in the litter and soil pools is decomposed and back into the atmosphere (Luo and Weng, 2011; Luo et al., 2022). Plants assimilate the nitrogen and phosphorus from mineral soils and then transfer following the same flow with organic matter in terrestrial ecosystems. Therefore, following the approach developed by Xia et al. (2013), the biogeochemical cycle processes can be mathematically represented by three matrix equations:

where X(t) = (X₁(t), X₂(t), …, X₉(t))^T , N(t) = (N₁(t), N₂(t), …, N₉(t))^T , and P(t) = (P₁(t), P₂(t), …, P₉(t))^T are 9 × 1 vector, which describes the C, N and P pool size of leaf, root, wood, metabolic litter, structural litter, coarse wood debris (CWD), fast soil, slow soil and passive soil pool at time t in CABLE model. B i = b 1 b 2 b 3 0 … 0 T i = C , N , P is a vector of allocation coefficients of C, N, and P among different plant pools. For the C allocation, the partitioning coefficients from NPP to root and wood carbon pools are equations of availably of light, nitrogen, and water, respectively. Then the rest of NPP goes to the leaf pool. For the N and P, the allocations of N and P uptake among different plant pools are calculated based on the proportion to each pool’s demand of N and P. The inputs of C, N, and P are represented by U i t i = C , N , P . U C is the fixed carbon by canopy photosynthesis, i.e., net primary production (NPP). U N and U P are the assimilated N and P via plant uptake from soil minerals. ξ t is a diagonal matrix, the diagonal components representing the environmental scalars (such as temperature and soil moisture) effects on carbon decomposition rate at time t. C is a 9 × 9 diagonal matrix with diagonal entries by 9 × 1 vectors c = c 1 c 2 … c 3 T . The diagonal elements indicate the C decay rate for each pool. A is a transfer coefficients matrix, which can quantify how much carbon can be transferred among different pools. Therefore, the first term on the right side of Equation (1), B i U i i = C , N , P , describes the C, N, and P inputs and allocation among different plant pools, and the second term on the right represents the transfer and exit rates (Xia et al., 2013; Luo et al., 2017).

By letting Equation (1) equal zero, we obtained the C, N, and P pool size at steady-state as the product of ecosystem residence time and influx:

where X s s , N s s , and P s s are vectors that includes all the organic pools at steady state. U C , U N , and U P are the ecosystem C, N, and P influx at steady state. U C represents NPP in this study, which can be further decomposed to gross primary production (GPP) and carbon use efficiency (CUE) based on some previous studies (Bradford and Crowther, 2013; Xia et al., 2017). The term A ξ C − 1 B i i = C , N , P in Equation (2) are vectors of the residence time of individual pools as:

The ecosystem-level residence time of C, N, or P was summed from all the individual pools. The ecosystem residence time can be determined by multiple ecological processes, such as allocation (i.e., the B vector), carbon transfer among different pools (i.e., the A matrix), decomposition rates (i.e., the C matrix), and the environmental scalars (i.e., ξ ). The scalar ξ usually consists of the temperature ( ξ T ) and water scalars ( ξ W ) as:

Generally, the parameters of B, A, and C matrices are preset in a specific model according to model structure, soil properties, and vegetation characteristics (Zhou et al., 2018). By rearranging Equation (3), we can further decompose the residence time to the environment scalar and the corresponding preset parameters:

The residence time can be further decomposed to environment scalars and baseline residence time vectors. The equation of baseline residence time can be expressed as:

The Australian Community Atmosphere Biosphere Land Exchange (CABLE) model version 2 is a global land surface model that can simulate biophysical and biogeochemical processes (Wang et al., 2010, 2011). It includes five submodels: (1) radiation, (2) canopy micrometeorology, (3) surface flux, (4) soil and snow, and (5) ecosystem respiration, and it also incorporates carbon (C), nitrogen (N), and phosphorus (P) cycles. CABLE has been widely evaluated by other global observations (Piao et al., 2015), eddy-flux measurements (Best et al., 2015), and manipulated field experiments (De Kauwe et al., 2014). This model can be applied to attribution analysis (Zhang et al., 2016) or plant feedback effects (Lei et al., 2019) by doing a series of simulation experiments. The default settings of LAI are prognostic in the CABLE model, but a switch can control them. When the switch is turned on, the prognostic LAI can be calculated as the product between specific leaf area (SLA) and leaf biomass. SLA and the phenology phases (used to determine the leaf growth) are prescribed for each plant functional type.

In this study, we turned the switch off and replaced the modeled LAI with data from three satellite-observed products (GIMMS LAI3g, GLASS, and GLOBMAP). Based on the traceability analysis approach, we performed three simulations to diagnose the source of uncertainty in the biogeochemical cycle caused by LAI. The CABLE model was first spun up with the C-N-P coupled schemes to the steady state in 1,900 using a semi-analytical solution (Xia et al., 2012). The forcing data (Zhang et al., 2016) used to spin up the model concludes seven 6-hourly meteorological forcing variables (i.e., temperature, precipitation, downward shortwave radiation, downward longwave radiation, specific humidity, pressure, and wind speed) from the CRUNCEP version 5 (New et al., 1999, 2000, 2002). Using spin-up results as an initial value, we run the model from 1901 to 1981. After that, we performed three simulations by replacing the modeled LAI with three satellite-based LAI products (GIMMS LAI3g, GLASS, GLOBMAP), respectively. Lastly, we spun up CABLE to a steady state forced by the satellite-derived LAI datasets and time-variant CO₂ concentration from 1982 to 2011 (Supplementary Table S1).

The estimated global mean LAI was 1.24 ± 0.18 m² m⁻² (mean ± SD) across the three data products. The simulated global stocks of C, N, and P were 11.0 ± 1.4 C kg C m⁻², 504.3 ± 68.1 g N m⁻², and 97.7 ± 13.6 g P m⁻², respectively. All three simulations showed the highest mean LAI in tropical regions among the eight biomes (Figure 2A). Furthermore, the variability of LAI across three simulations was also higher in the tropics than in any other area (Figure 2B). However, the simulated element storage (i.e., C, N, and P) showed a divergent spatial pattern in magnitude and variability compared with LAI (Figures 2C–H). Specifically, northern high-latitude regions showed the highest annual mean element storage (12.0 kg C m⁻² for the C storage, 531.3 g N m⁻² for the N storage, and 93.1 g P m⁻² for the P storage) compared with other climate regions. In contrast, tropical regions showed a relatively high C storage (10.0 kg C m⁻²) but low storage of N (400. 2 g N m⁻²) and P (78.5 g P m⁻²) (Figures 2C,E,G). A similar distribution was also found in the spatial variability of elements (Figures 2D,F,H). In addition, it is noted that the high disagreement in P storage across three simulations is mainly located in the regions covered by herbaceous vegetation, such as eastern Australia, southeastern South America, and southern Africa (Figure 2H).

The ecosystem element storage can be decomposed into the corresponding element residence time ( τ C , τ N , τ P ) and the influx rate [U_C (NPP), U_N, and U_P] based on the traceability framework according to Equations (2, 3). Deciduous needle leaf forest (DNF) had the largest ensemble annual mean results of ecosystem C (86.0 kg C m⁻²), N (424.3 g N m⁻²), and P (24.4 g P m⁻²) storage among the eight biomes, resulting from its longest residence time ( τ C , 154.3 years; τ N , 78.4 years; τ P , 65.5 years) and mediated element input (NPP, 557.4 g C m⁻² yr.⁻¹; U_N, 5.4 g N m⁻² yr.⁻¹; U_P, 0.37 g P m⁻² yr.⁻¹). However, the order of the size in C, N, and P storage was inconsistent in the remaining seven biomes (Supplementary Table S2). Shrubland had the lowest C storage (10.7 kg C m⁻²) as a result of the smallest NPP (257.6 g C m⁻² yr.⁻¹) and a moderate τ C (41.4 years). Evergreen needle leaf forest (ENF) had a relatively long τ N (50.4 years) and low U_N (3.1 g N m⁻² yr.⁻¹). While evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), and C3 grassland (C3G) had a short τ N (~25.0 years) and relatively high U_N (~9.0 g N m⁻² yr.⁻¹), leading to a moderate N storage. For P storage, EBF, DBF, and C3G had a relatively short τ P (~23.9 years) and relatively high U_P (~0.50 g P m⁻² yr.⁻¹), resulting in a moderate P storage. Although Tundra had a relatively long τ N (53.0 years) and τ P (50.9 years), it still has the lowest N (16.5 g N m⁻²) and P (1.4 g P m⁻²) storage as the result of the smallest U_N (0.31 g N m⁻² yr.⁻¹) and U_P (0.03 g P m⁻² yr.⁻¹). Additionally, the carbon influx via canopy photosynthesis (NPP) can be further decomposed into GPP and CUE (Xia et al., 2017). The results showed that EBF had the highest annual mean GPP (3312.2 g C m⁻² yr.⁻¹), followed by DBF (2103.1), C4 grassland (C4G, 1873.4), C3G (951.3), ENF (912.6), DNF (871.2), Shrubland (511.6) and Tundra (370.8) regions (Supplementary Figure S2 and Supplementary Table S2). The ranges of carbon use efficiency of all eight biomes from 0.37 to 0.71.

The ecosystem element storage derived from three satellite-based LAI differed among the eight biomes. For example, ENF, EBF, and DBF had similar C, N, and P storage for the three simulations (Figure 3). Although Shrubland and C3G had comparable values of element residence time ( τ C , τ N , and τ P ), their different element uptake rates (NPP, U_N, and U_P) led to variations in element storage across the three simulations. In addition, compared with the simulations derived from GIMMS LAI3g and GLASS, the simulation derived by GLOBMAP LAI had the smallest NPP, U_N, and U_P (Figure 3). Tundra and DNF had comparable C storage across three simulations. However, the N and P storage magnitude varied widely in these regions. The differences in N and P storage across the three simulations are mainly due to τ N and τ P (Supplementary Table S2).

Ecosystem residence time can be decomposed into baseline residence time and environmental scalars. The baseline residence time is determined by the carbon transfer coefficients matrix (A matrix), decomposition rates (C matrix), and the allocation coefficients (B vector) among different plant element pools based on Equation 6. Considering that the element cycles (C, N, and P) share the same A and C matrix, the difference in baseline element residence time is mainly caused by the element allocation coefficients ( B C , B N , B P ). The baseline C residence time ( τ C ′ ) is the longest compared with baseline N residence time ( τ N ′ ) and baseline P residence time ( τ P ′ ) (Supplementary Figure S3) in all the eight biomes. Deciduous needle leaf forest has a relatively long τ C ′ (21.8 years), which is almost three times that of τ N ′ (8.2 years) and τ P ′ (6.8 years). Baseline element residence times were similar in C3 grassland ( τ C ′ , 4.6 years; τ N ′ , 3.7 years; τ P ′ , 3.6 years), C4 grassland ( τ C ′ , 5.6 years; τ N ′ , 4.3 years; τ P ′ , 4.3 years), and shrubland ( τ C ′ , 8.7 years; τ N ′ , 7.1 years; τ P ′ , 7.1 years) regions. Additionally, τ C ′ were similar to each other in three LAI-derived simulations. While τ N ′ and τ P ′ differed substantially across three simulations, especially in ENF, DNF, and tundra regions (Figure 4).

Environmental scalar ( ξ ) can regulate ecosystem residence time by limiting the decomposition rates of litter and soil organic pools. By decomposing environmental scalar into temperature ( ξ T ) and water scalar ( ξ W ), we can found that the multi-year mean of ξ T ranges from 0.09 for deciduous needle leaf forest to 0.70 for C4 grassland, while ξ W ranges from 0.60 for tundra and 0.92 for C4 grassland (Figure 5). The global average of ξ T (0.38) is considerably lower than that of ξ W (0.82). In general, ξ T dominates the difference in ξ among eight biomes compared with ξ W .

We decomposed the variations of element storage into several traceable components on each vegetated grid by a traceability analysis approach. Figure 6 shows the dominant traceable component for each grid which explains the greatest contrition to the element storage variation. The results showed that the element influx rates (i.e., NPP, U_N, and U_P) are the primary uncertainty source in 67.6, 93.2, and 93.0% of the vegetated grid cell for the C, N, and P storage, respectively. By further tracing the modeled variation of NPP into GPP and CUE, we found that GPP and CUE explained 91.6 and 8.4% of the variation across simulations, respectively. The baseline residence time has a larger uncertainty contribution than the environmental scalars (Figure 6). Specifically, the contributions of baseline element residence time (i.e., τ C ′ , τ N ′ , and τ P ′ ) to the variation of ecosystem element residence time (i.e., τ C , τ N , and τ P ) is 75.3, 91.2, and 91.8%, respectively. While the contribution of ξ is 24.7% for the variation of τ C , 8.8% for τ N , and 8.2% for τ P . In addition, the contributions of ξ T and ξ W to the variation of element storage are relatively small compared with other contributors, i.e., ξ T and ξ W only contribute 3.8 and 4.2% for C storage, 0.32 and 0.28% for N storage, and 0.28 and 0.29% for P storage variation, respectively.