Global vegetation variability and its response to elevated CO2, global warming, and climate variability - a study using the offline SSiB4/TRIFFID model and satellite data

. The climate regime shift during the 1980s had a substantial impact on the terrestrial ecosystems and vegetation at different scales. However, the mechanisms driving vegetation changes, before and after the shift, remain unclear. In this study, we used a biophysical dynamic vegetation model to estimate large-scale trends in terms of carbon ﬁxation, vegetation growth, and expansion during the period 1958–2007, and to attribute these changes to environmental drivers including elevated atmospheric CO 2 concentration (hereafter eCO 2 ), global warming, and climate variability (hereafter CV). Simulated leaf area index (LAI) and gross primary production (GPP) were evaluated against observation-based data. Signiﬁcant spatial correlations are found (correla-tions > 0.87), along with regionally varying temporal correlations of 0.34–0.80 for LAI and 0.45–0.83 for GPP. More than 40 % of the global land area shows signiﬁcant positive (increase) or negative (decrease) trends in LAI and GPP during 1958–2007. Regions over the globe show different characteristics in terms of ecosystem trends before and after the 1980s. While 11.7 % and 19.3 % of land have had consistently positive LAI and GPP trends, respectively, since 1958, 17.1 % and 20.1 % of land saw LAI and GPP trends, respectively, reverse during the 1980s. Vegetation fraction cover (FRAC) trends, representing vegetation expansion and/or shrinking, are found at the edges of semi-arid areas and polar areas. Environmental drivers affect the change in ecosystem trend over different regions. Overall, eCO 2 consistently contributes to positive LAI and GPP trends in the tropics. Global warming mostly affects LAI, with positive effects in high latitudes and negative effects in subtropical semi-arid areas. CV is found to dominate the variability of FRAC, LAI, and GPP in the semi-humid and semi-arid areas. The eCO 2 and global warming effects increased after the 1980s, while the CV effect reversed during the 1980s. In addition, plant competition is shown to have played an important role in determining which driver dominated the regional trends. This paper presents new insight into ecosystem variability and changes in the varying climate since the 1950s.


10
Y. Liu et al.: Global vegetation variability and its response to elevated CO 2 , global warming, and climate variability neni et al., 1997;Piao et al., 2011;Ichii et al., 2013;Los, 2013) and leaf area index (LAI, defined as the one-side leaf area per ground area) products Zhu et al., 2016). The Earth's terrestrial vegetation has acted as an important carbon sink in the past 60 years (Ballantyne et al., 2012;Le Quéré et al., 2013), with a significantly strengthening carbon sink rate, about 0.06 Pg C yr −2 , after the 1980s , revealing growth in plant productivity (Nemani et al., 2003;Anav et al., 2015). In the meantime, vegetation fractional coverage (hereafter FRAC) has been changing, including some large-scale increases in total vegetation cover (Piao et al., 2005;Donohue et al., 2009;Mc-Dowell et al., 2015), and shifts in the spatial distributions of plants species, such as woody plants encroachment in the savanna area (Stevens et al., 2017) and shrubification in the tundra biome (Epstein et al., 2012;Mod and Luoto, 2016).
Many studies have attributed these large-scale ecosystem trends to climatic drivers and eCO 2 after applying statistical methods to satellite-based observations or the results from process-based land surface models (Myneni et al., 1997;Liu et al., 2006;Ichii et al., 2013;Mao et al., 2013;Piao et al., 2015;Schimel et al., 2015;Devaraju et al., 2016;Zhu et al., 2016;Smith et al., 2016). Statistical regression and cross-correlation have been applied to attribute the recent biosphere changes to precipitation, temperature, and solar radiation variability Myers-Smith et al., 2015). Results from these studies indicated that northern mid-to high-latitude NDVI anomalies were positively correlated with temperature and positively associated with precipitation in temperate to tropical semi-arid and arid regions . However, statistical methods rarely isolate the drivers' contribution to the interannual or decadal variability of the terrestrial ecosystem (Ahlbeck, 2002;Piao et al., 2015). Moreover, satellite products only cover the period after 1980 .
Process-based land surface models overcome these limitations and are also able to include atmospheric CO 2 as an external driver. Dynamic global vegetation models (DGVMs) simulated vegetation cover changes in response to changes in climate and atmospheric CO 2 , and update associated surface characteristics such as plant functional type (PFT) distribution and LAI (Claussen and Gayler, 1997;Bonan et al., 2002;Sitch et al., 2003;Woodward and Lomas, 2004;Krinner et al., 2005;Zeng et al., 2005;Zaehle and Friend, 2010;Lawrence et al., 2011;Zhang et al., 2015). By applying DGVMs in a model intercomparison project (called TRENDY), a general consensus has been reached that eCO 2 explains the greater part of the increasing trend of LAI and gross primary production (GPP) towards the end of the 1980s (Schimel et al., 2015;Zhu et al., 2016). Air temperature, precipitation, land use and land cover change, and nitrogen decomposition also play roles in the changing terrestrial biosphere (Cramer et al., 2001;Schimel et al., 2015;Zhu et al., 2016). However, DGVMs should be applied with caution. The Coupled Model Intercomparison Project phase 5 (CMIP5) reported that most DGVMs overestimated LAI in comparison to Global Inventory Modelling and Mapping Studies (GIMMS) data (Murray-Tortarolo et al., 2013;. In addition, large discrepancies between models were found when predicting ecosystem variability and trends Zhu et al., 2017). Unsurprisingly, the dominant factors obtained from different models are often significantly different (Beer et al., 2010;Huntzinger et al., 2017). Furthermore, DGVM simulations were sensitive to meteorological forcing data (Slevin et al., 2017;Wu et al., 2017). Therefore, a comprehensive evaluation of large-scale terrestrial ecosystem vegetation trends and potential drivers is crucial for improved DGVM application.
Most ecosystem trend detection and attribution studies have focused on the period after the 1980s when satellite data were available (Myneni et al., 1997;Schimel et al., 2015;Zhu et al., 2016). However, a climate regime shift, identified by abrupt shifts in temperature, precipitation, and other climate variables (e.g., wind speed and sea surface pressure), was observed during the 1980s (Gong and Ho, 2002;Lo and Hsu, 2010;Reid et al., 2016). The responses of vegetation to these climate shifts have not yet been comprehensively investigated, especially at the level of individual PFTs.
In this study, we investigate the effect of eCO 2 and climate drivers including global warming and climate variability (i.e., meteorological forcing excluding global warming, referred to as CV) on the trends of FRAC, LAI, and GPP during the period 1958-2007 by using the SSiB4/TRIFFID (Simplified Simple Biosphere model version 4/Top-down Representation of Interactive Foliage and Flora Including Dynamics) DGVM (Xue et al., 1991;Cox, 2001;Zhan et al., 2003;Zhang et al., 2015;Harper et al., 2016) at both grid and PFT levels, and using satellite products whenever they are available. Changes in the ecosystem trends are attributed to changes in eCO 2 and climate effects, focusing particularly on the climate regime shift during the 1980s. The key focuses of this paper are on (1) how the vegetation trends change before and after the 1980s, and (2) the effect of climate regime shifts during the 1980s on the vegetation trend before/after the 1980s.

