Wildfire‐Initiated Talik Development Exceeds Current Thaw Projections: Observations and Models From Alaska's Continuous Permafrost Zone

As the Arctic warms and wildfire occurrence increases, talik formation in permafrost regions is projected to expand and affect the cycling of water and carbon. Yet, few unified field and modeling studies have examined this process in detail, particularly in areas of continuous permafrost. We address this gap by presenting multimethod, multiseasonal geophysical measurements of permafrost and liquid‐water content that reveal substantial talik development in response to recent wildfire in continuous permafrost of boreal Alaska. Results from observation‐based cryohydrogeologic model simulations suggest that predisturbance subsurface conditions are key factors influencing thaw response to fire disturbance and air temperature warming. Our high‐resolution integrated study illustrates enhanced vulnerability of boreal continuous permafrost, with observed talik formation that exceeds coarse‐scale model projections by ~100 years even under the most extreme future emissions scenario. Results raise important scaling questions for representing extreme permafrost thaw phenomena of growing widespread importance in large‐scale predictive models.


Introduction
Permafrost degradation in northern high latitudes is ongoing and expected to continue with impacts on water resources, infrastructure, and ecosystems (Meredith et al., 2019). Thermal models at various scales and complexity (Delisle, 2007;Parazoo et al., 2018;Zhang et al., 2008) project a thickening active layer (ground that freezes and thaws seasonally) in response to warming trends in permafrost regions (e.g., Koven et al., 2013;McGuire et al., 2016;Slater & Lawrence, 2013). Recent studies underscore the growing importance of talik (perennially fully unfrozen ground) formation as the next phase of permafrost thaw following active-layer thickening (Connon et al., 2018;Devoie et al., 2019;Lamontagne-Halle et al., 2018;Streletskiy et al., 2015a;Walvoord et al., 2019). In this advanced thaw phase, taliks form as the depth to permafrost exceeds the depth of annual freezing. Talik formation is a precursor to accelerated permafrost thaw (Lawrence et al., 2008) and influences the magnitude, seasonality, and pathways of soil carbon released in the gas and solute phases (Commane et al., 2017;Vonk et al., 2019;Walvoord et al., 2019). A pan-Arctic-scale modeling study by Parazoo et al. (2018) projects widespread talik development across permafrost regions in the coming decades; consequently, we need to address major gaps in estimating talik formation, spatial distribution, and function. Measurements of liquid-water content from terrestrial taliks unassociated with surface-water features are typically nonexistent, despite being a key determinant for predicting thaw acceleration (Romanovsky & Osterkamp, 2000), advective energy transport (Kurylyk et al., 2016;Lamontagne-Halle et al., 2018;Zipper et al., 2018), and winter soil respiration (e.g., Natali et al., 2019).
The vulnerability of continuous permafrost to talik formation is relatively unknown as direct observations of terrestrial taliks are limited and have generally been confined to discontinuous to sporadic permafrost (Connon et al., 2018;Gibson et al., 2018). Field detection typically relies on frost probe and ice auger measurements (Connon et al., 2018;Gibson et al., 2018) or vertical temperature profile monitoring (Streletskiy et al., 2015b), all of which can provide challenges in talik characterization. For example, it is difficult to penetrate shallow frozen soils with manual techniques during winter months to detect unfrozen zones below. Additionally, using temperature data and a 0°C threshold for soil freezing can underestimate talik formation, as >25% liquid-water content can exist in soils at temperatures below 0°C (e.g., Watanabe & Mizoguchi, 2002;Yoshikawa & Overduin, 2005).
Fire activity in Northern Hemisphere permafrost regions is projected to increase in the coming decades (Meredith et al., 2019) and has been shown to enhance carbon respiration via deep soil carbon decomposition (Gibson et al., 2019). Wildfire promotes permafrost thaw through the combustion of insulative organic soil and reduced shading from vegetation, both of which enhance subsurface summer warming (e.g., Yoshikawa et al., 2002). However, field studies demonstrate a wide range in permafrost response to fire (e.g., Minsley et al., 2016) rendering thaw prediction across large heterogenous landscapes challenging. For example, a study by O'Donnell et al. (2011) from Alaska's discontinuous permafrost zone found no evidence of permafrost thaw >1 m 4 to 40 years after fire. In contrast, Nossov et al. (2013) found thaw depths >1 m across most of their study sites in Alaska's discontinuous permafrost zone, with some burns that resulted in talik formation. However, environments with high active-layer moisture or ice content, similar to continuous permafrost settings, are considered relatively resilient to wildfire (Jorgenson et al., 2010;Nossov et al., 2013).
This study integrates field observations and cryohydrogeologic modeling to assess the vulnerability of continuous permafrost in boreal Alaska to talik formation after wildfire. We present novel multimethod, multiseasonal geophysical data, including measurements of seasonal in situ liquid-water content to quantify subsurface response to fire disturbance. Data-informed model simulations enable examination of the role of prefire conditions in influencing the rate and magnitude of talik development. This unique data set shows that continuous permafrost is vulnerable to talik formation sooner than previously projected by coarse-scale modeling that does not account for disturbance, even under extreme warming scenarios.

