Identifying environmental controls on vegetation greenness phenology through model-data integration

16 Existing dynamic global vegetation models (DGVMs) have a limited ability in reproducing 17 phenology and decadal dynamics of vegetation greenness as observed by satellites. These 18 limitations in reproducing observations reflect a poor understanding and description of the 19 environmental controls on phenology, which strongly influence the ability to simulate longer 20 term vegetation dynamics, e.g. carbon allocation. Combining DGVMs with observational data 21 sets can potentially help to revise current modelling approaches and thus to enhance the 22 understanding of processes that control seasonal to long-term vegetation greenness dynamics. 23 Here we implemented a new phenology model within the LPJmL (Lund Potsdam Jena 24 managed lands) DGVM and integrated several observational data sets to improve the ability 25 of the model in reproducing satellite-derived time series of vegetation greenness. Specifically, 26 we optimized LPJmL parameters against observational time series of the fraction of absorbed 27 photosynthetic active radiation (FAPAR), albedo and gross primary production to identify the 28 Formatvorlagendefinition: Überschrift 3: Einzug: Links: 0,32 cm, Tabstopps: 1,59 cm, Listentabstopp + Nicht an 1,27 cm


Introduction
The greenness of the terrestrial vegetation is directly linked to plant productivity, surface roughness and albedo and thus affects the climate system (Richardson et al., 2013).Vegetation greenness can be quantified from satellite observations for example as Normalized Difference Vegetation Index (NDVI) (Tucker, 1979).NDVI is a remotely sensed proxy for structural plant properties like leaf area index (LAI) (Turner et al., 1999), green leaf biomass (Gamon et al., 1995) and plant productivity.In particular, NDVI of green vegetation has a linear relationship with the fraction of absorbed photosynthetic active radiation (FAPAR) (Fensholt et al., 2004;Gamon et al., 1995;Myneni et al., 1995Myneni et al., , 1997b;;Myneni and Williams, 1994).Satellite-derived FAPAR estimates are Published by Copernicus Publications on behalf of the European Geosciences Union.
often used to estimate terrestrial photosynthesis (Beer et al., 2010;Jung et al., 2008Jung et al., , 2011;;Potter et al., 1999).Decadal satellite observations of NDVI demonstrate widespread positive trends ("greening") especially in the high-latitude regions (Lucht et al., 2002;Myneni et al., 1997a;Xu et al., 2013) but also in the Sahel, southern Africa and southern Australia (Fensholt and Proud, 2012;de Jong et al., 2011de Jong et al., , 2013b)).Surprisingly, these trends are accompanied by negative trends ("browning") which were observed regionally in parts of the boreal forests of North America and Eurasia, and in parts of eastern Africa and South America (Baird and Verbyla, 2012;Bi et al., 2013;de Jong et al., 2013b).Regionally different causes have been identified for the observed greening and browning trends.The greening of the high latitudes is supposed to be mainly induced by rising air temperatures (Lucht et al., 2002;Myneni et al., 1997a;Xu et al., 2013).Browning trends in subtropical regions were related to changing drought conditions and land use change (Cook and Pau, 2013;van Leeuwen et al., 2013).On the other hand, the environmental controls on the browning of boreal forests have been intensively investigated but no concluding or general explanation has been found so far (Barichivich et al., 2014;Beck et al., 2011;Beck and Goetz, 2011;Bunn et al., 2007;Goetz et al., 2005;Piao et al., 2011;Wang et al., 2011).Trends in vegetation greenness are often related to changes in vegetation phenology like an earlier onset and an associated lengthening of the growing season in mid-and high-latitude regions (Atzberger et al., 2013;Høgda et al., 2001Høgda et al., , 2013;;Tucker et al., 2001;Zeng et al., 2011).Changes in vegetation greenness are linked to changes in primary production and thus affect atmospheric CO 2 concentrations and the terrestrial carbon cycle (Barichivich et al., 2013;Keeling et al., 1996;Myneni et al., 1997a).Additionally, vegetation greenness affects the climate system by influencing surface albedo.For example, greening trends in high latitudes are associated with decreasing surface albedo (Urban et al., 2013) which alters the surface radiation budget (Loranty et al., 2011).This can potentially further contribute to a warming of Arctic regions (Chapin et al., 2005).Thus, satellite observations of vegetation greenness demonstrate the recent interactions and changes between terrestrial vegetation dynamics and the climate system.
Dynamic global vegetation models (DGVM) or generally climate/carbon cycle models are used to analyse and project the response of the terrestrial vegetation to the past, recent and future climate variability (Prentice et al., 2007).DGVMs can be used to explain observed trends in vegetation greenness (Lucht et al., 2002) or to quantify the related terrestrial CO 2 uptake.While most global models simulate an increasing uptake of CO 2 by the terrestrial vegetation under future climate change scenarios, the magnitude of future changes in land carbon uptake largely differs among models (Friedlingstein et al., 2006;Sitch et al., 2008).The spread of land carbon uptake estimates among DGVMs might be partly related to insufficient representations of vegetation phenology and greenness (Richardson et al., 2012).Coupled climate-carbon cycle models and uncoupled DGVMs have been compared against 30-year satellite-derived time series of LAI (Anav et al., 2013;Murray-Tortarolo et al., 2013;Zhu et al., 2013).Models usually overestimate mean annual LAI in all biomes and have a too long growing season because of a delayed season end (Anav et al., 2013;Murray-Tortarolo et al., 2013;Zhu et al., 2013).Additionally, most DGVMs have more positive LAI trends than the satellite-derived LAI product, i.e. they underestimate browning trends in boreal forests while a few DGVMs do not reproduce the general greening of the high latitudes (Murray-Tortarolo et al., 2013).The limitations of DGVMs in reproducing observed LAI or FAPAR time series is mostly related to limited phenology routines that often miss environmental controls on seasonal leaf development (Kelley et al., 2013;Murray-Tortarolo et al., 2013;Richardson et al., 2012).In conclusion, with improved modelling approaches for vegetation phenology and greenness, DGVMs can potentially more accurately reproduce the recent, and project the future response of the terrestrial vegetation to climate variability.
Past studies successfully demonstrated the use of vegetation greenness observations to improve stand-alone phenology models or to optimize phenology and productivityrelated parameters in DGVMs.The growing season index (GSI) is an empirical phenology model that is used to estimate seasonal leaf developments (Jolly et al., 2005).Empirical parameters of GSI have been optimized against globally distributed 10-year FAPAR and LAI time series from MODIS to reanalyse climatic drivers for vegetation phenology (Stöckli et al., 2008(Stöckli et al., , 2011)).This optimization resulted in a good representation of temporal FAPAR and LAI dynamics in all major biomes except evergreen tropical forests (Stöckli et al., 2011).Model parameters of the Biome-BGC model were optimized against eddy covariance flux observations and NDVI time series from MODIS for poplar plantations in northern Italy which resulted in a more accurate representation of carbon fluxes and NDVI (Migliavacca et al., 2009).The BETHY-CCDAS model was optimized against FAPAR time series from MERIS for seven eddy covariance sites (Knorr et al., 2010) and later for 170 land grid cells using coarse 8 by 10 • spatial resolution (Kaminski et al., 2012).These studies demonstrated the improvements in simulated vegetation phenology by optimizing model parameters against observations of vegetation greenness.
Nevertheless, spatial patterns and temporal dynamics of vegetation greenness were not yet optimized in a DGVM globally at a higher spatial resolution (0.5 • ) and by using long-term (30 year) satellite-derived time series of vegetation greenness.Newly developed 30 year time series of LAI or FAPAR from the GIMMS3g data set (Global Inventory Modeling and Mapping Studies, third generation of data sets; Zhu et al., 2013) make it possible to improve DGVMs not only based on seasonal cycles of single years (i.e.phenology) but additionally against decadal time series properties including inter-variability and trends.By integrating the GIMMS3g FAPAR data set in a DGVM, we can potentially improve spatial patterns and seasonal to long-term temporal dynamics of vegetation greenness.We use the LPJmL DGVM (Lund-Potsdam-Jena managed lands).Similar to other DGVMs, LPJmL does not accurately reproduce the growing season onset and seasonal amplitude of observed LAI and FA-PAR time series, presumably because of a limited phenology model (Kelley et al., 2013;Murray-Tortarolo et al., 2013).Thus integrating long-term observations of FAPAR in the LPJmL DGVM potentially requires the development of an improved phenology scheme.
We aim to improve environmental controls on vegetation phenology and greenness in LPJmL by (1) developing a new phenology module for LPJmL, by (2) optimizing FAPAR, productivity and phenology-related parameters of LPJmL against 30-year satellite-derived time series of FA-PAR, against 10-year satellite-derived time series of vegetation albedo and against spatial patterns of mean annual gross primary production (GPP) from a data-oriented estimate and by (3) integrating further data streams into LPJmL to constrain land cover dynamics and disturbance effects on vegetation greenness in diagnostic model simulations.This modeldata integration approach for LPJmL (LPJmL-MDI) will be applied to identify the environmental controls on vegetation greenness phenology.