Model description
SSiB is a biophysically based model which simulates fluxes of radiation, momentum, sensible heat, and latent heat, as well as runoff, soil moisture, and surface temperature (Xue et al., 1991). A photosynthesis model (Collatz et al., 1991(Collatz et al., , 1992 has been implemented into SSiB to calculate carbon assimilation, forming SSiB2 (Zhan et al., 2003). The TRIF-FID DGVM (Cox, 2001) was coupled to SSiB version 4 (Xue et al., 2006) to calculate vegetation dynamics, including relevant land surface characteristics of vegetation cover and structure. Zhang et al. (2015) updated the competition dominance hierarchy from tree-shrub-grass (i.e., trees dominate shrubs and grasses, and shrubs dominate grasses) to tree-grass-shrub but still allowed shrubs and grasses to compete for sunshine and space. SSiB4 estimates net plant photosynthesis assimilation rate, autotrophic respiration, and other surface conditions such as canopy temperature and soil moisture for TRIFFID. TRIFFID updates the coverage of a PFT based on the net carbon available to it and the competition with other PFTs, which is controlled by the Lotka-Volterra equations. Vegetation is described by leaf, wood, and root with associating carbon pools. Leaf phenology is simulated as a function of canopy temperature and soil moisture. In addition, tundra was separated from the original single shrub category in order to better reflect the arctic biomes. Evergreen and deciduous broadleaf trees are also separated as different PFTs. To date, SSiB4/TRIFFID therefore includes seven PFTs: (1) evergreen broadleaf trees, (2) deciduous broadleaf trees, (3) needleleaf broadleaf trees, (4) C 3 grasses, (5) C 4 plants, (6) shrubs, and (7) tundra.

Experimental design
In this study, SSiB4/TRIFFID was used to simulate the global vegetation distribution and assess the sensitivity of ecosystem trends to climate and eCO 2 . Two sets of simulations were performed: (1) a 100-year quasi-equilibrium simulation driven by climatological forcing, and (2) sensitivity simulations driven by real forcing from 1948 to 2007 (Table 1). In the first set, SSiB4/TRIFFID was driven with the climatological forcing and 1948 CO 2 concentration to reach a steady state, which was used as the initial condition in the second set of simulations.
Using the quasi-equilibrium simulation results as the initial condition, the historical meteorological forcing and yearly updated atmospheric CO 2 concentration were used to drive SSiB4/TRIFFID from 1948 to 2007. In this control simulation, we evaluated the model performance in reproducing the climatology and variability of vegetation coverage, LAI, and GPP in comparison with multiple observationbased datasets. The long-term trends were diagnosed before and after the climate regime shift of the 1980s. Furthermore, three sets of experiments were conducted to quantify the effects of environmental drivers (climate and CO 2 ) and vegetation competition on the ecosystem trends. These experiments were designed as following: 1. For the Fixed-CO 2 experiment, the model was driven by the same meteorological forcing as the control experiment, but CO 2 concentration was fixed at the level of 1948 (310.33 ppm). The difference between control experiment and Fixed-CO 2 indicates the eCO 2 effect.
2. For the Detrend-Temp experiment, the mean warming trend over each 10 • of latitude, from 1948 to 2007, was subtracted in this experiment. Then, the detrended temperature along with other meteorological forcing and annually varying CO 2 concentration were used to drive the model. Subtraction of Detrend-Temp from the control experiment isolates the effect of global warming.
3. For the climate variability experiment, subtraction of both Fixed-CO 2 and Detrend-Temp from the control experiment was regarded as representing the effect of CV.

Data
A SSiB vegetation and soil map is used as the preliminary initial condition for the quasi-equilibrium simulation. Three-hourly meteorological forcing data from 1948 to 2007 (Sheffield et al., 2006) are used for this study. The observation-based LAI and GPP products 12 Y. Liu et al.: Global vegetation variability and its response to elevated CO 2 , global warming, and climate variability Xiao et al., 2014;Jung et al., 2009) are used to validate and calibrate the model to produce proper vegetation spatial distribution and temporal variability.

Initial condition for equilibrium simulation
There are different ways to initialize the surface condition for the quasi-equilibrium simulation. Based on our previous study (Zhang et al., 2015), we set up the initial condition using the SSiB vegetation map and SSiB vegetation table, which are based on ground survey and satellite-derived information (Dorman and Sellers, 1989;Xue et al., 2004;Zhang et al., 2015) with 100 % occupation at each grid point for the dominant PFT and zero for other PFTs.

Meteorological forcing data
The Princeton global meteorological dataset version 1 for land surface modeling (Sheffield et al., 2006) is used to drive SSiB4/TRIFFID for the period of 1948-2007. This dataset is constructed by combining a suite of global observation-based datasets with the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis starting from 1948 (http:// hydrology.princeton.edu/data/pgf/, last access: 3 July 2011). The spatial resolution is 1 • × 1 • and temporal interval is 3-hourly. This dataset includes surface air temperature (K), pressure (Pa), specific humidity (g kg −1 ), wind speed (m s −1 ), downward shortwave radiation flux (W m −2 ), downward longwave radiation flux (W m −2 ), and precipitation (mm day −1 ). Its 60-year mean climatology with 3 h intervals from 1 January to 31 December was generated to drive quasi-equilibrium simulation.

Observation-based data
Two satellite-derived global land cover maps are used to evaluate the vegetation distribution in both quasi-equilibrium and real-forcing simulations. The Global Land Cover (GLC) database for the year 2000 (Bartholome et al., 2002) is used and was derived from SPOT (Satellite Pour l'Observation de la Terre) at a spatial resolution of about 1 km. This dataset provides a global map with one consistent legend, as well as regional maps with separate legends containing more detail for certain regions. For instance, tundra is not included in the global legend but is included in the regional product for northern Eurasia (Bartalev et al., 2003). The regional land cover maps (download from http://forobs.jrc.ec. europa.eu/products/glc2000/glc2000.php, last access: 23 August 2017) are used to calculate the land cover fraction by counting the percentage of each PFT in a 1 • grid, then are merged to form a global land cover fraction map. Other than GLC2000, the land cover type climate modeling grid (CMG) product (MCD12C1), which is derived using the same algorithm that produces the V051 global 500 m land cover type product (MCD12Q1) from the observation input of Terra and Aqua on board the Moderate Resolution Imaging Spectroradiometer (MODIS), is also used as a reference (Friedl et al., 2010). The MODIS-MCD12C1 product provides land cover fraction at the spatial resolution of 0.05 • , which is then converted to 1 • resolution. SSiB4/TRIFFID only includes primary land cover types, while both GLC2000 and MODIS International Geosphere-Biosphere Programme (IGBP) have more sub-level classes. For easy comparison of the distribution of dominant vegetation types with different products, we hierarchically combine the GLC2000 and the MODIS IGBP classifications with the SSiB4/TRIFFID PFTs.
SSiB4/TRIFFID GPP was evaluated using the FLUXNET model tree ensemble (MTE) GPP, which was downloaded from https://www.bgc-jena.mpg.de/geodb/projects/Data.php (last access: 7 March 2017). The FLUXNET-MTE GPP was upscaled from FLUXNET observations of carbon dioxide flux to the global scale using the machine learning technique MTEs. This method was trained to predict site-level GPP based on remote sensing indices, climate and meteorological data, and information on land use (Jung et al., 2009). This dataset provides global monthly mean GPP at 0.5 • spatial resolution for the period from 1982 to 2011. The FLUXNET-MTE GPP was resampled to 1 • spatial and monthly temporal resolution.

Vegetation initial conditions
The DGVM's initial conditions for long-term simulations are obtained from a quasi-equilibrium solution in a long-term simulation using the climatological forcing as presented in Sect. 2.3.2. The effect of large-scale disturbance (LSD) on regulating tree fraction over the savanna areas is also investigated. The parameter for LSD is tuned to generate a reasonable tree cover distribution there.

Quasi-equilibrium simulation
DGVMs require initial conditions for a number of state variables. DGVMs normally take 50-1000 years' simulation under specified meteorological forcing to reach this steady state Zeng et al., 2008). Since our purpose was to generate the initial condition for the decadal simulations, we applied a shortcut to reach the quasi-equilibrium coexistence of PFTs under the climatological forcing. We started the model from a SSiB 1 • dominant vegetation map (Xue et al., 2004), with 100% occupation of the dominant PFT and zero for other PFTs at each grid point. The 1948-2007 averaged meteorological forcing along with 1948 CO 2 concentration were used to drive SSiB4/TRIFFID for 100 years. SSiB4/TRIFFID is a water, energy, and carbon balanced model. Plant expansion and biotic properties such as vegetation height and LAI are constrained by carbon allocation. FRAC links the carbon accumulation within plants and intra-species competition via a system of the Lotka-Volterra equations (Cox, 2001). In the first 10 years of the simulation, most PFTs' FRACs change rapidly but then take decades to reach a steady state (Fig. 1). To qualify the steady state, we define quasi-equilibrium status as occurring when the rate of change of vegetation fraction is less than 2 % of the mean vegetation fraction, over the last 10 years of simulation. The result shows 100-year spinup time is sufficient for our model. Overall, in the tropical areas (23.5 • S-23.5 • N), C 4 plants and evergreen broadleaf trees are of mixed dominance and coexistent with C 3 grasses, shrubs, and deciduous broadleaf trees. The subtropical areas (23.5-35 • S and 23.5-35 • N in both hemispheres) are dominated by C 3 grasses, C 4 plants, and shrubs with similar occupation for each (∼ 18 %), whereas 40 % of the subtropical areas are occupied by bare land. Needleleaf trees, C 3 grasses, shrubs, and deciduous broadleaf trees are the dominant mix in the temperate zones (35-66.5 • S and 35-66.5 • N, particularly in the Northern Hemisphere). Over the polar areas (66.5-90 • S and 66.5-90 • N in both hemispheres), shrubs and tundra are of mixed dominance (Fig. 1).

Effect of large-scale disturbance
Large-scale disturbance, such as fire, and insect outbreaks alter physical structure and/or arrangement of biotic elements with great effect. TRIFFID introduces only a global uniform and PFT-dependent parameters to represent the rate of vegetation loss caused by LSD (units: yr −1 ). The preliminary quasi-equilibrium run shows that, under this approach, trees extended into the South American and African savanna areas (Fig. 2a), where the climate acting alone would seem to favor tree growth (Bond et al., 2005). However, major eco-14 Y. Liu et al.: Global vegetation variability and its response to elevated CO 2 , global warming, and climate variability logical disturbances vary spatially and temporally (Giglio et al., 2006). We raised the LSD coefficient (largely representing fire at this scale) from 0.004 to 0.04 (yr −1 ) for tree PFTs that coexist with C 3 grasses and C 4 plants. With the updated setting, SSiB4/TRIFFID produced reasonable dominant tree coverage over the tropical rainforest areas (Fig. 2b) in comparison with the GLC2000 dataset (Fig. 2c). The global PFT distributions in the equilibrium run are close to the results using the real meteorological forcing.

Model evaluation of simulated vegetation distribution, LAI, and GPP
Using initial conditions derived from the equilibrium run for both biotic and abiotic variables such as FRAC, LAI, vegetation height, soil moisture, and temperature, the model is then driven with the historical meteorological forcing and yearly updated atmospheric CO 2 concentration from 1948 to 2007. In this section, the global (hereafter referring to the regions of 180 • W to 180 • E, 60 • S to 75 • N) distributions of the simulated FRAC, LAI, and GPP are compared to the observation-based datasets to ensure SSiB4/TRIFFID generates a reasonable spatial pattern and temporal variability of those variables, which provides a base for further assessment of the impact of external drivers on some surface variables.

Vegetation spatial distribution
Satellite-derived products are used to assess the modelsimulated mean FRAC averaged over the last 10 years (Fig. 2). Overall, the model-simulated vegetated land covers 79.6 % of the global land surface, which is less than GLC2000 estimation (80.8 %) and higher than the MODIS estimation (79.3 %) (Fig. S1 in the Supplement). There is no human activity included in the model simulation, as such an agricultural category is not included. Therefore, in the following vegetation coverage comparison with the GLC products, simulated PFT coverage in the grid boxes with agriculture is reduced based on the GLC agriculture fraction. On this basis, the total simulated tree cover (the sum of evergreen broadleaf trees, deciduous broadleaf trees, and needleleaf trees) is 28.8 %, close to 29.8 % in GLC2000. The evergreen broadleaf trees in the Amazon, central Africa, and southeast Asia, deciduous broadleaf trees in southeast North America, and needleleaf trees in the high to midlatitudes of North America and Eurasia are reasonably predicted. SSiB4/TRIFFID simulates 12.7 % C 3 grass occupation, which is slightly higher than 11.9 % in GLC2000, with reasonable simulation in the midlatitudes in both hemispheres such as the central US, Eurasian steppes, South America, south and east Africa, and east Australia. The model simulates 10.1 % natural C 4 plants, compared to 7.9 % in GLC2000. The discrepancy could be partially attributed to the absence of C 4 plants in some GLC2000 regional maps (such as southeast Asia). The global GLC2000 map is assembled from these regional maps. In fact, a satellite-based physiological model simulation estimated 13.9 % of C 4 plant coverage with no agriculture category (Still et al., 2003).
In the SSiB4/TRIFFID prediction without excluding agricultural land, the C 4 plants cover 13.5 %, close to Still et al.'s (2003) estimation. C 4 plants are primarily located in South American and African savanna areas, the Indian subcontinent, southeast Asia, the southeast US, and northern Australia. The model simulates 15.9 % shrubs and tundra occupation, which is close to 16.7 % in GLC2000, with shrubs primarily located in the semi-arid areas in both hemispheres and the pan-Arctic area, while tundra is located in the pan-Arctic area and Tibetan Plateau (Fig. 3, and also see vegetation fractional distribution in Fig. S2).

Leaf area index
This section discusses the spatial and temporal correlations between the simulations and observations and compares with other model results. Since other published studies on this subject have not excluded agricultural areas when evaluating LAI and GPP simulation, to make our results comparable with others, the agricultural areas are not subtracted. In fact, the difference between the results with and without the exclusion of agriculture area for our results is less than 0.01. SSiB4/TRIFFID produces a similar global LAI pattern compared to both the GIMMS and GLASS products, confirmed by global spatial correlation coefficients of 0.86 (GIMMS, p < 0.05) and 0.87 (GLASS, p < 0.05), and above 0.74 (P < 0.05) against both observations over the Northern Hemisphere (Fig. 4). Previous studies reported spatial correlation coefficients between models and GIMMS LAI over the globe and Northern Hemisphere in the range

Gross primary production
The spatial correlation coefficient between model and FLUXNET-MTE GPP is 0.93 (P < 0.05) (Fig. 5). Anav et al. (2015) reported less than 0.8 correlation against FLUXNET-MTE GPP for multi-model comparison. Over the globe, SSiB4/TRIFFID simulates 151 Pg C yr −1 , greater than the FLUXNET-MTE average of 122 Pg C yr −1 . However, our simulation was still within the range of 130-169 Pg C yr −1 reported by Anav et al. (2015) and 111-151 Pg C yr −1 reported by Piao et al. (2013). In addition to our model's deficiencies (such as lack of N limitation), the lack of CO 2 fertilization during the MTE model training may have contributed to an underestimation in the FLUXNET-MTE GPP. Model performance in predicting temporal variability has been less evaluated in previous studies on ecosystem trend detection and attribution (Ichii et al., 2013;Piao et al., 2013;Zhang et al., 2015). However, performance in estimating LAI and GPP trends and variability has been found to vary among models (Murray-Tortarolo et al., 2013;Piao et al., 2013;Anav et al., 2015;Zhu et al., 2017). To better assess model performance in this regard, we select 13 subregions associated with different regional climate and land cover conditions (Table 2). Although the MTE excludes CO 2 fertilization during its model training, FLUXNET-MTE GPP still incorporates variability at different scales associated with climate variability and is widely used by the community for model evaluation. Both annual LAI and GPP correlation coefficients are calculated over the period of 1982-2007. Globally, the correlations for annual mean LAI between the SSiB4/TRIFFID and the satellite-based products are 0.58 (P < 0.05) for GIMMS and 0.64 (P < 0.05) for GLASS. The correlation for annual GPP is 0.59 (P < 0.05) between the SSiB4/TRIFFID and FLUXNET-MTE GPP. Regionally, LAI correlations over west Africa are 0.79 (P < 0.05) with GIMMS and 0.77 (P < 0.05) with GLASS (Fig. 6), and GPP correlation is 0.80 (P < 0.05) with FLUXNET-MTE GPP in that region. For other semiarid areas in western North America, South American savanna areas, and east Africa, the simulated LAI significantly matches at least one of the two reference datasets with correlation in the range of 0.46-0.58 (P < 0.05). The simulated GPP correlations with FLUXNET-MTE GPP are in the range of 0.63-0.70 (P < 0.05). The LAI over the forested areas is better correlated to GLASS LAI, while GPP is only significantly corrected to FLUXNET-MTE GPP over the boreal forests. Over the cold regions (subarctic and Tibetan Plateau), SSiB4/TRIFFID matches the reference data in a varying range of 0.46-0.74 (P < 0.05) for LAI and 0.49-0.72 (P < 0.05) for GPP.
Distinct decadal variabilities are identified in most subregions. For instance, trends reverse sign in west Africa (from negative to positive) and western North America (from positive to negative) during the 1980s (Fig. 6). Areas with enhancement in trend slopes, such as the subarctic after the climate regime shift in the 1980s, will be discussed in detail in the next section.
Generally, SSiB4/TRIFFID simulates reasonable predictions of terrestrial ecosystem climatology and variability compared to the observation-based datasets (Table 3). Compared to other DGVMs, SSiB4/TRIFFID shows above average performance in reproducing the spatial distribution but with certain bias in absolute numbers. In particular, SSiB4/TRIFFID captures the ecosystem temporal variabilities over different regions across the world, which provides a basis for pursuing the ecosystem trend detection and attribution study presented in the next section.

Three major types of vegetation trend change since the 1950s
The climate regime shifted abruptly during the 1980s, giving rise to changes in ecosystem trends in many parts of the world. Here, we compare trends of FRAC, LAI, and GPP over two periods : 1958-1982 and 1982-2007. The model performance for the second period, for which satellite observations are available, is evaluated in Sect. 3.3.1. Spatial patterns of the trends are shown in Fig. 7. At the global scale, significant vegetation trends are only found in the simulations after the 1980s. During this period, FRAC increases at the rate of 0.032 yr −1 . LAI has a positive trend of 0.0029 yr −1 , which matches very well to GIMMS' results (0.0029 yr −1 ). GPP has a positive trend of 2.22 g C m −2 yr −2 , within the range of 1.60-4.69 g C m −2 yr −2 for GPP over similar periods Yue et al., 2015). In contrast to LAI and GPP, there are relatively few areas with a significant simulated FRAC trend (Fig. 7).
For the global land surface, over 40.2 % has a significant LAI trend from 1958 to 2007, and over 48.1 % has a significant GPP trend. In response to the climate regime shift     during the 1980s, the terrestrial ecosystem has three major trend changes in different parts of the world after the 1980s (Table 4). (1) There is trend sign reversal from negative to positive in the East Asian monsoon area, west Africa, central Asia, and the eastern US, over 14.2 % (LAI) and 11.4 % (GPP) of the land surface. In particular, west Africa experiences the largest vegetation deterioration in the world before the 1980s, associated with LAI and GPP reductions of 0.0258 yr −1 and 18.54 g C m −2 yr −2 , respectively -approximately 10 times the trends of the global average. After the 1980s, recovery is simulated at the rate of 0.0137 yr −1 and 8.02 g C m −2 yr −2 for LAI and GPP, respectively. (2) Trend Table 4. Linear trend of vegetation FRAC, LAI and GPP changes during P1 (1958) and P2 (1982-2007  sign reversal from positive to negative is found in western North America, South American savanna, and east Africa, which accounted for 2.9 % (LAI) and 2.7 % (GPP) of the land surface.
(3) There are consistent positive trends but substantially enhanced by at least 50 % of prior period trends after the 1980s, which are found in equatorial rainforest areas, boreal forest areas, south Africa, north Australia, subarctic areas, and the Tibetan Plateau, representing over 11.7 % (LAI) and 19.3 % (GPP) of the land surface. There are also areas with consistent positive trends but no substantial change during the entire period or other types of trend change over much smaller areas. The first three major trend changes as described above will be discussed in the following sections.

Attribution of three environmental drivers on ecosystem trends
3.4.1 Global overview of three simulated environmental drivers' effects on the ecosystem trends Sensitivity experiments were conducted to isolate the contributions of elevated atmospheric CO 2 concentration, global warming, and climate variability. The differences between the control experiment and Fixed-CO 2 show that eCO 2 stimulated vegetation growth mainly in the equatorial areas and eastern North America, western Europe, and eastern China in the midlatitudes. Substantially enhanced positive trends are found after the 1980s for both LAI and GPP over those areas (Fig. 8). eCO 2 promoted rainforest LAI increase only after the 1980s; however, its effect on GPP appeared during the entire period. GPP is directly linked to CO 2 through the photosynthesis process, while LAI, in addition to the photosynthesis process, is also affected by respiration and carbon allocation in plants, which are influenced by both climate and eCO 2 (O'Sullivan et al., 2017). The differences between the control experiment and Detrend-Temp show that global warming has minor effects on the trends of LAI and GPP before the 1980s (Fig. 9). After the 1980s, the rapidly enhanced warming contributes positive LAI trends at high latitudes, while the GPP change seems less substantial. Meanwhile, there are negative trends due to heat stress in low latitudes, particularly in the semiarid regions such as South American savanna, east Africa, and central Asia.
The differences between the control experiment and Fixed-CO 2 and Detrend-Temp show the CV effect, which has (a) LAI 1958(b) LAI 1982-2007(d) GPP 1982-2007 (c) GPP 1958GPP -1982   complex influences on the ecosystem. The CV in this study includes contribution of changes in surface pressure, precipitation, surface wind speed, downward longwave and shortwave radiation, and surface humidity, along with temperature that excludes the global warming trend. Precipitation is however found to play a dominant role. The correlation coefficients between the annual mean CV effect on LAI and GPP and annual mean precipitation at the grid points with significant CV effect are greater than 0.60 (P < 0.05). Overall, the CV effect alone can explain the total FRAC trends in the control experiment (Fig. 10). Before the 1980s, CV causes LAI to decrease in East Asian monsoon areas, eastern North America, west Africa, western Europe, central Asia, Siberia, and eastern Australia. The GPP also decreases in these areas except for eastern North America, western Europe, and Siberia. In contrast, the CV effect before the 1980s leads LAI and GPP increase in the Tibetan Plateau and south Asia, western North America, South American savanna areas, east and south Africa, and northern Australia. Due to the climate regime shift, CV has produced the opposite sign to the trends of LAI and GPP in East Asian monsoon areas, central Asia, west Africa, North America, South American savanna areas, and east Africa. In some areas, such as south Africa and northern Australia, persistent precipitation increase/decrease leads to sustained positive/negative trends from the 1950s.
Overall, after the 1980s, the effects of eCO 2 and global warming have been generally enhanced, but the CV effect has exhibited distinctly different regional features before/after the 1980s over many regions in the world. The enhanced or opposite contribution of the primary driver and the changes in their relative importance on the ecosystem trends occur during the 1980s, resulting in different ecosystem responses in many regions across the world.

Dominant factor in influencing trend reversal from negative to positive in west Africa and east Asia
CV is found to be the dominant driver of the ecosystem trends in west Africa, explaining most of the LAI and GPP trends and trend changes (Fig. 11a). Before the 1980s, CV causes C 4 plants' LAI, GPP, and FRAC over the region to decrease, followed by shrubs, whereas eCO 2 caused C 3 grasses' LAI and GPP to slightly increase. Global warming shows little effect during the entire period from 1958 to 2007. The ecosystem trends in west Africa reverse from decreasing to increasing when the precipitation trend changes to increasing after the 1980s, with the major increase in C 3 grasses and shrubs over the region. A previous study using satellite data also showed recent west Africa greening is highly correlated with the precipitation increase (Herrmann et al., 2005). eCO 2 has played a role in increasing C 3 grasses' coverage since the 1950s. However, the PFT competition outcomes reduce the C 4 plant coverage over the region, mainly after the 1980s when eCO 2 has a large impact. As such, the change in regional FRAC overall within west Africa is not significant and has been compromised by positive and negative contributions of the individual PFTs after the 1980s. Regional average trends reverse in the East Asian monsoon area because CV and eCO 2 dominate LAI and GPP trends, before and after the 1980s, respectively. Their combined effects cause trend reversal in the East Asian monsoon areas. CV contributes to decreasing trends of LAI and GPP before the 1980s but with minor effects after the 1980s, while eCO 2 has dominated the PFT LAI and GPP trends since the 1950s, which caused a significant increase in C 3 grasses and trees but significant decrease in C 4 plants (Fig. 11b). Meanwhile, enhanced global warming after the 1980s stimulates C 4 plant growth, but this effect is compromised by its detrimental effect on C 3 grasses after the 1980s. Overall, the CV and eCO 2 relative contribution change during the 1980s dominates the negative to positive trend shift in this area.

Dominant factor in influencing trend reversal from positive to negative in western North America
The eCO 2 effect has persistently caused LAI, GPP, and FRAC increase since the 1950s, while global warming reduced both LAI and GPP only after the 1980s. However, CV dominated the LAI and GPP trends and trend reversal in western North America by causing the dominant PFTs (C 3 grasses and shrubs) to increase/decrease before/after the 1980s (Fig. 11c and Table 5). The CV effect on FRAC change is more complex due to its different effects on LAI and GPP and FRAC expansion in C 3 and shrub PFTs after the 1980s: both C 3 and shrubs expand with LAI and GPP decrease. This discrepancy suggests that expansion might be coupled with carbon fixation less than with growth in the model. We conjecture that CV may promote vegetation expansion into some areas that are largely unvegetated, but this requires further investigation.
3.4.4 Dominant factor in influencing the enhanced positive trend in rainforest, boreal forest, subarctic, and Tibetan Plateau eCO 2 and CV have had persistent positive impacts on tropical rainforest growth in terms of LAI and GPP since the 1950s (Fig. 11d). eCO 2 dominates the LAI and GPP trends in both periods except for the LAI positive trend before the 1980s, which is dominated by CV. LAI trend en-

22
Y. Liu et al.: Global vegetation variability and its response to elevated CO 2 , global warming, and climate variability hancement after the 1980s is associated with increased CO 2 fertilization, while GPP trend enhancement is attributed to increase in both eCO 2 and CV effects. The importance of CO 2 and CV impacts on the rainforests is confirmed by previous analyses on the trends of LAI and NDVI (Hilker et al., 2014;Zhu et al., 2016). eCO 2 and global warming have increased LAI and GPP in North American boreal forest areas since the 1950s and caused significant positive trends over the Eurasian boreal forest area after the 1980s (Fig. 11e). However, due to CVinduced negative effects on tree LAI, no significant trend is found in regional average LAI in boreal areas before the 1980s. The LAI and GPP trend enhancement in the boreal forest areas can be attributed to the enhanced eCO 2 and global warming effects, accompanied by reduced CV negative effects after the 1980s.
North American subarctic areas have enhanced LAI and GPP positive trends after the 1980s, which were caused Table 5. Climate drivers and eCO 2 effect on the trends of FRAC, LAI, and GPP during P1 (1958) and P2 (1982-2007 in subregions with regard to their regional average 1,2 .

Regions
Var by the increase in eCO 2 and CV positive effects, while all three environmental drivers have effects on LAI and GPP positive trends in the Eurasian subarctic (Fig. 11f). Meanwhile, remarkable FRAC changes have been found since the 1950s. Our simulation suggests that global warming continually favored shrub invasion into tundra biomes, except in the Eurasian subarctic before the 1980s. After the 1980s, this shrubification is enhanced due to increase in the warming effect. In contrast, eCO 2 has promoted tundra expansion and shrub decline over subarctic areas since the 1950s, which mitigates the shrubification. Meanwhile, CV plays a role to help tree and C 3 grass expansion into subarctic areas and also alters the shrub and tundra competition.
Over the Tibetan Plateau (Fig. 11g), CV has dominated the positive LAI and GPP trends since the 1950s, except in the case of the GPP increase before the 1980s which is dominated by eCO 2 . The positive trend enhancements for LAI and GPP after the 1980s are caused by the impact of both eCO 2 and CV. Furthermore, our simulation also suggests that CV favors C 3 grasses but harms tundra biome expansion. However, eCO 2 has the opposite effects on those PFTs, in contrast to the CVs.

Discussion
The spatial distribution of the dominant PFT is closely related to large-scale climate (MacDonald, 2003), and DGVMs Y. Liu et al.: Global vegetation variability and its response to elevated CO 2 , global warming, and climate variability are designed to reproduce the observed ecosystem-climate relationship. There is diverse performance in reproducing the spatial distribution and temporal variability of the ecosystems (Murray-Tortarolo et al., 2013;Piao et al., 2013;Anav et al., 2015;Zhu et al., 2017), which resulted in large discrepancies between models in identifying attributed dominant drivers of changes (Beer et al., 2010;Huntzinger et al., 2017). It is therefore important to validate the model performance in reproducing the observed ecosystem variability and spatial distribution first before using DGVMs for attribution studies. As a matter of fact, it is challenging for DGVMs to reproduce PFT coexistence, particularly for the smaller PFTs in semi-humid and semi-arid areas, as they are fragile and sensitive to climate and vulnerable to competition (Fu et al., 2006). Because of that difficulty, Zeng et al. (2008) had to introduce a specific submodel to grow temperate shrubs in the spaces unoccupied by trees and grasses. SSiB4/TRIFFID allows smaller fractions of PFTs to coexist with full competition with other PFTs (Cox, 2001). After modifying the competition coefficients in the Lotka-Volterra equation (Zhang et al., 2015) and updating large-scale disturbance parameters, it produces a reasonable global distribution of temperate shrubs and high-latitude tundra (Fig. 3). During the validation process, some parameters in SSiB4/TRIFFID have also been calibrated (Zhang et al., 2015).
With all these efforts, the SSiB4/TRIFFID-produced LAI and GPP show higher temporal correlation with observation compared to the start-of-the-art offline models in the TRENDY intercomparison project Zhu et al., 2016). The improvement may mainly be due to better capturing of the interannual variability by SSiB4/TRIFFID in semi-arid areas, which has been considered as dominating global interannual variability (Poulter et al., 2014;Ahlstrom et al., 2015). Meanwhile, both simulated LAI and GPP are also well correlated with reference data over the Northern Hemisphere boreal forests. Our evaluation is based on satellite-based products, which are the only sources providing global distribution at long-term periods. Although these products showed a general consistency among them, large relative uncertainties were identified over some regions (Jiang et al., 2017), which contribute to large discrepancy of interannual correlations when the simulated LAI was compared to GIMMS and GLASS LAIs (Fig. 6). It should be pointed out that the GPP simulation over the rainforests exhibit inconsistency with the FLUXNET-MTE GPP. We consider that the missing CO 2 fertilization in FLUXNET-MTE GPP could be a predominant limitation to its GPP there. By and large, SSiB4/TRIFID's performance suggests this model is appropriate to be applied for the attribution study.
SSiB4/TRIFFID simulates increased LAI and GPP after the 1980s (Fig. 8), which is confirmed by observation and the TRENDY models' simulation (Anav et al., 2013 for GPP;Piao et al., 2013 andZhu et al., 2016 for LAI). These increases are considered to be responding to elevated atmospheric CO 2 concentrations and warming surface temperature in high latitudes . Some areas with decreased LAI and GPP are due to decrease in precipitation and/or increase in stress due to warming temperature in low latitudes Zhu et al., 2016). Our study, however, further estimates large-scale trends in all three aspects, i.e., carbon fixation (GPP) and vegetation growth (LAI), as well as expansion, instead of focusing only on one aspect, such as the LAI trend in Zhu et al. (2016). Our results also reveal different LAI and GPP responses to the environmental changes (Figs. 8-10), indicating LAI and GPP are involved in different processes as discussed in O'Sullivan et al. (2017). The results suggest that GPP is more directly linked to atmospheric CO 2 (Fig. 8).
In SSiB4/TRIFFID, net CO 2 assimilation is proportional to the gradient of atmosphere and leaf CO 2 concentration (Zhan et al., 2003; also see the Supplement). Hence, the elevated atmospheric CO 2 concentration leads to an increase in GPP, while LAI in TRIFFID is related to carbon allocation and competition between PFTs. As such, LAI is not only affected by the atmospheric carbon concentration but also other processes, such as phenological processes and the percentage collocated carbon for growth. Therefore, GPP is more sensitive to the change in atmospheric carbon concentration compared to LAI. Integrated analysis and observation with multiple variables, such as LAI and GPP, are required to improve the understanding of vegetation biochemical process and climate effect on ecosystem changes.
The competition between PFTs within a grid box contributes to the ecosystem trend discussed above. Our analysis with grid points has shown intensive interactions between PFTs. For instance, shrubs are found to expand into tundra ecosystems (Fig. 11f), which is linked to climate change (Myers-Smith et al., 2015), particularly to global warming (Tape et al., 2006;Elmendorf et al., 2012) and precipitation (Martin et al., 2017). The response to eCO 2 of a particular PFT not only depends on its own physiological and morphological characteristics but is also determined by the interactions that arise with other PFTs, competing for the same resources (Fig. 11a). Field experiments reported that the differential growth and competitiveness responses of C 3 grasses and C 4 plants to eCO 2 are complex and under debate (Leakey et al., 2009;Lee, 2011;Miri et al., 2012). In this paper, we have discussed the competition between C 3 and C 4 , and its contribution to the trend change. It seems under the elevated atmospheric CO 2 concentration scenario C 3 grasses show enhanced competitive ability over C 4 plants at regional scale (Fig. 11b). Different responses of the coexisting PFTs to the climate regime shift either enhance or mitigate the environmental drivers' contribution at grid-averaged scale.
Furthermore, our results show that the boundary between the Sahara and Sahel have experienced significant variation since the 1950s. Based on the observed precipitation data and the precipitation/NDVI correlation, Thomas and Nigam (2018) suggested a Sahara desert expansion starting from the 1950s. Our results are in agreement at large with the Thomas and Nigam (2018) study but also with substantial differences in the rate at two different climate regimes. A comprehensive discussion on this issue is out of the scope of this paper and will be addressed in a separate paper.

Conclusion
This work employs a biophysical dynamic vegetation model (SSiB4/TRIFFID) to explore the responses of the terrestrial ecosystem to climate variability, global warming, and elevated atmospheric CO 2 concentration during 1948-2007. SSiB4/TRIFFID is evaluated by available satellite data in simulating the land surface carbon fixation, and plant growth and competition. We have shown that the SSiB4/TRIFFID model can simulate the vegetation and temporal variability for the period 1982-2007. A series of sensitivity experiments is then conducted to detect the ecosystem trends and attribute the trends to elevated atmospheric CO 2 concentration (eCO 2 ), global warming, and CV.
The effects of the external drivers on the ecosystem trends manifest distinct spatial and temporal characteristics. For the global land surface, over 40.2 % has a significant LAI trend and over 48.1 % had a significant GPP trend during the 1958-2007 period. In responding to the climate regime shift during the 1980s, the terrestrial ecosystem had three major changes in different parts of the world after the 1980s. Over 14.2 % (LAI) and 11.4 % (GPP) of the land surface, primarily located in the East Asian monsoon area, west Africa, central Asia, and the eastern US, had trend sign reversal from negative to positive. In contrast, trend reversal from positive to negative is found in western North America, South American savanna, and east Africa, which accounted for 2.9 % (LAI) and 2.7 % (GPP) of the land surface. Meanwhile, there are consistent positive trends substantially enhanced in equatorial rainforest areas, boreal forest areas, south Africa, north Australia, subarctic areas, and the Tibetan Plateau, representing over 11.7 % (LAI) and 19.3 % (GPP) of the land surface, respectively.
In general, the major types of trend change are attributed to the changes in relative contributions of environmental drivers and, consequently, the changes in the dominant driver, or changes in the dominant driver's "direction" in its effect (enhancing or suppressing) on the ecosystem. The eCO 2 stimulates vegetation growth through fertilization effects mainly in the equatorial areas, as well as eastern North America, western Europe, and eastern China in the midlatitudes. The rapidly enhanced global warming after the 1980s contributes positive LAI trends at high latitudes, while the GPP change seems less substantial. Meanwhile, there are negative trends in LAI and GPP due to the heat stress in low latitudes. CV dominates the variability of FRAC, LAI, and GPP in the semi-humid and semi-arid areas. The overall effects on the ecosystem are the integrated contribution of all environmental drivers. Data availability. SSiB4/TRIFFID simulated vegetation fraction, LAI, and GPP are available at https://ucla.box.com/v/ssib4-offline (Liu and Xue, 2018).