The E ﬀ ect of Post-Fire Disturbances on a Seasonally Thawed Layer in the Permafrost Larch Forests of Central Siberia

: We examined and simulated the consequences of the degradation of the litter and the moss–lichen layer after ﬁre impact, which could a ﬀ ect the seasonal temperature of the soil and the depth of the seasonally thawed layer (STL) in the permafrost zone. According to the analysis of satellite imagery for 2000 to 2019, the ﬁre-disturbed area in the region of interest amounted to 20%. The main aims of the study included quantitative evaluation of the variation range of summer temperature anomalies at ﬁre-damaged plots, summarizing the statistical norm of the STL depending on natural conditions, and numerical simulation of the response of the STL. Using Terra and Aqua / MODIS imagery, we analyzed surface temperature (in bands of λ = 10.780–11.280 and 11.770–12.270 µ m) coupled with the normalized di ﬀ erence vegetation index (NDVI) for non-disturbed and ﬁre-damaged sites under the same natural conditions of larch forests in Central Siberia. Heat transfer, freezing and thawing processes were numerically simulated for two extreme cases of soil conditions: dry soil and water-saturated soil. The model was also applied to soil with non-homogeneous water content. As input parameters, we used data on the properties of cryogenic soils collected in larch forests ( Larix gmelinii ) in the ﬂat-mountainous taiga region of the Evenkia (Central Siberia). For post-ﬁre plots, surface temperature anomalies observed during summer months remained signiﬁcant for more than 15–20 years after ﬁre impact, while the NDVI values were restored to the statistical norm within 7–10 years of the ﬁre. According to the results of numerical simulation, the thickness of the STL could show a 30–50% increase compared to the statistical norm. In the ﬁrst approximation, we showed the annual soil temperature dynamics at various depths in disturbed and non-disturbed plots.


