Influence of Rainfall Changes on the Temperature Regime of Permafrost in Central Yakutia

Climate change effects, such as melting of glaciers and sea ice in response to rising temperatures, may lead to an increase in global water availability and thus in precipitation. In Central Yakutia, as one of the possible options for climate change, an increase in rainfall is possible, which makes up more than 60% of the annual precipitation. Rainfall is a highly variable meteorological parameter both spatially and temporally. In order to assess its effect on the ground temperature regime in Central Yakutia, we conducted manipulation and numerical experiments with increased rainfall. The manipulation experiment results suggest that a significant (three-fold) increase in rainfall can lower the mean annual ground temperatures locally. The long-term simulation predicts that a 50% increase in rainfall would have a warming effect on the ground thermal regime on a regional scale. For Central Yakutia, infiltration of increased precipitation has been shown to have both warming and cooling effect depending on the area affected.


Introduction
Heat transfer processes and ground thermal conditions near the surface are governed by the surface energy balance, near-surface air circulation, precipitation, and land cover type, as well as the type, properties, and temperature of underlying material [1]. Precipitation, both solid (snow) and liquid (rain), is important to the ground thermal regime; however, the role of the latter has not been adequately addressed.
The migration of water to the freezing front results in the release of the latent heat of fusion, which reduces the rate and depth of freezing [2,3]. The more water migrates to the freezing front, the shallower is frost penetration depth. Moisture migration in sandy soils is known to be less intensive than in finer-grained sediments, generally resulting in deeper freezing [4][5][6].
Saturation with water and ice has a significant effect on the ground thermal regime and thaw depths. Water percolating from the surface moves downwards under gravity effects. In the permafrost zone, infiltration (in the active layer) is associated with heat transfer from the surface to the frozen-thawed interface. Thus, heat conduction is accompanied by convective heat transport by water [7].
Many studies have assessed the influence of precipitation infiltration on the ground thermal regime and thaw depths, including those by V.V. Veselov [8], V.G. Goldtman [9][10][11], V.T. Balobaev [12], V.A. Kudryavtsev [13], G.Z. Perlshtein [14,15], G.M. Feldman [16], and Iijima et al. [17]. However, no experimental studies with rainfall manipulation have been conducted, except for permafrost pre-thawing by sprinkling for dredging operations. Understanding the contribution of precipitation infiltration to the thermal regime of soils is of considerable interest from a theoretical standpoint. It is also important to engineering and agricultural practice in Central Yakutia [18]. Global climatic change with associated increases in mean annual air temperature, melting of glaciers and rising of sea level may lead to an increase in global water availability and, consequently, in overall precipitation. In Central Yakutia, increases are anticipated both for snowfall and for rainfall, which comprises over 60% of the annual precipitation [19].
In this study, we attempt to determine the effect of rainfall infiltration on the ground thermal regime using theoretical models and experimental methods.

Manipulation Experiment
Two experimental plots were established at the Tuymaada Permafrost Research and Monitoring Station in Central Yakutia to quantify the effect of rainfall infiltration on the thermal regime of sands ( Figure 2). The plots, 3 × 3 m in size, were fenced with thin, solid metal sheets inserted to a depth of 150 cm to prevent lateral water movement.
The geological cross-section between the two plots is shown in Figure 3. The surface is covered with a 15 cm thick organic layer followed by sandy silts to depths of 70-80 cm. Below, fine-and medium-grained sands occur between depths of 70-80 cm and 10 m interbedded with thin layers of sandy silt and plant detritus. The distance between the plots is 15 m. Global climatic change with associated increases in mean annual air temperature, melting of glaciers and rising of sea level may lead to an increase in global water availability and, consequently, in overall precipitation. In Central Yakutia, increases are anticipated both for snowfall and for rainfall, which comprises over 60% of the annual precipitation [19].
In this study, we attempt to determine the effect of rainfall infiltration on the ground thermal regime using theoretical models and experimental methods.