Study Site
In interior Alaska, the continuous permafrost zone near the continuous/discontinuous transition has mean annual air temperatures of < −5°C (Jorgenson et al., 2008;Figure 1a). South of the Brooks Range, vegetation in unburned areas with underlying permafrost typically comprises open canopy black spruce forest and shrubs with lichens and an organic layer consisting of mosses (Calef et al., 2015;Turetsky et al., 2011). Burned area vegetation is often composed of deciduous trees, grasses, forbs, and shrubs that vary based on the time since fire and fire severity (e.g., Bernhardt et al., 2011;Brown et al., 2016); additionally, these areas commonly have an organic layer thickness that has been reduced from fire Genet et al., 2013).
Field data were collected near the West Fork of Dall Creek (WFD),~100 km southwest of Coldfoot, AK, along the Dalton Highway, south of the Brooks Range (Figure 1a) near the continuous/discontinuous permafrost transition. A 15,000 m 2 burn scar (Dall City Fire 2004) at WFD is bordered by a small ephemeral tributary to the north, the creek to the east, and unburned black spruce forest to the south and north (Figure 1b). Burn scar vegetation consists of a thin (~0-5 cm) layer of sphagnum moss with interspersed grasses and forbs. The unburned area surrounding the burn is extensive and representative of the larger watershed. In burned and unburned areas, sphagnum moss overlies a thin cryoturbated organic-mineral soil mixture that transitions to mineral soil characterized as clay loam . Previous soil observations show liquid-water saturation values at the field site (burned and unburned) to be~85% when thawed  and thus ice rich when frozen. Electrical resistivity tomography (ERT) data collected in September of 2016 showed shallow, continuous permafrost within unburned areas and deep thaw within the burn scar ( Figure 1c; Minsley et al., 2018). These data suggested a potential talik within the burn but only provided a single observation during maximum annual thaw.

Geophysical Measurements
Geophysical methods were used to capture subsurface talik characteristics, including liquid-water content, that are unattainable with typical physical measurements. Nuclear magnetic resonance (NMR) measures in situ liquid-water content (Kenyon, 1997) and has been shown to provide reliable measurements of liquid-water content in permafrost regions (Kass et al., 2017;Minsley et al., 2016). Ground-penetrating radar (GPR) is an established method to determine permafrost depth due to its sensitivity to liquid-ice transitions (Bradford et al., 2005;Hinkel et al., 2001). GPR reflections are caused by contrasts in relative dielectric permittivity (Annan, 2009), which is distinctly different for water (~81) and ice (~3-5). Consequently, GPR reflections at WFD are primarily a function of contrasts in water phase.
NMR data were acquired with a Dart down-borehole tool (Vista Clara, Inc., Mukilteo, Washington). The Dart measures a~1-2 mm cylindrical shell of soil with a diameter of~0.15 m centered on the tool and averages over a vertical interval of~0.23 m. Boreholes were augered to various depths up to~2.45 m, and NMR data were logged upward from the bottom of the borehole in 0.25 or 0.125 m intervals. Measured T2 decay time series were analyzed using proprietary software (Vista Clara, Inc.) that estimates volumetric liquid-water content by fitting a multiexponential function to the measured T2 decay-curve data (Walsh et al., 2013).
GPR data were collected in winter 2019 using a Pulse EKKO-PRO (Sensors and Software, Inc. Ontario, Canada) with 100 MHz central-frequency antennas. Common-offset surveys had 0.25 m spacing between traces and an antenna separation of 1 m. A common midpoint survey (CMP) was conducted using an initial offset of 1 m with 0.1 m step sizes moving outward. Trace locations were collected using an Arrow series EoS GNSS receiver (Eos Positioning Systems, Inc., Terraborne, QC, Canada). For detailed information on NMR and GPR collection and processing see supporting information.

