Environmental drivers and remote sensing proxies of post-fire thaw depth in eastern Siberian larch forests

Boreal fire regimes are intensifying because of climate change, and the northern parts of boreal forests are underlain by permafrost. Boreal fires combust vegetation and organic soils, which insulate permafrost, and as such deepen the seasonally thawed active layer and can lead to further carbon emissions to the atmosphere. Current understanding of the environmental drivers of post-fire thaw depth is limited but of critical importance. In addition, mapping thaw depth over fire scars may enable a better understanding of the spatial variability in post-fire responses of permafrost soils. We assessed the environmental drivers of post-fire thaw depth using field data from a fire scar in a larch-dominated forest in the continuous permafrost zone in eastern Siberia. Particularly, summer thaw depth was deeper in burned (mean=127.3 cm, standard deviation (SD) = 27.7 cm) than in unburned (98.1 cm, SD=26.9 cm) landscapes 1 year after the fire, yet the effect of fire was modulated by landscape and vegetation characteristics. We found deeper thaw in well-drained upland, in open and mature larch forest often intermixed with Scots pine, and in high-severity burns. The environmental drivers basal area, vegetation density, and burn depth explained 73.3 % of the measured thaw depth variability at the study sites. In addition, we evaluated the relationships between field-measured thaw depth and several remote sensing proxies. Albedo, the differenced normalized burn ratio (dNBR), and the pre-fire normalized difference vegetation index (NDVI) derived from Landsat 8 imagery together explained 66.3 % of the variability in field-measured thaw depth. Moreover, land surface temperature (LST) displayed particularly strong correlations with post-fire thaw depth (r=0.65, p<0.01). Based on these remote sensing proxies and multiple linear regression analysis, we estimated thaw depth over the entire fire scar. Our study reveals some of the governing processes of post-fire thaw depth development and shows the capability of Landsat imagery to estimate post-fire thaw depth at a landscape scale.

1 Introduction

Permafrost regions, covering Arctic and many boreal ecosystems, store approximately twice as much carbon as is currently in the atmosphere (Hugelius et al., 2014; Miner et al., 2022; Schuur et al., 2022). Permafrost ecosystems are rapidly changing because of climate warming (Natali et al., 2021; Schuur et al., 2022), and changes in fire regimes may accelerate permafrost thaw (Chen et al., 2021b; Descals et al., 2022; Gibson et al., 2018; Scholten et al., 2022). Fires are a major disturbance in permafrost landscapes (Descals et al., 2022; Holloway et al., 2020; Walker et al., 2019), influencing the carbon, energy, and water cycles (Jin et al., 2021; Larjavaara et al., 2017; Li et al., 2021; Liu et al., 2005; Minsley et al., 2016; Randerson et al., 2006).

