Implications of the steady-state assumption for the global vegetation carbon turnover

Vegetation carbon turnover time (τ) is a central ecosystem property to quantify the global vegetation carbon dynamics. However, our understanding of vegetation dynamics is hampered by the lack of long-term observations of the changes in vegetation biomass. Here we challenge the steady state assumption of τ by using annual changes in vegetation biomass that derived from remote-sensing observations. We evaluate the changes in magnitude, spatial patterns, and uncertainties in vegetation carbon turnover times from 1992 to 2016. We found the robustness in the steady state assumption for forest ecosystems at large spatial scales, contrasting with local larger differences at the grid cell level between τ under steady state and τ under non-steady state conditions. The observation that terrestrial ecosystems are not in a steady state locally is deemed crucial when studying vegetation dynamics and the potential response of biomass to disturbance and climatic changes.


Introduction
One of the largest uncertainties in Earth system models is in quantifying how the carbon uptake by terrestrial ecosystems will respond to changes in climate (Friedlingstein et al 2006, Friend et al 2014).As an emergent ecosystem property that partially determines carbon sequestration capacity, the vegetation biomass turnover times (τ ) have been used as a diagnostic metric to quantify the feedback between the carbon cycle and climate (Carvalhais et al 2014, Thurner et al 2016).However, there is a large uncertainty in the simulations of vegetation carbon stock as well as τ across earth system models, indicating different representations of the response of vegetation to future climate change (Friend et al 2014).Furthermore, our current understanding of τ and the vegetation dynamics is limited due to the lack of long-term observations of changes in vegetation.
As a result, the estimation of τ has relied so far on the assumption that the vegetation carbon in an ecosystem would eventually reach a steady state (steady state assumption, hereafter SSA) at which the net change of vegetation biomass becomes zero (∆C veg = 0), or so small compared to the total biomass that becomes negligible.The SSA has been shown to be a useful assumption at a large spatial scale.However, at local scales, an ecosystem is unlikely to maintain a steady state due to the influences from external factors such as disturbances and climate variability (Ge et al 2019).It is still unknown whether the SSA can hold at local spatial domains and how much the difference it can make to the τ estimation if one neglects the temporal changes in vegetation carbon.
In this study, we used estimates of annual changes in vegetation carbon derived from a multi-decadal dataset (Besnard et al 2021, Santoro et al 2022) and global estimations of gross primary productivity (GPP) that are driven by meteorological observations (Tramontana et al 2016, Jung et al 2020), for estimating and comparing τ estimates that are derived from SSA and non-steady-state assumption (hereafter NSSA), respectively, at local, biome and global scales.The validity of SSA was evaluated in different spatial domains to better quantify the effect of spatial scales on the patterns of carbon turnover times.

Data and methods
In this section, we first introduce the datasets used to estimate τ including above-ground vegetation biomass, below-ground biomass and GPP.We used a forest canopy cover dataset to examine the relationship between the changes in τ and tree canopy cover.Then the calculations of τ using three methods are introduced next with detailed explanations.

The multi-decadal estimates of AGB dataset
Annual above-ground biomass (AGB) estimates were derived from C-band satellite radar signals between 1992 and 2016 with a pixel size of 25 km (Santoro et al 2022).The very dense time series of observations by the European remote sensing WindScatterometer, the MetOp Advanced SCATterometer (ASCAT), and the Envisat advanced synthetic aperture radar (ASAR) were used to maximize the information content of forest structure in the signal, allowing for AGB estimates of higher accuracy compared to values obtained from a single observation (Santoro et al 2022).The annual estimation of AGB is obtained by synthesizing all daily observations of the radar backscatter at one location in a pixel (0.25 • × 0.25 • ), enabling the inference of a continuous time series of AGB estimation.By adapting the AGB retrieval method in time and space and computing a weighted average of individual AGB estimates, the annual AGB estimates were less impacted by data noise, instantaneous moisture conditions, precipitation, and snow cover (Santoro et al 2011).

Estimation of total vegetation carbon stock
The total vegetation biomass consists of AGB and below-ground biomass (BGB).We estimated BGB from the AGB time series by scaling with the rootshoot ratio (Huang et al 2021), R rs : Therefore the total vegetation biomass is: (1)

