CC BY
149
И. А. Судаков, Л. П. Бобылёв, С. А. Береснев
МОДЕЛИРОВАНИЕ ТЕРМИЧЕСКОГО РЕЖИМА ВЕЧНОЙ МЕРЗЛОТЫ ПРИ СОВРЕМЕННЫХ ИЗМЕНЕНИЯХ КЛИМАТА
Введение
В настоящее время наблюдаются глобальные климатические изменения, связанные с повышением среднегодовой температуры воздуха на планете (так называемое глобальное потепление). Прогнозируется, что в ближайшем будущем эти изменения будут только усиливаться и приведут к негативным последствиям для всей климатической системы. Особо уязвимым климатическим элементом природной среды является вечная мерзлота, которая занимает большую часть Северного полушария Земли.
Частичное или полное исчезновение вечной мерзлоты под воздействием глобальных климатических изменений может привести к серьезным экономическим и политическим проблемам в северных регионах планеты [1, 2]. Таяние вечной мерзлоты в условиях глобального потепления обусловливает дополнительную эмиссию парниковых газов, которые до этого времени были законсервированы в мерзлотной толще. Изменение состояния вечной мерзлоты должно повлиять на функционирование экосистем, которые расположены на ее территории. Кроме того, понимание физических процессов в мерзлоте и ее эволюции необходимо для исследования подобных природных образований на других планетах.
Исследованиям эволюции и устойчивости вечной мерзлоты в условиях изменяющегося климата посвящены многочисленные публикации (например, [3-6]). К ним в первую очередь относятся исследования, основанные на применении глобальных климатических моделей. Недавние оценки термического режима вечной мерзлоты по климатическим характеристикам, определенным на основе данных моделей, представлены в [5, 7]. Предлагаемые модели достаточно хорошо описывают динамику температурного поля в верхних слоях вечной мерзлоты (так называемые деятельные слои глубиной в несколько метров) с учетом сезонных климатических колебаний в течение небольших временных периодов (от 1 года до 5 лет). С другой стороны, актуальна разработка оценочных моделей эволюции термического состояния вечной мерзлоты, которые могли бы прогнозировать развитие процессов в толще мерзлых грунтов на длительных временных промежутках (до 100 лет).
Целью настоящего исследования является разработка достаточно простой физико-математической модели, описывающей изменение термического режима вечной мерзлоты, а также расчет и анализ эволюции состояния многолетнемерзлых грунтов п-ова Ямал в XXI в. на основе входных данных глобальной климатической модели ЕСНАМ 5 [8].
© И. А. Судаков, Л. П. Бобылёв, С. А. Береснев, 2011
Физико-математическая модель для определения термического режима вечной мерзлоты и схема ее численной реализации
Результаты исследований свидетельствуют о том, что многолетнемерзлые грунты представляют собой сложные гетерогенные многокомпонентные и многофазные физико-химические системы, формирующиеся в значительных интервалах вариаций термодинамических параметров их состояния. Взаимодействие компонентов и фаз в этих средах осуществляется посредством поверхностных сил различных видов, часто через сложные многослойные граничные переходные зоны, содержащие, в общем случае, жидкие, жидкообразные, твердообразные, твердые слои. Все перечисленные компоненты находятся в физико-химическом и механическом взаимодействии, интенсивность и формы которого зависят от температуры.
При моделировании многолетнемерзлых грунтов часто предполагают, что они состоят всего из двух фаз — твердой и жидкой, на границе которых происходит фазовый переход «промерзание—оттаивание» [9, 10]. Термический режим многолетнемерзлых грунтов описывается моделями на основе типа задачи Стефана. Такой подход к моделированию термического режима многолетнемерзлых грунтов оправдан, когда климатические условия меняются периодически (при сезонных колебаниях температуры воздуха) и когда имеются эффективные численные методы решения типа задачи Стефана.
С другой стороны, модель многолетнемерзлых грунтов можно представлять как модель гетерогенной, многокомпонентной и многофазной системы, описывающую основные физико-химические процессы между компонентами и фазами. Такая модель полезна для детального описания свойств вечной мерзлоты, но она мало пригодна и вычислительно не эффективна для описания процессов, происходящих на макроскопическом уровне за длительные промежутки времени.
Предложим физическую модель, описывающую термический режим многолетнемерзлых грунтов на достаточно длительном временном промежутке (от года до столетия), которая базируется на следующих основных допущениях:
1. Предположим, что многолетнемерзлый грунт — однофазная система с некоторыми эффективными теплофизическими параметрами. В связи с поставленными задачами рассматриваем изменение термического состояния многолетнемерзлых грунтов на протяжении длительного периода (до 100 лет). Для этого в качестве граничного условия выбираем среднегодовую температуру воздуха на поверхности. Такой подход позволяет избежать излишней детализации в описании процессов, но позволяет сохранить основную климатологическую цель исследований.
2. В модели рассматривается изменение температурного профиля в однофазной системе в течение длительного времени (больше одного года). Выберем критерий для оценки состояния вечной мерзлоты: если температура на линии температурного профиля выше 273 К (0°С), то происходят процессы оттаивания, ниже этого значения — процессы промерзания.
3. Температурное поле многолетнемерзлых грунтов нестационарно. Поэтому для описания термического режима используется нестационарное уравнение теплопроводности. Теплофизические характеристики грунтов полагаются постоянными во всей мерзлой толще (из-за небольшого числа геокриологических стационаров практически отсутствуют данные по теплофизическим свойствам мерзлых грунтов для многих географических районов).
4. Толща многолетнемерзлых грунтов является вполне определенной для конкрет-
ного географического региона. Из анализа геотермальных температур в толщах различных геологических платформ следует, что на достаточно больших глубинах эти температуры в основном определяются величиной геотермального потока из мантии. По этой причине будем задавать температуру на нижней границе слоя мерзлоты постоянной, определяемой либо экстраполяцией известного хода геотермальной температуры по мощности залегания мерзлых толщ, либо из малочисленных экспериментальных данных. Проблема определения температуры на нижней границе подробнее рассмотрена в [11].
Постановка задачи
На поверхности мерзлого грунта как однофазной системы по определенному закону происходит изменение температуры воздуха Та^г. Оно, благодаря кондуктивному теплообмену, влияет на многолетнемерзлый грунт, изменяя внутри него температурный режим. Теплообмен описывается нестационарным уравнением теплопроводности. Нижняя граница грунта, расположенная на существенном удалении от поверхности, имеет известную температуру Тд. Эта «глубинная» температура принимается постоянной в течение всего периода времени. Математическая формулировка задачи в одномерной постановке следующая:
дТ (х,€) д2Т (х,Ь)
уравнение теплопроводности: С --------= Л——^— и > 0, 0 < х < х ), (1)
от дхх
Здесь Л — коэффициент теплопроводности, С — объемная теплоемкость грунта.
Уравнение теплопроводности, сформулированное в такой постановке, является одномерным параболическим уравнением 2-го порядка с граничными условиями первого рода для температуры на поверхности.
Обратимся к методу аналитического решения подобной краевой задачи, основанной на теории возмущений. Для этого необходимо ввести некоторую неизвестную функцию ю(х,Ь), представляющую собой отклонение от известной функции и(х,Ь) и учитывающую нестационарность процесса:
Задача имеет решение, если известен вид функции и(х,Ь). В [9] рассмотрен пример с функцией и(х, £) для случая, когда граничные условия как на верхней, так и на нижней границе слоя постоянны. В рассматриваемой задаче граничное условие на верхней границе задается некоторой функцией времени ф(Ь), а вид и(х,Ь) в общем случае не известен. Это затрудняет прямое применение метода [9], который, по всей видимости, эффективен лишь для определенного класса функций в аналогичных проблемах.
В монографии [9] представлены аналитические решения краевой задачи для случая полуограниченной пластины с известным начальным распределением температуры. Для пластины конечной толщины аналитических решений такого рода, по-видимому, не имеется.
На сегодняшний день существуют эффективные методики численного решения подобных краевых задач. В частности, известны подходы, использующие различные формулировки метода конечных объемов.
Нестационарное уравнение теплопроводности решается методом Патанкара [12], который основан на комбинации различных численных методов. Для реализации метода
граничные условия: Tair = <p(0,t), Tg = ф(х* ,t) = const, начальные условия: T(0, 0) = Tair, T(x*, 0) = Tg.
(2)
T (x,t) = U (x,t) +v(x,t).
(3)
Патанкара была создана вычислительная программа, которая состоит из двух частей: неизменяемой и адаптируемой. Неизменяемая часть программы была заимствована из [13]. Преимущество неизменяемой части программы состоит в том, что вычислительный алгоритм и программный код могут быть использованы для решения любого уравнения переноса (теплопроводности, диффузии). Адаптируемая часть программы создана самостоятельно в соответствии с задачами настоящего исследования. Программный код реализовывался на языке программирования FORTRAN 90.
Численные эксперименты и обсуждение результатов
Предложенная ранее модель использовалась для расчета состояний термического режима и прогноза изменения температуры вечной мерзлоты п-ова Ямал. Район был выбран по трем причинам: во-первых, п-ов Ямал является типичной бореальной экосистемой, для которой характерны практически все типы мерзлотных грунтов и геокриологических условий; во-вторых, экспериментальные данные, полученные с геокриостационаров п-ова Ямал, часто используются для сравнения с расчетными данными различных моделей; в-третьих, в настоящее время идет активное освоение данной территории, что предполагает повышенное внимание к различным модельным оценкам характеристик вечной мерзлоты.
Для численной реализации модели потребовались данные по температуре поверхности для различных временных масштабов, по температуре на нижней границе мерзлоты, значения теплофизических параметров многолетнемерзлых грунтов.
Данные по температуре поверхности (принимаемой равной температуре воздуха) для XX в. были оценены по климатической модели ЕСНАМ 5 Метеорологического института им. Макса Планка (Германия), а для XXI в. — по сценарию МГЭИК А2 для выбросов парниковых газов из той же модели [8]. Этот сценарий предполагает негативное развитие событий и увеличение температуры воздуха за столетие на 4-5 градусов.
В файле данных по п-ову Ямал из модели ЕСНАМ 5 представлены сведения о температуре поверхности, пространственной сетке модели, приведены долготы и широты центров модельных ячеек. Из них были выбраны те, центры которых попадают в область, ограниченную широтами 67°N-75°N и долготами 65°Е-72°Е, причем это точки для суши по маске океан/суша (0/1) модели ЕСНАМ. Непосредственно для расчетов выбиралась ячейка с координатами центра: широта 69°94^ и долгота 67°50'Е. Этот район расположен достаточно близко к геокриологическому стационару Марре-Сале и соответствует морфологии многолетнемерзлых грунтов с типом «полигональная тундра». На основе этих данных был построен ход среднегодовой поверхностной температуры начиная с 80-х годов XX в. и до конца XXI в. Наблюдение за вечной мерзлотой п-ова Ямал началось только в 80-е годы прошлого века, поэтому для проверки модели ход поверхностной температуры рассматривается с этого же времени.
На протяжении XXI в. действительно наблюдается тенденция к резкому повышению температуры поверхности, что может негативно сказаться на развитии всей бореаль-ной экосистемы п-ов Ямал. Для удобства численного моделирования среднегодовой ход температуры поверхности был аппроксимирован квадратичным полиномом, подобранным программой 0^ш7.5 и именно в таком виде использован в расчетах.
Как уже отмечалось ранее, нет однозначных рекомендаций по выбору температуры на нижней границе исследуемого слоя вечной мерзлоты. Анализируя данные [14] и экстраполируя известный ход геотермальной температуры по толще залегания мерзлоты,
приходим к выводу, что нижнюю границу слоя залегания целесообразно оценить в 15 м, а значение температуры на ней как Т = 266К (—7°С).
Значения теплофизических параметров для полигональной тундры задавались согласно ГОСТу 25100-82—«Грунты. Классификация» и по данным [3], а именно: коэффициент теплопроводности Л = 1,7 Вт/(м • К), объемная теплоемкость грунта С = 2,4 Дж/(м3 • К).
Прежде всего, была рассчитана зависимость температуры мерзлых грунтов на глубинах 2 м и 10 м от времени за период 1980-2000 гг. Оказалось, что термический режим на глубине в 2 м зависит от изменения среднегодовой температуры поверхности. На глубине в 10 м температура слоя монотонно увеличивалась за весь анализируемый период времени.
Кроме того, были приведены расчетные зависимости температуры мерзлых грунтов на глубинах 2 м и 10 м в период 2001-2100 гг. На глубине 2 м наблюдается тенденция к повышению температуры в течение всего исследованного промежутка времени. С 2089 г. температура слоя становится равной 273 К (0°С), что свидетельствует о начале активного таянии, а к 2096 г. слой вечной мерзлоты выше 2 м полностью исчезает. Отметим, что к 2093 г. на глубинах 3-4 м температура приблизится к отметке 273 К (0°С), а к 2100 г. слой вечной мерзлоты выше 4 м, возможно, полностью исчезнет. Для глубины 10 м также характерно достаточно сильное повышение температуры вплоть до конца XXI в. К 2100 г. температура слоя может достичь величины 270,8 К (-2, 2°С), однако мерзлый грунт на этой глубине сохранится.
Как отмечалось ранее, для расчетов был выбран район вблизи геокриологического стационара Марре-Сале. Это единственный на Западном Ямале стационар, где за последние 30 лет были получены данные для температуры многолетнемерзлых грунтов на различных глубинах. Из них следует, что в целом наблюдается тенденция к повышению температуры верхних мерзлотных слоев, а также к увеличению глубины сезонного протаивания.
На рис. 1 показан временной ход наблюдаемой и расчетной температуры вечной мерзлоты за 1980-2000 гг. на глубине 10 м. Видно, что в начальный момент времени экспериментальные и расчетные температуры были достаточно близки к значениям 266К (-7°С) и
265К (-8°С) соответственно. Для расчетной температуры наблюдается ее монотонное увеличение. При этом расчетная температура не выходит за пределы известного наблюдательного интервала 264К (—9°С)-268К(-5°С) для вечномерзлых грунтов ис-
Т, К 270
269
268
267
266
265
264
263
262
.А.а-
А А
Аа
8 10 12 14
А Эксперимент
16 18 20 22 Годы (1980-2000)
■ Модельный расчет
Рис. 1. Сравнение экспериментальных данных [6] (треугольники — исходные данные, кривая 2 — их линейный тренд) и расчетов (кривая 1) для района Марре-Сале п-ова Ямал на глубине 10 м в 1980—2000 гг. Жирные линии сверху и снизу — экспериментальный интервал наблюдаемых температур по [15].
следуемого района [15]. По всей видимости, расхождение в 1,5 К может служить косвенной оценкой точности температурных предсказаний предлагаемой модели. Путем аппроксимации экспериментальных данных получен линейный тренд температуры, который в течение исследуемого периода увеличивается. Темп повышения экспериментальной температуры составил 0,018 К/год, темп повышения расчетной температуры — 0,010 К/год. Одинаковые поведения линейного тренда и расчетной температуры и близкие темпы повышения температур свидетельствуют о хорошей согласованности пред-
На рис. 2 представлено сравнение прогнозируемого изменения среднегодовой поверхностной температуры и температуры вечной мерзлоты на глубине 10 м в XXI в.
В [6] представлена обобщенная зависимость изменений поверхностной температуры и температуры мерзлых грунтов на глубине 1020 м в рамках регионального прогноза для арктических и субарктических районов России (врезка на рис. 2). Здесь же представлена аналогичная расчетная зависимость по предлагаемой модели для изучаемого района. Видно, что эти зависимости качественно согласуются: при повышении среднегодовой температуры поверхности существенно изменяется и термический режим вечной мерзлоты. Оцененный темп повышения температуры поверхности для исследуемой территории составляет 0,09 К/год, а темп повышения температуры многолетнемерзлых грунтов — 0,06 К/год. Отметим, что в [6] при прогнозировании использовался другой сценарий — А1 (МГЭИК). Темп повышения температуры поверхности для арктических и субарктических регионов составил 0,021 град/год, а темп повышения температуры вечной мерзлоты — 0,017 град/год.
В исследованиях термического состояния вечной мерзлоты для данного района в XXI в. важное значение имеет прогноз для глубины протаивания мерзлого грунта. Протаявшим считался слой вечной мерзлоты, для которого расчетная температура достигла значения в 273К (0°С). К 2100 г. мерзлота на глубине, возможно, 4 м полностью исчезнет. Это предельное значение глубины протаивания вечной мерзлоты на изучаемой территории, полученное в расчетах по используемой модели. Не исключено, что глубина сезонного протаивания может быть и больше, но предлагаемая модель не учитывает данных сезонных амплитуд. Отметим, что, по данным [4, 7] для субарктических регионов Западной Сибири, прогнозные оценки для глубины протаивания лежат в интервале 2-6 м.
ложенной модели с экспериментальными данными. Т, К
, „ Годы
-------1 - тренд поверхностной температуры
-------2- тренд температуры вечной мерзлоты на 10 м
Рис. 2. Прогнозируемое изменение среднегодовой температуры поверхности (кривая 1) и температуры вечной мерзлоты на глубине 10 м (кривая 2) для исследуемого района в 2001—2100 гг. На врезке: то же самое, но для перепада температур (в градусах) на глубине 10—20 м в арктических и субарктических районах России по [6].
На данном этапе предложенная модель адекватна задачам выполненного исследования. Однако возможно и дальнейшее ее развитие. В настоящей модели максимальное число итераций в вычислительной программе составляет 109. Этого вполне достаточно для прогнозных оценок в несколько веков при временном шаге в 1 год. При уменьшении шага (например, до полугода с целью учета сезонных изменений) потребуется большее число итераций, что допустимо к изменению в данной программе. Временной шаг в 1 год выбран в соответствии с климатологической задачей исследования, где не учитывались сезонные колебания температуры. Возможны уменьшение временного шага до полугода и выделение двух периодов: холодного и теплого. При этом температура поверхности должна аппроксимироваться соответствующей временной функцией.
В данной версии программы шаг по глубине в 1 м близок к оптимуму как с точки зрения характерных масштабов задачи, так и с точки зрения вычислительного процесса. Его увеличение или уменьшение не даст выигрыша в точности вычислений.
Рассматриваемый вариант модели не учитывает детального описания физических процессов на поверхности вечной мерзлоты (теплообмен со снежным покровом, органическим слоем и т.д.). Однако нет принципиальных ограничений модели для учета данных процессов при обоснованном уменьшении временного шага.
Открытым остается вопрос о температурных условиях на нижней границе мерзлоты. Такая проблема типична для многих аналогичных моделей. Однозначных рекомендаций по выбору температуры на нижней границе расчетной области на сегодняшний день, к сожалению, не существует.
Заключение
В данной статье предлагается простая оценочная модель описания термического режима вечной мерзлоты для больших промежутков времени, основанная на решении краевой задачи на базе нестационарного уравнения теплопроводности. Численная реализация предложенной модели осуществляется простым и достаточно эффективным методом контрольных объемов (методом Патанкара). Основываясь на той же вычислительной схеме Патанкара, предлагается также метод оценки фильтрационного вклада в потоки эмиссии парниковых газов из вечной мерзлоты при изменении ее термического режима.
На основании предложенной модели был проанализирован термический режим вечной мерзлоты п-ова Ямал в условиях современных климатических изменений. Изучение эволюции термического состояния вечной мерзлоты в XXI в. проводилось с использованием сценария выбросов парниковых газов А2 (МГЭИК). Расчетные данные по температуре вечной мерзлоты в конце XX в. в целом согласуются с экспериментальными данными, полученными на стационаре Марре-Сале. Установлено, что к концу XXI в. ожидается значительный рост температуры многолетнемерзлых грунтов, который может привести к их частичному оттаиванию. Слой вечной мерзлоты выше 2 м полностью исчезнет, а процессы протаивания распространятся на глубину 3-4 м. Следует отметить, что прогнозная оценка термического режима вечной мерзлоты и перестройки геокриологических условий существенно зависит от выбранного сценария ожидаемых в XXI в. климатических изменений. Кроме того, описаны дальнейшие пути развития предлагаемой модели. Возможно совершенствование как самой физикоматематической модели, так и вычислительной программы, ее реализующей.
1. The Physical Science Basis, Contribution of Working Group I to the Fourth Assessment Report of the IPCC. Cambridge University Press, 2007. 390 p.
2. Paepe R., Melnikov V. P. Permafrost Response on Economic Development, Environmental Security and Natural Resources. Kluwer Academic Publishers, 2001. 664 p.
3. Zhuang Q., Romanovsky V. E. McGuire A. D. Incorporation of a permafrost model into a large-scale ecosystem model: Evaluation of temporal and spatial scaling issues in simulating soil thermal dynamics // J. Geophys. Res. 2001. Vol. 106. P. 649-670.
4. Клименко В. В., Ершов Э. Д. и др. Изменения климата и динамика толщ многолетнемерзлых пород на Северо-Западе России в ближайшие 300 лет // Криосфера Земли. 2007. №3. С. 3-13.
5. Мачульская Е.Е., Лыкосов В. Н. Моделирование термодинамической реакции вечной мерзлоты на сезонные и межгодовые вариации атмосферных параметров // Известия АН. Физика атмосферы и океана. 2002. №1. С. 20-33.
6. Мельников В. П., Павлов А. В. Изменение климата на Севере России и геокриологические последствия // Возможности предотвращения изменения климата и его негативных последствий: материалы Совета-семинара при Президенте РАН / Отв. ред. Ю. А. Израэль. М.: Наука, 2006. C. 341-352.
7. Павлова Т. В., Катцов В. М., Надежина Е. Д. и др. Расчет эволюции криосферы в XX и XXI веках с использованием глобальных климатических моделей нового поколения // Криосфера Земли. 2007. №2. С. 3-14.
8. Official Website of Max Planck Institute for Meteorology. ECHAM — Global Climatic Model. [Electronic resource] / [Germany, Hamburg], [2008]. Mode of access: http://www.mpimet.mpg. de/en/wissenschaft/modelle/echam.html (дата обращения: 28.05.2010).
9. Карташов Э. М. Аналитические методы в теории теплопроводности твердых тел. М.: Высшая школа, 2001. 550 с.
10. Комаров И. А. Термодинамика и тепломассообмен в дисперсных мерзлых породах. М.: Научный мир, 2003. 608 с.
11. Ткаченко Е. И. Исследование тепловых задач для сред с изменением агрегатного состояния на основе новой формулировки нижнего граничного условия: Автореф. дис. . . . канд. тех. наук / Е. И. Ткаченко. Тюмень: Тюм. ун-т, 2006. 18 с.
12. Патанкар С. В. Численные методы решение задач теплообмена и динамики жидкости: Пер. с англ. М.: Энергоатомиздат, 1984. 152 с.
13. Патанкар С. В. Численное решение задач теплопроводности и конвективного теплообмена при течении в каналах / Пер. с англ. М.: Изд-во МЭИ, 2003. 312 с.
14. Дучков А. Д., Соколова Л. С., Балобаев В. Т. и др. Тепловой поток и геотемпературное поле Сибири // Геол. и геоф. 1997. №11. С. 1716-1729.
15. Природные условия Байдарацкой губы / Под ред. В. А. Совершаева и др. М.: ГЕОС, 1997. 432 с.
Статья поступила в редакцию 1 июня 2010 г.
