Региональная модель динамики льда. Часть 1. Описание модели, постановка численных экспериментов и современная динамика потока в окрестностях станции Конен
О. Рыбак1,4, Ф. Хёбрехтс1,2, Ф. Паттэн3, Д. Штайнхаге1
1Alfred-Wegener-Institut für Polar- und Meeresforschung, Bremerhaven, Germany; 2Vrije Universiteit Brüssel, Brüssels, Belgium; 3Université Libre de Bruxelles, Brussels, Belgium; 4Сочинский научно-исследовательский центр РАН
Статья поступила в редакцию 23 июня 2006 г. Предтавлена членом редколлегии В.Н. Голубевым
Рассмотрены методологические аспекты реконструкции нестационарной динамики потока льда в области 600x400 км на Земле Королевы Мод в Антарктиде с помощью региональной математической модели.
Введение
В последнюю четверть века математические модели стали важнейшим инструментом изучения макродинамики ледниковых щитов и связей их состояния с компонентами климатической системы Земли. Модели позволяют изучать как современные щиты, так и реконструировать исчезнувшие, а также прогнозировать состояние ледниковых щитов в будущем. Традиционно комплексные трехмерные модели ледниковых щитов [10, 20 и др.] строились на основе приближения уравнений движения мелкого льда SIA (Shallow Ice Approximation) [9]. В рамках SIA-моде-лей считается, что деформации льда обусловлены только касательными напряжениями в вертикальной плоскости. SIA-модели накладывают ограничения на соотношение между характерной толщиной щита и горизонтальным разрешением модели, приблизительно 1:10. Это ограничение не существенно для большей части Антарктиды, однако в таком случае нельзя с надлежащей точностью описать динамику
ледниковых потоков в окрестностях ледоразделов, линии налегания, в областях со сложной топографией ложа и др. Ограничение на пространственное разрешение не позволяет ассимилировать в модели данные полевых наблюдений, получаемые с большим пространственным разрешением. Собственно, уже для описания динамики шельфового льда и переходной зоны между шельфовым льдом и ледниковым щитом приходится отходить от 81А-моделей. Для описания в таких моделях полного тензора напряжений требуются вычислительные ресурсы на порядки большие, чем это необходимо при решении задач в рамках 81А-моделей. Альтернативой служит применение аппроксимации неполного второго порядка [2]. Физический смысл ее заключается в пренебрежении сопротивлением меняющемуся градиенту напряжения сдвига в направлении движения для вертикально проинтегрированного баланса сил. В [17] показано, что этот член на два порядка меньше по маг-нитуде, чем другие касательные и нормальные на-
- з-
Материалы гляциологических исследований, вып. 102
0
пряжения даже для тех районов, где происходит смена режима течения, например, в непосредственной близости от ледоразделов.
Чтобы оптимизировать расчеты деформаций льда и скоростей потока в модели, в которой используется более высокая, чем в SIA, степень аппроксимации, целесообразно разбить задачу на две части: 1) имитация динамики ледникового щита на комплексной SIA-модели с относительно грубым пространственным разрешением и 2) решение уравнений движения более высокого порядка аппроксимации в ограниченной области (встроенном домене) с большим пространственным разрешением. Этот подход был реализован при моделировании и прогнозе регионального климата. Модели такого типа в англоязычной литературе носят название «nested models», а мы будем называть их встроенными. Основная задача, которая должна быть решена при построении встроенной модели, заключается в организации обмена данными между двумя суб-моделями: большой (комплексная SIA-модель, БМ) и малой (модель встроенного домена с большим пространственным разрешением, ММ). Существует два пути организации обмена данными [4]: с динамической обратной связью между большой и малой моделями и без нее. В первом случае рассчитанные в малой модели поля транслируются обратно в большую, во втором нет. Первый подход чреват разного рода вычислительны-
40 80 120 160 200 240 280
Рис. 1. Положение домена 600x400 км, соответствующего области интегрирования региональной модели (выделен прямоугольником) на карте Антарктиды. Маркировка границ дана в соответствии с нумерацией узлов пространственной сетки с шагом 20 км. Южный полюс расположен в точке с координатами (141, 141) Fig. 1. Location of the nested domain of 600x400 km corresponding to the region of integration of the regional model (shown as grey square) on the map of Antarctica. Boundaries are marked in accordance with numbering of the spatial grid with 20 km step. South Pole is located in the point (141, 141)
ми сложностями в силу потенциальной несовместимости описания физических процессов в малой и большой моделях. Так, динамическая обратная связь была причиной численной неустойчивости в [22]. Обе суб-модели использовали одни и те же системы уравнений, только во встроенном домене проводилась десятикратная пространственная детализация. В большой модели генерировались краевые условия, которые интерполировались на границе встроенного домена. Малая модель возвращала в большую рассчитанную топографию поверхности щита в совпадающих узлах сетки. По мнению авторов [22], целесообразнее использовать подход без обратной связи, что было реализовано в модели Гренландского щита [7], где во встроенном домене решались уравнения более высокого порядка точности по сравнению с SIA, но без обратной связи с внешним доменом.
В настоящей работе мы рассматриваем встроенную модель, в которой во внутреннем домене применяется как пространственная детализация, так и более высокий порядок аппроксимации по сравнению с SIA. Глобальная динамика Антарктического ледникового щита имитируется трехмерной комплексной моделью (БМ), которая детально рассмотрена в [11, 12]. Модель состоит из взаимосвязанных компонентов, описывающих эволюцию собственно ледникового щита и шельфового льда, взаимодействие щита и подстилающих пород (гляциоизостатическое приспособление земной коры), климатический форсинг и баланс массы щита. Пространственное разрешение модели составляет 20 км по горизонтали и 30 слоев по вертикали. В малой модели использован неполный второй порядок аппроксимации [18]. Ее пространственное разрешение равно 2,5 км по горизонтали и 100 слоев по вертикали.
Моделирование эволюции Антарктического щита — одно из направлений программы EPICA (European Project for Ice Coring in Antarctica), в рамках которой в 1996—2005 гг. были пробурены две скважины в Восточной Антарктиде — на франко-итальянской станции Купол Конкордия (Dome C) и немецкой станции Конен (Kohnen) на Земле Королевы Мод (рис. 1). В течение нескольких лет в этом районе проводились полевые исследования, включающие закладку шурфов и бурение мелких скважин, радиоэхозондирование и т.д. В результате были построены карты топографии ложа и скорости аккумуляции с разрешением 2,5 км, что и определило пространственное разрешение малой модели. Эти данные были ассимилированы в модели в ходе численных экспериментов.
Станция Конен расположена на слабо выраженном гребне, и на этом же гребне находится станция Купол Фуджи, приблизительно в 1200 км выше по потоку (рис. 2а) в области со сложной топографией ложа (рис. 2б). Интерпретация результатов изучения ледяного керна требует детальных сведений о скоростях потока и деформациях льда в окрестности станции, что и послужило целью моделирования. Резуль-
9
4
ЭН
dt
■■ -V(vH) + M
(1)
ЭТ dt
k— Эz
- vVT + Ф
(2)
1 _Э_ рс^ дz
где Н — толщина льда, ? — время, М — баланс массы, V — трехмерный вектор скорости потока, V — вектор горизонтальной скорости, осредненный по глубине, Т — температура льда, Ф — приток тепла за счет внутреннего трения. Константы в уравнениях (1) и (2) и в последующих приведены в табл. 1. Уравнение сохранения количества движения может быть записано в форме:
О*
р— = Vo-dt
pg ,
(3)
где о — симметричный тензор напряжений Коши с компонентами о.у. После некоторых преобразований [18] уравнение (3) в скалярной форме имеет вид:
Рис. 2. Топография поверхности (а) и ложа (б) во встроенном домене (м над ур. моря) Fig. 2. Surface (a) and bedrock (b) topography in the nested domain (m a.s.l.)
таты, которые обсуждаются ниже, и методы, с помощью которых они получены, были кратко представлены ранее в виде докладов на конференциях [13, 14].
Описание модели
Баланс сил и уравнения движения льда. Модель динамики ледникового щита строится в предположении, что лед представляет собой несжимаемую вязкую неньютоновскую жидкость. Базовая стандартная система уравнений модели включает уравнения сохранения массы, энергии и момента количества движения. В общепринятой форме уравнения сохранения массы и энергии имеют вид:
Эо хх + Эо ху
Эх Эу Эz
Эо ух , Эо уу , Эо yz
Эх Эу Эz
^zz = pg .
эz
= 0
= 0
(4)
Интегрируя последнее уравнение из (4) от поверхности щита s(x,y) вниз до подстилающей поверхности и подставляя результат в первые два уравнения, получим:
д п — (2т Эх
хх + туу
_э_
Эу
(2тУУ + Тхх
Этху ,Э 'х = pg Эs
Эу эz Эх
Эт ух = pg Эs
Эу эz Эу
(5)
где
- 3 5j X k °kk
и Sj
символ Кронекера. Уравнения (5) описывают баланс сил в аппроксимации неполного второго порядка. Пренебрегая первыми и вторыми членами в (5), получим уравнения 8!Л-аппроксимации:
Таблица 1
Константы, входящие в уравнения (1)—(13)
Константа
Символ Значение Единица измерения
Р 910 кгм-3
Pw 1000 кгм-3
^ 1,8 ■ 10-10 Н-3год-1м8
T0 273,15 K
ср g 2009 Джкг1^-1
9,81
k 6,6210-7 Джм-1^год-1
m 1
Плотность льда Плотность воды Коэффициент трения Температура плавления льда Теплоемкость льда Гравитационная постоянная Коэффициент теплопроводности льда Коэффициент усиления
5 -
Материалы гляциологических исследований, вып. 102
дтхг дз —— = рё— дг дх
дт
уг _,
д5
д ' = Pg Т" дг ду
(6)
Закон Глена [16] связывает напряжения и скорости деформации: Ту = 2уёц , где
у = I А(Т *)-1/п е/1"«)/« ,
ее — второй инвариант тензора скоростей деформации, А(Т ) — функция температуры, исправленной на давление, п=3. Закон Глена определяет термомеханическое сопряжение в моделях динамики ледниковых щит
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.