Evaluating changes of biomass in global vegetation models: the role of turnover fluctuations and ENSO events

2.1.1. The global vegetation model simulations

2.1.2. Gridded observation-based annual above-ground biomass dataset.

2.1.3. Masking areas where biomass is affected by land use change

2.1.4. Multivariate ENSO Index (MEI)

2.2.1. Relative rate of change of biomass (RRB)

This study aims at corroborating biomass changes in the GVM simulations, using the AGB data as observational reference. However, the GVMs do not report AGB, but only on total biomass. Thus, a direct comparison would require first to reconstruct the AGB to below-ground biomass (BGB) relation, in order to transform AGB onto total biomass (or the other way around). Alternatively, changes in total biomass and AGB can be consistently compared in relative terms (i.e. using metrics quantifying changes per unit of biomass), provided the ratio of AGB to BGB is constant in time and space (see section S1.4 in SI). The plausibility of this assumption is supported by empirical evidence (Niklas 2005, Yang and Luo 2011, Cheng et al 2015). An additional advantage of applying a relative metric is that it enables one to compare the intensity of biomass changes across different ecosystems. Among the possible relative metrics that can be used to assess changes in biomass, we consider the relative rate of change of biomass (RRB), defined as:

$$\begin{equation} {\rm{RR}}{{\rm{B}}_{{t}}} = \frac{1}{{{{\rm{B}}_{{t}}}}}\frac{{{\rm{\Delta }}{{\rm{B}}_{{t}}}}}{{{{\Delta t}}}}. \end{equation} \tag{ 1 }$$

Here, B_t denotes the value of biomass (or AGB) at time t and ΔB_t = B_t − B_t − Δt the change occurred within the time window [t, t −δt] (hereafter Δt = 1 yr). We obtain global and regional sequences of aggregated RRB by feeding in equation (1) the time-series of the spatial sum of the gridded biomass (or AGB) values. We note that if the values of RRB are considerably small compared to Δt, the trajectory of biomass is described by the exponential

$$\begin{equation} {{\rm{B}}_{{t}}} \approx {{\rm{B}}_0}\,{\rm{exp}}\left( {{{\rm{K}}_{{t}}}} \right) \end{equation} \tag{ 2 }$$

with K_t denoting the sum ${{\rm{K}}_{{t}}} = \mathop \sum \limits_{{\rm{j}} = 1}^{\rm{m}} {\rm{RR}}{{\rm{B}}_{\rm{j}}}{{\Delta t}}$ and m the number of time steps after the beginning of the study interval, i.e. t = m Δt (for the historical period 1994–2010, m ∈ [1, 17]) (see derivation of equation (2) in S1.5). The exponential relation in equation (2) shows the pertinence of using RRB to evaluate biomass changes. According to equation (2), a systematic difference between independent estimates of K_t implies, in the long-term, an exponential deviation at the level of the corresponding biomass trajectories (as explained in section S1.5 and illustrated by figure S2 in SI).

Trends of biomass and AGB over the historical period were characterized in terms of the temporal mean of RRB, denoted by $\overline {{\rm{RRB}}}$ ̶ which in the short term (or as long as K_t ≪ 1) provides similar information as the mean percentage change of biomass with respect to the beginning of the historical period (see S1.5 in SI). The amplitude of annual variations of biomass were evaluated in terms of the temporal standard deviation of RRB, denoted as $\overline {{\rm{STD}}}$. Both, $\overline {{\rm{RRB}}}$ and $\overline {{\rm{STD}}}$ were obtained for each model and climate forcing combination. In addition, for every quantity under investigation we estimate the ensemble mean, i.e. the average over the whole set of simulations.

2.2.2. Multivariate ENSO Index-based reconstruction of biomass time series

We model the response of RRB to ENSO by assuming a heuristic linear relationship RRB_t = β + γMEI_t. After computing the intercept (β) and sensitivity (γ) parameters using least-squares regression, the reconstruction of biomass trajectories can be achieved by combining the RRB to MEI linear relation with equation (2), i.e.

$$\begin{equation} \rm B_{\it t} \approx B_0 \exp (m\beta \Delta {\it t} + \gamma \mathop \sum \limits_{j = 1}^m {\text{MEI}_j } \Delta {\it t}). \end{equation} \tag{ 3 }$$

2.2.3. The turnover rate of biomass and its relation to RRB

Biomass changes can be described by a simple balance equation (Friend et al 2014)

