НАУКИ О ЗЕМЛЕ
УДК 910.1 551.8
ЭЙЛЕРОВ И ЛАГРАНЖЕВ МЕТОДЫ РАСЧЕТА ВОЗРАСТА ЛЬДА В ПРОСТОЙ ЧИСЛЕННОЙ МОДЕЛИ ЛЕДНИКОВОГО ЩИТА
© 2008 г. О.О. Рыбак1
Для интерпретации косвенных данных (изотопного состава, концентрации примесей и др.), получаемых из ледяных кернов в результате глубокого бурения ледниковых щитов, необходимо их точное датирование. Независимые методы датирования реализованы в математических моделях растекания дна. В статье анализируются эйлеров и лагранжев методы датирования льда. В этих целях проведена серия численных экспериментов с простой моделью ледникового щита. Показано, что в целом точность лагранжева метода выше по сравнению с эйлеровым, особенно в нижних слоях ледникового щита.
Ключевые слова: ледниковый щит, математическая модель, возраст льда, уравнение адвекции.
ВВЕДЕНИЕ
Результаты исследования криосферы Земли год от года играют все более важную роль в понимании механизмов изменения климата планеты и связанных с этим вариаций уровня Мирового океана, атмосферной циркуляции, эволюции экосистем и т.д. Ледниковые щиты Антарктиды и Гренландии представляют собой своеобразные природные архивы, которые содержат косвенную информацию об изменениях климата на протяжении последних сотен тысяч лет. В последние три десятилетия были выполнены многочисленные национальные и международные программы глубокого бурения антарктического и гренландского ледниковых щитов. Цель этих проектов - не только восстановить историю колебаний климата в прошлом, но и оценить перспективы его изменения в будущем. Интерпретация данных о физико-химическом составе льда, получаемых в ходе бурения, требует максимально точного датирования. Непосредственный подсчет годовых слоев льда возможен лишь до определенного временного горизонта, характерный возраст которого 10-20 тыс. лет, ниже которого происходит стирание физических и химических свойств между слоями отдельных лет. Для определения возраста льда ниже этого горизонта прибегают к различного рода синхрониза-циям однородных событий в различных кернах (например, синхронизация вулканических извер-
1 Сочинский научно-исследовательский центр Российской академии наук, 354000, Краснодарский край, Сочи, ул. Театральная, 8; тел. (862-2) 62-59-86, e-mail: snic@sochi.ru.
жений) или к датировке с помощью изотопов, скорость накопления которых условно считается постоянной (10Ве), и т.д. Независимое датирование производится методами математического моделирования растекания льда.
Задача датирования решается аналитически лишь для ограниченных условий, в частности, когда горизонтальной адвекцией можно пренебречь, или же она может быть описана в очень упрощенном виде. Другие сложности связаны с пространственными вариациями баланса массы на поверхности ледникового щита, изменениями толщины и температуры льда вдоль потока и особенностями его реологии [1]. Существенное влияние на формирование поля возраста льда оказывают процессы базального таяния и базального скольжения. Базальное таяние начинается в областях ледникового щита, где температура на границе раздела льда и подстилающих пород под действием геотермического тепла достигает точки плавления. В этих же областях масса льда в дополнение к составляющей скорости, обусловленной внутренними деформациями, приобретает составляющую за счет скольжения относительно подстилающих пород (базальное скольжение).
В трехмерных моделях динамики ледникового щита задача расчета возраста льда сводится к решению уравнения адвекции [2]. В рамках такого подхода возраст рассматривается как консервативное свойство льда, которое переносится вместе с потоком. Задача расчета адвекции может решаться в рамках одного из двух базовых методов - эйлерова или лагранжева. Лагранжев метод в некоторых физических приложениях получил название РагТМе-т-СеП [3]. Каждый из этих
методов имеет достоинства и недостатки. В частности, реализация эйлерова метода в конечно-разностных схемах требует введения дополнительной «искусственной диффузии» для обеспечения устойчивости решения [2]. Основные сложности лагранжева метода связаны с процедурами интерполяции, которые решаются в рамках так называемого полулагранжева метода [4], который можно рассматривать как адаптированный эйлеров метод [5].
Все три метода (эйлеров, лагранжев и полула-гранжев) с недавних пор используются в гляциологических приложениях. Так, в [6] эйлеров и лагранжев методы применены в исследовании распределения возраста с глубиной в одной из глубоких скважин в гренландском ледниковом щите. Показано, в частности, что эйлеров метод недостаточно точен в нижних слоях ледникового щита. Лагранжев метод позволяет решить проблему дополнительной диффузии и провести более точное датирование в нижних слоях. Преимущество полулагранжева метода над эйлеровым заключается в том, что он позволяет обойти ограничение на величину шага по времени [4, 6, 7]. С помощью полулагранжева метода были расчитаны характеристики изотопного баланса Гренландского щита [7]. К недостаткам полулагранжева метода следует отнести потерю индивидуальности виртуальными частицами (трассерами) на каждом временном шаге интегрирования, и, следовательно, возникновение дополнительных сложностей при построении линий тока [4, 7].
Построению хронологий (датированию) конкретных ледяных кернов посвящены сотни статей. Цель настоящей работы - сравнение возможностей эйлерова и лагранжева методов датирования льда в простой численной модели ледникового щита. Простота модели позволяет сделать общие заключения о влиянии тех или иных факторов на распределение возраста льда в теле щита. Поскольку одной из исследуемых характеристик является искривление траекторий частиц льда под влиянием базального таяния и базаль-ного скольжения, в работе используется не полу-лагранжев метод, а метод Рагййе-т-СеП.
1. ЛАГРАНЖЕВО И ЭЙЛЕРОВО ОПИСАНИЕ ДВИЖЕНИЯ ЧАСТИЦ ЛЬДА
Рассмотрим декартову систему координат (х, I), где х - горизонтальная ось, 7 - вертикальная ось, и точку на плоскости, координаты которой в начальный момент времени были (х0, г0).
Пусть в момент времени t эта точка переместилась в точку (х, I) Назовем лагранжевым такой способ описания, когда перемещение является функцией начальных координат точки:
xt Х0 = UL (Х0,70,t)
7 - 70 = WL (Х0 ,70> О '
(1.1)
Для заданной пары начальных координат (х0, г0) линия тока материальной точки задается в параметрической форме хг = х0 + иь(^, = 10 + WL(t).
Будем называть эйлеровым такой способ описания, когда перемещение частицы льда выражается как функция конечных координат, т.е. в каждый новый момент времени рассматривается новая материальная точка:
■ Х0 = ^Е (Х£, , О 7 - 70 = ще(Х, 7,О '
(1.2)
Обозначим через А момент времени t, когда частица льда была отложена на поверхности ледникового щита. Естественно считать, что эта величина со временем не меняется. Таким образом, время отложения (возраст льда) есть консервативная характеристика, которую можно считать пассивным трассером и описать его эволюцию уравнением
dA( х, I, £) dt
= 0.
(1.3)
Лагранжев подход к решению (1.3) будет заключаться в интегрировании уравнений (1.1) с
начальным условием |0,10 = (Н + И)Хо % ^ отражающим то обстоятельство, что в начальный момент времени частица льда (трассер) находилась на поверхности щита.
В эйлеровой форме уравнение (1.3) запишется как
дА
+ цУА = 0,
(1.4)
где и - вектор скорости потока льда с компонентами и и w. Граничное условие на поверхности щита А(х, zH+h, ^ = t, где Н - толщина щита, к - абсолютная высота подстилающей поверхности. Граничное условие у подстилающей поверхности [2]
ЭА ЭА — = В--
д1
, ЭА' и 11=» Т"
дх
(1.5)
где В - скорость базального таяния (м/год в ледовом эквиваленте). Чтобы обеспечить устойчи-
вость численного решения (1.3) в конечных разностях, в [2, 7] было предложено добавить в правую часть (1.3) дополнительный диффузионный член. Таким образом, его можно переписать в виде
дА дА дА п д2А --+ и— = w— = иА —2-•
д? дх дz дz
(1.6)
На нижней границе
дА
дz
=0. В [8] оценена
2=к
оптимальная величина DA = 5 • 10-8 м2/с. Показано также, что поскольку характерное время адвекции на два порядка меньше, чем характерное время диффузии, то введение дополнительного диффузионного члена в (1.6) не приведет к искажениям поля возраста льда, кроме возраста в придонной области.
2. ТЕСТИРОВАНИЕ МЕТОДОВ
2.1. Задание профиля щита и компонент скорости
Для того чтобы провести сравнение методов, воспользуемся аналитически полученным полем скоростей внутри ледникового щита. В простейшем случае профиль ледникового щита, известный как профиль Ная-Вьялова (Nye-Vialov [9]), на ровной подстилающей поверхности задается уравнением
Н (х) = Но
т
п+1 / п
п /2(п+1)
(2.1)
где L - расстояние от ледораздела до края щита, Н0 - толщина щита на ледоразделе,
Но =
20 М
1/2( п+1)
vPg/
,1/2( п+1)
г1/2
(2.2)
М = 0,1 м/год - характерная величина баланса массы щита, F = 10-16 Па/год, р = 910 кг/м3 -плотность льда, п = 3, g = 9,181 м/с2 - ускорение свободного падения. Ось х направлена вдоль потока.
В первом приближении всеми напряжениями, определяющими деформацию и скорость течения льда, за исключением тангенциального напряжения в вертикальной плоскости, можно пренебречь [9]. В этом случае в изотермическом леднике горизонтальная компонента скорости будет определяться только пространственным
градиентом толщины льда%
и( X, 2) = ■
2£_ п + 1
дН
- ( Н - 2)
дх
п +1
п-1
дН
дх
[ Н
+1
(2.3)
Зададим вертикальную скорость в простом виде как
w( X, 2) =--
Н
^ + В + и( х, Н) —
дх
- В, (2.4)
где О - скорость снегонакопления (аккумуляции) на поверхности щита (м/год в ледовом эквиваленте), и(х, Н) дН - составляющая вертикаль-дх
ной скорости, обусловленная наклоном оси х по отношению к горизонтальной плоскости. Простой вид уравнения (2.4) позволяет делать аналитические расчеты вертикальной скорости на ледоразделе.
2.2. Постановка экспериментов
Зададим размер щита как 2L = 2000 км. Расчеты будем вести на пространственной сетке 101 х 101 узлов с равным расстоянием между ними по горизонтали. В этом случае расстояние между узлами с
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.