Environmental Controls on the Increasing Gpp of Terrestrial Vegetation across Northern Eurasia

Terrestrial ecosystems of northern Eurasia are demonstrating an increasing gross primary productivity (GPP), yet few studies have provided definitive attribution for the changes. While prior studies point to increasing temperatures as the principle environmental control, influences from moisture and other factors are less clear. We assess how changes in temperature, precipitation, cloudiness, and forest fires individually contribute to changes in GPP derived from satellite data across northern Eurasia using a light-use-efficiency-based model, for the period 1982–2010. We find that annual satellite-derived GPP is most sensitive to the temperature , precipitation and cloudiness of summer, which is the peak of the growing season and also the period of the year when the GPP trend is maximum. Considering the regional median, the summer temperature explains as much as 37.7 % of the variation in annual GPP, while precipitation and cloudiness explain 20.7 and 19.3 %. Warming over the period analysed, even without a sustained increase in precipitation, led to a significant positive impact on GPP for 61.7 % of the region. However, a significant negative impact on GPP was also found, for 2.4 % of the region, primarily the dryer grasslands in the southwest of the study area. For this region, precipitation positively correlates with GPP, as does cloudiness. This shows that the southwestern part of northern Eurasia is relatively more vulnerable to drought than other areas. While our results further advance the notion that air temperature is the dominant environmental control for recent GPP increases across northern Eurasia, the role of precipitation and cloudi-ness can not be ignored.


Introduction
Several analyses of normalized difference vegetation index (NDVI) data derived from satellite remote sensing have pointed to a positive trend in gross primary productivity (GPP) and leaf area index (LAI) of the northern high latitudes in the recent decades (Myneni et al., 1997;Carlson and Ripley, 1997;Zhou et al., 2001;Guay et al., 2014).Warming has also occurred over this time.Global mean surface air temperatures increased by 0.2 to 0.3 • C over the past 40 years, with warming greatest across northern land areas around 40-70 • N (Nicholls et al., 1996;Overpeck et al., 1997).Precipitation increases have also been observed over both North America and Eurasia over the past century (Nicholls et al., 1996;Groisman et al., 1991).Urban et al. (2014) describe the co-occurrence of these climatic and ecosystem changes.Here we investigate increasing GPP of terrestrial ecosystems of northern Eurasia and determine the relative attribution arising through changes in several geophysical quantities, hereinafter referred to as "environmental variables", as they potentially drive observed temporal changes in vegetation productivity.
GPP is a physical measure of the rate of photosynthesis, or the rate at which atmospheric CO 2 is fixed by autotrophic (generally green) plants to form carbohydrate molecules.Photosynthesis, being a biological process, is regulated by several environmental factors.Productivity is highest at the optimum temperature, though this optimum can be modified by cold or warm acclimation (Larcher, 1969(Larcher, , 2003)).Water availability also affects plant hydraulics and chemistry by overlaid with the spatial distribution of the 10 flux tower sites whose GPP (gross primary productivity) data were used to validate the GPP data derived from satellite NDVI (normalized difference vegetation index).For our statistical analysis, we show the distribution of two fundamental types of vegetation types: (i) herbaceous, i.e. without woody stems, which includes tundra in the north and grasslands (Eurasian Steppe) to the south, and (ii) wooded, i.e. plants with wood as its structural tissue, which includes the boreal forests appearing in the middle and extending from the western to the eastern boundary.This land cover map has been derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) Type 5 land cover product (Friedl et al., 2002).The details of the flux tower sites are listed in Table 2.
controlling the nutrient uptake through shoot transportation (Sharp et al., 2004;Stevens et al., 2004).Increasing atmospheric CO 2 concentration increases GPP by biochemical fertilization for C 3 plants and increasing water use efficiency for both C 3 and C 4 plants (Bowes, 1996;Rötter and Geijn, 1999).
There is both direct and indirect evidence of increasing productivity across the northern high latitudes.Flask-and aircraft-based measurements show that the seasonal amplitude of atmospheric CO 2 concentration across the Northern Hemisphere has increased since the 1950s, with the greatest increases occurring across the higher latitudes (Graven et al., 2013).This trend suggests a considerable role of northern boreal forests, consistent with the notion that warmer temperatures have promoted enhanced plant productivity during summer and respiration during winter (Graven et al., 2013;Kim et al., 2014;Myneni et al., 1997).Observed at eddy covariance sites, net ecosystem exchange (NEE), the inverse of net ecosystem productivity (NEP), is a strong function of mean annual temperature at mid-and high latitudes, up to the optimum temperature of approximately 16 • C, above which moisture availability overrides the temperature influence (Yi et al., 2010).Other studies have found vulnerabilities in ecosystems of North America as well as Eurasia from warming-related changes in hydrological patterns (Parida and Buermann, 2014;Buermann et al., 2014), thereby highlighting the importance of precipitation.With warming, lowtemperature constraints to productivity have relaxed (Nemani et al., 2003;Zhang et al., 2008;Yi et al., 2013).Tree-ring data suggest that black spruce forests have experienced drought stress during extreme warmth (Walker et al., 2015).Over northern Eurasia, precipitation trends have complicated the relationship between temperature and productivity, as the increasing moisture constraints have made northern Eurasia more drought-sensitive (Zhang et al., 2008;Yi et al., 2013).Increasing atmospheric CO 2 concentration is another factor, as CO 2 fertilization has been demonstrated through observations, models, and FACE (free-air CO 2 enrichment) experiments (Ainsworth and Long, 2005;Hickler et al., 2008;Graven et al., 2013).Cloudiness or shade can strongly influence vegetation productivity (Roderick et al., 2001), particularly over northern Eurasia (Nemani et al., 2003).Disturbances through forest fires also affect vegetation productivity by destroying existing vegetation and allowing for regeneration (Goetz et al., 2005;Amiro et al., 2000;Reich et al., 2001).
The role of temperature and precipitation in the positive trend of GPP of northern high latitudes, especially northern Eurasia, has not been firmly established.Few studies have examined the effect of CO 2 concentration, cloudiness, and forest fires.Of these environmental variables, CO 2 concentration is unlike the others, given its long atmospheric lifetime (∼ 100-300 years; Blasing, 2009).Thus, CO 2 concentration is assumed to be more spatially uniform.As a result, any statistical analysis using this variable will not be comparable with the other variables.We consequently do not analyse the influence of CO 2 concentration.While some studies have focused on terrestrial ecosystems of the pan-Arctic (Urban et al., 2014;Myneni et al., 1997;Guay et al., 2014;Kim et al., 2014) or the high latitudes of North America (Goetz et al., 2005;Buermann et al., 2013;Thompson et al., 2006), few studies have investigated the relative role of different environmental variables on increasing GPP of northern Eurasia.Therefore, we assess in this study how vegetation productivity trends in northern Eurasia are influenced by the environmental variables air temperature, precipitation, cloudiness, and forest fire.Objectives are to (1) calculate the long-term trend of both GPP and the environmental variables, (2) assess the magnitude of the effect of the environmental variables on GPP, (3) identify the seasonality of the variables, and (4) identify the regions of northern Eurasia where the variables boost or reduce GPP.Exploiting the availability of long-term time series observation-based data we perform a spatially explicit grid point statistical analysis to achieve the above objectives.