$$\begin{equation} \frac{{{\rm{\Delta }}{{\rm{B}}_{{t}}}}}{{{{\Delta t}}}} = {\rm{NP}}{{\rm{P}}_{{t}}} - {\rm{RL}}{{\rm{B}}_{{t}}}{{\rm{B}}_{{t}}}. \end{equation} \tag{ 4 }$$

Here, NPP_t is the average net primary production and RLB_t the rate of loss of biomass per B_t (in the models leaf and fine root turnover and mortality), during the time interval [t, t − Δt]. After dividing both sides of equation (4) by B_t, and with TRB_t defined as NPP_t/B_t one finds

$$\begin{equation} {\rm{RR}}{{\rm{B}}_{{t}}} = {\rm{TR}}{{\rm{B}}_{{t}}} - {\rm{RL}}{{\rm{B}}_{{t}}} \end{equation} \tag{ 5 }$$

which expresses the relation between the RRB_t, the equilibrium turnover rate of biomass, TRB_t, with respect to NPP input and the RLB_t (indistinctly referred in the text as the actual turnover rate). Regional scale NPP was obtained as the spatial sum of gridded values. The RLB_t was computed from the simulations, by first calculating TRB from NPP and biomass and then subtracting it from RRB, diagnosed from equation (1).

2.2.4. Relative sensitivity of RRB to NPP and RLB

2.2.5. Fraction of area affected by biomass decrease.

For a region R, we analyze the spatial share of biomass trends. For this purpose, we consider the fraction (f_A) of the area affected by a decrease of biomass as

$$\begin{equation} {{\rm{f}}_{{A}}} = \frac{{\mathop \sum \limits_{{i}} {\rm{\chi }}\left( { - \overline {{\rm{RR}}{{\rm{B}}_{{i}}}} } \right){{\rm{A}}_{{i}}}}}{{\mathop \sum \limits_{{i}} {{\rm{A}}_{{i}}}}}. \end{equation} \tag{ 6 }$$

With A_i the area of the grid-cell i in region R, and χ a step function equal to 1 if its argument is positive or 0 otherwise. Alternatively, one may consider the area fraction where $\overline {{\rm{RRB}}} \, > 0$. This would amount to estimating 1-f _A, given the virtual absence of areas where strictly $\overline {{\rm{RRB}}} = 0$.

The GVM simulations show positive significant trends of global biomass, steeper than the observed one. Whereas the observed $\overline {{\rm{RRB}}} = 0.03{\rm{\% }}$ yr⁻¹, the ensemble of simulations shows $\overline {{\rm{RRB}}} = 0.3{\rm{\% }}$ yr⁻¹ (figure 1(a)) (hereafter RRB results are reported in units of % of change per unit of biomass per year). The $\overline {{\rm{RRB}}}$ in the individual simulations range between 0.1 (in VEGAS- WFDEI) and 0.6 (in CARAIB-WFDEI). In the case of VISIT, LPJmL and CARAIB, the $\overline {{\rm{RRB}}}$ results strongly depend upon the climate forcing.

Zoom In Zoom Out Reset image size

The inter-annual variability of biomass is controlled by annual RRB fluctuations. The observations and the ensemble mean are characterized by a temporal standard deviation of RRB of $\overline {{\rm{STD}}} = 0.5$ and $\overline {{\rm{STD}}} = 0.2$, respectively. The smallest $\overline {{\rm{STD}}}$ equals 0.1 (ORCHIDEE-PGFv2) and the highest 0.9 (JULES-GSWP3) (table 3). The exceedingly large $\overline {{\rm{STD}}}$ in JULES-GSWP3 results from a single strong RRB fluctuation around 2002. Amid the models, similar values of $\overline {{\rm{STD}}}$ as in the observation are found for VEGAS-GSWP3 and CARAIB-WFDEI.

When looking at the sign changes of detrended RRB anomalies, the ensemble mean matches the observation in 11 out of 17 years. The best match between individual models and observation is found for VEGAS-GSWP3 (11 years) and the lowest one for DLEM-GSWP3 (6 years) (table 3).

Table 3. Summary statistics of the analysis of global RRB.