Overview
LPJmL is a dynamic global vegetation model that simulates ecosystem processes such as carbon and water fluxes, carbon allocation in plants and soils, permafrost dynamics, fire spread and behaviour and vegetation establishment and mor-tality.We used LPJmL version 3.5.This version is based on the original LPJ model (Sitch et al., 2003).The model has been extended for human land use (Bondeau et al., 2007), and agricultural water use (Rost et al., 2008).It includes a process-oriented fire model (Thonicke et al., 2010), an improved representation of surface albedo and snow coverage (Strengers et al., 2010) and a newly implemented soil hydrology scheme and permafrost module (Schaphoff et al., 2013).This study focuses on the natural vegetation plant functional types (PFTs) (Sitch et al., 2003), i.e. our model developments and optimizations were not applied for crop functional types (CFTs) (Bondeau et al., 2007) because crop phenology is highly driven by human practices.
We developed a model-data integration approach for the LPJmL DGVM (LPJmL-MDI, Fig. 1).LPJmL-MDI allows us (1) to directly insert land cover, tree cover and burnt area data sets in LPJmL for diagnostic model applications (Sect.2.4.1);(2) to optimize LPJmL model parameters against data sets (here FAPAR, GPP, albedo; Sect.2.4.2); and (3) to evaluate and benchmark LPJmL simulations against observations or observation-based data sets (Sect.2.4.3).Like in a prognostic mode, LPJmL was driven by climate forcing data.Additionally, observed burnt areas were directly inserted into LPJmL to consider observed fire dynamics in diagnostic model applications.For this, we directly replaced the simulated burnt area in the LPJmL-SPITFIRE fire module (Thonicke et al., 2010) by observed burnt areas using the approach of Lehsten et al. (2008).Thus, the timing and location of fire spread is constrained by observations whereas fire effects on vegetation are still simulated by LPJmL-SPITFIRE.We further prescribed observed land cover and tree cover fractions to control for vegetation dynamics in parameter optimization experiments.Observed FAPAR and albedo time series as well as observation-based mean annual spatial patterns of GPP were used in a joint cost function to optimize productivity, phenology, radiation, and albedo-related model parameters using a genetic optimization algorithm.

FAPAR
FAPAR is defined as the ratio between the photosynthetic active radiation absorbed by the green canopy (APAR) and the total incident photosynthetic active radiation (PAR).Thus, the total FAPAR of a grid cell is the sum of FAPAR that is distributed among the individual PFTs: where n is the number of established PFTs in a grid cell.
The FAPAR of a PFT depends on the annual maximum foliar projective cover (FPC), on the daily snow coverage in the green canopy (F snow, gv ), green-leaf albedo (β leaf ) and the daily phenology status (Phen): Thus, the temporal dynamic of FAPAR in LPJmL is affected on an annual time step by changes in foliar projective cover (FPC PFT ) and on daily time steps by changes in phenology (Phen PFT ) and snow coverage in the green canopy (F snow, gv, PFT ; Supplement Fig. S1).This approach extends the previous implementation of Sitch et al. (2003) in which FAPAR depends only on FPC and phenology but leaf albedo and snow effects on FAPAR were missing.FPC PFT expresses the land cover fraction of a PFT.It is the annual maximum fractional green canopy coverage of a PFT and is annually calculated from crown area (CA), population density (P ) and LAI (Sitch et al., 2003): The last term expresses the light extinction in the canopy which depends exponentially on LAI and the light extinction coefficient k of the Lambert-Beer law (Monsi and Saeki, 1953).The parameter k had a constant value of 0.5 for all PFTs in the original LPJmL formulation (Sitch et al., 2003).We changed k to a PFT-dependent parameter because it varies for different plant species as seen from field observations (Bolstad and Gower, 1990;Kira et al., 1969;Monsi and Saeki, 1953).Crown area and LAI are calculated based on allocation rules and are depending on the annual biomass increment (Sitch et al., 2003).Population density depends on establishment and mortality processes in LPJmL (Sitch et al., 2003).

Phenology
The daily phenology and green leaf status of a PFT (Phen PFT ) in LPJmL expresses the fractional cover of green leafs (from 0 = no leafs to 1 = full leaf cover).Thus, it represents the temporal dynamic of the canopy greenness.We explored two phenology models in this study: First, we were trying to optimize model parameters of the original phenology module in LPJmL (LPJmL-OP, Sitch et al., 2003).Secondly, we implemented a new phenology module based on the growing season index (GSI) concept (Jolly et al., 2005), hereinafter called LPJmL-GSI.
LPJmL-OP has three different routines for summergreen (i.e.temperature-driven deciduous), evergreen (no seasonal variation) and raingreen (i.e.water-driven deciduous) PFTs (details in Supplement 1.1).Obviously, LPJmL-OP misses important controls on phenology, like effects of light in all PFTs or effects of water in summergreen and herbaceous PFTs.Additionally, in herbaceous PFTs the end of the growing season is not controlled by environmental conditions but is defined based on fixed calendar dates.
Because of the obvious limitations of LPJmL-OP, we developed the alternative LPJmL-GSI phenology module.The growing season index (GSI) is an empirical phenology model that multiplies limiting effects of temperature, day length and vapour pressure deficit (VPD) to a common phenology status (Jolly et al., 2005).We modified the GSI concept for the specific use in LPJmL (LPJmL-GSI).We defined the phenology status as a function of cold temperature, short-wave radiation and water availability.Additionally to the original GSI model, we added a heat stress limiting function because it has been suggested that vegetation greenness is limited by temperature-induced heat stress in several ecosystems (Bunn et al., 2007;Verstraeten et al., 2006) and has been demonstrated that heat stress reduces plant productivity also without additional water stress (Jiang and Huang, 2001;Van Peer et al., 2004;Poirier et al., 2012).Thus, the daily phenology status of a PFT is the product of the daily cold temperature (f cold, PFT ), light (f light, PFT ), water (f water, PFT ) and heat stress (f heat, PFT ) limiting functions: Examples for the four functions are shown in Fig. 2. The cold temperature limiting function at a daily time step t is defined as where sl cold, PFT and base cold, PFT are PFT-dependent slope and inflection point parameters of a logistic function based on mean daily air temperature T ( • C).The parameter τ cold, PFT is the change rate parameter based on the difference between the actual predicted limiting function value and the previous-day cold temperature limiting function value.This parameter introduces a temporal autocorrelation in the phenology status and avoids abrupt phenological changes because of changing weather conditions.The light-limiting function was implemented accordingly: where sl light, PFT and base light, PFT are the PFT-dependent slope and inflection point parameters of a logistic function based on daily shortwave downward radiation SW (W m −2 ).The parameter τ light, PFT is the temporal change rate for the light-limiting function.
The water-limiting function f water, PFT depends on the daily water availability W ( %) in LPJmL: where sl water, PFT and base water, PFT are the PFT-dependent slope and inflection point parameters of a logistic function based on daily water availability.W is a ratio between water supply from soil moisture and atmospheric water demand (Supplement 1.2) (Gerten et al., 2004).The parameter τ water, PFT is the temporal change rate for the water-limiting function.
The heat stress limiting function is defined as the coldtemperature limiting function based on daily air temperature but with a negative slope parameter: where sl heat, PFT and base heat, PFT are the PFT-dependent slope and inflection point parameters of a logistic function based on T .The parameter τ heat, PFT is the temporal change rate for the heat limiting function.
Besides the additional use of the heat stress limiting function, LPJmL-GSI has important differences to the original GSI phenology model (Jolly et al., 2005).We made the water limiting function dependent on water availability.VPD has been used instead in the original GSI phenology model.Nevertheless, it has been shown that phenology is more driven by soil moisture and plant available water than by atmospheric water demand especially in Mediterranean and grassland ecosystems (Archibald and Scholes, 2007;Kramer et al., 2000;Liu et al., 2013;Yuan et al., 2007) and that GSI performed better when using a soil moisture limiting function instead of the VPD limiting function (Migliavacca et al., 2011).With the implementation of the water limiting function in LPJmL-GSI, phenology depends not only on atmospheric water demand as in the original GSI model but also on water supply from soil moisture.Additionally, the soil moisture can be modulated through seasonal freezing and thawing in permafrost soils according to the permafrost routines in LPJmL (Schaphoff et al., 2013).Another important difference to the original GSI phenology model is the use of logistic functions instead of stepwise linear functions with fixed thresholds because smooth functions are generally easier to optimize than functions with abrupt thresholds and potentially better represent biological processes.A moving average of 21 days has been used in the original GSI model to create smooth phenological cycles and to avoid abrupt phenology changes because of daily weather variability (Jolly et al., 2005).It has been shown that PFT-and limiting functiondependent time-averaging parameters are needed instead of one single time averaging parameter (Stöckli et al., 2011).We implemented change rate parameters τ cold , τ light , τ water and τ heat that are PFT-and limiting function-dependent instead of moving average window lengths because LPJmL cannot use the same running window time-averaging approach as a prognostic model.

