Salinity effect of permafrost temperature predictions by the example of the Yamal Peninsula

Nikitin Kirill Alekseevich Postgraduate student; Faculty of Geology; Lomonosov Moscow State University 119991, Russia, Moscow, Leninskie Gory str., 1 Nikitin.kirill@yandex.ru Komarov Il'ya Arkadevich Doctor of Geology and Mineralogy Professor; Faculty of Geology; Lomonosov Moscow State University 119991, Russia, Moscow, Leninskie Gory str., 1 ilya_komarov@mail.ru Mironenko Mikhail Viktorovich PhD in Geology and Mineralogy Senior Researcher; Vernadsky Institute of Geochemistry and Analytical Chemistry of the Russian Academy of Sciences 19 Kosygina str., Moscow, 119991, Russia mikhail_mironenk@mail.ru

Kiyashko Nadezhda Vladimirovna PhD in Geology and Mineralogy Researcher; Lomonosov Moscow State University 119991, Russia, Moscow, Leninskie Gory str., 1 nadin130187@mail.ru

Introduction

Saline permafrost soils are the basis for most of the designed and operated buildings and structures on the Yamal Peninsula. Under the conditions of climatic changes in the Yamalo-Nenets Autonomous Okrug, there is an increase in the number of deformations and destructions of foundations and superfundament structures of buildings and engineering structures designed primarily according to the first principle, that is, with preservation of the frozen state ^[1]. The predicted activation of cryogenic geological processes associated with climate change and the thermal humidity regime of the upper horizons of frozen rocks has an additional negative impact on the technical condition of infrastructure and business conditions in Yamal ^[2-5].

The salinity of frozen soils places additional requirements not only on the operational reliability and safety of buildings and structures erected on such soils, but also on the conduct of geocryological forecasting. Changes in thermophysical, water-physical, strength and deformation properties associated with the dynamics of the temperature regime of rocks are much more pronounced in the presence of easily soluble salts ^{[6, 7]}. Quantitative assessment of rock salinity is important when conducting predictive thermal engineering calculations, especially taking into account the effect of climate change.

A significant amount of diverse empirical information has been accumulated (field observations, experimental, laboratory studies) concerning the characteristics of saline rocks, as a result of which a certain gap has formed between the depth and detail of developments to clarify the physical nature of processes and models embedded in the schemes of geocryological forecasting ^[8].

The refinement of forecast estimates can occur by increasing the accuracy of the input empirical information on the characteristics of saline rocks, creating object-oriented databases of experimental data and regional regulatory documents that take into account the specifics of the composition, structure and properties of rocks.

In this regard, the relevance of developing a methodology and forecasting the temperature regime of saline frozen rocks, taking into account climate changes and the water-ion composition of pore solutions accompanied by chemical interactions, increases. The change in the water-ion composition is considered as a result of physico-chemical reactions of the transformation of dissolved gases, mineral and organic substances in solution as a result of the formation and melting of ice. The main purpose of the work is a quantitative study of the effect of salts in pore solutions of frozen rocks on the forecast estimates of their temperature regime by the middle of the century on the example of Yamal. The problem has been solved for three areas of the peninsula — northwestern and western Yamal, the lower reaches of the Ob, within which the rocks are characterized by a marine type of salinity.

Research area

The territory of the Yamal Peninsula is mainly a flat, terraced accumulative lowland plain located in the north of the West Siberian Plate. Long-term regional studies have established that saline frozen rocks north of the latitude of the mouth of the Baydaraty river—settlement. The new Port has a continuous distribution, salts are present throughout the section of the cryogenic strata ^[9]. South of the selected boundary, the sediments are unsalted or slightly saline from the surface to a depth of 100 m. The salt content in rocks depends on the facies conditions of their formation and age. The highest salinity values are observed in modern coastal areas — in Holocene Laide and alluvial estuarine sediments. The continental Late Pleistocene-Holocene sediments are unsalted or slightly saline. In marine clay sediments that retain primary sedimentation salinity, the salt content can reach 3%, exceeding the salt content in the sands several times, except for sandy horizons with cryopags ^{[10, 11]}.