Estimation of GPP
We used estimations of GPP from the FLUXCOM project in which different machine learning approaches were applied to upscale global energy and carbon fluxes from eddy covariance flux measurements (Tramontana et al 2016, Jung et al 2020).In this study, GPP annual estimates driven by meteorological observations and remote sensing observations at the spatial resolution of 0.5 • and the time period from 1992 to 2016 are used as carbon influx into the vegetation carbon pool.The dataset was resampled (bilinear interpolation) to the spatial resolution of 0.25 • to match the AGB dataset.

Forest tree canopy cover change
Tree canopy cover (vegetation that is greater than 5 meters in height) was derived from the advanced very high resolution radiometer remote-sensing measurements (Song et al 2018).The version 4 long term data record (LTDR) was used to generate the data on tree canopy coverage from 1982 to 2016.Daily LTDR surface reflectance data were used to compute the normalized difference vegetation index (NDVI) at each pixel (0.05 • × 0.05 • ).Maximum NDVI composition was then used to obtain adjusted annual phenological metrics, which were used as input to supervised regression tree models to generate the annual product of tree canopy coverage.

Estimation of τ under steady state
Changing C veg over time is determined by the uptake of carbon and turnover times: C veg is the vegetation carbon stock.Assuming that the vegetation carbon pool is in a steady state, i.e. the change in C veg over time (dC veg /dt) equals zero, then vegetation carbon turnover times can be calculated as the ratio between vegetation carbon stock and GPP: Here τ SSA is calculated pixel-wise by using annual mean C veg and GPP over the period of 1992-2016.

Estimation of τ under non-steady state
Compared with the estimations of τ under steadystate assumption, the changes in C veg over time are considered (dC veg /dt ̸ = 0) when estimating τ under non-steady state (τ NSSA ).To derive a robust estimation of τ NSSA at each grid cell, we calculated τ NSSA using three different methods to assess the uncertainty built in the τ estimations.

Method 1
We estimate ∆C veg by calculating the annual difference of C veg between year t and year t−1.Then, a τ estimate can be derived for each year by applying GPP and ∆C veg at year t.Finally, we derive the τ under a non-steady state for each year: N Fan et al

Method 2
In the second method, we estimated the mean ∆C veg using the trend of C veg in a certain period to avoid the influence of outliers on the results.In this way, τ can be inferred as: . (5) Here the ∆C veg, trend is inferred by applying a simple linear regression model (least-square robust fitting) between the response variable C veg and time (C veg ∼ T).The coefficient of T is, therefore, the average annual ∆C veg over the whole period.Thus, the τ under a non-steady state can be estimated with the annual mean values of C veg , GPP, and ∆C veg .

Method 3
In the third method, we infer τ from equation ( 2) by applying a linear regression model (least-square robust fitting) at each grid cell in which (GPP-∆C veg ) is the target variable while C veg is the predictor, enabling annual turnover time to be inferred from the coefficient of C veg (1/ τ NSSA): Here ∆C veg is the difference of C veg between year t and year t-1.GPP is the carbon input in year t and C veg is the total carbon density in year t-1, respectively.

The impact of disturbance and climate variability on vegetation carbon turnover
We used the bootstrap aggregating (bagging) ensemble learning method to improve generalizability and robustness of the predictions of τ NSSA over a single estimator.The original training dataset is bootstrapped to form a subset, i.e. samples are drawn with replacement.Bootstrapping promote diversity of the predictions by increasing ensemble members.The issue of overfitting can be reduced by validating prediction with out-of-bag subset.Finally predictions of different ensemble members are aggregated to generate the combined results.In this study, we used multiple features to represent disturbance, temperature, atmospheric and soil moisture conditions (table 1).All variables are harmonized into the same spatial (0.25 • ) and temporal (annual) resolution as the τ NSSA .To estimate the relative importance of individual features, we employed the Shapley value (Lundberg et al 2017) which measures the average marginal contribution of a player in a cooperative game.The advantage of Shapley value is that it provide the magnitude and direction (negative or positive) of the impact of a feature to the deviation from the average prediction.Therefore, Shapley value provide a more intuitive and informative way to compare the relative importance of different features in a tree ensemble model.