Data sets for parameter optimization: FAPAR, albedo and GPP
We used FAPAR, albedo and GPP data sets to optimize phenology, FAPAR, productivity and vegetation albedo-related parameters in LPJmL (Fig. 2).We require long-term FAPAR data sets to improve vegetation greenness in LPJmL on seasonal to decadal timescales.Two recently developed data sets provide 30-year time series of FAPAR.The Geoland2 BioPar (GEOV1) FAPAR data set (Baret et al., 2013) (hereinafter called GL2 FAPAR) and the GIMMS3g FAPAR (Zhu et al., 2013) data sets were used in this study.GL2 FAPAR is defined as the black-sky green canopy FA-PAR at 10:15 solar time and has been produced based on SPOT VGT (1999VGT ( -2012) ) and AVHRR (1981AVHRR ( -2000) ) observations (Baret et al., 2013).The GL2 FAPAR data set has a temporal resolution of 10 days and a spatial resolution of 0.05 • for the AVHRR-period and of 1/112 • for the SPOT VGT period.GIMMS3g FAPAR corresponds to black-sky FAPAR at 10:35 solar time and has been produced based on the GIMMS3g NDVI data set (Pinzon and Tucker, 2014;Zhu et al., 2013).GIMMS3g FAPAR has a 15-day temporal resolution and a 1/12 • spatial resolution and covers July 1981 to December 2011.We excluded in both FAPAR data sets observations that were flagged as contaminated by snow, aerosols or clouds.Additionally, we excluded FAPAR observations for months with temperatures < 0 • C to exclude potential remaining distortions of snow cover.Both data sets were aggregated to a 0.5 • spatial and monthly temporal resolution to be comparable with LPJmL simulations.We found that the GL2 AVHRR and GL2 VGT FAPAR data sets have not been well harmonized (Supplement 2.1).Thus, we did not use the combined GL2 VGT and AVHRR FAPAR data set for parameter optimization and for analyses of interannual variability and trends but only for analyses and evaluations of mean seasonal cycles and spatial patterns of FA-PAR.The GIMMS3g FAPAR data set has no uncertainty estimates.Uncertainty estimates are necessary in multiple data stream parameter optimization to weight single data streams in the total cost function.As a workaround we estimated the uncertainty based on monthly varying quantile regressions to the 0.95 quantile between FAPAR and the FAPAR uncertainty in the GL2 VGT data set.We applied the fitted regressions to the GIMMS3g data set to estimate FAPAR uncertainties (Supplement 2.2).The fit to the upper quantile provides conservative uncertainty estimates for the GIMMS3g FAPAR data set.
We used monthly shortwave white-sky albedo time series ranging from 2000-2010 from the MODIS C5 data set (Lucht et al., 2000;Schaaf et al., 2002) to constrain vegetation albedo parameters.Albedo observations in months with < 5 • C air temperature and above an albedo of 0.3 were excluded from optimization because we are optimizing only vegetation-related albedo parameters.High albedo values at low temperatures are probably affected by changing snow regimes which is not within our focus of model development and optimization.Thus we are only optimizing growing season albedo.
We used mean annual total GPP patterns from the dataoriented MTE (model tree ensemble) GPP estimate (Jung et al., 2011).This GPP estimate uses FLUXNET eddy covariance observations together with satellite observations and climate data to upscale GPP using a machine learning approach (Jung et al., 2011).This data set is not an observation but a result of an empirical model.Nevertheless, evaluation and cross-validation analyses have shown that this data set well represents the mean annual spatial patterns and mean seasonal cycles of GPP whereas it has a poor performance in representing temporal GPP anomalies (trends and extremes) (Jung et al., 2011).Thus, we are only using the mean annual total GPP from this data set for parameter optimization to constrain LPJmL within small biases of mean annual GPP.We used the mean seasonal cycle from the MTE GPP product as an independent benchmark for model evaluation.

Data sets for the prescription of land cover, tree cover and burnt area
The FAPAR, albedo and GPP data sets do not presumably contain enough information to constrain all processes that control FAPAR dynamics -especially processes like establishment, mortality, competition between PFTs, allocation and disturbances control FPC and thus FAPAR.The optimization of parameters of these processes against appropriate data streams is not feasible within this study.Thus, we directly prescribed land and tree cover fractions as well as burnt areas from observed data to control for some of these processes.
To prescribe land and tree cover in LPJmL, we combined several data sets to create observation-based maps of FPC (Supplement 3.1).Land cover maps from remote sensing products are not directly comparable with PFTs in global vegetation models due to differences in classification systems (Jung et al., 2006;Poulter et al., 2011a).PFTs in LPJmL are defined according to biome (tropical, temperate or boreal), leaf type (needle-leaved, broadleaved) and phenology type (summergreen, evergreen, rain green).We extracted the biome information from the Köppen-Geiger climate classification (Kottek et al., 2006) whereas leaf type and phenology were extracted from the SYNMAP land cover map (Jung et al., 2006).FPC was derived from MODIS tree cover (Townshend et al., 2011).Because LPJmL so far classified herbaceous vegetation according to their photosynthetic pathway (i.e.C 3 , temperate herbaceous and C 4 , tropical herbaceous), we further subdivided herbaceous PFTs according to biome and introduced a polar herbaceous PFT (PoH) based on the existing temperate herbaceous PFT (TeH) to differentiate tundra from temperate grasslands.
Burnt area data were prescribed directly in LPJmL by combining three data sets, the Global Fire Emissions Database (GFED) burnt area data set (Giglio et al., 2010), the Alaska Large Fire Database (ALFDB) (Frames, 2012;Kasischke et al., 2002) and the Canadian National Fire Database (CNFDB) (CFS, 2010;Stocks et al., 2002).GFED provides monthly burnt area estimates in 0.5 • resolution from 1996-2011.Burnt areas from the Alaska (ALFDB) and Canada (CNFDB) fire databases were used to extend burnt area time series before 1996 for boreal North America.Fire perimeter observations from 1979-1996 from ALFDB and CNFDB were aggregated to 0.5 • × 0.5 • gridded monthly burnt area time series.Observations before 1979 were excluded because fires were not reported for all provinces in Canada.Although the CNFDB contains only fire perimeters > 200ha, in both databases some fires are missing due to different mapping techniques, and fire perimeters do not agree with burned area, the integration of these data sets provides unique information about spatial-temporal patterns of disturbances especially in boreal ecosystems.It is necessary to simulate fire activity also during the model spin-up as fire influences the equilibrium between vegetation, soil and climate as well.Otherwise biomass would be overestimated at the beginning of the transient model run.For this purpose, we created artificial burnt area time series for the periods 1901-1978 (North America) and 1901-1995 (rest of the world).For this, observed annual total burnt areas from the periods 1979-2011 (North America) and 1996-2011 (rest of the world) were resampled according to temperature and precipitation conditions and assigned to the pre-data period in order to include fire regimes that agree with observed fire regimes in the spinup of LPJmL.This approach assumes that fire regimes in the pre-data period were not different than in the observation period.

Data sets for model evaluation
LPJmL was evaluated against data sets that are independent of the optimization and prescription data sets and against independent temporal or spatial scales of the optimization and prescription data sets.We compared LPJmL against mean annual patterns and mean seasonal cycles of ET from the MTE estimate (Jung et al., 2011).Further, we evaluated model results against spatial patterns of biomass.Ecosystem biomass estimates were taken from satellite-derived forest biomass maps for the tropics (Saatchi et al., 2011) and for the temperate and boreal forests (Thurner et al., 2014) including an estimation of herbaceous biomass (Carvalhais et al., 2014).Additionally, we evaluated LPJmL against independent temporal and spatial scales of the integration data (mean seasonal cycle of GPP, tree cover, inter-annual variability and trends of FAPAR).We were using tree cover from MODIS to evaluate LPJmL model runs with dynamic vegetation.

Climate forcing data and model spin-up
LPJmL was driven by observed monthly temperature and precipitation data from the CRU TS3.1 data set ranging from 1901-2011 (Harris et al., 2013) as well as by monthly shortwave downward radiation and long-wave net radiation reanalysis data from ERA-Interim (Dee et al., 2011).
LPJmL needs a model spin-up to establish PFTs and to bring vegetation and soil carbon pools into equilibrium.The spin-up was performed according to the standard LPJmL modelling protocol (Schaphoff et al., 2013;Thonicke et al., 2010): LPJmL was run for 5000 years by repeating the climate data from 1900-1930.After the spin-up model run, the transient model run was restarted from the spin-up conditions in 1901 and LPJmL was run for the period 1901-2011.Model results were analysed for the observation period .
For model optimization experiments we used a different spin-up scheme because the spin-up is computation time demanding and many model runs are needed during optimization experiments.As in the standard modelling protocol, we firstly spin-up the model for 5000 years by repeating the climate from 1901-1930.Secondly, a transient model run was restarted from the spin-up conditions in 1901 and was performed for the period 1901-1979.Thirdly, each optimization experiment was restarted from the conditions in 1979 and a second spin-up for 100 years by recycling the climate from 1979-1988 was performed.The transient model run was restarted from the conditions of the second spin-up and simulated for the period 1979-2011.This second spin-up is needed to bring the vegetation into a new equilibrium which can be caused by a new parameter combination during optimization.From visual analyses of model results, we found that a spin-up time of 100 years for the second spin-up was enough to eliminate trends in FAPAR and GPP that resulted from other equilibrium conditions.

