Attribution of divergent northern vegetation growth responses to lengthening non-frozen seasons using satellite optical-NIR and microwave remote sensing

The non-frozen (NF) season duration strongly influences the northern carbon cycle where frozen (FR) temperatures are a major constraint to biological processes. The landscape freeze-thaw (FT) signal from satellite microwave remote sensing provides a surrogate measure of FR temperature constraints to ecosystem productivity, trace gas exchange, and surface water mobility. We analysed a new global satellite data record of daily landscape FT dynamics derived from temporal classification of overlapping SMMR and SSM/I 37 GHz frequency brightness temperatures (Tb). The FT record was used to quantify regional patterns, annual variability, and trends in the NF season over northern (≥45°N) vegetated land areas. The ecological significance of these changes was evaluated against satellite normalized difference vegetation index (NDVI) anomalies, estimated moisture and temperature constraints to productivity determined from meteorological reanalysis, and atmospheric CO2 records. The FT record shows a lengthening (2.4 days decade−1; p < 0.005) mean annual NF season trend (1979–2010) for the high northern latitudes that is 26% larger than the Northern Hemisphere trend. The NDVI summer growth response to these changes is spatially complex and coincides with local dominance of cold temperature or moisture constraints to productivity. Longer NF seasons are predominantly enhancing productivity in cold temperature-constrained areas, whereas these effects are reduced or reversed in more moisture-constrained areas. Longer NF seasons also increase the atmospheric CO2 seasonal amplitude by enhancing both regional carbon uptake and emissions. We find that cold temperature constraints to northern growing seasons are relaxing, whereas potential benefits for productivity and carbon sink activity are becoming more dependent on the terrestrial water balance and supply of plant-available moisture needed to meet additional water use demands under a warming climate.


Introduction
The high northern latitude (HNL) land areas have experienced greater surface air temperature warming than the global average in recent decades (Liu et al. 2007; Serreze and Francis 2006;Winton 2006). Recent environmental changes attributed to HNL warming include earlier and longer potential growing seasons (Nemani et al. 2003;Kim et al. 2012), vegetation greening Bunn and Goetz 2006;Hudson and Henry 2009) and productivity increases ), a northward shift in vegetation biomes , and tundra shrub expansion (McManus et al. 2012;Tape et al. 2012). HNL vegetation growth is primarily constrained by seasonal cold temperatures (Nemani et al. 2003;Friedlingstein et al. 2006;Qian, Joseph, and Zeng 2010), but recent reports indicate that widespread drought and wildfire disturbances exacerbated by continued warming have resulted in the frequent occurrences of tree mortality and declines in boreal productivity (Girardin and Mudelsee 2008;Goetz et al. 2005;Peng et al. 2011). Spatially extensive patterns of drought-induced vegetation growth decline have also been reported across the HNL domain (Goetz et al. 2012;Kim et al. 2012;Schaphoff et al. 2006;Zhang et al. 2008), including interior Alaska (Baird, Verbyla, and Hollingsworth 2012;Verbyla 2008), Canada (Ma et al. 2012;Peng et al. 2011), and Eurasia (Park and Sohn 2010;Piao et al. 2011).
Arctic tundra and boreal forests represent the dominant HNL biomes and cover approximately 25% of the global vegetated land surface ). Both biomes store a considerable amount of carbon in seasonally frozen (FR) ground that favours soil organic carbon (SOC) accumulation over microbial decomposition. The status of HNL terrestrial carbon sink activity under continued climate warming remains unclear due to uncertainty regarding the potential productivity benefits of lengthening growing seasons and increasing vulnerability of vegetation and SOC stocks to disturbances from drought, wildfire and insects, and increasing SOC decomposition and respiration from permafrost melting (McGuire et al. 1995;Zimov, Schuur, and Chapin 2006;Kimball et al. 2007).
The normalized difference vegetation index (NDVI) derived from satellite optical and near-infrared (NIR) remote sensing has been widely used to document vegetation responses to climate variability and global warming (Lloyd, Bunn, and Berner 2011;Raynolds et al. 2012;Xu et al. 2013). The NDVI is sensitive to canopy photosynthetic activity (Field, Randerson, and Malmström 1995;Huete et al. 2002;Tucker, Townshend, and Goff 1985), and has been commonly used to analyse regional patterns and temporal variability and trends in HNL vegetation growth Goetz et al. 2005;Piao et al. 2011). However, seasonal reductions in solar illumination and persistent clouds, smoke, and other atmospheric effects can degrade satellite NDVI retrieval accuracy over many areas, especially at higher latitudes (Alcaraz-Segura et al. 2010;Nagol, Vermote, and Prince 2009;Shuai et al. 2008).
Recently, Kim et al. (2012) developed a consistent global record of daily landscape freeze-thaw (FT) dynamics by temporal change classification of satellite passive microwave brightness temperature retrievals that are relatively insensitive to atmosphere and solar illumination effects (Jones et al. 2010;Kim et al. 2011;Naeimi et al. 2012). The satellite microwave FT signal senses the predominant FR or non-frozen (NF) status of water in surface soil and vegetation canopy layers Zhao et al. 2011), and has been linked to FR temperature constraints to vegetation productivity, landscape water mobility, and land-atmosphere carbon exchange over the HNL domain McDonald et al. 2004;Zhang et al. 2011). Kim et al. (2012) applied the satellite FT record and documented significant Northern Hemisphere trends of lengthening NF seasons consistent with recent global warming from 1979 to 2008. The FT results were compared against satellite (MODIS) optical-NIR-derived summer NDVI growth anomalies over a 9 year (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) record. They showed that longer NF seasons have generally promoted widespread vegetation canopy growth, with the relative benefits increasing at higher latitudes, and reduced or degraded productivity effects at lower latitudes.
For this investigation we hypothesized that the regional patterns in the HNL canopy growth response to longer NF seasons are dependent on underlying cold temperature or moisture constraints to annual vegetation productivity. A positive canopy growth response to longer NF seasons is expected to occur in cold-temperature-dominant landscapes, whereas in regions where summer moisture supply plays a greater role in limiting vegetation productivity the growth response will be reduced or reversed due to more frequent summer droughts and stronger evaporative demands. Potential moisture and temperature constraints underpinning the regional growth changes were evaluated using ancillary moisture-and temperature-constraint metrics derived from global reanalysis daily surface meteorology records. To verify our hypothesis we analysed a 29 year satellite microwave FT record to quantify NF season variability and trends over the HNL (≥45ºN) domain. We evaluated the ecological significance of these changes by comparing the FT record with an independent and well-calibrated satellite NDVI record in order to clarify the effects of FT-derived NF season variability on vegetation canopy growth for dominant cold temperature-and moisture-constrained areas within the HNL domain. We also compared the FT record with regional atmospheric CO 2 observations to specify how potential growing season changes associated with FT cycle variability are influencing HNL atmospheric CO 2 seasonal cycles.

