Spatial variation in vegetation productivity trends, fire disturbance, and soil carbon across arctic-boreal permafrost ecosystems

In arctic tundra and boreal forest ecosystems vegetation structural and functional influences on the surface energy balance can strongly influence permafrost soil temperatures. As such, vegetation changes will likely play an important role in permafrost soil carbon dynamics and associated climate feedbacks. Processes that lead to changes in vegetation, such as wildfire or ecosystem responses to rising temperatures, are of critical importance to understanding the impacts of arctic and boreal ecosystems on future climate. Yet these processes vary within and between ecosystems and this variability has not been systematically characterized across the arctic-boreal region. Here we quantify the distribution of vegetation productivity trends, wildfire, and near-surface soil carbon, by vegetation type, across the zones of continuous and discontinuous permafrost. Siberian larch forests contain more than one quarter of permafrost soil carbon in areas of continuous permafrost. We observe pervasive positive trends in vegetation productivity in areas of continuous permafrost, whereas areas underlain by discontinuous permafrost have proportionally less positive productivity trends and an increase in areas exhibiting negative productivity trends. Fire affects a much smaller proportion of the total area and thus a smaller amount of permafrost soil carbon, with the vast majority occurring in deciduous needleleaf forests. Our results indicate that vegetation productivity trends may be linked to permafrost distribution, fire affects a relatively small proportion of permafrost soil carbon, and Siberian larch forests will play a crucial role in the strength of the permafrost carbon climate feedback.