Prescription of land and tree cover
Land cover is expressed as FPC in LPJmL.We used the observation-based FPC data set to prescribe land and tree cover in LPJmL (Sect.2.3.2,Supplement 3.1).The presence of a PFT in a grid cell depends on establishment and mortality in LPJmL (Sitch et al., 2003).A PFT establishes in a grid cell if the climate is within the bioclimatic limits of the PFT for establishment and survival.On the other hand, a PFT dies in a grid cell if the climate is no longer suitable for the PFT.Additionally, mortality occurs because of heat stress, low productivity, competition among PFTs for light, and because of fire disturbance (Sitch et al., 2003;Thonicke et al., 2010).
FPC is the major variable that contributes to interannual variability of FAPAR in LPJmL despite the daily www.biogeosciences.net/11/7025/2014/Biogeosciences, 11, 7025-7050, 2014 phenological status.Thus fixing FPC to the observed value is not a desired solution to prescribe land cover in LPJmL.
Fixing FPC would neglect mortality effects on land cover but would also permit the simulation of post-fire succession trajectories.Consequently, we prescribed land cover in LPJmL using a hybrid diagnostic-dynamic approach.In this approach we prescribed the annual maximum FPC in LPJmL similar to previous approaches (Poulter et al., 2011b).Firstly, we switched off the effects of bioclimatic limits on establishment and mortality.Only these PFTs were allowed to establish in a grid cell that occurred in the observed land cover data set.Vegetation growth depends on the annual biomass increment and allocation rules in LPJmL.This leads to an extension of FPC of each PFT.We limited a further expansion of FPC if the simulated FPC exceeded the observed FPC by replacing the simulated FPC with the observed FPC (prescribed maximum FPC).Consequently, the simulated FPC can be lower than the observed FPC because the PFT is still growing or because the FPC was reduced due to fire, heat stress or low productivity.For herbaceous PFTs we only reduced the FPC if the observed total fractional vegetation cover in a grid cell was exceeded.This allowed herbaceous PFTs to replace tree PFTs if the FPC of trees was reduced due to fire or other mortality effects in the model.With this approach a prescription of land cover can be achieved in LPJmL which well represents observed PFT distributions (Supplement 3.2) but still allows for main processes of dynamic vegetation.

Parameter optimization
Photosynthesis, albedo, FAPAR and phenology-related model parameters of LPJmL were optimized against observed FAPAR and albedo satellite observations and dataoriented estimates of GPP.A description of all parameters including parameter values is given in Supplement 4.1.The parameter α a is the most important parameter in LPJmL for photosynthesis (Zaehle et al., 2005).This parameter accounts for the amount of radiation that is absorbed at leaf level in comparison to the total canopy.Thus, this parameter is a replacement for a more enhanced model formulation for canopy structure and leaf clumping.We used this parameter to adjust biases in GPP.The PFT-dependent leaf, stem and litter albedo parameters (β leaf , β stem and β litter ) are mostly sensitive for model simulations of albedo.The parameter β leaf affects additionally the maximum FAPAR of a PFT.The light extinction coefficient k controls the FPC of a PFT and thus affects mainly land cover, maximum FAPAR and the available radiation for photosynthesis.All other parameters that were considered in optimization experiments are the parameters of the LPJmL-OP and LPJmL-GSI phenology modules.These parameters contribute mainly to seasonal variations in FAPAR.Some parameters were excluded from optimization experiments that were identified as insensitive to GPP and FAPAR simulations in PFTs.The temporal change rate pa-rameters τ cold , τ light , τ heat and τ water are insensitive in most PFTs because of the monthly temporal resolution of the climate forcing data used.The optimization of model parameters was performed by minimizing a cost function between model simulations and observations using a combined genetic and gradientbased optimization algorithm (GENOUD, genetic optimization using derivatives, Mebane and Sekhon, 2011, see Supplement 4.2 for details).The cost function J of LPJmL for a single model grid cell (gc) depends on the scaled model parameter vector d (d = proposed parameter value/prior parameter value) and is the sum of square error (SSE) between model simulation and observation weighted by the number of observations (nobs) for each data stream (DS): The SSE for a single data stream is calculated from the LPJmL simulation of this data stream (x LPJmL ) and the corresponding observed values (x obs ) weighted by the uncertainty of the observations (x obsunc ) for each time step t: where p 0 are LPJmL prior parameters.That means that the minimization of the cost function J is based on scalars of LPJmL parameters relative to the prior parameter values.Different model optimization experiments were performed for individual grid cell and for multiple grid cells of the same PFT for LPJmL-OP as well as for LPJmL-GSI (Table 1).In the grid-cell-based optimization experiments model parameters of the established target tree PFT and the established herbaceous PFT were optimized at the same time.The purpose of grid cell level optimization experiments was to explore the variability of parameters within different regions and PFTs.In the PFT level optimization experiments the cost of LPJmL was calculated as the sum of the cost for each grid cell weighted by the grid cell area A: For PFT level optimizations parameters of herbaceous PFTs were first optimized for grid cells where only the herbaceous PFT was dominant.In a second step, the optimized parameters of the herbaceous PFTs were used in the optimization of the target tree PFT (Fig. S9).The purpose of PFT level optimization experiments is to derive optimized parameter sets that can be used for one PFT in global model runs.
For grid cell as well as PFT level optimization experiments, we only used grid cells that are vegetated, dominated by one PFT and that are only marginally affected by agricultural use or fire disturbances.These grid cells are called candidate grid cells in the following.We randomly selected grid cells from the set of candidate grid cells to perform grid cell-or PFT level optimization experiments.Table 1 gives an overview of all optimization experiments for LPJmL-OP and LPJmL-GSI with the number of used grid cells.Grid cells that were selected for optimization experiments are also shown in Fig. 3.The PFT level optimization of LPJmL-OP (OP.pft) did not result in plausible posterior parameter sets because of structural limitations of the LPJmL-OP phenology model for herbaceous PFTs (i.e.no water effects, calendar day as end of growing season), raingreen PFT (i.e.binary phenology) and evergreen PFTs (i.e.constant phenology) and was therefore excluded from further analysis.
Posterior parameter sensitivities, uncertainties and correlations were explored by analysing the maximum likelihood and the posterior range of each parameter as derived from all parameter sets from the genetic optimization algorithm (Supplement 4.3).

Model evaluation and time series analysis
Global model runs of LPJmL were performed in order to evaluate model results against the integration data, against independent metrics of the integration data and against independent data streams.We evaluated results from LPJmL-OP with standard parameters (LPJmL-OP-prior), from LPJmL-OP with optimized productivity, albedo and FAPAR parameters from grid cell level optimization experiments (LPJmL-OP-gc) and from LPJmL-GSI with optimized parameters from PFT level optimization experiments (Table 2).We did not use optimized phenology parameters in the LPJmL-OPgc model run because we were not able to derive plausible phenology parameters in optimization experiments of LPJmL-OP.All model runs were performed with dynamic vegetation and prescribed burnt areas.
We aggregated monthly FAPAR time series to mean annual FAPAR to evaluate inter-annual variability and trends.Mean annual FAPAR time series were averaged from all monthly values with mean monthly air temperatures > 0 • C to exclude potential remaining effects of snow in the observed FAPAR time series.Trends in mean annual FAPAR time series and trend breakpoints were computed using the "greenbrown" package for the R software (Forkel et al., 2013).In this implementation, trends are computed by fitting piece-wise linear trends to the annual FAPAR time series using ordinary least squares regression.The significance of trends was computed using the Mann-Kendall trend test (Kendall, 1975;Mann, 1945).Parameters as in GSI.prior (Table S4) Table S5 (one optimized parameter set per PFT) 3 Results and discussion