The selected areas for forecasting the temperature regime of rocks differ in patterns in the distribution of salinity and mineralization of pore solutions. In the upper horizons of marine, coastal-marine Late Pleistocene and modern sediments composing watersheds in northwestern Yamal, the salt content in loams is 0.4-1.1%, in sands and sandy loams no more than 0.2% ^[12]. The composition of the solutions is mainly sodium chloride. The salt concentration in section I of the marine terrace reaches 100 g/l in clay rocks and 20 g/l in sandy rocks ^{[13, 14]}. On the territory of western Yamal, the salinity of marine Late Pleistocene loamy-clay rocks is in the range of 0.2—2.5%, sandy 0,05—0,2% ^{[15, 16]}. In the layer of annual temperature fluctuations, there is an increase in salt concentration by 1.5—2 times with depth. The type of soil salinization is predominantly chloride. In the lower reaches of the Ob River, the greatest salinity is observed in sections of lagoon-marine terraces. Salinity reaches 0.1—1.1% in loams and 0.05—0.1% in sands ^{[10, 17]}. The composition of the solutions is mainly chloride-sulfate-sodium.

Saline frozen soils containing lenses of unfrozen, negative-temperature brines — cryopags are widespread in the studied regions ^[18]. Such lenses have been observed not only in modern marine and alluvial-marine sediments, but also in marine Pleistocene deposits of terraces ^{[19, 20]}. The cryopags of the terraces lie in the section in the form of low-power lenses and are confined to areas with elevated soil temperatures. On the slides, the deposits containing cryopags are cooled. Cryopegs are characterized by the highest mineralization up to 150 g/l. In cryopegs, compared with pore solutions of frozen soils, there is an increased content of chlorine and sodium ions and a reduced content of sulfate ion ^[16].

Climate change with a warming trend is observed in the north of Western Siberia. The results of ground-based observations show an increase in the average annual air temperature by an average of 1.4 ° C relative to the climatic norm, together with an increase in the duration of the warm period. The height of the snow cover gradually increases at a rate of up to 2 cm/year, with significant interannual variability ^{[21, 22]}. These parameters significantly affect the temperature of the rocks. Long–term regime observations at hospitals in the typical tundra zone of Yamal have established that since the late 1970s, the average annual temperature of rocks has increased by 1.5—2.2 °C ^{[22, 23]}. Since 2010, the linear trend of increasing the maximum depth of seasonal thawing has changed from 0.8 cm/year (Bely Island) to 7.3 cm/year (Yerkuta) ^[24].

Materials and methods

A feature of the developed methodology for predicting the temperature regime of saline frozen rocks is to take into account changes in their phase and chemical composition during freezing—thawing against the background of climatic changes using thermodynamic and mathematical modeling methods performed in the "Freezbrine" and "QFrost" programs ^{[25, 26]}. For each of the districts, the calculation was performed for two types of model sections — sandy and loamy. Sections of elevated terraced plots were selected within each of the three regions. The calculation was performed for several values of the mineralization of the pore solution — 35, 90 and 150 g/l.

1. Thermodynamic modeling of the transformation of the water-ion composition of saline frozen rocks

The Freezbrine program is a continuation and improvement of the thermodynamic Frezchem model developed at The CRREL laboratory (The Cold Regions Research and Engineering Laboratory, USA) as part of the NASA (National Aeronautics and Space Administration) project to study the behavior of salt solutions on other planets. The Frezchem model is designed to calculate the chemical equilibrium between aqueous solutions of electrolytes, ice and salts, using the methods of statistical thermodynamics by K. Pitzer ^{[27, 28]}.

The database of the Freezbrine program includes 8 cations, 7 anions, 8 neutral particles, 8 gases and 56 solid components. The input data to the program are: the initial total mineralization and concentration of the main components of the chemical composition of the sample obtained on the basis of chemical analysis. The Gibbs free energy minimization method in multicomponent aqueous-salt solutions is used to calculate the phase composition and content of salts, ice, gas and liquid phases, osmotic coefficient, ionic strength and pH of solutions, chemical potential, etc. in a wide range of temperature values ^{[29, 30]}.

The following assumptions were used in predicting the water-ion regime: the thermophysical properties of frozen, thawed and cooled zones are piecewise constant, there is no chemical and adsorption interaction of pore moisture; ion transfer due to diffusion and convection is negligible; thermal effects associated with the precipitation of salts from solution into the solid phase and the precipitation (absorption) of the heat of dissolution, they are small compared to the heat of the water—ice phase transitions.