Observation/Model	$\overline {{\rm{RRB}}} \left( {{{10}^{ - 3}}{\rm{y}}{{\rm{r}}^{ - 1}}} \right)$	$\overline {{\rm{STD}}} \left( {{{10}^{ - 3}}{\rm{y}}{{\rm{r}}^{ - 1}}} \right)$	Fraction of years where the sign of detrended RRB agrees with the sign of RRB derived from VOD
VOD	0.4	4.6
VEGAS
GSWP3	1.3	4.5	0.7
WFDEI	1.2	4.3	0.5
VISIT
PGFv2	3.0	2.0	0.6
GSWP3	1.4	2.1	0.5
WFDEI	2.8	2.3	0.6
DLEM
PGFv2	2.1	2.0	0.4
GSWP3	2.1	2.1	0.4
WFDEI	2.1	2.3	0.4
LPJ-GUESS
PGFv2	3.5	3.0	0.5
GSWP3	3.6	3.3	0.5
WFDEI	4.0	3.3	0.6
ORCHIDEE
PGFv2	3.8	1.0	0.6
GSWP3	4.1	1.3	0.4
WFDEI	4.1	2.4	0.5
LPJmL
PGFv2	3.1	3.1	0.6
GSWP3	5.6	2.7	0.6
WFDEI	3.7	2.9	0.7
JULES-UoE
PGFv2	4.2	1.3	0.5
GSWP3	4.5	9.4	0.5
CARAIB
PGFv2	3.6	3.0	0.4
GSWP3	4.3	5.0	0.5
WFDEI	6.3	3.7	0.5
ENSEMBLE	3.4	1.8	0.7

The results of the regression of global RRB with the Multivariate ENSO Index (MEI) (equation (3)) show the importance of ENSO events in explaining the inter-annual variability of global biomass (figure 1(d)). The observed RRB to MEI regression shows a R² value of 0.5 (figure 1(b)), as opposed to a model ensemble mean R² of 0.2. The highest value R² = 0.5 in the simulations corresponds to ORCHIDEE-WFDEI (table 4).

In the observations, the sensitivity of RRB to MEI (figure 1(c)) (denoted by γ and hereupon reported in the text in units of % yr⁻¹ per MEI unit change) is negative and equals −0.32 ± 0.08. The ensemble mean does not show a significant sensitivity (at p-value ≤ 0.05). However, significant negative sensitivities were identified in VISIT-GSWP3 where γ = −0.12 ± 0.1, as well as in all the simulations of ORCHIDEE (from −0.05 ± 0.02 in PGFv2 to −0.17 ± 0.05 in GSWP3) and CARAIB (from −1.8 ± 0.7 in PGFv2 to −2.8 ± 1.1 in GSWP3). The intercept regression coefficients in the simulations can be practically assimilated to $\overline {{\rm{RRB}}}$ (see figure S3 in SI).

After restricting the analysis to the Northern Hemisphere extra-tropical region (23°N–85°N), the observed RRB to MEI regression shows a negative sensitivity with γ = −0.8 ± 0.2 and R² = 0.4 (see figure S4 in SI). In this region, the sensitivities of RRB in the GVM simulations are not significant (see figure S4 in SI). The opposite occurs in the Tropical region (23°S—23°N), where the sensitivity of RRB is not significant in the observation (see figure S5 in SI), but in the simulations of VISIT-GSWP3 (−0.2 ± 0.1), JULES-PGFv2 (−0.1 ± 0.03), CARAIB (from −0.3 ± 0.1 in WFDEI to −0.4 ± 0.1 in GSWP3) and of ORCHIDEE (from −0.05 ± 0.02 in PGFv2 to −0.2 ± 0.06 in GSWP3). The simulation CARAIB-PGFv2 shows the highest R² (0.5).

To gain understanding on the discrepancies in global biomass change, between the GVM simulations and the observation (figure 1(a)), the RRB analysis was applied to regions of natural forest with high biomass densities, namely the tropical forests of Amazonia (AMA) and Congo (CON) and the boreal forests of Western Siberia (WES), Central-Siberia (CES) and North-Western America (NWA)–the latter including also temperate forest (figure 2 and table S1 in SI). Altogether, the contributions of these regions sum up ca. 40% of the total AGB considered in the global analysis.

Zoom In Zoom Out Reset image size

