Consequences of changes in vegetation and snow cover for climate feedbacks in Alaska and northwest Canada

We performed model simulations for 1900–2100 at a 1 km² spatial resolution across Alaska and western Canada, and analyzed the outputs for the 90 years between 2010 and 2099. We used a transient landscape-level model of vegetation and fire dynamics, the Alaskan Frame-based Ecosystem Code (ALFRESCO), to estimate the changes in vegetation dynamics associated with: (1) post-fire distribution of the ages of the ecosystems, (2) treeline, and (3) shrubification of the tundra. We used the Terrestrial Ecosystem Model (TEM) to calculate changes in the length of the snow season, incoming solar irradiance and longwave radiation, and changes in vegetation and soil C. We then computed changes in atmospheric heating across the study region due to projected changes in vegetation dynamics, the length of the snow season, and C storage. Below, we provide further details of the study region, the ALFRESCO and TEM models, model input datasets, and calculations of atmospheric heating.

Our study region focuses on four Landscape Conservation Cooperatives (LCCs; Millard et al 2012) in Alaska and northwest Canada including the Western Alaska LCC and parts of the Arctic, Northwest Boreal, and North Pacific LCCs (figure 1(a)). For simplicity, we refer to these as the Western Alaska, Arctic, Northwest Boreal, and North Pacific LCCs throughout the rest of this paper. These four LCCs correspond to the major ecoregions in Alaska and northwest Canada. To simulate plant growth and treeline dynamics, vegetation composition and distribution, wildfire, and biogeochemistry across the landscape, we relied on a baseline vegetation map of our study region. This map included information from the North American Land Change Monitoring System  2005 dataset (for arctic and subarctic ecosystem types) and the National Land Cover Dataset  2001 (for coastal temperate rainforest ecosystem types), the only two presently available maps with region-wide data for this study.

This vegetation map permitted us to classify the major vegetation types in our study region by LCC, including black and white spruce forest, deciduous forests, maritime coastal temperate forests, heath tundra, shrub tundra, graminoid (tussock) tundra, and wetland (wet sedge) tundra. Although further details of this classification are available in Zhu and McGuire (2016, chapter 2), here we provide an overview of the vegetation composition in each LCC.

The Arctic LCC contains all four types of tundra, with the dominant tundra types relatively evenly distributed between shrub and graminoid tundra, 108 226 km² and 119 027 km² respectively (Zhu and McGuire 2016, chapter 2). There were minor components of deciduous forest, and white and black spruce forests, primarily distributed near the southern border of the Arctic LCC, which generally corresponds to latitudinal treeline. The Northwest Boreal LCC consists primarily of deciduous forest (305 051 km², where the 'deciduous' category consists of both pure deciduous and mixed deciduous-spruce forests), with substantial coverage of pure white spruce forest (120 395 km²) and black spruce forest (88 995 km²), and minor components of all tundra types. The North Pacific LCC is dominated by maritime coastal temperate forest (76 000 km²), and small components of deciduous, white spruce, and black spruce forests. The Western Alaska LCC is relatively evenly dominated by deciduous forest (130 904 km²) and shrub tundra (119 517 km²), with small components of graminoid tundra, and white and black spruce forests. Although we do not explicitly include wetlands, the simulations with TEM (described below) include separate calibrations for three upland and lowland forest types (black spruce, white spruce, and deciduous forests), with lowland vegetation in the black spruce boreal ecosystems generally corresponding to wetlands, while graminoid and wet sedge tundra correspond to wetlands in the tundra ecosystems.

We used a transient, biogeographic model, ALFRESCO (figure 1(b); e.g., Rupp et al 2000, Johnstone et al 2011, Gustine et al 2014, Hewitt et al 2015), to simulate changes in vegetation type due to fire disturbance, natural succession, treeline advance, and shrubification from 2010 to 2100. This model includes the dominant vegetation types of Alaska and northwest Canada as described above, including black spruce forests, white spruce forests, deciduous forests, heath (barren lichen-moss) tundra, shrub tundra, wetland (wet sedge) tundra, graminoid (tussock) tundra, and maritime coastal temperate forests. The heath tundra, wetland tundra, and maritime forests are currently treated as static vegetation states in the model.