Introduction
Vegetation in arctic and boreal ecosystems is changing in response to climatic drivers , Elmendorf et al 2012, e.g. Hagedorn et al 2014 and climate mediated fire regime changes (Johnstone et al 2010a, e.g. Barrett et al 2011). These changes have potentially important climate feedback implications that include effects of changes in biomass accumulation on global carbon cycling (Ma et al 2012, Pearson et al 2013, Abbott et al 2016 and albedo feedbacks on regional and global atmospheric temperatures (Chapin et al 2005, Beck et al 2011c, Pearson et al 2013. The effects of vegetation change on warming induced thaw and mineralization of permafrost soil carbon is another potentially large yet less studied positive climate feedback from arctic and boreal ecosystems. Permafrost soil temperatures under current and future climate scenarios are modulated by a series of interrelated ecological factors including vegetation cover, soil properties, surface wetness, and snow cover Jorgenson 2007, Jorgenson et al 2010). As such, the strength of the permafrost carbon-climate feedback will be strongly modulated by changes in ecosystem structure caused directly by vegetation responses to climate (Elmendorf et al 2012, Myers-Smith et al 2015 and indirectly by climate induced increases in wildfire extent (Yoshikawa et al 2003, Tchebakova et al 2009, O'Donnell et al 2011 and severity (Jafarov et al 2013, Nossov et al 2013 that affect vegetation structure (Johnstone et al 2010b).
Changes in vegetation productivity and composition can affect permafrost thermal dynamics through impacts on ground insulation, surface hydrology and shading. Widespread productivity increases, or greening trends, inferred from satellite observations  have been linked to a transition from graminoid to shrub dominance in tundra ecosystems (Frost andEpstein 2014, Myers-Smith et al 2015). In contrast, both greening and declining productivity, or browning trends, have been observed in boreal forests (Verbyla 2008, Beck and and linked with growth trends captured in tree rings (Beck et al 2011a). Expansion of woody shrubs in tundra ecosystems is associated with decreased summer soil temperature and permafrost thaw depth via canopy shading (Blok et al 2010, Myers-Smith andHik 2013). Forest cover also reduces summer soil temperatures and active layer depths (Roy-Léveillée et al 2014, Jean andPayette 2014b). Experimental removal of shrubs (Blok et al 2010) and trees (Iwahana et al 2005) has increased soil temperature and permafrost thaw depth via altered surface energy partitioning and associated changes in soil hydrology. Snow accumulation around shrubs (Sturm et al 2001) and trees (Jean and Payette 2014b) can warm soil in winter and affect surface energy dynamics of the subsequent growing season (Stiegler et al 2016). However, the insulating effects of snow cover can be reversed in larger patches (Roy-Léveillée et al 2014) or offset by other factors such as organic layer thickness (Jean and Payette 2014a) and so the large scale effects of vegetation change on winter soil temperatures remains unclear.
The extent and severity of wildfires in arctic and boreal ecosystems is increasing (Soja et al 2007, Turetsky et al 2011, but the effects on permafrost vary within and between boreal forest and tundra ecosystems (Jiang et al 2015). Combustion of vegetation and organic soil, and decreased surface albedo after fire typically lead to rapid increases in soil temperature and active layer thickness Shaver 2011, Jiang et al 2015). Degraded permafrost may recover during succession of vegetation and organic soils in the decades after fire (Rocha et al 2012, Jiang et al 2015. Fire severity influences successional pathways in ways that affect post-disturbance regrowth trajectories and plant species composition (Racine et al 2004, Johnstone et al 2010b, Jones et al 2013. Thus, fire influences on ecosystem succession may impact permafrost thermal dynamics in ways comparable to those associated with climate-induced changes in vegetation composition and distribution (Pearson et al 2013).
It is clear that ecosystem changes associated with vegetation productivity responses to climate and fire disturbance will exert strong influence over future permafrost soil temperatures, and therefore the permafrost soil carbon dynamics (Grosse et al 2011(Grosse et al , 2016). Yet despite the high degree heterogeneity in ecosystem structure, vegetation productivity trends, and fire regimes across the arctic-boreal region there has been no systematic analysis of how the spatial distribution of these properties and permafrost soil carbon co-vary in space. The objective of this study is to quantify the distribution of permafrost soil carbon in relation to vegetation type, vegetation productivity changes, and fire disturbance in order to understand the relative importance of these factors with respect to future permafrost soil carbon dynamics. In addition we seek to understand how vegetation productivity trends and fire disturbance vary with permafrost distribution and vegetation type. To accomplish these objectives we quantify the distribution of carbon in permafrost soils (C P ; defined as 0-1 m depth) (Hugelius et al 2014) in areas affected by changes in vegetation productivity and fire, by vegetation type underlain by continuous and discontinuous permafrost across the circumarctic-boreal region.

Study area delineation
The spatial domain for this study was defined using the Circum-Arctic Map of Permafrost and Ground Ice Conditions (CAMP; Brown et al 1998). This map covers the northern hemisphere poleward of 20°N latitude, and categorizes permafrost affected areas as continuous (90%-100% areal cover), discontinuous (50%-90% areal cover), sporadic (10%-50% areal cover), and intermittent (<10% areal cover) permafrost zones. Our analyses were restricted to the continuous and discontinuous permafrost zones to avoid including large areas unaffected by permafrost. Because one of our aims was to quantify the amount of permafrost soil carbon underlying areas affected by fire and vegetation change within different ecosystems we further restricted our study area using the Northern Circumpolar Soil Carbon Database v2 (NCSCD; Hugelius et al 2013). This effectively removed altitudinal permafrost at lower latitudes (i.e. in the Alps and Himalayan mountain ranges). The final study domain included approximately 11 million km 2 underlain by continuous permafrost and 3.5 million km 2 with discontinuous permafrost (figure 1, table 1).

Mapped variables
Soil organic carbon content to a depth of 1 m across the study region was quantified using the NCSCDv2 (NCSCDv2; Hugelius et al 2013). The NCSCDv2 represents a synthesis of multiple national and regional soil maps and pedon data, and is the source of the most recent and widely accepted estimates of organic carbon stored in permafrost soils. We chose to quantify organic carbon stocks in the upper 1 m because these near surface soils will be most affected by changes in ecosystem structure and disturbance regimes, and also because the 0-1 m data set is better constrained than deeper pools which are based upon fewer observations. The original dataset was in geographic coordinates on a 0.012°grid.
Variability in ecosystem types within the continuous and discontinuous permafrost zones were quantified using the Global Land Cover 2000 dataset (GLC2000; Bartholomé and Belward 2005). Percent tree canopy cover was mapped using the MODIS Vegetation Continuous Fields Collection 5.1 data product (MOD44B; Hansen et al 2003). We combined the 19 original GLC2000 land cover classes present in the study region into six aggregate classes for analysis. Aggregation of the land cover classes was informed by preliminary assessment of the spatial extent, percent canopy cover, other regional land cover products (Walker et al 2005, Sulla-Menashe et al 2011, and overall relevance to the study. Needleleaf deciduous and needleleaf evergreen classes were retained due to their extent and importance in the region, and all broadleaf or mixed leaf classes and burned forests were aggregated to an other forest class. These classes were aggregated because of their relatively small areal extents, and in the case of burned forests the dominant pre-fire vegetation type could not be determined. Mosaic forest and evergreen shrub were aggregated to the mosaic forest class. Evergreen shrub was included in this class because the two original classes were found in close proximity near latitudinal treeline (figure S1), and because the dominant shrub species in these systems are deciduous (e.g. Walker et al 2005). Moreover, the GLC2000 criteria for forest height is a minimum of 5 m, and it is likely that the evergreen shrub class consists of open woodlands with short slow-growing trees, and analysis of MODIS VCF data showed the evergreen shrub class to have a mean canopy cover of 16% (table 1), which is higher than the original mosaic forest class, and on par with needleleaf deciduous forests in the continuous permafrost zone. The remaining shrub and herbaceous classes were aggregated into the tundra class. Lastly, all non-vegetated classes that included any sort of management were aggregated into the other class. A summary of the classes is shown in table 1.
Vegetation productivity changes were quantified using a normalized difference vegetation index (NDVI) trend map produced by Guay et al (2014) that depicts linear changes in average annual June-August landscape greenness from 1982 to 2012. The trend map was derived from the Global Inventory Modeling and Mapping Studies 3rd generation (GIMMS 3g ) NDVI dataset, which is based on measurements from Advanced Very High Resolution Radiometer instruments carried by National Oceanic and Atmospheric Administration satellites (Pinzon and Tucker 2014). The per-pixel trend analysis involved pre-whitening the NDVI time series, followed by assessment of statistical significance and slope using the Mann-Kendall (Mann 1945) and Theil-Sen (Theil 1992) approaches, respectively, as implemented by the zyp package (Bronaugh and Werner 2012) in R (R Development Core Team 2014).
Fire extent and distribution across the study region for the 2000-2014 period was mapped using the Collection 5.1 MODIS Burned Area product (MCD45A1; Roy et al 2005Roy et al , 2008. Annual composites were derived from monthly files spanning the months of May through September. We excluded October through April because the study area is typically snow covered during this period and fires are highly unlikely. Monthly burned area files indicate the first day on which a pixel was determined to have burned, and annual composites were constructed by taking the minimum value for the five-month period (i.e. earliest date of burning). The MODIS product tends to underestimate burned area, and errors of omission are typically higher than errors of commission in boreal forests (Padilla et al 2014) therefore data were not screened using the embedded QA flags.

Resampling and data analysis
Prior to analyses all data were transformed from native resolution and projection to Lambert Azimuthal Equal Area Projection (LAEA) with MODIS 500 m resolution. The 500 m spatial resolution was chosen to avoid loss or creation of data associated with resampling the MODIS burned area data. The MODIS Reprojection Tool (https:// lpdaac.usgs.gov/tools/modis_reprojection_tool) was used to mosaic the monthly burned area product, reproject each file from sinusoidal to LAEA, and finally convert files to GeoTiff. These mosaics included data poleward of approximately 40°N latitude, and served as a template for resampling the remaining datasets. The CAMP, NCSCDv2, GLC2000, VCF, and GIMMS datasets were then reprojected and resampled to match the spatial extent and resolution of the burned area data set. The CAMP and GLC2000 include categorical variables and so were transformed using nearest neighbor resampling. The remaining maps included continuous variables and so were resampled using bilinear interpolation. All transformations were performed using the Raster package (Hijmans and van Etten 2013) in R. Once all data were resampled to the common projection and spatial resolution the CAMP map was reclassified to exclude the sporadic and intermittent permafrost classes, and then masked to further exclude areas not covered by the NCSCDv2. Similarly, the NCSCDv2 was then masked using the modified CAMP map so that both maps covered the same area. We then masked the remaining data sets to include only data within the delineated study area (figure 1). The zonal function in the Raster package was used to quantify the distribution and variability in ecosystem structural and functional properties within the continuous and discontinuous permafrost zones, and also ecosystem structural and functional properties within vegetation zones.

Distribution of ecosystem type and permafrost carbon
Areas of continuous permafrost store an estimated 274 Pg of C p (table 2). Within this zone, deciduous needleleaf (i.e. larch) forests contain 71.5 Pg or 26% of the total C p pool ( figure 2, table 2). An additional 141.5 Pg or 52% of the total C p pool occurs in tundra ecosystems, with the remaining 61 Pg distributed among four additional land cover classes (figure 2, table 2). Areas of discontinuous permafrost contain 65 Pg of C p , with 18.2 Pg or 28% occurring in tundra ecosystems. Evergreen needleleaf, deciduous needleleaf, and forest mosaic contain 12.6 Pg or 19%, 11.0 Pg or 17%, and 15.1 Pg or 23% of C p , respectively. Note that the mosaic forest class occurs primarily in North America because the Eurasian latitudinal treeline ecotone occurs in the deciduous needleleaf class (figure 1). Differences in total C p stocks between ecosystem types reflect, in large part, their spatial extent.

Vegetation productivity trends
Over the 30 year study period we observed greening trends in 41% and 23% of areas underlain by continuous and discontinuous permafrost respectively, while browning trends were observed in 5% of the continuous and 12% of the discontinuous zones. For all vegetation classes, the magnitude of greening trends in areas of continuous permafrost was greater than those in discontinuous permafrost (table 3). The amount of C p in ecosystems affected by productivity changes was 130 Pg in continuous and 22 Pg in discontinuous permafrost, and the majority, 117 Pg, occurred in areas of continuous permafrost exhibiting greening trends ( figure 3, table S4). Areas of increasing productivity in tundra and deciduous needleleaf forests of the continuous permafrost zone contain 58 Pg C P (49%) and 33 Pg C P (28%), respectively. By comparison, for all other vegetation classes, areas of greening or browning in either permafrost zone contained no more than 11 Pg C P .

Wildfire distribution
In contrast to decadal trends in vegetation productivity, wildfire has immediate effects on permafrost through combustion of vegetation and organic soil, and consequences of these changes on ground thaw. During the 2000-2014 period wildfire occurred over 181 000 km 2 (1.7%) of the continuous and 108 000 km 2 (3.1%) of the discontinuous permafrost zones ( figure 1(B)). Burned areas contained a total of 7.7 Pg of C P , of which half (3.8 Pg) occurred beneath deciduous needleleaf forests (figure 4). In the continuous permafrost zone, the total burned area of deciduous needleleaf forests (112 000 km 2 ; table S4) and C P contained therein (3.3 Pg) was greater than the total for all classes in the discontinuous permafrost zone (109 000 km 2 and 2.2 Pg). An additional 1.0 Pg of C P was located in tundra ecosystems affected by fire over the study period. In discontinuous permafrost deciduous needleleaf, evergreen needleleaf, and mosaic forests, burned areas contained 0.5 Pg, 0.4 Pg, and 0.5 Pg of C P respectively. Here we note that the evergreen needleleaf and mosaic classes represent North American boreal forests likely comprised of similar vegetation types that have been separated as an artifact of mapping methods, as described above. Thus, wildfire in North American boreal forests affect the largest proportion of C P in the discontinuous permafrost zone, whereas deciduous needleleaf forests in Eurasian affect the largest proportion of C P in areas of continuous permafrost. Overall the annual area burned in deciduous needleleaf forests was approximately an order of magnitude larger than any other vegetation class (figure 5), which reflects its broad spatial extent.

Discussion
Regarding the potential role of vegetation in modulating the strength of the permafrost-carbon climate feedback our results indicate that: (1) Siberian larch forests represent a substantial yet relatively understudied component of the arctic-boreal permafrost region, and (2) at present vegetation productivity trends over the past 30 years affect a much larger portion of the C P pool than fire. However, this does not imply that the contribution of a particular phenomenon (i.e. greening or fire), or ecosystem type to the permafrost carbon feedback will necessarily be proportional to its areal extent or the amount of permafrost carbon contained therein. Changes in ecosystem structure associated with productivity trends and fire will alter permafrost soil thermal dynamics in predictable ways via influences on surface energy dynamics. But these changes and their effects on permafrost will vary within and between ecosystems, and are also likely to alter soil moisture, carbon allocation, and other factors that influence permafrost soil temperatures in ways that are less well understood. Moreover, the spatial distribution of ecosystems, productivity trends, and fire are likely to change in response to continued climate warming. Nonetheless our results provide important context for understanding the role of vegetation in the permafrost carbon climate feedback.

Vegetation productivity-permafrost interactions
The structural effects of changing forest productivity and associated impacts on permafrost temperatures remain largely unstudied. While treeline advance is slower than shrub expansion, both have been observed in recent decades (Frost and Epstein 2014). In deciduous needleleaf forests, average canopy cover is just 17.4%, so productivity trends may be indicative of canopy infilling, understory vegetation responses to climate, or both. The structural consequences of widespread productivity increases in boreal larch  forests represent an important yet understudied process influencing the vulnerability of permafrost to thaw with climate warming. Unlike shrub removal in tundra (Nauta et al 2015), experimental tree removal in a larch forest has not resulted in sustained permafrost thaw, as reduced transpiration increases latent heat content of the soil (Iwahana et al 2005). Widespread increases in boreal canopy cover result in regional warming (   to canopy shading (Blok et al 2010, Myers-Smith andHik 2013). At larger scales regional albedo feedbacks (Loranty et al 2011) associated with widespread expansion of tall shrubs (e.g. >100 cm) may overwhelm local effects and lead to permafrost thaw (Lawrence and Swenson 2011). However, at the landscape scale, tall shrubs may be fairly uncommon, and their distribution is limited by a combination of climatic and environmental factors (Beck et al 2011b, Swanson 2015. Warming experiments show consistent declines in moss associated with shrub expansion (Elmendorf et al 2012), which would warm soils by removing the insulating layer mosses provide. On the other hand, shrub expansion will increase the amount of recalcitrant woody litter (Cornelissen et al 2007, Elmendorf et al 2012, changing rates of organic matter accumulation, which also affects permafrost temperatures (Jorgenson et al 2010, Grosse et al 2016. It is therefore still undetermined how shrub expansion will affect permafrost soil carbon dynamics. Increasing vegetation productivity in areas of continuous permafrost will not prevent gradual prolonged thaw and release of C P to the atmosphere that is likely to occur over the next decades to centuries . Our results suggest productivity changes, which are important yet unaccounted for in recent estimates of feedback strength (Koven et al 2015), are pervasive enough to influence the rate of C P release. Initially, greening may counteract the effects of continued climate warming on permafrost temperatures, although the effects will likely vary within and between tundra and boreal forest ecosystems. Greater carbon uptake with greening will partially offset C P release (Pries et al 2015). Changes in soil moisture may accompany vegetation changes (Iwahana et al 2005, Nauta et al 2015, which could have important influences on the fate of mineralized soil carbon (Schädel et al 2016). Thus, productivity changes will have important effects on the permafrost carbon climate feedback via a wide range of controls, both direct and indirect, on regional carbon dynamics.
More widespread and stronger greening trends, and less prevalent browning trends over our 30 year study period suggest that soil moisture retention in areas of continuous permafrost may buffer ecosystems from hydrologic stress associated with climate warming, which has been implicated as a cause for boreal productivity declines (Juday et al 2015, Walker et al 2015. Differences between continuous and discontinuous permafrost were associated with a substantial reduction (>50%) in the relative extent of greening, along with a three-fold increase in the relative extent of browning, suggesting that rising temperatures are less beneficial to tree growth in the discontinuous permafrost zone and, in some areas, could be detrimental due to increased drought-stress and more frequent disturbance by wildfire and insects. Several tree-ring and remote sensing studies have demonstrated recent reductions in tree growth in parts of northern Eurasia and North America that have been linked with temperature-induced droughtstress (Barber et al 2000, Beck et al 2011a, Berner et al 2011, Porter and Pisaric 2011, Juday et al 2015. These observed changes are congruent with model projections suggesting that higher temperatures and increased rates of disturbance will lead to a northward migration of the boreal forest biome over the 21st century (Lucht et al 2006, Tchebakova et al 2010, Pearson et al 2013. Alternatively, several recent studies in Canadian boreal forests have shown that inundation associated with permafrost degradation in poorly drained lowlands can result in productivity declines (Chasmer et   et al 2014, Frost and Epstein 2014), and that this phenomenon appears to be widespread (Helbig et al 2016). These results imply that browning trends may become more prevalent in areas where permafrost degradation results in increased surface wetness. However, productivity declines or forest loss associated with this type of hydrologic change has not been widely identified as a cause for observed browning trends. Regardless of the explanatory mechanism, increased prevalence of browning trends in areas of discontinuous permafrost will likely (i) translate into greater permafrost thaw associated with stronger coupling between air and permafrost temperatures, (ii) increase the strength of the permafrost carbon climate feedback as a result of diminishing offsets from vegetation productivity.

Fire disturbance-permafrost interactions
The vast majority of wildfire across the continuous and discontinuous permafrost zones occurs in deciduous needleleaf forests. The effects of wildfire on permafrost temperatures and C P vary with the degree of organic layer combustion, a common metric of burn severity that varies within and between vegetation types (Jiang et al 2015, Rogers et al 2015. Removal of the insulating soil organic layer leads to deeper seasonal thaw in the years to decades following fire (Viereck and Dyrness 1979, Swanson 1996, Viereck et al 2008, Jiang et al 2015, increasing the pool of thawed organic matter available to microbial decomposers. Fire effects on ground thaw also varies across landscape, soil, and vegetation classes, which affect organic layer depth, soil moisture, and ice content (Swanson 1996, Jorgenson et al 2010, Nossov et al 2013. Permafrost with high ice content may be more resistant to thawing following fire because of high latent heat content but, once thawed, ground subsidence leads to a shift in local topography (Mackay 1995) and mobilization of deeper soil carbon (Nossov et al 2013). At the same time, loss of labile carbon coupled with microbial community and abiotic changes may limit decomposition following fire (Taş et al 2014). Recovery of vegetation and soil accumulation after fire can promote permafrost aggradation enabling recovery of active layer depths to pre-fire conditions (Mackay 1995, Viereck et al 2008, Rocha et al 2012, Loranty et al 2014. In the discontinuous zone, permafrost is often thermally protected by denser vegetation canopies and deeper soil organic layers, meaning disturbance that removes these insulating layers may lead to greater thawing of permafrost (Jorgenson et al 2010).
On longer timescales, variability in fire severity may lead to shifts in ecosystem types (Mack et al 2008, Sofronov and Volokitina 2010, Beck et al 2011c, Jones et al 2013 with impacts on permafrost that are similar to the effects of productivity changes. This leads to shifts towards deciduous broadleaf dominance in evergreen conifer forests (Johnstone et al 2010b), shrub dominance in tundra (Jones et al 2013), and shifts in stand density in larch forests (Sofronov and Volokitina 2010). The associated impacts on permafrost temperatures may therefore be similar to those induced by structural changes accompanying warming induced greening trends. The influence of fire severity on post-fire succession is well documented in North American evergreen forests (Johnstone et al 2010b, Beck et al 2011c, however the total area and amount of C P affected by tundra fire is of a similar magnitude to that of North American boreal forests (evergreen needleleaf and mosaic classes), highlighting a need for improved understanding of tundra fire dynamics. There is also a relative paucity of information on fire severity effects in the extensive larch forests of Siberia where the majority of permafrost-affecting fires occur (Berner et al 2012).

Conclusions
Across the zones of continuous and discontinuous permafrost more than half of all soil carbon occurs in areas affected by productivity changes and wildfire that are likely to influence permafrost soil thermal dynamics via changes in ecosystem structure and function. Understanding the net feedback effects of such changes on climate and soil carbon vulnerability in permafrost regions is complex because processes such as greening, which tend to foster negative feedbacks to local permafrost temperatures, also have positive feedbacks to regional warming. Proportionally, areas of discontinuous permafrost exhibit less greening, more browning, and more fire relative to continuous permafrost, indicating the distribution of these processes will change over time as permafrost degrades. Our findings indicate that deciduous needleleaf boreal forests, and links between the distribution of permafrost and vegetation productivity trends are likely to influence the timing and magnitude arcticboreal terrestrial ecosystem feedbacks to climate.