The observation shows a major increase of biomass in WES (19% growth of AGB with respect to the beginning of the study period), that is characterized by $\overline {{\rm{RRB}}} = 0.8$ (figure 3(a)). This increment largely surpasses in magnitude the loss of biomass observed in NWA ($\overline {{\rm{RRB}}} = - 0.08$), CES ($\overline {{\rm{RRB}}} = - 0.07$), AMA ($\overline {{\rm{RRB}}} = - 0.07$) and CON ($\overline {{\rm{RRB}}} = - 0.01$). With the exception of CON, the decremental trends of biomass associated to the observed regional values of $\overline {{\rm{RRB}}}$ are significant. In contrast, the $\overline {{\rm{RRB}}}$ shown by the ensemble mean of the simulations is positive in all of the regions. According to the ensemble mean, the most intense growth of biomass occurs in CES ($\overline {{\rm{RRB}}} = 0.8$), followed by WES ($\overline {{\rm{RRB}}} = 0.6$), NWA ($\overline {{\rm{RRB}}} = 0.4$), CON ($\overline {{\rm{RRB}}} = 0.2$) and AMA ($\overline {{\rm{RRB}}} = 0.2$). The overestimation of biomass growth in the simulations is also expressed as an underestimation of the fraction f_A of the total area where $\overline {{\rm{RRB}}} < 0$ (see section 2.2.5). The least underestimation of this fraction in the boreal regions is found in NWA (ensemble mean: 32.4%; observation: 56.3%) and in the tropical regions in CON (ensemble mean: 27.3%; observation: 61.7%) (see summary in tables S2-S6 in SI).

Zoom In Zoom Out Reset image size

The individual simulations show positive values of $\overline {{\rm{RRB}}}$ in all of the regions, excerpt for LPJmL-PGFv2, VEGAS (forced by GSWP3 and WFDEI) and VISIT (GSWP3 and WFDEI) in CON (figure 3). Practically all of the trends of biomass associated to the $\overline {{\rm{RRB}}}$ values in the simulations are significant across the regions. The simulated increments of biomass show their best mutual agreement in AMA, with $\overline {{\rm{RRB}}}$ ranging between 0.06 (VEGAS-WFDEI) and 0.3 (JULES-GSWP3). On the contrary, the largest inter-simulation discrepancies in $\overline {{\rm{RRB}}}$ take place in CON, with values ranging from −0.4 (VISIT-WFDEI) to 1.4 (LPJmL-WFDEI). In the boreal regions, the largest spread of $\overline {{\rm{RRB}}}$ occurs in CES, with values ranging between 0.1 (VEGAS-GSWP3) and 1.4 (LPJmL-GSWP2). In WES, $\overline {{\rm{RRB}}}$ ranges from 0.3 (CARAIB-WFDEI) to 1.4 (ORCHIDEE-GSWP3) and in NWA from 0.08 (DLEM-PGFv2) to 0.9 (CARAIB-GSWP3) (see summary of $\overline {{\rm{RRB}}}$ results in tables S2–S6 in SI).

In those simulations showing the largest values of regional $\overline {{\rm{RRB}}}$, RRB is positive in at least 14 out of the 17 years covered by the study period. We remind that positive values of RRB occur whenever TRB, i.e. the ratio NPP/B, exceeds RLB (equation (5)). The simulations show strong TRB to NPP correlations (with a R² above 0.7 (0.8) in 100 (86) of the 110 cases comprised by the simulations across the regions). In all the regions the significant trends of NPP are positive (figure 4) (NPP trends are hereafter denoted as TNPP and reported in units of Tg C yr⁻²; significant refers to p ≤ 0.05)). Most of the simulations show significant TNPPs in CES, with slopes ranging from 2.2 ± 0.8 (CARAIB-GSWP3) to 10.9 ± 2.2 (ORCHIDEE-WFDEI). In WES and NWA, significant TNPPs are shown only by the simulations of ORCHIDEE. In WES, these trends range from 1.7 ± 0.5 (GSWP3) to 2.6 ± 0.5 (WFDEI) and in NWA from 1.5 ± 0.7 (GSWP3) to 1.8 ± 0.6 (WFDEI). In the tropical regions, significant TNPPs occur mostly in CON, in all of the simulations forced by PGFv2, in DLEM, LPJmL, ORCHIDEE and LPJ-GUESS forced by GSWP3, as well as in DLEM-WFDEI. The lowest significant TNPP in CON equals 1.2 ± 0.5 (ORCHIDEE-GSWP3) and the highest one 6.7 ± 1.3 ( LPJmL-PGFv2). In AMA, significant TNPPs are displayed only by DLEM, in the range between 2.2 ± 0.9 (PGFv2) and 3.6 ± 1.2 (WFDEI).

Zoom In Zoom Out Reset image size

