Is arctic greening consistent with the ecology of tundra? Lessons from an ecologically informed mass balance model

Climate change has been implicated in the widespread ‘greening’ of the arctic in recent decades. However, differences in arctic greening patterns among satellite platforms and recent reports of decreased rate of greening or of browning have made attributing arctic greening trends to a warming climate challenging. Here, we compared MODIS greening trends to those predicted by the coupled carbon and nitrogen model (CCaN); a mass balance carbon and nitrogen model that was driven by MODIS surface temperature and climate. CCaN was parameterized using model-data fusion, where model predictions were ecologically constrained with historical ecological ground and satellite based data. We found that, at long temporal and large spatial scales, MODIS greening trends were consistent with ecological and biogeochemical data from arctic tundra. However, at smaller spatial scales, observations and CCaN greening trends differed in the location, extent, and magnitude of greening. CCaN was unable to capture the high rates of MODIS greening in northern wetlands, and the patchy MODIS browning in the southern portion of the North Slope. This model-data disagreement points to disturbance and its legacy impacts on land cover as an important mechanism for understanding greening trends on the North Slope of Alaska.


Introduction
Satellite based observations of the normalized difference vegetation index (NDVI) have illustrated a 'greening of the arctic'. Because NDVI is a proxy for leaf area, several have hypothesized that these greening trends result from increased vegetation productivity as a result of the rapidly warming temperatures (Walker et al 2003, Goetz et al 2005, Beck and Goetz 2011, Forkel et al 2016. However, attributing greening trends to increased vegetation productivity is not straightforward (Alcaraz-Segura et al 2010). Satellites differ in the magnitude and spatial extent of greening, while some regions have exhibited insignificant greening trends or declines (i.e. 'browning') even as the arctic continues to warm (Stow et al 2007, Bhatt et al 2013, Ju and Masek 2016. Moreover, changes in background reflectance or satellite performance (e.g. orbital or sensor drift) may confound such trends resulting in a spurious greening signal without an increase in leaf area index (LAI) (Rocha and Shaver 2009, De Jong et al 2011, Wang et al 2012. Hence, in order to better understand the observed greening response, it must be determined whether the greening signal is expected given a priori knowledge of the ecophysiology and biogeography of arctic vegetation. Due to low diversity and low temperatures, panarctic CO 2 fluxes are thought to be functionally convergent and can be reasonably modeled with a single response to light, temperature, and canopy leaf area (Street et al 2006;Shaver et al 2007Shaver et al , 2013Loranty et al 2011;Street et al 2012;Mbufong et al 2014). This functional convergence is thought to arise from strong nitrogen (N) limitation of arctic vegetation, which controls both allocation to LAI and leaf photosynthetic rates Rastetter 1999, Shaver et al 2013). Low mineralization, deposition, and leaching rates cause a tight recycling of N among biomass, soil organic matter, and available N pools (Shaver et al 1992, Stieglitz et al 2006. Consequently, arctic greening is dependent on the tight coupling between C and N cycles. This conceptual framework was implemented into a parsimonious mass balance model of coupled C and N cycles (e.g. the coupled carbon and nitrogen (CCaN) model) (Wright and Rocha 2018). CCaN was parameterized with model-data fusion, where prior observed ecological information (i.e. eddy covariance, NDVI, and vegetation and soil data from the Arctic Long Term Ecological Research (LTER) program) provided realistic biogeochemical constraints on model parameters and predictions. CCaN performed well at capturing the seasonal and spatial variability in NDVI across arctic systems on the North Slope of Alaska (Wright and Rocha 2018). Consequently, it was hypothesized that CCaN along with a model-data fusion approach would provide insight into arctic greening trends.
The modern landscape temperature sensitivity of arctic productivity may place a strong constraint on arctic greening trends. Temperature and NDVI on the North Slope of Alaska are tightly linked and follow a general north to south gradient of increasing temperature and NDVI (Raynolds et al 2008;Wright and Rocha 2018). Such sensitivity is likely related to the long-term equilibrium response of arctic biogeochemistry to temperature; providing an upper limit to the expected shifts in ecosystem productivity with warmer temperatures. Here model-data fusion was used with CCaN and latitudinal MODIS NDVI observations on the North Slope to constrain the temperature sensitivity of LAI under tightly coupled C and N cycles. Model-data fusion is an iterative technique used to find the best set of parameters that fit a proposed model to data with quantified uncertainties (Keenan et al 2012). Parameters were constrained within a given range, given a priori information from historical LTER measurements (Wright and Rocha 2018). Model-data fusion provided a posteriori statistical likelihood parameter distribution based on model structure, data uncertainty, and parameter constraints. We hypothesize that such constraints will allow for attribution of the observed greening trends from 2000-2015 to increased temperatures in the region.

Region of interest
Analyses focused on the foothills and coastal plains of the North Slope of Alaska, but excluded the Brooks Range mountains (figure 1). The region of interest encompassed 151 427 km 2 and included a diverse array of topography, landscape age, soil types, and ecosystem communities. Elevation generally decreases with latitude from ∼1000 m above sea level (asl) in the South to <5 m asl in the North. Soils in the region vary widely in age from 14 kya to >900 kya and pH from <5.5 to >6.5 as a result of past glaciation events and climate (Walker et al 1998, Raynolds et al 2006. The region comprises nine out of fifteen arctic vegetation community types as defined in Raynolds et al (2006), and included erect-shrub tundra in the south, acidic and non-acidic graminoid tundra in the foothills of the North Slope, and wetlands in the north. The North Slope also is a region where satellite based NDVI trends have diverged across satellite platforms, indicating a need for further validation of arctic greening trends (Stow et al 2007, Ju andMasek 2016).

CCaN model and parameterization
Details on CCaN formulation and parameterization are described in Wright and Rocha (2018). Briefly, CCaN couples C and N cycles through the linear relationship between LAI and foliar N (Williams and Rastetter, 1999), a shared litter-fall rate for C and N, and plant N uptake limited by root biomass. The model is driven with daily light and temperature, and predicted daily net ecosystem exchange of CO 2 (NEE), LAI, and NDVI using the functional convergence model of Shaver et al (2007Shaver et al ( , 2013 and Loranty et al (2011). NDVI was obtained by inverting the exponential relationship between ground based LAI and NDVI observations (Shaver et al 2007, Rocha andShaver 2009). CCaN was constrained with field measured C and N cycling processes and their uncertainties, and was parameterized with model-data fusion (Richardson et al 2010, Keenan et al 2012).
The functional convergence model of Shaver et al (2007Shaver et al ( , 2013 and Loranty et al (2011) hypothesizes that NEE in all arctic tundra communities responds similarly to temperature, light, LAI, and NDVI throughout time. In accordance with this hypothesis, CCaN was originally parameterized with NEE and NDVI observations from a single site (Wright and Rocha 2018). With this parameterization, CCaN captured a majority of NDVI observations over 8 years at a single site and NDVI observations across a latitudinal gradient. The gradient was comprised of 32 tundra sites along a North-South transect that included all major vegetation communities found on the North Slope of Alaska (Wright and Rocha 2018). However, this original parameterization underestimated NDVI at sites with high LAI (Wright and Rocha 2018). Consequently, for this study we reparameterized CCaN using the MODIS NDVI from the 32 tundra sites previously mentioned, with a focus on constraining the temperature sensitivity of the proportion of N in leaves, which ultimately determined LAI and NDVI.
CCaN was reparameterized using a weighted leastsquares model-data fusion technique based on the approaches of Keenan et al (2012) and Richardson et al (2010). To determine starting values for CCaN pools, a model-spin up was run for thirty years using decadal averaged daily forcing data. This was done to ensure that the NDVI trends were not a result of the initial model conditions being out of equilibrium with the current climate. The pool starting values for the model spin-up were determined from previous publications (Jiang et al 2016). The model-data fusion routine was an adaptive Markov-chain Monte Carlo that used Metropolis-Hastings sampling (Metropolis et al 1953) and incorporated measurement uncertainty into the parameterization, allowing propagation of parameter uncertainty into model predictions. In the first step of the optimization, an individual cost function (ji) was calculated as the total uncertainty-weighted squared data-model disagreement for each data type (Richard- where N i was the number of observations for data type i, y i (t) was the data at time t, pi (t) was the value predicted by the model, and δ i (t) was the measurement uncertainty. NDVI uncertainty was quantified as a constant 0.07 based on analyses of ground and MODIS NDVI observations (Wright and Rocha, 2018). In the initial parameterization, data on soils and vegetation from the Toolik LTER were used to assign biologically valid maximum and minimum parameter values (Wright and Rocha 2018). For the reparameterization of the temperature sensitivity of the proportion of N in leaves, parameter values were randomly drawn from a normal distribution with a mean equal to the initial parameter (Wright and Rocha 2018) and standard deviation equal to a fraction of the initial parameter range that was adjusted to achieve an acceptance rate of 25%-30%, preventing the routine from getting stuck at local minima ( Step two was run until 1000 parameter sets were accepted, providing distributions and uncertainty estimates for each parameter.

CCaN model drivers and MODIS NDVI data
We obtained daily incoming shortwave radiation for the period of 2000-2015 from the North American Regional Reanalysis (NARR) for Arctic Research (http://www.esrl.noaa.gov/psd/data/narr/). The spatial resolution of NARR grid data was 0.3°by 0.3°. We downscaled the data to 1×1 km spatial resolution using bilinear interpolation and extracted the spatial domain of the North Slope as in Jiang et al (2016). Bilinear interpolation is a resampling method that uses the distance weighted average of the surrounding coarse scale pixels to estimate the value of the newly finer downscaled resolution (Zhang et al 2017). Shortwave radiation was multiplied by 0.1686 to translate the unit from W m −2 to mol PAR m −2 d −1 (Papaioannou et al 1993). We used 15 years (2000-2015) of 1 km, 8 d MODIS Land Surface Temperature and Emissivity (MOD11A2) to characterize land surface temperature (Wan et al 2004). We used surface temperature instead of air temperature for several reasons. First, permafrost creates a strong heat sink, and as a result, surface and atmospheric temperature rarely deviate from each other during the growing season in tundra ecosystems  (Eugster et al 2000, Rocha and. Second, surface temperature better reflects the thermal environment experienced by vegetation (Eugster et al 2000). Third, surface temperature can be reliably measured at large spatial resolution with the MODIS satellite, reducing potential errors from extrapolation of sparse meteorological measurements. We removed pixels with poor data quality and gap-filled with a modified Savitsky-Golay filter (Chen et al 2004). The Savitsky-Golay filter also was used to transform the 8 d surface temperature product into daily surface temperature product for the CCaN model. MODIS surface temperature was retrieved from the Land Processes Distributed Active Archive Center (LP-DAAC) using the R package 'ModisDownload' (r-gis.net).
We used 1 km, 8 d MODIS NDVI (MOD13A2) to characterize vegetation greenness across the North Slope of Alaska from 2000-2015. MODIS T s and NDVI products were downloaded from LP-DAAC using the R package 'ModisDownload' (r-gis.net) for the North Slope region. MODIS NDVI fill-value (poor quality) pixels from cloud or snow cover were removed and gap-filled with a modified Savitsky-Golay filter (Chen et al 2004). For CCaN parameterization, we extracted a 15 year NDVI time series from 32 sites (2×2 km in area) along a North-South transect extending from 65-70°N that ran along Alaska's Dalton Highway (Wright and Rocha 2018). This transect included all vegetation communities and soil types found on the North Slope of Alaska and ranged from erect-shrub tundra in the south, acidic and non-acidic graminoid tundra in the foothills of the North Slope, and wetlands in the far north. The area encompassed by the 32 sites represented less than 0.1% of the North Slope area. These data were used for CCaN parameterization and determination of MODIS NDVI uncertainty (Wright and Rocha 2018). Lastly, we aggregated MODIS NDVI into growing season (June-August) averages and then further aggregated to a 15 year average growing season NDVI map for the North Slope. This 15 year average growing season MODIS NDVI map was used to validate the CCaN NDVI map for the same period using the same growing season average.

Arctic vegetation map and water mask
Tundra vegetation was delineated using the Circumpolar Arctic Vegetation Map (CAVM). The CAVM includes the spatial distribution of 85 plant community descriptions across Alaska at a spatial resolution of ∼14 km (Raynolds et al 2006). Vegetation types on the North Slope were grouped according to the following main categories: Barrens (Class: B3, B4), Graminoid tundra (Class: G3, G4), Shrub tundra (Class: S1, S2), and Wetlands (Class: W1, W2, W3) (figure 1). Greening trends observed in MODIS and CCaN NDVI were compared using the difference between the average MODIS and CCaN greening trend for each vegetation category. A positive difference indicated that MODIS produced larger greening trends, while a negative difference indicated the CCaN produced larger greening trends.
To eliminate low NDVI values due to open water and exposed floodplain sediment, we removed MODIS and CCaN pixels that were located on or near water bodies. Analyses were carried out both with and without the water mask in place. Waterbody vector data for the North Slope was obtained from the National Hydrology Dataset (http://catalog.northslope.org/ or https://nhd.usgs.gov/). These data were used to locate large rivers (>100 miles in length) and lakes (>1 km 2 in area) that were masked out of subsequent NDVI and NDVI trend analyses with a 250-500 m buffer on each side.

Statistical analyses
NDVI trends (i.e. greening or browning trends) were determined with a modified Mann Kendall and a Theil Sen robust linear regression at the 95% confidence level. The modified Mann Kendall is a non-parametric test for a monotonic trend that also takes account of autocorrelation in the data (Hamed and Ramachandra Rao 1998). The test only determines the presence of a trend in the data without calculating a regression slope. The regression slope (i.e. the greening trend) was determined with a Theil Sen robust linear regression, which is also a non-parametric test that calculates the median of the slopes of all lines through pairs of points (Gilbert 1987). Non-parametric tests were chosen over traditional least squares regression methods to account for skewed or heteroscedastic data and outliers that may influence the estimation of trends (Christensen 2001). Greening trends that were below the 95% confidence level were not included in subsequent analyses.
Differences in MODIS and CCaN greening trends for each vegetation type were determined with a student's t-test. Both MODIS and CCaN data demonstrated a high degree of spatial autocorrelation, which violates the assumption of the t-test that values should be independent of one another. This spatial autocorrelation inflates the degrees of freedom of the statistical test and increases the likelihood of type II statistical errors (Christensen 2001, Zhang et al 2017. Consequently, the degrees of freedom was determined by the number of independent greening clusters in the MODIS and CCaN datasets using a density-based clustering analysis (self adjusting HDBSCAN method) in ArcGIS pro. The density-based clustering analysis identified spatially correlated groups of greening trends and delineated a number of unique clusters across the spatial domain. This removed the spatial autocorrelation in the data and allowed for independence among the samples. We delineated greening trend clusters for each vegetation type and used the number of clusters in each sample as the degrees of freedom. Statistical significance was determined at the 95% confidence level.

Results
CCaN captured the broad scale spatial patterns of average growing season NDVI from 2000-2015 (figure 2). Both MODIS and CCaN NDVI generally declined with latitude with distinct latitudinal bands across the North Slope (figures 2(A), (B)). MODIS NDVI latitudinal differences were more noticeable than those predicted by CCaN, especially in the western region of the North Slope with a spatial coefficient of variation for NDVI that was double that of CCaN (figure 2(C)). CCaN generally under predicted NDVI in the south western central region, and overestimated NDVI in the alluvial dominated Northeast region of the North Slope (figure 2(C)). MODIS NDVI had larger NDVI range than CCaN across the North Slope, ranging from 0.24 to 0.71 for MODIS and from 0.34 to 0.63 for CCaN ( figure 2(D)). The average latitudinal gradient of MODIS NDVI was well correlated with CCaN with an R 2 of 0.82 ( figure 2(D)). This compares an R 2 of 0.65 for the MODIS and CCaN NDVI data from the subset of 32 sites used to inform CCaN parameterization.
Both CCaN and MODIS indicated widespread NDVI changes across the North Slope from 2000-2015 (figure 3). CCaN indicated more significant NDVI trends than MODIS. MODIS indicated a statistically significant NDVI trend for 35% of the North Slope, whereas CCaN indicated a statistically significant NDVI trend for 71% of the North Slope (figures 3(A), (B)). A majority of the NDVI trends were positive, indicating greening, but 0.01% of MODIS and 0.001% of CCaN NDVI trends were negative, indicating browning. MODIS exhibited larger variation in greening trends than CCaN, ranging from 0.003 to −0.001 ΔNDVI y −1 for MODIS and 0.001 to −0.001 ΔNDVI y −1 for CCaN.
There were noticeable differences in the spatial extent of greening and browning trends for MODIS and CCaN (figure 3). At very fine scales (1 km resolution), spatial greening patterns were difficult to capture by the model with ∼30% of trends overlapping for both model and observation. CCaN browning NDVI trends were localized and formed several distinct clumps on the northern region of the North Slope. MODIS browning NDVI trends were scattered with a majority of patches occurring in the south central region of the North Slope. Strong MODIS greening trends occurred along the Northern coastal plains, whereas strong CCaN greening trends occurred along the southwestern region of the North Slope. This resulted in high positive model disagreement in the North, indicating larger MODIS NDVI trends, and high negative model disagreement in the Southwest, indicating larger CCaN NDVI trends ( figure 3(C)). NDVI greening trends were statistically similar for both CCaN and MODIS across the entire region (pvalue: 0.11) ( figure 3(D)). Although not statistically significant, MODIS produced a slightly larger average greening trend (0.004 ΔNDVI y −1 ) than CCaN (0.003 ΔNDVI y − ΔNDVI y −1 ) ( figure 3(D)).
There was significant variation in CCaN and MODIS greening trends with latitude and vegetation communities ( figure 4). At low latitudes and along the foothills of the North Slope (<69°N), averaged CCaN and MODIS greening trends did not significantly differ (p-value: 0.7). Differences between CCaN and MODIS greening trends were more pronounced at higher latitudes ( figure 4(A)). In general, greening MODIS and CCaN trend differences increased with increasing latitude with substantially higher greening trends reported for MODIS than CCaN. On average, greening trends were twice as large for MODIS than CCaN above 70°N. This difference was attributed to Wetland communities where MODIS trends were 44% higher than predicted by CCaN ( figure 4(B)). Differences in greening trends in wetland communities were significantly different (p-value <0.01), whereas differences in greening trends among other vegetation communities were non-significant

Discussion
Our analyses revealed that greening on the North Slope of Alaska from 2000-2015 was consistent with the landscape temperature sensitivity of arctic NDVI/ LAI. In general, increased temperatures increased NDVI/LAI, which is consistent with experimental warming manipulations with greenhouses and opentop chambers (Shaver et al 1992, Rusted et al 2001, Elmendorf et al 2011, Shaver et al 2014, Buruah et al 2017. Our approach of substituting the spatial temperature sensitivity of NDVI with the temporal sensitivity of NDVI was similar to Raynolds et al (2008).  ]. Error bars are 95% confidence intervals and stars indicate significant difference from zero at the 95% confidence level. Raynolds et al (2008) used a qualitative statistical model to determine the temperature response of NDVI, which precluded a mechanistic understanding of NDVI trends. However, our approach was mechanistic allowing for the prediction and attribution of NDVI trends to temperature by using a mass balance model that was constrained with uncertainty assessed field observations. Consequently, the model-data fusion approach tested our ecological knowledge of tundra systems and identified significant model-data mismatches (Rastetter 2017, Wright andRocha 2018) where NDVI is influenced by factors other than temperature.
It is clear that, at large spatial and long temporal scales, Arctic NDVI/LAI responses to temperature are tightly constrained. This is supported by CCaN's ability to capture the average latitudinal NDVI gradient and the range of greening trends with a single parameterization with information from <1% of the total region. These observations support a functionally convergent response of coupled carbon and nitrogen cycles to temperature and radiation across the arctic (Loranty et al 2011, Wright andRocha 2018). It was also surprising that CCaN was able to capture broad scale NDVI patterns without the inclusion of important arctic processes that control C and N cycling such as thaw depth and soil moisture. However, surface temperature, thaw depth and soil moisture covary through the energy budget, and these processes are likely represented in CCaN at longer decadal timescales and larger spatial scales through this covariation Rastetter 1999, Shaver et al 2013).
Although CCaN captured the broad scale arctic NDVI patterns, variability at finer spatial scales was more challenging to replicate. Some of the larger differences between MODIS and CCaN NDVI also coincided with spatial variability in pH and underlying bedrock. The large difference in NDVI in the northeast coincides with a patch of younger tertiary bedrock (Wilson et al 2015). A majority of the North Slope is comprised of older Mesozoic bedrock, and landscape age is a strong determinant of soil pH in this region (Tieszen 1978;Raynolds et al 2006). Soil pH is not included in CCaN, but is important because it influences arctic species diversity, nutrient availability, vegetation productivity, and land surface energy exchange (Walker et al 1998, Gough et al 2000, Walker et al 2003, Walker et al 2012. These large mismatches between model and observations highlight the importance of understanding soil impacts on vegetation and NDVI across the arctic.
On average, MODIS and CCaN NDVI trends agreed until 70°N and were similar for Barrens, Graminoid, and Shrub communities. However, there were areas where CCaN indicated either a greening or browning without congruent trends in MODIS. Since CCaN was largely driven by temperature, this lack of congruency indicated areas of warming, but no expected positive change in NDVI. This lack of tundra NDVI response to temperature has been shown elsewhere and supports a more nuanced explanation of attributing greening to temperature alone (Raynolds et al 2013, Lara et al 2018. It is likely that the mismatch among greening, no significant change, and browning patches across the North Slope highlights the importance of process other than temperature in determining NDVI trends. These processes are obviously occurring at much smaller scales than the large spatial surface temperature increases on the North Slope imply, due to the patchier distribution of MODIS NDVI trends as compared with CCaN. MODIS tundra browning was much patchier with a higher spatial standard deviation than predicted by CCaN. A majority of these small scale browning patches occurred in southern graminoid communities, and could be the result of several factors including increased litter production through either increased plant productivity, changes in vegetation composition (i.e. shrubification), or increased plant mortality. Litter decreases NDVI by brightening the surface and decreasing the amount of green leaf exposure through shading (Van Leeuwen, Huete 1996, Rocha et al 2008. Litter production in graminoid communities is high and can be the result of both proposed mechanisms. Tundra greening trends have decreased or browned in recent years even as the climate continues to warm (Bhatt et al 2013, Pheonix and Bjerke 2016, Vickers et al 2016. Early snow melt and frosting events also have been identified as an important process in tundra browning (Pheonix and Bjerke 2016). The ecological mechanisms regulating tundra browning remain elusive. However, our analyses highlights that these processes occur independently of large scale temperature increases.
CCaN assumed functional convergence in ecosystem cycles of C and N, which indicated that the parameters driving flux rates in all arctic community types are constant once light, temperature, and LAI are taken into account. In this model, all arctic vegetation communities function the same, and hence, individual parameterizations based on vegetation community type are unnecessary allowing for pan-arctic prediction of ecosystem function. CCaN also used the modern latitudinal temperature sensitivity of LAI and NDVI to provide an upper bound to their temporal temperature sensitivity. Such assumptions simplified the analyses and highlighted differences in the functioning among arctic vegetation community types. On average, most vegetation communities followed functional convergence in C and N fluxes provided the modern day latitudinal temperature sensitivity of NDVI and LAI, with the exception of northern wetland communities. These communities had higher temperature sensitivities than other community types on the North Slope because of the larger than expected (i.e. by the latitudinal NDVI temperature response in CCaN) northern wetland MODIS greening trends (figures 3 and 4). This suggests that either northern wetlands do not follow the functional convergence hypothesis and the intrinsic latitudinal temperature sensitivity of NDVI, which would be inconsistent with Shaver et al (2007Shaver et al ( , 2013 and Raynolds et al (2008), or points to a mechanism that allows these systems to more rapidly respond to temperature increases than anticipated from the latitudinal NDVI-temperature gradient.
Northern wetland areas have been subjected to large changes in hydrology as a result of permafrost thaw, resulting in decreased lake cover and surface moisture (Liljedahl et al 2016). These changes create disturbance and initiate succession of newly colonized vegetation that stimulates vegetation productivity (Odum 1969). Such dynamics accelerate temporal changes in NDVI through the initiation of succession and the stimulation of vegetation productivity during mid-successional stages (McMillan andGoulden 2008, Rocha et al 2012). The patchiness of greening trends also support a disturbance induced greening as disturbances tend to be smaller in scale than regional temperature increases. Successional dynamics also may explain why arctic greening has recently slowed even as temperatures continue to increase. Older disturbed patches will enter a period of senescence due to increased nutrient limitation (Rocha et al 2012). Permafrost thaw, subsidence, and thermokarst events are increasing in the arctic and impart successional land cover legacies that can decouple vegetation from year to year temperature change (Raynolds et al 2008, Liljedahl et al 2016, Pastick et al 2018. Other disturbances could include increased productivity and shrub cover in water tracks, which are likely more common in southern parts of the North Slope (Curasi et al 2016). The role of disturbance in explaining greening trends has gained recent attention (Lara et al 2018), and our results corroborate the importance of land cover legacies in understanding arctic greening trends.

Conclusions
Understanding links between climate change, arctic greening, and C cycling are critical in constraining the temperature sensitivity of arctic vegetation productivity and predicting future arctic C cycling (Schuur et al 2015, Abbott et al 2016, Fisher et al 2018. Because of its large spatial extent, remote sensing provides a valuable tool in understanding these links across the arctic, but remains limited because it includes satellite measurement artifacts and provides no mechanistic understanding due to the indirect measure of vegetation productivity through a spectral reflectance proxy. Consequently, there is a need to validate and verify vegetation changes measured with remote sensing with ground based observations. Here we demonstrate the utility of combining a parsimonious mass balance model with historical ground based measurements of arctic tundra function to develop a mechanistic understanding of arctic greening. We demonstrate that northern wetlands are greening at a much higher rate than expected, and that many areas on the North Slope are warming but not greening. We hypothesize that these patterns indicate that disturbance and successional processes are important for interpreting greening trends. Further testing of this hypothesis using such an approach would be useful for constraining the response of vegetation to warming and improving predictions of the arctic C cycle in a future warmer world.

Acknowledgments
This work was supported by a NASA Earth and Space Science Fellowship to Kelseyann Wright, NSF grants #1065587 and #1026843 to the Marine Biological Laboratory, and NSF grant #1556772 to the University of Notre Dame. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the US Government.