Estimation of uncertainty
The uncertainty of τ NSSA estimation comes from both the uncertainty in the AGB and GPP.The uncertainties of annual estimations of AGB is derived from the weighted sum of the individual daily AGB estimate and a covariance component that accounts for the correlation between errors (Santoro et al 2022).
The uncertainty of GPP derived from three machine learning method and two partitioning method (Jung et al 2020).Therefore, the propagated error (standard deviation) of τ NSSA estimation can be derived from equation (4): where And

Comparison of τ under NSSA and SSA at grid cell level and global scale
The τ values (figure 1) represent the turnover time of the entire forest living vegetation biomass, averaged over the whole period of the observations.The comparison between estimates of τ NSSA using three different methods and τ SSA shows a consistent pattern that carbon turnover processes are far from a steady state at grid cell level (figure 1) in spite of a high global correlation between the spatial patterns (R > 0.98, bottom off-diagonal plots in figure 1).
Although there are differences in the estimations of τ NSSA that derived from the three methods, a similar pattern emerges which shows the difference between τ NSSA and τ SSA is the largest in the temperate and boreal forest ecosystems whereas the differences are substantially smaller in the tropical forest ecosystems.
Our results show a high spatial variability of τ values ranging from 0 to 15 years.The longest turnover times are located in the northern boreal forest ecosystem, where part of the biome has τ values longer than ten years, whereas carbon in the temperate forest ecosystem turnovers over much faster where the τ values are mostly under five years.The assumption that vegetation biomass is in steady state results in an maximum overall bias of τ by 10% (90th percentile), compared to the τ estimates under a non-steady state at the grid cell level (figure 1).This finding indicates that the assumption of steady state imposes a maximum deviation in τ of 10% in global forest ecosystems, although the degree of deviation changes from one region to another.The discrepancies between τ SSA and τ NSSA are substantially higher in the boreal forest interquartile range of (8.62%) ecosystem than in the tropical forest interquartile range of (1.86%) ecosystems indicating that the forests in the tropics are closer to a steady state, whereas assuming SSA in the boreal forest may cause large bias (figure S1).Here we show that the forest biomass at the global scale is roughly in a steady state whereas the SSA can be locally largely violated at the grid cell level, especially in the northern boreal forest ecosystems where the τ values can be substantially underestimated or overestimated if assuming SSA.
In line with a previous study in which the SSAinduced biases are assessed at site level (Ge et al 2019), we show that SSA causes substantial underestimations of τ from 5% (median value) up to 40% (99th percentile) in China during the period of 2011-2016 (figure S2).However, our results show a high heterogeneity where SSA can also cause overestimation of τ .Further analysis shows that the pattern also changes across different periods of time.For instance, there is a contrasting pattern between 2005-2010 and 2011-2016 in which the former is characterized by overestimation of τ induced by SSA in a large part of the southwest China whereas there is a widespread underestimation of τ in the latter.

The effect of large-scale disturbances on carbon turnover times
The disturbances from natural causes or anthropogenic activities can make an ecosystem deviate from a steady state.By estimating carbon turnover times at different periods, we quantified the degree of deviation if disturbances, e.g.deforestation, happened in a forest ecosystem.Figure 2 shows that the pervasive deforestation in the 90 s primarily affected the carbon turnover times in the southeast part of the Amazon, which is known as the 'arc of deforestation' (hereafter AOD, Durieux et al 2003).Our results clearly show τ NSSA is underestimated ranging between 5% (25th percentile) and 24% (1th percentile) lower than τ SSA in the AOD region from 1993 to 1998, indicating anthropogenic activity (mostly deforestation) accelerated the carbon turnover rates to a large extent.Compared with the AOD, forests in the middle of Amazon, where there are less population and disturbances, are closer to a steady state, as shown by the much less difference between τ NSSA and τ SSA .Further analysis shows that tree canopy cover (figure 2, row 2) and C veg (figure 2, row 3) changes decreased mainly during the same period of 1993-1998, whereas the changes in GPP does not follow the trend in the arc of deforestation.These results indicate that the acceleration of turnover times during this period is directly caused by the large decrease in the vegetation biomass, which is intimately associated with a decrease in forest cover in this region.On the other hand, our findings show that the forest ecosystems started to recover during the 1999-2004 period as the vegetation biomass increased by 10% to 20%, in line with the increased tree canopy cover in the AOD region.As a result, the carbon turnover times increased by 10% to 30% during the same period.From 2011 to 2016, the magnitude of changes in τ , C veg and tree canopy cover significantly decreased, indicating the forest ecosystems are closer to a steady state due to less disturbances.These findings show how turnover times and the steady state of the forest ecosystem can be largely affected by anthropogenic activities.