Satellite microwave FT metrics
We used a global Earth System Data Record of daily landscape freeze-thaw status (FT-ESDR) derived from satellite microwave remote sensing to define FT metrics over all northern (≥45ºN) vegetated land areas . The FT-ESDR provides a daily (AM and PM) retrieval representing the predominant FR or NF status of the land surface within the satellite sensor field-of-view (approximately 25 km resolution). Overlapping 37 GHz, vertically polarized (V) brightness temperature (T b ) time series from the Scanning Multichannel Microwave Radiometer (SMMR) and Special Sensor Microwave Imager (SSM/I) have been integrated to produce a consistent and continuous long-term daily global FT data record currently extending from 1979 to 2010 ; the FT retrieval is derived on a grid-cell-wise basis using a temporal change classification of daily T b values in relation to FR and NF T b reference states using a seasonal threshold algorithm (STA) . Cross-sensor compatibility between the SMMR and SSM/I has been implemented through pixel-wise adjustment of the SMMR T b record based on empirical analyses of overlapping SMMR and SSM/I T b measurements (Kim et al. 2012). The FT-ESDR for this study is gridded to a global cylindrical Equal-Area Scalable Earth grid (EASE-Grid) projection with 25 km grid cell resolution (Brodzik and Knowles 2002). The FT-ESDR encompasses global vegetated land areas where seasonal FR temperatures significantly impact annual vegetation productivity . Four discrete FT classification levels are provided from daily composite conditions (CO): FR (AM and PM frozen), NF (AM and PM thawed), Transitional (AM FR and PM thawed), and Inverse Transitional (AM thawed and PM FR) statuses. The FT-ESDR showed mean annual FT spatial classification accuracies of 91.4 ±1.05 [inter-annual standard deviation (SD)] and 84.2 ±0.92 [inter-annual SD] percents for respective satellite ascending (PM) and descending (AM) overpasses in relation to daily surface air temperature records from the global WMO weather station network (Kim et al. 2012). Additional FT-ESDR quality assessment (QA) maps were developed for each annual data record, and provide more spatially explicit data quality information (Kim et al. , 2012. The FT-ESDR daily CO FT series was used in this investigation to assess HNL regional patterns, annual variability, and temporal trends in the timing and duration of the NF season, FR season, and transitional (TR) frost days for the 1979-2010 period. The NF and FR periods were derived from the daily CO FT record by summing the number of classified NF and FR days, respectively. Similarly, the number of TR days was determined as the sum of classified TR days for each grid cell within the domain. The day of primary seasonal thaw (T thaw ) for each cell was determined as the first day for which 12 out of 15 consecutive days (i.e. 80% rule) from January to June were classified as NF Kim et al. 2012;Bi et al. 2013).