Hence, only in the case of ORCHIDEE-GSWP3 in WES the highest regional $\overline {{\rm{RRB}}}$ coincides with a significant positive TNPP. In CES, the highest $\overline {{\rm{RRB}}}$ is associated with a negative significant trend in RLB (−0.16% yr⁻², in LPJmL-GSWP3) (figure 4) (RLB trends are subsequently denoted as TRLB and reported in units of % yr⁻²). Similarly in CON, the highest $\overline {{\rm{RRB}}}$ is associated to a negative significant TRLB (−0.095, in LPJmL-WFDEI). Other cases where biomass growth appears reinforced by negative significant TRLBs are: in WES, LPJmL-WFDEI (−0.05); in NWA, ORCHIDEE-WFDEI (−0.1) and the simulations of LPJmL (−0.08) and VISIT (−0.05) forced by PGFv2; in CON, in LPJ-GUESS-PGFv2 (−0.026), LPJ-GUESS-WFDEI (−0.03) and CARAIB-GSWP3 (−0.08). On the other hand, the increment of biomass is attenuated by positive significant TRLBs in CES, in all the simulations of VISIT (0.06–0.07) and ORCHIDEE (0.2–0.3); in WES, in VISIT-WFDEI (0.03); in CON, in VISIT-WFDEI (0.02) and in all the simulations of DLEM (0.004–0.006). Finally, in AMA, the positive significant TRLBs occur in DLEM (0.003, with forcings GSWP3 and WFDEI) and in LPJmL-PGFv2 (0.09) (summary of regional trends of NPP and RLB in tables S2–S6 in SI).

In the absence of significant TNPPs and TRLBs, as is the case of most of the simulations in NWA, differences in $\overline {{\rm{RRB}}}$ could be simply interpreted as discrepancies at the level of the background simulated values of TRB (NPP) and RLB during the historical period. However, important changes in biomass also occur during shorter subintervals. The negative $\overline {{\rm{RRB}}}$ values (figure 3) in CON, in the simulations of VEGAS, are mainly due to an overall decrease in TRB values within the periods 1995–1998 and 2003–2005, that co-occur along with important increments in RLB (see plots of regional trajectories of TRB and RLB in figures S6-S10 in SI). Similarly, the decrease of biomass described by LPJmL-PGFv2 in CON is mainly due to a sharp rise in RLB between 1994 and 1996, after which RLB remains above TRB until 2002. On the other hand, we notice that the high standard deviation of global RRB in JULES-GSWP3 (cf section 3.1) is related to large RLB variations occurring in CON, AMA and NWA, between years 2001 and 2004. Overall, in each region the simulations generated by different models and forcings show dissimilar trajectories of RLB and TRB. In spite of this, we notice that in AMA all the simulations describe a decrease in TRB between 1996–1998, followed by recovery until 1999.

These findings point to important differences regarding the relative influence of NPP (via TRB) and RLB on RRB (section 2.2.4). Altogether, in most of the simulations the annual variations in regional RRB are mostly influenced by NPP (see figure S11 in SI). Particular cases where the influence of RLB on RRB is dominant are: in WES, CES and NWA, the simulations of VEGAS, LPJ-GUESS and CARAIB; in AMA, in the simulations of VEGAS, LPJ-GUESS (forced by WFDEI and GSWP3), CARAIB (GSWP3 and PGFv2), JULES-GSWP3 and LPJmL-PGFv2 and, in CON, LPJ-GUESS, LPJmL-WFDEI and CARAIB-GSWP3. Finally, the influences of NPP and RLB are closely comparable in the simulations LPJ-GUESS-WFDEI in CES, CARAIB-PGFv2 in WES, JULES-GSWP3 in CON and CARAIB-WFDEI in AMA.

The main objective of this work is to analyze the ability of the GVMs to reproduce historical changes in biomass. Evaluating the physiological component of biomass changes is important as a preliminary step towards assessing the model representations of direct anthropogenic influences. Consequently, this study focused on natural terrestrial ecosystems.

We used the AGB maps provided in Liu et al (2015) as observational reference to verify the trends and inter-annual variability of biomass. Other similar datasets are available (cf Saatchi et al 2011), but the Liu et al (2015) dataset is the only product providing multi-year coverage. Since the GVMs do not report explicitly AGB but only total biomass, their comparison with observations can be achieved based on a relative metric (e.g. percentage change of biomass with respect to a baseline value). The underlying assumption of this approach is that the AGB to BGB ratio is constant in space and time (see explanation in S1.3 in SI). A linear AGB to BGB relation has been shown to be a valid approximation for non-woody plants (Niklas 2005), forest ecosystems (Yang and Luo 2011), as well as for aggregate biomass estimates involving different plant communities (Cheng et al 2015). The use of a relative metric has the additional advantage of enabling the comparison of the intensity of changes in biomass across different ecosystems.