The impacts of climate and disturbance on carbon turnover times
To identify the potential factors that control the temporal changes of τ NSSA and quantify the role of each factor, we investigated the link between τ NSSA and multiple variables that represent different perspectives of climate, fire and forest cover change (table 1).By constructing ensemble regression models with bagging decision trees, we take into account nonlinear interactions among different predictors and estimated the importance of individual features.Our results show that the ensemble model can explain a substantial amount of variance globally (median R 2 = 0.56, figure S3).The analysis of predictor importance shows that climate, including air temperature (mean, maximum and minimum), precipitation, atmospheric & soil moisture deficit and radiation dominate vast areas in the northern boreal forest, part of temperate forest in China and tropical forest in Amazon, Congo basin and Indonesia (figure S4(a)).In contrast, we show a dominant role of forest cover change for explaining the inter-annual changes of τ NSSA in forested regions where are subject to intensive anthropogenic land use/land change and deforestation, such as the AOD region (figure S4(b)).To better understand the role of forest cover change or disturbance on τ NSSA , Shapley values are calculated on top of the ensemble regression model to estimate the contribution of individual features to model predictions (see Methods).We show that the change of forest cover has a large effect on the inter-annual changes of carbon turnover, that is, the increase of forest cover lead to longer turnover times whereas decrease in forest cover leads to faster turnover of carbon (figure 3).This result is consistent with our previous analysis which showed that there is a strong co-varying spatial patterns between turnover times and forest cover change (figure 2).The comparison  of the Shapley values among all features indicate that forest cover change, which is an indicator of disturbance, is the main contributor to the changes of τ NSSA in the AOD region as the magnitude of Shapley value of forest cover is two to three times larger than the second ranking feature.Globally, we find similar effect of forest cover change in other forested regions (figure S5).These results show that, in addition to climate change, the effect of natural caused or anthropogenic induced disturbance have high impact on carbon turnover of vegetation in forested ecosystems.

The effect of spatial scale on the steady-state assumption
We further investigate the effect of spatial scale on the difference between τ NSSA and τ SSA in different biomes by gradually changing the spatial scale from 0.25 • (grid cell level) to 25 • (continental scale) via spatial aggregation as shown in figure 4.Here the difference between τ NSSA and τ SSA at each spatial scale is quantified by the 10, 25, 75 and 90 percentiles of the relative difference between τ NSSA and τ SSA (P10, P25, P75, P90, figure 4).We find that the difference between τ NSSA and τ SSA substantially decreases with increasing spatial scales.The P10 and P90 tropical forests decrease approximately 5%, whereas it decreases approximately 10% in temperate and boreal biomes when the spatial scale increases from grid cell to ecosystem scale.Globally, the difference between τ NSSA and τ SSA is approximately 3% at ecosystem scale, indicating that SSA will cause less errors in estimating carbon turnover times at larger spatial scales.