Performance of phenology models
The newly developed LPJmL-GSI phenology model resulted in significantly higher correlations with monthly GIMMS3g FAPAR than LPJmL-OP in all PFTs except in the tropical broadleaved evergreen (TrBE) and boreal broadleaved summergreen (BoBS) PFTs (Fig. 4).LPJmL-OP with prior parameters had high correlations with monthly GIMMS3g FAPAR in broadleaved summergreen PFTs (TeBS median r = 0.87, BoBS median r = 0.92) and medium correlations in boreal needle-leaved PFTs (BoNE median r = 0.53, BoNS median r = 0.6).In all other PFTs, LPJmL-OP had low correlations with monthly GIMMS3g FAPAR.The correlation against monthly GIMMS3g FAPAR did not significantly improve in all PFTs after grid cell level optimization experiments of LPJmL-OP (Fig. 4).The use of the newly developed LPJmL-GSI phenology model already significantly improved the correlation with monthly GIMMS3g FAPAR in all PFTs except in the temperate herbaceous (TeH) and BoBS PFTs.LPJmL-GSI had significantly higher correlations with monthly GIMMS3g FAPAR after grid cell level optimization experiments in the TrBR, TeNE, TeBS, TeH, BoBS and BoNS PFTs.After PFT level optimization experiments, LPJmL-GSI had median correlation coefficients > 0.5 in all PFTs except in broadleaved evergreen PFTs (TrBE, TeBE).These results prove that the raingreen, evergreen and herbaceous phenology schemes of LPJmL-OP were not able to reproduce temporal FAPAR dynamics despite the attempt of parameter optimization and that LPJmL-GSI can reproduce seasonal FAPAR dynamics in most PFTs.The low correlation coefficients between LPJmL-GSI and GIMMS3g FAPAR after optimization experiments in broadleaved evergreen PFTs (TrBE, TeBE) might be caused by the specific properties of the FAPAR data set in these PFTs.GIMMS3g FAPAR does not have a clear seasonal cycle but a high short-term variability in broadleaved evergreen forests.These regions are often covered by clouds that inhibit continuous optical satellite observations.The high short-term variability results ultimately in low correlation coefficients between both LPJmL versions (LPJmL-OP and LPJmL-GSI) and GIMMS3g FAPAR time series.In temperate broadleaved evergreen forests, the GIMMS3g FAPAR data set might have a wrong seasonality.In these regions, the mean seasonal FA-PAR cycles from the GIMMS3g and GL2 VGT FAPAR data sets are anti-correlated and FAPAR from LPJmL-GSI agrees   S3) LPJmL-GSI GSI-based phenology Parameters from the GSI.pft optimization experiment (Table S5) better with the GL2 VGT data set.Because of these reasons, we did not expect to improve seasonal FAPAR dynamics in broadleaved evergreen forests with the current model-data integration setup.
All optimization experiments of LPJmL-OP and LPJmL-GSI resulted in a significant reduction of the cost in comparison to the respective prior models (Supplement 4.4, Fig. S10).Nevertheless, the prior parameter set of LPJmL-GSI resulted already in a significant lower cost than the grid cell level optimized parameter sets of LPJmL-OP in tropical and polar herbaceous PFTs, and in temperate broadleaved summergreen and boreal needle-leaved summergreen PFTs.The reduction of the overall cost was in all model optimization experiments usually associated with a significant reduction of the annual GPP bias (Fig. S11).LPJmL-OP with prior parameters underestimated mean annual GPP in the tropical broadleaved evergreen PFT and overestimated mean annual GPP in all other PFTs.Grid cell level optimization experiments of LPJmL-OP resulted in a significant reduction of the GPP bias in all PFTs except in the polar herbaceous PFT (PoH).We were not able to remove the GPP bias and to reduce the cost of LPJmL-OP and of LPJmL-GSI in the PoH PFT in optimization experiments because of inconsistencies between the FAPAR and GPP data sets or in the LPJmL formulation.LPJmL was not able to sustain the relatively high peak FAPAR in tundra regions as seen in the GIMMS3g data set given the low mean annual GPP of the MTE data set (Supplement 4.4).These inconsistencies might be related to higher uncertainties of the GPP and FAPAR data sets in tundra regions where the MTE GPP data set is not covered by many eddy covariance measurement sites, and where satellite-based FAPAR observations are affected from high sun zenith angles (Tao et al., 2009;Walter-Shea et al., 1998).On the other hand, dominant tundra plant communities like mosses and lichen are not represented in LPJmL (Supplement 4.4).All model optimizations experiments kept growing season albedo within reasonable ranges in comparison to MODIS albedo (Fig. S12).These results demonstrate an improved performance of optimized model parameter sets over prior model parameter sets and of LPJmL-GSI over LPJmL-OP regarding a cost that is defined based on 30 years of monthly FAPAR, mean annual GPP and 10 years of monthly vegetation albedo.

Parameter sensitivities and uncertainties
The uncertainty of productivity and albedo-related parameters was reduced after optimization of LPJmL-GSI in most PFTs while the reduction of the uncertainty of phenologyrelated parameters depended often on plant functional type (Fig. 5).Prior and posterior parameter values from each optimization experiment are listed in the Supplement (Tables S2  to S5).
The parameter α a (absorption of light at leaf level in relation to canopy level) was sensitive within a narrow parameter range for all PFTs.The posterior α a parameter range was smaller than the uniform prior range in all PFTs.In all optimization experiments we found for the parameter α a a gradient from high values in tropical to low values in boreal PFTs (Fig. S13).This pattern reflects the initial overestimation of mean annual GPP in temperate and boreal PFTs and underestimation of GPP in tropical regions with the prior parameter set of LPJmL-OP.Thus, the low α a parameter values probably account for nitrogen limitation effects on productivity in boreal forests (Vitousek and Howarth, 1991) that are currently not considered in LPJmL.A future implementation of nitrogen limitation processes in LPJmL requires a re-optimization of the α a parameter.
The leaf albedo parameter β leaf was sensitive in all PFTs and the posterior β leaf parameter range was smaller than the prior parameter range in evergreen PFTs.In these evergreen PFTs the β leaf parameter was well constrained because albedo satellite observations are less affected by variations in background albedo (soil, snow) than in deciduous PFTs.In all other PFTs the β leaf posterior parameter range was equal the prior parameter range or the optimized parameter value was close to a boundary of the prior parameter range.This result indicates that the albedo routines in LPJmL should consider variations in background albedo caused by changes in soil properties, soil moisture, or snow conditions in order to accurately reproduce satellite-observed albedo time series (see discussion in Supplement 4.5).Nevertheless, the optimization of the leaf albedo parameter β leaf resulted in values that differed especially between broadleaved and needleleaved evergreen PFTs as well as herbaceous PFTs (Fig. 5, Fig. S14).Low leaf albedo parameters in needle-leaved evergreen PFTs (TeNE and BoNE) and high leaf albedo parameters in broadleaved summergreen and herbaceous PFTs agree well with the patterns reported by Cescatti et al. (2012).
The light extinction coefficient k was sensitive for all PFTs but the posterior parameter range was only in herbaceous PFTs and in the BoBS PFT smaller than the prior parameter range (Fig. 5).In all PFTs this parameter had a large spatial variability (Fig. S15).The parameter k affects mostly the FPC and thus the maximum FAPAR.Thus, this parameter cannot be well constrained for tree PFTs in the current optimization setup because the maximum FPC of trees was prescribed from the land and tree cover data set.On the other hand, the maximum FPC of herbaceous PFTs was not prescribed from observations which resulted in narrow k posterior parameter ranges for herbaceous PFTs.The parameter k was optimized towards a very high value in the BoNS PFT (k = 0.7) due to high tree mortality rates after lowproductivity years (Supplement 4.5).This parameter would result in an overestimated PFT coverage in model runs with dynamic vegetation.Thus, we performed a second optimization experiment for this PFT (blue in Fig. 5) where k BoNS was limited to 0.65.This optimization experiment resulted in similar posterior values for the other parameters.Although the k parameter was well constrained for the TrH, TeH and PoH PFTs, these parameters cannot be used in the final parameter set of LPJmL-GSI.In dynamic vegetation model runs, the relatively low k parameter values for the TrH and TeH PFTs and relatively high values for the PoH PFT would result in an underestimation of herbaceous coverage in temperate and tropical climates and an overestimation of herbaceous coverage in boreal and polar climates, respectively.Therefore, we performed three more optimization experiments for herbaceous PFTs where we fixed k at 0.5 (blue in Fig. 5).These optimization experiments resulted in similar α a parameters but different albedo parameters and phenology parameters in order to compensate for biases in FAPAR and albedo that were introduced by the fixed k parameter.Thus, the high spatial variability and the large uncertainty of the light extinction coefficient k require re-addressing this parameter in a model optimization setup with dynamic vegetation using tree and vegetation cover data or perhaps a replacement by a better representation of canopy architecture and radiative transfer.The sensitivity and posterior uncertainty of phenologyrelated model parameters depended often on plant functional type.The parameter base cold which controls the effect of cold temperature on phenology was sensitive in all PFTs except the TrBE and TrH PFTs.The posterior parameter range was smaller than the prior parameter range in temperate PFTs (TeNE, TeBS and TeH).The parameter base heat which controls the effect of heat stress on phenology was sensitive in TrBR, TrH, TeH, BoNE and BoNS PFTs while in other PFTs this parameter was only sensitive towards the boundaries of the prior parameter range.Nevertheless, the posterior parameter range was only smaller than the prior parameter range in TrBR and TrH PFTs.The parameter base light was sensitive in temperate and boreal PFTs.In tropical PFTs this pa-rameter is only sensitive above a certain threshold (i.e. 60 W m −2 for TrBE and 100 W m −2 for TrBR).The parameter base water was sensitive in all PFTs.The posterior parameter range of this parameter was smaller in all PFTs except in TeBS, BoNE, BoBS and BoNS PFTs.Although the parameter base water had a large variability among PFTs, it was generally optimized towards higher values in PFTs that are presumably water controlled (TrBR, TeBS, TrH, TeH) and optimized towards lower values in PFTs that are presumably less water controlled (TrBE, TeNE, BoNE, BoNS, PoH).This result indicates that FAPAR of water-controlled PFTs reacts already to small decreases in water availability whereas other PFTs react only to strong decreases in water availability.We phenology-related model parameters (Fig. S16) which indicates the ability to disentangle the relative effects of temperature, light and water on phenology.As the base water parameter was the only phenology parameter which was sensitive in all PFTs, this indicates that water availability is the only phenological control that acts in all PFTs.