Other Physical Observations
Vegetation, organic layer thickness, and permafrost depth were measured along the previously collected ERT transect (Figure 1c; Minsley et al., 2018). Two soil pits (burned and unburned) were excavated, and intact, thawed, mineral soil samples were collected from various depths. Soil samples were analyzed in the laboratory in unfrozen conditions to attain saturated soil-water content, field-saturated hydraulic conductivity, soil-water retention curves, volumetric specific heat capacity, and thermal conductivity . During soil pit excavation, a handheld time domain reflectometry (TDR) soil moisture probe was used to collect observations of liquid-water content later compared to adjacent NMR measurements.

Cryohydrologic Models
The USGS SUTRA code (Voss & Provost, 2002), with extended capabilities to account for freeze/thaw processes (SUTRA-ICE), was used to explore permafrost response to fire disturbance and to quantify the influence of prefire subsurface conditions on the rate and magnitude of talik development. SUTRA-ICE has been used in previous work (e.g., McKenzie & Voss, 2013;Evans & Ge, 2017;Kurylyk et al., 2016;Wellman et al., 2013;Zipper et al., 2018) and has been benchmarked against similar codes (Grenier et al., 2018). We used a 1-D variably saturated model since low field-saturated hydraulic conductivities (10 −6 to 10 −7 m s −1 ) at WFD  suggested negligible advection. The model was 100 m deep to minimize the influence from the bottom thermal boundary condition on the shallow areas of interest, and field observations were used to constrain model parameters (Tables S3). Pressure and temperature distributions from a steady-state model were used as initial conditions for a starting unburned model that had a 0.25 m thick organic layer overlying mineral soil. Top temperature boundary conditions were time varying (sinusoidal), and the temperature function was derived from mean annual air temperatures and annual temperature fluctuations calculated from 1970 through present for the three closest weather stations (Bettles, Coldfoot, and Seven-mile station; https://www.wcc.nrcs.usda.gov/ snow/ and https://wrcc.dri.edu/). Bottom temperature boundary conditions were constant and assumed a geothermal gradient of 2.5°C per 100 m.
An N-factor approximation (Lunardini, 1978) was used to account for surface-energy balance variability that leads to different mean annual air temperatures and mean annual ground surface temperatures in boreal forests (Fisher et al., 2016;Jorgenson et al., 2010). N-factors range from 0-1 where smaller values indicate a higher degree of thermal buffering. Different frozen (nF) and unfrozen (uF) N-factors were used to capture the strong thermal influence of seasonal snow cover at WFD, as previous studies identified significant differences between seasonal N-factors for Alaskan and Canadian soils (e.g., Jorgenson & Kreig, 1988;Karunaratne & Burn, 2004;Smith & Riseborough, 2002). Previously published uF values from black spruce forests in interior Alaska range between 0.5 and 0.7 (Jorgenson & Kreig, 1988;Smith & Riseborough, 1996). Consequently, a central value of uF = 0.6 was used for initial unburned model simulations to represent the open canopy black spruce forest observed at WFD. Frozen N-factors (nF) from previously published literature ranged between~0.12 and 0.44 (Jorgenson & Kreig, 1988;Karunaratne & Burn, 2003;Smith & Riseborough, 2002). Variation in nF from previous literature was primarily a function of mean annual air temperatures and snow depth (e.g., Smith & Riseborough, 2002). Smith and Riseborough (2002) suggest that sites with mean annual air temperatures between~4°C and 6°C (i.e., the range that inludes WFD air temperature observations from the three closest weather stations) and snow depths between~0.4 and 0.7 m (i.e., the range observed at WFD in the unburned black spruce forest during the 2017-2018 and 2018-2019 winters) would correspond to nF values between~0.15 and 0.30. Additionally, Karunaratne and Burn (2003) calculated nF of~0.12-0.44 for five sites across Canada spanning both continuous and discontinuous permafrost, and Jorgenson and Kreig (1988) calculated nF of~0.3 for closed canopy black spruce forest in interior Alaska. For our initial unburned models a nF of 0.3 was selected because it represented a central value from previously published literature and fit nF values from sites with similar vegetation and snow depths to those observed at WFD.
The unburned model was run for 55 years with a warming trend of~0.04°C per year, consistent with previously published data from the 55 years prior to the WFD burn (Stafford et al., 2000). Unburned model results were compared to field measurements from WFD, and outputs were used as initial conditions for burned model simulations. Fire disturbance was simulated by removing 0.20 of the 0.25 m organic layer (i.e., 0-5 cm observed after fire), nF were decreased to account for less snow interception, and uF were increased to account for decreased ground shading. Specficially, uF values were increased from 0.6 to 0.9 for burned model simulations to account for the lack of black spruce canopy to insulate the ground from incoming solar radiation in summer months. Burned model uF fall within the range of observed uF values for Alaskan and Canadian soils with thinner organic layers and sparse tree canopy observed from previous studies (Jorgenson & Kreig, 1988;Klene et al., 2001;Karunaratne & Burn, 2003). Greater snow depths were observed in the WFD burn scar (i.e.,~0.75-1 m) than surrounding forest, likely due to decreased canopy interception of falling snow. To capture this increase in snow depth, nF values were decreased from 0.3 to 0.2 and fall within the range of nF values from previous studies for site conditions. Unburned models with different N-factors generated different subsurface initial conditions for burned model simulations used to quantify the influence of prefire temperature and liquid-water saturation conditions on fire disturbance response. Three models representing cold (nF = 0.4, uF = 0.6), intermediate (nF = 0.35, uF = 0.5), and warm (nF = 0.2, uF = 0.6) conditions were used to capture a wide range of potential prefire surface-energy factors. N-factors for these three simulations were selected to represent surface-energy factors from black spruce forest that would create cold, intermediate, and warm subsurface conditions and fall within the range of previously published N-factor values. Burned models were run for 30 years, 15 years longer than the time since disturbance. Model parameters and N-factors can be found in the supporting information.

In Situ Liquid-Water Content Measurements Reveal Talik Formation in Response to Fire Disturbance
NMR data collected in burned areas during June 2018 and early September 2016 (near-maximal annual active-layer thickness) indicate fully unfrozen profiles of~50% volumetric water content (VWC) from the ground surface to the bottom of the boreholes (Figure 2a). In contrast, NMR profiles collected in burned areas during periods of near-maximal annual soil freezing (May 2017 and April 2019) indicate the presence of three key zones: (1) a frozen zone representing winter active-layer conditions (<25% VWC), (2) a partially thawed zone (~25-50% VWC), and (3) a talik (>50% VWC) with liquid-water content profiles nearly indistinguishable from those collected during periods of maximum thaw (Figure 2a). Currently, distributed latent-heat effects from the high liquid-water fraction observed within the talik and partially thawed zone is likely accelerating permafrost thaw and retarding the depth of seasonal frost at WFD. NMR profiles collected within the unburned black spruce forest show no evidence of a talik. September 2016 and June 2018 NMR profiles showed a distinct active layer with VWC of~60% during maximum thaw transitioning to permafrost with low VWC (~3-10%) at~0.7 m depth (Figure 2b). During maximally frozen conditions, VWC is low (<10% in May 2017 and <25% in April 2019) down to a depth of~1 m, indicating frozen active-layer conditions (Figure 2b). April 2019 data show an area of higher liquid-water content (~15-25% VWC) that suggests warmer permafrost between~0.5 and 0.8 m before decreasing to colder permafrost with lower VWC at~0.8-1 m (Figure 2b). The major contrast in talik occurrence between the unburned and burned black spruce forest highlights the sensitivity of continuous permafrost areas to fire disturbance.

Talik Footprint Extends Beyond NMR Boreholes
Talik depth exceeded our deepest NMR measurements; consequently, 2-D GPR transects were collected in April 2019 (Figure 1b) to locate the bottom of the talik and to determine the taliks' spatial extent. Small variations in permafrost depth and reflection amplitudes exist across the four collected profiles (e.g., Figure 2c); however, each transect displayed the presence of a talik, a strong permafrost table reflection, and minimal reflections from unburned areas. Consequently, the talik footprint likely encompasses the entire burned area. Burned area data show a consistent reflector at depth whose polarity indicates a transition from more polarizable (higher water content) to less polarizable (high ice content) material (Figure 2c). Minimal reflections from unburned areas indicate the presence of homogeneous frozen soil where limited GPR reflections suggest pore water is present in a similar phase (Figure 2c). This conclusion is consistent with previously collected ERT data (Figure 1c), our NMR data (Figures 2a and 2b), and the previously published soil pit data  from WFD.
Analysis from the CMP survey collected within the burn suggests a talik depth of~4 m, or~1.5 m deeper than the NMR boreholes. Multiple reflectors within the talik indicate liquid-water content variations from pockets of variably thawed material; this heterogeneity may explain variations in burned area NMR profiles and is consistent with NMR observations of partially thawed zones. NMR data suggested freezing-front 10.1029/2020GL087565

Geophysical Research Letters
propagation during the winters of 2016/2017 and 2018/2019 only extended 0.5-0.7 m below the surface. Consequently, GPR and NMR data suggest the combined thickness of the partially thawed zone and talik formed over the 15 years postfire to be~3.25-3.5 m in thickness.

Model Results Compare Favorably to Field Observations
Field data were compared to model outputs to ensure general temporal and spatial patterns of modeled liquid-water saturations were capturing observed system behavior. SUTRA-ICE models reproduced active-layer depths of~0.9-1 m similar to those observed via NMR and frost probe measurements from unburned areas (Figure 3a). Additionally, liquid-water saturation profiles obtained from NMR data closely matched those produced by unburned model simulations (Figure 3b). Model results for the burned case, 15 years postfire (i.e., present), shows the development of an envelope of high liquid-water content (talik) with a partially thawed zone below the talik, before transitioning to permafrost at depth (Figure 3c). Liquid-water saturation profiles obtained from the burned model captured the spatial and temporal patterns from NMR data collected in the burn scar (Figure 3d). Similarly, the bottom of the talik generally agreed with average depths to the permafrost table reflection observed in the GPR profiles (Figures 2c and 3c).

Observations Exceed Large-Scale Thaw Predictions
Soil temperatures from a recent coarse-scale modeling study that predicted ubiquitous future pan-Arctic talik formation (Parazoo et al., 2018) were extracted from the grid cell containing our study area for comparison. The model does not project widespread talik initiation for the study region until~2112 even under representative concentration pathway (RCP) 8.5, an extreme warming scenario (Figure 4a).  Figure 4a). However, for depths of 3.2 and 4.0 m minimum seasonal soil temperatures remained above 0°C after 2112, indicating talik formation (Figure 4a). In contrast, field observations from this study show a current talik thickness of~3.25-3.5 m (Figures 2a-2c), with modeled soil temperatures from this study remaining~5°C warmer than previous predictions during periods of maximum annual soil thaw and~12-20°C warmer than previous predictions during periods of maximum annual soil freezing (Figure 4b).
Coarse-scale pan-arctic models cannot be expected to incorporate the plethora of heterogeneity that exists at the field scale, and likewise high-resolution observational studies cannot fully cover the pan-arctic region. The intent of coarse-scale model comparison with our high-resolution observations is to highlight a disturbance response that may become increasingly more important at larger scales as thaw vulnerability expands into the continuous permafrost zone and wildfire occurrence increases. This study demonstrates the sensitivity of continuous permafrost to fire disturbance, a sensitivity that is projected to increase with future changes in the near-surface energy budget such as rising air temperatures, ecological shifts, or changing snow depth and timing (Meredith et al., 2019). Given the ubiquity of fire across the northern high lattitudes (e.g., Calef et al., 2015;Ponomarev et al., 2016) and the projected increase of wildfire occurrence in pan-arctic environments (Meredith et al., 2019), this study highlights the potential need for addressing scaling issues associated with representing extreme permafrost thaw phenomena in large-scale predictive models.