Thermodynamic modeling of the chemical transformation of salt solutions during freezing—thawing was performed on the basis of published data ^{[16, 17, 31-33]} and stock materials of engineering surveys of NTF Krios LLC. Information was used on the water-physical, thermophysical properties, the degree of salinity of rocks, the content of basic ions and mineralization of pore solutions of various stratigraphic-genetic complexes composing the upper horizons of the cryogenic strata of Yamal.

2. Mathematical modeling of the temperature regime of saline frozen rocks

Mathematical modeling of the dynamics of the temperature regime of saline rocks, taking into account the trend of increasing air temperature, was performed in the QFrost program. The predictive task is one—dimensional, that is, it is solved for a homogeneous half-space under constant boundary conditions, mass transfer characteristics of soil in thawed and frozen states and the dependence of the content of unfrozen water on temperature.

The mathematical formulation of the Stefan problem at the boundary of frozen and permafrost rocks has the following form:

where is the thermal conductivity of the soil in the freezing zone, W/(m*°C); — thermal conductivity of thawed soil, W/(m*°C); — the effective heat capacity of the soil in the freezing zone, J/(kg * °C); — the heat capacity of thawed soil, J/(kg * °C); t is the temperature, °C; ξ is the coordinate of the phase boundary; l is the depth of the calculated area, m; z is the coordinate; τ is the time, h.

The following conditions are met on the mobile boundary between thawed and frozen rocks:

where is the temperature of the phase transition, ° C; where is the latent heat of phase transitions per unit volume of soil, J/m3.

An explicit regularization solution scheme was applied, which made it possible to increase the time step several times and reduce the counting time by the same number of times. The enthalpy form of the problem and the balance method were used, which led to a visual physical interpretation of the results and significantly increased the accuracy of determining the position of the freezing—thawing boundaries. The dynamics of the temperature regime of rocks is estimated by the average annual temperature of rocks, determined at the bottom of the annual heat turnover layer.

The calculation in the QFrost program was performed in two independent iterations. In the first case, the program used constant heat exchange characteristics of soils for the forecast period. In the second case, these characteristics were calculated using the "Freezbrine" program with subsequent data processing according to the methodology proposed by I. A. Komarov ^{[34, 35]}.

Additionally, the problem of determining the time step at which a statistically noticeable change in the thermophysical properties of saline rocks was solved. In real natural conditions, such changes are noticeable with significant fluctuations in soil temperature, manifested on the scale of years or seasons of freezing—thawing. It was previously established that the concentration of pore solution of rocks in the layer of annual temperature fluctuations during the year can vary several times ^{[16, 36]}. The maximum increase in concentration is observed mainly at the end of the freezing period, and the minimum at the end of the thawing period. Thus, the heat exchange characteristics of the rocks were changed in the program 2 times a year. The boundary between these periods is determined by the sign (positive or negative) of the average monthly temperature at the upper boundary of the calculated area. A constant time step of 3 hours was adopted for all iterations.

3. Development of a scenario for changes in the average annual air temperature for the purposes of geocryological forecasting

Scenarios of changes in the average annual air temperature have been developed to predict the temperature regime of frozen rocks by the middle of the century. An empirical approach based on the analysis of the results of long-term observations at weather stations was used. The method is based on the identification of cyclicity, the allocation of multi-period oscillations of different genesis. A set of rhythms overlapping each other with different periods, amplitudes and phase shifts determine the course of the parameter under consideration.

Based on the harmonic analysis of meteorological observations of air temperature, L. N. Khrustalev proposed a method of author's retrospective analysis ^[37]. The method is based solely on the data of instrumental observations conducted with high frequency according to a single methodology. Periodic fluctuations in air temperature over the base interval are approximated by the trigonometric Fourier series. Long-term series of observations of the average monthly air temperature were used, available in the open database of the All—Russian Scientific Research Institute of Hydrometeorological Information - the World Data Center ^[21]. Based on this methodology, scenarios of changes in the average annual air temperature have been developed for three regions — Marre-Sale (western Yamal), Kharasaway (northwestern Yamal), Salekhard (lower Ob river).

Creation, verification and implementation of a thermal model

