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

This paper evaluates the ability of eight global vegetation models to reproduce recent trends and inter-annual variability of biomass in natural terrestrial ecosystems. For the purpose of this evaluation, the simulated trajectories of biomass are expressed in terms of the relative rate of change in biomass (RRB), defined as the deviation of the actual rate of biomass turnover from its equilibrium counterpart. Cumulative changes in RRB explain long-term changes in biomass pools. RRB simulated by the global vegetation models is compared with its observational equivalent, derived from vegetation optical depth reconstructions of above-ground biomass (AGB) over the period 1993–2010. According to the RRB analysis, the rate of global biomass growth described by the ensemble of simulations substantially exceeds the observation. The observed fluctuations of global RRB are significantly correlated with El Niño Southern Oscillation events (ENSO), but only some of the simulations reproduce this correlation. However, the ENSO sensitivity of RRB in the tropics is not significant in the observation, while it is in some of the simulations. This mismatch points to an important limitation of the observed AGB reconstruction to capture biomass variations in tropical forests. Important discrepancies in RRB were also identified at the regional scale, in the tropical forests of Amazonia and Central Africa, as well as in the boreal forests of north-western America, western and central Siberia. In each of these regions, the RRBs derived from the simulations were analyzed in connection with underlying differences in net primary productivity and biomass turnover rate ̶as a basis for exploring in how far differences in simulated changes in biomass are attributed to the response of the carbon uptake to CO2 increments, as well as to the model representation of factors affecting the rates of mortality and turnover of foliage and roots. Overall, our findings stress the usefulness of using RRB to evaluate complex vegetation models and highlight the importance of conducting further evaluations of both the actual rate of biomass turnover and its equilibrium counterpart, with special focus on their background values and sources of variation. In turn, this task would require the availability of more accurate multi-year observational data of biomass and net primary productivity for natural ecosystems, as well as detailed and updated information on land-cover classification.


Introduction
Large uncertainties of the global carbon budget (Le Quéré et al 2015) are associated with the terrestrial carbon cycle (Bloom et al 2016). Reducing those uncertainties requires an improved capacity to estimate the fluxes of carbon between the atmosphere and terrestrial ecosystems, along with gaining better understanding of the dynamics of the corresponding carbon pools and their sensitivity to climate. For this purpose, global vegetation models (GVMs) have been developed to describe the effects of climate, atmospheric CO 2 and other drivers on a large number of ecological processes, including photosynthesis, respiration, carbon allocation, phenology, and plant mortality which altogether control the dynamics of carbon pools. The possibility to evaluate the performance of such complex vegetation models is currently improving thanks to novel satellite-based observation methods.
Recently, the ensemble of GVMs available through the Inter-Sectoral Model Intercomparison Project Phase 2a (ISIMIP2a) has been evaluated with regards to the models' ability to reproduce the historical trends and variability of the gross primary production and net biome productivity (Ito et al 2017, Chang et al 2017, Chen et al 2017. Using the same ensemble of GVMs, this paper aims at corroborating recent trends and interannual variations of biomass stocks. Beyond comparing changes in biomass, derived from GVMs and observations, we propose to address the processes explaining these changes. As recently shown, the equilibrium turnover rate of biomass (TRB) Thurner et al (2017) and the rate of loss of biomass (RLB) (Friend et al 2014, Bloom et al 2016 are major sources of uncertainty in model representations of biomass dynamics. The RLB is the inverse of the residence time of carbon in vegetation, i.e. the time carbon remains in living biomass pools from fixation to removal. In natural ecosystems, biomass removal occurs through background mortality of stems, foliage and roots, and through losses from disturbances such as fire, storm or insect damage. The TRB is defined as the value of RLB when the biomass pool is in equilibrium, i.e. the inverse of the equilibrium residence time. The objective of this paper is to analyze changes in biomass stocks from the view point of the deviations of the biomass loss flux, from the loss in the virtual equilibrium state, i.e. the difference between TRB and RLB. This difference, here referred as the relative rate of change in biomass (RRB), is relevant for model evaluation since a positive/negative value in the temporal mean of RRB entails that the biomass system tends to move towards its equilibrium following an exponential increase/decrease in biomass. Observed RRB is diagnosed from annual maps of biomass retrieved from vegetation optical depth satellite observations (Liu et al 2015) during the period 1993-2010. We note that this product may be less accurate than other biomass datasets, but it is the only temporally resolved one allowing for comparison with model outputs on annual time-scale.
We focus on the variations of processes controlling biomass change in so called 'natural' terrestrial ecosystems, i.e. in ecosystems not being affected by land use change. For this purpose, areas covered with cropland and pastures were filtered out from the gridded data of the observation and GVM simulations, prior to the analysis of RRB. At the global scale, this work evaluates the performance of GVMs to reproduce the observed temporal mean of RRB. In addition, global biomass variations are evaluated with focus on the relation between RRB fluctuations and the El Niño Southern Oscillations (ENSO) that occurred since 1993.
The behavior of RRB in the models is further analyzed in specific regions representative of the main forest biomes (table S1 in supplementary information (SI) available at stacks.iop.org/ERL/13/075002/mmedia). The underlying differences in RRB in the GVM simulations are investigated in relation with the trends and variations of the net primary productivity (NPP) and RLB. In addition, the spatial share of regional biomass trends in the simulations is evaluated in terms of the land cover fraction where the temporal mean of RRB exhibits negative values -i.e. areas where the observation indicates loss of biomass.