Prefire Subsurface Conditions Affect System Response to Disturbance
Burned model simulations showed that predisturbance temperature and liquid-water saturation influenced the rate and thickness of talik formation. When unburned models with N-factors representing higher winter thermal buffering and lower summer thermal buffering were used as initial conditions, burned model simulations showed more rapid and spatially expansive thaw (Figures 4c-4e). Specifically, warm model simulations experienced~80% more thaw than cold simulations, and intermediate simulations showed~50% more thaw than cold model simulations (Figures 4c-4e). Additionally, warm model simulations took~7 years postfire to develop a 1 m thick talik, intermediate simulations took~12 years postfire to develop a 1 m thick talik, and cold model simulations did not experience a 1 m talik until~21 years postfire (Figures 4c-4e). Consequently, talik development rate may have implications for permafrost recovery postfire, as ground shading from vegetation reestablishment (Briggs et al., 2014) and the reaccumulation of surface organic material (e.g., Jorgenson et al., 2010) will promote permafrost aggradation. Faster thaw rates would promote deeper thaw prior to recovery and may release previously accumulated legacy carbon identified as a driving factor for the long-term arctic carbon balance, thus promoting enhanced climate warming via the permafrost carbon feedback (e.g., McGuire et al., 2012;Schuur et al., 2015). Deep thaw will be particularly prone to legacy carbon release if thaw affects depths with previously unfavorable conditions for organic carbon decomposition (e.g., Ping et al., 2015) or that escaped burning in previous fire cycles (e.g., Chapin et al., 2006;Walker et al., 2019). If areas with prefire subsurface conditions conducive to rapid talik development are affected by wildfire, thaw depths may become unrecoverable when combined with future warming trends despite the reestablishment of vegetation and mosses.