Boreal fires combust parts of the aboveground biomass and organic soils (Dieleman et al., 2020; Veraverbeke et al., 2021; Walker et al., 2018a, 2020b). When this insulating layer is consumed, charcoal is left on the surface, reducing the albedo immediately after fire in the snow-free season (Rocha et al., 2012; Yoshikawa et al., 2003). This increases the absorption of shortwave radiation and the accumulated heat results in ground warming. Moreover, the combustion of vegetation canopies increases the exposure of the ground to solar radiation and reduces the evaporative cooling effect (French et al., 2016; Johnstone et al., 2010; Li et al., 2021; Liu et al., 2018; Shur and Jorgenson, 2007; Yoshikawa et al., 2003; Zhao et al., 2021). The increased ground heat flux deepens the seasonally thawed active layer above the frozen ground (Holloway et al., 2020; Zhang et al., 2023, 2015). Active layer thickening can enhance soil respiration and consequent emissions of the greenhouse gases carbon dioxide and methane for several decades (Genet et al., 2013; Köster et al., 2018; O'Donnell et al., 2011; Taş et al., 2014; Veraverbeke et al., 2021). Fires thus accelerate permafrost thaw, and although they are a natural component of boreal ecosystems (Kharuk et al., 2021), ongoing changes to the severity and extent of boreal fires (Descals et al., 2022; Turetsky et al., 2011b; Zheng et al., 2023) in combination with climate warming are projected to irreversibly change ice-rich permafrost environments (Chen et al., 2021a; Gibson et al., 2018; Holloway et al., 2020; Zhang et al., 2015).

Previous work has highlighted fires as an important driver of the thickening of the active layer (Brown et al., 2016; Gibson et al., 2018; Jafarov et al., 2013; Zhang et al., 2015). In addition, the active layer thaw depth is strongly controlled by vegetation type and density, organic soil content, soil texture, soil moisture content, topography, and drainage (Bai et al., 2018; Duguay et al., 2013; Fisher et al., 2016; Loranty et al., 2018b; Park et al., 2013). Interactions between these diverse environmental variables lead to a large heterogeneity of the depth of the active layer in the landscape (Holloway et al., 2020; Liu et al., 2012; Shiklomanov and Nelson, 2002). A higher canopy density reduces the amount of solar radiation reaching the ground, which in turn reduces the energy available for soil thawing (Fisher et al., 2016; Stuenzi et al., 2021a, b). In addition, vegetation type and density may alter soil moisture through evapotranspiration, thereby affecting the thermal conditions of the permafrost active layer (Fedorov et al., 2017; Iwahana et al., 2005; Juszak et al., 2016; Marsh et al., 2010). Topography also influences the landscape's vulnerability to burning (Brown et al., 2016). For example, lowlands with limited drainage and wet soils tend to burn with lower severity (Benscoter et al., 2011; Dillon et al., 2011; Holloway et al., 2020; Turetsky et al., 2011a). Thus, the active layer depth is governed by complex interactions between topography, vegetation, and fires (Brown et al., 2016). However, quantitative understanding of these effects over space and time, as well as the strength of underlying drivers, is needed given the intensification of Arctic and boreal fire regimes (Gibson et al., 2018; Michaelides et al., 2019; Natali et al., 2021).

Remote sensing provides information on the ground surface and shallow soils yet cannot directly measure the state of the subsurface. Nevertheless, remote sensing acquires information on various environmental drivers with a known influence on active layer depth or that are indicative of changes in it, such as vegetation status (Brown et al., 2016; Li et al., 2019), ground subsidence (Chen et al., 2020; Michaelides et al., 2019; Molan et al., 2018), land surface temperature (LST) (Bai et al., 2018; Park et al., 2016; Ran et al., 2022; Wen et al., 2022; Yi et al., 2018; Zorigt et al., 2020), surface albedo (Liu et al., 2018; Rogers et al., 2014; Zhao et al., 2021), and fire severity (Delcourt et al., 2021; Lentile et al., 2006; De Santis and Chuvieco, 2007). LST has shown potential to map and model active layer depth over large areas given its relationship to the ground thermal regime (Batbaatar et al., 2020; Hachem et al., 2012; Langer et al., 2010; Obu et al., 2019; Westermann et al., 2017). Obu et al. (2021) produced an active layer thickness (ALT) dataset derived from a remote-sensing-driven model based on Moderate Resolution Imaging Spectroradiometer (MODIS) LST merged with reanalysis near-surface air temperature data. Liu et al. (2024) integrated several MODIS-derived environmental variables, including LST, in a machine learning approach to predict gridded ALT at a spatial resolution of 1 km. However, no studies have used Landsat thermal remote sensing data so far to assess and map thaw depth in a post-fire environment at a spatial resolution of 100 m or less, primarily because of the lower temporal resolution of Landsat data in comparison with, for example, MODIS. The higher spatial resolution from Landsat, however, may offer benefits to map and understand post-fire effects on permafrost soils with more spatial detail, thereby allowing comparisons with field measurements.

In recent years, eastern Siberian larch forests underlain by permafrost have experienced extreme fire activity (Scholten et al., 2022; Talucci et al., 2022a, b; Xu et al., 2022; Zheng et al., 2023). In this study, we investigated post-fire thaw depth in a fire scar in a larch-dominated forest underlain by continuous permafrost in eastern Siberia using field and remote sensing data. Our objectives are twofold. First, we assessed which environmental drivers influence post-fire thaw depth. Second, we investigated which remote sensing proxies relate to field-measured thaw depth. Using these relationships, we mapped post-fire thaw depth over the entire fire scar.

2 Materials and methods

2.1 Study area

Our study area is a fire scar near the Yert rural locality, Republic of Sakha (also known as Yakutia), Russia (Fig. 1). The fire occurred between 30 June and 21 July 2018 and burned approximately 97 000 ha. The area is approximately 200 km west of Yakutsk, the capital of Yakutia. The mean annual air temperature in Yakutsk is −10.2 °C, with means of −42.6 °C in January and 18.7 °C in July, and the mean annual precipitation is 234 mm (Fedorov et al., 2017). The fire scar is in continuous permafrost terrain dominated by Cajander larch (Larix cajanderi), sometimes intermixed with the presence of Scots pine (Pinus sylvestris) and silver birch (Betula pendula), as well as the tall shrubs alder (Alnus spp.) and willow species (Salix spp.). The surface vegetation includes cowberry (Vaccinium vitis-idaea), bog blueberry (Vaccinium uliginosum), crowberry (Empetrum nigrum), rhododendron (Rhododendron dauricum), dog rose (Rosa acicularis), spirea shrub (Spiraea spp.), juniper (Juniperus spp.), fireweed (Epilobium angustifolium), moss (as Ceratodon purpureus and Aulacomnium palustre), and lichens (Delcourt et al., 2021). Cajander larch trees are deciduous needleleaf trees that are adapted to grow on permafrost terrain. They are physiologically and morphologically adapted to withstand the low temperatures, short growing seasons, and cryogenic processes (Abaimov, 2010; Berner et al., 2012). Cajander larch and Scots pine have evolved and adapted to frequent low-severity fires. The fire investigated in this study included a wide range of severities, from low-severity surface fires to high-severity stand-replacing fires (Delcourt et al., 2021).

Figure 1 Map of the study area in the Republic of Sakha, Russia. The fire scar is shown in the reddish brown colors in a Landsat 8 Operational Land Imager false-color composite (RGB 754) from 23 July 2019. The locations of burned and unburned field plots are shown on the map.

2.2 Field data

We collected field data in 13 burned and 7 unburned plots (Fig. 1) between 30 July and 8 August 2019, approximately 1 year after the fire (Delcourt et al., 2024). Plot selection aimed to cover a gradient in vegetation composition, fire severity, and landscape position within accessibility constraints. The plots were represented by quadrats of 30 m by 30 m to match the spatial resolution of the Landsat Operational Land Imager (OLI) optical imagery. We aimed to sample plots in areas where fire severity, vegetation, and landscape characteristics were relatively homogeneous (Delcourt et al., 2021).

For each plot, center coordinates and elevation were collected using a high-precision Global Positioning System handheld device (Trimble Geo 7X, GeoExplorer) with 1 m horizontal and vertical accuracy. The geolocation was post-processed using data from the closest reference station (Scripps Orbit and Permanent Array Center (SOPAC), Seismic Station Yakutsk), resulting in decimeter accuracy. A clinometer was used to determine the slope. We also assigned the topographic position of the plot in relation to its surroundings (upland, midslope, lowland) and plot-level site moisture classes following Johnstone et al. (2008). This approach assesses site moisture based on local topographic drainage, thereby accounting for soil texture and permafrost presence. The resulting ordinal scale consists of six site moisture classes ranging between dry and wet (xeric, subxeric, subxeric to mesic, mesic, mesic to subhygric, and subhygric). The six classes represent the potential moisture available for plant growth and should not be confused with actual measurements of soil moisture. This site moisture classification has been used extensively for fire studies in boreal North America (Dieleman et al., 2020; Walker et al., 2018a, b, 2020a).

We measured vegetation characteristics (stand age, basal area, vegetation density, and larch proportion) in a 30 m by 2 m transect in the north–south direction of the plot and with the plot center as the centroid. To characterize vegetation type and density, we inventoried every live and standing dead tree along the transect. Particularly, we recorded tree species, measured diameter at breast height (DBH; 1.3 m), and estimated basal area and stem density.

To characterize fire severity, we used the Geometrically structured Composite Burn Index (GeoCBI) protocol (De Santis and Chuvieco, 2009). GeoCBI is a visual field-based fire severity index in a continuous numeric scale ranging from 0 to 3. The GeoCBI estimates the cumulative impacts of fire over the vertical structure of forest stands and soil surface. Its field protocol is based on a hierarchical and multi-layered sampling design, splitting the plot into five different strata: substrates (ground surface, litter, duff); herbs, low shrubs, and trees less than 1 m; tall shrubs and trees of 1 to 5 m; intermediate trees of 5 to 20 m; and trees higher than 20 m. The strata are divided into subcategories that are evaluated independently according to various visible fire effects. These subcategories, for example, include changes in the color and condition of the soil, consumption of fuels, and char height on tree boles. The subcategories are rated with decimal values between 0, meaning no fire effect, and 3, reflecting high fire severity. Then, the scores of each stratum are obtained by averaging the scores for all criteria and strata are weighted by their fraction of coverage to derive the plot-level GeoCBI value. Previous studies have used this index to assess fire severity in boreal forests (Delcourt et al., 2021; Dieleman et al., 2020; Rogers et al., 2014).

We estimated burn depth in the burned plots by measuring the vertical distance between the top of the residual soil layer and the uppermost adventitious root on larch trees, hereafter referred to as adventitious root height (ARH), following the methodology described in Delcourt et al. (2021). Adventitious roots are fine lateral roots that develop at upper soil horizons as trees grow, responding to unfavorable temperature and moisture conditions deeper in the soil (Rogers et al., 2014). These roots remain visible on the tree trunks several years after a fire (Kajimoto, 2010; Kajimoto et al., 2003). The position of adventitious roots in the soil column can be used to estimate the height of the pre-fire soil surface (Boby et al., 2010; Kasischke and Johnstone, 2005; Rogers et al., 2014). In unburned plots, we determined the location of the highest adventitious root beneath the moss layer as well as the soil organic layer depth (SOL) and derived a linear relationship between SOL and adventitious root height above the mineral soil. We used this relationship in burned plots to reconstruct pre-fire SOL from measurements of ARH and residual SOL depth. Burn depth was then estimated as the difference between reconstructed pre-fire SOL depth and residual SOL depth.

Thaw depth is defined as a measure of the seasonally unfrozen portion of the soil column (the active layer) that lies atop perennially frozen ground (permafrost) (Rocha et al., 2012). The thaw depth was measured using a standard frost probe (stainless-steel rod), which was inserted into the ground to the depth of resistance by the frozen ground. We performed measurements twice (with 1 m increments) every 7.5 m along the belt transect. In some measurements, we were uncertain whether the depth resistance stemmed from frozen ground or was caused by rocks. In the case of uncertainty, we excluded the measurement from the analysis. Our approach resulted in a maximum of 10 thaw depth measurements per plot, which were then averaged to derive the plot-level thaw depth. All our measurements were made between 30 July and 8 August 2019, and we therefore assume negligible seasonal influences on thaw depth between plots.

2.3 Remote sensing data

Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) Collection 2 Level-2 data were used in this study (Table 1). The OLI reflective and the TIRS thermal bands have spatial resolutions of 30 and 100 m, but the latter are resampled and provided at 30 m for consistency purposes. We acquired post-fire images for 23 July 2019, the closest date to the field campaign. For the pre-fire imagery, we used a cloud-free image from 7 July 2016, near the anniversary date of the post-fire imagery. No cloud-free summer images were available from 2017. The acquisition timing of the pre-fire image aligns with recommendations of Key and Benson (2006), since no other disturbance in the landscape interfered with the subject fire between the summers of 2016 and 2017. Two scenes were mosaicked to cover the entire fire perimeter for both pre- and post-fire imagery (Fig. 1). We used the Landsat data to calculate the pre-fire normalized difference vegetation index (NDVI), post-fire albedo, land surface temperature (LST), and differenced normalized burn ratio (dNBR).