Discussion
Our findings imply that the two different assumptions, i.e.SSA and NSSA, should be applied based on different ecological principles and spatial scales.The common approach of defining τ as the ratio between carbon stock and carbon influx based on SSA can be justified and properly applied when the changes in net carbon flux are negligible relative to the total carbon stock (Carvalhais et al 2014).Although disturbances from nature or human beings could cause nonsteady-state behavior, neglecting the changes, in some cases, only make a little difference to the quantification of the spatial pattern of τ , which does not hamper the understanding of the dynamics of the terrestrial ecosystem carbon cycle.However, at a grid cell level, neglecting the changes in vegetation carbon (assuming vegetation is in a steady state) may result in locallized large biases.Using three methods, we provide robust estimations of τ under a non-steady state.The comparisons between τ SSA and τ NSSA show high heterogeneity in both space and time.A study (Ge et al 2019) showed large SSA-induced biases on τ estimation in varied ecosystems of China by using the data at ten FLUXNET sites from 2005 to 2015 which is consistent with our results.However, we further show that the magnitude and the signs of the SSA-induced biases are characterized by high spatial heterogeneity and can change in time.This is mainly caused by the changes in vegetation biomass due to climate change or disturbances (figure 2).
By using ensemble regression models with multiple features that represent different perspectives of disturbance and climate variability, we identified the dominant role of local disturbance, which is represented by forest cover change, on the major forest ecosystems (figure S4(b)).We found similar patterns in these ecosystems where the impact of disturbance on vegetation carbon turnover is two to three times higher in magnitude than the climate factors (figure 3).For instance, we quantified the role of forest cover change in the arc of deforestation in Amazon forest where the well-documented anthropogenic-induced disturbance caused large changes in the biomass of rainforest.We show that the effect of large decrease in forest cover can significantly accelerate carbon turnover (decrease in τ ).Interestingly, the effect of increase in forest cover, probably due to recovery of forest from disturbance, have higher impact in turnover, but in the opposite direction, i.e. carbon turnover slow down due to recovery.
We have shown substantial heterogeneity in the degree of validity of the steady-state assumption across space.The comparison between τ SSA and τ NSSA quantitatively shows that most global forest ecosystems are not far from steady state, while the largest differences are in the temperate and boreal forests.However, disturbances could drive the forest ecosystem away from steady-state.Our results show that the arc of deforestation in Amazon Forest have large difference between τ SSA and τ NSSA caused by drastic changes in vegetation biomass (figure 2).These results indicate that applying SSA at the grid cell level can cause substantial bias, potentially leading to locally misleading conclusions based on poor estimation of carbon turnover times.
Furthermore, our study quantified the link between spatial scales and the validity of SSA.Our results imply that SSA is approximately valid at large spatial scales (>15 • or 1500 km), at which scale the differences are much lower (∼5%) at grid cell level.The current understanding of the temporal dynamics of the terrestrial carbon cycle nearly all relies on Earth system models in which the carbon turnover rates which have large discrepancies in carbon pools and turnover among different models (Todd-Brown et al 2013, Friend et al 2014).The estimation of τ under NSSA with observational long-term biomass data provides insights into better understanding and thus modeling turnover rate and its spatial patterns.

Conclusions
Our findings suggest that the SSA is robust at a global scale yet becomes much less realistic locally at the grid cell level as the difference between local τ SSA and τ NSSA can be as large as 20% in some regions.The usage of the SSA would result in a substantial bias of τ , especially in regions with a high degree of disturbance, either from anthropogenic activities or natural causes.We found the impact of disturbance on vegetation carbon turnover can be two to three times larger than climate.However, at a larger spatial scale, the differences in τ estimations at SSA and NSSA significantly decrease because the annual changes in vegetation biomass are small compared with the total amount of biomass.With the novel long-term observations of vegetation biomass, we revealed a detailed picture of the spatial distribution of carbon turnover times under different assumptions and its relationship with spatial scales, which can guide the proper application of the two assumptions on different conditions.

Figure 1 .
Figure 1.Comparison of τ under SSA and NSSA using different methods.The upper off-diagonal subplots show the relative difference between each pair of datasets (column/row).The bottom off-diagonal subplots show the density plots and major axis regression line between each pair of datasets (m: slope, b: intercept, r: correlation coefficient).The ranges of both of the color bars are between the 1st and the 99th percentiles of the data.

Figure 3 .
Figure 3.Estimated Shapley values in the AOD region show the contribution of individual features to the model output.The colorbar shows the scaled values of explanatory variables.Red/blue color represents positive/negative effect of a feature.The higher/lower of the Shapley value of individual feature indicate a higher/lower contribution to the changes in model output.

Figure 4 .
Figure 4. Effects of spatial scale on the difference between τ SSA and τ NSSA (Method 1).The x-axis represents the increase of spatial scales from grid cell level (0.25 • ) to continental level (25 • ).The y-axis represents the absolute value of 10th, 25th, 75th and 90th relative difference (%) between τ NSSA and τ SSA.

Table 1 .
Features that are used in the ensemble regression.