2.1 Simulations

CLM4.5 provides an accurate characterization of the physical and hydrological state of permafrost needed to evaluate permafrost vulnerability and identify key processes (Lawrence et al., 2008; Swenson et al., 2012). CLM4.5 includes a basic set of permafrost processes to allow projection of permafrost carbon–climate feedbacks, including snow schemes, vertically resolved SOM dynamics and soil hydrology, coupled hydraulic and thermal properties in frozen and unfrozen soils allowing realistic seasonal evolution of the active layer, and interaction with shallow (0–3 m) and deep (> 3 m) permafrost C (Swenson et al., 2012; Oleson, 2013; Koven et al., 2013, 2015; Lawrence et al., 2008). More abrupt thaw processes affecting permafrost C dynamics and talik formation such as thermokarst or other thaw-related landscape dynamics changes in wetland or lake distribution are not accounted for in CLM4.5 (see Riley et al., 2011, for more discussion).

CLM is spun up to C equilibrium for the year 1850 by repeatedly cycling through 20 years of pre-industrial climate forcing with CO₂ and N-deposition set at 1850 levels. C initialization is achieved via slow mixing by cryoturbation between the seasonally thawed active layers and deeper permafrost layers (Koven et al., 2009). Including vertically resolved processes leads to a sign change in the projected high-latitude C response to warming, from net C gains driven by increased vegetation productivity to net C losses from enhanced SOM decomposition (Koven et al., 2011). The soil grid includes 30 vertical levels that have a high-resolution exponential grid in the interval 0–0.5 m and fixed 20 cm layer thickness in the range of 0.5–3.5 m to maintain resolution through the base of the active layer and upper permafrost, and reverts to exponentially increasing layer thickness in the range 3.5–45 m to allow for large thermal inertia at depth. Soil C turnover in CLM4.5 is based on a vertical discretization of first-order multi-pool SOM dynamics (Koven et al., 2013; Oleson, 2013), where decomposition rates as a function of soil depth are controlled by a parameter Zτ (Koven et al., 2015; Lawrence et al., 2015). This depth control of decomposition represents the net impacts of unresolved depth-dependent processes. In this study, we utilize Zτ = 10 m, which yields a weak additional depth dependence of decomposition beyond the environmental controls and, as discussed and evaluated relative to Zτ = 1 m and Zτ = 0.5 m in Koven et al. (2015), results in CLM permafrost-domain soil C stocks that are in closest agreement (1582 Pg for Zτ = 0.5 m, 1331 Pg for Zτ = 1 m, and 1032 Pg for Zτ = 10 m) with observed estimates (1060 Pg C to 3 m depth; Hugelius et al., 2013). This reduction in initial C is due to higher decomposition rates at depth during the model initialization period. There is no C below 3.5 m, so additional thaw below 3 m has a small impact on the C cycle. We note that the relationship applied in CLM4.5, which implies multiplicative impacts of limitations to decomposition, is commonly applied in land biogeochemical models but is quite uncertain.

We use CLM4.5 configured as described in two recent permafrost studies (Lawrence et al., 2015; Koven et al., 2015) using time-varying meteorology, N deposition, CO₂ concentration, and land use change to capture physiological (i.e., CO₂ fertilization) and climate effects of increasing CO₂ over the period 2006–2300. We use an anomaly forcing method to repeatedly force CLM4.5 with observed meteorological from the CRUNCEP dataset for the period 1996–2005 (data available at https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm4.CRUNCEP.v4.html) and monthly anomalies added based on a single ensemble member from a CCSM4 Representative Concentration Pathway 8.5 (RCP8.5) simulation for the years 2006–2100 and Extended Concentration Pathway 8.5 (ECP8.5) for the years 2100–2300. The period from 1996 to 2015 represents a base climatological period used for calculating monthly anomalies, with a 20-year record chosen to minimize large anomalies in the first few years. This process is repeated for all variables and all times from 2006 to 2300 (constantly cycling through the same 1996–2005 observed data). Land air temperature for the period 2006–2300, shown in Fig. 1a, is projected to increase steadily over our simulation, with a slight decrease in the rate of warming.

We caution that we are using only a single ensemble member from CCSM4, and hence our results represent one realization from one model forced with one climate scenario. This results in uncertainties from the historical climate–weather forcing, the structure and parameterization of the model, and climate scenarios (both across models and across emissions scenarios).

Simulations are carried out on a global domain at a grid resolution of 1.25^∘ longitude × 0.9375^∘ latitude and saved as monthly averages. Simulation output is collected into decadal averages from 2011 to 2300 (e.g., 2011–2020 averages for the 2010s, 2021–2030 for the 2020s). Our method to link C balance changes to permafrost thermal state relies on identifying the timing of two key processes: (1) talik formation and (2) C source transition. Talik formation represents a critical threshold of permafrost thaw. The C source transition represents a shift of ecosystem C balance from a neutral or weak C sink to a long-term source as C balance shifts to increasing dominance of C source processes including permafrost thaw and fires (Koven, Lawrence and Riley, 2015). Using the hypothesis that talik formation triggers a transition to long-term C sources, we quantify the extent of talik formation and rate of transition to C source once talik has formed in permafrost-affected NHL ecosystems.

Table 1Site information for long-term borehole temperature measurements along the East Siberian Transect for the period 1957–1990. The nine sites reported in this table, presented in a north-to-south order, meet the criteria of at least 1 year of valid soil temperature data (≥ 10 months per layer, ≥ 55 months across five layers). Talik is observed in four of nine sites, two of which is observed in the first year of valid reported data. Site-specific thaw trends are provided for sites with at least 6 years of valid data. Regional trends are calculated from all available data for three regional locations.