Satellite NDVI
The NDVI from satellite (MODIS) optical-NIR remote sensing has been used to investigate the relative benefits of lengthening NF seasons for vegetation growth (Kim et al. 2012). For this investigation, bi-weekly 0.05º× 0.05º (5.6 km) resolution NDVI records provided by the University of Arizona Vegetation Index and Phenology lab (VIP) were used to identify HNL regional patterns and variability in mean summer (JJA) vegetation greenness anomalies over the longest available NDVI data record ; the VIP product represents a continuous, well-calibrated global NDVI record derived from multiple overlapping satellite optical-NIR sensors (Didan 2010). The VIP NDVI records for the 1982-1999 period were derived primarily from Advanced Very High Resolution Radiometer (AVHRR) daily surface reflectance measurements. For the period 2000-2010 the VIP record was derived primarily using MODIS (Moderate Resolution Imaging Spectroradiometer) spectral reflectances from the EOS Terra (starting from March 2000) and Aqua (starting from June 2002) satellites. Similar overlapping satellite spectral reflectance data from the SPOT4-VGT (vegetation instrument of the System Pour l'Observation de la Terre 4) sensor were used to bridge the 1999 AVHRR and MODIS records. The AVHRR NDVI was adjusted to the corresponding MODIS NDVI using spectral transformation equations for constructing longterm consistent NDVI records (Miura, Huete, and Yoshioka 2006). For this investigation we selected the VIP record because it represented the longest continuous global NDVI record available (1982-2010 period) with documented data quality information. Other satellite NDVI global data records are also available, including the Global Inventory Modeling and Mapping Studies (GIMMS; Tucker et al. 2005) product; however, regionally dissimilar trends and large discrepancies have been documented between GIMMS and other AVHRR NDVI records within the HNL domain (Alcaraz-Segura et al. 2010;Hall, Masek, and Collatz 2006;McCloy et al. 2005). The GIMMS NDVI record also only extended up to 2006 at the time of this investigation.
Cloudless and snow-free NDVI pixels determined from pixel-wise best data quality (QA) metrics were used for reprojection from the climate modelling grid (CMG) VIP product format to the 25 km global EASE-Grid of the FT-ESDR. We calculated the average summer NDVI conditions for each grid-cell (NDVI JJA ) by averaging June, July, and August NDVI values. The NDVI retrieval accuracy in spring is relatively more susceptible to degradation from artefacts, including snow cover and cloud effects, which are independent of canopy changes (Dye and Tucker 2003;Suzuki et al. 2011). We therefore used summer NDVI data as a surrogate measure of peak annual vegetation growth and to investigate how changes in FT status impact productivity in the HNL region as mediated by estimated cold temperature or moisture constraints.

Defining primary temperature-and moisture-constrained regions
A bio-climatological moisture stress index (MSI) and cold temperature stress index (TSI) were produced at 25 km spatial resolution for the 1982-2010 period and HNL domain following Zhang et al. (2008); the TSI and MSI were determined from surface daily minimum air temperature and vapour pressure deficit (VPD) parameters from global reanalysis data and used to define estimated environmental constraints to HNL annual vegetation net primary production (NPP). All daily surface meteorological variables used to define the MSI and TSI were derived from the National Centers for Environmental Prediction and National Center for Atmospheric Research (NNR; Kalnay et al. 1996;Kistler et al. 2001) global reanalysis following Zhang et al. (2008). Despite well-known deficiencies in the NNR product (Kanamitsu et al. 2002), NNR reanalysis data capture the major changes in surface climate anomalies (Simmons et al. 2004) and have been successfully used as meteorological inputs to estimate HNL and global vegetation productivity patterns and environmental trends (Nemani et al. 2003;Kimball, McDonald, and Zhao 2006;Zhang et al. 2008). To minimize influences of known errors in the NNR and other reanalysis data, the 29 year mean TSI and MSI values were used as climatological indicators of dominant cold temperature-and moisture-constrained grid cells, rather than attempting to represent more detailed MSI and TSI temporal trends. NNR reanalysis data at 1.9°(latitude) by 1.875°(longitude) were re-projected to the 25 km global EASE-Grid of the FT-ESDR. Daily VPD was used as a meteorological input for the MSI calculations rather than precipitation, because the atmosphere evaporative demand generally corresponds to soil wetness ) and large uncertainty in NNR precipitation has been documented (Drobot et al. 2006;Rawlins et al. 2006). Average daily air temperature and atmospheric vapour pressure were used to calculate daily VPD. The TSI was determined from minimum daily air temperature using a photosynthetic response curve that varies with plant functional type defined from a global land-cover classification ). An MSI value of 1 indicates effective photosynthesis cessation due to extreme moisture constraints to plant growth, whereas a value of 0 indicates that photosynthesis is not constrained by plant-available moisture limitations. The TSI defines the proportional loss of potential NPP due to cold temperatures, where the TSI ranges from 0 to 1 with increasing cold temperature constraints to NPP . We calculated annual summer mean TSI (TSI JJA ) and MSI (MSI JJA ) values for each grid cell within the HNL domain. These climate constraints were used as relative indicators of the moisture and temperature control factors influencing vegetation productivity. Dominant moisture (MSI)-and temperature (TSI)-constrained cells were determined from differences between MSI JJA and TSI JJA averaged over the 1982-2010 period for which VIP NDVI data are available; the resulting dominant moisture-and temperature-constrained cells were assigned as respective positive and negative values of these differences, respectively.