Moreover, we choose the RRB metric because: (i) biomass trajectories are determined by the exponential of the cumulative sum of RRB over time (see equation (2) in the text and section S1.4 in SI). This is relevant since a persistent bias of RRB in the models, with respect to observations, may imply a non-linear increase in the error of biomass trajectories in long term projections (see illustration in figure S2 in SI); (ii) RRB equals the deviation of the actual turnover rate of biomass RLB from its equilibrium value TRB (equation (5)). Gaining knowledge on the dynamics and climate sensitivity of the turnover rate of biomass is critical to understand possible future trajectories for the terrestrial carbon cycle (Bloom et al 2016). For instance, the behavior of RLB has been recognized as a substantial source of uncertainty in multi-model projections of biomass in a similar set of GVMs (Friend et al 2014). Moreover, a recent analysis of boreal forests by Thurner et al (2017) has shown that an underestimation of the average of TRB leads to important discrepancies in biomass between observation and simulations.

Recently, Luo et al (2017) introduced a new framework to analyze the transient dynamics of carbon pools in ecosystems from a systems perspective. In their formulation, changes in the ecosystem carbon storage are described by the carbon storage potential: the difference between the actual carbon storage and the carbon storage capacity (Olson 1963). The latter being defined as the maximum amount of carbon that an ecosystem can store, given the environmental condition at a point in time. Under stationary environmental conditions, the transient dynamics of the carbon pools is determined by the decrease of the carbon storage potential towards zero. The RRB bears similarity with the carbon storage potential concept. In absence of trends in environmental conditions, the biomass pool attains stationarity through the gradual balancing of RLB and TRB, occurring via the uptake and the loss of carbon.

Different factors may affect the reliability of the approach followed in this study to evaluate historical changes of biomass in the selected GVMs. For instance, intense drought episodes may compromise the consistency of the RRB comparison, between observed AGB and simulated total biomass, due to changes in carbon allocation patterns in vegetation (Doughty et al 2015). However, the AGB product in Liu et al (2015) also embodies significant uncertainties, that are partly related to the spatial extrapolation method (Mitchard et al 2014) based on Saatchi et al (2011). Moreover, this AGB dataset may be of limited use for evaluating impacts of climate extremes on tropical forests, in view of the possible underestimation of inter-annual variations of biomass by VOD signals, in comparison with other satellite indicators of vegetation (Liu et al 2011). In addition, we notice that the AGB time-series indicate a decrease of biomass in arid regions of Australia (figure 2), which still has to be confirmed by other techniques. It is also important to bear in mind the limitations of the land-mask applied to disentangle the anthropogenic and physiologic components of RRB. The land cover selection criterion applied for the regional analysis was especially strict, in the sense that only grid cells with a classification of 100% natural vegetation were considered (see comparison of surface areas in table S1 in SI). However, the land-cover classification was based on the available information from year 2005 and important changes in land-cover have occurred ever since Houghton et al (2012). Finally, the construction of the mask did not accounted for information on other human activities, such as wood harvest in managed and unmanaged forests, regrowth in abandoned agricultural areas (Pan et al 2011), or landscape fragmentation (Pütz et al 2014).

The global observed RRB indicates an increase in global biomass between 1994 and 2010 (figure 1(a)). As reported by Zhu et al (2016), recent increments of global biomass have been attributed to positive plant-physiological effects of increasing atmospheric CO₂, as well as to the prolongation of growing seasons in the Northern Hemisphere, as a result of global warming. The global increase of biomass described by the GVMs appears too optimistic, with the ensemble mean of the simulations showing a global $\overline {{\rm{RRB}}}$ one order of magnitude larger than the observed one (figure 1(a)).

The sensitivity of the global RRB variations to ENSO in the observation (figure 1(c)) is practically explained by the response of AGB in the Northern Hemisphere extra-tropic (see figure S4 in SI). The apparent absence of a significant ENSO sensitivity in the observation of RRB in the tropics (see figure S5 in SI) is striking, as it contrasts with evidence on the strong influence of El Niño events, on extreme heat and drought conditions, over the Amazonia forests (Jiménez-Muñoz et al 2016). This ambiguity suggests a limited capacity of VOD signals to adequately capture inter-annual variations of AGB in tropical forests. On the contrary, the simulations of 4 out of the 7 GVMs show significant sensitivities of RRB in the tropics (see figure S5 in SI), but none of them in the Northern Hemisphere (see figure S4 in SI). The underlying differences in RRB variations in the simulations are discussed in further detail below in section 4.3.

