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

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, pan-arctic CO₂ 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 2007, 2013; Loranty 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 (Williams and 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.

2.1. 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² 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 and Masek 2016).

Zoom In Zoom Out Reset image size

2.2. 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₂ (NEE), LAI, and NDVI using the functional convergence model of Shaver et al (2007, 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 and Shaver 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 (2007, 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 least-squares 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 (Richardson et al 2010) [equation (1)].

$$\begin{eqnarray}&&{j}_{i}=\,\left(\displaystyle \sum _{t-1}^{{N}_{i}}{\left(\displaystyle \frac{{y}_{i}(t)-{p}_{i}(t)}{{\delta }_{i}(t)}\right)}^{2}\right),\end{eqnarray} \tag{ 1 }$$

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 (Richardson et al 2010, Keenan et al 2012). The best parameter value was determined by minimizing the overall cost, which was the product of the individual cost functions (Richardson et al 2010). In step two, uncertainty-weighted data-model disagreements were normalized based on the variance of the j's obtained in step one, and parameter sets were only accepted if the cost function for each data stream (ji) passed a chi-square test at 90% confidence. (Franks et al 1999, Richardson et al 2010, Keenan et al 2012). Step two was run until 1000 parameter sets were accepted, providing distributions and uncertainty estimates for each parameter.

2.3. CCaN model drivers and MODIS NDVI data

2.4. Arctic vegetation map and water mask

2.5. Statistical analyses

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² of 0.82 (figure 2(D)). This compares an R² of 0.65 for the MODIS and CCaN NDVI data from the subset of 32 sites used to inform CCaN parameterization.

Zoom In Zoom Out Reset image size

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⁻¹ for MODIS and 0.001 to −0.001 ΔNDVI y⁻¹ for CCaN.

Zoom In Zoom Out Reset image size

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 (p-value: 0.11) (figure 3(D)). Although not statistically significant, MODIS produced a slightly larger average greening trend (0.004 ΔNDVI y⁻¹) than CCaN (0.003 ΔNDVI y⁻ΔNDVI y⁻¹) (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 (Barrens [p-value: 0.92]; Gramminoid [p-value: 0.13]; Shrubs [p-value: 0.7]).

Zoom In Zoom Out Reset image size

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 open-top 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). 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 and Rocha 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 and Rocha 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 (Williams and 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 (2007, 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 and Goulden 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.

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.

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.