Manipulation Experiment
Two experimental plots were established at the Tuymaada Permafrost Research and Monitoring Station in Central Yakutia to quantify the effect of rainfall infiltration on the thermal regime of sands ( Figure 2). The plots, 3 × 3 m in size, were fenced with thin, solid metal sheets inserted to a depth of 150 cm to prevent lateral water movement.
The geological cross-section between the two plots is shown in Figure 3. The surface is covered with a 15 cm thick organic layer followed by sandy silts to depths of 70-80 cm. Below, fine-and medium-grained sands occur between depths of 70-80 cm and 10 m interbedded with thin layers of sandy silt and plant detritus. The distance between the plots is 15 m.  Both plots are covered with grass vegetation (Figure 4). Two boreholes, 2, 5, and 10 m deep, were drilled at each plot and instrumented with automated (data loggers) and mechanical (thermistors) temperature measurement systems.
HOBO data loggers manufactured by Onset Computer Corporation (Bourne, MA, USA) and data loggers made at the Melnikov Permafrost Institute were used. These loggers have a measurement accuracy of ±0.1 °С and several recording channels (4 to 10) and are suited for research applications [22][23][24][25].  Both plots are covered with grass vegetation (Figure 4). Two boreholes, 2, 5, and 10 m deep, were drilled at each plot and instrumented with automated (data loggers) and mechanical (thermistors) temperature measurement systems.
HOBO data loggers manufactured by Onset Computer Corporation (Bourne, MA, USA) and data loggers made at the Melnikov Permafrost Institute were used. These loggers have a measurement accuracy of ±0.1 °С and several recording channels (4 to 10) and are suited for research applications [22][23][24][25]. Both plots are covered with grass vegetation (Figure 4). Two boreholes, 2, 5, and 10 m deep, were drilled at each plot and instrumented with automated (data loggers) and mechanical (thermistors) temperature measurement systems.
HOBO data loggers manufactured by Onset Computer Corporation (Bourne, MA, USA) and data loggers made at the Melnikov Permafrost Institute were used. These loggers have a measurement accuracy of ±0.1 • C and several recording channels (4 to 10) and are suited for research applications [22][23][24][25].  Mechanical instrumentation consisted of ground temperature cables made at the Melnikov Permafrost Institute with MMT-1 or MMT-4 type semiconductor thermistors. Temperatures were recorded by measuring thermistor resistance using a portable digital multimeter with a resolution of 0.1 Ohm. This method, with lead resistance accounted for, generally gives a temperature error of 0.05 °С, which is acceptable for the scale and status of the present study. The theory and methodology of thermistor application in ground thermal investigations are described in detail by Balobaev et al. [26,27].
The HOBO data loggers collected active-layer temperatures at depths of 0, 10, 20, 40, 70, 90, 120, and 160 cm at three hour intervals. Ground temperatures at depths of 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 m were measured at six hour intervals using the data loggers made in MPI and once every three days at 03:00 p.m. using thermistor cables also made in MPI. The thermistors were also used for cross-checking. The HOBO soil moisture sensors were installed to observe changes in the moisture regime.
A manipulation experiment was set up, with one plot receiving naturally occurring rainfall while the other receiving additional water. The purpose was to compare ground thermal dynamics between the plots in order to estimate the effect of rainfall infiltration on the ground thermal regime.
Prior to the experiment, observations under natural conditions were conducted in 2011-2013. Grain-size and mineral analyses of soils samples were performed. The results showed that the original thermal regimes were identical between the plots, warranting further experimentation.
The experiment was conducted during three summer seasons of 2014-2016. Added water, according to the principle per of 1 mm (rain) + 2 mm (additional water) = 3 mm (total), was watered immediately after rains. To do this, assuming that 10 L of water per area of 1 m 2 obtain a layer of water of 10 mm, we are can calculate how much we need to pour in the amount of water in liters per area of 9 m 2 . The temperature of the added water was adjusted to the air temperature at the time of rain. Climatic data (mean annual air temperature and precipitation) were analyzed for the experimental period (Table 1) [28,29].  Mechanical instrumentation consisted of ground temperature cables made at the Melnikov Permafrost Institute with MMT-1 or MMT-4 type semiconductor thermistors. Temperatures were recorded by measuring thermistor resistance using a portable digital multimeter with a resolution of 0.1 Ohm. This method, with lead resistance accounted for, generally gives a temperature error of 0.05 • C, which is acceptable for the scale and status of the present study. The theory and methodology of thermistor application in ground thermal investigations are described in detail by Balobaev et al. [26,27].
The HOBO data loggers collected active-layer temperatures at depths of 0, 10, 20, 40, 70, 90, 120, and 160 cm at three hour intervals. Ground temperatures at depths of 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 m were measured at six hour intervals using the data loggers made in MPI and once every three days at 03:00 p.m. using thermistor cables also made in MPI. The thermistors were also used for cross-checking. The HOBO soil moisture sensors were installed to observe changes in the moisture regime.
A manipulation experiment was set up, with one plot receiving naturally occurring rainfall while the other receiving additional water. The purpose was to compare ground thermal dynamics between the plots in order to estimate the effect of rainfall infiltration on the ground thermal regime.
Prior to the experiment, observations under natural conditions were conducted in 2011-2013. Grain-size and mineral analyses of soils samples were performed. The results showed that the original thermal regimes were identical between the plots, warranting further experimentation.
The experiment was conducted during three summer seasons of 2014-2016. Added water, according to the principle per of 1 mm (rain) + 2 mm (additional water) = 3 mm (total), was watered immediately after rains. To do this, assuming that 10 L of water per area of 1 m 2 obtain a layer of water of 10 mm, we are can calculate how much we need to pour in the amount of water in liters per area of 9 m 2 . The temperature of the added water was adjusted to the air temperature at the time of rain. Climatic data (mean annual air temperature and precipitation) were analyzed for the experimental period ( For the full-scale determination of the depth of seasonal thawing in the experimental plots, the method of a mechanical probing was used. The thawing depth was determined every ten days.
Soil properties are important parameters that affect heat transfer processes in soils. The primary control over heat transfer are thermal properties, along with density and moisture content.
For investigation of freeze-thaw processes and modeling of thermal interactions between the atmosphere and frozen ground, data on soil thermal conductivity λ [W/(m·K)], soil volumetric heat capacity Cγ [kJ/(m 3 ·K)], and soil thermal diffusivity a [m 2 /c] are needed. These parameters are related by the formula: Soil thermal conductivity was measured using KD2 and KD2 Pro thermal properties analyzers ( Figure 5). The devices have three interchangeable sensors which measure thermal diffusivity, specific heat, thermal conductivity, and thermal resistivity. Measurements can be made in the manual or auto-read mode with on-board and downloadable data storage. For the full-scale determination of the depth of seasonal thawing in the experimental plots, the method of a mechanical probing was used. The thawing depth was determined every ten days.
Soil properties are important parameters that affect heat transfer processes in soils. The primary control over heat transfer are thermal properties, along with density and moisture content.
For investigation of freeze-thaw processes and modeling of thermal interactions between the atmosphere and frozen ground, data on soil thermal conductivity λ [W/(m·K)], soil volumetric heat capacity Сγ [kJ/(m3·K)], and soil thermal diffusivity a [m2/c] are needed. These parameters are related by the formula: Soil thermal conductivity was measured using KD2 and KD2 Pro thermal properties analyzers ( Figure 5). The devices have three interchangeable sensors which measure thermal diffusivity, specific heat, thermal conductivity, and thermal resistivity. Measurements can be made in the manual or auto-read mode with on-board and downloadable data storage.  Needle temperature readings were taken at 1 s intervals during heating or cooling. Raw data were processed using the exponential function of a non-linear least squares procedure. Correction for linear temperature drift was performed during measurements to improve the accuracy of readings.
Control water contents were determined using the standard method. Soil samples for moisture measurement were collected from holes drilled into the active layer with a hand auger. The samples were weighed, oven dried at 105 °С for 24 h and reweighed. The soil moisture content was calculated as the difference between the initial and dry masses [30].

Numerical Experiment
There are two general approaches to mathematical modeling of heat and mass transfer in freezing and thawing soils. One assumes that the phase change is localized at the phase change interface (at a single temperature) and another considers the phase change spread over a significant spatial extent (over a range of temperatures). The latter approach was primarily developed from experimental studies by Tsytovich [31] and Nersesova [32], which demonstrated that these processes occur in an advancing front. Kolesnikov [33] and Martynov [34] proposed a mathematical formulation of the thermal problem for a soil Needle temperature readings were taken at 1 s intervals during heating or cooling. Raw data were processed using the exponential function of a non-linear least squares procedure. Correction for linear temperature drift was performed during measurements to improve the accuracy of readings.
Control water contents were determined using the standard method. Soil samples for moisture measurement were collected from holes drilled into the active layer with a hand auger. The samples were weighed, oven dried at 105 • C for 24 h and reweighed. The soil moisture content was calculated as the difference between the initial and dry masses [30].

Numerical Experiment
There are two general approaches to mathematical modeling of heat and mass transfer in freezing and thawing soils. One assumes that the phase change is localized at the phase change interface (at a single temperature) and another considers the phase change spread over a significant spatial extent (over a range of temperatures). The latter approach was primarily developed from experimental studies by Tsytovich [31] and Nersesova [32], which demonstrated that these processes occur in an advancing front. Kolesnikov [33] and Martynov [34] proposed a mathematical formulation of the thermal problem for a soil freezing over a temperature range. Later, this thermal problem was complemented with a mass transfer equation [1,[35][36][37].
Both approaches are widely used in research and practice and the choice between the two depends on unfrozen water content. Models based on the former approach are generally applied to coarse materials which contain little if any unfrozen water. The latter approach is used in models for fine-grained or saline soils that have appreciable unfrozen contents at temperatures below 0 • C.
These mathematical models differ in the way they consider latent heat and water (ice) and salt released "freeze" with soil moisture.
In this study, we used the model of the phase change over a temperature range, as it more adequately reflects the freezing-thawing process in layered fine-grained media. Under these conditions, free water freezes at T f = 273 • K, while bound water solidifies with decreasing medium temperature [31,32,38,39]. The change of state of the bound water occurs within some range of temperatures (T1 and T2), resulting in the development of a freezing zone.
The energy equation in the temperature range has the form [40]: Soil moisture migration may be written, using the migration model, in the form: Or using the equations of Richards [41] and Van Genuchten [42] for unsaturatedsaturated soils: The system of Equations (2)-(4) is closed by an equilibrium function for unfrozen water content: where c, c w -volumetric heat capacities of the soil and water, J/(m 3 ·K); T-temperature, K; τ-time, s; z-spatial coordinate, m; λ-thermal conductivity of the soil, W/(m·K); Lvolumetric latent heat, J/m 3 ; D-volumetric latent heat of condensation, J/m 3 ; ρ-the bulk density of the mineral skeleton, kg/m 3 ; ρ w -the density of water, kg/m 3 I c = ∂W c /∂τsource of condensation or sublimation; V-flow rate, m/s; k-diffusion coefficient, m 2 /s; θ = θ i + θ w -total volumetric water content (dimensionless), bulk ice θ i and water θ w ; k f -hydraulic conductivity, m/s; h = P − z-pressure, m; P-water suction pressure, m; W = W i + W w -total gravimetric soil moisture content (ice content W i and water content W w ), %. Equation (3) describes the diffusive pore water transport in unsaturated soils, while Richards's Equation (4) defines flow of water in saturated-unsaturated soils. The system of Equations (2)-(4) is non-linear and is solved numerically using an implicit finite difference scheme with iteration.
The model was tested using data from the manipulated experiment. An adapted grid technology for porous media (soils) was applied to estimate changes in temperature and moisture content, as well as thaw progression with account for rainfall infiltration. The problem was implemented with the Delphi 8 software and solved using the finite difference method.
Active layer thickness is an important variable characterizing the state and dynamics of permafrost environments. Seasonal thaw progression was simulated to check the validity of the model. A comparison of the simulated thaw depths with and without account for internal condensation with the measured values is presented in Figure 6. The proposed mathematical model more adequately represents freezing of the "locked' zone and shows a good agreement with measurements, and can reasonably describe the freezing-thawing process in the active layer. It can provide more accurate predictions of thaw depths and better describe seasonal thaw progression over the year.
We can conclude that, based on the field experiments, we developed an improved temperature and moisture model that accounts for rainfall infiltration. The model verification using field data showed that it simulates thawing-freezing processes quite well. At the end of the summer season (late September-early October), the wet thawed layer is "locked". The proposed mathematical model more adequately reflects the freezing process of the "locked" zone. The results demonstrate the adequacy of the improved model and its applicability in environmental change projections and engineering analyses.
A numerical experiment was implemented using the proposed model to predict long-term rainfall effects. Input data were selected to represent a forb meadow site with environmental conditions typical of Central Yakutia. The experiment accounted for heat and moisture transfer by precipitation and evaporation from the surface. Values for mean monthly precipitation, evaporation, air temperature and effective heat transfer coefficient were taken from weather websites, which use Roshydromet data [28,29], climate reference books [43], and field measurements [18].
It should be noted that, on an annual basis, evaporation exceeds precipitation in Central Yakutia [44].
The site lithology is shown in Figure 7. The thermal and mass transfer characteristics for different soil types as a function of temperature, total water content and ice content were set based on the field data obtained at the Tuymada Permafrost Research and Monitoring Station. The proposed mathematical model more adequately represents freezing of the "locked' zone and shows a good agreement with measurements, and can reasonably describe the freezing-thawing process in the active layer. It can provide more accurate predictions of thaw depths and better describe seasonal thaw progression over the year.
We can conclude that, based on the field experiments, we developed an improved temperature and moisture model that accounts for rainfall infiltration. The model verification using field data showed that it simulates thawing-freezing processes quite well. At the end of the summer season (late September-early October), the wet thawed layer is "locked". The proposed mathematical model more adequately reflects the freezing process of the "locked" zone. The results demonstrate the adequacy of the improved model and its applicability in environmental change projections and engineering analyses.
A numerical experiment was implemented using the proposed model to predict long-term rainfall effects. Input data were selected to represent a forb meadow site with environmental conditions typical of Central Yakutia. The experiment accounted for heat and moisture transfer by precipitation and evaporation from the surface. Values for mean monthly precipitation, evaporation, air temperature and effective heat transfer coefficient were taken from weather websites, which use Roshydromet data [28,29], climate reference books [43], and field measurements [18].
It should be noted that, on an annual basis, evaporation exceeds precipitation in Central Yakutia [44].
The site lithology is shown in Figure 7. The thermal and mass transfer characteristics for different soil types as a function of temperature, total water content and ice content were set based on the field data obtained at the Tuymada Permafrost Research and Monitoring Station. On the surface of the soil, the influence of air temperature and precipitation is taken into account: = ( ).
At the lower boundary (at z = H), the conditions of thermal insulation for temperature and impermeability for water content are fulfilled: = 0.
The initial distribution of temperature and total water content are specified as: where T-temperature, K; 0 -initial temperature distribution, K; с -surface temperature, K; τ-time, s; z-spatial coordinate, m; λ-thermal conductivity of the soil, W/(m‧K); k-diffusion coefficient, m 2 /s; = + -total gravimetric soil moisture content (ice content and water content ), %; 0 -initial moisture distribution, %; -effective heat transfer coefficient, W/m 2 ·ּ K; -precipitation infiltration flux, m/s. Model parameterization was done based on data from the experiment and climatic data for the respective period.  On the surface of the soil, the influence of air temperature and precipitation is taken into account: k ∂W w ∂z = q w (τ).
At the lower boundary (at z = H), the conditions of thermal insulation for temperature and impermeability for water content are fulfilled: The initial distribution of temperature and total water content are specified as: where T-temperature, K; T 0 -initial temperature distribution, K; T c -surface temperature, K; τ-time, s; z-spatial coordinate, m; λ-thermal conductivity of the soil, W/(m·K); k-diffusion coefficient, m 2 /s; W = W i + W w -total gravimetric soil moisture content (ice content W i and water content W w ), %; W 0 -initial moisture distribution, %; α e f f -effective heat transfer coefficient, W/m 2 ·K; q w -precipitation infiltration flux, m/s. Model parameterization was done based on data from the experiment and climatic data for the respective period.   The simulated data shows some departure from the field data; nevertheless, the model can be used for long-term projections.

Results of the Manipulation Experiment
Data from the manipulated experiment were analyzed to estimate the effect of rainfall infiltration on the summer thermal regime of sands in Central Yakutia (see Table 1). The differences in mean summer (May-September) ground temperature (∆t) at several depths between the manipulated and control plots are shown in Figure 9.  The simulated data shows some departure from the field data; nevertheless, the model can be used for long-term projections.

Results of the Manipulation Experiment
Data from the manipulated experiment were analyzed to estimate the effect of rainfall infiltration on the summer thermal regime of sands in Central Yakutia (see Table 1). The differences in mean summer (May-September) ground temperature (∆t) at several depths between the manipulated and control plots are shown in Figure 9.  The simulated data shows some departure from the field data; nevertheless, the model can be used for long-term projections.

Results of the Manipulation Experiment
Data from the manipulated experiment were analyzed to estimate the effect of rainfall infiltration on the summer thermal regime of sands in Central Yakutia (see Table 1). The differences in mean summer (May-September) ground temperature (∆t) at several depths between the manipulated and control plots are shown in Figure 9.  Four zones with distinctive ∆t patterns are evident in Figure 9. Zone I between the ground surface and 0.4 m depth is characterized by a ∆t increase in the range of 1.2 to 2.3 • C. This is partly due to the lithological inhomogeneity and occurrence of sandy silts at depths of 0.2 to 0.7 m. Sandy silts have lower hydraulic conductivity than sands; therefore, heat loss from precipitation is significant in this depth interval and the ground receives more heat, resulting in higher temperatures.
Zone II is identified between 0.4 and 1.4 m depths where the temperature difference is strongly reduced. This is due to the change in hydraulic conductivity of the soils in this interval.
Zone III in the depth interval from 1.4 to 2 m is a zone of transition to the quasi-steady regime with a relatively low gradient (∆t). The soil temperature difference decreases with depth by 0.7 to 0.3 • C. This is attributed to the proximity of the frozen-unfrozen interface and the relatively low temperature of infiltrating water at these depths.
Below the depth of 2 m is the zone with a quasi-steady thermal regime. It includes the layer of phase changes and frozen ground.
The data analysis (see Figure 9) suggests that precipitation manipulation increased the maximum temperature difference in the active layer to 1.5 • C. The ∆t peaks occur at depths between 0.4 and 0.9 m due to low infiltration properties of the soils in this interval and a sufficiently long time for heat exchange processes.
The mean annual ground temperature (MAGT) profiles for different years ( Figure 10) demonstrate the effect of rainfall infiltration on the ground thermal regime. The MAGT profiles for the year 2012/2013 (before the experiment) indicate that the temperature difference ∆t between the plots below 1 m was 1 • C and the temperatures were in the range of −1.5 to −2.5 • C. It should be noted that the interannual variation of the MAGT profiles at the natural plot was between −2 and −3 • C. Four zones with distinctive ∆t patterns are evident in Figure 9. Zone I between the ground surface and 0.4 m depth is characterized by a ∆t increase in the range of 1.2 to 2.3 °C. This is partly due to the lithological inhomogeneity and occurrence of sandy silts at depths of 0.2 to 0.7 m. Sandy silts have lower hydraulic conductivity than sands; therefore, heat loss from precipitation is significant in this depth interval and the ground receives more heat, resulting in higher temperatures.
Zone II is identified between 0.4 and 1.4 m depths where the temperature difference is strongly reduced. This is due to the change in hydraulic conductivity of the soils in this interval.
Zone III in the depth interval from 1.4 to 2 m is a zone of transition to the quasi-steady regime with a relatively low gradient (∆t). The soil temperature difference decreases with depth by 0.7 to 0.3 °C. This is attributed to the proximity of the frozen-unfrozen interface and the relatively low temperature of infiltrating water at these depths.
Below the depth of 2 m is the zone with a quasi-steady thermal regime. It includes the layer of phase changes and frozen ground.
The data analysis (see Figure 9) suggests that precipitation manipulation increased the maximum temperature difference in the active layer to 1.5 °C. The ∆t peaks occur at depths between 0.4 and 0.9 m due to low infiltration properties of the soils in this interval and a sufficiently long time for heat exchange processes.
The mean annual ground temperature (MAGT) profiles for different years ( Figure  10) demonstrate the effect of rainfall infiltration on the ground thermal regime. The MAGT profiles for the year 2012/2013 (before the experiment) indicate that the temperature difference ∆t between the plots below 1 m was 1 °С and the temperatures were in the range of −1.5 to −2.5 °С. It should be noted that the interannual variation of the MAGT profiles at the natural plot was between −2 and −3 °С.   The precipitation manipulation increased soil moisture contents and subsequently decreased the depths of thaw ( Figure 11) due to the increase in the layer of ice-saturated soils with higher latent heat. As a result, ground temperatures decreased. The assumption that precipitation falls uniformly over the entire surface and evaporates similarly is inappropriate, since this is influenced by differences in vegetation and lithology. Commonly, these parameters are highly variable and evaporation, as an output component and thermodynamic process, plays an important role in the thermal balance of landscapes.
Land 2021, 10, x FOR PEER REVIEW 11 of 20 The precipitation manipulation increased soil moisture contents and subsequently decreased the depths of thaw ( Figure 11) due to the increase in the layer of ice-saturated soils with higher latent heat. As a result, ground temperatures decreased. The assumption that precipitation falls uniformly over the entire surface and evaporates similarly is inappropriate, since this is influenced by differences in vegetation and lithology. Commonly, these parameters are highly variable and evaporation, as an output component and thermodynamic process, plays an important role in the thermal balance of landscapes. The pattern and depths of thaw penetration (see Figure 11) indicate that the manipulated plot thawed to a depth of 1.5 m faster in 2014 and 2015 (first two years of the experiment) than in 2013 (control year). However, thaw penetration (to a depth of 0.5 m) occurs slower in the early warm season of the second year (2015). During the third year, thaw penetration patterns are nearly identical with some lag in the upper 0.6 m of soil at the manipulated plot. This is attributed to increased saturation of the soils after additional water treatment. Some decrease in thaw depths since 2015 was due to higher natural rainfall in 2014 and relatively low rainfall quantity in 2015.
To obtain representative soil thermal properties, in-situ measurements of thermal conductivity, thermal diffusivity and volumetric heat capacity were performed during the summers of 2015 and 2017. In addition, soil samples were collected from the upper 1.5 m of soil for determination of moisture content, density and grain size at different levels ( Figure 12). The measurements were made once monthly from May to September (five times a year). The pattern and depths of thaw penetration (see Figure 11) indicate that the manipulated plot thawed to a depth of 1.5 m faster in 2014 and 2015 (first two years of the experiment) than in 2013 (control year). However, thaw penetration (to a depth of 0.5 m) occurs slower in the early warm season of the second year (2015). During the third year, thaw penetration patterns are nearly identical with some lag in the upper 0.6 m of soil at the manipulated plot. This is attributed to increased saturation of the soils after additional water treatment. Some decrease in thaw depths since 2015 was due to higher natural rainfall in 2014 and relatively low rainfall quantity in 2015.
To obtain representative soil thermal properties, in-situ measurements of thermal conductivity, thermal diffusivity and volumetric heat capacity were performed during the summers of 2015 and 2017. In addition, soil samples were collected from the upper 1.5 m of soil for determination of moisture content, density and grain size at different levels ( Figure 12). The measurements were made once monthly from May to September (five times a year). These measurements provided data on soil physical and thermal properties for geocryological and microclimatic conditions typical of Central Yakutia. The data obtained are presented in Table 2 for the soil types determined by grain-size analysis. The data indicate that for the range of dry density between 1640 and 1740 kg/m 3 and the range of moisture content between 5.3 and 14.5%, the thermal properties of the soils vary from 0.4 to 1.26 W/(m·K) for thermal conductivity, from 0.29·10 -6 to 0.84·10 -6 m 2 /s for thermal diffusivity and from 1280 to 1650 kJ/(m 3 ·K) for volumetric heat capacity.

Results of the Numerical Experiment
The proposed model was used to predict changes in MAGT with changing rainfall, using results of the manipulation experiment.
For this purpose, the climatic data for Central Yakutia, including rainfall amount, were analyzed. Data for Yakutsk were taken from climate reference books and open sources [28,29,43] (Figure 13). These measurements provided data on soil physical and thermal properties for geocryological and microclimatic conditions typical of Central Yakutia. The data obtained are presented in Table 2 for the soil types determined by grain-size analysis. The data indicate that for the range of dry density between 1640 and 1740 kg/m 3 and the range of moisture content between 5.3 and 14.5%, the thermal properties of the soils vary from 0.4 to 1.26 W/(m·K) for thermal conductivity, from 0.29·10 -6 to 0.84·10 -6 m 2 /s for thermal diffusivity and from 1280 to 1650 kJ/(m 3 ·K) for volumetric heat capacity.

Results of the Numerical Experiment
The proposed model was used to predict changes in MAGT with changing rainfall, using results of the manipulation experiment.
For this purpose, the climatic data for Central Yakutia, including rainfall amount, were analyzed. Data for Yakutsk were taken from climate reference books and open sources [28,29,43] (Figure 13). The results showed that interannual variability of rainfall was within the range of long-term normal (160 mm ± 50%), except for 2001 and 2006 when the amount of rainfall was anomalous (see Figure 13). For this reason, the following scenarios were used in modeling MAGT changes: normal +50%; normal +25%; control (normal); normal −25%; normal -50%. The only variable changing was rainfall, all others being equal. For correct calculation, long-term means were used for the initial and boundary conditions. Projection time was 20 years. The predicted ground temperatures at depths of 3, 6, and 10 m in 10 years are given in Figure 14.  The results showed that interannual variability of rainfall was within the range of long-term normal (160 mm ± 50%), except for 2001 and 2006 when the amount of rainfall was anomalous (see Figure 13). For this reason, the following scenarios were used in modeling MAGT changes: normal +50%; normal +25%; control (normal); normal −25%; normal -50%. The only variable changing was rainfall, all others being equal. For correct calculation, long-term means were used for the initial and boundary conditions. Projection time was 20 years. The predicted ground temperatures at depths of 3, 6, and 10 m in 10 years are given in Figure 14. The results showed that interannual variability of rainfall was within the range of long-term normal (160 mm ± 50%), except for 2001 and 2006 when the amount of rainfall was anomalous (see Figure 13). For this reason, the following scenarios were used in modeling MAGT changes: normal +50%; normal +25%; control (normal); normal −25%; normal -50%. The only variable changing was rainfall, all others being equal. For correct calculation, long-term means were used for the initial and boundary conditions. Projection time was 20 years. The predicted ground temperatures at depths of 3, 6, and 10 m in 10 years are given in Figure 14. The simulation results suggest that under the increased rainfall scenarios ground temperatures will warm, while the decreased scenarios lead to a cooling effect (see Figure  14). A decrease or increase in precipitation by 50% after 10 years will cause MAGT at 3 m depth to decrease by 0.37 °C and increase by 0.53 °C, respectively.
The range (amplitude) of MAGT changes tends to decrease with depth for all precipitation scenarios. For example, under a Normal ±50% scenario, the amplitude in 10 years is 0.9 °С at a depth of 3 m and 0.53 °C at a depth of 10 m. Under a Normal ±25% scenario, the amplitude in 10 years is 0.45 °С at 3 m and 0.26 °С at 10 m.
It should be noted that the warming effect of increased rainfall is greater than the cooling effect of reduced rainfall. This is attributed to the difference in thermal properties of water-saturated and dry soils, latent heat of phase changes, and evaporation rates which are directly dependent on soil moisture content [17].

Discussion
The MAGTs at the manipulated plot changed over a wide range (see Figure 10). A strong increase, up to -1 °С and higher, occurred at all depths during the first year of the experiment. This was attributed to additional heat input due to higher precipitation. It should be noted that the year was warmest over the observation period; increased heat input with increased precipitation was nevertheless evident. Considering that the difference in MAGTs between the plots was 0.5 to 1.0 °С prior to the experiment, the effect of rainfall infiltration (less the above difference) as a result of the experiment was 0.2 to 1.7 °С. These values are close to theoretical estimates derived by Kudryavtsev [13] and Varlamov et al. [45].
In subsequent years of the experiment, however, the mean annual temperatures at depth began to recover and lowered, on average, down to -2 °C Compared to MAGT at the control plot where it was directly related to climatic parameters (air temperature and precipitation regime), MAGT at the manipulated plot decreased each year (see Figures 10  and 15). It appears that MAGTs may decrease with increasing amounts of rainfall. In other words, other factors being equal, a long-term increase in rainfall in local areas will have a cooling effect on the ground temperature regime in Central Yakutia. A similar conclusion was reached by Zhi Wen et al. [46] for the Tibetan Plateau. The simulation results suggest that under the increased rainfall scenarios ground temperatures will warm, while the decreased scenarios lead to a cooling effect (see Figure 14). A decrease or increase in precipitation by 50% after 10 years will cause MAGT at 3 m depth to decrease by 0.37 • C and increase by 0.53 • C, respectively.
The range (amplitude) of MAGT changes tends to decrease with depth for all precipitation scenarios. For example, under a Normal ±50% scenario, the amplitude in 10 years is 0.9 • C at a depth of 3 m and 0.53 • C at a depth of 10 m. Under a Normal ±25% scenario, the amplitude in 10 years is 0.45 • C at 3 m and 0.26 • C at 10 m.
It should be noted that the warming effect of increased rainfall is greater than the cooling effect of reduced rainfall. This is attributed to the difference in thermal properties of water-saturated and dry soils, latent heat of phase changes, and evaporation rates which are directly dependent on soil moisture content [17].

Discussion
The MAGTs at the manipulated plot changed over a wide range (see Figure 10). A strong increase, up to -1 • C and higher, occurred at all depths during the first year of the experiment. This was attributed to additional heat input due to higher precipitation. It should be noted that the year was warmest over the observation period; increased heat input with increased precipitation was nevertheless evident. Considering that the difference in MAGTs between the plots was 0.5 to 1.0 • C prior to the experiment, the effect of rainfall infiltration (less the above difference) as a result of the experiment was 0.2 to 1.7 • C. These values are close to theoretical estimates derived by Kudryavtsev [13] and Varlamov et al. [45].
In subsequent years of the experiment, however, the mean annual temperatures at depth began to recover and lowered, on average, down to -2 • C Compared to MAGT at the control plot where it was directly related to climatic parameters (air temperature and precipitation regime), MAGT at the manipulated plot decreased each year (see Figures 10 and 15). It appears that MAGTs may decrease with increasing amounts of rainfall. In other words, other factors being equal, a long-term increase in rainfall in local areas will have a cooling effect on the ground temperature regime in Central Yakutia. A similar conclusion was reached by Zhi Wen et al. [46] for the Tibetan Plateau. In summary, increased rainfall treatment caused an increase in summer ground temperatures by up to 1.5 °С. The temperature difference between the manipulated and control plots was greatest in the depth interval between 0.4 and 0.9 m. A short-term increase in precipitation by 3 times increased ground temperatures in the layer of annual temperature fluctuations by 0.2-1.7 °С on an annual basis. A longer-term (3 or more years) increase in precipitation can locally have a cooling effect on permafrost temperatures in Central Yakutia. Figure 16 present the predicted MAGT profiles within the top 2 m-of soil and 10 m-of ground in 10 and 20 years, other factors being equal. In summary, increased rainfall treatment caused an increase in summer ground temperatures by up to 1.5 • C. The temperature difference between the manipulated and control plots was greatest in the depth interval between 0.4 and 0.9 m. A short-term increase in precipitation by 3 times increased ground temperatures in the layer of annual temperature fluctuations by 0.2-1.7 • C on an annual basis. A longer-term (3 or more years) increase in precipitation can locally have a cooling effect on permafrost temperatures in Central Yakutia. Figure 16 present the predicted MAGT profiles within the top 2 m-of soil and 10 mof ground in 10 and 20 years, other factors being equal.
The model predicts that a 50% reduction in rainfall will lower MAGT at the permafrost table (2 m depth) by 0.4 • C in 10 years and by 1.1 • C in 20 years. For 10 m depth, ground cooling will be by 0.2 and 0.8 • C, respectively. A 50% increase in precipitation will warm MAGT at 2 m depth by 0.5 • C in 10 years and by 1. The model predicts that a 50% reduction in rainfall will lower MAGT at the permafrost table (2 m depth) by 0.4 °С in 10 years and by 1.1 °С in 20 years. For 10 m depth, ground cooling will be by 0.2 and 0.8 °С, respectively. A 50% increase in precipitation will warm MAGT at 2 m depth by 0.5 °С in 10 years and by 1.2 °С in 20 years. MAGT at 10 m will warm by 0.3 °С in 10 years and 0.5 °C in 20 years. Hence, the ±50% change in precipitation relative its normal results in the range of MAGT between −1.7 °С and −4.0 °С at 2 m and between −1.8 °С and −3.0 °С at 10 m depth.
The simulation results show that an increase in summer precipitation by half of its long-term normal during 10 years will increase MAGT by −1.9 to −2.4 °С at various depths. After 20 years, these values can reach −1.0 to −1.8 °С. If precipitation decreases by 50%, ground temperatures can decrease by−2.5 to −3.4 °С after 10 years and by −3.0 to −4.3 °С The simulation results show that an increase in summer precipitation by half of its long-term normal during 10 years will increase MAGT by −1.9 to −2.4 • C at various depths. After 20 years, these values can reach −1.0 to −1.8 • C. If precipitation decreases by 50%, ground temperatures can decrease by−2.5 to −3.4 • C after 10 years and by −3.0 to −4.3 • C in 20 years. Thus, the threshold values of MAGT changes have been obtained for different rainfall change scenarios.
Thus, the cooling and warming effect of infiltration with an increase in the amount of liquid atmospheric precipitation in the conditions of Central Yakutia has been established, depending on the area of rainfall ( Figure 17). in 20 years. Thus, the threshold values of MAGT changes have been obtained for different rainfall change scenarios.
Thus, the cooling and warming effect of infiltration with an increase in the amount of liquid atmospheric precipitation in the conditions of Central Yakutia has been established, depending on the area of rainfall ( Figure 17). The differences between the observed and simulated effects can be attributed to different rainfall inputs. In the manipulation experiment, rainfall was increased threefold locally in a limited area (9 m 2 ). The water which infiltrated through the active layer froze completely due to the surrounding and underlying permafrost and winter cooling. The process repeated over 3 years of the experiment resulted in an ice-saturated layer, which thawed slower, decreasing the active-layer thickness. In this case, convective heat transfer by the rainfall was small and had no significant effect on the total heat balance, thus the cooling effect was observed. The numerical simulation assumed a 50% increase in rainfall, up to 240 mm in summer season (May to September), in a semi-infinite medium. In this case, a large quantity of water enters the system which introduces a large amount of heat, changing the total heat balance towards warming.
It can be concluded that, in the general case, increased rainfall would have a warming effect on the ground thermal regime. However, an increase in summer precipitation in limited, localized areas can lower ground temperatures, decreasing the active-layer thickness and restoring the transient ice-saturated layer. This may explain the stability of permafrost conditions observed in some areas despite the increasing mean annual air temperatures.

Conclusions
Analysis of the observational and experimental data from the Central Yakutia site with sandy grounds has shown that rainfall timing and amounts have a significant effect on the temperature regime of upper permafrost: - The timing of rainfall varies from year to year in the region. The effect of rainfall infiltration on the ground temperature regime is greatest if rainfall intensity is high in the fall. This increases moisture contents and, subsequently, ice contents of the soils. Thaw depths during the next summer are reduced due to the latent heat effects and hence ground temperatures are lowered. If precipitation amounts are higher in the early and middle summer, no significant changes to ground temperatures or thaw depths are observed, since a large portion of the rainfall evaporates. -A single short-term (one season) increase in rainfall increases upper permafrost temperatures by 0.2 to 1.7 °С on an annual basis. The differences between the observed and simulated effects can be attributed to different rainfall inputs. In the manipulation experiment, rainfall was increased threefold locally in a limited area (9 m 2 ). The water which infiltrated through the active layer froze completely due to the surrounding and underlying permafrost and winter cooling. The process repeated over 3 years of the experiment resulted in an ice-saturated layer, which thawed slower, decreasing the active-layer thickness. In this case, convective heat transfer by the rainfall was small and had no significant effect on the total heat balance, thus the cooling effect was observed. The numerical simulation assumed a 50% increase in rainfall, up to 240 mm in summer season (May to September), in a semi-infinite medium. In this case, a large quantity of water enters the system which introduces a large amount of heat, changing the total heat balance towards warming.
It can be concluded that, in the general case, increased rainfall would have a warming effect on the ground thermal regime. However, an increase in summer precipitation in limited, localized areas can lower ground temperatures, decreasing the active-layer thickness and restoring the transient ice-saturated layer. This may explain the stability of permafrost conditions observed in some areas despite the increasing mean annual air temperatures.

Conclusions
Analysis of the observational and experimental data from the Central Yakutia site with sandy grounds has shown that rainfall timing and amounts have a significant effect on the temperature regime of upper permafrost: - The timing of rainfall varies from year to year in the region. The effect of rainfall infiltration on the ground temperature regime is greatest if rainfall intensity is high in the fall. This increases moisture contents and, subsequently, ice contents of the soils. Thaw depths during the next summer are reduced due to the latent heat effects and hence ground temperatures are lowered. If precipitation amounts are higher in the early and middle summer, no significant changes to ground temperatures or thaw depths are observed, since a large portion of the rainfall evaporates. -A single short-term (one season) increase in rainfall increases upper permafrost temperatures by 0.2 to 1.7 • C on an annual basis. - The experimental results indicate that a longer-term (over 3 years) increase in rainfall, in Central Yakutia areas with sandy grounds, leads to water (ice) accumulation in the near-surface layers which has a cooling effect on the temperature regime of upper permafrost. This is related to the formation of moisture, and later ice, anomalies due to water saturation, changes in soil thermal properties, and latent heat effects. This explains the stability of permafrost in some local areas even under increasing air temperatures. The numerical simulation of future changes in the ground thermal regime suggests that increased rainfall will have a warming effect regionally.

Conflicts of Interest:
The authors declare no conflict of interest.