NOAA GLOBALVIEW-CO 2
We used the FT-ESDR as a surrogate indicator of FR temperature constraints to the HNL atmosphere seasonal CO 2 cycle. The FT metrics were compared against atmosphere CO 2 seasonal shape metrics defined from HNL monitoring stations represented by the NOAA ESRL GLOBALVIEW-CO 2 record (GLOBALVIEW-CO 2 2011). Because seasonal CO 2 cycles are strongly influenced by regional differences between ecosystem net carbon uptake (NPP) and carbon release from soil heterotrophic respiration, biomass burning, and anthropogenic emissions (Patra et al. 2005;Qian, Joseph, and Zeng 2010;Randerson et al. 1997;Sitch et al. 2007), it is necessary to address issues on the sparseness and temporal discontinuity in CO 2 observations when defining atmospheric CO 2 seasonal shape metrics (Vadrevu and Choi 2011;Xueref-Remy et al. 2010). Since the GLOBALVIEW-CO 2 record is derived from the integration of surface, tower, and aircraft measurements by an atmospheric transport model (Masarie and Tans 1995), this record is considered suitable for representing aggregated HNL atmospheric CO 2 variability (Higuchi et al. 2003;Murayama, Taguchi, and Higuchi 2004). The timing of atmospheric CO 2 drawdown by vegetation net photosynthesis in spring (C spr ) was computed from the weekly GLOBALVIEW-CO 2 record (>50°N) on an annual basis as the timing (day of year) of the downward 0-ppm crossing of the normalized (weeklyannual mean) atmospheric CO 2 seasonal record Randerson et al. 1999), and coinciding with the HNL growing season initiation Zhang et al. 2008). The annual minimum CO 2 concentration of the normalized seasonal record (C min ) was used as a relative measure of net photosynthetic carbon uptake Randerson et al. 1997). The magnitude of the normalized CO 2 concentration at the end of each calendar year (C end ) was used as an indicator of net ecosystem carbon release from soil heterotrophic respiration, biomass burning (Yi et al. 2009;Bond-Lamberty et al. 2007), and anthropogenic emissions (Le Quéré et al. 2009). In this study, the HNL aggregated mean annual NF season and T thaw anomalies from the FT-ESDR were compared against the corresponding atmosphere CO 2 seasonal metrics for the 1979-2010 period to investigate the potential impacts of annual variability in the NF season and primary spring thaw timing on the HNL atmospheric CO 2 seasonal cycle.

