Scientific journal
Modern high technologies
ISSN 1812-7320
"Перечень" ВАК
ИФ РИНЦ = 0,940

MODELING OF DEGRADATION OF PERMAFROST IN THE PROCESS OF CLIMATE WARMING IN CENTRAL YAKUTIA FOR THE COMING 300 YEARS

Ivanov V.A. 1 Rozhin I.I. 1
1 Institute of Oil and Gas Problems SB RAS
The unprecedentedly rapid climate warming in the northern regions is having a significant impact on many natural ecosystems. Probably one of its effects will be large-scale thawing of the permafrost layer. To predict the depth of thawing of the frozen ground, the 1D Stefan problem is solved with the help of a model developed for paleoreconstruction of the thermal regime of frozen rocks. The model used takes into account the dependences of the coefficients of thermal conductivity and heat capacity of the porous medium on the content of water and ice in it, the dependences of the thermal conductivity of water and the matrix of the porous medium on temperature, and also the dependence of the melting point of ice on pressure and salt content. Well thermometry, core research data, and data on the modern climate of Central Yakutia, including monthly average air temperatures, radiation balance of the earth’s surface, and heat transfer coefficient on the surface taking into account snow cover, were used as input data for numerical calculation. Accounting for annual fluctuations in air temperature and radiation balance allows us to more accurately predict the thermal state of shallow ground. Two scenarios of climate warming are considered – an increase in average annual temperatures at a rate of 2 °C/100 years and at a rate of 4 °C/100 years. At the output of the computational experiment, we obtained pictures of an unsteady thermal field in the soil, which show the gradual thawing of permafrost soils over time. The results will allow us to clarify our ideas about the spatio-temporal patterns of how permafrost react to the ongoing process of climate warming, which is necessary to develop adaptation and mitigation measures.
permafrost
Stefan problem
mathematical modeling
global warming

В последнее столетие наблюдается постепенный рост среднегодовых температур воздуха во всем мире, регулярно фиксируются температурные рекорды и климатические аномалии [1]. Глобальное потепление набрало особенно высокую скорость в последние десятилетия и в силу этого способно вызвать множество последствий, которые коренным образом могут поменять мир, в котором мы привыкли жить. Это и затопление прибрежных территорий Мировым океаном вследствие таяния ледников, и наступление пустынь в аридной зоне, и изменения в биосфере планеты, и многое другое. Одним из подобных последствий может оказаться масштабная деградация многолетней мерзлоты в северных широтах, которая усугубляется тем фактом, что в Арктике потепление климата выражено еще сильнее. Сложно оценить экологические, экономические и социальные риски, связанные с резким разрушением криолитозоны. Некоторыми из имеющих место процессов для детального изучения этих рисков могут служить возможная потеря устойчивости грунтов под инженерными сооружениями или эрозия почв на берегах водоемов [2].

Многолетнемерзлые породы (ММП) распространены в арктической и субарктической климатических зонах, а также в горных районах. Их толщина может достигать 500 м и более. В настоящее время температура мерзлых пород находится в пределах нескольких градусов ниже нуля, поэтому их оттаивание может начаться в скором времени.

В Центральной Якутии в настоящее время имеет место слабый отклик глубины сезонно-талого слоя и в целом теплового состояния ММП на наблюдаемое повышение температуры воздуха [3]. В основном это объясняется синхронным уменьшением высоты снежного покрова и, следовательно, уменьшением его теплоизолирующей способности. Однако на окраинах криолитозоны в других регионах происходит деградация ММП [4]. Например, на северо-востоке европейской части России и на севере Западной Сибири. В Печоро-Уральском регионе по результатам мониторинга отмечаются рост температуры мерзлых пород, рост глубины таликов и осадка грунтов [5]. На территории Уренгойского месторождения в Западной Сибири идет устойчивое повышение температуры верхних горизонтов пород, а в некоторых местах произошел отрыв кровли ММП и ее продвижение вниз идет со скоростью примерно 30 м/100 лет [6]. В приведенных работах указывается, что, за исключением Центральной Якутии, измеренное повышение температуры ММП отстает от повышения температуры воздуха в 2–7 раз.