The regional value of $\overline {{\rm{RRB}}}$ in WES is higher in the observation than in most of the simulations (figure 3). In the rest of the boreal regions analyzed, the sign of $\overline {{\rm{RRB}}}$ is negative in the observation, and positive in all of the simulations. A similar situation is found for AMA. The analysis of $\overline {{\rm{RRB}}}$ in CON shows no significant change of biomass in the observation, but increments in most of the simulations. These findings are relevant for the purpose of understanding the difference between the observed and simulated trends of global biomass (figure 1(a)), given the broad extent and high concentration of biomass in these regions (whose contributions sum up to ca. 40% of the global total AGB in this analysis).

The overestimation of the regional $\overline {{\rm{RRB}}}$ in the simulations may be related to the CO₂ sensitivity of carbon uptake. All of the regional significant TNPPs simulated by the GVMs are positive (figure 4). The significant TNPPs in the boreal regions are mostly localized in CES, whereas in the tropics these mainly occur in CON. A recent analysis of the CMIP5 earth system models, forced only by the increase of CO₂ (Smith et al 2016) over the period 1982–2011, has pointed to an overestimation of incremental NPP trends. The excessive CO₂ response of NPP was attributed to the lack of nutrient constraints in most of the CMIP5 models. This possibility could be debated in the case of the GVMs used here, where the degree of overestimation of trends of gross primary productivity (GPP) depends on the observational reference that is used as benchmark (Ito et al 2017). Although the presence of an excessive CO₂ fertilization effect in all of the GVM simulations cannot be excluded, its role in explaining differences in $\overline {{\rm{RRB}}}$ is not conclusive. Nitrogen limitation is included in VEGAS, DLEM and LPJ-GUESS (table 1) and the significant TNPPs displayed by VEGAS and LPJ-GUESS in CES are below the ensemble mean (figure 4). However, the $\overline {{\rm{RRB}}}$ values shown by the LPJ-GUESS simulations in this region (figure 3) are comparably high as those shown by the models where an excess in CO₂ fertilization could be expected (due to their lacking representation of nitrogen limitation). Overall, we find that the single case where the largest regional value of $\overline {{\rm{RRB}}}$ coincides with a positive significant TNPP corresponds to the simulation ORCHIDEE-GSWP3 in WES (see figures 3 and 4). Indeed, regardless how disproportionate the CO₂ response of simulated NPP may be, we recall that an increase of NPP at a point in time will lead to a greater increment of biomass provided it raises further the value of TRB above RLB (equation (5)). In other words, understanding discrepancies in biomass growth necessarily requires to assess NPP (or more specifically TRB) and RLB simultaneously.

In the Boreal regions, the variations of RRB are mainly influenced by RLB in the simulations of VEGAS, LPJ-GUESS and CARAIB, and are dominated by NPP (via TRB) in the rest of the models (see figure S11 in SI). The simulations do not reproduce the observed significant ENSO sensitivity of RRB in the Northern Hemisphere extra-tropics. The analysis of boreal regions by Thurner et al (2017) pointed out the failure of a very similar set of GVMs at reproducing phenomenological relations between climate variables and average TRB. Amidst the physiological processes involved in the turnover of biomass, mortality is highly complex, not thoroughly understood and, consequently, insufficiently described in global vegetation models (Steinkamp and Hickler 2015). In addition to the responses of TRB and RLB to climate and CO₂ concentrations, disturbances can constitute an important source of differences in biomass change. During the historical period covered by the RRB analysis, decrements of biomass occurred in NWA as a result of widespread and severe pine beetle outbreaks (Kurz et al 2008). The positive $\overline {{\rm{RRB}}}$ values shown by the simulations in NWA may thus be attributed to the lack of representation of insect infestation effects in the GVMs. In CES, incremental trends of fire intensity and frequency were reported also for the analyzed period (Ponomarev et al 2016, Schaphoff et al 2016). Consequently, significant increments of RLB in CES could be expected in the simulations of those models featuring fires, i.e. CARAIB, LPJmL, LPJ-GUESS, VISIT and VEGAS (table 2). Among them, the simulations of CARAIB show a high fraction of the total area of CES with $\overline {{\rm{RRB}}} < 0$, as in the observation (see table S5 in SI). However, significant positive TRLBs in CES are shown only by VISIT and ORCHIDEE (figure 4). Given the time average of TRB displayed by VISIT and ORCHIDEE in CES, an increase in the time average of RLB by at most 3.8% and 6.5%, respectively, would be required to match the observed $\overline {{\rm{RRB}}}$ in this region. It is important to mention that limitations of fire models in reproducing observed trends in burned area have been identified in a recent GVM intercomparison (Andela et al 2017). In WES, the drastic increase of the observed carbon sink can be partly attributed to the expansion of forested areas, following agricultural abandonment and reduced harvesting (Zhu et al 2016). The generally lower $\overline {{\rm{RRB}}}$ values shown by the simulations in this region, compared to the observation, could be explained by the lack of representation of these processes in the models (or equivalently, by the failure of the land-mask applied to exclude their influence from the RRB analysis). Nonetheless, the case of WES illustrates the importance of analyzing the spatial aspects of the NPP and RLB trends, as the observation also shows a decline of biomass over ca. 32% of the total area (see table S3) –occurring mainly in the North-Eastern part of this region (figure 2). We find that comparable area fractions in WES, as in the observation, are displayed by LPJmL-PGFv2 (32.7%) and CARAIB-GSWP3 (34.6%), whereas those shown in the rest of the simulations are lower. Overall, these findings suggest that further analysis of differences in average TRB (NPP) and RLB, as well as of the implementation of fire dynamics and disturbances in the GVMs is required.