Liquid-Water Content Is a Fundamental Measurement for Talik Detection and Prediction
This study found liquid-water content to be a fundamental measurement for determining talik presence. Using a set temperature threshold can alter talik interpretation in the absence of liquid-water saturation information (Figure 4a). Only warm initial condition scenarios met the cryotic definition (perennially >0°C) of talik presence within the first 15 years postfire (Figure 4b), despite >70% liquid-water saturation being present to depths approaching 3 m for all scenarios (Figures 4c-4e). Our current understanding of liquid-water saturation as a function of soil freezing is limited (Amiri et al., 2018;Kurylyk & Watanabe, 2013), making it difficult to define thresholds for a binary classification of soils as frozen or thawed based on temperature alone. Consequently, field observations or models that rely on temperature thresholds for permafrost table and freezing-front depths may provide a conservative estimate of talik ubiquity, particularly in regions where soils retain significant portions of liquid water below 0°C.

Conclusions
Field observations and site-based high-resolution cryohydrologic models show deep talik formation after wildfire in continuous permafrost where conditions are generally considered to be resilient to major near-term thaw. Observed liquid-water saturations and soil temperatures associated with talik formation likely provide the conditions necessary for cold season respiration of deep carbon that has recently shown to be a significant component of the Alaskan carbon budget. Deep thaw responses to wildfire, such as observed here, could promote the release of organic carbon sequestered in deep permafrost. In areas prone to lateral flow, similar talik features could enhance pathways for dissolved carbon export to rivers and streams. Cryohydrogeologic model analyses performed as part of this study suggest that prefire liquidwater saturation and temperature conditions substantially influence spatial and temporal thaw response to fire disturbance and therefore should be considered in permafrost model predictions. Prefire thermal priming will impact short-and long-term permafrost recovery postfire, particularly as projected changes in surface-energy conditions promote irrecoverable permafrost thaw. Taliks in cold terrestrial hillslopes may become increasingly important in the long-term arctic carbon balance as wildfire activity increases and subsurface conditions become increasingly more sensitive to disturbance.
Talik observations presented here exceed coarse-scale modeling projections of talik formation in the study region under high emission scenario RCP 8.5 by~100 years, with differences in winter soil temperatures often >15°C. Differences between field observations and coarse-scale model projections underscore the growing importance of supportive high-resolution field studies of extreme thaw (Farquharson et al., 2019). Given the biogeochemical significance of taliks in permafrost systems and the potential for aggregated localized thaw to evolve toward regional importance as widespread wildfire activity is expected to increase across circumboreal regions in the future, our findings raise critical issues regarding representation of extreme permafrost thaw phenomena in large-scale predictive earth system models.