В работе [7] проведено моделирование температурного режима грунтов за последние 3,5 млн лет на территории Уренгойского месторождения. В разработанной очень подробной модели учтены зависимости коэффициента теплопроводности породы от фазового состава порового флюида, его зависимость от температуры, зависимость температуры таяния льда от давления и солености, а также возможное наличие незамерзшей воды в мерзлых грунтах. Моделирование показало, что за рассмотренный период произошло как минимум 10 циклов похолодания, которые сопровождались появлением слоя ММП и зоны стабильности газовых гидратов. Выявлено, что температурное поле находится в сильной зависимости от физико-химических свойств грунтов, в особенности от коэффициента теплопроводности. В работе [8] построен прогноз изменений температуры грунтов в Центральной Якутии при двух сценариях потепления – умеренный со скоростью 0,04 °C/год и агрессивный со скоростью 0,08 °C/год. Согласно использованной модели даже при агрессивном сценарии температура приповерхностного грунта поднимется выше 0 °C только через 40 лет. При этом глубина сезонного оттаивания при тренде 0,04 °C/год практически не меняется, а при тренде 0,08 °C/год – увеличивается с 2,0 до 2,2 м. Также следует отметить, что в этой модели потепление климата на 1 °C вызывало повышение температуры грунтов на 0,33 °C.

В настоящей работе ставится цель провести анализ темпов деградации слоя ММП при возможном повышении атмосферной температуры с помощью математической модели, предложенной в работе [7].

Материалы и методы исследования

Согласно модели из [7], для массива горных пород решается одномерная задача Стефана по глубине

missing image file (1)

где Cv(T) – объемная теплоемкость пористой среды, зависящая от её состава, температуры и глубины залегания и в которую включена теплота фазового перехода «лед – вода»; T – температура горных пород; t – переменная времени; z – глубина; λ(T) – коэффициент теплопроводности пористой среды, зависящий от ее состава, температуры и глубины залегания; A – внутреннее тепловыделение горных пород, вызванное естественной радиоактивностью.

Объемная теплоемкость пористой среды вычисляется как взвешенная средняя арифметическая по ее составу:

missing image file (2)

где TL – температура плавления льда; Cv,m – объемная теплоемкость скелета пористой среды, которая находится из исследований керна; ϕ – пористость среды, которая зависит от глубины и также определяется из исследований керна; Cv,w – объемная теплоемкость воды, равная 4,187 МДж/(м3· °С); W(T) – доля незамерзшей поровой воды, которая зависит от температуры и типа горных пород; Cv,i – объемная теплоемкость льда, равная 1,926 МДж/(м3· °С); ρw – плотность воды, равная 1000 кг/м3; L – удельная теплота плавления льда, равная 334 кДж/кг.

Температура плавления льда зависит от гидростатического давления, содержания солей в поровой воде натрий- и калийхлоридного типа

TL = 0 °C – 0,073P – 0,064Cs, (3)

где TL – температура замерзания воды в °С, P – давление в МПа, Cs – концентрация солей в г/л.

Коэффициент теплопроводности пористой среды вычисляется как средняя геометрическая по ее составу:

missing image file (4)

где λm – коэффициент теплопроводности скелета пористой среды, который зависит от температуры; λw – коэффициент теплопроводности воды, тоже зависящий от температуры; λi –коэффициент теплопроводности льда, который принят равным 2,26 Вт/(м∙ °С).

Коэффициент теплопроводности скелета пористой среды находится по следующей зависимости от температуры:

missing image file (5)

где λm,20 – коэффициент теплопроводности скелета при 20 °С, который определяется из исследований керна; αT – температурный коэффициент.

Коэффициент теплопроводности воды также зависит от температуры:

missing image file (6)

На верхней границе задается условие 3-го рода, учитывающее конвективный и радиационный теплообмен

missing image file (7)

гд α – суммарный коэффициент теплообмена между атмосферным воздухом и грунтом, Tair – температура воздуха, qsurf – радиационный баланс земной поверхности.

На нижней границе массива горных пород задается условие 2-го рода с геотермическим потоком тепла

missing image file (8)

Для единственности решения системы уравнений (1)–(8) математическую модель необходимо дополнить начальным условием

missing image file (9)

