Soil Temperature in Disturbed Ecosystems of Central Siberia: Remote Sensing Data and Numerical Simulation

Abstract

We investigated changes in the temperature regime of post-fire and post-technogenic cryogenic soils of Central Siberia using remote sensing data and results of numerical simulation. We have selected the time series of satellite data for two variants of plots with disturbed vegetation and on-ground cover: natural ecosystems of post-fire plots and post-technogenic plots with reclamation as well as dumps without reclamation. Surface thermal anomalies and temperature in soil horizons were evaluated from remote data and numerical simulation and compared with summarized experimental data. We estimated the influence of soil profile disturbances on the temperature anomalies forming on the surface and in soil horizons based on the results of heat transfer modeling in the soil profile. According to remote sensing data, within 20 years, the thermal insulation properties of the vegetation cover restore in the post-fire areas, and the relative temperature anomaly reaches the level of background values. In post-technogenic plots, conditions are more “contrast” comparing to the background, and the process of the thermal regime restoration takes a longer time (>60 years). Forming “neo-technogenic ecosystems” are distinct in special thermal regimes of soils that differ from the background ones both in reclamated and in non-reclamated plots. An assumption was made of the changes in the moisture content regime as the main factor causing the long-term existence of thermal anomalies in the upper soil horizons of disturbed plots. In addition, we discussed the formation of transition zones (“ecotones”) along the periphery of the disturbed plots due to horizontal heat transfer.

1. Introduction

Regimes of soil functioning, dynamics, and rate of soil processes are subjects of physical characteristics of soils (temperature and moisture content). Analysis and modeling of temperature distribution in horizons of cryolithozone soils, including the seasonally thawed layer (STL) of soil, are urgent and widely discussed problems [1,2,3,4,5,6,7]. Central Siberia is an important subject of these studies due to that 50% of forest ecosystems (about 3 mln km² of 6 mln km² of forested area) is in the permafrost zone [8].

Changes in thermal regimes of soil reveal themselves in the root-inhabited soil layer (active layer of soil). Possible changes in the state of permafrost soils are a significant factor affecting eco-systems [7,9,10,11,12,13] and even on global climate [14,15,16].

The heat exchange in the soil profile is a subject of the temperature gradient, thermal radiation balance, turbulent heat exchange in air, and evaporation and soil moisture. Temperature gradients determine the heterogeneity of the soil cover functioning [17]. The vegetation cover causes the most significant influence on the soil thermal regime, moss–lichen layer, and organic soil layer [18,19,20,21]. The temperature of the upper soil layer (0–0.4 m of depth) depends on the thermal insulating properties of moss and lichen covers and on the thickness of the organic soil horizon [21,22,23], as well as on the tree stand crown canopy [24]. Heat exchange in soils depends also on the granulometric composition of mineral horizons [18].

Under conditions of loss of thermal insulating organic layer, soil heating causes changes in ecosystems functioning according to the seasonal variations of the active layer of soil comparing to the statistical norm [3,9,13]. Such a situation may arise as a result of a number of destructive natural or human-made factors (wildfires, logging, and activities of the mining complex).

Disturbances may have a larger and immediate impact on Low Arctic ecosystems of Siberia than temperature warming alone. While soil conditions strongly affect vegetation patterns, plants may also have strong feedback effects on soil thermal and hydrological properties [25].

Destructive impacts affect surface albedo, emissivity, moisture, and water regimes of the upper soil horizons. As a result, there is a significant change in the temperature regime of soils comparing to the background. It is estimated that the post-fire effects have an “accumulative” character [26,27] on the area up to 20% of boreal forests of Siberian and can stay significant for 15 years [11,27]. This factor is likely to increase in the future [28]. Anomalous heating of the surface is also observed for post-technogenic areas (open pit mining, quarries, overburden dumps, logging, etc.), characterized by the intense mechanical impact on the vegetation cover and soil [29,30]. This issue is relevant for the zone of the resource-mining complex of Siberia and, in particular, for the Arctic zone (>65° N).

In the context of the diversity of soil properties, the issue of mathematical modeling based on ground-based and remote sensing data is urgent [7,27,31,32,33,34,35]. Available models of heat and moisture transfer in soil under freezing and thawing conditions take into account the moisture conductivity of soil [36,37,38,39]. An obligatory element of models of moisture transfer in freezing soil is taking into account cryosuction, i.e., capillary force associated with the presence of a phase transition [40,41,42], natural convection in water that fills soil pores [43], and vapor transfer in water-free soil pores and snow cover [44,45,46]. After all, it is important to account correctly for solar radiation, which depends on the location and state of the atmosphere [47,48].

The main aim of this work was to assess the influence of the transformation of the soil profile on the formation of temperature anomalies in the surface and in the soil horizons in disturbed plots based on numerical simulation of heat transfer in the soil profile. The following aspects of the issue were considered:

• (1). remote sensing of the state of disturbed natural and technogenic ecosystems based on satellite monitoring data and subsequent numerical assessment of the influence of changes in the optical properties of the surface (spectral albedo and surface emissivity) and variability of soil moisture content on the formation of thermal anomalies in disturbed natural and technogenic ecosystems;

• (2). modeling of conditions for abnormal heating of soil horizons in disturbed areas during the growing and non-growing seasons of a year;

• (3). time-lag of recovery of the thermal state in disturbed ecosystems and the level of stabilization of the temperature characteristics of the surface under conditions of long-term post-pyrogenic and post-technogenic recovery;

• (4). analysis of heat transfer at the boundary of disturbed plots/background areas and the peculiarities of the temperature regime of soils in the formed transition zones (“ecotones”) along the periphery of the disturbed plots due to horizontal heat transfer.

In the current work, it will be demonstrated that the characteristic features of the temperature effect of the damaging impact, revealed during field observations, are reproducible in a mathematical model that takes into account the differences in the structure and moisture content of disturbed and undisturbed soil. We intend to verify, by means of the numerical simulation, that not only the change in the optical properties of the surface, but also the disturbance of the structure and properties of soil horizons are the factor providing the long-term existence of temperature anomalies.

2. Area of Interest and Research Objects

In the current research, we took into account the properties of cryogenic soils, which are typical for the taiga region of Central Siberia (59–66° N, 90–107° E). The forest stands of the region are dominated by larch (Larix sibirica Lebed., Larix gmelinii Rupr. Rupr.) and sparse larch stands, accounting for >50% of the total forested area [26,49,50]. Herb-shrub layer, mosses, and lichens represent the vegetation cover as well as litter loads [51].