Land cover
The study domain is the Northern Eurasia Earth Science Partnership Initiative (NEESPI) region (Groisman and Bartalev, Thus the effect of each environmental variable accounts only for changes in NDVI and does not track potential changes in land cover type.While the GPP products use the standard IGBP MODIS global land cover classification, for our statistical analysis we simplify the LC distribution into two fundamental types.One is "herbaceous", without woody stems, found in the tundra to the north and grasslands to the south, one of the driest biomes of northern Eurasia.The second is "woody vegetation", plants with woody stems, located within the area of boreal forests extending from west to east across much of the centre of the domain (Fig. 1).

Vegetation productivity -long-term data
GPP represents the total amount of carbon fixed per unit area by plants in an ecosystem utilizing the physiological process of photosynthesis (Watson et al., 2000).GPP is one of the key metrics useful in assessments of changes in vegetation productivity.It is also a standard output of processbased vegetation models.The GPP fields used in this study represent model estimates driven by satellite data.The GPP model used is based on a light use efficiency (LUE) model that prescribes theoretical maximum photosynthetic conversion efficiency for different land cover classes.LUE is reduced from potential (LUE max ) rates for suboptimal environmental conditions determined as the product of daily environmental control factors defined for the different land cover types using daily surface meteorological inputs from ERA-Interim reanalysis data.Daily surface meteorology inputs to the model include incident solar radiation (SW rad ), minimum and average daily air temperatures (T min and T avg ), and atmospheric vapour pressure deficit (VPD).GPP is derived on a daily basis as (Running et al., 2004;Zhang et al., 2008) where ε is a LUE parameter (g C MJ −1 ) for the conversion of photosynthetically active radiation (PAR, MJ m −2 ) to GPP.FPAR is estimated from NDVI using biome-specific empirical relationships emphasizing northern ecosystems (Yi et al., 2013).Several studies demonstrated the linear relationship between NDVI and FPAR through field measurements and theoretical analysis (Fensholt et al., 2004;Myneni and Williams, 1994;Ruimy et al., 1994;Sellers, 1985).Two sets of NDVI records are obtained for this study and used to derive alternative FPAR and GPP simulations: (i) the third generation Global Inventory Modeling and Mapping Studies (GIMMS3g; Zhu et al., 2013;Pinzon and Tucker, 2010), downloaded from https://nex.nasa.gov/nex/(referred to as GIMMS-GPP), and (ii) the Vegetation Index and Phenology (VIP) database (Didan, 2010;Barreto-Munoz, 2013), downloaded from http://phenology.arizona.edu/(University of Arizona's Vegetation Index and Phenology Lab; referred to as VIP-GPP).The 16-day NDVI records are first interpolated to a daily time step using temporal linear interpolation to estimate daily FPAR following previously established methods (Yi et al., 2013).The use of daily NDVI and FPAR inputs rather than coarser (8-day or 16-day) temporal composites reduces potentially abrupt step changes in the model calculations due to temporal shifts in the coarser time series canopy inputs.Moreover, the daily interpolation was found to improve simulations of GPP seasonality es-

P. Dass et al.: Vegetation greening across northern Eurasia
pecially during spring and autumn transitional periods over northern land areas (Yi et al., 2013).PAR is estimated as a constant proportion (0.45) of incident shortwave solar radiation (SW rad ).ε max is the potential maximum ε under optimal environmental conditions.T f and VPD f are scalars that define suboptimal temperature and moisture conditions represented by respective daily T min and VPD inputs.T f and VPD f are defined using a linear ramp function (Yi et al., 2013;Heinsch et al., 2006), as well as minimum and maximum environmental constraints defined for different biome types (T mn_min and T mn_max , VPD min and VPD max ).Table 1 summarizes the biome property look-up table (BPLUT) used to define the environmental response characteristics in the model.These GPP data sets are currently available through a public FTP directory (ftp://ftp.ntsg.umt.edu/pub/data/HNL_monthly_GPP_NPP/).The GPP data are derived at a daily time step and have been aggregated to a monthly time step for this study.Spatial resolution is 25 km, with a temporal range from 1982 to 2010, restricted to the northern high latitudes (> 45 • N).In many of the statistical analyses to follow we use the ensemble mean of the two satellite-derived GPP data sets, henceforth denoted as "GPPsat".Winter is characterized by extremely low productivity, and technical constraints of optical-IR remote sensing due to low solar illumination and persistent cloud cover make for a particular challenge in estimating vegetation indices and consequently computing GPP across the high latitudes (Pettorelli et al., 2005).Given the limited confidence in GPP data over winter (driven mainly by the uncertainty in winter NDVI) we focus on the remainder of the year in our analysis.
Accuracy of the GIMMS-NDVI data set has been examined in several recent studies.Analysing trends in growingseason start over the Tibetan Plateau, Zhang et al. (2013) found that GIMMS NDVI differed substantially over the period 2001-2006 from SPOT-VGT and MODIS-NDVIs, indicating significant uncertainty among NDVI retrievals from different satellite sensors and data records.The GIMMS3g data set is based on the NOAA-AVHRR (Advanced Very High Resolution Radiometer) long-term time series record, which is comprised of AVHRR2 and AVHRR3 sensors on board the NOAA-7 through to NOAA-19 satellites spanning multiple overlapping time periods; this leads to potential artifacts from cross-sensor differences and inter-calibration effects influencing long-term trends in the AVHRR NDVI time series (Pinzon and Tucker, 2014).The Vegetation Index and Phenology (VIP) NDVI data set applies a different data processing scheme from that of GIMMS3g (Fensholt et al., 2015), and involves an integration and calibration of overlapping AVHRR, SPOT, and MODIS sensor records for generating consistent NDVI (Didan, 2010).The ensemble mean and variance of alternate GPP calculations derived using the GIMMS3g and VIP NDVI records was used as a metric of uncertainty in the regional productivity trends and underlying satellite observation records.

Flux tower data
To verify the satellite-based GPP estimates we use gap-filled daily tower GPP data at 10 flux tower sites distributed across northern Eurasia, available for different periods of time.Details of the individual towers are provided in Table 2.The data, generated using the eddy covariance measurements acquired by the FLUXNET community, were collected from http://www.fluxdata.org/for the "free fair-use" data subset.The spatial distribution of the flux towers used in this study is shown in Fig. 1.Unless otherwise noted, we use seasonal totals of the daily gap-filled tower GPP data.Monthly and seasonal values were aggregated from the daily data.
We also use monthly GPP data computed using FLUXNET observations of carbon dioxide, water and energy fluxes upscaled to the global scale for additional verification of the satellite-derived GPP record for the entire study area, on a per grid cell basis.Upscaling of the FLUXNET observations was performed using a machine learning technique and model tree ensembles (MTE) approach from the Max Planck Institute of Biogeochemistry, Jena, Germany, and available online at https://www.bgc-jena.mpg.de/geodb/projects/Data. php.Description and benchmarking of this data set can be found in Jung et al. (2009) and Jung et al. (2011).Of the two versions available, we use the one which incorporates flux partitioning based on Reichstein et al. (2005).

Temperature, precipitation, and cloudiness
Monthly values of 2 m air temperature (in • C), precipitation (in mm), and cloudiness (in %) are taken from monthly observations from meteorological stations, extending over the global land surface and interpolated onto a 0.5 • grid (Mitchell and Jones, 2005).The data set, CRU TS 3.21, is produced by the Climatic Research Unit of the University of East Anglia in conjunction with the Hadley Centre (at the UK Met Office) and is available at http://iridl.ldeo.columbia.edu/SOURCES/.UEA/.CRU/.TS3p21/.monthly/(Jones and Harris, 2013).
Although the LUE-based GPP model does not use precipitation as an input, we assume that precipitation is a useful metric of water supply to vegetation and thus analyse it as one of the environmental variables affecting GPP.Here we use monthly values of temperature, precipitation, and cloudiness for the period of 1982 to 2010, since this is the common period for which both GPPsat and the environmental variable data are available.Seasonal means for spring (March, April, May), summer (June, July, August), and autumn (September, October, November) are derived from the monthly values.As explained in Sect.2.1.2,lower reliability and availability of satellite NDVI observations and associated GPP data for the winter months lead us to focus on the spring, summer, and autumn seasons.

Fire
Fire is represented by proportional burnt area (% of each grid cell) estimates from the Global Fire Emissions Database (GFED) Monthly Burned Area Data Set Version 3.1 released in April 2010.This product was developed on a global scale at a 0.5 • spatial resolution and covers the period from 1997 to 2011.The GFED is an ensemble product of burn areas derived from multiple satellite sensors, though primarily emphasizing MODIS surface reflectance imagery (Giglio et al., 2010).

Spatial interpolation
Data not on a 0.5 • grid were interpolated to that resolution using spherical version of Shepard's traditional algorithm (Shepard, 1968;Willmott et al., 1985).This method takes into account (i) distances of the data points to the grid location, (ii) the directional distribution of stations in order to avoid overweighting of clustered stations, and (iii) spatial gradients within the data field in the grid point environment.

Verification
The GIMMS-GPP and VIP-GPP simulations are evaluated against co-located tower-based GPP observations for model grid cells corresponding to each of the ten regional flux tower locations (Table 2).The evaluation is carried out using five different approaches: 1. Pearson's product moment correlation, which is a measure of the linear dependence between simulated (GIMMS-GPP and VIP-GPP) and observed (towerbased GPP) values and its value ranges from −1 to +1, where 0 is no correlation and −1/ + 1 is total negative or positive correlation respectively.
2. Percent bias, which measures the average tendency of the simulated values to be larger or smaller than the corresponding observations.The optimal value is 0.0 with low-magnitude values indicating accurate model simulations.Positive values indicate overestimations and vice versa (Yapo et al., 1996;Sorooshian et al., 1993).
3. The Nash-Sutcliffe efficiency (NSE) coefficient, which is a normalized statistic that determines the relative magnitude of the residual variance compared to the measured data variance (Nash and Sutcliffe, 1970).The statistic indicates how well the plot of observed vs. simulated data fits the 1 : 1 line.Nash-Sutcliffe efficiencies range from −∞ to 1.An efficiency of 1 corresponds to a perfect match of model-simulated GPP to the observed www.biogeosciences.net/13/45/2016/Biogeosciences, 13, 45-62, 2016 Table 3. Validation of GIMMS3g and VIP-GPP data sets along with their ensemble mean using flux tower GPP from 10 flux tower sites across northern Eurasia.The spatial distribution of the flux tower sites is shown in Fig. 1.Validation was carried out using the following approaches.
(1) Pearson's product moment correlation, which is a measure of the linear dependence between the simulated and observed GPP and its value ranges from −1 to +1, where 0 is no correlation and −1/ + 1 is total negative or positive correlation.
(2) Percent bias, which measures the average tendency of the simulated values to be larger or smaller than their observed ones.The optimal value is 0.0, with low-magnitude values indicating accurate model simulations.Positive values indicate overestimations and vice versa (Yapo et al., 1996;Sorooshian et al., 1993).( 3) Nash-Sutcliffe model efficiency coefficient (Nash and Sutcliffe, 1970), values of which range from −∞ to 1.An efficiency of 1 corresponds to a perfect match of model-simulated GPP to the observed data.An efficiency of 0 indicates that the model predictions are as accurate as the mean of the observed data, whereas an efficiency less than zero occurs when the observed mean is a better predictor than the model or, in other words, when the residual variance (between modelled and observed values) is larger than the data variance (between observed values and the observed mean).Essentially, the closer the model efficiency is to 1, the more accurate the model is.data.An efficiency of 0 indicates that the model predictions are as accurate as the mean of the observed data, whereas an efficiency less than zero occurs when the observed mean is a better predictor than the model or, in other words, when the residual variance (between modelled and observed values) is larger than the data variance (between observed values and the observed mean).
Essentially, the closer the model efficiency is to 1, the more accurate the model is.
4. A scatter plot, which demonstrates using Cartesian coordinates the correlation between satellite-derived GPP and tower-derived GPP at the respective sites for the respective time periods.This along with the line of best fit helps determine how well the two data sets agree with each other.
5. Spatially explicit, pixel-by-pixel validation using the upscaled GPP data from FLUXNET observations (described in Sect.2.1.3)using correlation and difference maps for the entire period.

Trend analysis
Temporal changes for each environmental variable are determined using linear regression.Both annual and seasonal time integrations are examined.Trends are deemed statistically significant at the 95 % level.For each variable, we compute the trend per decade (10 yr −1 ) from the monthly values (month −1 ).Other studies have implemented a similar methodology to identify trends (Piao et al., 2011;de Jong et al., 2011;Forkel et al., 2013;Goetz et al., 2005).In order to determine whether the temporal rate of change differs for different periods of the study period we plot the percentage difference of the annual means (of the regional average) from that of the first 5-year mean.
For the entire period of study, a few of the variables assessed show strong trends.Moreover, we assume the variables to be linearly associated.This introduces the issue of collinearity, as a consequence of which the study of influence of one variable on another becomes less precise.Therefore, in order to make accurate assessments of correlation between two variables, correlation analysis has only been carried out after long-term trends (for period 1982-2010) have been removed and consequently only the interannual variability is preserved.

Correlation
We use the Pearson product-moment correlation coefficient (represented as R), one of the more popular measures of dependence between two variables and which is sensitive only to a linear relationship between two variables.This metric is defined from +1 (perfect increasing linear relationship) to −1 (perfect decreasing linear relationship or "inverse correlation") and as the value approaches zero, the relationship becomes uncorrelated (Dowdy and Wearden, 1983).When a single variable is affected by more than one independent factor, simple correlation is inappropriate.We perform partial correlation to better assess the relationship between two variables after eliminating the influence of other variables.

Attribution
The primary objective of this study is to determine the magnitude and spatio-temporal variations in trends for environmental conditions (variables) which have contributed to the increase in GPP of northern Eurasia indicated from the satellite records.Ideally one would study the direct influence of one condition on another in experiments in which all other possible causes of variation are eliminated.However, since this study involves only large-scale observational data and not process-based models or laboratory-based experiments, there is no control over the causes of variation.Investigations into the structure and function of terrestrial ecosystems, like those for many elements of the biological sciences, involve quantities which are often correlated.In some cases, the derived relationship may be spurious.The coefficient of determination (represented as R 2 ) is a common measure to estimate the degree to which one variable can be explained by another (percentage; Wright, 1921), while correlation anal-ysis (R) can explain this dependence of one variable on another keeping the sign of the relationship (±) intact (Aldrich, 1995).

Verification of satellite-derived GPP
The GIMMS-GPP and VIP-GPP, as well as their ensemble mean (GPPsat), are individually verified against the fluxtower-based GPP data using Pearson's correlation coefficient, percent bias, and the Nash-Sutcliffe normalized statistic.Scatter plots (Fig. 2) show that GPP derived from the satellite NDVI records is generally higher than the towerbased GPP at the flux tower sites that have comparatively lower productivity (and vice versa).Moreover, the agreement is stronger at lower-productivity sites than at higherproductivity sites.Though Table 3 lists all of the verification statistics, we focus primarily on the annual GPPsat results for the rest of the study.The correlation coefficients are all positive and high (0.7 for annual GPPsat); percent bias is predominantly negative (18.3 %); and since all the values of the Nash-Sutcliffe efficiencies are above zero (0.33), we conclude that the satellite NDVI-derived values are a more accurate estimate of GPP than the observed mean for the respective flux tower sites.Spatially explicit verification of GPPsat reveals that the correlation is high and statistically significant for almost the entire study area (Fig. 3a).GPPsat shows a general underestimation in the boreal forests of the western parts of northern Eurasia and overestimation in the Eurasian steppes to the south of the study area (Fig. 3b).
Satellite-derived vegetation indices have been evaluated using a variety of techniques.Using tree-ring width measurements as a proxy for productivity, Berner et al. (2011) examined its relationship with NDVI from AVHRR instruments and found the correlation to be highly variable across the sites, though consistently positive.Remarkably strong corre- shows yearly change in the regional average GPP for the data sets derived from the GIMMS3g (red) and VIP (blue) NDVI data sets.The interannual variation is smoothed using a smoothing spline using a smoothing parameter of 0.8.lations were observed in comparisons of GIMMS3g NDVI to aboveground phytomass at the peak of summer at two representative zonal sites along two trans-Arctic transects in North America and Eurasia (Raynolds et al., 2012).From comparison of production efficiency model-derived NPP (Zhang et al., 2008) to the stand level observations of boreal aspen growth for the 72 CIPHA (Climate Impacts on Productivity and Health of Aspen) sites, the correlation was found to be positive.LUE algorithms similar to the one used in this study for the generation of GPP data sets from satellite NDVI produce favourable GPP results relative to daily tower observations, with a strong positive correlation (Yi et al., 2013;Yuan et al., 2007;Schubert et al., 2010).Evaluating the uncertainties in the estimated carbon fluxes computed using a similar LUE-based GPP model, Yi et al. (2013) concluded that the uncertainty in LUE (ε) characterization is the main source of simulated GPP uncertainty.GPP simulation errors under dry conditions are increased by an insufficient model vapour pressure deficit (VPD) representation of soil water deficit constraints on canopy stomatal conductance and ε (Leuning et al., 2005;Schaefer et al., 2012).It was also found that the GPP model does not consider the response of ε to diffuse light due to canopy clumping (Chen et al., 2012) and shaded leaves (Gu et al., 2002).

Temporal changes in GPP
Across the study domain, regionally averaged GPPsat exhibits a trend of 2.2 (±1.4) g C m −2 month −1 decade −1 .Figure 4a displays the annual GPP trend map.Increases are noted across most of the region except for a small area in the north-central part of the region, just east of the Yenisey River.The largest increases are located in the western and south-eastern part of the region.Over half (69.1 %) of the study area exhibits a statistically significant positive trend (95 % significance level), while 0.01 % of the area has a statistically significant negative trend.Uncertainty in the ensemble mean GPP is illustrated by the coefficient of variation map (Fig. 4b).The highest uncertainty is noted in the north-central and the south-western part of the region.The yearly increase in annual GPP for both GIMMS-GPP (red) and VIP-GPP (blue; Fig. 4c) reveal the difference between the two data sets, which is highest at the beginning of the study period.The nature of increase in GPP is also different for the two data sets, with the rise in one being more linear than the other.A possible explanation for the differences in the two data sets is discussed in Sect.2.1.2.Examining the seasonality of GPP trends (of GPPsat; Fig. 5), we find that the summer trend is greatest among all other seasons.This implies that the response of GPP to environmental changes is greatest at the peak of the growing season.While the productivity of the region is predominantly increasing, there are clearly certain areas each season with decreasing productivity.
The GPP increase described here is consistent with the results of Sitch et al. (2007), who also noted considerable interannual and spatial variability, with many areas demonstrating decreased greenness and lower productivity.Using a processbased model (LPJ-DGVM) to perform a retrospective analysis for the period of 1982-1998, Lucht et al. (2002) ) found, after accounting for the carbon loss due to autotrophic respiration, that boreal zone NPP increased by 34.6 g C m −2 yr −1 , which is comparable to our estimate.The higher GPP trend in summer (Fig. 5), especially over the northern Eurasia portion of the domain, suggests that the vegetation of this region is predominantly cold-constrained, a finding described in other recent studies (Yi et al., 2013;Kim et al., 2014).

Temporal changes in the environmental variables
The regionally averaged air temperature increase is nearly monotonic and the distributions displayed in Fig. 6a show that the region has a predominantly positive trend for all parts of the growing season.Warming is highest in autumn.A statistically significant increase in temperature is noted for approximately half of the region.The greatest increases are found in the north-eastern and south-western parts of the region (maps not shown).Unlike temperature, precipitation does not exhibit a sustained increase over the study period.While the regional median trend for precipitation is highest for spring (Fig. 6b), the range of trends for this region, from minimum to maximum, is highest for summer.The fraction of the region experiencing significant increases in annual precipitation is about 3 times the area experiencing significant decreases.The significant positive trends are located in the  north-eastern and western parts (mainly boreal forests) of the domain, while significant negative trends are located in the west-central (boreal forests) and south-eastern (steppes) parts of the region (maps not shown).Along with the regional averages of other environmental variables, Table 4 reveals the regional average of cloudiness, which shows a negative trend.However, similar to precipitation, the spatial standard deviation is very high, implying a high spatial variability in cloudiness trends across the region.Unlike precipitation, a greater fraction of the region is experiencing significant decreasing cloudiness or a significant clear-sky trend (Fig. 6c).
Compared to the rest of the region, annual cloudiness shows higher negative trends in the southern parts of the study area (maps not shown).Burnt area exhibits significant trends, both positive and negative, over only 1 % of the region, with the total yearly burnt area for the study area increasing from 15.9 to 17.1 million hectares from 1997 to 2010.The negative trend of the regional mean (Table 4; Fig. 6d) is not significant.
Recent studies have reported similar changes in these environmental variables.For the period of 1979period of to 2005period of , Trenberth et al. (2007) ) found temperature trends over the region range from 0.3 to 0.7 • C decade −1 , and for most regions of the higher latitudes, especially from 30 to 85    (Harris et al., 2014).Burnt area data from the Global Fire Emissions Database (GFED; Giglio et al., 2010).
positive precipitation trends have occurred.Contrary to the cloud cover trend we find here, studies reported in AR4 suggest an increase in total cloud cover since the middle of the last century over many continental regions, including the former USSR and western Europe (Sun et al., 2001;Sun and Groisman, 2000).The large spatial variability in the gridded cloud cover trends (Table 4) may explain the disagreement.Burnt area data, representing fire disturbance, is dissimilar from the other environmental variables in that it spans only 14 years of the 29-year study record, and it is spatially nonuniform, involving only a fraction of the total study area.This limitation makes it difficult to assess impacts on vegetation productivity (Balshi et al., 2007).While the model used to generate the satellite NDVI-derived GPP data does not account for CO 2 fertilization directly, the fertilization effect may be partially represented through associated changes in NDVI.As stated in Sect. 1, we do not analyse atmospheric CO 2 concentration due to its spatial homogeneity.