Поставленная задача решается методом конечных элементов с использованием пакета прикладных программ для технических вычислений.

Объектом исследования является неоднородный массив горных пород глубиной 3000 м, из которых в условиях Центральной Якутии примерно на 500 верхних метрах расположены ММП, а на оставшуюся часть приходятся талые горные породы. Для сокращения краевых эффектов используется неравномерная сетка разбиения, шаг которой постепенно уменьшается при приближении к верхней границе. Начальное распределение температуры массива горных пород по данным термометрии представлено на рис. 1. Видно, что температурное распределение в большей части мерзлой зоны практически равномерное с близкими к нулю отрицательными температурами. Это может свидетельствовать о том, что в настоящий момент ММП находятся на пороге начала оттаивания.

missing image file

Рис. 1. Начальное распределение температуры горных пород по глубине

В расчетах содержание незамерзшей воды в мерзлых песчаных породах, а также содержание солей в поровой воде принимаются равными нулю. Также для условий Центральной Якутии приняты усредненные значения A = 1,5·10–6 Вт/м3, αT = 0,002 °С–1, qgeo = 0,057 Вт/м2. Физические параметры многослойного массива горных пород представлены в табл. 1. Данные по среднемесячной температуре воздуха и по радиационному балансу земной поверхности, приведенные в табл. 2, взяты из работы [9]. Значения суммарного коэффициента теплообмена (табл. 2) рассчитаны с учетом снежного покрова, который препятствует поступлению холода и, соответственно, понижению температуры верхних горизонтов ММП в зимний период.

Таблица 1

Характеристики слоев массива горных пород

Интервал, м

Плотность скелета пористой среды, кг/м3

Объемная теплоемкость скелета пористой среды, МДж/( °С∙м3)

Коэффициент теплопроводности скелета пористой среды, Вт/(м∙ °С)

Пористость, д.е.

0–86

2000

2,13

2,24

0,21

86–488

2000

2,27

2,12

0,21

488–539

2120

2,37

1,77

0,19

539–980

2300

2,18

2,39

0,13

980–1831

2350

2,19

2,72

0,12

1831–2561

2380

2,19

2,86

0,12

2561–3000

2330

2,19

2,51

0,12

Для проверки адекватности численной модели проведено сравнение среднемесячных температур грунта на глубине 3,2 м с наблюдениями на площадке с естественным покровом на стационаре «Туймаада» ИМЗ СО РАН [10] (рис. 2). В этом расчете для верхнего граничного условия (7) использованы архивные данные по погоде с сайта rp5.ru, а также данные по радиационному балансу луговой площадки из табл. 2.

missing image file

Рис. 2. Начальное распределение температуры горных пород по глубине

Расчетные и наблюдаемые на стационаре кривые совпадают в достаточной мере.

Таблица 2

Среднемесячные значения данных атмосферы на начальный момент расчета

Месяц

I

II

III

IV

V

VI

VII

VIII

IX

X

XI

XII

Средняя месячная температура воздуха, °С

–37,5

–32,7

–21,4

–5

7,3

16,6

19,8

15,4

5,8

–7,4

–27,6

–38

Радиационный баланс, Вт/м2

–9,65

–4,63

–5,02

45,14

113,04

126,16

119,6

79,48

35,49

0

–14,27

–9,65

Коэффициент теплообмена между атмосферным воздухом и грунтом, Вт/(м2· °С)

0,39

0,40

0,43

9,28

27,87

27,46

26,19

25,30

24,38

8,89

0,54

0,42

Рассматриваются две скорости будущего потепления климата: 1) рост температуры воздуха с темпом 4 °С/100 лет; 2) рост с темпом 2 °С/100 лет. Первая скорость соответствует повышению температуры поверхности Земли при отсутствии дополнительной деятельности по уменьшению выбросов CO2, согласно докладу МГЭИК [1]. Вторая скорость соответствует целевому уровню увеличения средней температуры планеты, согласованному на конференции по климату в Париже в 2015 г. [11]. В течение каждого года температура воздуха задается циклически изменяющейся функцией, аппроксимирующей его среднемесячные значения. Так как моделирование проводится на срок, составляющий 300 лет, то повышение температуры воздуха в ХХI в. линейно экстраполируется на два последующих столетия.