The global vegetation model simulations
The simulations were performed by eight GVMs, namely, CARAIB (Warnant et al 1994, Gérard et al 1999, Dury et al 2011, DLEM (Tian et al 2015), JULES-UoE (Best et al 2011, Clark et al 2011, Harper et al 2016, LPJ-GUESS (Smith et al 2014), LPJmL (Bondeau et al 2007), ORCHIDEE (Krinner et al 2005), VEGAS (Zeng et al 2005), and VISIT (Ito and Inatomi 2012). These models embody different degrees of detail in the representation of biomass dynamics reflecting, productivity, allocation, maintenance respiration, vegetation dynamics and losses from mortality and disturbances (see tables 1 and 2). Each simulation was forced by a time-series of global atmospheric CO 2 concentration (Keeling and Whorf 2005) and historical reconstructions of climate on a 0.5 • global grid in daily temporal resolution. In addition, simulations of all models, except CARAIB, were constrained with historical land-use maps. The landuse dataset describes cropland and pasture areas for the year 2000 (Ramankutty et al 2008). Changes in cropland and pasture covers were extrapolated into past and future using the trends reported in HYDE3 (Goldewijk et al 2011). A further detailed classification of crops and irrigated land was achieved using the MIRCA2000 crop dataset (Portmann et al 2010). Table 1. Summary of key processes represented in the global vegetation models used in this study: dynamical vegetation (DYNV), nitrogen limitation (NITL); CO 2 fertilization effect (CO2E); light interception (LIGI); light use efficiency (LIGU); phenology (PHEN); number of plant functional types in natural vegetation (#PFTS); water stress on photosynthesis (WASP); heat stress on photosynthesis (HEAP); evapo-transpiration (EVAT); differences in root depth (DROD); root distribution over depth (RDOD); closed energy balance (CBAL); latent heat (LATH); sensible heat (SENH); information non-specified (N/S) (for further details see ISIMIP2a, https://esg.pik-potsdam.de/projects/isimip2a/).

MODELS DYNV NITL CO2E LIGI LIGU PHEN #PFTS WASP HEAP EVAT DROD RDOD CBAL LATH SENH
--All simulations were accomplished following a two-step process: (1) a 'spin-up' run in order to ensure that the carbon pools in the models reach an equilibrium, according to the climate, CO 2 and land-use conditions at the beginning of the 20th century ; (2) starting from year 1901, models were forced away from equilibrium by historical climate, CO 2 concentrations and land-use data, up to 2010 (www.isimip.org/gettingstarted/#simulationprotocol). The model output used in this work covers the period 1971-2010.
One original feature of this study is the application of three different climate forcing datasets, instead of only one forcing usually used in other model inter-comparison projects (Sitch et al 2015, Huntzinger et al 2013. We consider the following climate forcing datasets: PGFv2 (http://hydrology. princeton.edu/data.pgf.php) (not available for VEGAS), GSWP3 (http://hydro.iis.u-tokyo.ac.jp/ GSWP3) and WFDEI (combined with WATCH) (www.waterandclimatechange.eu/about/watch-forcin g-data-20th-century) (not available for JULES); The data pre-processing is described in section S1.1 in SI.

Gridded observation-based annual aboveground biomass dataset.
The above-ground biomass (AGB) inter-annual dataset used for model evaluation was produced by Liu et al (2015) by regressing satellite observations of vegetation optical depth (VOD) with the map of tropical AGB from Saatchi et al (2011), VOD is an indicator of the water content of woody and leaf vegetation tissues, related to biomass. The VOD-AGB regression coefficients were used to extrapolate AGB on a grid of 0.25 • for each year during the period 1993-2012, based on annual maps of VOD. We re-gridded the Liu et al (2015) biomass dataset to the 0.5 • resolution of the GVMs, using a conservative remapping algorithm (see details in S1.2 in SI).

Masking areas where biomass is affected by land use change
Since changes in land-use and land-cover drastically affect the turnover rate of carbon in vegetation (Erb et al 2016) and thereby biomass change estimates, we focus on so-called 'natural' ecosystems (those not affected by land use). For this purpose, a land-mask was applied to both the AGB and GVM output to exclude: 1) gridcells where the cover fraction occupied by pastures and cropland in 2005 (the last available time-step in the land-use pattern used in the simulations) exceeds 5%; 2) grid-cells where VOD is known to be affected by the presence of open water bodies, snow and ice cover (see figure S1 in SI). Excluding the effect of land use changes is especially critical when comparing GVMs with observations at the regional scale. Therefore, for the regional scale analysis we consider only grid cell values where the land-use pattern indicates a natural vegetation cover of 100% (see table S1 in SI). Finally, it is important to consider that managed forests have a faster turnover rate and a lower biomass than primary forests simulated by the GVMs. However, our mask does not exclude forested areas under management even though forest management is not represented in the GVMs. This constitutes a source of systematic error in the model evaluation, as the effect of management should be captured in the AGB observations.

Multivariate ENSO Index (MEI)
The time-series of multivariate El Niño Southern Oscillation index (MEI) has been derived on the basis of six observed climate variables (Wolter and Timlin 2011) over the 1993-2010 period of interest: sea level pressure, sea surface temperature, zonal and meridional components of the surface wind, surface air temperature and cloudiness (www.esrl.noaa.gov/psd/ enso/mei/).

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: Here, B denotes the value of biomass (or AGB) at time t and ΔB = B − B −Δ 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 Δ , the trajectory of biomass is described by the exponential  (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 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 RRB -which in the short term (or as long as K ≪ 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 STD. Both, RRB and 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.

Multivariate ENSO Index-based reconstruction of biomass time series
We model the response of RRB to ENSO by assuming a heuristic linear relationship RRB = + MEI . 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. (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) Here, NPP is the average net primary production and RLB the rate of loss of biomass per B (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 , and with TRB defined as NPP ∕B one finds which expresses the relation between the RRB , the equilibrium turnover rate of biomass, TRB , with respect to NPP input and the RLB (indistinctly referred in the text as the actual turnover rate). Regional scale NPP was obtained as the spatial sum of gridded values. The RLB was computed from the simulations, by first calculating TRB from NPP and biomass and then subtracting it from RRB, diagnosed from equation (1).

Relative sensitivity of RRB to NPP and RLB
The relative influences of NPP and RLB on RRB were evaluated in terms of the coefficients NPP and , obtained from the least-squares regression of the linear equationRRB = NPPÑ PP + RLBR LB + c, with symbol ∼ denoting standardized variables (i.e. for each variable, subtracting the temporal mean over the historical period and then dividing by the corresponding standard deviation). Here, NPP > RLB ( NPP < RLB ) indicates that the variations of RRB are mostly influenced by NPP (RLB).

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 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 RRB > 0. This would amount to estimating 1-f A , given the virtual absence of areas where strictly RRB = 0.

Global RRB and its sensitivity to ENSO
The GVM simulations show positive significant trends of global biomass, steeper than the observed one. Whereas the observed RRB = 0.03% yr −1 , the ensemble of simulations shows RRB = 0.3%yr −1 (figure 1(a)) (hereafter RRB results are reported in units of % of change per unit of biomass per year). The 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 RRB results strongly depend upon the climate forcing. 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 STD = 0.5 and STD = 0.2, respectively. The smallest STD equals 0.1 (ORCHIDEE-PGFv2) and the highest 0.9 (JULES-GSWP3) (table 3). The exceedingly large STD in JULES-GSWP3 results from a single strong RRB fluctuation around 2002. Amid the models, similar values of 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).
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 interannual variability of global biomass ( figure 1(d)). The observed RRB to MEI regression shows a R 2 value of 0.5 ( figure 1(b)), as opposed to a model ensemble mean R 2 of 0.2. The highest value R 2 = 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 −1 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 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 2 = 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 The simulation CARAIB-PGFv2 shows the highest R 2 (0.5).

The regional RRB and underlying differences in NPP and RLB in the global vegetation models.
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.
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 RRB = 0.8 ( figure 3(a)). This increment largely surpasses in magnitude the loss of biomass observed in NWA (RRB = −0.08), CES (RRB = −0.07), AMA (RRB = −0.07) and CON (RRB = −0.01). With the exception of CON, the decremental trends of biomass associated to the observed regional values of RRB are significant. In contrast, the 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 (RRB = 0.8), followed by WES (RRB = 0.6), NWA (RRB = 0.4), CON (RRB = 0.2) and AMA (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 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).
The individual simulations show positive values of 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 RRB values in the simulations are significant across the regions. The simulated increments of biomass show their best mutual agreement in AMA, with RRB ranging between 0.06 (VEGAS-WFDEI) and 0.3 (JULES-GSWP3). On the contrary, the largest inter-simulation discrepancies in 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 RRB occurs in CES, with values ranging between 0.1 (VEGAS-GSWP3) and 1.4 (LPJmL-GSWP2). In WES, 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 RRB results 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 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 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 andNWA, between years 2001 and2004. 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).

The RRB concept, its interpretation and caveats
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).

Global aspects of RRB
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 plantphysiological effects of increasing atmospheric CO 2 , 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 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.

Regional aspects of RRB
The regional value of 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 RRB is negative in the observation, and positive in all of the simulations. A similar situation is found for AMA. The analysis of 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 RRB in the simulations may be related to the CO 2 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 2 (Smith et al 2016) over the period 1982-2011, has pointed to an overestimation of incremental NPP trends. The excessive CO 2 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 2 fertilization effect in all of the GVM simulations cannot be excluded, its role in explaining differences in 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 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 2 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 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 2 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 2 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 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 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 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 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 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 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 RRB in AMA, compared to the rest of the regions. On the contrary, the largest spread of simulated RRB occurs in CON. The positive and negative extreme values of 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).

Conclusions
The historical changes of biomass in natural terrestrial ecosystems described by the simulations of an ensemble of eight GVMs was evaluated against time-series of AGB, reconstructed from the signals of Vegetation Optical Depth retrieved from satellite passive microwave sensors over the period 1993-2010. For this purpose, the model output and the AGB observation were compared in terms of the relative rate of change of biomass RRB. In particular, the study focused on: the temporal mean and standard deviation of global RRB and its relation to ENSO events, as well as on the regional aspects of the temporal mean and annual variations of RRB, in relation to changes in NPP and RLB. The main findings of the analysis are: (a) the temporal mean of RRB in the models is one order of magnitude higher than in the observed AGB records; (b) the observed RRB shows a significant sensitivity of AGB in the Northern Hemisphere to ENSO events, but not in the tropics; (c) on the contrary, some of the GVMs show a significant sensitivity of RRB to ENSO in the tropics, while none of them reproduces the observed sensitivity in the Northern Hemisphere; (d) in most of the models the annual variations of RRB are mostly influenced by NPP, in comparison to RLB. Overall, these findings underline the importance of conducting a further detailed analysis of both RLB and its equilibrium counterpart TRB and more specifically on their background values and sources of variation, including the effect of CO 2 , climate extremes and disturbances. In turn, this may require improving our current capacity to accurately quantify annual changes in biomass and NPP. In addition, a detailed and updated land-cover classification is necessary to effectively disentangle the physiological and anthropogenic components of TRB and RLB.