Statistical methods
Temporal trend analysis was used to quantify regional trends in NF season and NDVI summer growth changes over the HNL domain and for the 1982-2010 period. Monthly trends of HNL FT metrics derived from the longer (1979-2010) FT-ESDR record were used to examine seasonal changes of individual FT metrics and their relative contributions to the annual FT trends. The temporal trends in FT metrics and NDVI JJA growth anomalies were determined using a non-parametric Theil-Sen slope estimator approach, defined as the median of all pairs of slopes (Sen 1968). This approach is more accurate than simple linear regression when dealing with skewed and heteroscedastic data, and provides similar performance to the simple least squares method for normally distributed data (Akritas, Murphy, and LaValley 1995;Sen 1968). Pre-whitened series removed autocorrelation was used to compute the trends and their significance using the ZYP package in R statistics (Yue, Pilon, and Cavadias 2002) to mitigate autocorrelation effects and avoid trend test failure in the satellite time series (de Beurs and Henebry 2005). This method is insensitive to outliers, which are a challenge to all remote-sensing-based datasets (Alvera-Azcárate et al. 2012). Similar non-parametric approaches have been successfully used for environmental trend studies on a grid-cell-wise basis (Alcaraz-Segura et al. 2010;de Jong et al. 2011;Pouliot, Latifovic, and Olthof 2009).
We first calculated temporal anomalies of the FT metrics, NDVI JJA , and atmospheric CO 2 seasonal shape metrics with respect to average values defined from the period of record in case of having non-significant trends over the grid cells. If a significant (p ≤ 0.1) trend was determined, the temporal anomalies were derived as differences from the longterm detrended mean. When the trends and correlations were analysed, an outlier detection approach was applied to remove anomalous values with respect to the adjacent time sequence, and largely due to remaining optical-NIR cloud and aerosol contamination, and sensor discontinuity (Alvera-Azcárate et al. 2012). Outliers were identified on a grid-cellwise basis as values exceeding ±2 times the SD of the long-term mean (Moore 2006) and screened from the analysis. All statistical analyses were evaluated using a 90% (p ≤ 0.1) minimum threshold of significance.
We evaluated grid-cell-wise temporal correspondence between VIP NDVI JJA and NF season anomalies from the FT-ESDR record over the HNL domain. To evaluate the effects of NF season variability on HNL vegetation canopy growth within dominant temperature (TSI)-and moisture (MSI)-constrained areas we defined the NF period from January to August (NF JaAu ), coinciding with the period affecting NDVI JJA (Kim et al. 2012). The Pearson correlation coefficient (r) between NDVI JJA (independent variable) and NF JaAu (dependent variable) anomalies was used to assess the sign and strength of these relationships. Tundra, boreal forest, grassland, and temperate forest biomes within the HNL domain were categorized by a global terrestrial biome map (Olson et al. 2001), and the NDVI JJA and NF JaAu relationships were summarized for these individual biomes.

FT-ESDR QA
The FT-ESDR QA map (Figure 1(a)), averaged for the 1979-2010 period, provides a discrete, qualitative indicator of FT product quality for each grid cell within the HNL domain . The QA geographical map shows regions of relative high Figure 1. The mean annual FT-ESDR data QA map for the 1979-2010 record and HNL domain, aggregated by relative low (estimated spatial classification accuracy < 70%) to best (estimated accuracy > 90%) quality categories (a); the dimensionless QA values range from 0 to 1 for relative low to best quality categories and are responsive to landscape heterogeneity, with generally lower FT-ESDR data quality under higher fractional (% of grid cell) open water cover (b) and upper elevations (c).
to low product quality in relation to potential negative impacts on FT classification accuracy from temporal gaps in sensor T b data time series, active precipitation, open water, terrain and land-cover heterogeneity effects, and uncertainty associated with the use of global reanalysis temperature data to define per grid-cell FR and NF T b reference state thresholds for the STA-based FT classifications. Stepwise linear regression was used to define the independent variables influencing (R 2~0 .44) FT classification accuracy and the resulting QA values derived from the SMMR and SSM/I T b records. The QA values provide a dimensionless relative indicator of FT-ESDR data quality and range from 0 to 1 with increasing quality of the FT classification. The resulting HNL QA map ranges from low (estimated mean annual FT spatial classification accuracy < 70%) to moderate (70% < accuracy < 80%), good (80% < accuracy < 90%), and best (accuracy > 90%) quality categories. The mean proportions of the four QA categories encompass 30.6% (best), 59.5% (good), 9.2% (moderate), and 0.7% (low) of the HNL domain for the 1979-2010 FT-ESDR data record. The mean QA values are generally proportional to regional gradients in open water fraction and elevation over the domain (Figures 1(b) and (c)), with generally greater QA values at lower elevations and lower open water fractional coverage.