Attributing GPP changes to environmental variables and assessing seasonality
Annual GPP is affected by more than one environmental variable.To study the impact of an individual environmental variable, we eliminate the impact of other variables by performing partial correlations.With the temporal range of the fire data (GFED) being a fraction of that of the other environmental variables, it is not possible to compute the partial correlation.Consequently, we are unable to assess the effects of only fire by eliminating the effects of the other variables.Moreover, fires have been found to be significantly correlated with annual GPP (GPPsat) for only a small fraction (1.7 to 3.4 % depending on season) of the entire study area.The impact of fires on annual GPP for the region is therefore ignored in this study.The regional median partial coefficient of determination (R 2 ) for significant values (Table 5) suggests that the summer values of the environmental variables have the highest influence on annual GPPsat.The contrast between summer and the other seasons is strongest for temperature, highlighting the importance of summer temperatures to annual productivity.Figure 7 reveals that the relationships between annual GPP and the environmental variables are not completely explained by simple correlation (R 2 ), as the distributions of partial correlations provide more information about the interaction.Considering only significant correlations (Fig. 7), we find that increasing temperatures predominantly increase GPP.The relationship between precipitation or cloudiness and GPP, on the other hand, leads to a predominantly bimodal distribution, with both positive and negative effects.Other than spring, areas demonstrating significant negative partial correlations appear to be larger than the areas of significant positive partial correlations.Among the environmental variables assessed, temperature has the highest partial coefficient of determination (Table 5).Moreover, unlike precipitation and cloudiness, temperature has a predominantly Bean plots of the multi-modal distribution for significant (95 % significance) partial correlation between annual de-trended GPP (GPPsat) and the values of each de-trended environmental variable after eliminating the influence of the other variables.A bean plot is an alternative to the box plot and is fundamentally a one-dimensional scatter plot.Here it is preferred over a box plot as it helps to show a multi-modal distribution.The thickness of a "bean" is a function of the frequency of the specific value -that is, the thicker a "bean" is for a value, the relatively higher the number of grid points having that value.The values shown are the Pearson's correlation coefficients which are based on the linear least-squares trend fit.Correlation values range from −1 to +1.Values closer to −1 or +1 indicate strong correlation, while those closer to 0 indicate weak correlation.The colour of the box indicates the season of the environmental variable being investigated (annual: grey; spring: green; summer: red; autumn: amber).The short horizontal black lines for each "bean" is the median of that distribution.
positive relationship with annual GPP.These relationships imply that, over recent decades, low temperatures have been the major constraint for GPP in northern Eurasia.
Similar results were reported by Yi et al. (2014), who concluded that satellite-derived vegetation indices show an overall benefit for summer photosynthetic activity from regional warming and only a limited impact from spring precipitation.The dominant constraint of temperature was described by Zhang et al. (2008), who found the same constraint to be decreasing.However, our results contrast with those of Piao et al. (2011), who concluded that at the continental scale of Eurasia, vegetation indices in summer are more strongly regulated by precipitation, while temperature is a relatively stronger regulator in spring and autumn.Regarding the dominance of temperature as a regulator, Yi et al. (2013) concluded that, over the last decade, Eurasia has been more drought-sensitive than other high-latitude areas.
Table 5. Medians of the distributions of the relative partial significant contribution (R 2 -95 % significance) of each de-trended environmental variable (except fire) of each season to the interannual variability in de-trended annual GPP (GPPsat).In each case the total contribution may not add up to 100 %.In these cases the factors behind the unexplained attribution are not identified.Since GPP trends are highest in summer (Fig. 5), the peak of the growing season, we are interested more in the impact of the environmental variables during summer on annual GPP since the terrestrial vegetation is likely to be more responsive to variations in summer environmental conditions relative to other seasons.Spatial analysis helps to elaborate on the results shown in Table 5 and Fig. 7. Assessing the partial significant correlation of annual GPP and summer temperature (Fig. 8a; Table 6), we find that areas with a positive correlation (62 % of the area) are concentrated to the north and east of the region, which include both tundra and boreal forest areas.Negative correlations occur across 2 % of the region, largely in the south within the Eurasian steppes.For other parts of the year (maps not shown for spring and autumn correlations but distributions represented in Fig. 7), significant negative correlations become more spatially disperse, while significant positive correlations are limited to the centre and west of the region for spring, becoming more disperse in autumn.Determining the partial correlation between annual GPP and summer precipitation, Fig. 8b reveals that the areas of significant positive correlations (4 % of area) are scattered over the southern part of the study area (steppes vegetation), while the significant negative correlations (16 % of area) are scattered across the north (tundra and boreal).Correlations for spring precipitation with annual GPP (maps not shown) are predominantly positive, while that for autumn precipitation is predominantly negative.The spatial correlations for summer cloudiness and summer precipitation are similar (Fig. 8c), though the area under significant correlation is comparatively less.Negative correlation areas are about 9 times more extensive than positive correlation areas (Table 6).Compared to summer, the area under significant positive correlation is higher for spring, while the area under negative correlation is higher for autumn (maps not shown).