Effects on seasonal and inter-annual greenness dynamics
LPJmL-GSI represents better the observed spatial patterns and seasonal-to-decadal temporal dynamics of vegetation greenness (FAPAR) than LPJmL-OP (Fig. 6, Supplement 5.3).Whereas LPJmL-OP overestimated mean annual FAPAR in many high-latitude and semi-arid regions, LPJmL-GSI was closer to both data sets and within the uncertainty of the GL2 VGT FAPAR data set in most regions and under most climate conditions (Fig. S22).LPJmL-GSI still overestimated mean annual FAPAR in temperate dry regions, but this overestimation was reduced in comparison to LPJmL-OP.
We further observe a substantial improvement in LPJmL-GSI regarding the seasonal cycles, monthly and annual dynamics of FAPAR as retrieved from the GIMMS3g and GL2 VGT FAPAR data sets (Fig. 6, Figs.S23-S25).Monthly FA-PAR time series from LPJmL-GSI were significantly (p ≤ 0.05) higher correlated with GIMMS3g than from LPJmL-OP in boreal forests of eastern Siberia, in the North American tundra, in temperate and tropical grasslands of central Asia, North America, Australia and especially, in the Sahel (Fig. 6a).This is because of an improved representation of spring onset and the end of the growing season in temper-ate and boreal forests and in herbaceous PFTs (Fig. S24).The highest differences between simulated and observed mean seasonal FAPAR cycles were observed in the temperate broadleaved evergreen PFT, where both model versions had opposite, although insignificant, relationships to the GIMMS3g data sets.For this PFT, the observational constraints are particularly problematic, where a weak agreement and opposite relationship is observed between the two data sets (r = −0.48).
Globally, LPJmL-GSI describes better the inter-annual dynamics of GIMMS3g FAPAR when compared to the previous model versions (Fig. 6b).In 20 % of the land the difference to other model versions is statistically significant, and in 40 % does not detract from the previous model versions.This improvement in inter-annual variability is especially seen in temperate and tropical dry regions, with sparse tree cover and grassland dominated ecosystems (western United States, central Asia, the Sahel, southern Africa, and Australia) (Fig. S25).In the Arctic, boreal and temperate climates LPJmL usually shows a higher correlation with the GIMMS3g data set than the correlation observed between both data sets (GIMMS3g and GL2 VGT).These results demonstrate that LPJmL-GSI can explain the interannual variability of the GIMMS3g FAPAR data set especially in temperate and boreal forests and temperate and tropical grasslands.
Overall, the global spatial representation of phenological dynamics in LPJmL-GSI improves significantly over the previous model versions from seasonal to inter-annual timescales.Given the inclusion of water controls on phenological development, these results emphasize the importance of water availability in explaining the mean spatial patterns of vegetation greenness, but also the seasonal phenology as well as inter-annual dynamics in vegetation development.

Effects on trends in vegetation greenness
The role of different climate drivers underlying the greening and browning trends in vegetation activity is still highly debated and the dominant factors show a strong spatial variability (de Jong et al., 2013a).The consideration of different environmental controls on the phenological development in LPJmL shows a significant improvement in representing such dynamics when compared to the previous model formulations (Fig. 7).
Both LPJmL-OP and LPJmL-GSI reproduced the observed greening trends in tundra regions and in boreal forests of Siberia.In both model versions this greening is mostly driven by annual changes in foliar projective cover and effects of temperature on spring phenology.This agrees with observational studies that identified temperature increases as drivers for an increasing shrub cover in tundra ecosystems (Blok et al., 2011;Forbes et al., 2010;Myers-Smith et al., 2011;Raynolds et al., 2013;Sturm et al., 2001) and that found a positive association between warming, increasing tree ring widths and NDVI greening in boreal forests of eastern Siberia (Berner et al., 2011(Berner et al., , 2013)).
Parts of the boreal forests in North America showed significant browning trends in the GIMMS3g data set but a tendency to positive trends in the GL2 data set.The simulation results from LPJmL-GSI are in agreement with the GIMMS3g-based browning trends, rather than greening trends.However, these model-based browning trends were not as strong as in the GIMMS3g data set.In LPJmL-GSI these browning trends are caused by the effects of seasonal light and water effects on phenology, and by fire activity.In the GIMMS3g data set these browning trends were related to several environmental factors like fire activity (Goetz et al., 2005), temperature-induced drought stress (Beck et al., 2011 Bunn and Goetz, 2006) and to snow-regulated changes in soil water availability (Barichivich et al., 2014).
The Sahel had widespread greening trends in the GIMMS3g FAPAR data set.Whereas LPJmL-OP simulated browning trends, the implementation of water availability effects on phenology enabled LPJmL-GSI to reproduce the observed greening trends.Increases in precipitation and rainuse efficiency were also identified in observational studies as the main drivers of positive trends in vegetation greenness in the Sahel (Fensholt et al., 2013).
Overall, we observed that both LPJmL-OP and LPJmL-GSI reproduced the greening trends in tundra, boreal and temperate forests, although LPJmL-GSI showed a wider agreement in the extent of browning trends in the boreal forests of North America.Further, in the Sahel region, the greening trends can only be reproduced through the inclusion of water availability controls on the phenology development.These results demonstrate that environmental controls like light, heat stress and water availability contribute to a better description of regional greening and browning trends in very different bioclimatic regions of the globe.Hence, a comprehensive characterization of the different environmental controls on phenological development is essential in performing model-based analysis of long-term trends in vegetation activity.

Effects on carbon fluxes and stocks
LPJmL-GSI and LPJmL-OP-gc with optimized parameters represented better the global patterns and mean seasonal cycles of gross primary production and biomass than LPJmL with original phenology and prior parameters (LPJmL-OPprior) (Fig. 8).LPJmL-OP-prior overestimated mean annual GPP and biomass in most polar, boreal and temperate regions.LPJmL-OP-prior underestimated mean annual GPP but overestimated mean annual biomass in tropical regions around the equator.These biases were reduced in LPJmL-OP-gc and LPJmL-GSI.LPJmL generally overestimated GPP also in arid regions but these biases were reduced after optimization in LPJmL-OP-gc and LPJmL-GSI (Fig. S18).We also found that the mean seasonal cycle of GPP from LPJmL-GSI agreed better with the mean seasonal GPP cycle from the MTE estimate especially in temperate forests and in tropical, temperate and polar grasslands (Supplement 5.1, Fig. S17) although no information about the seasonality of GPP was included in optimization experiments.LPJmL-GSI still overestimated biomass in some tropical regions (African Savannas, southeast Brazil, south and southeast Asia) (Fig. S19).These regions were mainly simulated as agricultural lands in LPJmL, i.e. as different crop functional types (CFTs).The LPJmL-GSI phenology module was not applied or optimized for agricultural regions, where the seasonal phenological development is prescribed according to the CFTs parameterizations from Bondeau et al. (2007).Generally, LPJmL-GSI performed substantially better than LPJmL-OP-prior and LPJmL-OP-gc when comparing the global total carbon fluxes and stocks to the data-oriented estimates (Supplement 5.1, Table S6).These results demonstrate that in addition to the optimization of productivity parameters in LPJmL, the implementation of the new GSI-based phenology improved estimates of spatial patterns, seasonal dynamics, and global totals of gross primary production and biomass.

Effects on forest distribution
LPJmL-GSI with dynamic vegetation better represented the observed tree cover in high-latitude regions than LPJmL-OPprior and LPJmL-OP-gc (Fig. 8d).LPJmL-OP-prior highly overestimated tree cover in boreal and Arctic regions and simulated a too northern Arctic tree line in comparison with tree cover from MODIS observations.Although this overestimation was reduced after optimization, LPJmL-OP-gc still highly overestimated tree cover in boreal and temperate regions.The occurrence of trees was shifted southwards in LPJmL-GSI.Although LPJmL-GSI still overestimated tree cover in boreal regions, this overestimation was much lower than in LPJmL-OP-gc.LPJmL-OP-prior and LPJmL-OPgc slightly underestimated tree cover in temperate regions around 45 • N but this was well reproduced by LPJmL-GSI.We found no differences in tree cover between LPJmL-OP and LPJmL-GSI in other parts of the world where tree cover is highly affected from agricultural land use and thus implicitly prescribed to LPJmL.These results demonstrate that additional to the optimization of productivity parameters in LPJmL-OP-gc, the newly developed GSI-based phenology model and the optimized model parameters contribute to a better representation of tree cover in high-latitude regions.

Effects on evapotranspiration processes
Evapotranspiration from LPJmL agreed well with the dataoriented MTE estimates (Fig. 8b).The implementation and optimization of the new GSI-based phenology did not affect ET much.ET increased only in tropical rainforests around the equator in LPJmL-GSI and LPJmL-OP-gc in comparison to LPJmL-OP-prior because of the increased GPP in these regions.In other regions ET remained almost unchanged.But this does not imply that the structural improvements in LPJmL-GSI did not affect the transpiration processes (Figs.S20, S21).Indeed, LPJmL-GSI had lower interception losses than LPJmL-OP in boreal forests because of the reduced tree cover.On the other hand this implies that simulated soil evaporation was increased.Furthermore, interception and soil evaporation had slightly shifted seasonal cycles in LPJmL-GSI compared to LPJmL-OP due to the seasonal differences in timing of leaf development and senescence stages (Fig. S21).Consequently, small differences in total evapotranspiration result from opposite and compensatory changes in interception and soil evaporation and slight changes in transpiration fluxes in LPJmL-GSI.

