В последние годы возросло количество исследований, направленных на оценку баланса массы ледниковых покровов и подледного таяния. Вода в ледниках и под ними занимает сравнительно небольшой объем относительно льда (1–2 %), но именно она влияет на динамику изменения климата. Согласно современным представлениям, антарктический ледниковый покров скрывает около пятисот подледных водоемов. Изучение жизненного цикла подледниковых водоемов, включающего их формирование и последующее развитие вплоть до прорыва, либо замерзания, становится особенно актуальным в свете наблюдаемых глобальных климатических изменений и повышения уровня моря. Кроме того, исследование субгляциальных гидрологических процессов дает понимание не только состояния и динамики ледника, но и формирования современного подледного рельефа. Наряду с данными полевых исследований все большее значение приобретают результаты численного моделирования.
Компьютерное моделирование является одним из способов определения вероятных местоположений подледниковых озер. Чаще всего в качестве математической модели используются либо уравнение Навье-Стокса, либо его приближения [1, 2]. Они позволяют описывать движение ледника как текучей непрерывной среды, но при этом не учитывают наличие фазового перехода лед–вода. Для учета последнего можно использовать математическую модель на основе краевой задачи Стефана с подвижной границей – задачи определения температурного распределения в каждой фазе вещества и положения межфазной границы.
В данной работе представлена численная модель донного таяния антарктического ледника на основе одномерной мультифазной задачи Стефана, которая описывает появление, эволюцию и вырождение под- и надледниковой жидкой фазы. Модель учитывает ряд конкретных условий протекания процесса. В частности, учитывалась аккумуляция, поверхностный температурный режим, вертикальное движение ледника.
Трехмерное представление ледяного щита Антарктиды было построено на основе топографических данных из базы международного проекта Bedmap2 [3]. Данные содержат три топографические характеристики: рельеф ложа ледника, рельеф поверхности ледника и толщина льда. Характеристики определены на равномерной сетке с шагом 1 км в декартовой системе координат и измеряются в метрах над уровнем моря. В каждом узле данной сетки независимо друг от друга решались одномерные трехфазные задачи Стефана вдоль вертикального направления, где первая и третья фазы вещества – подледниковая и надледниковая вода соответственно, а вторая фаза – ледник.
В данной работе не учитывались горизонтальное распространение тепла, перераспределение талой воды, взаимодействие с океаном, растекание ледника вследствие его аморфности.
Математическая модель. Одномерная задача Стефана формулируется следующим образом: для j-той фазы вещества с распределением тепла, занимающей объем пространства , необходимо определить в каждый момент времени и координаты границ и фазы в каждый момент времени.
Ось X направлена от ложа ледника вверх. Фазы вещества расположены вдоль оси X в следующем порядке: подледниковая вода ледник надледниковая вода.
Для каждой фазы вещества записывается краевая задача вида
(1)
где – удельная теплоемкость, – плотность, – коэффициент теплопроводности воды.
Для подледниковой воды – геотермический поток от ложа горных пород, – линейная зависимость температуры плавления льда от толщины льда [4].
Для ледника – температура плавления льда.
Для надледниковой воды – температурный режим на поверхности ледника.
Движение границ фаз определяется условием Стефана [5]:
(2)
(3)где – удельная теплота плавления льда, – скорость движения границы без учета проседания ледника. Одна и та же масса льда занимает больше объема, чем та же самая масса воды, и это влияет на движение поверхности ледника. Скорректированная скорость дается выражением:
(4)
где ζ – поправка, связанная с аккумуляцией осадков на поверхности ледника.
Уравнение для границы s3:
(5)
Перед дискретизацией задачи была проведена процедура обезразмеривания, заключающаяся в переходе от размерных координат , температуре и времени к безразмерным координатам , и следующими преобразованиями:.. В качестве характерного размера выбрана начальная толщина льда в данной точке континента. В качестве характерной температуры используется температура плавления льда Характерное время выбрано таким, чтобы максимально упростить вид уравнения теплопроводности:
Существенную сложность для численного решения представляет подвижность границ. Для решения этой проблемы был применен метод выпрямления фронта [6]. Он заключается в переходе от нестационарной области к стационарной путем применения преобразования:
(6)
Дискретизация задачи. После перехода к новым координатам задачу можно решать любыми сеточными методами. В данной работе для дискретизации задачи применялась неявная конечно-разностная схема на равномерной и неравномерной сетке. Введение неоднородной сетки, сгущающейся у фазовых границ, позволяет уменьшить количество узлов, что существенно снижает время расчетов. Для неоднородной сетки использовалась неявная разностная схема с первым порядком точности по времени и по координате, сводящая задачу Стефана к системе линейных уравнений с трехдиагональной матрицей.
Вырождение и появление фаз вещества. Математическая формулировка задач Стефана никак не учитывает возможность появления/исчезновения фаз. Границы разных фаз вещества не знают друг о друге, и при совпадении их координат в условии Стефана появляется деление на ноль. Поэтому, данная ситуация требует отдельной обработки на уровне численного решения.
В случае трехфазной задачи Стефана реализована возможность исчезновения и появления водных фаз. Вырождение фазы льда не предусмотрено, поскольку характерная толщина ледника значительно больше уровня талой воды под или над ним.
Механизм вырождения фаз следующий: если в ходе эволюции координат оказывается, что , то граница устанавливается на координату , и заменяются краевые условия. Если исчезла верхняя фаза, то ее краевое условие с границы приписывается границе верхней кромки ледника . Аналогично в случае исчезновения нижней водной фазы.
Механизм появления фазы был реализован на основе закона сохранения энергии. При формальном вычислении распределения температуры никак не обрабатывается состояние переохлажденной воды или перегретого льда. Поэтому, возможна такая ситуация, когда часть узлов фазы льда, находящихся у границы фазы, имеет температуру выше температуры плавления льда. В этом случае происходит вычисление количества «лишней» энергии путем численного интегрирования методом трапеций, и эта энергия направляется на зарождение новой фазы.
Условия и результаты моделирования. На поверхности задан синусоидальный режим изменения температуры:
(7)
где– время, – среднегодовая температура, – годовая амплитуда температуры. Две последние величины получены из [7].
Данные по аккумуляции осадков в Антарктиде и по геотермическому потоку взяты из [8] и [9] соответственно.
Топографические характеристики ледника (рельеф ложа, рельеф поверхности и толщина льда) получены из базы Bedmap2. Моделирование производилось на разреженной сетке с шагом 10 км. Также были отброшены точки с ненулевым расстоянием между горными породами и нижней кромкой ледника. Под данный критерий попали все шельфовые ледники и подледниковое озеро Восток. В итоге, расчет был проведен для 118 453 точек континента.
Параметры моделирования: шаг по времени месяца, лет, на фазу льда 5000 узлов расчетной сетки, на водные фазы по 500 узлов. В качестве начального распределения температуры в толще ледника выбрана функция [4].
Моделирование было выполнено в пакете MATLAB. Расчет проводился в паралельном асинхронном режиме для всех узлов сетки Антарктиды.
Расчеты проводились на ресурсах РЦ ВЦ СПбГУ, использовался вычислительный узел с процессором Intel Xeon E5-2690 v4 с 14 ядрами (28 потоков) и 250 GB RAM. Расчет занял 12 ч 55 мин. После оптимизации кода удалось сократить это время до 1 ч 42 мин.
Результаты моделирования представлены на рис. 1. Наблюдается скорость донного таяния в пределах от
0 до 15 мм/год (рис. 1а) с наибольшей интенсивностью в районе Западной Антарктиды и Антарктического полуострова. Наблюдается проседание ледника, связанное с ненулевой скоростью донного таяния и различными плотностями воды и льда (рис. 1б).
Обсуждение результатов. В двух районах Антарктиды (Антарктический полуостров и Западная Антарктида) наблюдается интенсивное донное таяние (порядка 10–
15 мм/год). В обоих случаях это объяснимо крайне высоким тепловым потоком порядка 120–150 мВт/м² при норме для остального континента около 50 мВт/м². На полуострове толщина льда не превышает 500 метров, что также могло повлиять на результат.
В реальности, как правило, вода не скапливается под ледником, а движется по системе каналов, которая в свою очередь также меняется со временем. Полноценный учет перераспределения воды представляет из себя отдельную непростую задачу.
В пересчете на массу талой воды в единицу времени для всей Антарктиды скорость таяния составила
29 Гт/год (усредненное значение для континента
2,5 мм/год) Для сравнения, в [10] приведены скорости донного таяния для различных моделей, где значения достигают 4000 Гт/год. Столь существенное расхождение в первую очередь связано с шельфовыми ледниками, которые не учитываются в данной работе. Например, в [11] приводятся данные по средней скорости таяния на шельфовом леднике Туэйтса, которая составляет 10–30 м/год. Такое большое число обусловлено откалыванием айсбергов от шельфовых ледников и обменом тепла с океаном.
Выводы. Получена средняя скорость таяния для ледяного щита Антарктиды в 29 Гт/год. Для сравнения, в работе [2] приводится значение 65 Гт/год, причем авторы не рассматривают обмен тепла с океаном и откалывание шельфов, но в их модель включена динимика льда как аморфного целого. В работе [10] сравниваются результаты приблизительно десятка моделей, где скорости донного таяния варьируются от 100 до 4000 Гт/год. В эти модели включены шельфовые ледники с учетом взаимодействия льда с океаном и откалыванием айсбергов, чем объясняется расхождение с настоящей работой на два порядка.