The continental climate of the permafrost zone and characteristics of parent rocks determine the genetic specificity of the soils in the region. The soils are characterized by low thickness (0.20–0.50 m) and weak differentiation of the mineral part of the soil profile. The permafrost level (<1.0 m) depends significantly on the topography [11]. Permafrost is widespread on flat areas and is not found in the soil profile on drained areas of slopes and valleys. Periodic freezing and thawing of the soils determine the soil structure, the low rate of humus accumulation in the soil, the peculiarities of moisture transport, and waterlogging. In addition, freezing and thawing processes affect slightly acidic or neutral reaction of the soil, as well as loamy or clayey granulometric composition with a high gravel content in the lower part of the soil profile [49,50]. Finally, cryoturbation, solifluction, cryogenic heaving, thermokarst, etc. characterize such soils.

According to the Russian soil classification of 2004 [52], the soil cover is represented by cryozems (In World Reference Base for Soil Resources (WRB, 2014 [53]) equivalents is Histic Cryosols), podburs (WRB equivalents is Cambic Crysols), and peat soils (WRB equivalents is Cryic Histosols). Cryozems are formed in the lower parts of slopes and microdepressions on rocks of heavy granulometric composition. The soil profile contains a peat litter (thickness of 0.05–0.10 m), a thin humus or coarse humus horizon, and a uniformly mixed cryoturbated mineral horizon with inclusions of plant residues of varying degrees of decomposition. A high (10–14%) humus content characterizes such soils, as the consequence of the low rate of organic matter decomposition. The lower supra-permafrost part of the soil profile is saturated with moisture, structureless, and is more dense and homogeneous. Podburs are formed on fine-earth-clastic products of destruction of metamorphic rocks, which leads to good drainage and to large thickness of the active soil layer. The profile consists of a peat litter, under which lies the Al-Fe-humus horizon, which gradually turns into the parent rock. Peat-cryozems are present in combination with cryozems, located in low relief areas, and are diagnosed by the presence of peat (T) and cryoturbated (CR) horizons.

Background cryogenic soils typically are full-profile, and according to the modern classification [50] have the formula T(Oao)-O(F+H)-CR-C_┴ (Figure 1).

Wildfire impact changes the structural organization of the soil profile. After a wildfire the soil loses its upper horizon, the level of permafrost occurrence decreases, and the soil formula transforms to Opir-O(F+H)pir-CR-C. Under the conditions of high-intensity burning, the profile structure can transform more strongly. Depending on the degree of litter burnout the soil formula changes to O(F+H)pir-CRpir-C or O(H)pir-CRpir-C. Locally, complete burnout of the upper organic horizons of the soil is possible [54]. Upper part of the soil profile, as a rule, have a thin humus-accumulative (underdeveloped) horizon W or a fragmented coarse humus horizon Oao both in case of post-fire mineralization of cryozems and in case of technogenic transforming of soils. Such soils formula is Oao(W)-C.

Evolution of soils after destructive impacts on thermal insulating covers includes an increase in thawing depth, an increase in soil moisture due to thawing of schlieren and wedge ice, an increase in gleying, and cryoturbation processes, and it proceeds under strong influence of succession.

For instance, the effect of post-fire rapid increase in the thawing depth (by 30–100% in the first year after the fire) [34] leads to a unique path of soil evolution, when the underlying gleyed rock layers are included in soil formation for a significant period [55]. This type of evolution is possible only in the permafrost zone. Moreover, an increase in the STL thickness leads to the thawing of schlieren (sometimes wedge ice). An excessive amount of moisture under conditions of difficult drainage stimulates the gleying process [55].

The objects of modeling were (i) native variant, which is undisturbed natural state of cryozem, and (ii) variants of soils with the topsoil layer disturbed as a result of a destructive factor.

The set of input parameters (Table 1 and Table 2) for the numerical simulation of the soil temperature regime was summarized from the data of field studies of the morphological and physical properties of soils carried out in the region [49,50,56,57,58,59].

According to field experiments in the permafrost zone [50], “thermal inversion” is possible in different layers of cryogenic soils of disturbed plots. Thus, there is stronger heating of the surface in summer (due to lower thermal resistance). Next, there is a higher temperature in the soil horizons during the non-growing period (early spring, winter, autumn) (due to the higher thermal resistance) in comparison with the distribution in the background plots. During the growing season (summer), such differences are explained by the state (or loss) of disturbed upper soil horizons and the reduced albedo value. The relative change in thermal resistance during the non-growing season can only be explained by a higher ice content in the soils of non-damaged plots, if the snow cover is equal.

As the model used did not consider the transfer of moisture within the soil, values of the volumetric water/ice content were tabulated as constants for the summer and winter seasons (Table 1 and Table 2).

3. Materials and Methods

3.1. Multispectral and Infrared Survey of Research Objects