Table 1 Metadata for Landsat 8 imagery used in this study.

We calculated the pre-fire NDVI as follows:

$$\begin{array}{}\text{(1)}& \text{NDVI}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{5}}-{\mathit{\rho }}_{\mathrm{4}}}{{\mathit{\rho }}_{\mathrm{5}}+{\mathit{\rho }}_{\mathrm{4}}}},\end{array}$$

where ρ represents the reflectance of the Landsat 8 OLI bands 4 (0.64–0.67 µm) and 5 (0.85–0.88 µm). We included the pre-fire NDVI as a proxy for fuel availability.

To retrieve albedo (α), we used the narrow-to-broadband conversion formula from Liang (2001) adapted for Landsat 8 as follows (Naegeli et al., 2017):

$$\begin{array}{}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle }\mathit{\alpha }=& {\displaystyle }\mathrm{0.356}{\mathit{\rho }}_{\mathrm{2}}+\mathrm{0.130}{\mathit{\rho }}_{\mathrm{4}}+\mathrm{0.373}{\mathit{\rho }}_{\mathrm{5}}+\mathrm{0.085}{\mathit{\rho }}_{\mathrm{6}}\\ {\displaystyle }& {\displaystyle }+\mathrm{0.072}{\mathit{\rho }}_{\mathrm{7}}-\mathrm{0.0018},\end{array}\end{array}$$