The fire module of ALFRESCO uses a cellular automata approach with separate subroutines for cell ignition and spread to simulate annual fire season activity, where additional details on the fire module are provided in Zhu and McGuire (2016, chapter 2). Following a wildfire, burned spruce forest (white or black) transitions into early successional deciduous forest, while burned deciduous forest self-replaces (figure 1(b)). Transitions in tundra are due to succession or colonization and infilling. These tundra transitions are influenced by climate and fire history, which in turn are influenced by proximity to seed source, seedling establishment and growth conditions. For the transition from tundra to forest at treeline, seed dispersal occurs within a 1 km² neighborhood. White spruce colonization/infilling may occur in both graminoid and shrub tundra with transition rates to spruce forest regulated by climate effects and basal area growth accrual. Drought stress may also influence treeline and vegetation shifts, and is modeled in ALFRESCO based on temperature and precipitation thresholds. Vegetation succession from graminoid to shrub tundra is modeled probabilistically, with a greater likelihood of transition to shrub tundra post-fire. When fires occur in the tundra, shrub tundra transitions to graminoid tundra, and while graminoid tundra is self-replacing (figure 1(b)). Vegetation transitions, defined as at least one shift in vegetation type during the projection period (2010–2099), were calculated as a percent of total area for each of the four LCC subregions. Although a grid cell generally cannot contain a fraction of a vegetation type, there may be an accumulation of white spruce basal area in tundra cells prior to transition to white spruce once the basal areal threshold has been reached. Furthermore, even though ALFRESCO uses an input vegetation map, the model is robust to initialization conditions. The map is used only at the beginning of a 1000 year model spin-up, with random age classes distributed over the landscape, across 200 model simulation replicates. This results in a rapidly converging ratio of vegetation types on the landscape by present time, with additional details provided in Mann et al (2012). The original 'mixed' class was reclassified as deciduous in the initialization vegetation map and the model only considers 'Deciduous or 'Black' or 'White' Spruce. Age is the key factor in determining the transition between 'Deciduous' or 'Black' or 'White' Spruce, and stands that are on the border between age thresholds could be considered 'Mixed', but this is presently not represented in the model. Depending on the age class of a given pixel, the associated probability of burning (cell ignition) is adjusted. The probability of a cell igniting and fire spreading to neighboring cells is also related to temperature and precipitation as well as vegetation type.

We implemented a linear coupling between ALFRESCO and TEM, allowing for the exchange of information between the models to occur in series. Here, the spatially-explicit time series of fire occurrence simulated by ALFRESCO is used to drive the process-based TEM with a Dynamic Organic Soil Module (DOS-TEM; Raich et al 1991, Yi et al 2010, Genet et al 2013). DOS-TEM then simulates the effects of wildfire, warming, and historical logging on carbon pools and aerobic carbon processes. DOS-TEM has been extensively used in arctic and boreal regions, with applications to the soil environment (e.g., Zhung et al 2001, 2003, Yi et al 2010), fire (e.g., Balshi et al 2009), and interactions between fire and carbon storage (Genet et al 2013). Although DOS-TEM has been extensively described in these studies, and most recently in Zhu and McGuire (2016), we briefly describe the model here. DOS-TEM belongs to the TEM family of process-based ecosystem models that are designed to simulate carbon and nitrogen pools of the vegetation and the soil, and carbon and nitrogen fluxes among vegetation, soil, and the atmosphere (Raich et al 1991, McGuire et al 1992). DOS-TEM is comprised of four modules: an environmental module, an ecological module, a disturbance module and a dynamic organic soil module. The environmental module simulates the dynamics of biophysical processes in the soil and the atmosphere, including soil temperature (including permafrost dynamics) and moisture conditions for multiple layers within various soil horizons including moss, fibric and humic organic horizons, and mineral horizons. The ecological module simulates carbon and nitrogen dynamics in the atmosphere, vegetation, and soil, while the dynamic organic soil layer module calculates thickness of the fibric and the humic organic layers, after soil carbon pools are altered by ecological processes (litterfall, decomposition and burial) and fire disturbance. The disturbance module simulates the effects of logging and wildfire on soil and vegetation carbon and nitrogen pools, taking into account the C lost to combustion during wildfires as described in Yuan et al (2012) and Genet et al (2013). For wildfire, the module computes combustion emissions to the atmosphere, the fate of un-combusted C and N pools, and the flux of N from the atmosphere to the soil via biological N fixation in the years following fire. In boreal forest, the amount of soil C combusted during a wildfire is determined using input data on topography, drainage, and vegetation, as well as soil (moisture and temperature) and atmospheric (evapotranspiration) environmental data (Genet et al 2013). In tundra, the rates of combustion are based on estimates from the 2007 Anaktuvuk River Fire (Mack et al 2011).