When creating a mathematical model in the QFrost program, boundary conditions were set. A type II condition has been established at the lower boundary — a time-constant heat flux of 0.06 W/m2 ^[38]. There is no heat flow on the side faces. In accordance with the condition of the third kind, the heat transfer coefficient is set at the upper boundary, having the form:

where is the coefficient of convective heat exchange of the soil surface with air, depending on the wind speed, W /(m2 * ° C), is determined by the formula ^[39]:

where is the average monthly air velocity, m/s;

— total thermal resistance of snow in winter, (m2 * °C)/W, is from the expression:

where δ is the thickness of the layer, m.

For the forecast period, the average annual values of snow height for 1970-2020 were used ^[21]. This approach is justified by the fact that the problem considered sections of elevated flat watersheds, where the predicted increase in snow height can be minimized due to wind redistribution. The influence of vegetation on elevated areas was not taken into account. The thermal conductivity of the snow column is determined by the empirical dependence of B. V. Proskuryakov ^[38]:

where is the density of the snow cover, g/cm3.

The organomineral layer (modern soil), which has a small capacity, is excluded from the calculation area. For the forecast period, the possible precipitation of thawing rocks and the development of cryogenic processes are not taken into account. The size of the calculated area is 20 m, which corresponds to or exceeds the depth of the sole of the annual heat turnover layer.

Calibration was carried out to obtain a model of a watered rock mass corresponding to real natural conditions. Its main goal is to find boundary conditions that form a temperature field corresponding to the natural one for the entire computational domain. The initial temperature field is set based on the results of thermometric observations made by NTF Krios LLC mainly in the 1990s.

Monthly values of temperature and heat transfer coefficient averaged over the period from 1970 to the year of thermometric observations were used at the upper boundary. 1970 was chosen as the beginning of the averaging period, since since the early 1970s of the XX century, according to meteorological observations, there has been a directed increase in the average annual air temperature in Yamal. The resulting calculation result is accepted as the initial conditions of the main task.

Results

Based on the numerical solution of the Stefan problem, the values of the average annual temperature of frozen rocks up to and including 2050 for three regions of Yamal were determined. The calculation was performed for two types of model sections with three values of mineralization of pore solutions. In the calculation process, climate warming was taken into account based on the developed scenarios for changes in the average annual air temperature.

For northwestern Yamal, by the middle of the century, the predicted value of the linear trend of air temperature increase is 0.3 °C/10 years based on the processing of long-term observations at the Kharasaway weather station. The predicted temperature of sandy rocks will change from -5.0...-6.8 °C to -4.8...-6.0 °C, loamy rocks from -5.3...-6.9 °C to -5.2...-6.2 °C. The linear trend of change reaches 0.2 °C/10 years (Fig. 1).

Figure 1. Temperature change of saline frozen rocks of northwestern Yamal.

1 and 2 are the results of solving the problem with and without taking into account the transformation of the phase and chemical composition of pore solutions, respectively.

For western Yamal, the predicted value of the linear trend of air temperature increase is 0.2 °C/10 years based on the processing of long-term observations at the Marre-Sale weather station. The results of solving the problem show an increase in the average annual temperature in sands from -5.5...-6.9 °C to -3.5...-4.7 °C, in loams from -6.1...-7.0 °C to -4.0...-4.9 °C. The projected rate of temperature change is 0.4 °C/10 years (Fig. 2).

Figure 2. Temperature changes in saline frozen rocks of western Yamal.

1 and 2 are the results of solving the problem with and without taking into account the transformation of the phase and chemical composition of pore solutions, respectively.

In the lower reaches of the Ob, the predicted value of the linear trend of air temperature increase is 0.4 °C/10 years based on the processing of long-term observations at the Salekhard weather station. The results of mathematical modeling show a gradual increase in temperature in sands from -2.4...-3.6 °C to -1.5...-2.5 °C, in loams from -2.7...-3.7 °C to -1.6...-2.6 °C. The value of the linear trend of increasing ground temperature is 0.2 °C/10 years (Fig. 3).

Figure 3. Temperature change of saline frozen rocks of the lower Ob river.

1 and 2 are the results of solving the problem with and without taking into account the transformation of the phase and chemical composition of pore solutions, respectively.

For all the studied regions and types of sections, there is a difference between the results of modeling with constant and changing heat exchange characteristics of rocks. The difference between iterations gradually increases both with increasing mineralization of pore solutions and with the forecast period.