Applicability and challenges of the LPJmL-GSI phenology module
The LPJmL-GSI phenology module is part of a DGVM that is applied for climate impact studies.In order to assess how well the model performs under different climate conditions, we additionally tested how the model performance changes in grid cells that were not used during optimization (Fig. S26).We found no general decrease in model performance with distance to the nearest grid cell used for optimization, or under different temperature conditions.Especially, no significantly lower correlations (p ≤ 0.05, Wilcoxon rank-sum test, Fig. S26) were found between simulated and observed FAPAR time series in grid cells that were 3 to 5 • C warmer than the closest optimization grid cell.
From a typical perspective of space for time substitution, this could indicate that the confidence in the simulation of FAPAR dynamics should not detract under climate warming scenarios of 0.3 to 4.8 • C (IPCC, 2014).Nevertheless, model optimization experiments and model evaluation indicated further needs for improvement in future studies -for example, simulations of surface albedo could improve through the implementation of time-varying effects of snow conditions and surface moisture on albedo.Also, an enhanced representation of canopy architecture and canopy radiative transfer could reduce the large spatial variability and parameter uncertainty found for the light extinction coefficient and hence improve the simulation of tree coverage and peak FAPAR.In addition to temperature, light and water availability, phenology also depends on other factors that are not considered in LPJmL-GSI.Phenology is also driven by leaf age (Caldararu et al., 2012(Caldararu et al., , 2014) ) and nutrient availability (Wright, 1996).These effects are neither considered in the original GSI phenology model (Jolly et al., 2005;Stöckli et al., 2011) nor in the LPJmL-GSI or other traditional formulations.Here, the lower posterior values found for the parameter α a may be compensating for missing nitrogen limitation effects on productivity in boreal forests (Vitousek and Howarth, 1991).Thus a future implementation of nitrogen limitation processes in LPJmL requires a re-optimization of the α a parameter.Additionally, the current implementation of phenology in LPJmL affects photosynthesis only through changes in APAR.In future model developments a stronger coupling between phenology and ecosystem carbon cycle dynamics could be explored.For example, the LPJmL-GSI phenology module could demand carbon for leaf development from photosynthesis or additional storage pools on the one hand and could trigger carbon turnover through litterfall on the other hand.In this case a phenology module could partly regulate an optimal carbon gain for a canopy similar to the approach of Caldararu et al. (2014).Nevertheless, such an analysis needs to go beyond the approach of Caldararu et al. (2014) and demands for additional observational constraints on ecosystem carbon fluxes, leaf area, biomass and litterfall.In order to better understand couplings between leaf phenology, changes in carbon allocation and photosynthesis it will be of benefit to use site level eddy covariance measurements from the FLUXNET network (Baldocchi et al., 2001) together with ancillary data in ecosystem-scale model optimization experiments (Carvalhais et al., 2010;Kuppel et al., 2012;Williams et al., 2009).Thus the LPJmL-GSI phenology module and the LPJmL model-data integration approach can serve as a framework to further explore hypotheses of ecosystem processes and vegetation dynamics.
We demonstrated the improved performance of LPJmL-GSI over LPJmL-OP in representing observed carbon fluxes and stocks, forest cover and seasonal to decadal dynamics of vegetation greenness.Thus, similar approaches to the LPJmL-GSI phenology module can be applied in other DGVMs to improve model simulations in comparison with observations.However, the adaptation of current results to other models should be cautionary because the phenology scheme of LPJmL-GSI is an empirical approach with PFTdependent parameters that need to be estimated.This estimation is model-specific because (1) different models do not necessarily use the same definition and set of PFTs; (2) our parameterizations depend on model structure, e.g.different models often use different hydrology routines; and (3) our posterior parameters for phenology were also constrained by using albedo and GPP data.Thus LPJmL-GSI model parameters cannot be easily transferred to other models.It might be possible to use the parameters of the temperature and light limiting functions in other models because these functions depend uniquely on the forcing data.On the other hand, the parameters for the water availability limiting function might need to be re-optimized because of differences in soil moisture computations.However, depending on the co-variability between forcing variables and simulated water availability by other models, the best parameterizations may differ from the ones presented here.Consequently, we acknowledge the potential need to optimize parameters of the LPJmL-GSI phenology model in order to obtain plausible results in other modelling structures.However, it is likely that the LPJmL-GSI phenology model can be easily applied to other models of the LPJ group of models (Prentice et al., 2011;Smith et al., 2001) that are using the hydrology routines of Gerten et al. (2004) while probably additional parameter optimization exercise are required to adapt the model to other types of DGVMs or ecosystem models.

Environmental controls on vegetation greenness phenology
As the newly developed GSI-based phenology model of LPJmL can reproduce the seasonality and monthly dynamics of observed FAPAR in most biomes, it can be used to identify environmental controls on vegetation greenness phenology.The importance of phenological controls differed by climate regions, ecosystems and season (Fig. 9).We identified environmental controls on seasonal FAPAR dynamics by analysing the mean seasonal cycles of FAPAR, of the cold temperature, light, water availability and heat stress limiting functions for phenology from the LPJmL-GSI model run.This analysis is comparable to previous investigations of limiting factors for vegetation phenology (e.g.Jolly et al., 2005;Caldararu et al., 2014).FAPAR seasonality in high-latitude regions (tundra, boreal forests) was mainly controlled by cold temperature (entire year) and light (October to February).We also found an important control by water availability in February to April in the tundra and in boreal forests of North America and eastern Siberia.This water limitation in early spring was due to the seasonal freezing of the upper permafrost layer in LPJmL.FAPAR seasonality in temperate grasslands in western North America and central Asia was controlled from a mixture of cold temperature (January to April), of water availability (May to November) and light (November to January).FAPAR seasonality in temperate forests in Europe was mainly limited by cold temperature in spring and by a combination of cold temperature and light in autumn.Additionally, heat stress and water availability contributed to a small reduction in summer FAPAR in temperate and boreal forests.The FAPAR seasonality in savannas (Sahel) was limited by water availability in the entire year and additionally by heat stress before the beginning of the rain season.The FAPAR seasonality of temperate regions in South America was limited by water availability in the entire year.Cold temperature was additionally limiting between May and September.Thus, water availability was the only environmental factor in LPJmL-GSI that controlled phenology globally from tropical to Arctic biomes.The implementation of the water limiting function on phenology in LPJmL-GSI resulted in unique patterns of phenological controls that were different from results reported in similar analyses (Caldararu et al., 2014;Jolly et al., 2005).LPJmL-GSI showed water limitation on phenology in many subtropical and dry temperate regions (especially Mediterranean, Pampas and Patagonia in South America, Mongolia, and northern Great Plains).The original GSI model showed mainly temperature and light limitation in these regions.In contrast to the original GSI, our implementation considers water limitations on phenology based on plant available water and not on VPD (Jolly et al., 2005).As considered by Caldararu et al. (2014), soil water availability exerts a more direct control on phenology development, which has been demonstrated for Mediterranean ecosystems (Kramer et al., 2000;Richardson et al., 2013) and in dry temperate grasslands (Liu et al., 2013;Yuan et al., 2007).
Additionally, we identify water availability as an important limiting function for spring phenology in boreal and Arctic regions in LPJmL-GSI because of the seasonal freezing of the upper active layer in permafrost soils.Although no relationships between active layer depth and vegetation greenness were found so far (Mcmichael et al., 1997), frozen grounds limit the seasonal tree growth in boreal forests because of limited water supply and nutrient uptake (Benninghoff, 1952;Jarvis and Linder, 2000).As the seasonal freezing and thawing of permafrost soils is to a large extent driven by changes in air temperature, one might argue that air temperature is enough to explain phenology dynamics in boreal and Arctic regions.Nevertheless, we found weak correlations between posterior model parameters for the cold temperature and water limiting function for phenology in PFTs that experience strong permafrost dynamics (BoNS r = 0.2, PoH r = −0.28)(Fig. S16).This indicates that the water and cold temperature limiting functions in boreal and Arctic regions are only weakly correlated.Indeed, we did not find a completely synchronized temporal dynamic of the cold temperature and water limiting functions for phenology (Fig. 9).These results emphasize the ability to disentangle effects of seasonal air temperature and soil moisture on phenology in boreal and Arctic regions.Air temperature and soil thawing are not completely synchronized because soil temperature depends also on topography, substrate and the insulating effects of the snow, litter and vegetation cover (Jorgenson et al., 2010;Shur and Jorgenson, 2007;Zhang, 2005).Soils might be still frozen if air temperature is already positive or vice versa.Also experimental studies highlighted the role of permafrost-regulated soil moisture on phenology and productivity in boreal and Arctic ecosystems (Natali et al., 2012; Schuur et al., 2007).It also has been observed that the seasonal freezing and thawing in permafrost regions regulates ecosystem evapotranspiration (Ohta et al., 2008) and fire activity (Forkel et al., 2012) especially during extreme dry years.Thus, although temperature might be enough to explain average spatial patterns of phenology in boreal and Arctic regions we acknowledge that variations in snow or vegetation cover that affects soil temperature and thus moisture might be important factors in explaining inter-annual variations of land surface phenology.
The heat stress limiting function was additionally introduced in LPJmL-GSI.Heat stress had no importance for seasonal FAPAR dynamics in most regions except in temperate and tropical grasslands.The heat stress function was highly correlated with the water availability function in temperate grasslands.This suggests that summer FAPAR is both regulated by water-induced and temperature-induced drought conditions in temperate grasslands.In tropical grasslands, heat stress and water availability were driving the temporal dynamics of seasonal FAPAR but asynchronously (in www.biogeosciences.net/11/7025/2014/Biogeosciences, 11, 7025-7050, 2014 the Sahel).These results suggest that soil moisture needs to be considered in observational data analyses and in other ecosystem models as a controlling factor for vegetation phenology in all biomes.Interestingly, Caldararu et al. (2014) identify leaf age as the dominant factor for phenology development in many permanent moist subtropical and tropical forests, but also in several water-limited regions which were here identified as seasonally controlled by water availability.We cannot identify a dominant control on seasonal FAPAR dynamics in evergreen forests, as leaf age is not explicitly simulated in LPJmL-GSI.We acknowledge that the consideration of leaf age effects on phenology could further enhance the representation of ecosystem processes.However, the seasonal co-variation between LAI or FAPAR and environmental controls on phenology complicates the ability to disentangle the leaf aging signal from a temperature, light or water availability-driven signal, especially in seasonally deciduous vegetation types, where climate-driven models explain a significant fraction of seasonal variability and the realized age of leafs is shorter than a year.In addition, cloud cover contamination over moist tropical or subtropical forests pertain usually a weak seasonal signal and a high short-term variability, hinging on the reliability of the seasonal signal.In particular, Morton et al. (2014) show that seasonal changes in MODIS LAI in the Amazon forests are linked to insufficient corrections of the sun-sensor geometry, which challenge the representation of vegetation phenology.However, in these tropical moist regions, where we find no environmental seasonal controls, and the realized age of oldest leafs are higher than a year, leaf age may be an important contributor for further consideration regarding the above-seasonal frequency of phenology.Hence, grasping the relevance of leaf longevity, especially in tropical perennial systems, would necessarily require ground observations of leaf development and litter fall to constrain leaf age parameters, as well as measurements of soil water content to address the appropriateness of soil moisture effects.