We used satellite images of average spatial resolution (15–30 m) from Landsat-5/7/8 for 1975–2020, which are freely available in The United States Geological Survey (USGS) database (https://earthexplorer.usgs.gov/, accessed on 21 May 2021). The surface temperature was evaluated from the calibrated В6 channel (λ = 10.4–12.5 μm, Landsat-5/TM—Thematic Mapper), B6/1 channel (Landsat-7/ETM—Enhanced Thematic Mapper), and B10 channel (λ = 10.6–11.9 μm, Landsat-8/OLI—Operational Land Imager). We implemented radiometric correction method to the initial data using calibration constants from the metadata files [60,61]. First, we analyzed the availability of the imagery for selected experimental post-technogenic and post-fire plots to obtain time-series for each plot for 1975–2020. The only criteria were the percentage of cloudiness, the binding of the survey dates to the middle of the growing season (July), and the regularity of survey covering of selected objects. Daytime imagery was used only. Thus, we did not implement any special procedure for the selection of Landsat images.

Additionally, we used low spatial resolution survey data (250–1000 m) from Terra/MODIS (Moderate Resolution Imaging Spectroradiometer) for 2002–2020. Standard products L2G and L3 used (free USGS database, https://lpdaac.usgs.gov/dataset_discovery/modis, accessed on 21 May 2021). We operated with reflectances (albedo) measured in MODIS band #1 (λ = 0.620–0.670 μm) and band #2 (λ = 0.841–0.876 μm) (product MOD09GQ), and to analyze surface temperature anomalies we used the MOD11A1 standard product. We used 10-day averaged data across the entire set of initial data, with a special focus on the disturbed ecosystem’s recovery succession stages. A monthly data average procedure, as well as a procedure of spatial averaging within each natural and technogenic disturbed ecosystem plots, were applied.

We have selected the time series of satellite data for two variants of plots with disturbed vegetation and on-ground cover: (1) natural ecosystems of post-fire plots (PF) for 1996–2018 and (2) post-technogenic plots for 1975–2019.

Using time series of satellite data, we performed indirect assessments of the state of disturbed ecosystems (disturbed plots) at different periods of recovery successions based on changes in the anomalies of spectral features comparing to characteristics of the background areas. We analyzed data for 30 post-fire plots recorded in Siberia in 1996, 2006, and 2017 (Figure 2a), and the dynamics of spectral anomalies for two technogenic plots: Borodinsky Coal Mine (BCM) for 1975–2018 and Olimpiada Mining Plant (OMP) for 2000, 2011, 2019 (Figure 2b,c).

We analyzed the dynamics of albedo and temperature anomalies during 1, 5, 10 (or 12), and 20 years of post-fire vegetation and on-ground cover restoration. For technogenic disturbances, we analyzed albedo and surface temperatures during 1, 10, 20, 40, and 60 years after impact a destructive factor on soil covers.

To determine the relative anomalies of the surface temperature of disturbed plots to the characteristics of background territory (ΔT/T_bg (°C/°C), in %), we compared averaging over 10 measurements for each period of recovery:

ΔT/T_bg = 100%·(T_tg—T_bg)/T_bg,(1)

where T_tg is surface temperature of target (disturbed plot), °C, and T_bg is surface temperature of background (non-disturbed) plot, °C.

3.2. Mathematical Model of Heat Transfer in the Soil

Mathematical models of heat transfer in soils taking into account water-ice phase transformation are quite diverse [62]. The choice of the model depends on the specifics of the problem being solved. Climatological studies often impose severe constraints on the computational resources that can be assigned to calculations of soil heat transfer, therefore analytical and semi-analytical methods are widely used [63,64]. These methods are based on constant or piecewise constant approximations of the main thermophysical characteristics (density, specific heat capacity, and thermal conductivity coefficients) in their dependence on the depth and on the assumption that the temperature depends on the vertical coordinate only.

The formulation of the issue of heat transfer in disturbed cryogenic soil took into account the following factors.

• (1). The multidimensional character of the temperature distribution, being a result of vicinity of disturbed and non-disturbed plots and, accordingly, the presence of a significant difference in the state of the permafrost layer. Under these conditions, significant heat transfer in the horizontal direction can arise.

• (2). Inhomogeneity of the thermophysical properties of the soil, due to the vertical structure of the soil profile and its differences in disturbed and background plots.

• (3). The presence of a transition zone with a mixed phase state of water and ice. Such areas inevitably arise due to local heterogeneities of the soil, its composition, and moisture content. For the used scale of grid sampling in the considered problem, the transition zone is described in terms of the average content of the liquid and solid phases of water.

In a mathematical model of heat transfer in soil, either the temperature field or the specific enthalpy field can be considered as a sought field [65]. The choice of temperature as a sought field is more expedient for the considered problem. This facilitates greatly the process of composition of the thermophysical properties of soil. Thus, the following formulation for the heat transport equation was chosen:

(2)$c_p\left(x,T\right)\rho \left(x\right)\frac{\partial T}{\partial t}-\frac{\partial }{\partial x_i}\left(\lambda \left(x,T\right)\frac{\partial T}{\partial x_i}\right)=S,$

where T is the temperature (K), $c_p\left(T\right)$ is the thermal capacity (J/(kg∙K)), ρ is the density (kg/m³), $\lambda \left(T\right)$ is the thermal conductivity coefficient (W/(m·K)), and S is the heat source associated with the phase change latent heat (W/m³).

Formulating the model of the heat source S, we consider the liquid and solid phases as separate components of the medium. They can present at the same place (i.e., in the same finite volume of the discrete mathematical model) simultaneously. This approach allows taking into account a fine structure of intermitting water and ice that is typical for freezing and thawing soils [66] but cannot be resolved in a discrete model due to the restrictions for a grid dimension. The mass sources of solid and liquid phases are opposite and proportional to the temperature on the Celsius scale [66]. In accordance with the chosen approach, the source term S in Equation (2) has the following form:

(3)$S=-y_w\left(x,T\right){\rho }_{w}L\frac{\partial {\alpha }_l}{\partial t},$

where $y_w$is the soil volumetric content of water regardless to the state (m³/m³), ${\alpha }_l$ is the liquid water fraction, and L is the latent heat (J/kg).

The phase transition rate is determined by the heat transfer flux to the contacts between the phases. This flux is proportional to the difference between the mean temperature of the soil and the phase change temperature, in other words, to the mean temperature of the soil measured in Celsius degrees. Considering the phase transition process as the exponential decay of the diminishing phase, we can write the following equation for the liquid water fraction:

(4)$\frac{\partial {\alpha }_l}{\partial t}=k(\max \left(T \left[°C\right], 0\right)\left(1-{\alpha }_l\right)+\min \left(T\left[°C\right], 0\right){\alpha }_l),$

where $k$ is a coefficient connecting the characteristic timescale of the exponential decay of the diminishing phase and temperature.

Equation (2) was solved in a two-dimensional formulation, in which the vector of coordinate x has vertical and horizontal components (to take into account the effects of the transition region between disturbed and undisturbed soil areas). Large-scale variability of the soil structure and its moisture content, as well as the presence of dynamically changing depth of snow cover, are taken into account in the expressions for the coefficients of the Equation (2): $c_p\left(x,T\right), \rho \left(x\right), \lambda \left(x,T\right), y_w\left(x,T\right)$.

Equation (2) was solved with the use of the finite volume method in that the equation domain is divided into a large number of cells that have four faces as the case is two-dimensional. The sought temperature fields are determined as a set of temperature values at the centers of the cells and the differential Equation (2) is replaced by a system of linear algebraic equations (grid equations) connecting the values at the centers of the cells [67]. Integration of Equation (2) along the cell volumes leads to the algebraic grid equations in that the sought variables are the temperature values in the cell centers [68].

Available models [1,4,5,69,70] were used to derive expressions for the thermal conductivity depending on the vertical coordinate. Empirical models by Pavlov [69], Anisimov et al. [70] are results of the generalization of the great amount of field measurements. In these models, the thermal conductivity coefficient is a function of the dry soil density and the volume fraction of water or ice. The models proposed in the works [1,4,5] are based on the supposition [71] that the thermal conductivity of a partially water-saturated soil correlates linearly with the thermal conductivities of a dry and a water-saturated soil:

(5)$\lambda =Ke{\lambda }_{sat}+\left(1-Ke\right){\lambda }_d,$

where ${\lambda }_{sat}$ and ${\lambda }_d$ are the thermal conductivity coefficients of a soil in the water-saturated and dry states (W/(m·K)), and Ke is Kersten number.

We used the model by Porada et al. [5] for the moss and lichen layer (Oao horizon). For the horizon О(F+H) that consists mainly of medium and well decomposed organic residues, the model by Anisimov et al. was used [70]. For the following two mineral horizons CR and C two types of models were considered [4,69] (Table 3). In calculations relating to periods with stable snow cover the thermal conductivity of snow was admitted to be a constant, ${\lambda }_{snow}$ = 0.23 W/(m·K).

Other thermophysical properties of snow (density, heat capacity, and albedo) were also assumed to be constant, but the depth of the snow cover was considered as alternating according to meteorological data.

The boundary conditions at the upper boundary of a soil or snow take into account the presence of solar (shortwave) and atmospheric (longwave) radiation and convective heat transfer by air [72,73]:

(6)$-\lambda \left(\frac{\partial T}{\partial x}\right)=q_c+\left(1-{\alpha }_s\right)·q_{Sr}+q_{Lr},$

where $q_c$ is the convective heat flux, $q_{Sr}$ and $q_{Lr}$ are the shortwave and longwave radiation heat fluxes, and α_s is the surface albedo.

The surface convective heat flux is calculated with the use of the following expression [69]:

(7)$\begin{array}{c}{q}_c=a_s·\left(T_{atm}-T_s\right),\\ {\alpha }_s=u_1^{0.5}·\left(6+3.1·\Delta T/u_1^2\right)\end{array}$

where ${\alpha }_s$ is the heat transfer coefficient, W/(K∙m²); $u_1$ is the wind velocity at the altitude 1 m; and $\Delta T$ is the temperature difference between the solid surface and atmosphere.

The shortwave radiation heat flux divides into the direct solar radiation $\left(I_{Si}^{Dir}\right)$ and the scattered one $\left(I_{Si}^{Dif}\right)$:

(8)$q_{Sr}=q_{Sr}^{Dir}+q_{Sr}^{Dif},$

The direct solar heat flux on a horizontal surface is found as [47]

(9)$q_{Sr}^{Dir}=q_{Sr0}·\cos Z·{{\displaystyle \prod }}_i{\tau }_i,$

where Z is the reduced incidence angle, ${\tau }_i$are the atmosphere transmittance factors that characterize absorption by aerosols, water vapor, ozone, other gases, and Rayleigh scattering; $q_{Sr0}$ (W/m²), is the Solar heat flux reaching the Earth atmosphere. The value of $q_{Sr0}$ depends on the position of the Earth relative to the Sun [73]. The factors ${\tau }_i$ in (9) are calculated according to [47].

Contributions to the scattered radiation flux are made by the Rayleigh and aerosol scattering and multiple reflection of the direct solar radiation. The meteorological data for cloud cover were used to calculate the change in relative insolation.

Longwave radiation heat exchange between the solid surface and the atmosphere is calculated as

(10)$q_{Lr}={\epsilon }_s·\left({\epsilon }_{atm}\sigma ·T_{atm}^4-\sigma ·T_s^4\right),$

where${\epsilon }_s$ is the surface emissivity and ${\epsilon }_{atm}$ is the atmosphere emissivity. For the surface of the damaged and undamaged soil, the values of the integral albedo were taken equal to 0.13 and 0.18 (a decrease by 40% for the case of a significant change in the state of the ground cover and the upper organic horizon), the emissivity values are 0.9 and 0.98, respectively [74]. The atmosphere emissivity ${\epsilon }_{atm}$ was calculated as a function of cloud cover and humidity [75,76].

3.3. Meteorological Data

We used long-term series of meteorological information from open data banks for the period 1990–2020: Climatic Research Unit (http://www.cru.uea.ac.uk/, accessed on 25 March 2021), Weather archive (http://rp5.ru, accessed on 25 March 2021), and the National Climatic Data Center (NCDC) (http://www7.ncdc.noaa.gov/CDO/cdo, accessed on 25 March 2021).

To calculate solar heat flux (4), we used the geographic coordinates of Tura (64° N, 100° E) (Evenkia, Krasnoyarsk region, Siberia, Russia). The data on the depth of snow cover (mm), air temperature (°C), wind speed (m/s), and cloudiness (oktas) were taken from the annual series of meteorological data. The observations have a regularity of 4 h (6 times a day). Numerical simulations are performed for 2018 data (Table S1, Supplementary Materials).

3.4. Main Features and Assumptions of the Mathematical Model of Heat Transfer in the Soil

The governing equations of the mathematical model are (1) and (4), that describe heat transfer in soil in presence of water-ice phase transition neglecting the moisture transfer processes. The main thermophysical properties, in addition to thermal conductivity, were set constant for different soil horizons, but different for temperatures above and below 0 °C. The calculations were performed for a time interval up to three years required to establish a quasiperiodic solution.

The result of the numerical solution of the model equations was the dynamics of two-dimensional temperature fields and phase transition boundaries.

In the mathematical model, the discretization of the computational domain was built based on a Cartesian (orthogonal, structured) grid. The vertical grid step was 0.01 m, the horizontal one was from 0.25 m, and the time step was 15 min.

In the mathematical model, the boundaries of the soil horizons were assumed constant and dynamic change of the snow cover depth according to the data of meteorological observations was taken into account. A linear temperature distribution changing from the air temperature at the surface to the specified temperature at the lower boundary (Z_end) (see Table 1 and Table 2) of the equations domain was assumed as the initial temperature distribution in the soil.

4. Results

4.1. Surface Temperature Anomalies of Disturbed Plots on Satellite Data in IR Range

Time series of satellite imagery reflect the change in spectral characteristics and temperature anomalies of disturbed areas caused by recovery of vegetation and of a ground cover (Figure 3). For PF plots it was evaluated for the period of 20 years, for both BCM and OMP plots it was analyzed for the period of >60 years. The peculiarity of plots after fire disturbances is that natural recovery processes occur there. The peculiarity of post-technogenic plots is that there are two possible options for dynamics of the ecosystem: (1) restoration after reclamation or (2) long-term condition in the format of non-reclaimed lands (dumps, mineralized surfaces). Each of these options is of interest from the point of view of the formation of surface temperature anomalies and its effect on the properties of soil horizons.

The long-term dynamics of spectral albedo (α) and surface temperature anomalies in disturbed plots of natural and technogenic landscapes in comparison with the background data are summarized in Table 4.

4.2. The Results of Numerical Simulation

The direct results of numerical simulation are dynamics of fields of temperature and fraction of liquid and solid states of water, radiation, and convective heat flux on the surface (Table S2, Supplementary Materials). The phase transition boundary and the magnitude of the relative surface temperature anomaly for disturbed areas in comparison with the background ones are the indicators most important for analysis.

Figure 4a shows the formation of a periodic solution during a three-year calculation period. Figure 4b demonstrates the differences in the dynamics of the phase transition boundary for damaged and undamaged soils. Figure 4c shows the dynamics of the surface temperature anomaly in the summer months through the values of the relative temperature deviation (ΔT/T_bg, °C/°C) from background values. Figure 4d reveals the influence of the season and optical properties of the surface (soil or snow) on the annual dynamics of the total radiation flux.

The influence of changes in the optical properties and moisture content of disturbed soil areas on the formation of thermal anomalies in the numerical simulation results are shown in Table 5. Several variants of the set of soil conditions were examined. One of them was considered as the basic set of conditions (Table 1 and Table 2) and others differed from it in a certain parameter value in order to reveal the influence of this parameter. We considered the following variants: (i) basic set of conditions, (ii) no change in the moisture content for the horizon O(F+H), (iii) equal moisture content in the upper soil horizons, and (iv) moisture content similar to var. 3 with greater change in albedo. In the cases considered, ΔT/T_bg varied from 5 to 25%.

5. Discussion

5.1. Long-Term Surface Temperature Anomalies

Temperature anomalies in the post-fire areas are typical for permafrost conditions of Siberia [11,12,56].

A decrease in the spectral and broadband albedo (α, %) of the surface in post-fire areas, due to partial or complete combustion of the ground cover, provokes excessive heating of the surface. As shown by the remote measurements (Table 4), the spectral albedo in the short-wavelength bands (MODIS band #1 and #2) is reduced by 20–48% relative to the background immediately after the fire, by 15–25% after 10 years of recovery, and by 3–12% after 20 years of vegetation regeneration. For the variant of post-technogenic plot BCM (with reclamation), the albedo values correspond to the dynamics in PF, while for OMP (dumps without reclamation) the level of the albedo reduction remained 45–60% throughout the observation period (Table 4).

Under conditions of positive air temperatures in summer, as well as in spring and early-autumn, the effect of significant surface temperature anomalies (ΔT/T_bg, %) was recorded in all variants of the disturbed plots (PF, BCM, OMP) in comparison with background surface temperature.

In PF plot ΔT/T_bg on average reaches values of 33 ± 6% with maxima of ~46%, which were recorded immediately after fire impact. During 10–12 years, on average, the value of the temperature anomaly decreases to ~12.5 ± 1% (sporadic maxima ~20%). During at least 20 years (the stabilization time lag T_stab, Figure 5a), the anomaly amplitudes decrease to the level of background values and the measurement error of 3 ± 1% (Table 4, Figure 5a). Thus, during the period under consideration, stabilization of temperature parameters is observed in PF plot, which is determined by the thermal insulation properties of the surface layer recovering to the pre-fire state under conditions of a successful vegetation regeneration [11,34,56,77].

The considered post-technogenic plots were characterized by a higher level of initial relative surface temperature anomalies (up to 100% of the background). Further, an exponential decrease was observed, similar to PF plot (R² = 0.87–0.92) (Figure 5). However, the stabilization time lag (T_stab) is 40–60 years, which is characterized with ΔT/T_bg of 18–20% under the conditions of recovery (Figure 5b).

A significantly high level of ΔT/T_bg in technogenic areas remains twice as long as in PF plot. In disturbed ecosystems, the morphological and physical properties of the upper soil horizons are not restored to the background level for a long time [56,78]. Thermal regimes remain significantly anomalous for up to 60–80 years (Figure 5b). Even after the time of stabilization (T_stab > 60 years), the level of relative anomaly overestimated at least by 15–20% in relation to the background values. Probably, we can talk about “neo-technogenic ecosystems” forming, which are characterized by special thermal regimes of soils that differ from the background ones both for reclamated (BCM) and for non-reclamated (OMP) plots. The main reason for such changes is a significant mechanical effect on the soil and/or complete destruction of its structure, which restoration requires much longer time.

These results can be explained by the fact that mechanical reclamation cannot provide an increase in the rate of restoration of the state of the litter and upper soil horizons to the level of the background territories. Although, the restoration of vegetation (as a factor that changes the albedo) proceeds faster in reclaimed areas, where successions acquire a character closer to natural restoration processes.

The formation of thermal anomalies in disturbed plots can be caused by a sharp change in three factors (i) surface emissivity, (ii) albedo, and (iii) moisture content in the upper soil horizons.

Surface emissivity determines the heat loss due to the re-emission of absorbed short-wave solar radiation into the atmosphere. The coefficient of surface emissivity decreases in the disturbed plots [74,77]. However, heat losses from the surface with damaged vegetation and soil cover at significant ΔT/T_bg should be higher than from the undisturbed surface, due to the strong nonlinear dependence of thermal radiation on temperature (according to Stefan–Boltzmann law). Consequently, heat losses contribute to a decrease in ΔT/T_bg values. Thus, the surface emissivity does not fully determine the formation of thermal anomalies.

According to field observations conducted in the study area [7,50], changes in albedo and moisture content conditions can also contribute to the formation of thermal anomalies.

The spectral albedo of disturbed areas in all variants (PF, BCM, OMP) significantly differs from the background values (Table 4). At the initial stages (during the first years after the destructive factor impact), this factor significantly affects the absorption of solar radiation. However, it is numerically shown (Table 5) that the albedo reduction by 45–55% leads to a decrease in the proportion of absorbed solar radiation by <10%, and only by 6% for the model of “native” variant (see Table 5).

In addition to changing the optical properties in disturbed areas, conditions are formed for changing the water content in the near-surface soil horizons [35]. Such changes in moisture content directly depend on the degree of disturbance of the vegetation cover and soil, including the disturbance of its structural organization. Probably, it is this parameter that determines thermal anomalies and their presence at later stages of recovery processes (10–20 years in PF plots or >60 years in BCM and OMP), when the spectral albedo and the surface emissivity are leveled as a result of the restoration of vegetation and litter cover. It was previously shown [27,79] that vegetation signs are restored no later than 7 years after a wildfire. Nevertheless, a change in the soil moisture regime in disturbed plots provokes a decrease in the thermal conductivity of the upper soil horizons. We believe that this is the main reason for the existence of significant temperature anomalies even after the normalization of the optical properties (albedo, surface emissivity) of the disturbed ecosystems.

5.2. Results of Numerical Simulation of the Annual Dynamics of Heat Transfer

As the results of the numerical simulation show, the initial temperature distribution is changed by a periodically time-dependent one starting from the second year. At a depth of 0.1 m, this change occurs after five months (Figure 4a).

At the beginning of the negative air temperatures period a thawed cavity is formed, surrounded by frozen soil in the area above the permafrost zone (Figure 4b). It decreases on both upper and lower sides. This phenomenon exists for 14 days in the non-disturbed soil, and up to 21 days in the soil of the damaged plot. Freezing in the intact soil starts two weeks later than in the soil of the damaged plot.

The maximum inflow and losses of heat due to longwave and shortwave radiation occur in summer and amount to about 350 and 100 W/m², respectively (Figure 4d). Heat loss in summer occurs mostly during nights. In winter, with a high snow albedo, the soil loses 30–50 W/m² due to radiation (Figure 4d). Abrupt changes in the profile of the radiation flux are associated with snowfalls and melting of snow, which affects the surface albedo sharply.

The numerical solution shows a positive temperature anomaly ΔT/T_bg = 15–25% in the damaged areas in summer (Table 5, var. 1), which correlates with the data of satellite observations (Table 4, Figure 5), in particular for the post-fire areas (PF). Such values of ΔT/T_bg are achieved in the case of a simultaneous significant decrease in moisture content and a change in the optical characteristics of damaged areas.

The simulation results demonstrate that a significant role in the formation of the temperature anomaly belongs to the change in moisture content. When the influence of the moisture content factor is excluded, the effect of changing in the radiant flux does not exceed ΔT/T_bg = 5% (Table 5, var. 3) that is 4–5 times less than data from remote sensing. With a decrease in moisture content for the upper soil horizon to the level of the second soil horizon of the non-disturbed soil, ΔT/T_bg becomes two times higher but remains below the observed values (Table 5, var. 2). Only with further 2.5 times decrease in the moisture content (Table 5, var. 1), the simulation results correlate with the data of satellite observations. To reproduce the observed thermal anomaly in the mathematical model changing only the amount of the absorbed radiation, the albedo in the disturbed area must be several times less than in the undisturbed areas (Table 5, var. 4).

5.3. Thermal Inversion

Thermal inversion in soil temperature values in comparison of disturbed and undisturbed areas was observed both in numerical simulation (Figure 6a) and in field measurements (Figure 6b) [50]. The results of numerical simulation allow us to state that the condition for this effect is a higher ice content in the soils of the background areas, in comparison with the variants of disturbed territories.

Due to the higher ice content, despite the presence of the heat-insulating soil horizon Oao, the temperature in the soil of non-disturbed plots decreases faster (Table 3 [5,70]). In the summer, the opposite situation is realized. Greater heating of soils in disturbed plots causes an increase in the depth of the active layer by ΔZ = 29% according to the results of numerical modeling (Figure 4b) in comparison with the background areas, both due to the absence of an additional thermal insulating horizon and due to lower albedo and emissivity values. The calculated increase in thawing depth for disturbed soil is consistent with observational data, which show the magnitude of the relative increase in thawing depth in the range of 25–40% [34].

Comparison of the results of modeling and field observations (Figure 6a,b) shows the similarity of heat transfer processes: thermal inversion, similar nature of temperature dynamics at a depth of 0.1 m for soils of disturbed and background areas. The most significant differences are observed in the summer period (in general, the maxima are higher by about 5 °C), which can be explained by differences in meteorological conditions and model assumptions of soil characteristics (Table 1 and Table 2) that differ from natural ones. In particular, in the presented results for the mineral horizon CR clayey granulometric composition, the model from [69] was used, as the model from the works in [1,4] led to a significant overestimation of the depth of soil thawing [34]. The necessity to refine the model coefficients in accordance with the characteristics of natural soils is obvious.

5.4. Transition Zone

The problem of heat transfer simulation is two-dimensional due to boundary effects at the periphery of the damaged plots. The calculated two-dimensional distribution of the temperature field in the summer and autumn months (Figure 7) indicates that at a depth of up to 0.5 m (within the root-inhabited layer), the horizontal propagation of thermal disturbances from plots with damaged soil cover to a non-damaged area can reach ~5.5 m.

Transferring the modeling results to post-fire and post-technogenic disturbed ecosystems, it should be taken into account that along the boundary of disturbed plots, transitional zones (or “ecotones”) (in terms of thermal and moisture conditions) are formed. Their characteristic width is up to ~5 m and area can be up to 6% of the area of the disturbed plots (Figure 8).

The zones formed along the border of disturbed areas should also be considered as transitional zones or “ecotones” characterized by “edge effect”, where conditions favorable for vegetation and biota can develop in the root-inhabited layer of the soil. This factor can be favorable for conditions of vegetation and biota earlier exit from the winter state due to increased temperatures in the root-inhabited layer, which is confirmed for Siberia by field studies [9,10,13,56]. As an ecological factor, this can determine an increase in the density of some populations and be a zone of additional introduction of plant, microbial, and animal communities into the disturbed plots territory [80].

6. Conclusions

Under conditions of the same insolation, plots with disturbances in the upper soil horizons and ground cover are accompanied by the formation of long-term surface temperature anomalies. Similar processes were recorded for both natural (post-fire) and post-technogenic landscapes. Within 20 years, the thermal insulation properties of the vegetation recover in the post-fire areas. Thus, the relative temperature anomaly reaches the level of the mathematical error of measurement. In post-technogenic plots, conditions are more “contrast” compared to the background, and the processes of restoration of the thermal regime take a longer period (>60 years). “Neo-technogenic ecosystems” are formed, characterized with special thermal regimes of soils that differ from the background ones both for reclamated (BCM) and for non-reclamated (OMP) plots.

Along with a change in the optical properties of the surface (the spectral and the broadband albedo, the emissivity coefficient), the main reason for the increased surface and soil temperatures of disturbed plots is a decrease in moisture content in the upper soil horizons. Probably, it is this parameter that determines thermal anomalies in disturbed plots as well as their presence at later stages of recovery processes for 10–15 years after a fires impact and for >60 years in post-technogenic ecosystems of Siberia. Thus, monitoring of surface thermal anomalies can be used as a diagnostic criterion of post-technogenic ecosystems state in the context of territories development after technogenic impacts. The data on the moisture content change in the damaged areas are crucially important for the numerical studies of the temperature profiles evolution. However, these data cannot be obtained directly by the remote observation. Further work on this topic involves not only field measurements of moisture content, but also taking into account moisture transfer in a mathematical model.

The destruction of the upper soil layer and the complete degradation of the moss–lichen cover leads not only to increased surface temperatures in summer, but also causes thermal inversion in soil temperature in winter when the thermal resistance in the upper organic layer of the disturbed plots becomes less than in the upper organic horizons of the background soil. This factor can be favorable for the formation of conditions for an earlier exit from the winter state of vegetation and biota due to the increased temperatures of the root-inhabited layer, which is confirmed for Siberia by field studies [9,10,13,56].

The simulation results show the formation of transition zones (“ecotones”) along the periphery of the disturbed plots due to horizontal heat transfer. The typical size of such a zone is ~5.5 m, and the area is up to 6% of the area of disturbed plots. Probably, zones of rapid recovery successions can form along the boundaries of disturbed areas, which makes further study of their properties urgent.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/f12080994/s1, Table S1: weather data for 2018 in Tura (Evenkia region, Central Siberia, Russia); Table S2: results of numerical simulation.

Author Contributions

Conceptualization, K.Y.L., T.V.P., E.I.P. and K.A.F.; methodology, K.Y.L. and E.I.P.; software, K.Y.L., K.A.F., N.D.Y. and A.V.S.; validation, T.V.P., E.I.P. and K.Y.L.; resources, E.I.P. and T.V.P.; data curation, T.V.P.; writing—original draft preparation, E.I.P., T.V.P. and K.Y.L.; writing—review and editing, E.I.P., K.A.F.; visualization, E.I.P., N.D.Y., A.V.S. and K.Y.L.; supervision, T.V.P., E.I.P.; funding acquisition, T.V.P. All authors have read and agreed to the published version of the manuscript.

Funding

This work was performed using the subject of project of IF SB RAS no. 0287-2021-0010. The study of heat transfer in soil was performed using the subject of project of IT SB RAS no. 0257-2021-0001. This research was partly funded by the Russian Foundation for Basic Research (RFBR) and Government of the Krasnoyarsk krai, and Krasnoyarsk krai Foundation for Research and Development Support, no. 20-44-242002 (“Instrumental monitoring of physical properties and numerical modeling of the state of technogenically disturbed soils in Siberia”), and by Siberian Federal University and Government of the Krasnoyarsk krai, and Krasnoyarsk krai Foundation for Research and Development Support, 2020, no. KF-782 49/20 (“Long-term consequences of extreme fires in the permafrost zone of Siberia by the materials of satellite monitoring”).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://ladsweb.modaps.eosdis.nasa.gov/ (accessed on 25 March 2021), The United States Geological Survey (USGS) database (https://earthexplorer.usgs.gov/, accessed on 21 May 2021), Climatic Research Unit (http://www.cru.uea.ac.uk, accessed on 25 March 2021), Weather archive (http://rp5.ru) (accessed on 25 March 2021), and the National Climatic Data Center (NCDC) (http://www7.ncdc.noaa.gov/CDO/cdo, accessed on 25 March 2021), and https://worldview.earthdata.nasa.gov (accessed on 13 April 2021).

Acknowledgments

The satellite data-receiving equipment used was provided by the Centre of Collective Usage of Federal Research Center “Krasnoyarsk Science Center, Siberian Branch of Russian Academy of Sciences”, Krasnoyarsk, Russia.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Models of soil profiles of non-disturbed cryozems (a), post-pyrogenic cryozems (b,c), and technogenic soils (d), which are typical for the study area in the permafrost zone of Siberia with marks of peat (T) and cryoturbated (CR) horizons, permafrost level (C┴), postfire litter O(F+H)pir and mineral horizon (CRpir), coarse humus horizon (Oao), vegetation cover (VC), and underdeveloped organic horizon (W).

Figure 2. Initial images on natural and technogenic disturbed plots: (a) an example of a post-fire plot (burned 2 July–24 July 2018), image from Landsat-8 OLI/TIRS, 12 July 2019; (b) technogenic disturbances at the Borodinsky Coal Mine (BCM), Landsat-8 OLI/TIRS, 02 July 2018; and (c) technogenic disturbances at the Olimpiada Mining Plant (OMP), 12 July 2019. Plots of different ages (from 1970th up to 2018) are marked as Z1–Z5 polygons.

Figure 3. Time series of satellite data on surface temperature in July for (1) post-fire plot in 1997 after 1-years recovery (a), 10-years recovery in 2006 (b), and 20-years recovery in 2016 (c); (2) OMP plot with the identification of zones of different ages (since 2000) formation of a technogenic ecosystem on dumps in (d) 2000, (e) 2011, and (f) 2019; and (3) BCM plot with the identification of zones of different ages (since 1960–1970th) formation of a technogenic ecosystem with reclamation in (g) 2000, (h) 2010, and (i) 2018.

Figure 4. Results of numerical simulation of the heat transfer in soils: (a) temperature at the depth of 0.1 m, (b) the vertical coordinate of the phase change boundary, (c) relative temperature anomaly at the soil surface during summer ΔT/Tbg, °C/°C (%), and (d) annual cycle of the total heat flux on the damaged soil surface (W/m2).

Figure 5. Long-term dynamics of surface temperature anomalies (ΔT/Tbg, %) marked with the period of intensive recovery (Trec) and stabilization time lag (Tstab) averaged for post-fire plots (a) and for post-technogenic plots (b).

Figure 6. Annual temperature dynamics at a depth of 0.1 m: (a) simulation results for the basic variant of model parameters set (see Table 5) and (b) data of field observations in the research area (summarized experimental data from the work in [50]).

Figure 7. Results of numerical modeling of a two-dimensional temperature field (°C) for 200th day (a) and for 294th day (b).

Figure 8. Areas and percentage of the emerging transitional zone (“ecotone”) versus the total area of the disturbed plots. In the calculations, we used the approximation of the circular shape of the disturbed plots. 1 is the percentage of ecotone; 2 is the area of ecotone.

Table 1

Input parameters used for numerical simulation of the temperature regime in the soil profile. Physical characteristics of non-disturbed soils, summarized from the works in [49,50,56,57,58,59].

Horizon	Specific Heat, J/(kg·K)	Volumetric Water/Ice Content (Qw) *	Porosity	$\mathbf{Bulk}\mathbf{Density}\mathbf{of}\mathbf{Dry}\mathbf{Soil},{\mathit{\rho }}_{\mathit{b}},\mathbf{kg}/{\mathbf{m}}^3$	Depth, m
Oao	1500	0.3/0.8	0.91	60	0.0–0.12
O(F+H)	1880	0.35/0.79	0.81	70	0.12–0.17
CR	750	0.275/0.3	0.48	1100	0.17–0.35
C	920	0.3/0.33	0.36	1500	0.35–Zend **

* Constant values of the volumetric water/ice content. ** Z_end—maximum depth of computational domains, m.

Table 2

Input parameters used for numerical simulation of the temperature regime in the soil profile. Physical characteristics of disturbed soils (post-fire and technogenic impact), summarized from the works in [49,50,56,57,58,59].

Horizon	Specific Heat, J/(kg·K)	Volumetric Water/Ice Content *	Porosity	$\mathbf{Bulk}\mathbf{Density}\mathbf{of}\mathbf{Dry}\mathbf{Soil},{\mathit{\rho }}_{\mathit{b}},\mathbf{kg}/{\mathbf{m}}^3$	Depth, m
O(F+H)pir	1880	0.35/0.79	0.81	70	0.12–0.17
CR (CRpir)	750	0.275/0.3	0.48	1100	0.17–0.35
C	920	0.3/0.33	0.36	1500	0.35–Zend **

* Constant values of the volumetric water/ice content. ** Z_end—maximum depth of computational domains, m.

Table 3

Models of thermal conductivity coefficients used for different soil horizons.

Horizon	Models of Thermal Conductivity Coefficients	Parameters Used	Author, Reference
Oao	Equation (5) with ${\lambda }_d=0.05$ and ${\lambda }_{sat}={\lambda }_o{}^{\left(1-\left(w_f+w_s\right)\right)}·{\lambda }_l{}^{w_l}·{\lambda }_s{}^{w_s}$ ${\lambda }_o=0.25;{\lambda }_l=0.56;{\lambda }_s=2.25$	${\lambda }_d$ is the thermal conductivity of dry organic substance, ${\lambda }_o$ is the thermal conductivity of the organic substance in the water-saturated condition, w is the volumetric fraction of a certain state of water, s is the solid state of water, l is the liquid state of water.	Porada et al., 2016 [5]
O(F+H)	${\lambda }_t=0,08·e^{3.88w_f}$ ${\lambda }_{fr}=\left(615w_s+22.2\right)·{10}^{-3}$	$t$ is the thawed soil, fr is the frozen soil	Anisimov et al., 2012 [70]
CR	$\lambda =k\left(0.001·{\rho }_b+{10}^4·w/{\rho }_b-1.1\right)-1.16·{10}^4·w/{\rho }_b$Equation (5) was used, where ${\lambda }_{sat}={\lambda }_o{}^{\left(1-\left(w_f+w_s\right)\right)}·{\lambda }_f{}^{w_f}·{\lambda }_i{}^{w_i}$ ${\lambda }_d=\left(A·{\rho }_b+B\right)/\left({\rho }_{soil}+C·{\rho }_b\right)$Coefficients A, B, С, D summarized from [1,4]	${\rho }_b$ is the bulk density of dry soil;${\rho }_{soil}$ is the density of the material of soil particles; w is the volumetric water/ice fraction; k is the empiric coefficient specific for the type and condition of soil	Pavlov, 1979 [69], Ekici et al., 2014 [4], Tarnawski et al., 2000 [1]

Table 4

Changes in the spectral characteristics of disturbed areas in % to background values for the options PF, BUR (reclamation), and OMP (without reclamation). Mean value for July ± SD, p < 0.05.

Recovery Time, Years	Albedo Decreasing (%)λ = 0.620–0.670 μm		Albedo Decreasing (%)λ = 0.841–0.876 μm		ΔT/Tbg (%)λ = 10–12 μm
	PF/BCM (Reclamation)	OMP (without Reclamation)	PF/BCM (Reclamation)	OMP (without Reclamation)	PF	BCM (Reclamation)/OMP (without Reclamation)
1	17.5 ± 4.8	53.1 ± 11.3	48.5 ± 1.5	46.8 ± 5.2	28.5 ± 3.4	77.9 ± 10.4
5	16.9 ± 5.6	–	27.5 ± 3.6	–	15.0 ± 2.5	–
10–12	15.4 ± 5.5	50.6 ± 8.1	25.0 ± 3.7	51.4 ± 5.6	12.5 ± 1.1	55.0 ± 8.3
22	3.1 ± 0.2	62.9 ± 13.9	12.5 ± 1.1	50.2 ± 5.5	3.6 ± 0.6	43.0 ± 6.2
30	–	–	–	–	–	32.4 ± 4.5
>40	–	–	–	–	–	18.7 ± 0.3

The data obtained for the two considered channels of satellite equipment (two spectral ranges) (Table 4) correspond to the change in the integral albedo of the disturbed plots [74] incorporated in the model of current research as one of the input parameter.

Table 5

Moisture content and albedo in the variants of model problem statement and calculated values of ΔT/T_bg for summer. Moisture content and albedo are given for disturbed/non-disturbed plots.

Horizon	Volumetric Water Content	Albedo for Disturbed Plot/Background	ΔT/Tbg, %
Var. 1. Basic set of conditions
Oao	–/0.6	0.13/0.18	15–25
O(F+H)	0.35/0.14
Var. 2. No change in the moisture content for the horizon O(F+H)
Oao	–/0.6	0.13/0.18	7–12
O(F+H)	0.35/0.35
Var. 3. Equal moisture content in the upper soil horizons
Oao	–/0.6	0.13/0.18	<5
O(F+H)	0.35/0.6
Var. 4. Moisture content similar to Var. 3 with greater change in albedo
Oao	–/0.6	0.07/0.28	10–25 (mean ~15%)
O(F+H)	0.35/0.6