Thelon-Kazan Peatland (TKP) was selected as a site for an earlier modeling study using a climatic gradient through Western-Central Canada (Treat et al., 2021). In that study, we discovered divergent peat profiles at TKP depending on the parameterization for the hydrological landscape setting. Here, we are exploring those different patterns to develop a more general hypothesis or possible explanation about peat/permafrost dynamics in contrasting hydrological and landscape settings, rather than trying to recreate the C accumulation dynamics at TKP. TKP is located in the Canadian Arctic, approximately 400 km north of treeline in the continuous permafrost zone (66° 27.07′ N, 104° 50.08′ W), where the mean annual air temperature was −11.3°C and mean annual precipitation was 272 mm (Station Baker Lake A 1980–2010; https://climate.weather.gc.ca/climate_normals/; accessed 12 December 2019). The TKP peatland is a fen with high-centered polygonal features and an active layer thickness of 30–40 cm, where peat cores were described by Vardy et al. (2005); Vardy et al. (2000). The manually digitized data from TKP were included in an earlier synthesis (Treat et al., 2016a) and used in an earlier modeling effort (Treat et al., 2021). Three cores were taken at TKP using monolith extraction in the surface soil and motor-driven CRREL Permafrost augers at some point prior to 1998; the exact timing of core collection is unknown. Here we focus on core TK1P2 (also called TKP-2). The TK1P2 core was collected from the center of a high-center polygon and accumulated 193 cm of peat since 6,620 cal. BP (basal age: 5820 ± 70 ¹⁴C BP; WAT-3066). The peat core (TK1P2) was mainly composed of Cyperaceae with mosses, lichens, and dwarf shrubs in the surface 35 cm (Vardy et al., 2000; Vardy et al., 2005). Data available from the TKP cores includes bulk density, carbon or organic matter content, radiocarbon dates, and plant macrofossil analysis. The data from TKP are available for download (doi: 10.1594/PANGAEA.863697; Treat et al., 2016b); the dataset key (Variable: Auth. Site.CoreID) for TKP is VAR-TKP-01.

We used HPM-Arctic (Treat et al., 2021), which couples a carbon-hydrologic model for peatlands (Frolking et al., 2010) with the Geophysical Institute Permafrost Lab soil thermal model GIPL2 (Marchenko et al., 2008). HPM-Arctic simulates peat formation since initiation to present and into the future, also generating time series of active layer thickness and age of the peat at the bottom of the active layer. HPM-Arctic has been tested and implemented for a transect of peatland sites across western Canada, including TKP (Treat et al., 2021). The model code, parameters, and climate driver data are available from: https://doi.org/10.5281/zenodo.4647666; derivation of parameters and climate drivers is described below.

Briefly, HPM simulates the development of a peat profile over millennia, from initiation, using an annual litter cohort approach with three major plant types (shrubs, herbaceous taxa, and mosses). In HPM-Arctic, the changes in vegetation and NPP are driven primarily by changes in water table depth, following the original HPM formulation (Frolking et al., 2010). Changes in water table can be driven either by external climate forcings (precipitation vs. ET), or by internal peat dynamics, including peat accumulation. NPP temperature sensitivity was modeled as a Q₁₀ function, with a Q₁₀ value of 1.8, based on an empirical relationship between mean annual air temperatures and above-ground net primary productivity for mosses, vascular plants, and trees that was developed for a transect of peatland sites in boreal Manitoba, Canada (Camill et al., 2001). CO₂ concentration impacts on productivity are not accounted for in HPM.

Each year’s annual litter input is aggregated into a litter cohort and tracked as it becomes buried and decomposed through the millennia of the simulation. HPM tracks the fraction of total plant litter inputs (above- and below-ground) remaining in annual peat cohorts through the entire profile. This fraction modifies the initial peat litter quality, so that under the same environmental conditions, peat that has lost half its mass to decomposition will decompose at half the rate of fresh peat (Clymo, 1992; Frolking et al., 2001). Decomposition rates are dynamic and dependent on plant litter types (e.g., species-specific decomposition rates), water table or water content, and soil temperature. HPM simulates the water-filled pore space of the unsaturated peat water filled pore space, based on peat bulk density and distance above the water table. In that zone, there is an optimum WFPS for decomposition, with rates falling for drier and wetter peat. Full anoxia occurs some distance below the simulated water table for several reasons: the actual water table is more dynamic than monthly time step simulations, oxygen and other terminal electron acceptors can be delivered to the peat in groundwater flow (much greater for fen than bog), and precipitation water, and also diffuse into the deeper peat from the underlying subsoil or the atmosphere. HPM simplifies these complex phenomena by specifying an exponential decline in the anoxia multiplier on decomposition from its value at the water table to a minimum value, following the original HPM formulation (Frolking et al., 2010). This decline is parameterized to happen much more abruptly in bog than fen conditions.

When decomposing peat becomes frozen into permafrost, it stops decomposing until permafrost thaw, preserving the quality of the peat at freeze-up, which affects the rate of decomposition at thaw. Active layer thickness, updated annually, is determined by identifying the soil depth immediately above where the temperature remains below 0°C for 2 years continuously, in accordance to the definition of permafrost (Harris et al., 1988). HPM-Arctic uses a simple “old-new” carbon tracking algorithm, whereby after a specified year all moss, sedge, and shrub plant litter gets labeled as “new,” so that its accumulation as peat and loss through decomposition can be tracked separately from the older peat derived from plant litter inputs prior to the specified year. The year 2015 C.E. is defined in the model as both present-day and the boundary between new (future) and old (past) carbon inputs.

Some site-specific calibrations were done for several model parameters to capture variability (often not quantified) related to individual watershed and site characteristics (Table 1). Peat initiation often occurs in a local topographic low, receiving run-on from the surrounding watershed; as the peat accumulates and the peat surface rises, it can shift to a local topographic high point, and shed water (run off) rather than receive it (Charman, 2002). The site-specific model parameters include the accumulating peat height at which this shift from run-on to run-off occurred (H_run-on/off) along with a constant scaling factor, and the peat height when initial fen-type vegetation transitioned to bog-type vegetation (H_FBT). When the peat height exceeds the site-specific H_FBT, the peatland transitions from a fen to a bog, which involves a decrease in annual NPP to a varying degree (Rydin and Jeglum, 2006) modeled with a site specific fractional parameter (F_NPP-bog). With greater lateral hydrological inflow, and therefore a shorter water residence time in the saturated zone, fen conditions are assigned a longer scale length (persistence of partial anoxia below the water table) to a full anoxia impact on decomposition rate. This is modeled with an anoxia scale length parameter (Frolking et al., 2010), which controls how far/quickly reduced electron acceptors are replenished below the water table. This affects peat decomposition rates below the water table. Frozen peat decomposition rates are set to zero (also for seasonal winter frost), while decomposition persists in the seasonally-thawed surface active layer.

Site-specific parameter values were determined from a combination of observations and/or optimization routines (Table 1; Treat et al., 2021). The parameter H_FBT was determined by trial and error from the final peat height and observations of the height of the fen-to-bog transitions (H_FBT) in the site core profile; the model parameter H_FBT generally was higher than in the observations. H_run-on/off was determined by trial and error when the difference between the final modeled peat height and observed peat height reached <10 cm, as well as the macrofossil composition compared with the dominant modeled PFTs, which indicated relative water table position over time (e.g., dry or wet). The other three parameters (F_NPP-bog, and anoxia scale lengths for fen and bog) were determined from minimizing the root mean squared error between the observed and modeled age-depth profiles, where the age of the peat surface was assumed to be the year of sampling. At TKP, the parameterization used results from earlier model runs at Ennadai Lake, a polygonal permafrost peatland site in NW Territories (c.f. Treat et al., 2021) and generalized parameters for temperate and boreal peatlands (Frolking et al., 2010) due to the limited number of radiocarbon dates (e.g., Figure 2), which limited the number of tunable model parameters (under constraint problem). However, the study goal was not to find the closest agreement to a particular peat core at TKP, but to test a hypothesis about how local hydrology (wet/dry) in a permafrost setting can impact peat accumulation rates and carbon loss upon thaw.

For calculating organic matter stocks and C stocks, model output (mass of peat) was multiplied by mean values determined from a synthesis of over 10,000 peat layers spanning the pan-Arctic permafrost region (Treat et al., 2016a). The conversion factor from model peat mass to organic matter (OM) was 0.924 g OM g⁻¹ peat. The conversion factor from organic matter to carbon was 0.495 gC g⁻¹ OM. While this assumption does introduce some uncertainty into the model results for predictions of future carbon losses, these are likely not the major sources of uncertainty in predicting future C loss as other processes, including changes in landscape hydrology in permafrost regions, are less well known and have larger impacts (Quinton et al., 2010).

The two landscape scenarios at TKP used identical climate drivers, initiation times, vegetation parameters, and most site-specific model parameters (Table 1). We created two model scenarios, dry and wet landscape position, by varying a single site-specific model parameter, H_run-on/off (Table 1), which is the height of the peat when it begins to shed water, exerting a strong control on the water table level and soil moisture. The dry landscape scenario used H_run-on/off = 0.5 m, which meant that the peat height quickly crossed the threshold to generate runoff, i.e., peat accumulation caused the site to become a local high point in the landscape, and so the peat profile became relatively dry. The wet landscape scenario used H_run-on/off = 4.5 m, a peat height that was not reached during the simulation, ensuring a relatively wet peat profile throughout the simulation. While HPM-Arctic does not simulate a surface energy balance, both scenarios are driven by identical climate forcings and modeled active layer depths are consistent with recent meta-analysis results (Clayton et al., 2021). Secondly, an earlier study showed relatively small differences in net radiation among different tundra sites without standing water, relative to sites with standing water, within the same region (Rouse et al., 2000); mean annual water table in these two landscape scenarios is below the peat surface, meaning that there is no standing water and that the net radiation absorbed by the sites is similar.

We utilize the TraCE-21ka transient paleoclimate simulations (Liu et al., 2009: https: //www. earthsystemgrid. org/project/trace.html) to drive HPM-Arctic with monthly temperature and precipitation from 8,000 B.P. to 1990 CE (8 Kyr timeseries). The continuous TraCE-21ka simulations include time-dependent changes in atmospheric greenhouse gases, insolation, and paleogeography as sea level rises from the melting of the large northern hemisphere ice sheets. Monthly temperature and precipitation time series for TKP were extracted from the 3.75° × ∼3.75° climate model grid by bilinear interpolation and then bias-corrected using contemporary regional climate data. To obtain timeseries of temperature for input to HPM-Arctic, we computed anomalies of the 8 Kyr time series to the TraCE-21ka 1950–1990 CE mean and applied the anomalies additively to a modern observed gridded data set (CRU TS v3.32; Harris et al., 2014) using the same base period. The 8 kyr precipitation timeseries was derived by applying the ratio of the 8 Kyr precipitation to the 1950–1990 CE mean as a scaler to the CRU data. As with the TraCE-21ka data, the CRU grid (0.5° × 0.5°) was interpolated to our site locations.

To continue the 8 Kyr timeseries into the future, we adopt the CCSM4 RCP8.5 simulation from the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2012). The RCP8.5 scenario was chosen as an end-member or upper bound to capture the greatest projected changes in temperature. CCSM4 (Gent et al., 2011) is the successor to the CCSM3 model used in the TraCE-21ka simulations, and so was chosen for model consistency. While focusing on one model projection is a limitation, CCSM4 also has an transient climate response (the response to a 1% yr⁻¹ doubling of atmospheric CO₂) of 1.8°C, which is the same as the CMIP5 ensemble multi-model mean of 1.8°C (Flato et al., 2013), giving us confidence the projected temperature changes at our sites are not unreasonable relative to the full CMIP5 ensemble. To extract projected time series at the peatland site locations, we use the same process of bilinear interpolation from the CCSM4 grid (0.9° × 1.25°) and applying CRU bias correction using the 1950–1990 CE climatology period. In our analysis and results below, HPM-Arctic is driven by TraCE-21ka output prior to 1990 CE and CCSM4 afterwards. The RCP8.5 simulations span 2005 through 2100 CE. To capture the effects of new equilibrium climate conditions at 2100 on the relatively slow processes of peat accumulation and loss, we extended RCP8.5 from 2100 for an additional 200 years by randomly sampling full years from the final 25 years of the simulation (2076–2100 CE) using the R command “sample” (replace = T). The modeled climate drivers are shown in Supplementary Figure S1.

All statistical analyses were conducted using R statistical software (Team, 2008). Errors represent standard deviation. In order to compare modeled C losses with observations, we used the peat C stock in 2015 CE, and the results from old/new C tracker in 2100 CE and 2300 CE. Net C loss (or gain) was calculated from the difference in total C stocks in 2015 and 2100 CE or 2300 CE; new C additions and old C losses were calculated relative to the stocks in 2015.

The model of peat development for TKP under wet and dry landscape settings produced divergent peat water table and active layer depths (Figure 1), and vegetation species composition (Table 2) despite using the same climatic drivers (Supplementary Figure S1). Simulated peat depths between the wet and dry scenarios in the present day show that 1.9 m accumulated in the wet scenario and 2.8 m in the dry scenario (Figure 1; Table 2). The water table in the wet scenario remained near the surface (mean: 3 ± 2 cm depth) throughout the simulation (Figure 1A), while the water table in the dry scenario was deeper (mean: 31 ± 9 cm depth, Figure 1B). The timing of peat and permafrost development also differed between the wet and dry scenarios. Modeled active layer depths in the present day were <50 cm in the dry scenario and >80 cm in the wet scenario (Table 2). Present-day vegetation differed between the two scenarios, with shrubs having the highest NPP in the dry scenario and herbaceous species (e.g., sedges) having the highest NPP in the wet scenario (Table 2).

The RCP 8.5 driven HPM-Arctic simulations at TKP under wet and dry landscape settings for the future resulted in substantially different outcomes at 2100 CE for peat temperatures and permafrost, hydrology, vegetation, and carbon cycle, despite having the same climatic drivers. By 2100 CE, transient warming projected under RCP 8.5 increased mean annual air temperature from −11.3°C to −4.4°C. Despite the strong warming, permafrost remained in the peat profile (Figure 3). Warmer temperatures increased evapotranspiration, which combined with gradual thawing of permafrost increased the depths to the perched water table (Table 2). The changes in water table shifted the modeled species composition, while the warmer air temperature increased modeled NPP (Figure 4). Along with warming temperatures, the changes in water table, active layer thickness, increased decomposition as well as net primary productivity altered the peat profiles and C balance by 2100 (Figures 5, 6).

Under the assumption that 21st century warming stabilized between 2100 and 2300 CE, the peat temperatures at TKP differed little between 2100 and 2300 CE in the two scenarios, and were still cold enough to maintain permafrost (Figure 3B). However, active layer thickness decreased by roughly 10 cm between 2100 and 2300 in both scenarios (Table 2), as the decrease in peat height (5–30 cm) due to net peat loss (decomposition) exceeded the gradual thawing at the top of the permafrost peat (Figure 5). Still, an additional ∼10 cm of permafrost peat thawed in both scenarios between 2100 and 2300. Net C losses by 2300 were larger than in 2100 as the increase in NPP and net peat accumulation were smaller than the decomposition of the old carbon (Figure 6). These C losses represented a small to moderate fraction of total peat C stocks which is in rough agreement with other high Arctic cold and wet sites (Elberling et al., 2013). Overall, these findings agree with earlier studies that projected larger permafrost C losses at 2300 than 2100, and show that these net losses are dependent on both new C sequestration in the soil and old permafrost C losses (Figure 6).

In this modeling study, an important result was that permafrost formation and stability was strongly influenced by the parameter that represents the regional hydrology and the wetland connectivity to the watershed, H_run-on/off (Figure 1). Specifically, H_run-on/off represents the height of the peat when the peatland stops receiving run-on from the surrounding landscape and begins to shed water via run-off, as occurs in domed bogs (Charman, 2002). While the H_run-on/off parameter was important for maintaining a shallow water table level, for example representing wetlands that are still connected to the surrounding watershed (Woo and Young, 2006), it also controlled rates of peat accumulation through feedbacks between the water table position and the peat thermal properties, which influenced the thickness of the active layer and ultimately decomposition rates. Peat accumulation results from rates of productivity exceeding rates of decomposition. During the period of exceptional high peat accumulation rates between peat initiation in 5700 BP to 4550 BP (Figure 1D), NPP in the wet landscape was on average 50% greater than the dry landscape NPP (Figures 4A,C), yet Figure 1D shows that peat accumulation rates are >2.5 times larger for this dry landscape scenario. Since the differences in peat accumulation that can be attributed to differences in productivity would favor the wet scenario, the observed difference must be due to slower rates of decomposition.

In peatlands, hydrology controls the decomposition rates through two mechanisms: oxygen availability (the standard peat situation of saturation and slow decomposition) and peat thermodynamics. Peat thermodynamics controls soil temperature and therefore, permafrost formation and active layer thickness. Permafrost formation and active layer thickness, in turn, controls decomposition rates, which are very low to negligible in permafrost (Schädel et al., 2016). In the dry landscape scenario during rapid peat accumulation, these effects on thermodynamics (i.e., shallower active layer) are stronger than the effects of substantially higher aerobic decomposition rates occurring in the peat acrotelm.

In this study, vertical peat accumulation was ultimately governed by the residence time of peat in the active layer, or the length of time when peat can be decomposed before it is frozen into the permafrost (Figure 1). In these permafrost systems, the active layer depth rather than acrotelm controlled decomposition in dry peats so that decomposition was limited not by the normal (in peatlands) anaerobic conditions but rather temperature limited by permafrost presence. The duration of time for decomposition (e.g., number of years with unfrozen conditions) in the active layer controlled the degree of peat decomposition, which has both a direct effect (loss of mass/material) and an indirect effect (increasing bulk density with increased humification). This effect is illustrated by contrasting patterns of peat accumulation and permafrost formation, with a higher degree of decomposition prior to permafrost preservation (Figure 1) and lower carbon loss with subsequent thaw in wetter conditions (Figure 6).

The higher degree of decomposition in the wet landscape scenario is related to the shallow water table, which results in a larger thermal mass, enhanced heat transfer to deeper peat and subsequently a deeper active layer (Figures 1, 5). In an analysis of active layer depth and soil water content data from sites across Alaska, Clayton et al. (2021) found that while active layer depth generally increased with increasing active layer soil water content, as water content in the top 12 cm of the soil increased the active layer also increased, and active layers also tended to be deeper as soils became saturated. These exceptions reflect the conditions of the TKP wet scenario. With this deeper active layer, the residence time for a peat cohort within the active layer (transit time from the peat surface to the permafrost surface) was about 1,000 ± 100 years or more (Supplementary Figure S2). By then the peat was highly decomposed, had lost most of its original litter mass (mean: ∼73% mass loss), and the organic matter density was near the maximum, resulting in less peat thickness per unit mass of peat and a shallower peat profile (Figure 2B), and relatively low rates of peat accumulation (overall: 2.9 mm decade⁻¹; Figure 1C).

In contrast, when the peatland began to shed water earlier in its development, the drier peat, deeper water table, and resultingly shallow active layer led ultimately to deeper peat with lower organic matter density (Figures 1, 2). These results show that higher peat accumulation and preservation can occur with a deep water table and dry peat surface, which are usually associated with low rates of peat accumulation such as those found in permafrost palsas and bogs (Treat et al., 2016a). Notably, the rapid peat accumulation rates simulated during c.5700–4550 BP were not due to years of very high accumulation (i.e., >2 mm y⁻¹ or 20 mm decade⁻¹), but rather no years of decomposition that exceeded litter input. In contrast, negative accumulation years were common during other periods during the dry scenario and the entire wet scenario (Figures 1C,D). These results illustrate an important role for the decomposition of older peat (roughly 200–800 years in this case) found in the deeper active layer and overall C balance of the peatlands more generally.

Contrasts between wet and dry landscape scenarios (Figure 1) present a possible explanation for how rapid rates of peat accumulation could occur as a result of shallow and relatively rapid permafrost formation in dry peat. This rapid permafrost aggradation coupled with a shallow active layer provides an alternative, or additional, explanation to the hypothesis that high rates of peat accumulation are driven primarily by enhanced seasonality attributed to solar insolation in the early Holocene (Lavoie and Richard, 2000; Jones and Yu, 2010). In a core, the rapid peat accumulation in permafrost would appear as well-preserved, relatively undecomposed peat at depth; unfortunately, the profile description from TKP is not detailed enough to determine whether this is observed in-situ but this seems unlikely given the relatively good agreement between the TKP core and the wet scenario (Figure 2; Table 1). Future work could focus on testing this hypothesis.

Other Arctic peat deposits support the hypothesis for rapid peat accumulation under dry permafrost conditions. Peat cores collected in a polar desert (1980–2010 MAAT: −17.1°C; Canadian Climate Normals, Environment Canada, Alert Station) on Northern Ellesmere Island, Canada at southern Piper Pass show extensive (>3 m thick) peat deposits that formed on the edges of steep banks during the mid-Holocene (LaFarge-England et al., 1991). High peat accumulation rates of 12–38 mm decade⁻¹ are associated with well-preserved plant macrofossils, indicating the rapid incorporation of peat into permafrost before it decomposes (syngenetic permafrost). The species assemblage indicates periods of dominant xeric moss species followed by soligeneous fen mosses and vascular plant species. These dry, well-drained conditions supported peat formation primarily through snowmelt water inputs, which enabled vegetation and moss growth, and the rapid incorporation of this material into permafrost prior to decomposition, similar to the dry model scenario (Figures 1B, 5B).

Permafrost C vulnerability and resilience to warming in these simulations was determined by both the water table position (related to the landscape setting and regional hydrology) and the legacy of the peat profiles themselves. In the future scenario in the dry landscape setting, warming soil temperatures (Figure 3) released the temperature control on decomposition rates, which were multiplied by aerobic conditions in the peat acrotelm (Figure 5). The legacy effect of previous decomposition also affected the vulnerability: the lower degree of decomposition in peat currently residing at the top of the permafrost in the dry landscape setting (Figure 5), when combined with a strong increase in peat temperature (Figure 3), increases the peat vulnerability to decomposition with permafrost thaw (Figure 6). Similar results have been shown in incubations of permafrost peat, where poorly decomposed peat in the permafrost is decomposed at similar rates to peats from the active layer (Treat et al., 2014). The dry conditions also decrease the resilience of these peatlands through lower rates of productivity (Figure 4), which are unable to offset losses of peat C (Figure 6), ultimately leading to substantially higher predicted net C loss, e.g., higher C vulnerability with permafrost thaw than in the wet landscape scenario (Figure 6). Therefore, while soil temperatures also warmed, anaerobic conditions associated with hydrology (Figure 5, high WT) and higher NPP (Figure 4; driven by less drought stress) led to lower net C losses (Figure 6).

Under wetter conditions, more humified peat accumulated more slowly due to a longer residence time in the active layer (Figure 1, Supplementary Figure S2), but then was more resilient against C loss in future scenarios (Figure 6). This was due to increases in vegetation productivity enabled by sufficient moisture (Figure 4) and smaller decomposition peat losses (Figure 5).

We note that the model simulations presented here have some limitations to assessing permafrost vulnerability: they currently do not include abrupt permafrost thaw, which will likely affect ice-rich permafrost peatlands (Olefeldt et al., 2016) and strongly alter ecosystem function and C exchange (Turetsky et al., 2020). If abrupt thaw results in a water table closer to the surface, this could potentially enhance vegetation productivity similar to the outcome in the wet landscape scenario (e.g., Camill et al., 2001: Figure 4) which might offset more of the peat C losses (Figure 6). Changes to the regional hydrology in the future due to permafrost thaw have a key role in the future vulnerability but are very difficult to predict (Quinton et al., 2010).

Some field studies in Alaska have shown substantial differences in the observed peat C dynamics after permafrost thaw that could be explained by timing of permafrost aggradation and peat accumulation demonstrated by the modeling scenarios. In Southwestern Alaska, permafrost peatland sites at Innoko and Koyukuk showed evidence of syngenetic permafrost aggradation within the peat (Jones et al., 2017; O'Donnell et al., 2012). At the APEX Bog site in interior Alaska, permafrost aggraded millennia after peat accumulation (e.g., epigenetic permafrost; Manies et al., 2021). Following permafrost thaw, substantially less carbon was lost from the epigenetic permafrost site than the syngenetic sites that accumulated peat and permafrost simultaneously (Jones et al., 2017; Manies et al., 2021). In the modeled scenarios, rapid permafrost accumulation in the dry scenario (e.g., syngenetic formation) showed carbon losses with warming that were nearly twice that of the wet scenario (Figure 6), where the longer residence time in a non-frozen state allowed for greater decomposition and a less labile carbon source upon thaw (Figure 5). While the context differs from the modeled scenario presented here both in terms of permafrost history (significantly longer permafrost-free conditions at the epigenetic site in the chronosequence) and permafrost thaw (abrupt thaw instead of active layer deepening), these model results provide a plausible explanation for these differing vulnerabilities of these permafrost peatlands.

The main result shown here is that wet and dry landscape settings can result in different peat accumulation histories and peat qualities, even in the same climate setting (Figure 1). If the active layer depth is shallow (dry landscape scenario), peat is less degraded before freezing into the permafrost and so is more susceptible to enhanced decomposition upon permafrost thaw. Because the wet scenario experienced greater decomposition prior to making it into the permafrost, less of it decomposes following thaw. Both peat accumulation history and peat quality have implications for peat C loss in future warming scenarios (Figures 5, 6) due to both permafrost thaw and also legacy effects. To simulate this in a model requires permafrost dynamics (common to many models) and peat stratigraphy, i.e., resolving the peat quality (degree of humification or decomposition) down the peat profile into the permafrost. This probably does not necessarily require annual litter cohorts as used in HPM-Arctic, but a single or just a few peat pools would probably not characterize the peat quality accurately enough.

Very few models simulate peat accumulation over millennia in the permafrost zone. This is a complex task for most earth system models as only a few incorporate peatlands into their land surface schemes (Kleinen et al., 2012; Largeron et al., 2018; Qiu et al., 2018; Chaudhary et al., 2020), and millennial length earth system model simulations are computationally costly. Often, the peatlands are initialized using datasets for timing and peatland extent. Here, we show that key parameters for determining carbon vulnerability and feedbacks, such as peat thickness, bulk density, and humification, must be considered and may vary depending on permafrost history. Some of these important uncertainty factors may be quantifiable from the global peat core databases (Loisel et al., 2014; Treat et al., 2016a; Treat et al., 2016b), but this may require details that many researchers may not have included their analyses. A first step would be to determine whether the hypothesis presented here (dry landscape setting, presence of permafrost, and shallow active-layer depth coincident with rapid peat accumulation) can be evaluated with existing peat core data, and if so, to determine whether these conditions can be associated with landscape topography and climate conditions in a way that provides some predictive power for peat initialization.

Furthermore, it is important to keep in mind that future peat vulnerability to climate change is related to several factors: thaw rate, subsequent water table dynamics, and nature of thawing peat (i.e., susceptibility to further decomposition, which may depend on both the inherent tissue lability of the plants comprising the peat and the degree of humification upon freezing into the permafrost. While this study only addresses this final point, we show that decomposition processes (i.e., degree of decomposition and litter quality) are quite important for both peat accumulation rates and post-thaw carbon losses.

Ultimately, the differences in past C accumulation rates and future C losses between these wet and dry landscape setting scenarios were driven by the model parameterization representing the regional hydrology, including the water inputs from the surrounding watershed and the runoff to downstream areas. By demonstrating the differences in peat profiles between the wet and dry scenarios using the same climate drivers, as well as an inter-site comparison between two sites with contrasting hydrology (TKP and southern Piper Pass), these results suggest that permafrost behavior and peat profiles are strongly dependent on landscape characteristics controlling hydrology (Woo and Young, 2006). The regional hydrology and wetland basin characteristics subsequently determined both peat accumulation history and future changes in peat organic matter stocks with permafrost thaw and on-going warming. These findings point to the importance of considering regional hydrology and wetland basin characteristics due to their controls on wetland hydrology, which will likely be a major control on the response of permafrost carbon in peatlands to warming and permafrost thaw. This study highlights the limitations for accurately predicting changes in C balance in response to permafrost thaw without careful hydrologic assessments. This presents a significant challenge for global-scale models and future emphasis would be better placed on improving hydrologic assessments at local to regional scales.