In the case of tropical forests, the approach followed in this study may not be reliable to verify the behavior of RRB in the GVMs. The positive sign of $\overline {{\rm{RRB}}}$ in the simulations of AMA is in qualitative agreement with previous observations of biomass change in intact forest (Pan et al 2011). However, recent analyses of in-situ data of Amazonia have shown a significant positive trend in mortality during the last decades (Brienen et al 2015) and subsequent biomass loss and carbon release after droughts in 2005 and 2010 (Lewis et al 2011). The role of drought as a driver of tree mortality is corroborated by experimental results in this region (Rowland et al 2015, Feldpausch et al 2016). Yet the analysis of Brienen et al (2015) asserts that, in spite of the impact of severe droughts, mature forests in Amazonia remained accumulating biomass during the 1983–2010 period. In this view, in addition to droughts, the negative sign of the AGB derivation of $\overline {{\rm{RRB}}}$ in AMA can be related to the impacts of human activities, such as clearing, landscape fragmentation and selective logging (Morton et al 2011, Houghton et al 2012) –not being excluded by the land-mask. A similar situation holds for the CON region, which has also experienced degradation (Zhuravleva et al 2013). Moreover, the presence of intense droughts in Amazonia may cast doubts on the validity of the AGB-to-BGB isometry over time, due to changes in carbon allocation patterns in vegetation (Doughty et al 2015). Moreover, the AGB product in Liu et al (2015) may underestimate the effects of climate extremes related to ENSO on biomass variations in tropical forests.

The GVM simulations show the best mutual agreement in $\overline {{\rm{RRB}}}$ in AMA, compared to the rest of the regions. On the contrary, the largest spread of simulated $\overline {{\rm{RRB}}}$ occurs in CON. The positive and negative extreme values of $\overline {{\rm{RRB}}}$ in CON are associated, respectively, with negative and positive significant TRLBs. In relation to the ENSO sensitivities of RRB exhibited by some of the simulations in the tropics (see figure S5 in SI), the analysis of AMA and CON shows that RRB variations are mostly influenced by NPP in JULES-PGFv2, VISIT-GSWP3, ORCHIDEE-GSWP3; whereas these are dominated by RLB in CARAIB-GSWP3. Altogether, in models where the description of mortality and/or leaf turnover is directly modulated by climate (see table 2), the degree of influence of RLB on RRB varies across regions and climate forcings (see figure S11 in SI). Conversely, in VISIT and DLEM, where the turnover of biomass is simply modulated by one factor (fire and tree age, respectively), the variations of RRB are consistently dominated by NPP. In spite of these differences, in AMA all of the simulations show a drop of NPP between 1996 and 1998, that leads to a decrease in TRB (see figure S7 in SI). This can be interpreted as a response to the extreme El Niño event that occurred during this period. For further information on the response of the carbon uptake to climate variations in the GVMs analyzed here we refer to Ito et al (2017), Chang et al (2017). Finally, the RRB results at the global (figure 1(a)) and regional (figure 2) scales are in general sensitive to the climate forcing dataset. This can be partly attributed to differences in the carbon uptake flux induced by the discrepancy in solar radiation among the forcing datasets (Ito et al 2017).