where ρ is the reflectance of the Landsat 8 OLI bands 2 (0.45–0.51 µm), 4, 5, 6 (1.57–1.65 µm), and 7 (2.11–2.29 µm).

The dNBR was obtained from the bitemporal difference of the normalized burn ratio (NBR) (Eq. 3) between pre- and post-fire scenes (Eq. 4) (Epting et al., 2005; García and Caselles, 1991; Key and Benson, 2006).

$$\begin{array}{}\text{(3)}& {\displaystyle }\text{NBR}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{5}}-{\mathit{\rho }}_{\mathrm{7}}}{{\mathit{\rho }}_{\mathrm{5}}+{\mathit{\rho }}_{\mathrm{7}}}}\text{(4)}& {\displaystyle }\text{dNBR}={\text{NBR}}_{\text{pre-fire}}-{\text{NBR}}_{\text{post-fire}}\end{array}$$

The dNBR is an often-used fire severity index and theoretically can range between −2 and 2. In practice, however, most pixel values range between 0 and 1. A dNBR value of zero denotes an unburned pixel, and increasing dNBR values represent higher fire severity (Allen and Sorbel, 2008; Key and Benson, 2006).

For the LST retrieval, we opted for an approach based on the direct inverse solution of the radiative transfer equation (RTE) (Price, 1983). This is a physically based method that is often used with a single thermal band (Jiménez-Muñoz et al., 2009). Inverting and simplifying the RTE to a thermal band center at wavelength λ, the LST can be calculated as

$$\begin{array}{}\text{(5)}& \text{LST}={\displaystyle \frac{K_{\mathrm{2}}}{\ln \left({\displaystyle \frac{K_{\mathrm{1}}}{\frac{L_{\mathit{\lambda }}^{\text{sen}}-L_{\mathit{\lambda }}^{\uparrow }-{\mathit{\tau }}_{\mathit{\lambda }}\left(\mathrm{1}-{\mathit{\epsilon }}_{\mathit{\lambda }}\right)L_{\mathit{\lambda }}^{\downarrow }}{{\mathit{\tau }}_{\mathit{\lambda }}{\mathit{\epsilon }}_{\mathit{\lambda }}}}}+\mathrm{1}\right)}},\end{array}$$

where $L_{\mathit{\lambda }}^{\text{sen}}$ refers to the at-sensor spectral radiance of the corresponding thermal band. In this study, we used the Landsat 8 TIRS band 10 (10.60–11.19 µm). ε_λ represents the land surface emissivity (LSE). K₁ and K₂ are calibration constants for band 10 of the Landsat 8 TIRS (Ihlen and Zanter, 2019). $L_{\mathit{\lambda }}^{\downarrow }$ and $L_{\mathit{\lambda }}^{\uparrow }$ are the downwelling and upwelling atmospheric radiances and τ_λ is the atmospheric transmittance. These atmospheric parameters were derived from the Goddard Earth Observing System Version 5 (GEOS-5) Forward Processing for Instrument Teams (FP-IT) data, which are available as atmospheric auxiliary data for the Landsat 8 products (EROS, 2021).

LSE information is essential for the LST retrieval. LSE determines the efficiency of the surface to convert heat energy into radiant energy. We adopted the NDVI threshold method (NDVI^THM) of Sobrino et al. (2004) to estimate LSE, but we adapted the method to the specific characteristics of a fire scar. For that, we first computed the NDVI and the fractional vegetation cover (P_V) (Carlson and Ripley, 1997):

$$\begin{array}{}\text{(6)}& P_{\mathrm{V}}={\left({\displaystyle \frac{\text{NDVI}-{\text{NDVI}}_{\text{min}}}{{\text{NDVI}}_{\text{max}}-{\text{NDVI}}_{\text{min}}}}\right)}^{\mathrm{2}},\end{array}$$

where NDVI_min=0.2 and NDVI_max=0.5 (Sobrino et al., 2004, 2008; Sobrino and Raissouni, 2000). Then, we estimated LSE as follows. The spectral notation (λ) is omitted, since a single TIRS band was used here.

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}{\displaystyle }\mathit{\epsilon }& {\displaystyle }=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\mathrm{0.962}& \mathrm{0.0}\le \text{NDVI}<\mathrm{0.2}\\ \mathrm{0.990}{P}_{\mathrm{V}}+\mathrm{0.962}(\mathrm{1}-P_{\mathrm{V}})+d\mathit{\epsilon }& \mathrm{0.2}\le \text{NDVI}\le \mathrm{0.5}\\ \mathrm{0.990}& \text{NDVI}>\mathrm{0.5}\\ \mathrm{0.993}& \text{NDVI}<\mathrm{0.0}\end{array}\right.\end{array}\text{(8)}& {\displaystyle }d\mathit{\epsilon }=(\mathrm{1}-\mathrm{0.962})\mathrm{0.990}{F}^{\prime }(\mathrm{1}-P_{\mathrm{V}})\end{array}$$

In this approach, pixels with NDVI values between 0 and 0.2 were considered non-vegetated (P_V=0). Since our study focuses on a fire scar, an LSE value of 0.962 was assigned to these pixels. This value was derived from charcoal spectra from Veraverbeke et al. (2012, 2014) integrated over the Landsat 8 TIRS band 10 spectral response function. If the pixel had an NDVI value between 0.2 and 0.5, it was a mixture of charcoal and vegetation and LSE depended on the P_V value calculated according to Eq. (7). dε accounts for the cavity effect due to surface roughness (i.e., dε=0 for flat surfaces) and was approximated by Eq. (8), with an assumed shape factor F^′ of 0.55 following Sobrino et al. (1990). Pixels with NDVI values higher than 0.5 were fully vegetated (P_V=1), and we assumed an emissivity of 0.990. Finally, when the pixel had a negative NDVI value corresponding to water surface, an LSE value of 0.993 was assigned (Vanhellemont, 2020).