Environmental variable
The negative correlations for temperature and positive correlations for precipitation and cloudiness in the southern grasslands (Eurasian steppes) are not surprising, as these grasslands are relatively dry compared to other biomes in the broader region.In this part of the study area, increasing temperatures in summer may lead to greater water stress (Gates, 1964;Wiegand and Namken, 1966;Jackson et al., 1981)  over, increasing cloud cover would tend to lead to a higher probability of rain (Richards and Arkin, 1981), thus relieving water stress induced by warming in this relatively dry area.The cause of the negative correlations in the north is unclear.The relationship may be attributable to the predominantly positive relationship between cloud cover (equivalent to inverse of sunshine duration) and precipitation (Sect.3.5).
In the light-limited and relatively colder north, an increase in cloud cover could, on the one hand, cause a decrease in direct radiation and increase in diffuse radiation, which may increase GPP through higher LUE (Alton et al., 2007;Gu et al., 2002;Williams et al., 2014;Roderick et al., 2001).However, an increase in cloud cover could decrease total solar radiation and, in turn, productivity (Nemani et al., 2003;Shim et al., 2014).
Recent studies have shown similar relationships to those found here.Zhang et al. (2008) showed that, across the pan-Arctic basin, while productivity increased with warming, increasing drought stress can offset some of the potential benefits.However, Yi et al. (2013) concluded that while GPP was significantly higher during warm years for the pan-Arctic, the same was not true for the Eurasian boreal forests, which showed greater drought sensitivity.Positive impacts of warming on GPP have been suggested in warming experiments (Natali et al., 2013).However, decreasing growing- season forest productivity, represented as a decline in "greenness" across northern Eurasia, may be a reflection of continued summer warming in the absence of sustained increases in precipitation (Buermann et al., 2014;Zhou et al., 2001).

