Environmental controls on the greening of terrestrial vegetation across northern Eurasia

Terrestrial ecosystems of northern Eurasia are greening, 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 Gross Primary Productivity (GPP) derived from satellite data across northern Eurasia for 5 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 regional median, 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 analyzed, even without a sustained increase in precipi10 tation, 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 south-west of the study area. For this region, precipitation positively correlates with GPP, as does cloudiness. This shows that the south-western part of northern Eurasia is relatively more vulnerable to drought than other areas. While our results further advance the notion that air temperature is 15 the dominant environmental control for recent GPP increases across northern Eurasia, the role of precipitation and cloudiness can not be ignored.


Introduction
Several analyses of Normalized Difference Vegetation Indices (NDVI) data derived from satellite remote sensing have pointed to greening, i.e. 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; Ove).Precipitation increases have also been observed over both North America and Eurasia over the past century (Nicholls et al., 1996;Gro).Urban et al. (2014) describe the co-occurrence of these climatic and ecosystem changes.Here we investigate greening 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.
Gross Primary 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 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).Highlighting the importance of precipitation, 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).With warming, low temperature constraints to productivity have re-laxed (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 greening 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 it's 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 analyze 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 greening of northern Eurasia.Therefore, we assess in this study how vegetation 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 variables on GPP; (3) identify the seasonality of the variables; (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, 2007), defined as the area between 15 • E longitude in the west, the Pacific coast in the east, 45 The product provides global land cover at 1 km spatial resolution, produced from several classification systems, principally that of the International Geosphere-Biosphere Program (IGBP).Friedl et al. (2002) describe the supervised classification methodology which leveraged a global database of training sites interpreted from high-resolution 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.So 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 center of the domain (Fig. 1).

Vegetation productivity -long term data
Gross Primary Production (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 process-based 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 vapor pressure deficit (VPD).GPP is derived on a daily basis as (Running et al., 2004;Zhang et al., 2008): where ε is a light use efficiency 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).The NDVI is from two datasets: (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) Vegetation Index and Phenology (VIP) (Didan, 2010;Barreto-Munoz, 2013) (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 sub-optimal temperature and moisture conditions represented by respective daily T min and VPD inputs.T f and VPD f are defined using simple linear ramp function (Yi et al., 2013;Heinsch et al., 2006), and minimum and maximum environmental constraints defined for different biome types (T mn_min and T mn_max , VPD min and VPD max ).extremely low productivity and technical problems of remote sensing make for a particular challenge in estimating 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 dataset has been examined in several recent studies.Analyzing trends in growing season 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.This deviation suggests data quality issues.The GIMMS dataset is based on the NOAA-AVHRR long term time series record, which is comprised of two different sensors for two different time periods, which leads to the difficulty of among-instrument inter-calibration and complicates the use of AVHRR NDVI data in long-term time series trend analysis (Pinzon and Tucker, 2014).This issue generates the need for a seamless and consistent sensor-independent record of land surface vegetation index and phenology.Used here in parallel with the GIMMS3g dataset, the Vegetation Index and Phenology (VIP) dataset, whose data processing applies a different correction scheme from that of GIMMS3g (Fensholt et al., 2015), is one such product.It is produced through continuity algorithms which enable the seamless conversion of land surface vegetation records from AVHRR, MODIS and VIIRS (Didan, 2010).The limited use and validation of the latter product necessitates the use of the ensemble mean of the two datasets.

Flux tower data
To validate the satellite based GPP estimates we use gap-filled daily tower GPP data at ten flux tower sites distributed across northern Eurasia, available for different periods of time.Details of the individual towers are provided in Table 4.The data, generated using the eddy covariance measurements acquired by the FLUXNET community, has been 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 to validate the satellite derived GPP for the entire study area, on a grid-by-grid basis.Upscaling of the FLUXNET observations was performed using the machine learning technique, model tree ensembles (MTE), a dataset provided by 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 dataset 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 & cloudiness
Monthly values of 2 m air temperature (in

Fire
Fire is represented by the 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.
Primary satellite-based data source is the MODIS surface reflectance imagery (Giglio et al., 2010).

Spatial interpolation
Data not on a 0.5 degree 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.

Validation
The GIMMS-GPP and VIP-GPP are evaluated against the tower-based data after grids and temporal ranges of the satellite data derived GPP, corresponding to the respective ten flux towers (Table 4) have been extracted.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 and observed values and it's 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 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 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 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 modeled 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) 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 datasets agree with each other; (5) Spatially explicit, pixel-by-pixel validation using the upscaled GPP data from FLUXNET observations (described in SubSect.2.1.3)using correlation and difference maps for the entire period.

Trend analysis
Temporal changes in each 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 (10yr −1 ) from the monthly values (mon −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 preserving only the inter-annual variability.

Correlation
We use the Karl Pearson's 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 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 change 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 analysis (R) can explain this dependence of one variable on another keeping the sign of the relationship (+/−) intact (Aldrich, 1995).

Validation of satellite derived GPP
The GIMMS-GPP and VIP-GPP, as well as their ensemble mean (GPPsat) are individually validated against the flux tower-based GPP data using Pearson's correlation coefficient, percent bias as well as the Nash-Sutcliffe normalized statistic.Scatter plots (Figure 2) show that GPP derived from satellite NDVI are generally higher than the tower-based GPP at the flux tower sites having comparatively lower productivity (and vice versa).Moreover, the agreement is higher at lower productivity sites than at higher productivity sites.Though  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 advanced very high resolution radiometers (AVHRR) and found the correlation to be highly variable across the sites, though consistently positive.Remarkably strong correlations 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).Comparing 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, correlation was found to be positive.LUE algorithms similar to the one used in this study for the generation of GPP datasets from satellite NDVI, produce favorable GPP results relative to daily tower observations, with a strong positive correlation (Yi et al., 2013;Yuan et al., 2007;Schubert et al., 2010).

Temporal changes in GPP
Across the study domain, regionally averaged GPPsat GPP exhibits a trend of 2.2 (±1.4) g C m −2 month −1 decade −1 .)) reveal the difference between the two datasets, which is highest at the beginning of the study period.The nature of increase in GPP is also different for the two datasets with the rise in one being more linear than the other.A possible explanation for the differences in the two datasets is discussed in Section 2.1.2.Examining the seasonality of GPP trends (of GPPsat) (Fig. 5), we find that among the seasons the summer trend is greatest.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 having a 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 process based 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 at the high latitudes of this region suggests that 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. 6(a.) shows 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 southwestern 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. 6(b.)), the range of trends for this region, from minimum to maximum, in highest for summer.The fraction of the region experiencing significant increases in annual precipitation is about three times the area experiencing significant decreases.The significant positive trends are located in the north-eastern and western parts (mainly boreal forests) while the 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.The latter is found to have a negative trend.However, similar to precipitation, the spatial standard deviation is very high, implying a high spatial variability and that the regional mean does not represent the trends of the whole region.Unlike precipitation, a greater fraction of the region sees a significant negative cloudiness trend or a significant clear sky trend (Fig. 6(c.)).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. 6(d.)), is not significant.
Recent studies have reported similar changes in these environmental variables.For the period of 1979 to 2005, Trenberth et al. (2007) found temperature trends over the region range from 0.3-0.7 • C decade −1 and for most regions of the higher latitudes, especially from 30 to 85 • N, significant positive precipitation trends have occurred.Contrary to the cloud cover trend we find here, studies reported in AR4 suggest an increased 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, is dissimilar from the other environmental variables in that it spans only 14 years of the 29 years, and is spatially non-uniform, 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 analyze 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.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 to annual GPP for only 1.7% to 3.4% (depending on season) a small fraction 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 GPP.The contrast between summer and the other seasons is strongest for temperature, highlighting the importance of summer temperatures to annual productivity.Fig. 7 reveals that the relationships between annual GPP and 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 bi-modal 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 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.
Since GPP trends are highest in summer (Fig. 5), the peak of the growing season, we are interested more in the impact of the summer values of environmental variables on the annual GPP since the terrestrial vegetation is likely to be more responsive to variations of environmental variables in summer than 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. 8(a.);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 the tundra as well as parts of the boreal forests.
Negative correlations occur across 2% of the region, largely in the south, where Eurasian steppes can be found.For other parts of the year (maps not shown for spring and fall correlations but distributions represented in Fig. 7) significant negative correlations become more spatially disperse, while significant positive correlations are limited to the center and west of the region for spring, becoming more disperse in autumn.Determining the partial relationship between annual GPP and summer precipitation, Fig. 8(b.) 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 fall precipitation is predominantly negative.The spatial correlations for summer cloudiness and summer precipitation are similar, (Fig. 8(c.)) though the area under significant correlation is comparatively less.The area of negative correlation is about nine times that of the positive correlation (Table 6).
Compared to summer, area under significant positive correlation is higher for spring while area under negative correlation is higher for fall (maps not shown).
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 might lead to greater water stress (Gates, 1964;Wiegand and Namken, 1966;Jackson et al., 1981).Decreasing precipitation would increase water stress.Moreover, 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 one hand, cause a decrease in direct radiation and increase the diffuse radiation.This would increase GPP through higher light use efficiency (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.On the contrary, 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 "greenness" across northern Eurasia, is possibly a reflection of continued summer warming in the absence of sustained increases in precipitation (Buermann et al., 2014;Zhou et al., 2001).Wang et al. (2014) documented a positive relationship between sunshine duration (equivalent to the inverse of cloud cover) and vegetation greenness.

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 temperature-precipitation 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 should be leading to increasing water stress, evidence of which is noted in a subset of the region.Indeed, some 2.4 % of the area in the southern parts of the study area (Fig. 8(a.) shows significant negative partial correlation between annual GPP 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 correlations between precipitation and cloud cover are predominately positive, with the significant correlations spread out across the region.As described in Section 3.4, the correlations between precipitation and cloud cover helps to explain why spatial distributions of the correlation coefficients of precipitation and cloud cover with GPP are similar.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.But Figure 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 correlation analysis is performed after detrending (removing long term 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 southwestern 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 datasets derived from GIMMS3g and VIP NDVI data shows 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 datasets.The regional mean trend for the ensemble mean GPP is 2.2 (±1.4) g C m −2 month −1 decade −1 .The grid analysis is consistent with results of prior studies which have suggested that air temperature is the dominant environmental variable influencing the increases in productivity across the northern high latitudes.Examining partial coefficient 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 s, 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 southwest, 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 any significant correlation, and consequently an impact of fire on GPP was ignored from 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, with warming having occurred with precipitation failing to keep up with the increased demand.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    et al., 1996;Sorooshian et al., 1993); (3) Nash-Sutcliffe model efficiency coefficient (Nash and Sutcliffe, 1970).Its values 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 modeled 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.columns list the fraction of the region with significant (95 % significance level) positive trends and negative trends respectively.The 3rd column is the regional mean trend of the variables per decade.
The 4th column is the coefficient of variation, estimated as the distribution mean divided by the standard deviation.Cover Product (Friedl et al., 2002).The details of the flux tower sites are listed in Table 4.  Unit (CRU TS 3.21) dataset (Harris et al., 2014).Burnt area data from the Global Fire Emissions Database (GFED) (Giglio et al., 2010).
Partial correlation (R) served mean for the respective flux tower sites.Spatially explicit validation of GPPsat reveals that the correlation is high and statistically significant for almost the entire study area (Fig.3 (a.)).GPPsat underestimates in the boreal forests of the western parts of northern Eurasia and overestimates in the Eurasian steppes to the south of the study area(Fig.3 (b.)).

Figure 4
Figure 4(a.)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 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.
the flux towers whose GPP data have been used to validate the satellite NDVI based GPP data.The spatial distribution of these flux towers are shown in Figure 1.845 Table 3: Validation of GIMMS3g and VIP GPP datasets along with their ensemble mean using flux tower GPP from ten flux tower sites across northern Eurasia.The spatial distribution of the flux tower sites is shown in Fig. 1.Validation was carried out using: (1) Pearson's product moment correlation.It is a measure of the linear dependence between the simulated and observed GPP and it's 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

Figure 1 :
Figure 1: Simplified land cover for northern Eurasia for year 2007 overlaid with the spatial distribution of the ten flux tower sites whose GPP (Gross Primary Productivity) data was 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; (ii) wooded, i.e. plants with wood as it's structural tissue, which includes the boreal forests appearing in the middle and extend from the western to the eastern boundary.This land cover map has been derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) Type 5 Land

)Figure 2 :
Figure 2: 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 deviate from the 1:1 perfect relationship.

Figure 3 :Figure 4 :)Figure 5 :)Figure 6 :
Figure 3: Spatially explicit validation of GPPsat using upscaled FLUXNET observations.(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.(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 7 :Figure 8 :
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 and one dimensional scatter plot.Here it is preferred over a box plot as it helps to visualize 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 is 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 indicates strong correlation while those closer to 0 indicates weak correlation.The color 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.
• N latitude in the south and the Arctic Ocean coast in the north.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 Land Cover Product Type 5 Land Cover Product for the year 2007, available online at https://lpdaac.usgs.gov/data_access/data_poolfrom Land Processes Distributed Active Archive Center (LP DAAC), Sioux Falls, South Dakota, USA.

Table 1
(Jones and Harris, 2013)5)) 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 dataset, 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/(JonesandHarris,2013).Although the LUE based GPP model does not use precipitation as an input, we assume that precipitation is a good measure of water supply to the vegetation and thus analyze 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 GPP 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 SubSect.2.1.2,lower reliability and availability of satellite NDVI observations and associated GPP data for the winter months leads us to focus on spring, summer, and autumn seasons.
Table 3 lists all the validation statistics, we focus only on those of annual GPPsat since this is what we investigate for the rest of the study.The correlation coefficients are all positive and high (0.7 for annual GPPsat), percent bias is predominantly nega- tive (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 ob-

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 dataset are Evergreen Needleleaf Forest (ENF), Evergreen Broadleaf Forest (EBF),

Table 2 :
Details of

Table 4 :
Trend statistics for annual monthly averages of environmental variables.The 1st and 2nd

Table 5 :
Medians of the distributions of the relative partial significant contribution (R 2 -95% sig-In these cases the factors behind the unexplained attribution are not identified. nificance) of each de-trended environmental variable (except fire) of each season to the inter-annual variability in de-trended annual GPP (GPPsat).In each case the total contribution may not add up to 100 %.

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.