We calculated the average pixel values for all remote sensing metrics within a three-by-three pixel window centered on a given plot's centroid coordinates, which minimized the effects of potential satellite misregistration when comparing these remote sensing metrics with field data (Ahern et al., 1991). We also compared NDVI, albedo, dNBR, and LST values of burned pixels from within the fire perimeter with values of these variables of unburned pixels within a 2 km buffer from the fire perimeter.

2.4 Statistical approach

Our statistical analysis consisted of two major components. First, we assessed the environmental drivers of post-fire thaw depth using our field observations. Second, we evaluated how remote sensing variables can be used to predict field-measured thaw depth. For both analyses, we used an ordinary least-squares (OLS) multiple linear regression (MLR) model with thaw depth as the response variable using the “statsmodels” Python module (Seabold and Perktold, 2010). We hypothesized that post-fire thaw depth relates to variables associated with topography (topographic position, site moisture, slope, and elevation), vegetation (stand age, basal area, vegetation density, and larch proportion), and fire severity (GeoCBI and burn depth). We used the original data (i.e., non-standardized) and ranked the categorical variables topographic position and site moisture on an ordinal scale to include them in the linear analysis. The selection of the environmental variables for the MLR model was based on a stepwise regression approach testing both forward selection and backward elimination. We calculated the relative importance of the retained variables to the MLR model explanatory power using the “pingouin” Python module (Vallat, 2018), which is based on the R package “relaimpo” (Grömping, 2006). The R² from the MLR model is partitioned by averaging over orderings among regressors using simple unweighted averages as proposed by Lindeman et al. (1980), providing a decomposition of the explained model variance into non-negative contributions. Finally, we developed an MLR model using the remote sensing metrics (post-fire albedo, dNBR, LST, and pre-fire NDVI) as potential predictor variables to assess their capacity as proxies for post-fire thaw depth. The resulting MLR was then spatially extrapolated to derive a continuous map of post-fire thaw depth over the study area.

In addition to these statistical analyses, which included all field plots, we also focused separately on two burned–unburned plot pairs. While our study area captured large variability in landscape position, vegetation characteristics, and fire severity, these plots had the advantage of being spatially adjacent, with a road acting as the fire barrier. Because of this, the landscape and vegetation characteristics of these plots were very similar. As a result, these plot pairs represent interesting study cases in which the fire effect on thaw depth can be unequivocally separated from the landscape and vegetation influences.

3 Results

3.1 Environmental drivers of post-fire thaw depth

On average, summer thaw depth was deeper in burned (mean=127.3 cm, standard deviation (SD) = 27.7 cm) than in unburned (98.1 cm, SD=26.9 cm) plots (Fig. 2). This difference was statistically significant at p<0.05. Thaw depth was deeper in uplands (Figs. 3, B1a and d) and well-drained (Figs. 3 and B1b) landscape positions; in open (Figs. 3 and 4b), mature (Figs. 3 and B1e), and mixed forests, which included tree species other than Cajander larch (Figs. 3 and B1f); and in high-severity burns (Figs. 3, 4c, and B1g).

Figure 2 The 1-year post-fire summer thaw depth was deeper in burned plots compared to unburned plots. Each box ranges from the first to the third quartile. Whiskers extend to points that lie within 1.5 times the interquartile range. The median is indicated by the horizontal line and the mean by the black triangle.

Figure 3 Correlation matrix between field-measured environmental variables representing landscape, vegetation and fire severity characteristics, and thaw depth. GeoCBI is the Geometrically structured Composite Burn Index.

Figure 4 Relationships between field-measured thaw depth and (a) basal area, (b) vegetation density, and (c) burn depth. The shading around the regression lines indicates the 95 % confidence interval, while r is the Pearson correlation coefficient.

Based on the results of the stepwise regression for both forward selection and backward elimination, we retained basal area, vegetation density, and burn depth as predictor variables for the MLR (Fig. 4). The multiple linear regression results indicated that the selected environmental drivers together explained 73.3 % of the thaw depth variability. Burn depth was the main contributor to the explanatory power of the MLR model (Table 2).

Table 2 Summary of the multiple linear regression for the environmental drivers of thaw depth. SE is the standard error of the corresponding variable's coefficient. The R² column displays the coefficient of determination of each predictor variable separately as well as for the multiple linear regression model. The relative importance indicates the contribution of each independent variable to the multiple linear regression model R².

3.2 Remote sensing proxies of thaw depth

Within the fire perimeter, there was considerable spatial variability in albedo, dNBR, LST, and pre-fire NDVI (Figs. 5 and 6). Albedo values inside the fire perimeter had a mean of 0.10 (SD=0.02) compared to mean of 0.12 (SD=0.02) outside the perimeter (Figs. 5a and 6a). The dNBR values of our 13 burned field plots ranged between 0.37 and 0.80, with a mean of 0.59 (SD=0.13), while for the entire fire perimeter the mean was 0.41 (SD=0.27, Figs. 5b and 6b). In addition, LST values inside the fire scar were on average 5.24 K higher than the LST values outside the fire scar (Figs. 5c and 6c). However, for some parts of the fire scar, LST values exceeded the surrounding unburned pixels by approximately 20 K. For the pre-fire NDVI, the distributions of values inside and outside the fire perimeter were similar (Figs. 5d and 6d), with a mean NDVI value of 0.78 (SD=0.08) inside the perimeter and a mean NDVI value of 0.79 (SD=0.07) outside the perimeter.