Download Print Version | Download XLSX

Following Koven et al. (2015), we define the timing of C source transition from net annual sink to net source as the first decade when annual net biome production (NBP) decreases below −25 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{y}}^{-\mathrm{1}}$ and remains a source (NBP < 0 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{y}}^{-\mathrm{1}}$) through 2300. Here, we use the sign convention of NBP < 0 to represent net C flux from land to atmosphere (e.g., source). The timing of talik formation is defined as the first decade when soil temperature (T_s) for any layer between 0 and 40 m exceeds −0.5 ^∘C for all months in a calendar year (January–December), assuming that soils start off as permafrost at the beginning of our simulations in 2006. We use a negative freezing point threshold to account for availability of liquid water below 0 ^∘C due to freezing point depression. We note that the real threshold temperature at which liquid water remains available varies depending on the soil salinity or mineral content, the latter effect of which is included in the actual respiration calculations used by CLM. Here we use −0.5 ^∘C as the freeze–thaw cutoff and examine cutoffs at 0.5 ^∘C increments from 0 to −2.0 ^∘C.

We introduce the thawed volume–time integral, or “thaw volume”, as a metric to better understand thaw dynamics and help identify thaw instability thresholds. We integrate permafrost in both time (month of year) and depth (soil layer from the surface to 40 m) into a logical function that is one for thawed layers (T_s > −0.5 ^∘C), zero for frozen layers, and multiply each thawed layer by layer thickness to convert to units of meter months. This conversion accounts for nonuniform layer thicknesses, providing a consistent metric for comparing simulated and observed thaw.