Relationships among individual environmental variables
Environmental variables are not independent of one another.We examine correlations among the de-trended individual variables to better understand their interactions.Figure 9 shows distributions of the correlations.The temperatureprecipitation correlation is predominantly negative, indicating that increases in precipitation did not accompany recent warming.Significant negative trends are located in the southern parts of the study area (steppes) as well as the boreal forests at the western and eastern ends of the region.These changes may be leading to increasing water stress, evidence of which is noted in a subset of the region.Indeed, approximately 2.4 % of the area in the southern parts of the study area (Fig. 8a) shows significant negative partial correlation between annual GPP (GPPsat) and summer temperature.The relationship between temperature and cloud cover is similarly predominantly negative.Spatially, however, the significant negative correlations are located in the central and western parts of the region.Grid-cell-wise correlations between precipitation and cloud cover are predominately positive, with the significant correlations spread out across the region.
As described in Sect.3.4, the correlations between precipitation and cloud cover help to explain why spatial distributions of the correlation coefficients of precipitation and cloud cover with GPP are similar.Wang et al. (2014) documented a positive relationship between sunshine duration (equivalent to the inverse of cloud cover) and vegetation greenness.
While increasing cloud cover leads to an increased probability of precipitation, and thus reduces water stress, it also reduces the sunshine duration and hence GPP.According to Table 4, regional mean precipitation has a positive trend, while cloudiness has a negative trend.However, Fig. 9 reveals the predominantly positive correlation between these two variables.This apparent contradiction is because the long-term trends are calculated for the actual values, while the correla-tion analysis is performed after de-trending (removing longterm trends) the variables.Consistent with our results, Thompson et al. (2006) found that, in the boreal and tundra regions of Alaska, NPP decreased when it was warmer and dryer and increased when it was warmer and wetter.They also described how colder and wetter conditions also increased NPP.Yi et al. (2013) concluded that while, globally, annual GPP for boreal forests is significantly higher in warmer years, the relationship does not hold true for Eurasian boreal forests, which they identify to be more drought-sensitive.For this reason, regional GPP variations are more consistent with regional wetting and drying anomalies, as we note for the south-western part of the study region.In this study we assessed only GPP.Other carbon cycle processes such as autotrophic and heterotrophic respiration and disturbances may not be responding in a similar manner.Additional studies are required before extrapolating these results to other carbon cycle components.