Figure 5 Maps showing (a) albedo, (b) differenced normalized burn ratio, (c) land surface temperature, and (d) pre-fire normalized difference vegetation index derived from Landsat 8 imagery. In (b) and (d), the northwestern tip with “no data” indicates a part of the study area that was not covered by the Landsat 8 scenes.

Figure 6 Relative frequency distributions of remote sensing proxies derived from Landsat 8 imagery inside and outside the fire perimeters: (a) albedo, (b) differenced normalized burn ratio, (c) land surface temperature, and (d) pre-fire normalized difference vegetation index. The pixels outside the fire perimeter were within 2 km of the fire perimeter.

Thaw tended to be deeper in the plots with lower summer albedo (Fig. 7a), higher dNBR, (Fig. 7b), higher LST (Fig. 7c), and lower pre-fire NDVI (Fig. 7d). The correlation with LST was the strongest, while the correlations between thaw depth and dNBR, albedo, and pre-fire NDVI were not statistically significant at p<0.05 (Fig. C1).

Figure 7 Scatter plots and linear regression lines between field-measured thaw depth and remotely sensed (a) albedo, (b) differenced normalized burn ratio, (c) land surface temperature, and (d) pre-fire normalized difference vegetation index. The shading indicates the 95 % confidence interval for the regressions and r is the Pearson correlation coefficient.

Stepwise forward selection indicated that LST was the only variable to be selected for the linear regression, capturing 42.9 % of the variability in field-measured thaw depth. By contrast, stepwise backward elimination pointed dNBR, albedo, and pre-fire NDVI as the variables to be retained in the MLR model. Taken together, these three remote sensing proxies explained 66.3 % of field-measured thaw depth variability. dNBR was the most important predictor in this model configuration (Table 3).

Table 3 Summary of the multiple linear regression for the remote sensing proxies of thaw depth. SE is the standard error of the corresponding variable's coefficient. The R² column displays the coefficient of determination of each predictor variable separately as well as for the multiple linear regression model. The relative importance indicates the contribution of each independent variable to the multiple linear regression model R².

The MLR model for thaw depth based on dNBR, albedo, and pre-fire NDVI enabled spatially continuous thaw depth estimates over the study area, including uncertainty estimates (Fig. 8). We decided on this MLR configuration due to the higher R² when compared to LST alone. Post-fire thaw depth varied largely within the fire scar, for example between 83.49 cm (5th percentile) and 228.39 cm (95th percentile), with a mean of 139.72 cm (SD=43.76 cm). In contrast, for the pixels outside the fire perimeter, the estimated thaw depth varied between 66.59 cm (5th percentile) and 237.99 cm (95th percentile), with a mean of 125.08 cm (SD=52.15 cm) (Fig. 8c). On average, burned areas thus experienced 14.64 cm deeper summer thaw depth 1 year after the fire than unburned areas.

Figure 8 Maps of (a) estimated thaw depth, (b) its standard error, and the (c) relative frequency distribution for the pixels inside and outside the fire perimeter. The area outside the fire scar considers a 2 km buffer from the fire perimeter. The northwestern portion of the fire scar was not covered by the pre-fire Landsat imagery, so there is no estimate of thaw depth in this area. In (a) and (b), post-fire Landsat 8 Operational Land Imager (OLI) band 5 is used as the background image.

3.3 Differences between neighboring burned and unburned plots

Our field dataset contained two pairs of burned–unburned plots, which were separated by a road (Fig. 9). The landscape and vegetation characteristics were very similar between adjacent burned and unburned plots and represented mature larch-dominated forest in a subxeric to mesic environment (Table 4). The burned plots experienced medium to high fire severity with GeoCBI values of approximately 2.5 and burn depth from 9.89 to 10.52 cm. The thaw depth in the two burned plots was 53 and 73 cm deeper compared to the paired unburned control plots. In addition, LST values in the burned plots were between 3 and 5 K higher than in the unburned plots, and the burned plots experienced a clear drop in albedo.

Figure 9 Location of the two pairs of burned–unburned plots in the study area (a), including field photographs, which contain the plot identification labels (b–e).

Table 4 Comparison of environmental variables and remote sensing proxies of the adjacent burned–unburned plot pairs.

4 Discussion

4.1 Environmental drivers of post-fire thaw depth

Fire accelerates permafrost active layer thickening (Gibson et al., 2018). However, this process is also modulated by topographic and vegetation characteristics (Brown et al., 2015, 2016; Fisher et al., 2016). We found that well-drained terrain positions generally experience deeper thaw. Previous studies have reported contrasting findings concerning moisture effects. For example, Fisher et al. (2016) found deeper thaw with increased soil surface moisture. Other studies, in contrast, found thinner active layers for wet soils (Clayton et al., 2021; Zhang et al., 2005) because these wet soils require more energy for all ice to melt (i.e., latent heat of fusion), resulting in shallower active layers (Clayton et al., 2021). We found that forest density was a thaw-limiting factor, in agreement with Fisher et al. (2016), who identified canopy shading and evapotranspiration as mechanisms to explain this relationship between forest density and active layer depth. With denser canopies, less solar radiation can reach the ground, be absorbed, and be converted into heat, resulting in a shallower thaw (Fedorov et al., 2017; Iwahana et al., 2005; Juszak et al., 2014, 2016). Denser vegetation may also transpire more, thereby reducing the thermal conductivity in the soil (Fisher et al., 2016; Iwahana et al., 2005). Moreover, the presence of Cajander larch trees is closely coupled with permafrost (Herzschuh, 2020). In turn, permafrost helps to provide sufficient water to the trees by preventing it from draining away quickly (Zhang et al., 2011). Our results reflect this vegetation–permafrost interaction, where plots with fewer larch trees showed more pronounced thaw. Our field observations also showed that larch trees were prevalent and tend to dominate in mesic and hydric parts of the landscape, including in lowlands with shallow permafrost thaw. Conversely, Scots pine was more prevalent in the drier upland areas, which demonstrated deeper thaw (Fig. 3). Scots pine trees in our study region tend to occur on moderately warm and dry landscape positions such as ridges and hilltops. These locations were characterized by sandy soils and deeper thaw. In contrast, larch trees are dominant in lowland depressions of watersheds, characterized by cool and often waterlogged soils and shallower permafrost thaw (Isaev et al., 2010). Eastern Siberian larch trees are one of the only tree species that can successfully grow on permafrost soils with very shallow thaw depth. This is due to their ability to develop an adventitious rooting system (Herzschuh, 2020; Kajimoto, 2010). Larch and also Scots pine have evolved under periodic fire conditions, with the capacity to regenerate and grow after fire (Kharuk et al., 2021). These competitive advantages have led to the establishment of larch dominance in Siberian permafrost environments, and fire occurrence has been a prerequisite for this evolutionary process. However, increased permafrost degradation could lead to a shift in dominant tree species, with larch trees being replaced by pine trees and other tree species (Zhang et al., 2011).