Introduction
In the permafrost zone of Siberia, significant parts of larch forests have been damaged by wildfires over the last two decades. According to satellite monitoring, fires damage at least 1% of the forested area per year, which determines a significant "cumulative" effect in terms of ecosystem disturbance [1]. There is an increasing trend in fire frequency [2], which will result in a significant increase in the disturbance of Siberian ecosystems in the future.
Along with the influence of climate change, the impact of massive wildfires has caused abnormal processes in permafrost soils in terms of the temperature balance, annual thermal regime and seasonal dynamics of the seasonally thawed layer (STL) [3][4][5]. Due to the destruction of vegetation cover in fire-damaged areas, the temperature regime of the surface and soil changes significantly. These changes determine the dynamics of STL, the thickness of which is strongly linked to the condition of vegetation cover, including the moss-lichen layer and litter [6][7][8].
Modern estimates allow us to state the significance of processes in the permafrost soils at the level of the global climate system [6,[9][10][11]. Locally, permafrost degradation may threaten the stability of the existing ecosystem [12,13]. At the same time, positive feedback has also been recorded in the permafrost area, namely, the radial tree growth coinciding with the increase in STL [12,14]. While field measurements of the temperature regime of permafrost soils are scarce and broadly dispersed due to the remoteness of the territory, a numerical simulation based on indirect features and generalizations of the available data can predict some general trends of the thermal regime and its dynamics [13,15,16].
The impacts of post-fire disturbances on the STL have been well studied in the permafrost zones of Alaska [17,18], Canada [19][20][21] and Siberia [3][4][5]. However, for the larch forests of Siberia, this issue has become more relevant with the increase in the number of fires and burned areas in recent decades.
The aims of this study were (i) to perform a quantitative assessment of summer temperature anomalies in fire-damaged plots of the permafrost zone using remote sensing data; (ii) to summarize the average rate of change in the STL based on field data and publications regarding the area of interest (Evenkia, Krasnoyarsk region, Central Siberia), and (iii) to simulate the annual dynamics of soil temperature and STL thickness so as to evaluate the rate of abnormal thawing after a fire impact in larch forests in the permafrost area.

Study Area
The area of interest is the Central Siberian flat-mountainous taiga region of the boreal taiga zone (60- The study area is located in a zone of continuous permafrost [22]. This region is characterized by a cold continental climate with a mean annual temperature of −9.5 • C and a mean monthly temperature ranging from −36 • C in January to 16 • C in July. The region's average annual precipitation is 300-350 mm, about 30-40% of which falls as snow, which averages 40-50 cm per year. The growing season lasts for 69-80 days [23]. Larch (Larix gmelinii Rupr. Rupr.) is the dominant species in the forest stands. Mosses, lichens and vascular species represent the vegetation cover [24].
The soil cover is represented by several types of Cryosols [25], and the permafrost occurs at a depth of <1 m from the surface. Such soils are characterized by cryoturbation and contain a large amount of poorly decomposed organic material in the upper soil layers [26]. Soils were characterized by their shallow (20-50 cm) depths, slight or neutral acidity, coarse texture (high gravel content) and medium to high clay contents. The soil-developing processes (cryoturbation, freezing and thawing) were variable, and their contributions to STL dynamics could be highly uncertain. Non-disturbed soils were characterized by the following soil profile formula: T (Oao)-O(F+H)-CR-C ⊥ [27,28].
As a result of fire disturbance, the soil profiles lost their upper areas, and the formula of the soil profile was transformed to Opir-O(F+H)pir-CR-C [27,28]. Permafrost soil layer characteristics for control plots and postpyrogenic plots were summarized in Table 1 from [28][29][30][31][32][33]. More details of the permafrost soil layers are discussed in Section 2.3.
As a result of fire disturbance, the soil profiles lost their upper areas, and the formula of the soil profile was transformed to Opir-O(F+H)pir-CR-C [27,28]. Permafrost soil layer characteristics for control plots and postpyrogenic plots were summarized in Table 1 from [28][29][30][31][32][33]. More details of the permafrost soil layers are discussed in Section 2.3. (c) Types of permafrost: 1-continuous (90-100% of territory); 2-discontinuous (50-90%); 3-sporadic (10-50%); 4-isolated patches (0-10%); 5-study area. We used the multispectral data from Terra and Aqua/MODIS for 2000-2019 to assess wildfire distribution over the region, and to select post-fire plots of different ages (Figure 2a). We used our own (V.N. Sukachev Institute of Forest, Krasnoyarsk, Russia) wildfire database for Siberia created using remote sensing data obtained in 1996-2019. The information on wildfires was stored as a geographic information system (GIS) polygon layer, with attributive data on coordinates, dates of fires and daily burn dynamics, as well as area burnt [34].
This wildfire database was generated using a multistep process, including: (1) contextual active fire detection; (2) the creation of fire polygons from adjacent fire pixels; and (3) the correction of resulting polygons. The processing chain for MODIS data was based on the approach of Giglio [35], incorporating several adjustments in background characterization and detection probability estimation [36]. The algorithm used by the Institute of Forest differs mainly in its approach of aggregating fire detections into fire polygons and subsequently using a correction procedure. This correction procedure was based on comparisons between the burned areas measured using low-resolution data (MODIS) and moderate resolution data (Landsat) for the test sample of fires (~5% of all annually detected fires) and obtaining the linear regression equations for several fire size classes. Then, these equations were used to correct fire polygons for the rest of the fires detected using low-resolution data [34]. High-resolution imagery (15-30 m) from Landsat/ETM/OLI (Enhanced Thematic Mapper/Operational Land Imager) and Sentinel-2 were used for the precise estimation of post-fire disturbance geometry. We used high-resolution imagery pre-processing data to outline fire scars using the geoinformation technique (GIS). We selected large-scale post-fire plots only, using the threshold S > 50,000 ha. These polygons were overlaid on the Terra/MODIS data, and time series were obtained for the post-fire pixels within the disturbed plots, as well as data on the territories nearest to each plot, which were used as "control". Finally, 50 post-fire plots of different successional stages were selected. The coordinates and dates of fires, as well as their linear characteristics (fire polygons), were used to further link the post-fire sites from the satellite images of different years.
We analyzed the surface temperature and normalized difference vegetation index (NDVI) for a series of post-fire plots and their neighboring undisturbed areas. To analyze NDVI, we operated with reflectances measured in MODIS bands 1 (λ 1 = 0.620-0.670 µm) and 2 (λ 2 = 0.841-0.876 µm) (product MOD09GQ [37]), and to analyze surface temperature anomalies we used the MOD11A1 product [38]. We used 10-day averaged data across the entire set of initial data, with a special focus on the ecosystem's recovery succession stages. A monthly data average procedure, as well as a procedure of spatial averaging within a post-fire plot, were applied.
For the selected post-fire plots, we collected datasets on the temperature of the surface and NDVI dynamics for the growing period (June-September) according to the satellite imagery of different years. Next, data for post-fire plots that were 1, 2, 5 and 10 years old were compared with values from the nearest non-disturbed areas, which were used as controls. We used the deviation (in percentage) of the temperature and NDVI values for the post-fire plots relative to the control plots' values as a measure of fire-caused anomalies. Finally, according to satellite data pre-processing, we quantified the changes in NDVI and thermal characteristics, and we also estimated the time lag and restoration rate of the NDVI and surface temperature after a fire disturbance. However, we did not calculate the broadband albedo and the emissivity coefficient, since a limited set of spectral channels was used. For the numerical simulation stage, we adopted these coefficients as described in [39]. The broadband albedo and the emissivity of the post-fire plots (Figure 2b) were set to 0.13 and 0.9, respectively, while these values were set to 0.18 and 0.98 for the non-damaged larch forests. For the selected post-fire plots, we collected datasets on the temperature of the surface and NDVI dynamics for the growing period (June-September) according to the satellite imagery of different years. Next, data for post-fire plots that were 1, 2, 5 and 10 years old were compared with values from the nearest non-disturbed areas, which were used as controls. We used the deviation (in percentage)

Data on Soils and STL Influenced by Natural Conditions
To characterize STL variations under the natural conditions of the region (Figure 2c), we collected all the available published studies conducted in Tura (Evenkia) over the period 1976-2015. We found 213 publications, and 51 sources (1062 observations, Supplementary Materials Table S1) were applied for the analysis. For further analysis, we used the data on STL measurements, tree stand age (0-241 years) and soil organic layer (litter) thickness, and we included the monthly variation of the STL. A multi-factor linear regression model was used for the assessment of the influence of different factors on the STL. In addition, the moss-lichen layer restoration versus the post-fire stage was evaluated, which was compared to the dynamics of vegetation index restoration obtained from the remote sensing data.
Ecosystem soils are characterized by the presence of well-developed organogenic horizons (peat or raw humus), with thicknesses varying from 5 to 20 cm [28,29]. The bulk density of organogenic horizons varies in the range of 60 to 200 kg/m 3 , depending on the degree of decomposition of vegetation residues, while the bulk density is 1100-1600 kg/m 3 in mineral horizons [28,29]. Post-fire soils are characterized by the upper layer of the organogenic horizon (Oao) being burned down, and unburned charred particles penetrating the underlying organic layer (O(F+H)), which is relatively loose. Thus, the bulk density of the upper horizons increases, and the porosity decreases after wildfire impact. In addition, the presence of charred particles leads to a change in the specific heat of the pyrogenic horizon [30], which was summarized in Table 1.
Porosity (P) was calculated by the formula : where D b is the bulk density and D is the density of the solid phase of the soil. These parameters, and specific heat (J/(kg·K)), were averaged according to the data of [28,[30][31][32][33], satisfying the soil conditions of the study area. Post-fire disturbance affects the STL via the increased exposure of land surfaces to solar radiation, which leads to the heating and drying of the land surface, as well as via the reduced thicknesses of the moss and surface litter layers, which provokes an increased conductance of heat from the surface into the soil profile. Naturally, these factors remained significant for many years after the fire. However, to simulate STL anomalies immediately after the fire, we used the approximation of the complete degradation of the moss-lichen cover, and the assumption of 100% mortality of the forest stand ( Figure 2b). Such a scenario is highly feasible in the fire-damaged larch forests of the permafrost zone of Siberia.
In our approximation, we considered that the differences between the fire-damaged and non-disturbed sites were significant for the first two organic soil horizons only. This is unique to the pyrogenic nature of the damage to the soil. This is how the soil is damaged by wildfire.
The model used the data from Table 1 for the non-damaged soils. At the same time, the complete absence of the upper organic layer (Opir) was assumed in the model for fire-damaged soils. Other features of the pyrogenic transformation of soil horizons were not taken into account in the employed version of the model.

Numerical Simulation of The STL's Thickness under the Conditions of the Post-Fire Plots
Next, we simulated the annual dynamics of soil temperature and STL depth in order to evaluate the rate of abnormality after a fire disturbance in permafrost larch forests. As noted above, the simulation was carried out under the assumption that the plots were damaged by fire recently, with complete degradation of the moss-lichen cover and decay of the forest stand.
The heat transport equation defined for the field of temperature, or enthalpy, is the general fundamental formulation used in mathematical models for heat transfer in soils [40]. In the present study, the temperature was chosen as a sought variable. This makes it much easier to calculate the thermophysical properties of the soil based on its known compounds. We chose the method of freezing and melting, in which the liquid and solid states of water are considered as separate components of Forests 2020, 11, 790 7 of 18 the medium. The general method for describing the phase transitions in dispersed multiphase media (see, for instance, [41]) is to relate the phase transition rate to the rate of the transfer of heat to the fragments of decaying phase. Considering that the heat transfer rate is proportional to the difference between the mean temperatures of the medium, the equation for the mass fraction of the liquid phase of water, α l , is written as follows: where k is the constant (1/(s · K)) characterizing the intensity of heat transfer, T is the temperature in Kelvin (K), and T 0 is the phase transition temperature, T 0 = 273.15 K. With the above relation for the phase transition rate, the heat transfer equation takes the following form where c p (T) is heat capacity (J/(kg·K)), ρ is density (kg/m 3 ), λ (T) is thermal conductivity (W/(m·K)), y w is the volumetric fraction of total soil moisture (dimensionless), α l is water mass per unit of soil moisture mass (dimensionless), and L is latent heat (J/kg). The boundary conditions for Equation (3), on the lower boundary as well as on the side boundaries, were all zero (the value of heat flux). The correctness of this condition was ensured by the location of the boundaries: the lower one is deep enough and the side ones are far enough from the damaged area, in the region where the condition of the soil is a function of the vertical coordinate only. The upper boundary conditions for (3) were determined by the presence of short-wave (solar) and long-wave (atmospheric) radiation, and convective heat exchanges with the air [42,43]. Solar radiation includes the direct and diffused components. The model described in the work of [43] was used to calculate the direct solar radiation.
The convective heat exchange rate between the atmosphere and solid surface is determined as where T air (t) is the air temperature, T sfc is the temperature of the surface of the solid material (ground or snow), and h sfc is the heat transfer coefficient, which depends on the wind speed and the type of terrain. The necessary expressions were taken from studies [44,45]. The annual data for the air temperature, wind speed and snow depth were taken from the observations of the weather station near Tura (Evenkia, Central Siberia). We used a long-term series of meteorological data from the free-access data banks Climatic Research Unit [46], Weather Archive [47], and NNDC Climate Data [48]. The frequency of updating the data on temperature was 1 h; for calculating the state in winter, the data on the depth of snow cover were averaged over months.
The thermal conductivity (λ (T)) for each soil horizon, depending on the water and ice content, was calculated using the models described in works [44,[49][50][51][52]. In some studies [44,50], a large amount of observational data was generalized to formulate the dependencies of the thermal conductivity coefficients on the dry soil density, and a volumetric fraction of water and ice. Other studies [49,51,52] used the expressions of the thermal conductivity coefficient, which were based on the model of Johansen [53] that proposes a mixing rule for these coefficients' values in dry (y w~0 .1) and water-saturated (ice-saturated) (y w~0 .9) soil: Forests 2020, 11, 790 8 of 18 where K is the Kersten number, which depends on the water content in the soil. The thermal conductivity of saturated soil depends on the liquid water mass fraction, changing within the limits set by the thermal conductivity of water-saturated and ice-saturated soil.
In the modeling, we considered two extreme conditions of water content in soil-dry soil and water-saturated soil, wherein the water occupies 10% and 90% of pore volume, correspondingly. For the types of soils with non-homogeneous water contents in the model, the water contents for each soil horizon have been varied (Table 2). Table 2. Volumetric water content (Q w ) in summer and ice content (Q i ) in winter for soil with non-homogenous water content (NHS), summarized from [45,51]. In the modeling, we considered three model variants of moisture distribution in the soil. In the case of dry soil (DS), where water fills 10% of the soil pore volume (see porosity (Q) in the Table 1), we used volumetric water content Qw = 0.1·Q, which was similar in all soil horizons. In cases of water-saturated soil (WS), where 90% of the soil pores are filled, we used Qw = 0.9·Q, which was not differentiated between soil horizons. The boundary cases DS and WS were considered to study the behavior of the mathematical model.

Soil
The case of non-homogeneous soil water content (NHS) was considered as a realistic case, accounted for by the field observations of the authors. In the case of NHS, we summarized volumetric water content according to limits of Pavlov's [45] and Ekici et al.'s [51] models ( Table 2).
The generalized characteristics of permafrost soils from the area of interest were used in the model. For the typical four horizons (see soil formula above), we assessed such data as soil texture and soil density, porosity, and specific heat (see Table 1 [49] were used for the next two loamy soil horizons (see Table 1).  Figure 1b,c). Finally, we calculated the total forested areas damaged by wildfires, which were 1.14 Mha and 9.89 Mha for 2000-2009 and 2010-2019, respectively. According to that data, the value of the relative burned area (RBA) index, which is the ratio of total burned area per year to the total forested area of the region, increased almost 10-fold during the last decade (from 0.2% up to 1.9% per year). Disturbances in the form of fire-scars are well-observed features of Terra/MODIS images, which are classified over a long time (see Figure 2a).

Characteristics of Fire-Disturbed Areas from Remote Sensing Data
For the 1-year-old burnt areas, the NDVI was 40-50% lower than for the non-disturbed plots ( Table 3). The deviation of NDVI abnormality (21 ± 7%) decreased more than twofold after a 5-year Forests 2020, 11, 790 9 of 18 post-fire recovery period. Finally, the 10-year-old burnt larch plots did not differ significantly from the non-disturbed areas in terms of NDVI, due to the restoration of on-ground vegetation cover and partial reforestation of disturbed plots. For these cases, the mean deviation did not exceed 9%, with a significant dispersion of σ = 5%. However, while anomalies in vegetation cover decreased over the 5-10 years after burning, the process of complete tree stand restoration was much slower, and took up to 50 years or longer [12,55]. The surface temperature anomalies in larch forests were significant during the 10-15 years following the wildfire (Table 3, Figure 3). According to the 10-day averaged Terra/MODIS data for 1-year-old burnt areas, the mid-summer maximum deviation of surface temperature was 7.0 ± 1.5 • C, which was 20-40% higher than those of non-disturbed plots. NDVI anomalies correlated with the data for mid-summer temperature anomalies (Figure 3a). According to Figure 3b, the rate of loss of temperature anomalies is 2.5 times lower than the rate of vegetation restoration in terms of NDVI (exponential approximation coefficients are −0.08 and −0.2, respectively). The speed of loss of NDVI deviation, as well as the temperature deviation, should be mediated by fire severity, drainage, stand density, etc. However, in the current study, we did not attempt to quantify such relations.
Forests 2020, 11, x FOR PEER REVIEW 9 of 18 5-10 years after burning, the process of complete tree stand restoration was much slower, and took up to 50 years or longer [12,55]. The surface temperature anomalies in larch forests were significant during the 10-15 years following the wildfire (Table 3, Figure 3). According to the 10-day averaged Terra/MODIS data for 1year-old burnt areas, the mid-summer maximum deviation of surface temperature was 7.0 ± 1.5 °C, which was 20-40% higher than those of non-disturbed plots. NDVI anomalies correlated with the data for mid-summer temperature anomalies (Figure 3a). According to Figure 3b, the rate of loss of temperature anomalies is 2.5 times lower than the rate of vegetation restoration in terms of NDVI (exponential approximation coefficients are −0.08 and −0.2, respectively). The speed of loss of NDVI deviation, as well as the temperature deviation, should be mediated by fire severity, drainage, stand density, etc. However, in the current study, we did not attempt to quantify such relations.
(a) (b) Figure 3. NDVI and temperature from 50 selected fire-damaged plots relative to non-disturbed areas averaged for June-August: (a) dependence between anomalies of temperature and NDVI; (b) dynamics of NDVI anomaly (1) and surface temperature anomalies (2) for 10 years after a fire.

Regression Model of STL
All available measurements of STL versus tree stand age (see Table S1) are summarized on a graph (Figure 4a). The linear regression coefficients varied, since the depth of the STL varies naturally during the growing season (from 46 ± 30 to 115 ± 50 cm (monthly mean ± SD)). We received a family of regression lines for each month in the June-September period (Equation (5)). Here, burnt areas are assumed to be larch stands less than 30 years old, and non-disturbed larch stands are assumed to be older than 60 years (Figure 4a). . NDVI and temperature from 50 selected fire-damaged plots relative to non-disturbed areas averaged for June-August: (a) dependence between anomalies of temperature and NDVI; (b) dynamics of NDVI anomaly (1) and surface temperature anomalies (2) for 10 years after a fire.

Regression Model of STL
All available measurements of STL versus tree stand age (see Table S1) are summarized on a graph (Figure 4a). The linear regression coefficients varied, since the depth of the STL varies naturally during the growing season (from 46 ± 30 to 115 ± 50 cm (monthly mean ± SD)). We received a family of regression lines for each month in the June-September period (Equation (5)). Here, burnt areas are assumed to be larch stands less than 30 years old, and non-disturbed larch stands are assumed to be older than 60 years (Figure 4a).  Tree stand age was a proxy for stand density and moss-lichen layer, as well as for litter layer depth. Thus, STL thickness decreased according to the tree stand age. Besides, the same dataset allowed the quantifying of the rate of moss-lichen layer recovery at the disturbed plots after fire impact (Figure 4b).
where A is tree stand age (years), L is soil organic layer (litter) (cm), t is time step (month of the growing period), and a, b, c and d are coefficients of regression (Table 4).  Tree stand age was a proxy for stand density and moss-lichen layer, as well as for litter layer depth. Thus, STL thickness decreased according to the tree stand age. Besides, the same dataset allowed the quantifying of the rate of moss-lichen layer recovery at the disturbed plots after fire impact (Figure 4b).

Simulated Annual Dynamics of Soil Temperature and STL Depth
The dynamics of STL thickness were assessed using the two above-mentioned extreme cases of soil conditions. The models were also applied for the more reliable case of permafrost soil conditions with non-homogeneous water content in its horizons (see Tables 1 and 2).
The soil profile thawing simulation proceeded via four phases (Figure 5c): A-active thawing phase; B-phase of the decreasing of the rate of thawing, the extreme of STL thickness; C-phase of upper soil horizon freezing at the transition of mean air temperature from a positive to negative value, when a part of the STL is still not frozen; D-phase of the transition state (freezing process) in a wide layer of the soil profile. Phase D cannot be registered in absolutely dry soil (Figure 5a,b), while its duration in over-wetted soils is determined by the environmental conditions and may last 20 days (see mark "1" in Figure 5d). The dynamics of STL thickness were assessed using the two above-mentioned extreme cases of soil conditions. The models were also applied for the more reliable case of permafrost soil conditions with non-homogeneous water content in its horizons (see Tables 1 and 2).
The soil profile thawing simulation proceeded via four phases (Figure 5c): А-active thawing phase; B-phase of the decreasing of the rate of thawing, the extreme of STL thickness; C-phase of upper soil horizon freezing at the transition of mean air temperature from a positive to negative value, when a part of the STL is still not frozen; D-phase of the transition state (freezing process) in a wide layer of the soil profile. Phase D cannot be registered in absolutely dry soil (Figure 5a,b), while its duration in over-wetted soils is determined by the environmental conditions and may last 20 days (see mark "1" in Figure 5d).  Pavlov (2008) [45]; (e) Soil with non-homogenous water content in the model developed by Ekici et al. (2014) [51]; (f) Soil with non-homogenous water content in the model built by Pavlov (2008) [45]. Results are from a model for the 365th day of the year, as the first 365 countings were used for system debugging and self-calibration. Figure 6 shows the results of the modeling of annual soil temperature dynamics at the depths of the root layer and below.
(2008) [45]; (c) Water-saturated soil (Qw = 0.9·Q, Qw not differentiated for soil horizons) and model developed by Ekici et al. (2014) [51]; (d) Water-saturated soil in the model built by Pavlov (2008) [45]; (e) Soil with non-homogenous water content in the model developed by Ekici et al. (2014) [51]; (f) Soil with non-homogenous water content in the model built by Pavlov (2008) [45]. Results are from a model for the 365th day of the year, as the first 365 countings were used for system debugging and self-calibration. Figure 6 shows the results of the modeling of annual soil temperature dynamics at the depths of the root layer and below.

Fire Disturbances and Post-Fire Effects
There is a significant trend of larch forests disturbance by wildfire, which is about 20% for the last two decades. Compared to the RBA for the Evenkia region, which increased almost 10-fold over the last decade (from 0.2% up to 1.9% per year), the RBA for the total Siberian region was estimated as 0.51-1.19% for different years [56], and the RBA for the forests of western Canada was 0.56% [57]. Thus, a significant "cumulative" effect of the fires provokes the large-scale anomalies in vegetation cover and the thermal balance of the soil surface (Table 3).
It was found that changes in the thermal regime of post-fire plots under the conditions of permafrost are accompanied by an abnormal increase in the average temperature of the surface and upper soil layer in the summer period. It was shown that NDVI values are restored within 7-10 years of the fire (Figure 3b), which is correlated with the rate of post-fire lichen layer recovery (Figure 4b). However, the partial regeneration of vegetation cover did not compensate for the changes, which lead to long-term surface temperature abnormalities. The thermal balance of post-fire sites with disturbed vegetation cover remains affected for more than 15-20 years, which is consistent with

Fire Disturbances and Post-Fire Effects
There is a significant trend of larch forests disturbance by wildfire, which is about 20% for the last two decades. Compared to the RBA for the Evenkia region, which increased almost 10-fold over the last decade (from 0.2% up to 1.9% per year), the RBA for the total Siberian region was estimated as 0.51-1.19% for different years [56], and the RBA for the forests of western Canada was 0.56% [57]. Thus, a significant "cumulative" effect of the fires provokes the large-scale anomalies in vegetation cover and the thermal balance of the soil surface (Table 3).
It was found that changes in the thermal regime of post-fire plots under the conditions of permafrost are accompanied by an abnormal increase in the average temperature of the surface and upper soil layer in the summer period. It was shown that NDVI values are restored within 7-10 years of the fire (Figure 3b), which is correlated with the rate of post-fire lichen layer recovery (Figure 4b). However, the partial regeneration of vegetation cover did not compensate for the changes, which lead to long-term surface temperature abnormalities. The thermal balance of post-fire sites with disturbed vegetation cover remains affected for more than 15-20 years, which is consistent with similar studies for different test plots [58,59]. A similar effect was observed, in particular, in Alaska [60]. The obtained results are also consistent with the data for Russian forests [10,16]. Thus, surface temperature anomalies are amongst the most significant post-fire effects, which determine the thermal balance of the soils and STLs of disturbed forests under the conditions of the permafrost zone of Central Siberia. Considering the predicted climate changes, the toughening of the fire regime and the increased fire activity in the forests of the permafrost zone [61], the integral effect of post-fire thermal anomalies will likely grow increasingly.

Statistical Data on STL and Regression Model
Numerical simulations and field measurements [3,29] have fixed the average statistical norm of STL in the range of 0.6-2.0 m, for conditions similar to those of the region of interest ( Figure 4a). As shown by the numerical simulation, the excessive heat flux on the surface caused an increase in depth thawing by, on average, 30% (up to 50% in some extreme cases), relative to the average statistical norm, which corresponds to estimates discussed earlier [1]. After a fire, the soil profile thaws for an extra 0.5 m, especially under the conditions of a stable anticyclone, a high level of insolation and high air temperatures, which were observed annually in Central Siberia over 30-50 days in mid-summer [62]. Although the maximum thermal anomaly was observed in mid-summer, the maximum thawing depth occurred in August-September due to the delay in temperature dynamics in the layers of the soil profile [15,16].
The dataset developed on the basis of published data (51 sources, 1062 observations, Table S1 Supplementary Materials) for the Evenkia region demonstrated the high variability of STL values even amongst larch stands of the same age (Figure 4a). Such variation can be caused by differences in soil type, microclimates, vegetation, ecological conditions and landscape heterogeneity [8]. Different disturbance agents (wildfires, permafrost landslides, clear-cuts, etc.) also contribute to the high variability of STL. These heterogeneities restrict model predictions of STL. The obtained dataset from the published research shows the good seasonal relationships between the STL and the soil organic layer (litter) in fire-damaged and non-disturbed larch stands. For the Arctic tundra, the modeling study suggested that STL was strongly linked to moss cover, soil temperature and carbon balance, and may affect the Arctic due to climate change [8].

Numeric Simulation of Soil Temperature and STL Depth
The simulation results depended strongly on the volumetric water/ice content in the loamy horizons of the soil (CR, C). The thermal conductivity function depending on this parameter was different in the models developed by Pavlov [45] and Ekici et al. [51] (Model 1 and Model 2). The thermal conductivity value in Model 2 was lower than in Model 1 for the low (y w < 0.1) values of the volumetric water/ice content in the loamy horizons of the soil. On the contrary, the thermal conductivity was higher in Model 2 when the volumetric water/ice content value was high. This was well reflected in the simulation results for the damaged soil of post-fire plots ( Figure 5). In the case of relatively dry soil (DS, the water volume content was significantly lower, at y w < 0.1), the thawing depth calculated with the use of Model 2 was 30% less than that of Model 1. On the contrary, the thawing depth was significantly shallower in the solution using Model 1 in the case of water-saturated soil (WS), and taking into account the layered distribution of water/ice content. However, the indicated models for the thermal conductivity coefficient of the undamaged section did not significantly affect the thawing depth due to the presence of the upper organic layer (Table 1), which worked as an additional heat-insulating layer [8,12,13].
The results of STL dynamics modeling were compared between post-fire plots and non-disturbed plots in parallel modes of modeling. This was done via the simultaneous observation of comparative changes resulting from the differences in the thermal-isolating moss-lichen cover and the soil organic layer. First, at fire-disturbed plots, the thawing process began earlier (∆T~20 days) (Figure 5b), which did not differ for the types of conditions or the modeling type choice. Second, for fire-disturbed plots, the STL change amplitude (∆Z) varied by 30-50% compared to the non-disturbed plots. The observed range goes well with the published data [3,55] from field observations at the study site (Table S1), and agrees with the evaluations based on the multivariate regression (Figure 4a), taking into account standard deviation and sporadic extremes.
According to results of the numeric simulation for non-disturbed plots (Figure 6a), the temperature of the soil at a depth of 0.05 m exceeded zero, with a one month delay in springtime (T1) compared to the post-fire plots, which is similar to the differences observed for the soils from south-and north-aspected slopes in the area of interest [13]. The same effect was found for the next depths of 0.10 m and 0.20 m (the root layer) (Figure 6b-c). The next peculiarity in the post-fire plots is that in the soil layers of 0.05, 0.10 and 0.20 m, the temperature was higher than the norm for undamaged areas during the late autumn. Long-term inversion was found in the temperature curves (see T 2 , T 4 ). Finally, thawing periods (see T 3 , T 5, Figure 6c-d) were observed in post-fire soils at depths of 0.10 and 0.20 m, while undisturbed soils had negative temperatures, or were close to 0 • C, at the corresponding depths.
In the middle of the growing season, the soil temperature difference at a depth of 0.05, 0.10 and 0.20 m (Figure 6a-c) equaled 5.3-8.0 • C, which was well related to the results from the satellite monitoring of surface temperature for plots in the first year after fire impact ( Table 3). The relative deviation of values for post-fire and non-disturbed plots varied up to 60%. Thus, numerical modeling implemented for the root layer (0.05-0.20 m) confirmed the significant change of the thermal regime [13], which stimulated the radial growth of tree rings (e.g., Larix) [12]. In fire-damaged plots, the absolute values of soil temperatures at a depth of 1.0 m modeled for the mid-growing season differed from the non-disturbed plots by 1.0-1.5 • C; as such, within post-fire plots with anomalies of surface temperature, the thawing period T 5 was predicted even at this depth ( Figure 6d). However, it was not present at the curve of non-disturbed plots. For the lowest soil horizons (e.g., 3 m), the model showed relatively low differences between post-fire and non-damaged cases (Figure 6e). Here, the influences of post-fire changes on the moss layer and soil organic layer (litter) were negligible over the course of the year.
The results of our numeric simulation experiments showed that fire-damaged soils lose their thermo-isolation properties, which directly influences the STL's dynamics. In non-disturbed soils, a lag-period in temperature dynamics (or thermal 'inversion') was observed at the beginning and the end of the growing season ( Figure 6, T 2 , T 4 marks). This may be the reason for the significant post-fire changes in the seasonal dynamics of many processes in the root-inhabited layer (microbial activity, gas exchange, etc.) of permafrost ecosystems [6,12,14,24,29]. Our findings showed that thermal anomalies influenced the thermal regime not just in the growing period (summer), but in the winter season as well. Especially, this was shown for springtime, when we observed higher (and close to 0 • C) soil temperatures at the 0.05-m soil horizon of fire-damaged plots compared to non-disturbed sites (Figure 6a, T 1 mark).

Limitations
The limitations of our work were as follows.
(1) The speed of loss of NDVI deviation, as well as the temperature deviation, should be mediated by fire severity, drainage, stand density, etc. However, in the current study, we did not attempt to quantify such relations. (2) Accounting for the movement of water in the soil is the primary purpose of further developments of the mathematical model. The vertical movement of moisture associated with precipitation, snow melting and the evaporation of moisture from the surface can make a significant contribution to the heat transfer in the soil. Moisture transfer can be significantly different depending on the presence of fire damage, since the damaged and intact organic layer differs in its ability to accumulate moisture. (3) The simulation was carried out under the assumption that plots were damaged by fire recently, in cases of the complete degradation of the moss-lichen cover and disturbance of the forest stand. The complete absence of the upper organic layer (Opir) was assumed in the model for fire-damaged soils. All these remain a topic for future research, and further validations of simulation results are also required.

Conclusions
According to the analysis of satellite imagery for 2000-2019, the cumulative fire disturbance in the territory of interest amounted to 20%. The partial or complete degradation of the litter and the moss-lichen layer after fire impact affect the seasonal temperature of the soil and the depth of the STL under permafrost conditions. In the current study, to simulate STL anomalies immediately after the fire, we used the approximation of the complete degradation of the moss-lichen cover and an assumption of the 100% mortality of the forest stand. This is a highly feasible scenario for the fire-damaged larch forests of the permafrost zone of Siberia.
Changes in the thermal regime of the post-fire areas of the permafrost larch forests of Central Siberia were accompanied by anomalies in the average surface temperature of ∆T = 7.2 ± 1.3 • C. This was up to 50% higher than the temperature norm for non-disturbed sites. While NDVI values were restored 7 to 10 years after the fire, the surface temperature anomalies in larch forests stayed significant over the 15-20 years after burning.
Seasonal relationships between STL and stand age in fire-damaged and non-disturbed larch stand, on the basis of published research in the permafrost area (Tura, Evenkia region), were found.
The mathematical model applied for two extreme cases of soil conditions (dry soil and water-saturated soil), and for soil with non-homogeneous water content, determined that STL thickness could undergo a 30-50% increase compared to the statistical norm for differently aged permafrost larch forests in Central Siberia. The results of our numerical simulation experiments showed that fire-damaged soils lose their thermo-isolation properties, which directly affects the STL dynamics. In the soils of non-disturbed plots, a lag-period in temperature dynamic (or thermal 'inversion') was observed at the beginning and the end of the growing season. Our findings presented the all-year-round impact of thermal anomalies on the thermal regime, including the winter period as well. Thus, at spring, fire-damaged soils were warmer, and close to 0 • C at the 0.05-m soil horizon, compared to non-disturbed soils.
The slow rate of the loss of surface temperature anomalies led us to consider this factor as one of the most significant for the state of the permafrost's soils, and the STL in particular. This appeared to be a significant factor in the stability of Siberian permafrost ecosystems, requiring long-term monitoring. Considering the predicted climate change, the worsening of the fire regime and the increased fire activity in the forests of the permafrost zone, the integral effects of post-fire disturbances will likely grow increasingly. Further observation is needed in order to assess how these processes contribute to the global situation under ongoing climate change.