FT-ESDR trend analysis
The HNL FT-ESDR results show a strong latitudinal NF season gradient with shorter NF seasons at more northerly latitudes and a mean annual NF period of 177.6 ± 47.6 [Spatial SD] days over the HNL domain (Figure 2(a)). The FT results show a generally longer NF period and larger NF season variability along coastal margins relative to inland areas (Figures 2(a) and (b)). The FT record also shows a strong, positive HNL trend (2.4 days decade −1 ; p < 0.005) in the mean annual NF period that is largely driven by an NF season increase from March to June (Figure 3); the positive NF days' trend in spring is partially offset by a decreasing NF days' trend in the fall (SON), resulting in a smaller increasing annual NF season trend than might otherwise occur. A positive trend in the number of annual TR frost days (2.4 days decade −1 ; p < 0.1) results from a strong TR increase from October to November. A significant decrease in the mean annual FR season (−3.0 days decade −1 ; p < 0.001) is largely driven by a decreasing FR trend from March to May. An increase in the spring NF season is consistent with previous reports documenting warming trends in Arctic spring temperatures derived from model reanalysis (Screen, Deser, and Simmonds 2012). In fall, the decreasing NF days' trend coincides with a positive trend in the number of TR days, implying that the fall TR season is beginning earlier and extending over a longer period. The seasonal decrease in NF days and the corresponding increase in TR days imply regional cooling in the fall; the mechanisms for this apparent HNL cooling are uncertain but are consistent with reported increases in HNL moisture and snow cover in fall attributed to summer warming and sea ice decline (Cohen et al. 2012a(Cohen et al. , 2012bPark et al. 2013). The recent increase in fall snow cover may promote regional cooling due to an enhanced snow-albedo feedback (Déry and Brown 2007).
Mean annual trends of NDVI JJA and NF JaAu over the HNL domain are summarized in Table 1. The mean annual NDVI JJA records show general increasing trends of 5.6% Figure 3. The HNL trends in regional monthly mean FT metrics for the 1979-2010 record, including the NF season, transitional (TR) frost days, and FR period.

NF season effects on vegetation canopy growth
The regional pattern (Figure 4) of correlations between the NDVI JJA and NF JaAu anomalies shows generally positive correspondence above approximately 50º N (Figure 4(c)). There is a large drop in mean correlation at the highest latitudes due to the relatively lower quality (QA) of the FT classification and the smaller number of high-quality NDVI retrievals. Widespread negative correlations occur at lower latitudes over the HNL domain with a more heterogeneous pattern of positive and negative correlation areas. Areal proportions of positive and negative correlations between NF JaAu and NDVI JJA anomalies are summarized for the four HNL biomes in Table 2. For the 29 year (1982-2010) NDVI record, tundra and boreal biomes show greater regional proportions of positive correspondence between NF JaAu and NDVI JJA anomalies, whereas grassland and temperate forest biomes show similar areal proportions of positive and negative correlations. The NF JaAu and NDVI JJA anomalies show widespread positive correlations over 68.1% of the HNL domain, including NA (70.4%) and Eurasian (64.2%) portions of the domain. Positive correspondence indicates that years with a longer NF JaAu season relative to the long-term record coincide with greater summer vegetation canopy growth, whereas a negative correlation indicates that longer NF JaAu seasons promote less vegetation growth. The correlation analysis indicates that NF JaAu variability corresponds positively with NDVI JJA growth anomalies for a majority (>80%) of HNL boreal and tundra biomes, whereas correspondence is reduced or reversed for other areas.

Vegetation growth patterns within moisture-and temperature-constrained areas
The geographic distribution of estimated temperature (TSI) and moisture (MSI) control factors influencing vegetation activity over the HNL domain is presented in Figure 5. Dominant moisture-constraint cells are defined by positive difference (Diff MSI-TSI > 0) values, whereas negative differences indicate cells with dominant temperature constraint (Diff MSI-TSI < 0). The regional differences between MSI JJA and TSI JJA metrics averaged for the 1982-2010 period indicate that cold temperature constraints most strongly influence vegetation activity over 66.1% of the HNL domain and primarily within boreal and tundra biomes, whereas the temperature constraint to vegetation productivity is generally weaker in grassland and temperate forest regions, especially where plant-available moisture supply indicated by the MSI JJA is a greater limitation to annual productivity ( Figure 5). Areal proportions of dominant temperature-constrained areas over NA and Eurasia are 77.6% and 61.1%, respectively. Dominant temperature-constrained areas also occur along coastal margins, due to relatively more cloudiness and cooler growing seasons than other HNL areas. Dominant moisture-constrained areas are generally located in southern and inland continental regions of the HNL domain. These results are consistent with HNL continental interior areas showing vegetation browning associated with recent drought stress Buermann et al. 2013;Goetz et al. 2005;Zhang et al. 2008). Relatively lower quality FT classification results and smaller number of high-quality NDVI retrievals lead to a large drop in mean correlations at the highest latitudes. Error bars (grey shading) in (c) denote one SD range around the correlation means.
The regional pattern of NF JaAu and NDVI JJA correlations generally follows the pattern of estimated MSI and TSI dominance (Figure 6). A longer NF JaAu season generally enhances NDVI JJA productivity in dominant cold temperature-constrained areas, while  Table 2. Spatial extent of positive and negative correlations between NF JaAu and NDVI JJA anomalies expressed as a proportion (%) of the regional biomes and the larger HNL domain; Four r-value categories are classified by significance level, including significant positive (p ≤ 0.1), positive (p > 0.1), negative (p > 0.1), and significant negative (p ≤ 0.1) levels.