Conclusions
The ensemble mean of the GPP data sets derived from GIMMS3g and VIP NDVI data indicates that vegetation productivity generally increased across northern Eurasia over the period 1982 to 2010, with a significant increase for as much as 69.1 % of the region.A significant decrease in GPP occurred across only 0.01 % of the region.We note some disagreement in the nature and magnitude of the increasing GPP among the two data sets.The regional mean trend for the ensemble mean GPP is 2.2 (±1.4) g C m −2 month −1 decade −1 .The regional analysis is consistent with results of prior studies which have suggested that air temperature is the dominant environmental variable influencing productivity increases across the northern high latitudes.Examining partial coefficients of determination (R 2 ), we find that the summer values of temperature, precipitation, and cloudiness have the highest influence on annual GPP.Considering the regional median of partial significant R 2 values, summer air temperature explains as much as 37.7 % of the variation in annual GPP.In contrast, precipitation and cloudiness explain 20.7 and 19.3 % respectively.A significant positive partial correlation between summer air temperature and annual GPP is noted for 61.7 % of the region.For 2.4 % of the area, specifically the dryer grasslands in the south-west, temperature and GPP are inversely correlated.Precipitation and cloudiness during summer also impart a significant influence, showing areas with both positive and negative significant partial correlation with annual GPP.Fire has a very small effect, with only up to 3.4 % of the region showing significant correlation, and consequently the impact of fire on GPP was ignored for the subsequent analysis.The spatial analysis reveals that the statistical relationships are not spatially homogeneous.While warming likely contributed to increasing productivity across much of the north of the region, the relationship reverses in the southern grasslands, which are relatively dry.That region exhibits increasing GPP, but with warming accompanying increased moisture deficits potentially restricting continued productivity increase.This result demonstrates that vegetation has been resilient to drought stress, which may be increasing over time.
We recommended that this study be followed up with experiments conducted using process-based models in which a single forcing variable independent of the others is manipulated.If feasible, multiple models should be used in order to quantify the uncertainty due to differences in model parameterization.Depending on emissions, population, and other forcing scenarios, rates of change in the environmental drivers such as air temperature and precipitation may be different than those found in this study.Thus it is critical to examine future scenarios of change across the region to better understand terrestrial vegetation dynamics under the respective model simulations.Environmental drivers influence other elements of the carbon cycle beyond the individual plant.In order to determine how terrestrial carbon stocks and fluxes have changed in the recent past, or may change in the near future, all aspects of the carbon cycle should be investigated in the context of changes in overarching climate influences.