Расчеты выполнены в коммерческом пакете прикладных программ для решения задач технических вычислений.

Результаты исследования и их обсуждение

На рис. 3 представлены распределения температуры по глубине в верхней толще мерзлых пород для двух рассмотренных сценариев потепления климата. После 100-го года при тренде увеличения температуры воздуха 2 °С/100 лет (рис. 3, a) и после 50-го года при тренде 4 °С/100 лет (рис. 3, b) заметны повышения температур пород, находящихся ниже сезонно-талого слоя глубиной до 3 м. Потепление климата провоцирует деградацию многолетней мерзлоты и увеличение глубины сезонного оттаивания, что может приводить к развитию термокарста при наличии сильнольдистых ММП, термоэрозии и термоденудации. Видно, что профили температуры в талой области имеют практически линейный характер. Следовательно, при потеплении климата тепловая энергия, подводимая на земную поверхность, в основном расходуется на плавление порового льда в верхних горизон- тах ММП.

missing image file

Рис. 3. Динамика распределения температуры по глубине в осенний период при темпе роста температуры воздуха: (a) 2 °С/100 лет (b) 4 °С/100 лет

missing image file

Рис. 4. (a) Прогнозная динамика изменения глубины кровли ММП в зависимости от сценария потепления климата. (b) Прогнозная динамика изменения глубины сезонного оттаивания ММП в зависимости от сценария потепления климата. Пунктирные линии – глубина кровли ММП после отрыва от поверхности

На рис. 4, a, представлены прогнозные глубины кровли ММП на период до 300 лет при действии темпов потепления климата на 2 °С/100 лет и 4 °С/100 лет. При этом в течение первых 60 и 120 лет, соответственно темпам потепления, растепление грунтов находится в пределах сезонно-талого слоя глубиной 3 м (рис. 4, b). Через 50 лет глубина сезонно-талого слоя увеличится на 30 см (13 %) при темпе потепления климата 2 °С/100 лет и на 50 см (23 %) при темпе потепления 4 °С/100 лет. Этот результат близок к значениям Росгидромета [12], согласно которому в районе г. Мирный к середине XXI в. ожидается увеличение глубины сезонного оттаивания на 20 % – с 1,95 м (2000–2010 гг.) до 2,34 м (2040–2050 гг.).

В дальнейшем происходит отрыв кровли ММП и ниже сезонно-талого слоя начинает формироваться таликовая зона, которая со временем увеличивается по линейному закону (рис. 4, a). Эта зона растет со скоростью примерно 12 м/100 лет при темпе потепления климата 2 °С/100 лет и примерно 15 м/100 лет при темпе потепления климата 4 °С/100 лет.

Погрешности результатов данного моделирования в наибольшей степени зависят от упрощающих допущений модели и качества исходных данных. Во-первых, сильное влияние на термическое состояние ММП оказывает растительность (возможно уменьшение мощности сезонно-талого слоя до 50 %) [13], вариации которой в данной работе не учитываются. Во-вторых, выходные данные МОЦ несут в себе определенную неточность, которая может быть довольно существенной для территории Якутии [14]. И в-третьих, теплофизические свойства грунтов могут значительно меняться в зависимости от местности.

Заключение

Проведен вычислительный эксперимент по определению масштабов оттаивания многолетней мерзлоты при двух сценариях потепления климата для условий Центральной Якутии. Использована численная модель, разработанная для палеореконструкции теплового режима горных пород. При темпе роста температуры воздуха 4 °С/100 лет, что соответствует сценарию с отсутствием дополнительной деятельности по уменьшению выбросов CO2, через 60 лет начинается деградация многолетней мерзлоты со скоростью примерно 15 м/100 лет. В случае если темп потепления климата удержится на отметке 2 °С/100 лет, то деградация многолетней мерзлоты начнется через 120 лет, а ее скорость составит примерно 12 м/100 лет. Эти прогнозные оценки получены при изменениях температуры в верхних горизонтах ММП, обусловленных только повышением температуры воздуха. Дополнительные исследования необходимы для оценки влияния многолетней динамики снежного покрова, солнечной радиации, альбедо земной поверхности, скорости ветра и других региональных характеристик климата на термическое состояние криолитозоны.