Sig. positive (%)
Positive ( the relative benefit of a longer NF season to vegetation growth is reduced or reversed in dominant moisture-constrained regions of the HNL domain. The sign of the NF JaAu and NDVI JJA relationship is generally positive and negative for respective cold temperatureand moisture-dominant areas, whereas NF JaAu and NDVI JJA correlation strength is proportional to the magnitude of TSI or MSI dominance. The relative strengths and regional patterns of the estimated temperature and moisture control factors may include substantial uncertainty associated with the coarse NNR surface meteorology used to derive the metrics (Buizza et al. 2005;Ma et al. 2008).

FT impacts on atmospheric CO 2 seasonal cycle
The HNL mean annual anomalies in FT-ESDR-defined T thaw corresponded significantly (p < 0.1) with C spr defined from GLOBALVIEW-CO 2 station observations over the 1979-2010 record (Figure 7). Temporal T thaw and C spr correspondence is predominantly positive (r = 0.409), indicating that relatively earlier (later) onset of the NF season coincides with earlier (later) growing season and C spr drawdown. The HNL T thaw anomaly for 1987 is approximately three-fold larger than other years of record and coincides with an extended period of TR frost days leading to a T thaw delay for this year. The T thaw anomaly also coincides with greater FT classification uncertainty and a lower annual FT-ESDR QA ranking for this year; the lower quality level was attributed to an extended gap in the satellite microwave T b record following the inauguration of the SSM/I record in 1987 (Kim et al. 2012). A relatively early C spr anomaly in 1990 is attributed to a shorter FR period in February and March relative to other years, and prior to the T thaw event, whereas an anomalous C spr delay for 1996 coincides with a longer FR period from January to April for this year. However, excluding outliers (1987, 1990, and 1996) from the correlation analysis resulted in similar T thaw and C spr correspondence (r = 0.407; p < 0.1). Correlations between HNL mean annual NF season variability and atmospheric CO 2 shape metrics anomalies for the 1979-2010 record were also analysed. We also examined these relationships over the shorter  record represented by the VIP NDVI record and found no difference in these relationships. Although the regional correlation between annual mean NF JaAu and NDVI JJA anomalies indicates a general direct benefit of longer NF seasons to vegetation productivity (r = 0.436, p < 0.05 without outliers; not shown), the spatial correlation pattern (Figure 4) implies that the relative benefit of a longer NF period to productivity and net ecosystem CO 2 uptake is spatially variable, leading to a weak HNL mean NF JaAu and C min relationship (r = -0.078, p > 0.1; not shown). A positive correlation between HNL mean annual NF season and C end anomalies (r = 0.708; p < 0.001) also accounts for the weak HNL NF and C min relationship because longer (shorter) HNL NF seasons promote elevated (lower) seasonal CO 2 maximums for the region. These results imply that longer NF seasons enhance both photosynthetic carbon (CO 2 ) uptake and terrestrial carbon emissions, resulting in a net annual increase in atmospheric CO 2 , whereas shorter NF seasons promote the opposite response. Potential sources of terrestrial carbon emissions associated with this response include soil heterotrophic respiration and disturbance, including wildfire and insect-related sources (Kurz et al. 2008;Randerson et al. 2006;Yi et al. 2013) and greater tree mortality Peng et al. 2011). However, the atmospheric CO 2 shape metrics may contain substantial uncertainty associated with strong annual CO 2 increase from fossil fuel emissions (Giorgi 2008).

