The impacts of recent permafrost thaw on land–atmosphere greenhouse gas exchange

Permafrost thaw and the subsequent mobilization of carbon (C) stored in previously frozen soil organic matter (SOM) have the potential to be a strong positive feedback to climate. As the northern permafrost region experiences as much as a doubling of the rate of warming as the rest of the Earth, the vast amount of C in permafrost soils is vulnerable to thaw, decomposition and release as atmospheric greenhouse gases. Diagnostic and predictive estimates of high-latitude terrestrial C fluxes vary widely among different models depending on how dynamics in permafrost, and the seasonally thawed ‘active layer’ above it, are represented. Here, we employ a process-based model simulation experiment to assess the net effect of active layer dynamics on this ‘permafrost carbon feedback’ in recent decades, from 1970 to 2006, over the circumpolar domain of continuous and discontinuous permafrost. Over this time period, the model estimates a mean increase of 6.8 cm in active layer thickness across the domain, which exposes a total of 11.6 Pg C of thawed SOM to decomposition. According to our simulation experiment, mobilization of this previously frozen C results in an estimated cumulative net source of 3.7 Pg C to the atmosphere since 1970 directly tied to active layer dynamics. Enhanced decomposition from the newly exposed SOM accounts for the release of both CO2 (4.0 Pg C) and CH4 (0.03 Pg C), but is partially compensated by CO2 uptake (0.3 Pg C) associated with enhanced net primary production of vegetation. This estimated net C transfer to the atmosphere from permafrost thaw represents a significant factor in the overall ecosystem carbon budget of the Pan-Arctic, and a non-trivial additional contribution on top of the combined fossil fuel emissions from the eight Arctic nations over this time period.