Figure 1 .
Figure1.Simplified land cover for northern Eurasia for year 2007 overlaid with the spatial distribution of the 10 flux tower sites whose GPP (gross primary productivity) data were used to validate the GPP data derived from satellite NDVI (normalized difference vegetation index).For our statistical analysis, we show the distribution of two fundamental types of vegetation types: (i) herbaceous, i.e. without woody stems, which includes tundra in the north and grasslands (Eurasian Steppe) to the south, and (ii) wooded, i.e. plants with wood as its structural tissue, which includes the boreal forests appearing in the middle and extending from the western to the eastern boundary.This land cover map has been derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) Type 5 land cover product(Friedl et al., 2002).The details of the flux tower sites are listed in Table2.

)Figure 2 .
Figure2.Relationship between the annual GPP recorded at the flux tower sites and the corresponding values of the satellite-derived GPP.The black solid line is the line of best fit and helps better understand the relationship between the two.The dashed line is the 1 : 1 line and demonstrates how much the relationship between the two sets of values deviates from the 1 : 1 perfect relationship.

Figure 3 .
Figure 3. Spatially explicit validation of GPPsat using upscaled FLUXNET observations.Panel (a) is the correlation map and displays the statistically significant (95 % level) correlations between the two sets of values of annual GPP for the period of 1982-2010.Panel (b) is the difference between the 29-year mean of GPPsat and upscaled FLUXNET database, with negative values demonstrating an underestimation of GPPsat and vice versa.