Fire results in a rise in surface and soil temperature, which consequently increases thaw depth (Gibson et al., 2018; Jiang et al., 2015; Li et al., 2019; Nossov et al., 2013; Smith et al., 2015). Our results are in agreement with the overall pattern of rapid active layer thickening in the years immediately after the fire (Holloway et al., 2020). We found a mean difference of about 30 cm in field-measured summer thaw depth between burned and unburned plots 1 year after the fire, and our study thereby also fills a data gap for the effects of fires on active layer thickening in Siberia (Petrov et al., 2022; Ponomarev et al., 2020). Our case studies of two burned–unburned plot pairs divided by a fire barrier, though limited in sample size, helped us to separate the impact of fire on the active layer thaw from topographic and vegetation influences.

Increasing fire severity tends to increase thaw depth (Alexander et al., 2018; Holloway et al., 2020; Jafarov et al., 2013; Jiang et al., 2015; Li et al., 2019). The positive correlations between burn depth and thaw depth and between GeoCBI and thaw depth found in this study are consistent with these previous results. In addition, fire severity is often influenced by vegetation and topographic conditions. In our study, for example, the well-drained plots are on upland landscape positions with deeper thaw depths, while lowland areas are usually wetter with limited drainage and shallow thaw (Figs. 3 and B1). These wet lowland ecosystems also tend to burn with lower severity (Benscoter et al., 2011; Dillon et al., 2011; Holloway et al., 2020; Turetsky et al., 2011a). The thaw measurements in burned plots, as measured in the field, are conservative estimates as we did not account for the combustion of the organic layer in these measurements. We found a mean burn depth of 10.2 cm (SD=1.0 cm; Fig. 4c) in our burned plots. As a consequence, the thaw depth in reference to the pre-fire ground level is even deeper than compared to the post-fire ground level. This further elucidates the importance of fire severity, or burn depth, as a driver of post-fire permafrost thaw.

4.2 Remote sensing proxies of thaw depth

Charcoal residues after a fire are dark in color, decreasing the summer albedo in the fire scar. More solar radiation is absorbed by the darker surface, leading to a thickening of the thawed layer. Liu et al. (2018) and Zhao et al. (2021) assessed post-fire summer albedo change in the boreal forests of eastern Siberia and North America. They reported small immediate (1 year after the fire) albedo declines: a maximum decline of 0.02 in Siberia and a mean decline of 0.01 in North America. We found a mean drop of 0.02, equaling a relative decrease of 17 %, in summer albedo due to fire. The studies by Liu et al. (2018) and Zhao et al. (2021) used MODIS products with coarser spatial resolution, which may partly explain the differences from our study, which used 30 m Landsat data. While the individual relationship between albedo and thaw depth was not strong, albedo may still be a useful metric to assess fire influences on permafrost soils, in part because albedo is found to be a good proxy for aboveground fire severity (Fig. C1; Rogers et al., 2014). In line with previous studies (Barrett et al., 2010; Delcourt et al., 2021; Rogers et al., 2014; Veraverbeke et al., 2015), the dNBR was also strongly correlated with burn depth and as a result reasonably strongly related to thaw depth (r=0.41, p=0.07). Somewhat surprisingly, pre-fire NDVI did not perform well as a proxy for fuel density as it showed a low correlation with the field-measured pre-fire vegetation density (Fig. C1). This may be because the NDVI in Siberian larch forest is strongly influenced by understory vegetation, thereby confounding its value as a proxy for forest density (Bendavid et al., 2023; Loranty et al., 2018a).

We found that LST was the strongest stand-alone remotely sensed predictor of post-fire thaw depth. This can be explained by the fact that it captures variations in topography, vegetation, and fire severity, which all showed important influences on post-fire thaw depth. LST has been used before to model ALT (Obu et al., 2019; Ran et al., 2022; Wen et al., 2022) and to study the impacts of fire on biophysical processes (Liu et al., 2018; Zhao et al., 2021). Although these studies used LST with a different approach, using time series at a lower spatial resolution, we have built upon these studies by relating LST to post-fire thaw depth in a predictive model at landscape scale.