Discussion and conclusions
The 1979-2010 FT-ESDR results show increasing HNL trends in the mean annual NF period (2.4 days decade −1 ; p < 0.005) and NF JaAu season (3.5 days decade −1 ; p < 0.001) that are approximately 26% larger than the corresponding mean NF season trends for the larger Northern Hemisphere domain (Kim et al. 2012). An increasing NF JaAu season trend occurs over 89% of the HNL domain, with 53% of the domain showing significant change relative to characteristic large annual NF season variability. These results Figure 7. Temporal correlation between mean spring drawdown of northern atmospheric CO 2 (C spr ) concentrations from the GLOBALVIEW CO 2 record and FT-ESDR-derived annual mean primary spring thaw timing (T thaw ) for the HNL domain.
generally coincide with the areal proportion of decreasing snow-cover duration estimated from Northern Hemisphere meteorological stations for the 1980-2006 record (Peng et al. 2013). Our results indicate that each day of increase in the NF JaAu season produces a relatively strong NDVI JJA increase (0.15 ± 0.2%) in dominant cold temperature-constrained areas, whereas the same NF increase produces a much smaller and insignificant growth increase (0.06 ± 0.1%) in dominant moisture-constrained HNL areas. A given daily increase in the NF season also produces a larger NDVI JJA increase in tundra (0.25 ± 0.2%) and boreal biomes (0.13 ± 0.2%), and relatively weak growth response in grassland (0.05 ± 0.1%) and temperate forest areas (0.06 ± 0.2%). Our FT and CO 2 correlation results confirm previous reports documenting a strong link between the timing of the primary spring thaw event defined from satellite microwave remote sensing and the onset of the vegetation growing season in northern biomes (Goetz et al. 2007;Kimball et al. 2004;McDonald et al. 2004). Our findings also support previous findings that a lengthening NF season coupled with increasing disturbance is likely to promote respiration and carbon emissions over productivity, diminishing HNL carbon sink strength (Angert et al. 2005;Yi et al. 2013).
The present geographic distribution of dominant moisture-and temperature-constrained regions indicates that a lengthening NF season is enhancing vegetation growth and productivity for a majority (82.8%) of HNL tundra and boreal biomes. These results are consistent with previous studies documenting satellite-derived HNL greening trends (Bhatt et al. 2010;Goetz et al. 2005;Xu et al. 2013). However, a longer NF season promotes reduced vegetation productivity benefits or adverse growth conditions for 51.9% of HNL grasslands and temperate forests due to summer water supply (MSI) restrictions to plant growth. Strengthening evaporative demands and more frequent summer droughts associated with earlier snowmelt and warmer summer temperatures may lead to more severe negative impacts on vegetation growth in more moisture-constrained areas (Barnett, Adam, and Lettenmaier 2005;Jepsen et al. 2012). These results are similar to previous studies indicating widespread drought-induced vegetation productivity declines in northern temperate and boreal forests Piao et al. 2011;Ma et al. 2012).
The positive correspondence between HNL mean T thaw and C spr anomalies indicates that earlier spring thawing promotes generally earlier ecosystem carbon (CO 2 ) uptake. Abundant spring solar radiation coupled with earlier seasonal thawing promotes deeper soil active layer development and release of plant-available moisture and nutrients, leading to rapid vegetation carbon uptake and earlier C spr (Bergh and Linder 1999;Tanja et al. 2003;Wang et al. 2011). However, our results indicate that longer NF seasons also promote increased amplitude of the atmospheric CO 2 seasonal cycle; thus the relative increase in CO 2 uptake by ecosystem productivity is being offset by enhanced carbon emissions, likely due to commensurate increases in regional disturbance and respiration processes (Angert et al. 2005;Buermann et al. 2013;Kasischke et al. 2005).
The results of this study document a consistent increase in the annual NF season over northern land areas from 1979 to 2010 as derived from the satellite passive microwave remote-sensing record. The lengthening NF season is consistent with global warming and appears to be largely benefitting HNL vegetation productivity, inferred from satellite NDVI summer growth changes, by promoting earlier and longer growing seasons. However, the relative benefits for vegetation productivity do not appear to be increasing terrestrial carbon sink activity indicated from the regional atmospheric CO 2 record, as carbon emissions from enhanced disturbance and respiration processes appear to be offsetting CO 2 uptake from additional vegetation growth. Our results also indicate that the relative benefits of a lengthening NF season for vegetation growth are spatially heterogeneous and are reduced or reversed in dominant moisture (MSI)-constrained areas, including grassland and northern temperate forest areas. A lengthening NF season is expected to continue to benefit productivity in dominant cold temperature (TSI)constrained regions, including tundra, although these constraints will likely relax with associated contraction of dominant TSI-constrained areas under continued HNL warming. In contrast, dominant moisture (MSI)-constrained areas are expected to become more widespread, with more persistent and severe summer drought constraints (Goetz et al. 2005;Ma et al. 2012). These changes are likely to promote a net decrease in HNL vegetation productivity and the terrestrial carbon sink without commensurate increases in regional precipitation and terrestrial water storage required to offset expected increases in evaporative demands under a warming climate (Hayes et al. 2011;Parmentier et al. 2011;Schuur et al. 2008).

Funding
This work was conducted at the University of Montana under contract to the National Aeronautics and Space Administration, with funding support provided by the NASA Terrestrial Ecology and Hydrology programmes. The FT-ESDR and VIP NDVI data records used for this investigation were provided by the National Snow and Ice Data Center (NSIDC) and the University of Arizona, with funding support provided by the NASA measures (Making Earth System Data Records for Use in Research Environments) programme. The GLOBALVIEW-CO 2 data record used for this study was provided by the NOAA ESRL (Earth System Research Laboratory).