Thermal Modelling Of Post-Fire Permafrost Change Under A Warming Coastal Subarctic Climate, Eastern Canada

Department of Geography, Environment and Geomatics, University of Ottawa; Northern Environmental Geoscience Laboratory, Department of Geography and Planning, Queen’s University (corresponding author). Email: yifeng.wang@queensu.ca Department of Geography, Environment and Geomatics, University of Ottawa. Email: alewkowi@uottawa.ca Department of Geography, Environment and Geomatics, University of Ottawa. Email: jhollowa@uottawa.ca Northern Environmental Geoscience Laboratory, Department of Geography and Planning, Queen’s University. Email: robert.way@queensu.ca


INTRODUCTION
Forest fires are a widespread natural disturbance in many permafrost regions. The post-fire loss of organic layer and canopy cover and related changes in snow distribution have long-lasting effects on the underlying frozen ground (Jorgenson et al. 2010). Immediate and short-term increases in ground temperatures have been reported post-fire, resulting in active layer thickening and suprapermafrost talik formation (Gibson et al. 2018). However, there is a paucity of research on longterm post-fire permafrost response, particularly within the context of a warming climate.
Thermal modelling is a useful approach for examining the response of frozen ground to surface disturbance and is essential for predicting change. One-dimensional thermal modelling has been used to study post-fire permafrost response within the North American boreal forest, but almost exclusively for subarctic sites in Alaska and the Northwest Territories (e.g., Nossov et al. 2013, Brown et al. 2015, Zhang et al. 2015. Here, we employ thermal modelling to examine the impacts of forest fire on frozen ground near Nain, Nunatsiavut ( Figure 1). This region is wetter than northwest North America, resulting in longer fire recurrence intervals (Foster 1983, Coops et al. 2018, faster surface organic layer growth (Williams and Flanagan 1998), and greater snow accumulation (Way and Lewkowicz 2016), all of which affect the response of permafrost to fire.

STUDY AREA
Coastal Labrador comprises the easternmost part of the Canadian Shield and is mostly characterized by igneous and metamorphic bedrock with thin layers of medium-to fine-grained materials deposited during and following the retreat of the Laurentide Ice Sheet (12-6 k years BP; Dyke 2004, Roberts et al. 2006, Bell et al. 2011. The subarctic climate is influenced by the Labrador Current which carries cold arctic waters southward, resulting in long, cold winters and short, cool summers (Banfield and Jacobs 1998). The mean annual air temperature at Nain is -2.5°C, average total rainfall is 450 mm, and average total snowfall is 475 cm (1981-2010 climate normal; Environment Canada 2020). Nain is located in the sporadic discontinuous permafrost zone (Heginbottom et al. 1995, Way andLewkowicz 2016), but permafrost has only recently been observed at nearby forested sites where it is thought to be climate-driven, ecosystem-protected, following Shur and Jorgenson (2007).
Forests around Nain are composed of black spruce (Picea mariana), white spruce (Picea glauca), balsam fir (Abies balsamea), and eastern larch (Larix laricina; Roberts et al. 2006). The estimated fire return interval of 1501-5000 years for this region (Coops et al. 2018) greatly exceeds the estimated 75-200 year intervals for the western North American boreal forest (Foster 1983).
Our simulations relate to ground conditions at a burned site located on a north-facing slope along Webb Bay (WB; 56.706°N, 62.209°W), ~25 km northwest of Nain (Figure 1-C). The WB burn was a small (1 km 2 ) low severity fire that occurred in summer 2004 (Brehaut and Brown 2020). The site is located on a series of raised beaches, although erratics are also present. Unburned stands are of mixed composition, with an average tree age of 227 years (Brehaut and Brown 2020).

MODELLING PARAMETERS AND METHODOLOGY
Modelling was undertaken with TEMP/W (Geoslope), a finite element geothermal modelling software. TEMP/W is commonly applied in infrastructure-related studies (e.g., Smith and Riseborough 2010), but this is the first to use it to investigate the effects of fire on frozen ground.

Profile and Material
Properties. An idealised material profile was developed for WB (Table 2). Materials were assumed to be saturated (McClymont et al. 2013). The material stratigraphy and boundaries were estimated from field measurements of organic layer thickness, laboratory-based analyses of near-surface sediment samples, and interpretations of sub-surface geophysical investigations (Wang 2020), supplemented by descriptions of regional geomorphology (Natural Resources Canada 1957, Bell et al. 2011. The base of the model was set at a depth of 100 m (Brown et al. 2015), and mesh elements measured 0.01 m (surface-0.   (Wang 2020). b bedrock was assumed to be gabbro (Natural Resources Canada 2017). c thermal conductivities were calculated using a geometric mean (Farouki 1981). d heat capacities were calculated using a weighted average (Johnston 1981 Hope et al. 2016). Projections were compiled for the current century under moderate (RCP4.5) and high (RCP8.5) warming scenarios (Figure 2-A). Projected daily temperatures were statistically downscaled to WB using the delta method (Ramirez-Villegas and Jarvis 2010), with anomalies calculated from a 2010-2017 overlap period. Ground surface boundary temperatures were generated by applying n-factor transfer functions to the air temperature record. The freezing n-factor (nf) and thawing n-factor (nt) summarize the influences of snow cover and vegetation shading on ground surface temperatures (Karunaratne and Burn 2004). An nf of 0.36 and an nt of 0.73, measured in the unburned region from 2017 to 2019, were applied throughout the unburned scenarios. An nf of 0.12 and an nt of 0.84, measured in the burned region for the same period, were used to inform and calibrate postfire changes in the burned scenarios. Nf and nt were adjusted in five year increments according to a fifty-year linear two-segment model, in which values recovered to 66% of the pre-fire conditions after the first 25 years following fire, before fully returning to pre-fire conditions after another 25 years (Zhang et al. 2015; Figure 2-B). Low geothermal heat fluxes (22 mW/m 2 ) were measured 50 km south of the study sites (Mareschal et al. 2000). To avoid producing unrealistically thick permafrost, the geothermal heat flux was set to 220 mW/m 2 , which is still below the values used to model permafrost mounds in nearby southern Labrador (Way et al. 2018).

Figure 2. A) Mean annual air temperatures estimated for WB (see text for data sources); B)
N-factors for burned scenarios at WB. The minimum freezing n-factor was set to 0.05. Thermal Modelling with TEMP/W. Six basic simulations were run to evaluate the response of frozen ground under moderate (RCP4.5) and high warming (RCP8.5) with no disturbance and low and high severity fire disturbance. Six additional simulations were run to test the model sensitivity to organic layer accumulation by excluding post-fire organic material regeneration. Each time step (30 minutes) was iterated up to 900 times to allow convergence within 0.005 K. The thermal model was initialized using the 1950-1979 daily average climate normal until equilibrium was achieved (node variability <0.01 K per 100 years; Zhang et al. 2015). The model was then run forward in time from 1965 using the daily air temperature composite for WB, comprised of BESTP (1965-2017) and RCP4.5 or RCP8.5 projections (2018-2099; Figure 2-A). Fire disturbance was modelled as a decrease in organic layer thickness and nf and an increase in nt, applied as of August 1, 2004 (Figure 2-B). Low and high burn severities were represented by a modelled reduction in the organic layer by 5 or 15 cm, respectively. Following this loss, the organic layer continued to accumulate in the model at a rate of 3.6 mm/year, based on in situ measurements (Wang 2020). Gradual adjustments to nf and nt were applied during the post-fire recovery period (Figure 2-B).

Model Output Processing and Validation.
Model outputs were processed in R (version 3.6.2; R Core Team 2019, Wang 2020) with results portrayed as frozen (dark blue; <-0.025°C; <11.39% UMC), partially frozen (light blue; between -0.025°C and 0°C; 11.40-99.99% UMC), or unfrozen ground (red; >0°C; 100% UMC). Model validation used ground surface temperatures, as boreholes were not established due to permitting limitations. Ground surface temperatures were measured at four-hour intervals using Thermochron DS1922L iButtons (Maxim; ±0.5°C accuracy), buried 5 cm below the ground surface in both burned and unburned locations at WB.

Model Validation.
Model outputs were compared to ground surface temperatures measured from July 2017 to July 2018. The modelled annual surface temperature averaged 0.41°C higher in the unburned region (daily RMSE: 2.53°C; Figure 3) and 0.13°C higher in the burned region (daily RMSE: 2.20°C; not shown). Discrepancies occurred following the transition from the historical to projected air temperature inputs in January 2018 and in the spring when observations showed a slower melt out than the modelled values due to late-lying snow.   (Table 2), mainly from the base up. A talik does not develop and there are no notable changes in active layer thickness. In this scenario, the temperature at DZAA increases at a linear rate of ~0.004°C/year from 2018 to 2099. Under RCP8.5, a climate-initiated talik develops in the 2060s (Figure 4-B). After talik initiation, the active layer, now defined by the depth of seasonal freezing, thins. The talik deepens throughout the remainder of the century, but a very thin layer of permafrost persists in 2099 (Table 2). Fire disturbance triggers talik development in all burned scenarios (Figure 4). Under RCP4.5, the thaw front advances rapidly into the permafrost at first, but gradually slows until the talik reaches a depth of 5-6 m. The active layer, defined by the depth of seasonal freezing, initially thins by an average of 35-46 cm when compared to the unburned scenarios over the same period. This thinning trend reverses once the talik reaches its maximum thickness and begins to refreeze. The temperature at the DZAA decreases and permafrost re-aggrades, causing the talik to close after three decades in the low burn severity simulation and four decades in the high burn severity simulation (Table 2). However, much of the re-aggraded permafrost is only partially refrozen, with temperatures between -0.005 and 0°C. Following talik closure, active layer thicknesses in the burned scenarios average 8-28 cm greater than those in unburned scenarios. By the mid-2090s, the temperature at the DZAA is the same in both burned and unburned scenarios, and 2-4 m of permafrost remains at the end of the century (Figure 4-C and 4-E; Table 2). Under RCP8.5, low severity fire disturbance results in the development of two taliks and the complete degradation of permafrost by 2080. The first talik is initiated immediately following fire, and the thaw front advances for approximately two decades. This is followed by a 20-year period of stability until refreezing occurs from the active layer downwards (Figure 4-D). The ground layers comprising the first talik only partially refreeze and then thaw out during the development of the second talik. In the high burn severity simulation, the fire-initiated talik does not refreeze, and permafrost degrades completely by 2069 (Figure 4-F; Table 2). Impact and Importance of Organic Layer Accumulation. Sensitivity analyses with no organic layer accumulation result in thaw from the surface downwards in all unburned scenarios. Under RCP4.5, a climate-initiated talik forms in 2082 and deepens so that only 4.9 m of permafrost remains in 2099 (Table 2). Under RCP8.5, a climate-initiated talik begins to form in 2059, 4 years before the talik in the comparable simulation with organic layer accumulation (Figure 4-B), and permafrost is lost by 2095 (Table 2). In all burned scenarios, omitting surface organic layer regeneration results in the total loss of permafrost by the end of the century. Under RCP4.5, permafrost is lost by 2082 and 2066 for the low and high burn severity simulations, respectively (Table 2). Under RCP8.5, the fire-initiated talik does not refreeze in either low or high burn severity simulation, and permafrost thaws by 2069 or 2060, respectively (Table 2).

DISCUSSION
Post-Fire Permafrost Response. In our simulations, permafrost thickness is relatively stable and varies by only 0.5 m (3%) during the first 40 model years. A local maximum in permafrost thickness and a corresponding minimum in the temperature at DZAA occurs in 1995 (Figure 4), reflecting the end of a period of regional climatic cooling (Way and Viau 2015; Figure 2-A).
Our simulations present a shallower active layer for the first two decades immediately following fire. This is because a supra-permafrost talik forms in the first year and the base of the active layer is then defined by the maximum depth of seasonal freezing rather than the maximum depth of seasonal thaw. Most field studies describe thickening of the active layer following fire (e.g., Zhang et al. 2015, Gibson et al. 2018, Holloway et al. 2020), but late-summer probing cannot be used to differentiate between the active layer and a talik. Our results do not contradict previous studies, but they emphasize the need for ground temperature monitoring in the field to distinguish between these two layers. The overall understanding of post-fire permafrost response would also benefit from a clear differentiation between an active layer and a talik in modelling studies.
Fire disturbance can trigger thawing of permafrost from the base, as has been modelled for forested regions in the Northwest Territories (Zhang et al. 2015). Despite the relatively high geothermal heat flux that was applied to the base of our model domain, the post-fire persistence of permafrost at the WB site speaks to the resiliency of perennially frozen ground in the region. Omission of surface organic layer accumulation caused rapid permafrost degradation due to the impacts of rising air temperatures. The organic layer is largely responsible for the thermal offset, which enables permafrost persistence following disturbance or in unfavourable climates (Smith and Riseborough 2002). Our results support previous findings of the importance of the organic layer in post-fire permafrost stability (Fisher et al. 2016) and of including organic layer regeneration in post-fire permafrost modelling (Zhang et al. 2015).
Permafrost Resiliency. The resiliency of the permafrost at WB can be assessed in relation to the climate-ecosystem-permafrost types defined by Shur and Jorgenson (2007). Prior to fire disturbance, permafrost is interpreted to be climate-driven, ecosystem-protected. Following fire, permafrost re-aggrades via ecosystem-driven processes and is re-classified as ecosystem-protected permafrost. Under RCP8.5, climate warming overcomes the ecosystem's ability to protect the permafrost, resulting in the development of a climate-initiated talik and the loss of permafrost.
Our simulations agree with previous studies that have demonstrated that fire acts to accelerate the loss of permafrost (e.g., Zhang et al. 2015, Gibson et al. 2018. Under the same climate scenario, burned simulations result in thinner and warmer permafrost that thaws out more quickly than in unburned simulations. Under RCP4.5, by 2100, low and high severity fire disturbances at WB result in respective permafrost thicknesses that are 56% and 73% thinner than the ~8 m of permafrost remaining in the undisturbed simulation. Loss of permafrost due to fire disturbance in the Northwest Territories is modelled to occur an average of 8 years earlier when compared to undisturbed sites (Zhang et al. 2015). At WB, loss of permafrost following fire is modelled to occur at least 24 years earlier when compared to undisturbed simulations, at a rate that is three to four times faster than that reported for the West. This accelerated degradation is likely due to the warmer, less ice-rich permafrost (i.e., with lower latent heat) at our site compared to conditions in the central Northwest Territories (Smith et al. 2015) and may also reflect the slower post-fire forest regeneration processes reported for WB (Brehaut and Brown 2020) compared to the pulse recruitment that has been documented following fire in the West (Johnstone et al. 2004).
The development of climate-initiated taliks in the latter half of this century under RCP8.5 shows the sensitivity of permafrost at WB to the impacts of climate warming. Taliks produced by warming deepen more linearly than those initiated by fire disturbance. Fire-initiated taliks develop in the immediate post-fire period due to the sudden effects of increased heat penetration in summer and decreased heat loss in winter (Gibson et al. 2018). However, post-fire talik development eventually slows as the ecosystem recovers from disturbance (Gibson et al. 2018). By contrast, permafrost degradation due to climate warming occurs more gradually, and climate-initiated taliks do not re-freeze (O'Neill et al. 2020).
Modelling Limitations. The modelling involves several generalizations that could have affected the individual predictions but are not thought to have altered the overall results. First, n-factors were held constant for 5-year periods (Figure 2-B), but in the field, they can vary over short distances and time periods (Karunaratne and Burn 2004). Nf is sensitive to inter-annual differences in snow cover (Karunaratne and Burn 2004) and can be impacted by permafrost presence or absence (Riseborough and Smith 1998). A reduction in nf following talik development would delay talik recovery and possibly accelerate permafrost loss. No adjustment method is currently available for this purpose. Second, post-fire variations in soil moisture, which typically relate to the remaining organic layer thickness (Kasischke and Johnstone 2005), were not considered due to limited field-based validation. Post-fire soil moisture is generally expected to increase due to thickening of the supra-permafrost layer, poor drainage, and reduced evapotranspiration (Nossov et al. 2013). Soil moisture is eventually expected to return to near pre-fire conditions as the site recovers from the fire disturbance (Smith et al. 2015, Gibson et al. 2018. Future simulations incorporating dynamic soil moisture and the lateral and/or vertical flow of water would provide insight into how changes to the post-fire hydrological regime might influence soil freezing properties and characteristics. Finally, equal-weighted averaging was performed to assemble the composite record for projected air temperatures (2018-2099) under RCP4.5 and RCP8.5 ( Figure  2-A). The resulting composite records mask inter-model variability. Individual climate scenarios would have produced more variability in the timing and rate of post-fire permafrost change.

CONCLUSION
The Webb Bay study site near Nain, Nunatsiavut is wetter, accumulates a deeper winter snowpack, and experiences slower and more continuous post-fire forest regeneration processes compared to sites in western North America where most research on forest fire and permafrost has been carried out. Despite these differences, the results of our thermal modelling agree with the prevailing view that fire disturbance accelerates permafrost loss due to climate warming. Following fire and within the context of a warming climate, permafrost at this coastal boreal site is modelled to persist to at least 2060 and, in many cases, through to the end of this century. Simulations highlight the role of the organic layer in permafrost persistence under otherwise unfavourable climatic conditions and the importance of differentiating between the active layer and a talik, if present, in the suprapermafrost layer. The results extend geographic coverage and improve our understanding of active layer change and talik formation as part of the overall long-term post-fire permafrost response.