Our analysis focuses on NHL grid points within the ABZ north of 55^∘ N. We analyze talik formation and C source transitions in the context of the simulated initial state of SOM as well as published maps of permafrost conditions from NSIDC (https://nsidc.org/data/docs/fgdc/ggd318_map_circumarctic/) and described in Brown et al. (2001). Permafrost extent is classified as continuous (90–100 %), discontinuous (50–90 %), and sporadic (10–50 %).

2.2 Observations

We compare simulated patterns of active layer dynamics and soil thaw to patterns observed from contemporary and historical borehole measurements of permafrost temperature profiles. We focus on sites in western North America and eastern Siberia with daily continuous observations year-round (January–December) over multiple consecutive years. The primary focus of data in North America (2004–2013) is to evaluate seasonal progression of soil thaw and talik formation near the surface (0–3 m). Siberian data, which have a longer record on average (1950–1994), are used to evaluate long-term trends in soil thaw at 0.0–3.6 m depth. Site locations are shown in Fig. 1.

Table 2Site information for borehole temperature measurements at three sites along a north-to-south transect in North America for the period 2004–2012. Climatological soil thermal states presented on a site-to-site basis in Fig. 4 are based on all available valid monthly data for each site, with valid data requiring at least 20 days of reported data for each layer. Layer of Deepest thaw represents the deepest layer in which mean soil temperature exceeds freezing (> −0.5 ^∘C) in at least 1 month. Month of latest thaw represents the latest month in which mean soil temperature exceeds freezing. Here, we define May as the earliest possible month and April as the latest possible month.

Download Print Version | Download XLSX

Siberian data are based on measurements along the East Siberian Transect (EST) (https://arcticdata.io/metacat/metacat/doi:10.5065/D6Z036BQ/default). The EST consists of 13 sites that cover a southwest-to-northeast transect in east Siberia (60.7^∘ N, 114.9^∘ E to 68.3^∘ N, 145^∘ E) during the period 1882–1994 (Romanovsky et al., 2007). For this study, we focus on the nine sites which report measurements as monthly averages at regular depths of 0.2, 0.4, 0.8, 1.6, and 3.2 m. Unfortunately, data gaps of years to decades exist on a site-by-site basis, and many years do not report the full annual cycle over multiple layers. To assess observed thaw trends from 1955 to 1990, we analyze individual sites which report at least 10 months yr⁻¹ of reported monthly mean soil temperature at each layer, and 55 months across the 5 layers (out of 60 possible layer months per year). Based on these requirements, we find that six of nine sites yield at least 6 years of data over multiple decades and are well suited for examining historical thaw trends. For comparison of observed trends to historical and projected trends from 1950 to 2300, we analyze clusters of sites by combining the nine sites into three groups based on approximate locations and calculate observed trends using the inter-site average at each location. We use two sites in northern Siberia (67^∘ N, 144^∘ E), six sites in southwest Siberia (61^∘ N, 115^∘ E), and one site in southeast Siberia (59^∘ N, 131^∘ E). Site information is shown in more detail in Table 1.

North American transect data are taken from the global terrestrial network for permafrost (GTNP) borehole database (http://gtnpdatabase.org/boreholes): (1) Borehole 1108 at Mould Bay in Canada (119^∘ W, 76^∘ N) from 2004 to 2012, (2) Borehole 33 in Barrow along the northern coast of Alaska (156^∘ W, 71.3^∘ N) from 2006 to 2013; and (3) Borehole 848 in Gakona in southeastern Alaska (145^∘ W, 62.39^∘ N) from 2009 to 2013. Mould Bay is a continuous permafrost tundra site with measurements at 63 depths from 0–3 m. Mould Bay has almost no organic layer (about 2 cm) and then sandy silt with high thermal conductivity. Barrow is a continuous permafrost tundra site with measurements at 35 depths from 0 to 15 m. The soil at Barrow is represented by silt with a bit of mix with some organics and almost no organic layer on top. Conductivity of the upper layer is ∼ 1 W mK⁻¹ for unfrozen and ≥ 2 W mK⁻¹ for frozen soil. Gakona is a continuous permafrost forest tundra site with measurements at 36 depths from 0 to 30 m. Gakona has a thick organic layer of moss (0 to 5 cm), dead moss (from 5 to 13 cm), and peat (from 13 to 50 cm), then silty clay at depth.

All North American transect datasets are reported as daily averages. For each site, we aggregate from daily to monthly averages requiring at least 20 days month⁻¹ at each layer and for each year. Measurements are reported at multiple depths and high vertical resolution (up to 0.1 m in shallow layers) but are generally nonuniform in depth (multiple layers missing, different layers reported for each site). Given these inconsistencies and records ≤ 8 years, we use these data for qualitative analysis of seasonal and vertical patterns in permafrost thaw. Site information and soil characteristics are summarized in Table 2.

3.1 Simulated talik onset in the 21st century

Our simulations show widespread talik formation throughout Siberia and northern North America over the period 2010–2300 (Fig. 1b), impacting ∼ 14.5 million km² of land in NHLs (55–80^∘ N) assuming a freeze–thaw threshold of −0.5 ^∘C and RCP8.5 and ECP8.5 warming scenarios. In Europe, southwest Asia, and North America (below 60^∘ N), 10.6 million km² of land either formed talik prior to the start of our simulation in 2010 in regions already experiencing degraded permafrost (e.g., Fig. 1d, permafrost extent < 90 % in southwest Siberia and southern North America) or did not have permafrost to begin with. A small amount of land along northern coastal regions (∼ 1.6 million km²) shows no talik formation prior to 2300.

The long-term trend and decadal variability of talik formation are quantitatively and qualitatively similar for freeze–thaw thresholds at or below −0.5 ^∘C (Fig. 1a). Peak formation generally occurs over the period 2050–2150, accelerating rapidly early in the 21st century, and leveling off in the late 22nd century. The timing and location of talik formation correlates with the annual mean temperature of permafrost at 3 m (T_soil-3 m) (Fig. 1c) and observed permafrost state (Fig. 1d, from Brown et al., 2001) at the start of our simulation; we see earlier talik formations in sub-Arctic regions (< 66^∘ N) with warm simulated permafrost (T_soil-3 m > 0 ^∘C) and permafrost extent less than 90 % and later formation in northern regions with cold permafrost (T_soil-3 m < 0 ^∘C) and continuous permafrost. Talik formation progresses northward from the sub-Arctic to the Arctic over time, starting in the warm, discontinuous permafrost zone in the 21st century then to the cold, continuous permafrost zone the 22nd century. This suggests a shift in permafrost state across the pan-Arctic from continuous to discontinuous over the next 2 centuries.

Figure 2Patterns showing the progression of soil thaw in the decades surrounding talik onset. Individual lines represent averages across the subset of talik-forming regions for each decade from the 2050s (darkest red) to the 2250s (darkest blue). (a) Integrated soil thaw volume, where the vertical solid line represents the mean timing of initial thaw at depth and late into the cold season (January–April). Note that the upper limit to the thaw volume metric in panel (a) is an artifact of the arbitrary maximum soil depth of 45.1 m in CLM4.5. Other panels show (b) date of spring surface thaw in the uppermost layer, (c) annual maximum active layer thickness, and (d) annual subsurface drainage (solid) and volumetric soil moisture averaged over the soil column (dashed). Grey shaded areas show the standard deviation of results for individual talik formation decades. Mean behavior exhibits a characteristic pattern: gradual increase in thaw volume and active layer depth prior to talik onset, abrupt shift in thaw volume, and active layer depth, followed by stabilization to constant thaw volume as soil drying and subsurface drainage increases.

Download

Our simulations demonstrate consistent patterns of changing thaw volume leading up to and following initial talik formation, independent of the decade of talik onset. Time series of thaw volume as a function of decade relative to talik onset (Fig. 2a) show a steady rise in thaw volume of 1–2 m month yr⁻¹ in the decades prior to talik formation, with thaw limited primarily to shallow soils (< 1.5 m) and summer–early fall. Thaw volume accelerates to 10–20 m month yr⁻¹ within 1–4 decades of talik onset, coinciding with thaw penetration at depth (∼ 2 m on average, Fig. 2b) and deeper into the cold season (∼ January–April). Thaw penetration into the January–April period occurs for the first time at 2.6 ± 0.9 decades prior to talik onset (vertical grey lines in Fig. 2a). At talik onset, thaw volume jumps from mean values of 60 ± 10.7 to 377 ± 44 m month yr⁻¹ at a mean depth of 4.1 m. Thaw volume levels out within 1 decade following initial talik formation and accelerated thaw of all soil layers; this leveling is an artifact of the maximum depth of soils in CLM4.5 (equal to 45.1 m) and represents the complete transition from permafrost to seasonally frozen ground in the model. The transition to deep cold-season thaw and rapidly increasing thaw volume represent key threshold signaling imminent talik onset.

Figure 3Evolution of simulated decadal thermal and hydrological state as functions of month and depth averaged across talik-forming regions in the 2090s. Each panel presents decadal average seasonal profiles in the decades surrounding talik onset from the 2050s (a) to the 2130s (i). Contours are soil temperature in 0.5 ^∘C intervals, with solid (dashed) lines denoting temperature above (below) a freeze–thaw threshold of −0.5 ^∘C. Stars indicate “thaw” months where soil temperature exceeds −0.5 ^∘C. Color shading is volumetric soil moisture anomalies relative to the 2040s, where red indicates drying. Note that soil depth on the y axis is plotted on a nonlinear scale. The soil thaw profile exhibits a shift from predominantly frozen and wet to perpetually thawed and drying conditions at depth while remaining seasonally frozen near the surface.

Download

Onset of surface thaw in the uppermost soils during the spring freeze–thaw transition provides another reliable predictor for talik onset. In particular, we find consistent dates and trends of spring thaw in the surface soil layer in the decades leading up to talik onset (Fig. 2c), shifting by about 1 week over 4 decades from day of year (DOY) 134 ± 2.8 (∼ mid-May) to DOY 127 ± 3.5 during talik formation (∼ early May).

Changes in total column soil water and subsurface drainage following talik onset may provide clues a posteriori that talik is already present. Lawrence et al. (2015) show that deepening of the active layer and thawing of permafrost allows water to drain deeper into the soil column, which dries out near surface soils. Our simulations show a similar, but very slight, drying pattern in shallow layers in the 4 decades prior to talik onset (1.3 % loss of soil moisture over 0–1 m depth; Fig. 2d), accounting for about half of total water storage loss in the column. More significant changes in water balance occur following talik onset, including more rapid drying in shallow layers (∼ 10 % over 4 decades) and in the column (∼ 16 %), and a substantial increase in subsurface drainage, as discussed below.

Figure 4Observed and simulated early 21st century soil thermal state as a function of month and depth for the North American Transect boreholes (black circles, Fig. 1). (a–c) Observed multi-year means for Mould Bay, Canada (2004–2012), Barrow, Alaska (2006–2013), and Gakona, Alaska (2009–2013). The color scale shows the mean temperature and the stars mark the months when each layer is thawed (T > −0.5 ^∘C). Simulated soil thermal state from 2006 to 2010 for borehole locations (d–f) and regions with 3 m permafrost temperature within 0.5 ^∘C of observed (g–i) show similar north-to-south spatial gradient to observations, especially for similar permafrost temperature. Note that the thaw state at Gakona, Alaska, persists at depths of 1–3 m into the deep cold season (January–February), perhaps signaling the threshold for rapid talik formation (see Fig. 3d).

Download

The time evolution of soil vertical thermal and hydrological structure for the subset of grid cells that form talik in the 2090s is shown in more detail in Fig. 3. Here, we have subtracted the thermal and hydrological profiles in the 2040s to show relative change. The 4 decades prior to talik onset are shown in Fig. 3a–d (2050s–2080s), the decade of talik onset in Fig. 3e (2090s), and the 4 decades following talik onset in Fig. 3f–i (2100s–2130s). CLM4.5 represents the process of soil thawing as passage of a “thaw front” in space and time through soil layers, penetrating and warming colder, deeper layers, and bringing the frozen soil environment at depth closer to thermodynamic equilibrium with the warming atmosphere. At 4 decades prior to talik onset (Fig. 3a), our simulated thawed layer exhibits a tilted time–depth profile with earlier thaw and longer thaw duration (∼ 4–5 months) in the near surface (< 1 m) compared to later thaw and reduced thaw duration (1–2 months) at maximum thaw depth (∼ 2 m). In the 3 decades leading up to talik onset, we find gradual deepening of the thawed layer to 3–4 m and penetration of thaw period into January–February.

Our simulations indicate an increased rate of heat transfer and thawing at depth following talik onset, leading to rapid subsequent thawing, drying, and decrease in the thickness of the seasonally frozen layer above talik (Fig. 3e–i). This rapid thawing is depicted in Fig. 2a as the large jump in thaw volume, and in Fig. 2d as enhanced drying and drainage, with drying peaking at 3.5–4.5 m depth. In our simulations, talik onset effectively pulls the “bath plug” that was the ice-filled pore space at depth, with year-round ice-free conditions allowing soil water to percolate and be diverted to subsurface drainage (Lawrence et al., 2015). We note that bedrock soil is not hydrologically active in CLM4.5, and thus the rate of thawing and drainage in response to permafrost thaw may be overestimated in deeper CLM4.5 layers near bedrock due to reduced heat capacity.

Figure 5Soil thaw observation time series from borehole measurements of soil temperature at sites along the East Siberian Transect over various periods from 1957 to 1990. Site coordinates are provided in the legend and plotted as crosses on the map provided in Fig. 1. Thaw trends are derived from estimates of thawed volume over a depth of 3.2 m for sites with > 55 months of data over multiple decades: Drughina, Lensk, Macha, Uchur, and Chaingda. Trend values are reported in Table 1. Vertical dashed lines mark the onset of talik formation at Lensk (1957), Macha (1970), Uchur (1974), and Chaingda (1989). Sites in southern Siberia show significant negative thaw volume trends over the 20th century, representing net increases in soil thaw. The trend at Drughina is not statistically significant, indicating that soil thaw is unchanged in northern Siberia.

Download

Our simulated pattern of phase lag for heat transfer to depth mimics observed thaw profiles in North America (Fig. 4), which are sensitive to latitude and ecosystem, but with more “vertical” time–depth tilt in CLM4.5 compared to observations. Borehole data show shallow (∼ 0.5 m) and seasonally short (∼ 3–4 months from June to September) thaw at the northernmost tundra site in the Canadian Archipelago (Fig. 4a; 76^∘ N, Mould Bay), shallow but longer thaw (5 months from June to October) moving slightly south to Alaska North Slope (Fig. 4b; 71.3^∘ N, Barrow), and deep (∼ 3 m) and seasonally long (May–February) thaw at the low-latitude continental boreal site in southeastern Alaska (Fig. 4c; 62.4^∘ N, Gakona). CLM4.5 shows reduced depth and seasonal duration of thaw when sampled at these specific geographical points, although the north–south gradient of increasing thaw moving south is preserved (Fig. 4d–f). Given the challenging task of comparing point locations with grid cell means, we also examine the mean behavior of CLM4.5 at locations where soil temperature at depth is similar to that observed. Accounting for permafrost temperature at 3 m (by sampling all locations with T_soil-3 m within 0.5 ^∘C of the observed temperature) better reproduces thaw depth, but with reduced seasonal duration throughout the soil column (Fig. 4g–i). These results suggest the current ensemble CLM4.5 run overestimates the rate of soil refreeze in early fall.

Figure 6Comparison of observed soil thaw to historical and future simulations at sites along the East Siberian Transect (crosses in Fig. 1). Observed thaw (filled circles) from 1955 to 1990 is based on soil thaw data in Fig. 5 and on the inter-site average at three locations: northern Siberia (blue), southwestern Siberia (yellow), and southeastern Siberia (brown). Simulated thaw from 1950 to 2200 is derived from CLM4.5 and sampled at the nearest grid cell of three above locations. Asterisks show simulated talik onset. Observed and simulated thaw trends are derived from soil thaw volume and estimated over the same period 1955–1990. We note a key discrepancy between observed and simulated thaw volume: simulated thaw volume is integrated over depths from 0 to 40 m; observed thaw volume is integrated from 0 to 3.6 m. The effect of this selection bias is a potential low bias in observed thaw volume. In general, soil thaw is projected to remain stable in northern Siberia but become increasingly unstable in southern Siberia.

Download

Based on the pattern of January and February freeze–thaw dynamics observed at Gakona in the 2010s and the time lag of 1–3 decades from this occurrence to talik onset in our simulations, we project that Gakona will form talik as early as the 2020s, assuming the atmosphere continues to warm as prescribed in CLM4.5. Talik onset in CLM4.5 is variable in the region containing Gakona (southeastern Alaska) with earliest onset by mid-century (∼ 2050s, Fig. 1a); however, our comparison to borehole temperatures at Gakona suggests that simulated thaw rates in southwestern Alaska and across pan-Arctic regions with similar permafrost temperatures are underestimated and that earliest onset may occur sooner than predicted. Overall, we find that simulated patterns of permafrost thermal state change are consistent with available observations but that the exact thaw rates are uncertain. Although there are many possible explanations for differences in observed and simulated thaw rates, we can attribute high observed thaw rates in part to a combination of (1) relatively dry upper soil at Gakona and Mould Bay and (2) low surface organic layer and high conductivity of the Barrow and Mould Bay soils. We keep these uncertainties in mind as we examine patterns of change and talik formation simulated into 2300.

3.2 Evaluation of simulated thaw rates and talik onset against Siberian borehole data

The Siberian borehole locations have similar permafrost extent (> 50 %) to the North American locations according to the Circumpolar Permafrost Map (Brown, 2001) and similar mean annual air temperature (∼ −13.6 ^∘C) in the 2000s according to CLM4.5. However, air temperature is more seasonal in Siberia, including colder winters (4 ^∘C colder) and warmer summers (6 ^∘C warmer). Spring thaw for the Siberian sites occurs 2 weeks earlier on average than for the North American sites in the 2000s, but follows the same pattern of later thaw date moving north along the borehole transect.

Next we examine thaw trends observed from borehole soil temperature data in Siberia in the 20th century and evaluate patterns of CLM4.5 projected trends in the 21st century. We note several caveats in these comparisons: (1) model simulations are based on only one realization (i.e., model ensemble member) of historic and future warming and projected permafrost thaw, (2) availability and access of long-term records in Siberia is limited, and (3) there is significant variability in space and time in simulated and observed thaw rates, making direct comparisons challenging. These comparisons thus serve primarily as a first benchmark for future model analysis and development.

We focus first on site-specific long-term historical trends by analyzing the six Siberian borehole sites which recorded at least 55 months and 5 years of temperature data spanning multiple decades: Drughina, Lensk, Macha, Oimyakon, Uchur, and Chaingda. Records at these locations show an increase in thaw volume with an average positive trend of 0.19 m month yr⁻¹ from 1955 to 1990 (Table 1, Fig. 5). All sites except Drughina show positive trends, with larger trends in southern locations, ranging from 0.51 m month yr⁻¹ from 1957 to 1990 at Chaingda in southern Siberia, to a statistically insignificant trend of −0.083 months yr⁻¹ from 1969 to 1990 at Drughina in northeastern Siberia, suggesting a more or less constant thermal state at this site. Further examination indicates that active layer thickness at Drughina actually decreased to 0.8 m from 1989 to 1990 compared to 1.2 m in the 1970s (data not shown). Drughina also shows smaller average thaw volume magnitude compared to other sites, consistent with shallower thaw. Together, these findings indicate that active layer thickness is decreasing at Drughina.

There is considerable spatial variability in thaw volume and trends, but in general thaw trends increase from west (0.18 m month yr⁻¹) to east (0.51 m month yr⁻¹). Talik forms at several sites, at different times between 1957 and 1990 (shown by vertical dashed lines on Fig. 5), with earlier talik to the west consistent with higher mean initial thaw volumes. We acknowledge the difficulty in identifying talik onset due to discontinuities in the dataset and limited vertical information; however, we note that the 15–30-year gap between talik formation in the western site cluster vs. Chaingda 15^∘ east is geographically consistent with model simulations of later talik formation in eastern Siberia in the 21st century (Fig. 1b) and thus may represent a gradual expansion of warming into the east. In general, permafrost appears to be degrading more rapidly at the southern locations compared to the northern location.

We recompute observed thaw trends at regional clusters using combined records at the two sites in northern Siberia (blue), six sites in southwest Siberia (yellow), and one site in southeast Siberia (brown, Table 1) and compare to historical and projected thaw volume trends in CLM4.5 (Fig. 6). Northern locations show a consistent pattern of low thaw volume (< 10 m month yr⁻¹) and negligible thaw trend (∼ 0 m month yr⁻¹) in the historical simulations and observed record from 1950 to 2000. Thaw projections in northern Siberia indicate an unchanged trend and continued stability of permafrost through the early 22st century, followed by a shift to accelerated soil thaw in the early 2120, marked by onset of deep soil thaw late in the cold season.

Southern locations show a systematic underestimate of mean thaw volume (< 20 m month yr⁻¹) compared to observations (∼ 40 m month yr⁻¹) from 1950 to 2000. Simulated thaw trends are negligible prior to 2000, but these likely represent an underestimate given low simulated thaw volumes and significant positive observed trends in both southeast and southwest Siberia beginning in the 1960s following talik onset (Fig. 5). Thaw projections show more abrupt shifts in thaw volume in the early 21st century in the southwest (∼ 2025) and in the mid-21st century (∼ 2050) in the southeast. The strong discrepancy between observed and simulated thaw and talik onset in southern Siberia warrants close monitoring and continued investigation of this region through sustained borehole measurements and additional model realizations of potential future warming.

3.3 Carbon cycle responses to changing ground thermal regime

Figure 7Projected decade when permafrost regions shift to long-term C sources over the period 2010–2300, and relation to talik onset, soil C, and fire emissions. (a) Map of the decade of transition to C source, reflected in the color code, showing earlier transitions in cold northern permafrost. (b) The area of land that transitions peaks in the late 21st century and is driven by regions where the C source leads talik onset (dashed). (c) The decadal time lag from talik onset to C source transition shows positive lags in warm southern permafrost (C source lags talik) and negative lags in cold norther permafrost (C source leads talik). (d) Histogram shows trimodal distribution of permafrost area as a function of decadal time lag, with negative lags related to high soil organic matter (green bars and map in e) and large positive lags related to fires (red bars and map in f) but delayed by high productivity. See text for details. These results assume a Representative Pathway 8.5 warming scenario through 2100 and an Extended Concentration Pathway 8.5 through 2300.

Download

Figure 7a plots the decade in which NHL ecosystems are projected to transition to long-term C sources over the next 3 centuries (2010–2300). A total of 6.8 million km² of land is projected to transition, peaking in the late 21st century, with most regions transitioning prior to 2150 (4.8 million km² or 70 %; Fig. 7b, solid black). C source transitions which occur in the permafrost zone, accounting for 6.2 million km² of land (91 % of all C source transitions), also form talik at some time from 2006 to 2300 (Fig. 7c). The remaining C source transitions (0.6 million km², or 9 %) occur outside the permafrost zone, primarily in eastern Europe.

Figure 8Cumulative net biome production (NBP) over northern high-latitude (NHL) regions (> 55^∘ N) from 2010 to 2300. NBP < 0 represents a net C source. NHL regions are divided into the following categories: all NHL land (diamonds), NHL land regions which form talik from 2010 to 2300 (crosses), and regions which transition to long-term C sources from 2010 to 2300 (black solid). C source transition regions are further broken down based on the lag relationship between talik onset and C source transition as follows: regions where the C source transition lags talik onset (dotted), leads talik onset (dashed), and occurs in the absence of talik (dashed dotted). C source transition regions also divided by soil C content and fire activity: regions where soil organic matter (SOM) exceeds 100 kg C m⁻² (green), fire emissions exceed 25 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (red), and SOM and fires do not exceed these thresholds (blue). Regions which transition to C sources prior to talik formation make up half of the total C source area but account for most of the cumulative C source (∼ 80 %) in large part due to high soil C.

Download

Net C emissions from C source transition regions are a substantial fraction of the total NHL C budget over the next 3 centuries (Fig. 8). The cumulative pan-Arctic C source increases slowly over the 21st century, reaching 10 Pg C by 2100 with RCP8.5 warming, then increases more rapidly to 70 Pg C by 2200 and 120 Pg by 2300 with sustained ECP8.5 warming (Fig. 8, solid black). This pan-Arctic source represents 86 % of cumulative emissions in 2300 from the larger NHL talik region (crosses), despite the 2-fold smaller land area, and exceeds the talik region through 2200 due to mitigating widespread vegetation C gains (Koven et al., 2015). Cumulative emissions over all NHL land regions (diamonds, > 55^∘ N) increase in similar fashion to the talik region, reaching 120 Pg C by 2200 and 220 by 2300, with no sign of slowing.

The geographic pattern of C sink-to-source transition date is reversed compared to that of talik formation, with earlier transitions at higher latitudes (the processes driving these patterns are discussed in detail below). Overall, the lag relationship between talik onset and C source transition exhibits a trimodal distribution (Fig. 7d), with peaks at negative time lag (C source leads talik onset, median lag = −5 to −6 decades), neutral time lag (C source synchronized with talik onset; median lag = −2 to 1 decade), and positive time lag (C source lags talik; median lag = 12 decades; red shading in Fig. 7c), each of which is associated with a distinct process based on soil C and fire emissions as discussed below. Roughly half of these regions (3.2 million km²) show neutral or positive time lag (lag ≥ 0). This pattern, characteristic of the sub-Arctic (< 65^∘ N), represents the vast majority of C source transitions after 2150 (Fig. 7b, dotted), but only accounts for 17 % of cumulative emissions (20 Pg C by 2300; Fig. 8, dotted). The remaining regions (3.0 million km²) in the Arctic and high Arctic (> 65^∘ N) show negative time lag and account for most of late 21st century sources and cumulative emissions (95 Pg C by 2300, or 79 %; Fig. 8, dashed). C sources in regions not identified as talik (0.63 million km²) either show talik presence at the start of our simulation or are projected to transition in the absence of permafrost or in regions of severely degraded permafrost (Fig. 7c, dash dotted). This region contributes only 5 Pg C (4 %) of cumulative C emissions in 2300.

Figure 9Net biome production (NBP) as a function of thaw volume. Symbols represent NBP and thaw volume values averaged over regions which transition to long-term C source from 2060 to 2140, binned into regions where the decade of C source transition (a) leads talik onset, (b) lags talik onset, and (c) lags talik onset AND where fires exceed 25 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. Colors indicate decade relative to C source transition, denoted by the large green marker, which occurs when NBP exceeds −25 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (grey horizontal dashed line). The grey square marker indicates the mean NBP and thaw volume values during talik onset. Cases where C source leads talik (a) show small thaw volumes during C source transition and amplified C sources during talik onset. Cases where C source lags talik (b)–(c) show large thaw volumes during C source transition and C sinks during talik onset.

Download

Here, we investigate biological and soil thermal processes driving these relationships, focusing first on regions where C source transition leads talik onset (blue shading in Fig. 7c). In these regions, thaw volume is low (< 50 m month yr⁻¹) and shows a weak relationship to NBP (NBP decreases much faster than thaw volume) prior to C source onset (indicated by large green circle in Fig. 9a). By the time thaw volume reaches 300 m month yr⁻¹ and talik formation occurs, these regions are already very strong sources (NBP > 150 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). This suggests that C sources in these regions are not driven by respiration of old C from deep soil thaw, and thus alternative explanations are needed.

Figure 10Evolution of simulated soil thermal and hydrological state, plotted as a function of month and depth, for regions which transition to long-term C sources in the 2060s but do not form talik for another 3 decades (≥ 2090s). This represents cases where C source leads talik (e.g., Fig. 9b). Each panel presents decadal average seasonal profiles in the decades leading up to C source transition. Shading and contour details are explained in Fig. 3. These profiles exhibit shifts in thaw period (October), depth (> 1.5 m), and soil moisture (drying) in the transition decade.

Download

Closer examination of thermal and moisture dynamics in shallow soils reveals three potential indicators of C source transition: (1) seasonal duration of thaw, (2) depth of thaw, and (3) soil drying. For example, vertical profiles of soil temperature and moisture (Fig. 10) in regions which transition to C sources in the 2090s show deeper seasonal penetration of soil thaw, a jump in active layer growth, and enhanced year-round soil drying during the C source transition decade (Fig. 10d). A broader analysis of soil thaw statistics over all regions and periods indicates that most C source transitions (∼ 2.3 million km², or 77 % of land where C source leads talik) occur at active layer depths below 3 m and thaw season penetration into November.

Further examination of ecosystem biogeochemistry also shows high initial C stocks in these regions (red shading in Fig. 7e). The median initial state of SOM, 109 kg C m⁻², is nearly a factor of 2 larger than the median value in regions where C source lags talik onset (SOM = 59 kg C m⁻²). These regions also show 40 % less gross primary production (median GPP = 755 vs. 1296 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) and higher over saturation prior to C source onset (water-filled pore space at 0.5 m depth at 10, 5, and 2 decades prior = 0.63, 0.59, and 0.57 mm³ mm⁻³ for cold permafrost vs. a near constant value of 0.57 mm⁻³ in warm permafrost). The total area of land in which SOM exceeds 100 kg C m⁻² represents two-thirds of all land where C sources lead talik onset (2.0 million km²) and peaks at a negative time lag of −5 to −6 decades (Fig. 7d, green bars), which perfectly aligns with the peak distribution of negative time lags. Cumulative C emissions from regions of SOM > 100 kg C m⁻² are also two-thirds of total C emissions (80 Pg C; Fig. 8, green). These results indicate peat-like conditions characterized by saturated soils, high C stocks, and low annual productivity, which allow low thaw volumes (active layer depth < 2 m and peak thaw month of October, on average) and rapid soil drying to produce early C losses in colder environments in the absence of talik.

In regions where C source transitions lag talik onset (red shading in Fig. 7c), NBP is strongly sensitive to changes in thaw volume until C source onset occurs (Fig. 9b), and talik formation occurs when these regions are weak sinks (NBP > 0 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). In general, C source onset under high thaw volume indicates these regions are more sensitive to C emissions from deep soil thaw. However, as noted above, neutral and positive time lags show a bimodal distribution peaking near 0 and 15 decades, and thus additional explanations are needed. Further examination shows high fire activity in these regions at the time of C source onset (red shading in Fig. 7f). The regions where fire C emissions exceed 25 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, representing our threshold for C source transition, are exclusively boreal ecosystems, account for one-third of all land with negative lags (∼ 1.1 million km²) and align perfectly with the peak distribution of positive time lags (Fig. 7d, red bars) and cumulative C emissions (20 Pg C in 2300, Fig. 8, red). NBP is less sensitive to thaw volume in regions where fire dominates the C balance, which are strong C sinks at talik onset (Fig. 9c), where soil C respiration is 13 % less than non-fire regions (median SOMHR = 331 vs. 382 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), and productivity is 25 % more (median GPP = 1548 vs. 1216 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). Fire regions are also 28 % drier on average in the surface layer than non-fire regions (volumetric soil moisture = 0.28 vs. 0.39 mm⁻³ in summer (May–September) in the upper 10 cm of soil). These results suggest that soil thermal processes and talik formation are significant factors driving C source transition in regions with reduced productivity, but fire activity, spurred by soil drying, drives C source transition in higher productivity regions.

The decadal time lag between talik onset and C source transition is more normally distributed in the remaining region, represented by the residual grey bars visible in Fig. 7d, which occurs predominantly in cold northern permafrost in northwestern Siberia, where low SOM (< 100 kg m⁻²) and fire emission (< 25 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) prevail. This region has a mean lag of 1 decade from talik onset to C source, with high standard deviation of lags (±8 decades) reflecting a skewed distribution of GPP; low productivity in cold permafrost (GPP = 385 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) increases the likelihood that soil thaw will lead to C source transition prior to talik onset, and high productivity in warm permafrost (GPP = 1111 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) increases the likelihood of a transition after talik onset. Cumulative C emissions from this region are on the low end (27 Pg C by 2300; Fig. 8, blue) due to low soil C (SOM = 59 kg C m⁻²).

Figure 11Time series of ecosystem C fluxes showing seasonal and decadal patterns during C source transition. This presents mean and standard deviations over the period 2040–2270 for (a–b) gross primary production (GPP), (c–d) sum of respiration from soils (SOMHR) and litter (LITHR), (e–f) difference of respiration from soils and litter, and (g–h) net biome production (NBP) where NBP < 0 indicates source. The left columns show seasonal fluxes during the decade of C source transition. The right column shows the evolution of decadal mean fluxes in the 3 decades preceding and following C source transition. Regions where C source transition leads talik (blue) show similar patterns to regions where transition lags talik (red), most notably a jump in soil vs. litter respiration during C source transition (f) corresponding in time and magnitude to decreasing NBP (h). The primary difference between regions is the seasonal distribution of SOMHR vs. LITHR (e), which shows a large soil respiration source throughout the cold season in cases where C sources lag talik. This indicates an annual source of deep old C.

Download

Independent of the presence of talik, a key effect of an increasing number of thaw months is an increasing rate of respiration from soil C pools. Warming and CO₂ fertilization increase the rate of photosynthetic C uptake, increasing soil respiration mainly from younger near-surface C pools; whereas deeper thawing affects both young and old C pools, so that the depth of thaw dictates the timing and dominant C age of the net respiration flux. Figure 11 illustrates this with a comparison of decadal respiration trends for SOM (SOMHR) and litter (LITHR) C pools for C source transitions in the mid-21st century, for scenarios where C source leads talik onset (blue line, cold permafrost) and lags talik (red lines, warm permafrost). Here, we examine combined respiration (SOMHR + LITHR) and respiration difference (SOMHR − LITHR) from soil and litter C pools.

GPP and combined respiration increase by ∼ 15 % decade⁻¹ for each permafrost regime surrounding the decade of C source transition with peak fluxes in the growing season (Fig. 11a–d). Combined respiration in cold permafrost is systematically larger than in warm permafrost in the growing season (May–September) and smaller in the cold season (October–April). In particular, combined respiration is effectively zero for the late cold season (January–April) in cold permafrost and significantly positive in warm permafrost over the same period. The respiration difference also increases surrounding the C source transition (Fig. 11e–f), but with two key differences from combined respiration: (1) the decadal increase is exponential, starting from a value near zero just 3 decades prior to C source transition, and (2) peak respiration difference occurs in late summer and early fall. Because litter respiration in the model is mainly drawing from C pools with short turnover times, the litter respiration flux equilibrates rapidly to changes in productivity and thus its change primarily reflects changes to inputs rather than decomposition rates. Conversely, soil C pools, which have much longer turnover times, equilibrate much more slowly to the productivity changes and thus primarily reflect changes to the turnover times.

The trend in the respiration difference in warm and cold permafrost, which increase by similar amounts (∼ 100 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), thus reflects an increasing dominance of respiration from younger and older soil C pools, respectively. These trends are identical to the corresponding NBP trends, which decrease by 100 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ over the same period from neutral to net source (Fig. 11g–h), such that the differences between GPP and respiration driving the NBP trends are explained almost entirely by the increasing fraction of soil vs. litter respiration. Furthermore, warm permafrost shows sustained dominance of soil respiration during the entire cold season. These results are consistent with an increasing thaw effect on C budgets during C source transitions, but where shallow thaw of young soil C dominates in cold permafrost and where talik formation and deep thaw of old soil C dominate warm permafrost.

These results suggest that where talik forms, soil respiration increases throughout the year as talik and perennial thaw mobilize deeper old soil C to respiration. In the absence of talik in colder environments, soil respiration increases primarily in the non-frozen season due to increased availability of thawed shallow soil C. The lower GPP in colder regions suggests that increased availability of substrate for respiration due to plant growth and soil C accumulation has less impact on C source transition in our simulations than soil thaw dynamics and the initial state of soil C. Thus, cold permafrost locations become C sources due only to thaw-season dynamics while warmer permafrost locations transition to C sources due largely to changes in cold-season dynamics.