The simulated trends in rock temperature for northwestern and western Yamal are quite close. For northwestern Yamal, the smallest difference between iterations of 0.1—0.2 ° C is typical for the first 5-10 years of the forecast at a mineralization of 35 g/l. By 2050, the difference reaches 0.3—0.4 ° C at a mineralization of 150 g/l. For Western Yamal, the simulation results differ by 0.2—0.3 °C in the first years of the forecast with minimal mineralization. By the middle of the century, the difference reaches 0.4—0.5 ° C with the mineralization of a pore solution of 150 g/l. In loams, the discrepancies between independent temperature calculations turned out to be 0.2—0.3 °C greater than in sands.

For the lower Ob river, the greatest difference between iterations is noted. The minimum difference is 0.2—0.3 ° C at a mineralization of 35 g/l, the maximum is 0.5—0.6 ° C at a mineralization of 150 g/l. The average discrepancy between the calculation results is 0.3 °C for sands and 0.5 °C for loams.

Discussion

The combined use of thermodynamic and mathematical modeling methods revealed a significant difference between the calculation results. Failure to take into account the salinity of frozen rocks and their properties for the period of the geocryological forecast leads to a calculation error of up to 20%. The values of the average annual temperature of frozen rocks obtained with constant heat exchange characteristics turn out to be underestimated compared with the results of the task with time-varying characteristics. The difference between the results is not constant, at the initial stage of forecasting (the first 5-10 years) it is 8-10%. In the future, the difference between calculations increases, reaching a maximum in 25-30 years. With an increase in the mineralization of the solution from 35 to 150 g/l, the difference between calculations increases by 15-20%.

The curves of changes in the predicted temperature for northwestern and western Yamal are close to equidistant. This makes it possible to use correction coefficients, if necessary, to reduce the possible error that occurs when solving a thermal problem without taking into account changes in the water-ion composition of pore solutions. However, this approach is justified only for short-term forecasts, since after 10-15 years this error increases significantly. For the lower Ob river, the use of coefficients is not possible — the nature of the temperature change in rocks is more complex. Thus, in order to reduce the uncertainty of forecast estimates of the temperature of frozen rocks at the regional level with climate change and economic development of the territory, the numerical solution of the heat transfer problem, taking into account changes in the phase and chemical composition of pore solutions, seems preferable to introducing various coefficients into the calculation scheme.

Conclusion

The widespread distribution of saline permafrost soils in Yamal and their use as foundations of buildings and structures place additional requirements on the prediction of the thermal state of frozen rocks. The use in predictive estimates of the results of observations, experiments and laboratory studies aimed at clarifying the physical nature of processes in saline rocks, taking into account the transformation of the water-ion composition of pore solutions, can increase the accuracy of the results obtained. At the same time, there is uncertainty as to how the results of a physically more reasonable modeling of temperature fields may differ from the results of solving the "classical" thermal problem with constant heat exchange characteristics for the forecast period. The work is aimed at a quantitative study of the influence of salinity on the forecast estimates of the temperature regime of frozen rocks for three regions of Yamal.

One of the limitations of this technique is the semi-automatic calculation mode and the formulation of the heat transfer problem exclusively in a conductive form. The semi-automatic operation mode is due to the fact that two independent programs are used to solve the problem, which are necessary to assess changes in the water-ion composition of salt solutions and to find temperature fields. In the future, this restriction can be removed by embedding the "Freezbrine" module as a subroutine for an application for solving thermal problems and fully automatic determination of the thermophysical properties of rocks at each time step with appropriate recalculation of thermal fields.

The physical picture of freezing—thawing of rocks saturated with salt solutions is much more complex than unsalted rocks. The basic equations of the boundary value problem in the presence of phase transitions and chemical transformations should contain differential equations of thermal conductivity, moisture conductivity, salt transfer, and a complex interrelated type of initial and boundary conditions. The improvement of existing or the development of new applications, taking into account the system of these equations, will allow a more reasonable approach to assessing the future thermal state of saline frozen rocks.

Thanks

The authors thank PhD V. A. Kondakov and PhD A. B. Osokin for the opportunity to use stock materials of NTF Krios LLC on the characteristics of saline rocks of the Yamal Peninsula.