We demonstrated the potential of Landsat imagery to predict post-fire thaw depth, resulting in a spatially continuous map of post-fire thaw depth at landscape scale. The MLR model using the Landsat derived proxies was able to explain 66.3 % of the field-measured thaw depth variability. This performance is only slightly lower than the 73.3 % explained variability based on only field data. Despite this achievement, some limitations are worth mentioning. The Landsat 8 TIRS images have a native spatial resolution of 100 m but are distributed at 30 m after a non-reversible cubic convolution resampling. Even with data collected at a lower native spatial resolution than the optical bands, LST was still our best stand-alone remote sensing proxy for post-fire thaw depth. The availability of LST data at a higher spatial resolution would therefore likely be beneficial for the model performance. In addition, our field dataset was somewhat limited in size. Because of this, the MLR models should be seen as first-order statistical models of the relationships between field-measured and spaceborne environmental variables with thaw depth. As more measurements of post-fire permafrost thaw may become available over time (Loranty et al., 2021), more advanced statistical approaches, including machine learning techniques, could be explored to better understand and map post-fire permafrost thaw over larger spatiotemporal scales.

We used a snapshot of LST that temporally overlapped with the field data collection, using only one image. Alternatively, one could use all available images within the growing season and gap-fill with, e.g., coarse-spatial-resolution but high-temporal-resolution MODIS LST imagery. It might thereby be possible to apply the concept of growing or thawing degree days, which is often used to approximate thaw depth evolution and ALT (Hinkel and Nelson, 2003; Romanovsky and Osterkamp, 1997; Streletskiy et al., 2015, 2008). In addition, we measured and estimated thaw depth in the middle of the growing season. However, the most common variable of interest for intercomparison is ALT, i.e., the maximum depth at the end of the growing season (Zhang et al., 2021; Michaelides et al., 2019; WMO, 2016). The rate at which the active layer deepens decreases towards the end of the growing season. The thaw depth in August has been found to be approximately 6 % less than the ALT (Boike et al., 1998; Clayton et al., 2021; Hinkel et al., 2001; Michaelides et al., 2019; Zhang and Stamnes, 1998). Streletskiy et al. (2008) found that 95 % to 99 % of maximum thaw propagation is achieved by mid-August and that landscape-specific thaw depth patterns are related to landform elements that show spatial regularity at the landscape scale. We thus assumed that the spatial variability of the mid-season thaw depth is a reliable indicator of the spatial variability of ALT at the end of the growing season in the study region.

For future research, in addition to optical and thermal data, synthetic aperture radar (SAR) datasets can be used to study spatiotemporal variations of fire effects on permafrost landscapes. SAR data enable the measurement of various physical properties of fire scars, including surface roughness and near-surface soil moisture. Moreover, the interferometric SAR (InSAR) technique allows precise measurements of ground movement, with millimeter vertical accuracy (Strozzi et al., 2018). The effectiveness of this approach has been demonstrated in previous studies conducted in permafrost environments in Alaska and Siberia (Michaelides et al., 2019; Yanagiya et al., 2023; Yanagiya and Furuya, 2020). Future work could integrate fire-induced permafrost subsidence from InSAR with optical and thermal data to predict post-fire thaw depth based on multi-sensor satellite products.

5 Conclusions

Boreal fires combust parts or all of the aboveground vegetation and soil organic layer. The resulting loss of insulation of the permafrost soil leads to an increased ground heat flux, promoting permafrost thaw. With the escalating extent and severity of fires in Siberia, there is a growing need to better understand the processes that govern fire-induced permafrost thaw, including spatial estimations of thaw depth. Here, we evaluated the impact of fire, landscape, and vegetation characteristics on post-fire thaw depth. In addition, we assessed the capacity of remote sensing proxies to estimate post-fire thaw depth at a landscape scale. Using field data from a recent fire scar in eastern Siberian larch-dominated forests, we confirmed that increasing fire severity indeed induces deeper thaw. This process was, however, mediated by landscape and vegetation characteristics. Thaw depth was deeper in well-drained uplands, characterized by open and mature forests, which included mixed stands of Cajander larch and Scots pine. The environmental drivers basal area, vegetation density, and burn depth together explained 73.3 % of the thaw depth variability. We explored the use of Landsat-8-derived metrics as proxies for post-fire thaw depth. A multiple linear regression model ingesting albedo, dNBR, and pre-fire NDVI explained 66.3 % of the field-measured variability in thaw depth. LST individually explained 42.9 % of this variability. The correlation between LST and thaw depth was particularly strong, highlighting the potential of LST as a proxy for post-fire thaw depth. We estimated post-fire thaw depth at a landscape scale by spatially extrapolating the regression model of the Landsat data. Future work could further investigate the consistency of our findings across space and time to better understand the impact of fire on permafrost ecosystems across the Arctic–boreal region.

Appendix A

Table A1 Summary of the remote sensing metrics used in the study. ρ represents the reflectance of the Landsat 8 Operational Land Imager (OLI) bands. $L_{\mathit{\lambda }}^{\text{sen}}$ refers to the Landsat 8 Thermal Infrared Sensor (TIRS) band 10 radiance. K₁ and K₂ are calibration constants. $L_{\mathit{\lambda }}^{\downarrow }$ and $L_{\mathit{\lambda }}^{\uparrow }$ are the downwelling and upwelling atmospheric radiances and τ_λ is the atmospheric transmittance.

Appendix B

Figure B1 Relationships between thaw depth and (a) topographic position, (b) site moisture, (c) slope, (d) elevation, (e) stand age, (f) larch proportion, and (g) Geometrically structured Composite Burn Index. In (a) and (b), each box ranges from the first to the third quartile. Whiskers extend to points that lie within 1.5 times the interquartile range. The median is indicated by the horizontal line and the mean by the black triangle. For site moisture, classes were grouped together for better visualization. In the regression lines, shading indicates the 95 % confidence interval and r is the Pearson correlation coefficient.

Appendix C

Figure C1 Correlation matrix of selected environmental drivers, remote sensing proxies, and thaw depth. GeoCBI is the Geometrically structured Composite Burn Index. LST is the land surface temperature. dNBR refers to the differenced normalized burn ratio and NDVI to the normalized difference vegetation index.