DOS-TEM takes as inputs information pertaining to climate (monthly mean air temperature, total precipitation, net incoming shortwave radiation and vapor pressure, which is converted to vapor pressure deficit in the environmental module), drainage, soil texture, elevation, and vegetation type. The model is calibrated by vegetation type for rate-limiting parameters represented in the model based on target values of C and N pools and fluxes of representative mature ecosystems. These parameters were 'tuned' until the model reached target values of the main C and N pools and fluxes (Clein et al 2002). For boreal forest communities, an existing set of target values for vegetation and soil carbon and nitrogen pools and fluxes was assembled using data collected in the Bonanza Creek Long Term Ecological Research (LTER) program (Yuan et al 2012). For the tundra communities, we used data collected at the Toolik Field Station (Shaver and Chapin 1991, Van Wijk et al 2003, Sullivan et al 2007, Euskirchen et al 2012, Gough et al 2012, Sistla et al 2013). Finally, for the maritime upland forest, we used data collected from a long-term carbon flux study in the North American Carbon Program (D'Amore et al 2012). The target values for maritime alder shrubland were assembled from Binkley (1982).

As described in Zhu and McGuire (2016, chapter 6), the model has been validated by comparing simulated soil and vegetation carbon stocks with field observations, including soil observations from the National Soil Carbon Network database (Johnson et al 2011) and vegetation biomass estimates from the Arctic LTER database (http://toolik.alaska.edu) for shrub, tussock, wet sedge, and heath tundra as well as forest inventory data collected by the Cooperative Alaska Forest Inventory for upland and lowland black spruce, white spruce and deciduous forest (Malone et al 2009). In terms of the validation of the other modules of DOS-TEM, regional applications of this model in northern high latitudes have included validation at seasonal to century scales for processes related to soil thermal activity (Zhuang et al 2001, 2003a, 2003b), snow cover (Euskirchen et al 2006, 2007), and fire (Balshi et al 2009, Yuan et al 2012, Genet et al 2013).

In this application, the biogeochemical outputs of interest include vegetation C stock estimates, derived from the sum of the above- and belowground living biomass, and soil C stocks. Soil C stocks are composed of carbon stored in the dead woody debris fallen to the ground, moss horizon, organic soil horizons, and mineral soil horizons. We computed regional area-weighted estimates of the vegetation and soil C pools in two ways: (1) by taking into account changes in the area of a given vegetation type, as derived in ALFRESCO, and (2) by assuming the area of each vegetation type remained constant to that in the first decade of the analysis, thereby permitting us to generally understand the impact of changes in vegetation on changes in C pools.

We use the environmental module of DOS-TEM to estimate the changes in the timing of snowmelt, snow return, and the duration of the snow season with methodology previously implemented in Euskirchen et al (2006, 2007, 2009a, 2009b). Snow pack accumulates whenever mean monthly temperature is below −1 °C, and snow melt occurs at or above −1 °C. At elevations of 500 m or less, the model removes the entire snow pack, plus any new snow by the end of the first month with temperatures above −1 °C. At elevations above 500 m, the melting process requires 2 months above −1 °C, with half of the first month's snow pack retained to melt during the second month. We incorporate an algorithm to estimate the date of snowmelt (or snow return) from the monthly estimates of snowpack. This algorithm uses a 'ramp' between monthly temperatures (Euskirchen et al 2007). Linear interpolations of data for monthly air temperature and the month(s) preceding snowmelt (or snow return), the month of snowmelt (or snow return), and the month following snowmelt (or snow return) are performed. For example, to calculate the date of snowmelt when all snow has disappeared by April, approximately 30 points are interpolated between mean monthly March and April air temperature to determine the 15 points for the first half of April and approximately 30 points are interpolated between mean monthly air temperature in April and May to determine the 15 points for the second half of April. The length of the snow-free season is calculated by subtracting the Julian date of snowmelt from the Julian date of snow return. To examine changes in the timing of snowmelt, snow return and snow cover duration between 2010 and 2099, we calculated anomalies based on the mean day of snow melt, snow return and number of days of snow cover duration for the period from 1980 to 2009. We calculate the anomaly for a given decade based on the difference between the day of snowmelt/return or total snow duration for a given decade from 2010 to 2099 minus mean day of snowmelt/return or total duration in the 30 years preceding the prognostic simulations, from 1980 to 2009. As discussed in Euskirchen et al (2006, 2007), the model estimates of the spatial extent and the temporal dynamics of snow cover are in agreement with those of Dye (2002).

The baseline period for both the ALFRESCO and TEM simulations was defined as the period from 1950–2009. This period corresponds to the beginning of reliable wildfire observations and statistics for this region (1950) and the end of contemporary downscaled historic climate observations (2009). The projection period was defined as the period from 2010 to 2099. Although the climate inputs to ALFRESCO include air temperature and precipitation, those for TEM include air temperature and precipitation as well as downwelling surface radiation and vapor pressure. Below, we describe the climate datasets and downscaling methods.

Historical (Climate Research Unit—CRU TS 3.1; Harris et al 2014) and projected (CMIP3) output variables of surface temperature (°C), precipitation (mm), and cloudiness (%) were downscaled via the delta method (Hay et al 2000, Hayhoe 2010) using PRISM 1961–1990 2 km resolution climate normals as baseline climate (Daly et al 2008). The delta method was implemented by calculating climate anomalies applied as differences for temperature or quotients for precipitation between monthly future CMIP3 data and calculated CRU climate normals for 1961–1990 (see http://www.snap.uaf.edu for data and details). These coarse resolution anomalies were then interpolated to PRISM spatial resolution via a spline technique, and then added to (temperature, downwelling radiation) or multiplied by (precipitation) the PRISM climate normals. The downscaled climate data were then bias corrected to account for mismatches between the historical and future data, and interpolated to a 1 km² resolution for modeling purposes.

Monthly surface downwelling shortwave radiation (R_Sin, MJ m⁻² d⁻¹) is calculated for use in TEM and for calculating changes in atmospheric heating. The algorithm takes into account the downscaled cloudiness data (as described above) and the baseline climatology global irradiance (GIRR, monthly, MJ m⁻² d⁻¹) to calculate R_Sin based on an equation from Allen et al (1998) at a monthly time step and 1 km² spatial scale:

$$\begin{eqnarray}{{R}}_{\mathrm{Sin}} & = & {\rm{GIRR}}\,* \,\{0.251+(0.509\\ & & * (1.0-\mathrm{cloudiness}/100))\}\mathrm{.}\end{eqnarray} \tag{ 1 }$$

Values of longwave radiation (W m⁻²) are calculated, including incoming longwave radiation (R_Lin), outgoing longwave radiation (R_Lout), and net longwave radiation (R_Lnet), based on Idso and Jackson (1969), and include a cloud correction factor (Parkinson and Washington 1979). The incoming longwave radiation is calculated as:

$$\begin{eqnarray}&&{{R}}_{{\rm{Lin}}}=\sigma {{{T}}_{{\rm{a}}}}^{4}\{1\,-\,0.26\exp [-7.77\times {10}^{-4}{(273{-T}_{{\rm{a}}})}^{2}]\}.\end{eqnarray} \tag{ 2 }$$

This is modified by the cloud correction factor of 1 + cn, where c is the fractional cloud cover value and n is an empirical factor set at 0.275 (Parkinson and Washington 1979). The Stefan–Boltzmann constant (σ) = 5.67 × 10⁻⁸ W m⁻² K⁻⁴. T is air surface temperature in K.

The outgoing longwave radiation is calculated taking into account the vegetation type and the emissivity values (Ɛ) for the given snow and non-snow covered vegetation type (Geiger et al 2009), in addition to the ground surface temperature (T_g) and the Stefan–Boltzmann law

$$\begin{eqnarray}&&{{R}}_{{\rm{Lout}}}={\rm{\varepsilon }}\sigma {{{T}}_{{\rm{g}}}}^{4}+{{R}}_{{\rm{Lin}}}(1-{\rm{\varepsilon }}).\end{eqnarray} \tag{ 3 }$$

The net longwave radiation is then calculated by subtracting R_Lin from R_Lout. This value is typically negative because the thermal radiative temperature of the atmosphere viewed from the air is cooler than the ground surface temperature (Betts and Ball 1997).

We selected two general circulations models (GCMs) of the best performing subset models for Alaska (Walsh et al 2008), which bound the climate scenarios from most warming and fire activity (the version 5 of the European Center Hamburg Model from the Max Planck Institute ECHAM5; Roeckner et al 2003, 2004) to least warming and more moderate fire activity (the Canadian Center for Climate Modeling General Circulation Model, CCCMA; McFarlane et al 1992, Flato 2005). Our analysis focused on a moderate forcing emission scenario (A1B) for the two GCMs from the Intergovernmental Panel on Climate Change's Special Report on Emission Scenarios (IPCC-SRES; Nakicenovic and Swart 2000). We chose only this one scenario due to the computational complexity and associated long run times with DOS-TEM. The AIB emissions scenario is a conservative choice given that global emissions more closely track the A2 scenario (Peters et al 2012, IPCC 2013), but we did not want to overestimate the amount of land area burned. Further, since we are using just this one emission scenario, we chose to bracket the variability across two bounding GCMs, choosing the CCCMA and ECHAM from a suite of GCMs, after ALFRESCO simulations indicated maximum (ECHAM) and minimum (CCCMA) area burned across the study domain with these two GCMs under the A1B scenario. We analyzed the climate data for the periods from: September–October (the fall shoulder season), November–March (the cold season), April–May (the spring shoulder season), and June–August (summer).

We assessed potential changes in atmospheric heating and feedbacks to climate due to both changes in vegetation following fire, treeline advance, and shrubification and changes in the snow season. Our calculations of atmospheric heating follow the methodology of Chapin et al (2005) and Euskirchen et al (2007, 2009a, 2009b), and are described below.

We first compiled information on albedo (α), latent (LE) and sensible heat (H) fluxes for the dominant vegetation types in our region (table 1). These parameterization sites were eddy covariance sites that collected information on the energy fluxes as well as meteorological variables, including both albedo and snow depth. We distinguished between four distinct periods based on the state of the snow cover as seen in the data: pre-snowmelt (six weeks prior to the final snowmelt, generally from March through mid-April for the boreal ecosystems and from early May to early June for the tundra ecosystems), post-snowmelt (mid-April through late August for the boreal ecosystems and from early June through mid-August for the tundra ecosystems), pre-snow return (six weeks prior to consistently snow covered ground, late August–mid October for the boreal ecosystems and mid-August–mid-September for the tundra ecosystems), and post-snow return (mid-October–early March for the boreal ecosystems and mid-September–early May for the tundra ecosytems).

Based on each period of snow cover and each value of albedo obtained from the field data and R_Sin obtained from (equation (1)), we calculated outgoing shortwave solar radiation (R_Sout) as:

$$\begin{eqnarray}&&{R}_{{\rm{Sout}}}={R}_{\mathrm{Sin}}* \alpha .\end{eqnarray} \tag{ 4 }$$

Shortwave net radiation (R_Snet) is then calculated as the difference between R_Sin and R_Sout:

$$\begin{eqnarray}&&{R}_{{\rm{Snet}}}={R}_{\mathrm{Sin}}-{R}_{{\rm{Sout}}}.\end{eqnarray} \tag{ 5 }$$

Seasonal atmospheric heating is then computed by multiplying incoming shortwave and longwave radiation by the proportion of incoming shortwave that is absorbed by the land surface times the proportion of net shortwave that is transferred to the atmosphere minus the absorbed and emitted longwave radiation (${R}_{{\rm{L}}{\rm{n}}{\rm{e}}{\rm{t}}}$):

$$\begin{eqnarray}&&\begin{array}{c}\{{R}_{\mathrm{Sin}}\,+\,{R}_{{\rm{Lin}}}\}\times \left\{\displaystyle \frac{{R}_{{\rm{Snet}}}\,\,}{{R}_{\mathrm{Sin}}\,+\,{R}_{{\rm{Lin}}}}\right\}\\ \,\times \left\{\displaystyle \frac{H\,+\,{\rm{L}}{\rm{E}}\,}{{R}_{{\rm{Snet}}}}\right\}-{R}_{{\rm{Lnet}}}\mathrm{.}\end{array}\end{eqnarray} \tag{ 6 }$$

Note that the ratios in this calculation are key, in that the terms described collapse into ({H + LE} − R_Lnet), but we have used the ratios to derive H + LE. We estimated changes in atmospheric heating due to changes in the length of the snow season, in terms of both snowmelt in the spring and snow return in the fall. We use DOS-TEM-derived values of R_Sin and literature-derived values of heat fluxes (table 1) and calculate atmospheric heating for pre- and post-snowmelt based on equation (6). We then compare pre- and post-snowmelt and pre- and post-snow return energy budgets to estimate the changes in snowmelt and snow return on atmospheric heating: {[Daily atmospheric heating post-snowmelt (or pre-snow return)] – [Daily atmospheric heating pre-snowmelt (or post-snow return)]} × [Change in snow cover duration], where the daily atmospheric heating is in units of MJ m⁻² d⁻¹ and the change in snow cover duration is in days yr⁻¹. The estimates of heating are then averaged over the length of the snow-free season, as calculated in the water balance model of DOS-TEM. Estimates of atmospheric heating are presented at the decadal time scale in W m⁻² based on the mean annual changes in the surface energy flux. Due to the computational complexity and large data files associated with the 1 km² model simulations in ALFRESCO and DOS-TEM, the vegetation transitions and snow cover are summarized annually across each LCC prior to computing the decadal atmospheric heating feedbacks. Atmospheric heating feedbacks from each LCC are then area weighted by the proportion of the LCC to the total region such that the total regional feedback is the sum of the feedbacks from each LCC.

Across all regions, both GCMs showed trends towards warmer and wetter conditions between 2010 and 2099, although these increases in temperature and precipitation were consistently greater in the ECHAM scenario than the CCCMA scenario (figure 2). The Western Alaska LCC showed the smallest increases in temperature. The Arctic LCC showed the smallest increases in the amount of precipitation in both the ECHAM and CCCMA scenarios, although as a percentage this increase in precipitation could be as large as 70% more precipitation in the future between November and March compared to the historical record during these same months. Increases in precipitation were generally greatest during the cold season from November to March, followed by the fall shoulder season from September to October (figure 2(a)). The fall shoulder season generally showed the greatest warming across each of the four LCCs, followed by the cold season, from November to March (figure 2(c)). There was a general trend towards a decrease in incoming shortwave radiation across both models (figure 2(b)).

Zoom In Zoom Out Reset image size

These differences in the amount of increase in temperature and precipitation across the two models resulted in differences in changes in the timing of snowmelt, snow return, and the duration of the snow season from 2010 to 2099 (figure 3). The snowmelt anomalies were greatest in the Western Alaska and North Pacific LCCs for the ECHAM scenario (figures 3(c) and (d)). In the Northwest Boreal LCC, the CCCMA scenario generally showed a later snowmelt anomaly from 2010 to 2099 (figure 3(b)). There was a greater change in the timing of snow return anomaly (figures (e)–(h)) compared to the change in the timing of snowmelt, with the change in snow return being smallest for the Northwest Boreal LCC under both the CCCMA and ECHAM scenarios (figure 3(f)). Overall, the change in the duration of the snow season was greatest in the Western Alaska LCC under the ECHAM scenario (an anomaly of 66 days shorter in the last decade, 2090–2099 compared to the mean snow duration from 1980 to 2009; figure 3(l)), with little change in the Northwest Boreal LCC under both the ECHAM and CCCMA scenarios (figure 3(j)).

Zoom In Zoom Out Reset image size

Across the region between 2010 and 2099, there was a shift towards a decrease in later successional (≥59 yrs) black and white spruce forests, and corresponding increases in early and mid-successional deciduous forests (0 to 59 yrs; figure 4(a); supplementary appendix 1;  figures 1–3). No change was projected in land cover for maritime coastal temperate forest, barren lichen-moss (heath) tundra, or wetland tundra due to the static nature by which ALFRESCO treats these ecosystem types. Consequently, the projected land cover changes were smaller in the North Pacific LCC compared to the other LCCs since the dominant vegetation type in the North Pacific LCC is maritime coastal temperate forest. The decreases in the black and white spruce forests were due to both drought stress and burning. Early, and eventually, mid-successional deciduous forests generally replaced the black and white spruce forests. The loss of graminoid tundra was due to treeline advance, fire, and shrubification. Below, we further describe these vegetation dynamics.

Zoom In Zoom Out Reset image size

In the Northwest Boreal LCC, there were similar losses of black and white spruce forest due to fire between the two scenarios (figure 4(b)). However, there was a much larger overall loss of black and white spruce forest in the ECHAM scenario compared to the CCCMA scenario (figure 4(c)), indicating that much of the boreal forest that burned in the CCCMA scenario eventually self-replaced by black or white spruce forest, while the forests in the ECHAM scenario burned more frequently and remained young deciduous forest. Although there was almost no fire in the Pacific LCC (figure 4(b)), there were shifts in vegetation due to natural successional processes, resulting in a decrease in black and white spruce forests and an overall increase in deciduous forests across both the CCCMA and ECHAM scenarios (figure 4(c)).

Across the Arctic, Northwest Boreal, and Western Alaska LCC, the shrub tundra was prone to fire, with greatest burn area in the Northwest Boreal and Western Alaska LCC under the ECHAM scenario (figure 4(b)). However, due to shrubification, these losses in shrub tundra from fire were alleviated in the Arctic and Northwest Boreal LCCs in both the ECHAM and CCCMA scenarios (figure 4(c)), resulting in an overall increase in the area of shrub tundra. This was not the case in the Western Alaska LCC, where the amount of shrubification was less than area of shrub tundra that burned under both the CCCMA and ECHAM scenarios, resulting in an overall loss of shrub tundra.

There were also losses in graminoid tundra across all scenarios due to fire, shrubification and treeline advance, with treeline advance, and the corresponding replacement of graminoid tundra with white spruce forest, accounting for the most of the loss. Regionally, this replacement totaled approximately 31 300 km² under the ECHAM scenario and 54 100 km² under the CCCMA scenario. The Northwest Boreal LCC comprised the greatest area where graminoid tundra either burned (figure 4(b)) or was lost due to treeline advance. Here, treeline advance accounted for a loss of approximately 22 000 km² for the ECHAM scenario and 35 000 km² for the CCCMA scenario in this LCC.

From 2010 to 2099, vegetation carbon pools increased (figure 5; supplementary appendix 2), even while the area covered by a given vegetation type may have decreased or remained static. There were two exceptions, both in the CCCMA scenario: a slight decrease in the Arctic LCC deciduous forests (−54 g C m⁻²) and in the black spruce forests of the Western Alaska LCC (−68 g C m⁻²). Those vegetation types that remained static in cover, barren lichen-moss (heath) tundra, wetland tundra, and coastal temperate forests, all showed increases in vegetation C (supplementary appendix 2). Out of all the vegetation types, the coastal temperate forests in the Pacific LCC showed the greatest increases in vegetation C, by +10 282 g C m⁻² in the CCCMA scenario and by +11 746 g C m⁻² in the ECHAM scenario, from 2010 to 2099.

Zoom In Zoom Out Reset image size

Regionally, soil C also increased from 2010 to 2099, although some vegetation types showed overall decreases in soil C. This included both the CCCMA and ECHAM scenarios for the shrub and heath tundra in the Arctic LCC, the black spruce forests in the Northwest Boreal LCC under the ECHAM scenario, and the deciduous forests in the Pacific LCC under the ECHAM scenario. Notably, in the ECHAM scenario, the Western Alaska LCC showed decreases in soil C across all vegetation types, except for the graminoid tundra and white spruce forests, and in the CCCMA scenario, soil C in this LCC decreased in the black spruce forests and heath tundra.

Regionally, from 2010 to 2099, when considering shifts across the vegetation types, the total change in vegetation and soil C was similar between the two scenarios: +1654 g C m⁻² for the CCCMA scenario and +1829 g C m⁻² for the ECHAM scenario. If the vegetation remained constant, then the estimates differed by 720 g C m⁻², with a change of +2674 g C m⁻² in the CCCMA scenario and +1954 g C m⁻² in the ECHAM scenario (figure 5).

Vegetation shifts resulted in a regional negative feedback to climate of −0.42 to 0.52 W m⁻² decade⁻¹ for the ECHAM and CCCMA scenarios, respectively (figure 6). The largest negative feedback occurred in the Western Alaska LCC under the CCCMA scenario (−0.34 W m⁻² decade⁻¹; figure 6(a)), and was due to the negative feedback associated with the decrease in shrub tundra following fire, and replacement by graminoid tundra. This negative feedback was large enough to cancel the positive feedback due to changes in snow cover in the Western Alaska LCC under the CCCMA scenario. The second largest negative feedback occurred under the ECHAM scenario in the Northwest Boreal LCC (−0.22 W m⁻² decade⁻¹; figure 6(b)). This coincided with the shift from white and black spruce forests to deciduous forests under an intensified fire regime. The only scenario where vegetation shifts resulted in a positive climate feedback was in the Arctic LCC under the ECHAM scenario. This small positive feedback (0.10 W m⁻² decade⁻¹; figure 6(b)) was due to both an advancement of treeline and shrubification, and a fire regime of low intensity.

Zoom In Zoom Out Reset image size

There was an overall positive feedback to atmospheric heating when changes in snow cover were considered in concert with changes in vegetation (figure 6). This positive feedback was larger in the ECHAM scenario across all LCCs (2.7 W m⁻² decade⁻¹ over all LCCs) compared to the CCCMA scenario (0.76 W m⁻² decade⁻¹ over all LCCs) due to the greater reduction in the length of the snow season in the ECHAM scenario (figure 3). The positive feedback was greatest in the Arctic and Western Alaska LCCs under the ECHAM scenario (∼0.85 W m⁻² decade⁻¹ from each LCC; figure 6(b)).

While we did not explicitly estimate the biogeochemical feedbacks due to changes in carbon storage in the vegetation and soils, the increases in C storage indicated that there would be an overall negative feedback to climate in all LCCs except the Western Alaska LCC, where losses in soil C were larger than increases in vegetation C (figure 5). While the biogeophysical feedbacks were small in the Western Alaska LCC (figure 6), large increases in vegetation and soil C in both the CCCMA and ECHAM scenarios by the maritime coastal temperate forests in the North Pacific LCC (figure 5) indicated that this region may be important to consider in terms of climate feedbacks related to changes in C cycling.

Across Alaska and northwest Canada, we simulated changes in vegetation composition, snow cover, and vegetation and soil C stocks under one emissions scenario (A1B) and the climate outputs from two general circulation models, CCCMA and ECHAM. Our simulations occurred at a 1 km² spatial scale, which is an unprecedented amount of detail for this region. We found that changes in snow cover duration, including both the timing of snowmelt in the spring and snow return in the fall, provided the dominant positive biogeophysical feedback to climate across all LCCs, while changes in vegetation composition generally resulted in negative biogeophysical feedbacks to climate. We expect that there would be a negative biogeochemical feedback to climate since vegetation C increased across all LCCs, and soil C increased across all LCCs except the Arctic LCC under the ECHAM scenario and the Western Alaska LCC under both scenarios. Although other studies have simulated some of these types of feedbacks in high latitudes (e.g., Groisman et al 1994, Betts 2000, Randerson et al 2006, Loranty et al 2011, 2014, de Wit et al 2014), this study provides a comprehensive comparative approach across the ecoregions in Alaska and northwest Canada. Below, we discuss additional details related to these feedbacks and consider other forcing agents that may be important to consider in future studies.

Our results emphasize the importance of considering both tundra and boreal forest wildfire. In the tundra, it was particularly important to consider the burning of shrub tundra. Although studies have suggested that an increase in shrub tundra will act as a positive feedback to climate warming since these shrubs have a lower albedo and greater sensible heat flux than graminoid tundra, previous work has not considered these feedbacks in association with increased burning of shrub tundra (Sturm et al 2005, Euskirchen et al 2009b, Bonfils et al 2012, Zhang et al 2014). Here, we found that while significant amounts of shrub tundra burned in the Arctic, Boreal, and Western Alaska LCCs (figure 4(b)), in the Arctic and Boreal LCCs, there were also increases in the amount of shrub tundra (figure 4(a)) large enough to counteract the greater albedo associated with the burned shrub tundra. In the Western Alaska LCC, there was an overall decrease in shrub tundra, resulting in a negative feedback to climate warming in the CCCMA scenario (figure 6). It is also important to note that shrub height, relative to snow cover thickness, has a key influence over feedbacks to climate, where research shows that shrubs of 0.5 m height having a significantly smaller influence on atmospheric heating than shrubs of 2 m in height (Bonfils et al 2012). Here, our parameterizations of the radiation and energy fluxes for shrub tundra (table 1) were based on a site with shrubs of ∼1 m in height, with shrubs partially exposed all winter. Thus, the positive biogeophysical climate feedbacks associated with shrub expansion would likely be less for the shorter shrubs ≤0.5 m in height, but greater than the taller shrubs >1 m in height.

Wildfire was also important in the Northwest Boreal LCC, with both the CCCMA and ECHAM scenarios indicating large losses of black and white spruce forests that were replaced by young, early successional deciduous forests. This loss of black and white spruce forest was greater in the warmer ECHAM scenario compared to the CCCMA scenario (figure 4(b)). These results support the concept that large expanses of black and white spruce forests in this region may be converted to grasslands or aspen/birch parklands in the coming century under warmer temperatures and an intensified fire regime (Chapin et al 2004, Johnstone et al 2010, Kasischke 2010). The grasslands and deciduous trees (aspen, birch) that establish following fire have a summer and winter albedo that is greater than the coniferous ecosystems (Amiro et al 2006, Liu and Randerson 2008; table 1), thereby resulting in a negative feedback to climate (Beck et al 2011, Rogers et al 2013).

Although fire proved a critical component of this analysis, the biogeophysical climate feedback due to the advance in treeline was negligible. The Northwest Boreal LCC had the greatest increase in the advancement of treeline, but the estimate of this positive feedback was miniscule (≪1 W m⁻² for the 90 years from 2010 to 2099). This differs from some other studies in other regions that have found stronger localized positive feedbacks to climate warming associated with changes in treeline (e.g., Pearson et al 2013, deWit et al 2014), and may be due to differences in forest species or faster estimated rates of treeline advance. It is also important to note that in our study area, the Brooks Range provides a natural barrier to treeline expansion (Rupp et al 2001).

It is now well known that changes in snow cover, particularly in spring, can have a strong impact on regional energy budgets due to the considerable amount of incoming solar energy reaching the snow surface (Groisman et al 1994, Déry and Brown 2007, Euskirchen et al 2007). We previously estimated that a pan-Arctic reduction in snow cover of 0.22 day yr⁻¹ during the 1970–2000 period of warming contributed an additional 0.8 W m⁻² per decade of energy (Euskirchen et al 2007), a result which is similar to that seen here in the CCCMA scenario (figure 6(a)) for this study region in Alaska and northwest Canada. When comparing climate feedbacks associated with changes in shrub tundra growth to changes in snow cover in arctic Alaska (Euskirchen et al 2009b), we found that greater aboveground biomass from increased shrub NPP produced a decrease in summer albedo, greater regional heat absorption (0.34 ± 0.23 W m⁻² decade⁻¹), and a net positive feedback to climate warming. However, this positive feedback was much smaller than the positive feedback due to a shorter snow season (3.3 ± 1.24 W m⁻² decade⁻¹; Euskirchen et al 2009b).

Our simulations with TEM indicated overall increases in carbon storage in the vegetation and soils across the study region, which would result in an estimated negative feedback to climate. Previous simulations with TEM have shown that increases in C storage are associated with the fertilization effect of a rising CO₂ concentration (e.g., Hayes et al 2011). Thus, with greater productivity under higher CO₂ concentrations, C inputs from litterfall increase, resulting in soil C storage. Historical estimates of simulated soil C stocks and heterotrophic respiration in TEM fall within the range of simulations with other biogeochemistry models across Alaska (Fischer et al 2014). One way to quantify the biogeochemical feedbacks to climate from carbon storage is to employ a fully coupled model of both biogeochemistry and climate that includes the high-level of detail in the arctic and boreal vegetation dynamics we incorporate here, while other methodology can also be implemented (e.g., Frolking et al 2006, Ward et al 2012, Bright et al 2016). It is important to recognize that biogeochemical feedbacks are a global phenomenon, while energy budget feedbacks are predominately local to regional.

The Western Alaska LCC showed loss of soil C, but there is greater uncertainty in this remote region due to a lack of observational data pertaining to climate, vegetation, soils, and disturbance. In future work in the four LCCs considered here, it may be also important to examine the losses of soil C associated directly with regional thawing permafrost and an intensified fire regime. Furthermore, it is important to consider differences in the drier uplands versus the wet lowlands. Simulations with TEM show that the 12% of the study region that is considered lowlands are sources of C, while the uplands are sinks of C (Zhu and McGuire 2016, chapters 6 and 7). We did not include changes in methane (CH₄) consumption in our estimate of feedbacks to climate. However, the wetlands in our study region are estimated as increasing sources of CH₄, and were the primary reason that the wetlands in the study region were estimated as sources of C (Zhu and McGuire 2016, chapter 7). Consequently, given the smaller proportion of lowlands than uplands, in this present regional assessment, the area remained a C sink, but an expansion of wetlands may result in a greater area that acts as a C source.

Since the version of TEM used here has been developed specifically and optimized for this region, other biogeochemical models with less specific development for our study area may show different results. For example, recent work shows large differences among land surface and biogeochemical models simulated in Arctic Alaska (McGuire et al 2012, Fisher et al 2014). Much uncertainty across models might be associated with differences in the representation of permafrost and to winter respiration (Belshe et al 2013).

Other forcing agents that we did not take into account, but which are present in our study area, include thermokarst formation, insect outbreaks, and shifts in the maritime coastal temperate forests. Increases in thermokarst formation (ground subsidence following permafrost thaw) may result in either a negative or positive feedback to climate. In boreal regions, areas that were previously forested, typically by spruce, may become bogs or fens (Jorgenson et al 2001), with a greater albedo and act as a negative biogeophysical feedback to climate. For example, data gathered from a thermokarst collapse scar bog (50 years old) and a moderate rich fen both show year round higher values of albedo and less atmospheric heating compared to spruce and deciduous forests (table 1). As a biogeochemical feedback to climate these thermokarst bogs and fens may sequester either more or less carbon (O'Donnell et al 2012, Jones et al 2013, Euskirchen et al 2014), making their net effect unknown. In arctic regions, typical thermokarst landforms include thermokarst lakes, collapsed pingos, sinkholes, and pits. A key question, both in arctic tundra and boreal forests, is if the number of newly formed thermokarst features would cover a large enough area to induce a significant change in land surface albedo and energy fluxes. Again, a related question is how the creation of these landforms alters the overall carbon storage capacity of the region (Walter Anthony et al 2014, Abbott and Jones 2015), and how this impacts the overall feedback to climate.

Across the forested portions our study area, insect damage can play an important role in tree mortality (FS-R10-FHP 2015). Under warming temperatures in recent decades, a reduction of the beetle life cycle from 2 years to 1 year has in part led to severe spruce beetle outbreaks, with widespread white spruce mortality (Werner et al 2006). Increased disturbance from insect outbreaks will also shift the forest to a younger age classes (Boucher and Mead 2006, Kurz et al 2008) increasing the surface albedo, possibly leading to a negative feedback on atmospheric heating. However, this negative feedback could be counterbalanced by reduced carbon storage capacity in these mature forests with insect outbreaks (Kurz et al 2008).

In future studies that include the North Pacific LCC, it may also be important to include shifts in the maritime coastal temperate forests, which have experienced documented mortality related to climate change (Hennon et al 2012). For example, mortality in yellow-cedar in this region has been related to earlier snowmelt and later freeze-up, where over 200 000 ha of forests have died in the past 100 years (Hennon et al 2006). Moreover, these forests can attract and hold cloud cover, which reflect heat and sunlight (Lawford et al 1996). This potentially strong negative feedback to climate, combined with these forest's carbon storage, are important factors to consider in light of net feedbacks to climate.