Figure 4 .
Figure 4. Change in annual GPP for GPPsat over the period 1982-2010.Panel (a) is the trend map for GPPsat, i.e. the ensemble mean (of two GPP data sets).Shades of green represent a positive trend and shades of red represent a negative trend.The trends have been derived from a linear least-squares fit to the GPP time series for GIMMS3g and VIP data sets.Trend values represent the rate of change of productivity per decade (g C(Carbon) m −2 month −1 10 yr −1 ).Panel (b) is the uncertainty map (uncertainty due to the use of two GPP data sets) represented by computing the coefficient of variation (CV).Darker values represent higher uncertainty and vice versa.Panel (c)shows yearly change in the regional average GPP for the data sets derived from the GIMMS3g (red) and VIP (blue) NDVI data sets.The interannual variation is smoothed using a smoothing spline using a smoothing parameter of 0.8.

Figure 5 .
Figure5.Box plot showing grid distributions of seasonal GPP trends for GPPsat.The GPP trends are in g C m −2 month −1 10 yr −1 .The black band and middle notch represent the 2nd quartile or median; box extents mark the 25th (1st quartile) and 75th (3rd quartile) percentiles.Whiskers extend from the smallest non-outlier value to the largest non-outlier value.The colours, green, red, orange, and grey represent spring, summer, autumn, and annual seasonal trends respectively.As described in Sect.2.1.2,GPP trends for winter have not been assessed in this study.

)Figure 6 .
Figure 6.Change in the environmental variables over the period of study represented by seasonal trends.Panels (a-c) show distribution of 2 m air temperature, precipitation, and cloud cover respectively for the period 1982-2010, and panel (d) illustrates seasonal trends of total burnt area for the period 1997-2011.The temperature, precipitation, and cloud cover data are taken from the Climatic Research Unit (CRU TS 3.21) data set(Harris et al., 2014).Burnt area data from the Global Fire Emissions Database (GFED;Giglio et al., 2010).
Figure 7. Bean plots of the multi-modal distribution for significant (95 % significance) partial correlation between annual de-trended GPP (GPPsat) and the values of each de-trended environmental variable after eliminating the influence of the other variables.A bean plot is an alternative to the box plot and is fundamentally a one-dimensional scatter plot.Here it is preferred over a box plot as it helps to show a multi-modal distribution.The thickness of a "bean" is a function of the frequency of the specific value -that is, the thicker a "bean" is for a value, the relatively higher the number of grid points having that value.The values shown are the Pearson's correlation coefficients which are based on the linear least-squares trend fit.Correlation values range from −1 to +1.Values closer to

Figure 8 .Figure 9 .
Figure 8. Spatial distribution of statistically significant (95 % significance level) partial correlation between de-trended annual GPP (GPPsat) and de-trended summer values of environmental variables (a) temperature, (b) precipitation, and (c) cloud cover.Negative correlations are shown with shades of red and positive correlations are shown in shades of blue.

Table 1 .
Biome property look-up table (BPLUT) for GPP algorithm with ERA-Interim and NDVI as inputs.The full names for the University of Maryland land cover classes (UMD_VEG_LC) in the MOD12Q1 data set are evergreen needleleaf forest (ENF), evergreen broadleaf forest (EBF), deciduous needleleaf forest (DNF), deciduous broadleaf forest (DBF), mixed forests (MF), closed shrublands (CS), open shrublands (OS), woody savannas (WS), savannas (SVN), grassland (GRS), and croplands (Crop).defined as the area between 15 • E longitude in the west, the Pacific coast in the east, 45 • N latitude in the south, and the Arctic Ocean coast in the north.The total area of this region is 22.4 million km 2 .Land cover distribution for the region is drawn from the Moderate Resolution Imaging Spectroradiometer (MODIS) MCD12Q1 Type 5 land cover product for the year 2007, available online at https://lpdaac.
Friedl et al. (2002)/data_pool from Land Processes Distributed Active Archive Center (LP DAAC), Sioux Falls, South Dakota, USA.The product provides global land cover at 1 km spatial resolution, produced from several classification systems, principally that of the International Geosphere-Biosphere Programme (IGBP).Friedl et al. (2002)describe the supervised classification methodology which leveraged a global database of training sites interpreted from highresolution imagery.The GPP products used in this study (described below) use a static land cover (LC) classification to define biome response characteristics over the study record.

Table 2 .
Details of the flux towers whose GPP data have been used to validate the satellite NDVI-based GPP data.The spatial distribution of these flux towers is shown in Fig. 1.

Table 4 .
Trend statistics for annual monthly averages of environmental variables.The first and second columns list the fraction of the region with significant (95 % significance level) positive trends and negative trends respectively.The third column is the regional mean trend of the variables per decade.The fourth column is the coefficient of variation, estimated as the distribution mean divided by the standard deviation.

Table 6 .
Connection between annual GPP of northern Eurasia (GPPsat) and summer values of environmental variables shown as percentage of the study area with statistically significant (95 % significance level) positive and negative partial correlation coefficients.