Introduction
Nearly 25% of the northern hemisphere land surface is underlain by perennially frozen ground (Zhang et al 1999). The northern high-latitude permafrost region has the potential for substantial influence on the global climate system through ecosystem feedbacks in energy, water and carbon (C) cycling . Climate warming over the high latitudes has been pronounced in recent decades (Serreze and Francis 2006) and is expected to increase at a greater rate than the global average over the next century (AMAP 2011). The response of the ecosystem C cycle to future warming is of particular concern, as the large stocks of C in permafrost soils  are vulnerable to thaw (Harden et al 2012) and release as atmospheric greenhouse gases (GHG) through decomposition. The release of C as CO 2 and CH 4 from thawing permafrost is a potentially important, but highly uncertain, positive feedback to climate warming .
Many coupled climate-carbon models project that the northern high latitudes will serve as a substantial land C sink during the 21st century because both climate warming and elevated global atmospheric CO 2 concentration favor increased vegetation productivity and C uptake in the region (Qian et al 2010). However, these models do not wellrepresent, or have not accounted for, important ecosystem dynamics associated with warming-driven permafrost thaw (Koven et al 2012). Back-of-the-envelope calculations suggest potentially substantial GHG forcing as a result of C emissions from permafrost thaw over the next century (Schuur et al 2013). Where coupled and uncoupled land surface models have incorporated the impacts of permafrost thaw in their simulations, they project that this region will become a net C source to the atmosphere of 12 to well over 100 Pg C over the next 100 years (e.g., MacDougall et al 2012, Zhuang et al 2006, Koven et al 2011, Schaefer et al 2011. These models simulate permafrost thaw as a thickening of the active layer, but estimates of the associated C emissions do not distinguish between C in the active layer and C from the permafrost below. Ultimately, the net impact of the permafrost carbon feedback will depend on the balance between future C uptake by arctic and boreal vegetation and the quantity (and form) of soil C emissions under a warming climate.
The influence of permafrost thaw on the high-latitude C cycle is complex (figure 1). The increase in active layer thickness (ALT) exposes more previously frozen soil organic matter (SOM) to thawing and decomposition, thus accelerating C and nutrient cycling throughout the system. Decomposition of this SOM results in the release of CO 2 from heterotrophic respiration (HR) under oxic conditions and the release of CH 4 from methanogenesis ( f CH 4 ) under anoxic conditions. The enhanced CO 2 and CH 4 release associated with decomposition of previously frozen C represents a positive feedback to climate warming. If decomposition is incomplete, some of the additional C may instead be lost from the land system to neighboring river networks as dissolved organic carbon (DOC). Enhanced decomposition also increases nitrogen (N) mineralization (N min ), releasing more available N (N avail ) from previously frozen SOM to be taken up by vegetation (N up ) to increase net primary production (NPP). This enhanced uptake of CO 2 by vegetation from N-stimulated NPP leads to a negative feedback to warming. The effect of permafrost thaw on net GHG depends on the relative influence of these competing feedbacks (Schuur et al 2008).
Recent estimates of the northern high-latitude C budget in studies by McGuire et al (2010) are based on a modeling approach that incorporates these complex mechanisms described above. Using the same approach, the study by Hayes et al (2011) suggests that the Pan-Arctic land-based C sink is currently weakening, which is due partly to temperature-driven increases in the decomposition rates of unfrozen SOM and partly to the mobilization of previously frozen SOM. Here, we disentangle these two mechanisms through a modeling experiment in order to assess the direct impact of permafrost thaw on the Pan-Arctic C budget. The experiment used a process-based ecosystem biogeochemistry model that includes an explicit representation of climate-driven dynamics in ALT, which directly influences the quantity of SOM that is available for decomposition as well as the mineralization and immobilization of N tied to decomposition (Hayes et al 2011). We report the results of an experiment where we compared a fully transient simulation against a reference simulation where ALT was held constant through the analysis period. The difference between the transient and reference simulations represents a quantitative estimate of the direct effect of active layer dynamics on the Pan-Arctic scale C budget and its component fluxes.

Study domain
The domain of northern permafrost covers 12.1 million km 2 and includes areas of mostly tundra (57%), and boreal forest (39%), each of which include wetland areas (14%), underlain by the continuous (76%) and discontinuous (24%) permafrost zones (Brown et al 2002) representing the highlatitude land areas in Asia (66%), North America (33%), Europe (<1%) and the North Atlantic Islands (<1%). The simulations described here do not include land below 45 • N, which thus excludes some permafrost area associated with high elevation (e.g. the Tibetan Plateau) from this analysis. Surface air temperature (Tair) anomalies, based on climate reanalysis data sets (Mitchell and Jones 2005), demonstrate a warming trend for all seasons across the domain over both the continuous and discontinuous permafrost zones (table 1). On average annually, the Pan-Arctic domain was about 0.8 • C warmer in the 1990s and 1.3 • C warmer between 2000 and 2006, relative to the 1960-1989 reference period. The greatest positive surface air temperature anomalies occurred over the spring months during the 1990s, whereas the pattern shifts to a noticeable increase in warming over the fall and winter months during the 2000-2006 time period.

Model description
The model experiments in this study build on previous model results for high-latitude ecosystems using the updated version 6 of the Terrestrial Ecosystem Model (TEM6; McGuire et al 2010). Several modifications from previous versions  have been implemented in TEM6; notably for this study, these include changes in the representation of C stored in SOM and its availability for decomposition, as well as in the simulated influence of multiple inorganic N pools, temperature and soil freeze-thaw state on plant gross primary productivity and maintenance respiration (Hayes et al 2011) TEM6 also considers DOC production as a function of incomplete decomposition, and DOC export ( f DOC) associated with runoff (Kicklighter et al 2013). The Methane Dynamics Module of the TEM (MDM-TEM) explicitly simulates the processes of CH 4 production and CH 4 oxidation as well as the transport of the gas between the soil and the atmosphere to estimate net biogenic CH 4 emissions (Zhuang et al 2007). The modeling framework also incorporates sub-grid heterogeneity in unique vegetation communities and disturbance histories (Hayes et al 2011).
In previous versions of the TEM, the potential influence of permafrost on the availability of soil organic matter to decomposition has not been considered. Estimates of heterotrophic respiration had assumed that the entire soil organic carbon (SOC) pool is susceptible to decomposition, controlled by soil temperature and moisture within the maximum rooting zone of the dominant vegetation (McGuire et al 1997). In permafrost regions, however, the thickness of the active layer is often less than the maximum rooting depth (0.3-2.9 m) assumed for these biomes (Vörösmarty et al 1989) indicating that some of the soil organic carbon stored within the rooting zone is protected from the seasonal or inter-annual variation in temperature or soil moisture, but may become more susceptible to decomposition with warming-induced permafrost thaw.
To specifically improve the simulation of C-N cycling in TEM6 for areas underlain by permafrost, the total thickness of unfrozen soil varies seasonally as a function of modeled below-ground thermal dynamics. In referring to the unfrozen soil layers, here we use the terminology as defined by Nelson and Hinkel (2003), where 'thaw depth' is the total thickness of thawed soil at a given point in time (monthly in TEM6) within a particular year. The term 'active layer' refers to the layer of ground between the surface and the permafrost that undergoes seasonal freezing and thawing; active layer thickness (ALT), then, is considered here to be equal to the maximum thaw depth as simulated by the model for a given year. In TEM6, ALT is first set to the vegetation-specific rooting depth (which also varies by soil texture); then, the soil thermal module (Zhuang et al 2003(Zhuang et al , 2002 is run to determine which portions of the profile are frozen versus unfrozen, for each cohort in each month over the simulation period. The model allows up to two layers of unfrozen soil between a frozen layer and the top of the permafrost or the rooting depth (whichever is more restrictive); thaw depth is then the sum of the thickness of the unfrozen layer(s) between the soil surface and the top of permafrost (or rooting depth). As with other land surface models that simulate permafrost C dynamics (e.g., Schaefer et al 2011, Koven et al 2011, TEM6 does not explicitly 'define' permafrost occurrence or characteristics, but rather infers its depth for each cohort in a grid cell based on the zero degree centigrade isotherm in the soil temperature profile of that cohort. Thaw depth determines the relative amount of soil organic carbon (SOC) available for decomposition along with Idealized representation of the distribution of total soil organic carbon (SOC) (dotted line) and reactive SOC (solid line) over a soil profile and its relationship to permafrost dynamics when the lability of soil organic matter decreases with depth. Note that as the active layer increases (top layer of permafrost changes from dashed line to dash-dot-dot line), more SOC becomes exposed to decomposition and that a larger proportion of reactive SOC in the profile may be exposed earlier than indicated by the distribution of total SOC. the amount of inorganic nitrogen available for uptake by plants and microbes. The depth to the frozen layer also represents the bottom impermeable soil layer for calculating water storage and associated runoff in the model.
To account for variations in SOC content and lability with depth TEM6 uses hyperbolic functions with biome-specific horizon information for describing the distribution of total SOC in the soil profile TEM6 only considers the quantity of soil organic matter (SOM) found within the maximum rooting depth as prescribed for each vegetation type, and assumes that the amount SOC available for decomposition is dependent upon the total thickness of the unfrozen layer(s) (thaw depth) relative to the rooting depth (figure 2). Because most of the SOM is located towards the top of the soil profile, the additional SOM exposed to decomposition at deeper depths in the profile will be less per unit depth than that exposed at shallower depths. This representation of SOC distribution is based on data for mineral soils (Jobbágy and Jackson 2000), but is a simplifying assumption that is used for all soils in TEM6. This exponential decrease in SOC and root density with depth is supported by data from mineral soils in boreal forests (Yi et al 2009), but is not always representative of organic and/or cryoturbated soils in boreal and tundra ecosystems. TEM uses this depth profile relative to the active layer to distinguish between 'frozen' and 'unfrozen' SOC in the active layer, but does not otherwise explicitly represent permafrost or 'ancient' frozen carbon in its initialization or spin-up.
TEM6 assumes that the SOM at shallower depths contains more labile components so that it is more 'reactive' to decomposition than SOM found at deeper depths. TEM6 thus represents the distribution of reactive SOC to mimic the distribution of total SOC through the soil profile, with the proportion of reactive SOC to total SOC varying by biome (Hayes et al 2011). The depth profile for the ratio of reactive to total SOC shows a wide range of variation across ecosystems and soil types (e.g., Lavoie et al 2011, Dutta et al 2006, and this assumption represents a 'null hypothesis approach'. In effect, the ratio can decrease (e.g., from fibric to humic to older unfrozen minerals soils), remain the same (e.g., Yedoma deposits where freshly frozen plant material was rapidly buried deep into permafrost), or increase (e.g., from modern, near-surface organic horizons to Yedoma and older frozen soils) with depth. The constant ratio implementation in TEM6 characterizes the intermediate case between these various possibilities.
The TEM was calibrated to site-specific vegetation parameters (Euskirchen et al 2006) and extrapolated across the Pan-Arctic study area based on spatially explicit time series data organized on a 0.5 • latitude by 0.5 • longitude grid, as defined by the climate driver data set. The model was driven by spatially referenced information on atmospheric chemistry, climate, elevation, soils, vegetation, disturbances and land cover to estimate monthly terrestrial C, N, and water fluxes and pool sizes. The climate driver data sets were obtained from the Climate Research Unit (CRU) corrected reanalysis data (Mitchell and Jones 2005) for monthly air temperature ( • C), precipitation (mm), and incident short-wave solar radiation (W m −2 ). The ozone (O 3 ) pollution data set, represented by a measure of the accumulated hourly ozone levels above a threshold of 40 ppb (AOT40 index), is based on Felzer et al (2005). The atmospheric N deposition data are based on Van Drecht et al (2003). The initial sub-grid cohort data set incorporated upland and wetland vegetation community types based on the 1 km Global Land Cover Characterization data set (Loveland et al 2000) and a 1 • grid fraction inundated database (Matthews and Fung 1987). Cohort dynamics were driven by fire disturbance based on mapped area burned data (Balshi et al 2007, van der Werf et al 2010) as well as modeled wood harvest and land use transitions (Hurtt et al 2006).

The simulation experiment
The simulations were run over all land grid cells above 45 • N, but the analysis presented here includes only the subset of those grid cells that are intersected by the continuous and discontinuous permafrost zones. As TEM does not explicitly define permafrost presence/absence spatially within and across grid cells, the Circum-Arctic Map of Permafrost and Ground-Ice Conditions (Brown et al 2002) was used to prescribe the continuous and discontinuous zonation used to summarize the simulation results in this analysis. To quantify and analyze the direct impact of increasing ALT on carbon flux, we compared results from a reference active layer simulation, in which thaw depth is held constant through the analysis period, to the simulation results reported previously by Hayes et al (2011), which includes explicit simulations of transient, climate-driven active layer dynamics. The transient active layer simulation, referred to in this study as 'S T ', was driven by temporally varying atmospheric CO 2 concentration, tropospheric ozone levels, N deposition, climate (Tair, precipitation and incident short-wave radiation), disturbance (fire and forest harvest) and agricultural (croplands and pasture) establishment and abandonment. The reference simulation conducted for this study, referred to as 'S R ', was identical to the S T except that thaw depth was held at a constant value for each month in each cohort over the analysis period. That constant value was calculated as the average thaw depth for each cohort in each month between 1960 and 1970, which represents the period before the onset of increasing warming trends over most of the Pan-Arctic region (Serreze andFrancis 2006, Euskirchen et al 2007).
Prior to 1970, the S T and S R were identical, using the same 'spin-up' run to reach dynamic equilibrium of the carbon pools and soil temperatures. From this same 1970 model state, modeling forward to 2006 we forced the constant monthly value of thaw depth for each cohort in the S R rather than have it be determined by the soil thermal module as in the S T ; otherwise, the simulation set up and the driver data used were duplicate between the two simulations. This creates an apparent contradiction where the soil temperatures themselves are the same between the two simulations, but we can have the 'frozen' layer in the S R be above 0 • C. This contradiction does not affect our analysis, however, since: (1) TEM does not explicitly define permafrost, as explained above, and (2) in this study we are only interested in the difference between simulation results for C and N fluxes as controlled by varying the amount of thawed C available for decomposition, which is determined by the dynamic (S T ) versus static (S R ) active layer. Active layer dynamics will also determine water storage and drainage; and, although they are not foci of this analysis, these dynamics will drive indirect effects on C and N cycling. In setting up the modeling experiment this way, the direct impact of active layer dynamics on ecosystem C and N fluxes could be determined by subtracting the simulated S R result for a particular TEM output variable from the associated simulated S T result for the same variable. We analyzed these 'S T minus S R ' results for the time period following the holding of thaw depth constant (in the S R simulation), i.e. 1970-2006. We calculated the amount of SOC in the 'frozen' soil pool using the inverse of the calculation described above for determining the amount of SOC within the active layer (figure 2). Since the frozen SOC pool is held constant in the S R simulation, the transfer of C from the frozen to unfrozen C pool ( f Thaw) over the analysis period was determined directly from the S T simulation. For both model experiments, we produced simulated, monthly estimates of land-atmosphere CO 2 exchange and lateral C transfer (e.g., f NPP, f HR and f DOC) using TEM6 and CH 4 exchange ( f CH 4 ) using MDM-TEM, employing the full C budget account approach as described by McGuire et al (2010). The MDM-TEM description and parametrizations for both upland and wetland ecosystems are documented in previous studies (Zhuang et al , 2006. The MDM-TEM used the same climate and land cover data sets described above for the TEM6 simulations as inputs, and was driven by simulated LAI and NPP output data from TEM6 (McGuire et al 2010).

Results
Over the 1970-2006 time period, the simulated maximum annual active layer thickness (ALT) showed greater change in the discontinuous than in the continuous permafrost zone ( figure 3). For the transient active layer simulation (S T ), ALT in the continuous permafrost zone increased by 5.8 cm from 1970 to 2006 (57-63 cm, respectively) versus the 9.5 cm increase over this time period in the discontinuous zone (70-79 cm). Although ALT shows a greater increase in the discontinuous zone, the model estimates a slightly larger cumulative amount of thawed carbon ( f Thaw) added to the unfrozen pool over the 1970-2006 time period in the continuous (5.9 Pg C) than the discontinuous zone (5.7 Pg C), as a function of the large area and soil carbon stocks there. Correspondingly, the results show a larger impact of active layer dynamics (i.e. S T -S R ) on ecosystem fluxes in the continuous zone, where the source effect on the vertical, land-atmosphere net carbon exchange of CO 2 and CH 4 ( f NCE; 2.6 Pg C) is more than twice as much as that in the discontinuous zone (1.1 Pg C). In both zones, the imbalance is a result of the sink effect of active layer dynamics on CO 2 uptake ( f NPP; 0.20 and 0.076 Pg C in the continuous and discontinuous zones, respectively) being outweighed by the source effect on CO 2 ( f HR; 2.9 and 1.2 Pg C) and CH 4 ( f CH 4 ; 0.017 and 0.013 Pg C) emissions tied to thawing of the frozen soil C pool from 1970 to 2006.
The S T results in a greater increasing annual trend in ALT over the discontinuous versus continuous zone (figure 4(a)). ALT increased from an average of 57.6 cm over the 1970s decade to 62.1 cm over the 2000s across the continuous permafrost zone, with a significant ( p < 0.01) increasing annual trend of 0.18 cm yr −1 from 1970 to 2006. Over the discontinuous zone, the trend, at 0.37 cm yr −1 , is more than twice that of the continuous zone, increasing from an average of 72.8 cm over the 1970s decade to 78.5 cm over the 2000s. In addition to a deeper ALT over time, the simulations show a strong bimodal pattern in the changes of monthly Thaw depth over the seasonal cycle for both permafrost zones ( figure 4(b)) where the greatest increases in Thaw depth occur in spring (May in the discontinuous zone and June in the continuous) and fall (October in both zones). For both zones, Thaw depth increased most in the spring during the 1980s but by the end of the simulation (2000s) the seasonal pattern is trending toward a greater increase during the fall season.
Over recent decades, our results show that active layer dynamics associated with permafrost thaw are causing the Pan-Arctic domain to become an increasing C source in both the continuous and discontinuous permafrost zones (figure 5). To compare the component fluxes of land-atmosphere exchange, each is shown in units of CO 2 equivalent representing the global warming potential of each greenhouse gas over a 100-year time scale, where one g of CH 4 (0.75 gC) is equivalent to 25 g of CO 2 (6.8 gC) (IPCC 2007). A positive value represents a source effect on C exchange from the land to the atmosphere whereas a negative value represents a net terrestrial C sink. The effect of increasing ALT on net greenhouse gas forcing ( f GHG) over the 1970-2006 time period is dominated by increasing CO 2 release from heterotrophic respiration (HR), whereas the contribution from methane flux ( f CH 4 ) was relatively small and without a significant trend in recent decades. The increasing effect on HR is apparent in the simulations since the 1970s and is much 0 Figure 3. Results from the S T and S R simulations comparing cumulative carbon stock changes and fluxes (Pg C) over the 1970-2006 time period for the (a) continuous and (b) discontinuous permafrost zones. Changes in the Live Biomass C stock are balanced by inputs from net primary production ( f NPP) and outputs from the transfer of dead organic matter to the unfrozen soil pool ( f Litter)-plus combustion in fire as well as harvested product removal (not shown). Changes in the Unfrozen Soil C stock balance the input from the live biomass pool against atmospheric emissions from heterotrophic respiration ( f HR) and methane fluxes ( f CH 4 ), along with lateral losses in the form of dissolved organic carbon export ( f DOC)-plus combustion in fire (not shown). We also show the input to the unfrozen soil pool from the thawing of the Frozen Soil C stock ( f Thaw) related to increasing annual thickness of the active layer (ALT, cm), which is shown for 1970 and 2006. The net C exchange ( f NCE) between land and atmosphere is uptake by vegetation ( f NPP) minus emissions from decomposition ( f HR plus f CH 4 ) and fire (not shown). In this diagram, aggregated fluxes and stock changes do not always balance due to our omission, for the purposes of clarity, fluxes from fire and harvest. stronger overall than the would-be compensating effect on net primary productivity (NPP), which shows a lag that begins after 1990. In the model experiment, the increase in HR is a result of the direct effect of adding additional C from the frozen to the unfrozen soil, or 'active', C pool.
The increase in NPP in the latter decades of the S T simulation is an indirect effect of increasing SOM decomposition resulting in greater mineralization of N that is then available for plant uptake. Analysis of these changing seasonal dynamics in C-N cycling by decade ( figure 6) shows an increase in gross N mineralization from the 1980s to 2000s largely occurring during the spring and fall months, which corresponds to the seasonal pattern by decade of S T -S R HR. The available N pool is the balance of additions from net N mineralization (gross mineralization minus immobilization) and reductions from plant N uptake. The model experiment shows an effect of greater available N accumulation over the winter months by decade, with a larger drawdown during the growing season that supports the increase in NPP.

Discussion
The impact of permafrost carbon release on global climate forcing is a function of (a) the rate of permafrost thaw, (b) the magnitude and rate of permafrost C mobilization, and (c) the form of thawed C release through decomposition as CO 2 or CH 4 . Our model simulation experiment allows an analysis of the sensitivity of each of these components to active layer dynamics. Averaged across the full domain, the simulated increase in annual maximum active layer thickness (ALT), as a proxy for permafrost thaw, was 6.8 cm over the time period of analysis . The increase in ALT transferred a total of 11.6 Pg C from the frozen soil C pool to the unfrozen soil C pool over this time period. We are able to make a quantitative estimate of the direct impact that adding this thawed C to the active layer pool had on ecosystem fluxes by comparing against a reference simulation where ALT was held constant. The result is a domain-wide, cumulative release of 4.0 Pg C, or 34% of the total thawed C, through decomposition ( f HR plus f CH 4 ) over the 37 years of the model experiment. Based on the 66.2 Pg C of cumulative decomposition flux from the transient simulation, the model experiment estimates that 6.1% is contributed by thawed permafrost C. According to the simulations, almost all of the thawed C release was in the form of CO 2 (from f HR), with only 1.2% as CH 4 . The reference simulation estimates a cumulative net sink of 4.1 Pg C with a constant active layer; the increasing ALT in the transient  (1980s, 1990s, and 2000-2006). Both plots summarize the results over the Pan-Arctic domain by continuous and discontinuous permafrost zones. simulation resulted in a 92% reduction in accumulation, or only a 0.3 Pg C sink over the time period of analysis.
Accounting for the greater global warming potential of CH 4 over a 100-year time scale means that f CH 4 represents a larger fraction (5.9%) of the cumulative estimated greenhouse gas forcing ( f GHG) from permafrost C emissions (15.6 Pg CO 2 eq) than when just the relative mass of C is considered. The forcing from emissions is offset by 6.4% from the sink effect of permafrost thaw that results increased f NPP. The net f GHG as a direct result of permafrost thaw is estimated in this study as 14.6 Pg CO 2 eq total since 1970. From 1990 to 2006, the average annual net forcing is 0.53 Pg CO 2 eq yr −1 , which represents a significant factor in the estimated 0.64 Pg CO 2 eq yr −1 pan-arctic GHG forcing estimated from full budget accounting (McGuire et al 2010). The net GHG forcing from active layer dynamics estimated in this study represents an additional 6.9% contribution on top of the combined 7.8 Pg CO 2 eq yr −1 fossil fuel emissions from the eight Arctic nations (Boden et al 2012), over the same time period.

Rate of permafrost thaw
There are several recent model studies that project rates of permafrost thaw under different climate scenarios (e.g., Schneider von Deimling et al 2012, Schaefer et al 2011, Burke et al 2012, Koven et al 2011 but few include retrospective analysis. Projections of ALT increase over the 21st century range from 5 to 30 cm per decade, based on the studies synthesized by Schaefer et al (2011). In their modeling study, Koven et al (2011) report a mean rate of ALT increase north of 60 • N for the 1990s and 2000s at 5 cm per decade, which is greater than the average 3 cm per decade estimate from this study over the same time period for the continuous and discontinuous zones combined. In this study, the greater change and increasing trend in simulated ALT over the discontinuous than in the continuous permafrost zone is to be expected since permafrost temperatures tend to be closer to the thaw point in the former (Romanovsky et al 2010). These simulated changes in ALT vary at the regional scale with the largest changes occurring in Europe and the smallest changes occurring in North America.
A comparison of simulated ALT with field-based measurements at sites across the Pan-Arctic collected by the Circumpolar Active Layer Monitoring (CALM) network  since the 1990s indicates very good agreement with the mean and spatial variability of the observational data across all North America (Alaska and Canada) sites, as well as most in the Asia (Siberian Russia) sub-region (figure 7). In North America, ALT has been relatively stable in high Arctic areas  with increasing trends restricted primarily to sites in the Alaskan interior (Jorgenson et al 2001). Increasing trends have also been observed for several sites in Scandinavia (Åkerman and Johansson 2008) and Russia (Romanovsky et al 2010). In Asian Russia, there are a few sites that show deeper ALT (greater than 80 cm) than the range of variability simulated by TEM6 in both tundra and taiga vegetation types. Based on the relatively few CALM sites in Europe (not shown), simulated ALT is generally an overestimate in the discontinuous zone and an underestimate compared to specific sites in Fennoscandia and western Russia that show deeper ALT than the rest of the CALM sites. Otherwise, both modeled and observed ALT over North America and Asia range from approximately 30-60 cm in tundra types and about 40-100 cm in forested (taiga) types for the majority of sites over the 1990-2006 time period.

Mobilization of permafrost carbon
The magnitude and rate of permafrost C mobilization depends on the amount, distribution and quality of permafrost C stocks. The northern permafrost region stores approximately 1700 Pg of soil C according to current estimates that include deep deposits . Our model estimates a total of 560 Pg C in total soil organic C at the start of the analysis period (i.e., 1970). This is much less than the total circumpolar estimate, but our estimate is for the continuous and discontinuous zones only and covers 12 million km 2 versus the 18 million km 2 in the study by Tarnocai et al (2009). Soil C stocks are determined in TEM6 according to the depth of the rooting zone, which varies by ecosystem type and soils but can be generally considered to represent an approximately 2 m profile. As such, our estimate is more comparable to the 496 Pg C to 1 m depth and 1024 Pg C to 3 m to depth, as reported by Tarnocai et al (2009). Soil C stocks in TEM6 are spun up over thousands of years according to ecosystem-specific calibrations, but otherwise are not initialized to capture special cases of high C content soils such as in peatlands, deltas or Yedoma.
The distribution of C in the soil profile is represented in TEM6 according to a generalized relationship for mineral soils (figure 2), which can better incorporate the range in conditions for permafrost soils by incorporating cryoturbation processes (Koven et al 2009), dynamic organic layers (Yi et al 2010) and parameters based on depth profiles from regionally specific soil pedon data sets (Harden et al 2012). The rate of mineralization of thawed soil C (turnover rate) is controlled in part by its quality (decomposability). For the full study domain over the 37-year time period in this analysis, the model has an average annual turnover rate (k equals f HR over the pool size) of 3.2% of the accumulating S T -S R unfrozen soil C pool, equivalent to a turnover time (1/k) of 31 years. This turnover time is integrated across the 'single-box' soil model in TEM6, which does not distinguish between C pools of varying lability. As a general comparison, though, the turnover time in this analysis is longer than the 5-15 years of the 'slow turnover pool' as determined by Schädel et al (2014) in their synthesis of permafrost soil incubation studies.

Permafrost carbon fluxes
The net impacts of warming-driven ALT increase on ecosystem C fluxes demonstrated in this study at the inter-annual and decadal-scales are a reflection of the changing seasonal dynamics of thaw and associated C-N cycling. In the 1990s, the positive temperature anomaly was disproportionately found in the spring months, and the results of the model experiment show a corresponding increase in Thaw depth and f HR. The recent trend in later-season Thaw depth increase and associated f HR is an apparent reflection of the shift in the 2000s toward warmer temperatures in the fall season. While there is a positive trend in the impact of increasing ALT on Figure 6. Plots subtracting the S R from the S T simulation results in this study to illustrate the impact of a dynamic active layer on the seasonal cycle of coupled carbon and nitrogen dynamics by decade (1980s, 1990s, and 2000s). Shown are the decadal average seasonal cycles of (a) HR and (b) NPP (gC m −2 month −1 ), which are coupled to (c) N mineralization rate and (d) change in the stock of plant available N (gN m −2 month −1 ). All flux effects are summarized over the Pan-Arctic domain by the continuous (solid lines) and discontinuous (dashed lines) permafrost zones.
f HR from the start of the analysis period (1970), the increasing effect on NPP is not as strong overall and shows a delay in the timing starting more in the latter decades (1990s). In the model experiment, gross N mineralization follows f HR with an increasing effect in spring and fall. In field-based experiments, a limited NPP response to nutrient addition in tundra ecosystems is most likely linked to a strong microbial immobilization response (Shaver andChapin 1980, Chapin et al 1986). Our results show immobilization increases along with gross mineralization in spring and fall by successive decade, though the plant-available N pool only eventually shows a net accumulation over the winter months in the 1990s-2000s. The model experiment suggests, then, that it is this accumulating available pool that allows greater N uptake by plants during the growing season. The details of the seasonal timing and complex interactions between mineralization, N storage and plant uptake are difficult to disentangle, however, both in the modeling environment as well as in the field (Natali et al 2012).
Observational data suggest high variability in CO 2 and CH 4 flux over time and space, therefore making it difficult to determine whether tundra ecosystems are currently acting as a net source or sink for atmospheric GHGs (McGuire et al 2009). We compared regionally aggregated flux estimates from our study to a recent synthesis of observations at tundra sites across the Pan-Arctic (McGuire et al 2012). We aggregated our simulation results by calculating the mean and standard deviation (as a measure of spatial variability) in flux estimates across all tundra cohorts in each sub-region, which is then comparable to the mean and standard deviation in observed fluxes across all sites in a sub-region. Over the 1990-2006 time period, the transient simulation (S T ) estimates CO 2 flux to be −7.2 ± 31.5 gC m −2 yr −1 (negative value denotes a sink) over all tundra areas as compared to −11.3 ± 21.0 gC m −2 yr −1 based on observations (figure 8). The reference simulation without ALT dynamics (S R ) results in CO 2 flux estimates that tend to overestimate the sink (−17.3 ± 39.2 gC m −2 yr −1 overall) and diverge further from the observations. The S T generally underestimates the regional-scale CH 4 source (1.5 ± 3.1 gC m −2 yr −1 ) for all tundra vegetation as compared to the observations (6.9 ± 6.8 gC m −2 yr −1 ), although a larger source (2.5 ± 4.6 gC m −2 yr −1 ) is estimated for wet tundra types only. The larger source in the observations could be due in part to the data being mostly typical of wetland sites, which will have higher CH 4 emissions than surrounding uplands. At the Pan-Arctic scale, CH 4 emissions estimates from the TEM Figure 7. Model-data comparison of TEM S T (solid lines) against measured annual maximum active layer thickness (ALT) from CALM network sites (points) in (a) North America, continuous zone, dwarf shrub, heath, herbaceous, and graminoid tundra types; (b) Asian Russia, continuous zone, dwarf shrub, heath, herbaceous, and graminoid tundra types; (c) North America, discontinuous zone, spruce forest and tall shrub tundra types ('taiga'); and (d) Russia, continuous zone, larch forest and tall shrub tundra types ('taiga'). The bars on mean TEM estimates for each year represent one sigma spatial variation across all modeled cohorts in each category. The legend associated with the measurements refers to the site code identifying individual sites in the network, as labeled on the CALM website (www.gwu.edu/∼calm/data/). have been shown to be in good agreement with atmospheric inverse estimates (McGuire et al 2010).

Conclusions
The state and flux of each component of the permafrost carbon feedback vary over time and space across the high-latitude region, reflecting the complex interaction of the key system drivers including climate, geomorphology, hydrology, vegetation dynamics and the physical properties of the permafrost regime. Here we describe the results of a model experiment on active layer dynamics suggesting that climate-driven permafrost thaw is already having a significant impact on ecosystem carbon fluxes across the high latitudes. Although not discussed here, on top of this gradual, widespread 'press' disturbance of permafrost thaw are the impacts of more rapid, localized 'pulse' disturbances such as fire and thermokarst that affect landscape-vegetation-hydrology dynamics, accelerate soil carbon mobilization, and determine the form of atmospheric greenhouse gas release as CO 2 or CH 4 (Grosse et al 2011). Ultimately, the impact of the permafrost carbon feedback on future climate change is largely determined by how these components interact as the net exchange of greenhouse gases between land and the atmosphere. Our results show that this feedback may have had a non-trivial impact on the Pan-Arctic carbon budget in recent decades, and studies agree on the likelihood of an increasing trend in permafrost C emissions for future decades (Koven et al 2011, Schaefer et al 2011. Taken together, this evidence points to the critical need for prognostic models to account for, and improve, their representation of permafrost carbon feedback components in order to make more reliable projections of future global climate change.