Conclusions
We have demonstrated a major improvement of the LPJmL dynamic global vegetation model by implementing a new set of phenological controls on vegetation greenness and by integrating multiple decadal satellite observations.We have proven that the original phenology model in LPJmL is unable to explain temporal dynamics of FAPAR.As an alternative we implemented a new phenology model (LPJmL-GSI) which considers effects of cold temperature, heat stress, light, and water availability on vegetation phenology.We developed a model-data integration approach for LPJmL (LPJmL-MDI) to (1) constrain model parameters against observations, (2) to directly integrate observed land cover fractions and burnt area time series and (3) to evaluate LPJmL against independent data streams.Specifically, phenology, produc-tivity, and albedo-related model parameters of LPJmL-GSI were optimized jointly against 30-year time series of satellite observations of FAPAR, against 10-year time series of vegetation albedo and against mean annual patterns of gross primary production using a genetic optimization algorithm.
The new phenology model and the parameter optimization clearly improved LPJmL model simulations.LPJmL-GSI better reproduces observed spatial patterns of gross primary production, tree cover, biomass and FAPAR than the original model.LPJmL-GSI simulates global total carbon stocks and fluxes that are closer to independent estimates than from the original model.LPJmL-GSI better represents observed seasonal, monthly, inter-annual and decadal FAPAR dynamics than the original model.The improvements of LPJmL in representing observed patterns and temporal dynamics of vegetation greenness allows assessing environmental controls on vegetation phenology and greenness.Contrasting to previous studies (Jolly et al., 2005;Stöckli et al., 2011), our results indicate that soil water availability is a major control of seasonal FAPAR dynamics not only in water-limited biomes but also in boreal forests and the Arctic tundra where water availability is regulated through seasonal thawing and freezing of the active permafrost layer.Until now the phenology of these ecosystems was mostly considered as temperature-limited.The consideration of the effect of soil water availability on phenology in LPJmL improved model simulations of greening trends in the Sahel and of browning trends in parts of the boreal forests of North America.Our results demonstrate that improved phenology models that consider seasonal effects of water availability on a continuous canopy development are needed in order to correctly explain seasonal to long-term dynamics in vegetation greenness.
The Supplement related to this article is available online at doi:10.5194/bg-11-7025-2014-supplement.Edited by: J. Kesselmeier

Figure 2 .
Figure 2. Examples of the cold temperature, heat stress, light and water limiting functions for phenology in LPJmL-GSI.Depending on the chosen parameters the functions have different shapes for each PFT.

Figure 3 .
Figure 3. Map of the dominant PFT in each grid cell as derived from SYNMAP, Köppen-Geiger climate zones and MODIS VCF.Grid cells that were used in any of the optimization experiments are shown as black crosses.Some grid cells were used in multiple optimization experiments.Grid cells that are dominated by agriculture were not used for optimization (TrML, tropical managed lands and TeML, temperate managed lands).

Figure 4 .
Figure 4. Distribution of the correlation coefficient between monthly LPJmL and GIMMS3g FAPAR (1982-2011) for several grid cells in prior model runs and optimization experiments grouped by plant functional types and biomes.Correlation coefficient for LPJmL-OP with default parameters (a, OP.prior), after grid cell level optimizations (b, OP.gc); cost for LPJmL-GSI with prior parameters (c, GSI.prior), after grid cell level optimizations (d, GSI.gc) and after PFT level optimizations (e, GSI.pft).Biomes are Tr (tropical), Te (temperate) and Bo (boreal/polar).Each distribution is plotted according to usual boxplot statistics.The point symbols indicate the plant functional type.The significance flag on top of each distribution shows if a distribution is significant different (p ≤ 0.01) to the corresponding distribution of the same PFT in another optimization experiment.The significance is based on the Wilcoxon rank-sum test.For example "acd" indicates a significant difference to the main categories (a) (OP.prior),(c) (GSI.prior) and (d) (GSI.gc)but no significant difference to (b) (OP.gc) and (e) (GSI.pft).

Figure 5 .
Figure 5. Uncertainty and sensitivity of LPJmL-GSI parameters derived from all individuals of genetic optimizations at PFT level.Shown is the relationship between parameter values and the likelihood of the corresponding parameter vector.The likelihood is normalized with the likelihood of the optimum parameter set.Only individuals with dAIC < 2 are shown.Grey areas indicate the uniform prior parameter range.Red crosses indicate the optimum parameter value.The optimum parameter value is indicated as text in a plot if it is outside of the plotting range.Results from two independent optimization experiments are shown for the BoNS, TrH, TeH and PoH PFTs (black and blue colours, respectively) but not all parameters were included in both experiments.The parameter ALBEDO_LITTER in the TrBE and TeBE PFTs was not considered in optimization experiments.

Figure 6 .
Figure 6.Best LPJmL model runs for (a) monthly FAPAR dynamics (1982-2011, n = 360 months) and (b) time series of mean annual FAPAR (1982-2011, n = 30 years).The best LPJmL model run has the highest correlation coefficient between simulated LPJmL FAPAR and GIMMS3g FAPAR.If one model run is shown the correlation coefficient of this best model is significantly higher than that of the second best model run (p ≤ 0.05, Fisher z transformation on difference in correlation).If two model runs are shown the correlation coefficients of the first and second best model runs are not significantly different from each other (p > 0.05).

Figure 7 .
Figure 7.Comparison of trends in mean annual FAPAR from LPJmL and from satellite data sets.Trends were computed between 1982 and 2011 as linear trends.The significance of a trend was determined using the Mann-Kendall trend test.Only significant trends slopes (p ≤ 0.05) are displayed in each map.Spatial correlations of trend slopes (Spearman coefficient) between LPJmL and the GIMMS3g data set are given in the map titles.Time series are showing mean annual FAPAR time series and trends spatially averaged for the regions as indicated in the first map.The blue area in time series represents the uncertainty of the GL2 VGT FAPAR data set.Numbers in the time series plot are correlation coefficients between mean annual FAPAR time series from GIMMS3g and from GL2 or LPJmL model runs, respectively.The significance of a trend and of the correlation is indicated as point symbol: * * * p ≤ 0.001, * * p ≤ 0.01, * p ≤ 0.05.p ≤ 0.1.

Figure 8 .
Figure 8. Latitudinal gradients of (a) gross primary production (GPP), (b) evapotranspiration, (c) biomass and (d) tree cover from data-oriented estimates and from LPJmL model simulations.Gradients were spatially averaged (median) from all 0.5 • grid cells for latitudinal bands of 1 • width.Grey areas represent uncertainty estimates for the data-oriented estimates.

Figure 9 .
Figure 9. Phenological controls on seasonal FAPAR dynamics.The maps are red-green-blue composites of the mean monthly values for the water (red), light (green) and cold temperature (blue) phenology limiting function values from the LPJmL-GSI model run.White regions in the maps are without vegetation or dominated by croplands for which the LPJmL-GSI phenology model was not applied.Time series represent the mean seasonal cycles (January to December) (averaged over 1982-2011) of simulated and observed FAPAR and phenology limiting function values averaged for different regions as indicated in the first map.Phenology limiting function values close to 0 indicate a strong control by phenology limiting functions whereas values close to 1 indicate no phenological control.The correlation coefficients of each time series with the simulated FAPAR time series are shown in each time series plot.The significance of the correlation is indicated as point symbol (see Fig. 7 for an explanation of significance symbols).

-
Sassan Saatchi for the tropical forest biomass map -Markus Kottek, Franz Rubel and the University of Veterinary Medicine Vienna for providing the Koeppen-Geiger climate map -Louis Giglio and Guido Van der Werf for providing the GFED database -Jennifer L. Northway and Gary Schmunk for compiling the Alaska Large Fire Database -Natural Resources Canada for providing the Canada National Fire Database -Phil Jones, Ian Harris and the CRU research unit for providing the CRU climate data set -ECMWF for ERA-Interim climate reananlysis data.M. Forkel received funding from the Max Planck Institute for Biogeochemistry and from the European Commission's 7th Framework Programme project CARBONES (grant agreement 242316).M. Forkel conducted this work under the International Max Planck Research School for Global Biogeochemical Cycles.The service charges for this open access publication have been covered by the Max Planck Society.

Table 1 .
Overview of optimization experiments with information sources for prior and posterior parameter sets.Parameter values and prior parameter ranges for each parameter set are listed in the Supplement 4.1.

Table 2 .
Overview of global model runs that were used in this study for model evaluation.