Impact of changes in GRACE derived terrestrial water storage on vegetation growth in Eurasia

We use GRACE-derived terrestrial water storage (TWS) and ERA-interim air temperature, as proxy for available water and temperature constraints on vegetation productivity, inferred from MODIS satellite normalized difference vegetation index (NDVI), in Northern Eurasia during 2002–2011. We investigate how changes in TWS affect the correlation between NDVI and temperature during the non-frozen season. We find that vegetation growth exhibits significant spatial and temporal variability associated with varying trend in TWS and temperature. The largest NDVI gains occur over boreal forests associated with warming and wetting. The largest NDVI losses occur over grasslands in the Southwestern Ob associated with regional drying and cooling, with dominant constraint from TWS. Over grasslands and temperate forests in the Southeast Ob and South Yenisei, wetting and cooling lead to a dominant temperature constraint due to the relaxation of TWS constraints. Overall, we find significant monthly correlation of NDVI with TWS and temperature over 35% and 50% of the domain, respectively. These results indicate that water availability (TWS) plays a major role in modulating Eurasia vegetation response to temperature changes.


Introduction
High latitude boreal and Arctic biomes are changing in response to recent climatic warming (Overland et al 2008). Recent atmospheric CO 2 anomalies, satellite normalized difference vegetation index (NDVI) records, satellite data driven vegetation productivity estimates, and ground-based observations indicate regional declines in Northern ecosystem productivity since the mid-1990s that may be a biophysical response to temperature induced drought stress , Piao et al 2008, Ma et al 2012, Yi et al 2013. A diverse response of Northern ecosystems to regional warming and drought has been reported from limited field experiments (Lloyd and Bunn 2007, Schwalm et al 2010, Yi et al 2010, Peng et al 2011, Tei et al 2015, while extension of these findings to the pan-boreal/Arctic region is constrained by the limited extent of these studies and a sparse regional measurement network. During the past 30 years, published studies have shown positive correlation between vegetation productivity and temperature changes during the growing season in the Arctic (e.g. Zhou et al 2001, Lucht et al 2002, however, more recent studies suggest a weakening of this relationship (e.g. Buermann et al 2014, Piao et al 2014, indicating complex patterns of change in vegetation growth in relation to regional temperature trends, including predominant greening in tundra areas, and widespread productivity declines in boreal forests attributed to water stress (Zhang et al 2008, Bhatt et al 2010, Beck and Goetz 2011, Ma et al 2012, Berner et al 2013, Guay et al 2014, Walker and Johnstone 2014. It is difficult to evaluate trends in vegetation response to water constraints at Northern high latitudes, because of uncertainty in regional precipitation and other water budget indicators due to intrinsic bias in available satellite and reanalysis data, and sparse ground measurements (Serreze et al 2006, Zhang et al 2009, Rawlins et al 2010. Satellite-derived soil moisture estimates have large uncertainties in the Northern regions (Mladenova et al 2014). In fact satellite soil moisture retrieval errors are large due to the presence of snow, ice, surface water and dense canopies (Naeimi et al 2012, Högström et al 2014. In previous studies, local and regional responses of vegetation to moisture variability have typically been inferred indirectly from correlations with model based temperature and precipitation, or atmosphere vapor pressure deficit (e.g. Nemani et al 2003, Zhang et al 2008, or from reanalysis based moisture indices (e.g. Berner et al 2013).
Since 2002 GRACE satellite measurements have provided regional estimates of monthly changes in terrestrial water storage (TWS) (Swenson et al 2006, Reager and Famiglietti 2009, Famiglietti et al 2011, i.e. of water content both on and below the land surface. GRACE data provide a unique tool for observing trends in water availability (Rodell et al 2009, Velicogna et al 2012. These data have been used to investigate regional climate variability, including drought related impacts (Chen et al 2013, Thomas et al 2014, and to study water constraints on vegetation growth in Australia . In Eurasia, previous studies have shown that GRACE TWS changes are consistent with observed changes in net precipitation (Landerer et al 2010, Velicogna et al 2012. Here we investigate the vegetation response to temperature changes and water storage availability from 2002 to 2011 in the Ob, Yenisei and Lena basins, which constitute 70% of all riverine freshwater inputs to the Arctic Ocean (Gordeev et al 1996). To study the relationship between vegetation growth and water storage and temperature changes, we use MODIS NDVI, GRACE TWS estimates, and ECMWF (European Centre for Medium-Range Weather Forecasts) ERA-Interim 2 m air temperature (Dee et al 2011). We discuss the rationale of using GRACE TWS as a proxy for plant available water supply. We study the spatial pattern of trends in NDVI, TWS, and temperature (T) and evaluate the trend and inter-annual relation of NDVI changes with temperature and water availability during 2002-2011. We conclude on discussing the role of water storage, as a proxy for plant-available water supply, in modulating the impact of temperature changes on regional vegetation productivity in Northern ecosystems.

Data and methodology
We use 113 monthly GRACE solutions from the Center for Space Research at the University of Texas (Tapley et al 2004), between August 2002 and December 2011. Each gravity solution consists of spherical harmonic coefficients, C lm and S lm , up to degree, l, and order, m, 60. We use monthly values of C 20 coefficients from satellite laser ranging (Cheng et al 2013) and include degree-1 coefficients calculated as described by Swenson et al (2008). The GRACE data directly measure monthly TWS changes, because this is the largest source of mass change within our area of interest; other mass changes such as glacial isostatic adjustment (GIA) are of much lower magnitude. TWS anomalies are calculated relative to the period August 2002-December 2011. To separate seasonal variability from inter-annual variability and trends we apply a 13 month smoothing to the GRACE monthly Stokes coefficients. To do this for each 13 month window, we simultaneously solve for annual, semiannual, 3 month period signals, a constant and a linear trend. We assign the sum of the constant term and the linear trend to the center point (i.e., the 7th month) of each time window. This moving average scheme yields a smoothed time series where all seasonal variations are removed. This smoothing procedure is not applied to the first and last 6 months of record because a one-year cycle is required to extract the seasonal signal. The GRACE data are also corrected for the GIA signal following Paulson et al (2007). To reduce the random error components, we apply a Gaussian smoothing with a 400 km radius (Wahr et al 1998) and then we generate monthly evenly spaced latitude-longitude grid. Linear trends are calculated using these monthly maps.
Water in seasonal snow cover is not directly available for vegetation growth. We therefore remove estimates of snow storage from the GRACE monthly TWS time series. We estimate the snow cover contribution using 25 km EASE-Grid monthly satellite derived snow water equivalent (SWE) data from the Advanced Microwave Scanning Radiometer on EOS Aqua, AMSR-E, (http://nsidc.org/data/amsre/) (Derksen et al 2003). The SWE retrieval accuracy is generally higher over flatter land areas with less vegetation cover. This is the case of the analyzed region that has relatively flat topography and is mainly covered by grassland, boreal taiga and tundra vegetation cover. Velicogna et al (2012) found good agreement between changes in snow mass from the AMSR-E SWE record and station-based cold season precipitation trends for the region. We use these data to remove the SWE signal from the GRACE TWS estimates. We only remove the SWE of the non-melted snow; hence the residual TWS includes the water storage from melted snow. In the remainder of the paper, TWS refers to the signal corrected for SWE.
In principle, the SWE-corrected GRACE TWS includes components such as groundwater that is not directly within the vegetation root zone. These other components vary on longer time scales, whereas short time scale (e.g. monthly) variations in water storage are related to soil water changes linked to vegetation growth (Vicente-Serrano et al 2010. In addition, previous studies have shown that in Eurasia, variations in SWE-corrected GRACE TWS are strongly correlated with net precipitation at the monthly and inter-annual time scale (Landerer et al 2010, Velicogna et al 2012, which has also been used as a proxy for water available for plant growth. Therefore we use the monthly SWE-corrected GRACE TWS estimates as a surrogate for plant-available water supply. We use ERA-Interim monthly 2 m height temperature (T) (Dee et al 2011) and MODIS derived NDVI records. The NDVI from satellite optical-IR remote sensing is sensitive to vegetation canopy greenness and has been used as a surrogate indicator of photosynthetic canopy growth changes over Northern land areas (Bogaert et al 2002. We use a monthly composited MODIS-MOD13Q1 NDVI dataset with 1°×1°spatial resolution (Huete et al 2011) to define vegetation canopy changes over the study domain. All datasets are processed as the GRACE data, i.e. anomalies are calculated relative to the common period (August 2002-December 2011), converted to spherical harmonics, truncated to degree 60, spatially averaged using a 400 km smoothing and converted back to the spatial domain on a regular 1°l ongitude-latitude grid. We calculate TWS, T, and NDVI trends. Trends are considered statistically significant when they exceed the respective 1−σ error (figure 1). We use the global Earth System Data Record of daily landscape freeze-thaw status derived from satellite microwave remote sensing to define daily landscape freeze-thaw status over all Northern vegetated land areas (Kim et al 2011). For each year, we define the non-frozen (NF) season at each grid cell as described in Kim et al (2014).
To investigate the impact of temperature and water storage changes on vegetation growth, we examine the partial correlation between NDVI and T (NDVI-T correlation) and NDVI and TWS (NDVI-TWS correlation) for the NF season. The partial correlation removes the inter-dependency between T and TWS. This study focuses on the relationship between trends and inter-annual NDVI variations with T and TWS. The trends and inter-annual changes are small compared to the seasonal variability, hence before calculating the partial correlations, we remove the seasonal signal from each time series as described earlier; otherwise it would dominate the correlation analysis. In the remainder of the paper, we use R NDVI-T and R NDVI-TWS to indicate the NF season partial correlation between NDVI and T and between NDVI and TWS, respectively. To identify whether trend or interannual variability is driving the correlations we calculated the correlations for de-trended time series and compare them with the one for the original time series.
Due to the spatial resolution of the GRACE data, TWS estimates can only be used to study vegetation growth changes at the regional scale, i.e. 400 km. To facilitate the interpretation of the results, we divide the study domain into a set of four regions. Each region corresponds to an area with significant NDVI-T or NDVI-TWS correlation and with a consistent trend in NDVI (figures 1(c)-(e)). Within each region, we analyze the spatially averaged correlation coefficients R NDVI-T and R NDVI-TWS to determine the primary moisture or temperature constraints and we discuss the results as a function of vegetation type within the region.
Change in fire disturbance regime also influence vegetation growth patterns and may affect our results based on correspondence between TWS, T and NDVI. To evaluate the impact of fire on vegetation growth, we use area burned data from the fourth-generation global fire emissions database (Giglio et al 2013) and for each vegetation type we compare NDVI trends on burned and unburned areas. show R NDVI-T and R NDVI-TWS calculated for the entire analyzed region; however, in our analysis we consider only areas where the partial correlation coefficient is statistically significant at more than 98% confidence level (p<0.02). We identify climatic constraints on vegetation activity for 68% of the analyzed region. In the remaining 32% of the region, we do not find significant correlation between NDVI and T or NDVI and TWS.

Results
Within the analyzed domain, NDVI correlates well (p<0.02) with TWS and T over 35% and 50% of the area, respectively. R NDVI-TWS is positive (p<0.02) in W and N Lena, SE Yenisei and SW and SE Ob basins (figure 1(e)) and is not significant elsewhere. We find a positive R NDVI-T in the E Lena (with R NDVI-T 35% larger than R NDVI-TWS ), N and SE Ob (with R NDVI-T 29% larger than R NDVI-TWS in SE Ob), and S Yenisei basins ( figure 1(d)). By comparing partial correlation results for the de-trended and original time series, we find that over the entire domain, the NDVI-T correlation is mostly driven by the inter-annual variability. In contrast, for the NDVI-TWS correlation, we find that in the Lena basin, the correlation is driven by the trend, whereas in the SW Ob basin, is driven by both the trend and the inter-annual variability, and in the SE Ob by the inter-annual variability.
Spatially averaged NDVI-T and NDVI-TWS correlations (p<0.02) and the land cover partitioning for the four identified regions are summarized in table 1. The identified dominant water and temperature constraints are summarized in figure 2. The time series of NDVI, T and TWS are shown in figure S2. We then use these results to calculate the areal proportion of relative water and temperature constraints for each vegetation type in the entire analyzed domain (table 2).
Region 1, the SW Ob basin, is cooling and drying, with a significant decrease in NDVI (figures 1(a)-(c). The land cover in this region predominantly consists of relatively arid grassland (97%) (table 1). R NDVI-TWS is positive over the entire region (average R 0.48, NDVI TWS -= p<0.02), while R NDVI-T is not significant. These results are consistent with a regional dominance of TWS constraints on vegetation productivity.
Region 2 covers the SE Ob and S Yenisei basins, and has a mixed land cover of 44% grassland, 24% temperate forest, 23% boreal forest and 8% tundra. The region accounts for 83% of the temperate forest and 53% of the grassland biomes of the entire study domain. This area is cooling and wetting, with a small negative NDVI trend significant only in the South (figures 1(a)-(c). R NDVI-T is positive (average R 0.53, T NDVI-= p<0.02) over the entire region,  R NDVI-TWS is significant in the temperate forest area in SE Ob (R NDVI-TWS =0.35, p<0.02). Here, T is the primary constraint on vegetation productivity. Region 3 (figure 1) includes almost the entire Lena basin and small portions of the NE and SE Yenisei. These two small areas are included here rather than in region 2 as they show an increase in NDVI consistent with the trend in region 3. Here the land cover is composed of boreal forest (83%), tundra (13%) and grassland (3%). The entire region shows TWS wetting, but with a spatially variable temperature trend: the E Lena basin is warming, the NE and SE Yenisei are cooling, and there is no significant T trend in the rest of the region (figures 1(a)-(b). In the E Lena basin, the average R NDVI-TWS and R NDVI-T are 0.46 and 0.61, respectively. R NDVI-T is 33% larger than R NDVI-TWS , indicating, in agreement with Piao et al (2011), that T has a generally stronger impact on vegetation in this region. In the remaining portion of the basin and in the NE Yenisei, however, R NDVI-TWS is positive (average R NDVI-TWS in these areas is 0.46, p<0.02), and the NDVI-T correlation is not significant. Here we observe an NDVI greening response to TWS wetting, indicating a dominant TWS constraint to vegetation activity.
The SE-Yenisei consists mostly of boreal forest (69%) and grassland (17%). This area is cooling and wetting, and features a small NDVI increase (1.2%/ decade on average, p<0.32). R NDVI-T is positive (average R NDVI-T =0.53, p<0.02) and 29% larger than the R NDVI-TWS (average R NDVI-TWS =0.41, p<0.02). Here temperature is the main constraint on vegetation growth.
Region 4 includes a portion of the NW Ob basin and a small area in the N Yenisei. Boreal forest and tundra cover 76% and 24% of the area, respectively (table 1). This region features a decreasing T trend, a small increase in TWS and no significant NDVI trend (figures 1(a)-(c)). R NDVI-T is positive (average R NDVI-T =0.51, p<0.02) over the entire region, while the TWS-NDVI correlation is not significant (table 1), indicating an overall T constraint.
Overall the analyzed domain includes 63% boreal forests, 20% grasslands, 10% tundra and 7% temperate forests. Within the tundra ( figure 1(f)), NDVI is significantly (p<0.02) correlated with T in 74% of the region and with TWS in only 8% (p<0.02) (table 2), suggesting cold temperature as the main constraint on vegetation growth. We find no significant correlation in 18% of tundra. Boreal forests show a diverse response to climatic variations; significant correlation of NDVI with T and TWS occur over 40% and 17% of boreal areas, respectively, and no correlation elsewhere (table 2). Overall temperature-constrained boreal forests are distributed over the NW Ob, SE Yenisei and E Lena, while TWS-constrained biomes  are found in the central and W Lena. Grasslands also show mixed climatic controls. TWS is the dominant constraint on NDVI growth in the SW Ob basin where we observe cooling and drying; this area accounts for 32% of the grasslands. While in the SE Ob and S Yenisei, where we observe cooling and wetting, we find significant correlation with T; these regions account for 58% of the grasslands. The analyzed domain also includes a relatively small area of temperate forests (S Ob and S Yenisei); 83% of this area shows significant NDVI-T correspondence, indicating a primary cold temperature constraint in temperate forests. When we evaluate the impact of fire disturbance on vegetation growth, we find mean annual burned area of 1.9% in grassland, and less than 0.5% for the other land covers in the analyzed domain. Over the grassland, we find similar NDVI trends (∼ −2%/decade) in both burned and unburned areas, indicating that fire at the regional scale is not a dominant factor controlling vegetation growth, consistent with previous studies in Eurasia (Buermann et al 2014).

Discussion
During the analyzed period, the greatest NDVI gains occur where we observe warming and wetting trends, and boreal forest is the locally dominant biome (Lena basin, region 3, average NDVI trend=2.4%/decade). The greatest NDVI losses are associated with cooling and drying trends, and occur in regions dominated by grassland (74%) (SW Ob basin, region 1, average NDVI trend=−1.9%/decade). NDVI changes are small or not significant in regions that are cooling and wetting. Our study domain does not sample warming and drying condition, but it is likely that warming would exacerbate drying-induced TWS constraints and reduce NDVI. The spatial pattern of these changes provides an indication of how regional vegetation growth may respond to future climate conditions. Regional climate model (IPCC AR5) projections indicate warmer climate conditions over Eurasia (e.g. Miao et al 2014), and an overall increase in precipitation with possible drying trends in the Southern regions (Monier et al 2013, Sillmann et al 2013. Our results imply that projected warmer, wetter trends could promote widespread productivity increases and conditions more suited to temperate forest than colder boreal and tundra areas; however, where warmer, drier conditions occur we may expect to see productivity declines and conditions better suited to drought tolerant grassland over tundra and forested areas (Wang and Overland 2004, Lenton et al 2008. Our results indicate that fire, at the regional scale, is not a dominant factor controlling vegetation growth, consistent with previous studies in Eurasia . However, area burned and fire severity have been increasing in Eurasia for the last decade Most of the Lena and a large portion of the Yenisei basins are underlain by permafrost (Brown et al 1997. Permafrost degradation may play a significant role in modulating the vegetation response to ongoing climate change in this region. Model projections predict shrinkage of the areal extent of permafrost and deepening of the active layer which has the potential to significantly impact plant-available soil moisture (Lawrence and Slater 2005, Anisimov 2007, Saito et al 2007, Lawrence et al 2008. The interaction between permafrost degradation, climate change and vegetation growth is a topic of ongoing research; the mechanism, timing and extent of this interaction is still uncertain (e.g. Tchebakova et al 2009, Walker et al 2009, Jorgenson et al 2010, Nauta et al 2014. Long-term monitoring is needed to determine if deepening of the active layer along with increasing drainage of root-zone soil water is leading to shifts in vegetation patterns and climatic control factors influencing productivity (e.g. Slater 2005, Jorgenson et al 2013).
During 2002-2011, T exerted a major constraint on NDVI over a large portion of Eurasian grassland and temperate forest areas, mostly in the SE Ob and S Yenisei basins, while previous studies, based on longer, multi-decadal observations, indicate that these biomes are mainly water limited (e.g. Piao et al 2011, Kim et al 2014. We attribute this difference to regional wetting and cooling trends, observed during the study period, that have been associated to a decadal and quasi-decadal climatic variability , Bhatt et al 2013, rather than a longer term climate change response (Cohen et al 2009, Piao et al 2011, Kim et al 2014. Our results indicate that the regional response of NDVI to temperature is influenced by TWS changes and associated water supply controls to vegetation growth. An observed increasing TWS trend is consistent with an apparent reduction in water supply constraints, so that T cooling represents the dominant constraint driving vegetation growth reductions in the region. These results are consistent with recent findings that soil moisture, as another water supply proxy, modulates temperature-productivity correlations and ecosystem net carbon uptake in Northern high latitude regions (Piao et al 2014, Yi et al 2014, while drought-induced water stress has been the dominant cause of observed reduction in vegetation growth for Canadian boreal forests (e.g. Ma et al 2012) despite regional warming and associated reductions in cold temperature constraints to productivity.
Summer cloud cover can also impact vegetation growth by limiting incident solar radiation on vegetation (Nemani et al 2003). Recent studies have linked changes in summer cloud to regional variability in temperature trends and vegetation growth (Tang andLeng 2012, Bhatt et al 2013). Piao et al (2014) have shown that, in the study domain, the NDVI correlations with shortwave radiation are weaker than those with temperature or precipitation. Still increase in cloud cover, observed during the analyzed period, may have contributed to the negative NDVI-TWS correlation observed in the Southern and Northern end of the study domain ( figure 1(e)) (e.g. Tang and Leng 2012). Future studies will incorporate additional observations to further investigate the relation between changes in TWS and summer cloud cover, and the related impact on vegetation growth.
Seasonal snow cover may affect vegetation growth as winter snow accumulation provides thermal insulation of soil and melting snow water provides water supply to vegetation growth in the NF seasons  . In this region the NDVI-TWS correlation is driven by the trends, and we find a positive correlation (R=0.81, p<0.01) between annual maximum TWS and early NF season NDVI. This suggests that seasonal snow pack and water supply from snowmelt have an important influence on vegetation activity. Further study is needed to investigate this correlation.
Overall, we find that cold temperature constraints are the main influence on vegetation activity in 50% of the study domain, while water supply constraints are the dominant control on vegetation activity in 18% of the domain (table 2). GRACE data provide a regional unbiased estimate of TWS (e.g. Velicogna et al 2012) and an observation-based proxy for analyzing changing water supply controls to vegetation productivity. Longer TWS time series will be obtained from continuing GRACE and GRACE follow-on mission (scheduled launch August 2017) observations, enabling improved diagnosis of regional climate trends from large characteristic climate variability. New soil moisture observations from the NASA Soil Moisture Active Passive (SMAP) mission and continuing improvements to regional snow cover (SWE) datasets may enable further partitioning of GRACE based TWS into surface moisture changes directly affecting vegetation growth versus deeper groundwater that is less accessible to plants.

Conclusions
In this study we analyze satellite observations of monthly TWS changes from GRACE, and compare the results with MODIS NDVI and ERA-Interim T to clarify regional changes in vegetation productivity patterns and underlying water supply and temperature constraints over the three largest basins in Eurasia during 2002-2011. In agreement with previous studies, we observe large spatial and temporal heterogeneity in patterns of vegetation growth and underlying environmental controls (e.g. Zhang et al 2008, Beck and Goetz 2011. We find both greening and browning trends for tundra in Northern Eurasia, driven by temperature warming and cooling, respectively. In boreal forests, we identify an overall greening trend, caused by relaxation of cold temperature and water supply constraints. In temperate forest and grassland regions, cold temperature is the dominant climatic control due to the relaxation of water constraints. Overall, we find that recent changes in plant-available water (either decrease or increase in available water supply) are a critical factor modulating the sign and magnitude of the regional vegetation growth response to temperature changes. Despite the short length of the satellite record, we find complex patterns of NDVI growth response to T and TWS changes that follow large regional climate and land cover gradients and provide for a space-for-time based assessment of potential vegetation responses to future climatic conditions. Our results indicate that future vegetation changes in the region will increasingly depend on having an adequate moisture supply to support potential vegetation growth enabled by regional warming and associated relaxation of cold temperature constraints. Warmer and wetter conditions will likely promote an expansion of temperate forest areas and loss of tundra, while warmer and drier conditions may promote grassland expansion and loss of forests. However ongoing increase in atmospheric CO 2 will likely modify the radiation budget and alter presentday plant physiology, and future studies are needed to understand its impact on vegetation growth and the associated climatic constraints. Longer satellite observation records from GRACE and follow-on missions in combination with planned soil moisture observations from SMAP and satellite SWE will expand and improve our assessment of plant-available moisture changes and environmental trends.