автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Математическая модель динамики медленнотекущих объектов

005001526

па правах рукописи

СЗЗ^

Рыбак Олег Олегович

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ДИНАМИКИ МЕДЛЕННОТЕКУЩИХ ОБЪЕКТОВ (НА ПРИМЕРЕ ВНУТРЕННИХ ОБЛАСТЕЙ ЛЕДНИКОВЫХ

ЩИТОВ)

05.13.18 - Математическое моделирование, численные методы и комплексы программ

1 О НОЯ 2011

АВТОРЕФЕРАТ диссертации на соискание ученой степени доктора физико-математических наук

Краснодар - 2011

005001526

Работа выполнена в лаборатории математического моделирования Сочинского научно-исследовательского центра Российской академии наук

Научный консультант: доктор физико-математических наук,

профессор Семенчин Евгений Андреевич

Официальные оппоненты: академик РАН, доктор физико-

математических наук Глико Александр Олегович

доктор физико-математических наук, профессор Зотов Владимир Александрович

доктор физико-математических наук, профессор

Толпаев Владимир Александрович

Ведущая организация: Южный федеральный университет (г.

Ростов-на-Дону)

Защита состоится « &» декабря 2011 г. в час. на заседании совета по

защите докторских и кандидатских диссертаций Д 212.101.17 при Кубанском государственном университете по адресу: 350040, г. Краснодар, ул. Ставропольская, 149, ауд. 231.

С диссертацией можно ознакомиться в научной библиотеке Кубанского государственного университета

Автореферат разослан (.Л^у, е^^г^^и^ 2011

Ученый секретарь совета по защите

докторских и кандидатских диссертаций Д 212.101.17

кандидат физико-математических наук,

доцент

Барсукова В.Ю.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ Актуальность проблем, решаемых в диссертационной работе, обусловлена как фундаментальными, так и прикладными аспектами, возникающими при описании динамики медленнотекущих объектов, таких, в частности, как покровные, горные, шельфовые ледники, вулканические лавы, некоторые виды глин и так далее. Математическое моделирование является эффективным инструментом и во многих случаях единственным средством решения научных и прикладных задач исследования таких природных объектов.

В диссертационной работе математическая модель разрабатывается для описания течения льда во внутренней области ледникового щита. Выбор объекта исследования не случаен. В последнее время во всем мире растет интерес к всестороннему изучению полярных регионов. В их исследовании участвуют как национальные организации, так и международные исследовательские консорциумы. В последние три десятилетия был выполнен ряд международных программ глубокого бурения антарктического и гренландского ледниковых щитов, призванных восстановить историю изменений климата на протяжении последних нескольких сотен тысяч лет, понять механизмы, определяющие закономерности его вариаций, и оценить перспективы изменений в будущем. Интерпретация данных, получаемых в ходе бурения, требует применения методов математического моделирования и новых методов обработки данных численных экспериментов для их сопоставления с данными, получаемыми иными методами.

Помимо чисто научного интереса к исследованию полярных районов, необходимо принять во внимание то, что в настоящее время в мире нарастает тенденция борьбы за минеральные ресурсы, значительная часть которых находится именно в труднодоступных полярных областях. Необходимо также учитывать, что в ледниковых щитах Антарктиды и Гренландии аккумулированы основные запасы пресной воды на Земле, и, в условиях нарастающего ее дефицита в отдельных регионах мира, неизбежно будет подниматься вопрос о том, как в той или иной форме использовать их. С этой точки зрения важно обеспечить постоянное присутствие нашей страны в полярных регионах, в Антарктиде в частности. Перспективное использование еще не разведанных полезных ископаемых и водных ресурсов полярных регионов, какую форму бы оно ни приняло в будущем, требует тщательного анализа динамических характеристик континентального оледенения в том числе и методами математического моделирования. Таким образом, актуальность диссертационного исследования наряду с чисто научными, фундаментальными, обусловлена практическими аспектами освоения полярных регионов.

Объектом исследования работы является лед как медленнотекущая природная жидкость, предметом - процесс течение льда в ограниченной пространственной области ледникового щита.

Цель исследования - разработать эффективный численный метод решения трехмерной задачи расчета течения изотропной несжимаемой вязкой неньютоновской жидкости в ограниченной пространственной области с

большим пространственным разрешением. На его основе предложить математическую модель эволюции заданного района ледникового щита в меняющихся внешних условиях. Основные задачи исследования

1. На основе алгоритмов решения трехмерной задачи эволюции ледникового щита с упрощенным описанием течения льда и алгоритмов решения стационарной задачи течения льда с высокой степенью приближения уравнений динамики льда (базовые субмодели), разработать метод решения задачи эволюции течения льда в ограниченной пространственной области с большим разрешением, разработать метод объединения двух базовых субмоделей, алгоритмы и программные коды рабочего интерфейса между ними, методы постановки численных экспериментов для моделирования эволюции ледникового щита. Это позволяет:

- включить в численные эксперименты новейшие подробные данные наблюдений о топографии подстилающей поверхности ледникового шита и скоростях аккумуляции и новейшие данные об изменениях окружающей среды в прошлом;

- оценить параметры течения льда в настоящем и прошлом (скорость, направление, стабильность, характер деформаций);

- восстановить топографию поверхности ледникового щита и ложа, поля приземной температуры воздуха и скорости аккумуляции в прошлом;

2. Разработать методику пост-экспериментальной обработки данных численных экспериментов, в рамках которой:

- построить алгоритмы и программные коды для расчета траекторий движения частиц льда в потоке;

- проанализировать различия в эйлеровом и лагранжевом подходах к расчетам возраста и изотопного состава льда в модели ледникового щита;

- уточнить современные скорости аккумуляции в области с недостаточным количеством прямых измерений;

- рассчитать функцию динамического сжатия льда и ее трансформацию вдоль траекторий движения частиц льда;

- рассчитать неклиматическое смещение приземной температуры воздуха, рассчитанной по косвенным изотопным данным.

- уточнить соотношение между температурой воздуха и скоростью аккумуляции в прошлом.

3. Разработать метод модельной датировки льда, с помощью которого:

- построить хронологическую шкалу ледяного керна;

- оценить влияние величины потока геотермического тепла на возраст льда и неклиматическое смещение в нижних слоях ледникового щита.

Методы исследования. Для решения поставленных задач в работе использовались методы механики сплошной среды, реологии, теории пластичности и ползучести, теоретические основы построения математических моделей, численные методы решения дифференциальных уравнений в частных производных, а также вычислительные эксперименты с использованием программных средств.

Научная новизна

1. Построена и успешно использована в геофизических приложениях трехмерная модель течения льда, в которой в комплексную трехмерную термомеханическую субмодель Антарктического ледникового щита с упрощенным описанием динамики и разрешением 20 км по горизонтали была встроена модель, основанная на более высокой степени приближения уравнений течения льда с разрешением 2,5 км.

2. Разработана методика пост-экспериментальной обработки данных численных экспериментов. Эта методика позволила рассчитать возраст льда в ледяном керне, оценить неклиматическое смещение в рядах реконструированной приземной температуры воздуха, рассчитать функцию динамического сжатия.

3. Предложен метод реконструкции скорости палеоаккумуляции на основе модельных функций динамического сжатия и с использованием данных радиозондирования.

4. В ходе численных экспериментов были исследованы закономерности зависимости возраста придонного льда в Антарктиде от величины потока геотермического тепла.

5. Модельными методами были определены перспективные области для глубокого бурения Антарктического ледникового щита для получения льда, имеющего максимально большой возраст.

Практическая значимость. Был предложен метод, в соответствии с которым была построена эффективная математическая модель динамики медленнотекущей жидкости, для чего была разработана схема объединения двух типов субмоделей с различным уровнем приближения уравнений течения льда и разным пространственным разрешением. Подобный подход может использоваться в дальнейшем для изучения региональных особенностей динамики ледниковых объектов. Разработанная модель может быть использована для прогноза полей скорости и деформаций льда в районах континентального, в том числе и горного, оледенения при хозяйственном освоении этих территорий в условиях реализации того или иного сценария климатических изменений. Алгоритмы расчетов могут быть включены как составная часть в комплексную модель горных ледников и использоваться для прогноза объемов стока. После незначительной модификации модель может быть использована для описания динамики иных природных медленнотекущих жидкостей (например, вулканической лавы).

Достоверность и обоснованность научных результатов и выводов обеспечивается использованием апробированного математического аппарата, адекватности математических моделей и подтверждена данными наблюдений за течением ледников и физико-химического анализа ледяных кернов. На защиту выносятся следующие основные положения:

• Трехмерная математическая модель динамики медленнотекущей несжимаемой неньютоновской жидкости;

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

неньютоновской жидкости и различное пространственное разрешение (на примере природных ледниковых объектов);

• Метод обратного отслеживания для обработки данных численных экспериментов;

• Модельный метод построения хронологической шкалы льда в глубоких скважинах в ледниковых щитах;

• Модельный метод расчета неклиматического (топографического) смещения в изотопно-температурных рядах;

• Метод восстановления скоростей палеоаккумуляции, основанный на модельном анализе данных полевых исследований;

• Модельный метод исследования возраста льда в придонных частях ледниковых щитов в условиях неопределенности потока геотермического тепла.

Апробация работы

Основные результаты диссертации неоднократно обсуждались

• В секциях Генеральных ассамблей Европейского Союза наук о Земле (в Ницце, Франция, в 2002-2004 гг., в Вене, Австрия, в 2005-2008, и 2010-2011 гг.);

• на симпозиумах Международного гляциологического общества:

По физическим и механическим процессам во льду и их связи с моделированинем ледников и ледниковых щитов в Шамони-Мон-Блан, Франция, в 2002 г.;

По гляциологии Антарктиды в Милане, Италия, в 2003 г.; По радиогляциологии и ее приложениям в Мадриде, Испания, в 2008 г.;

• на открытых научных конференциях SCAR (Scientific Committee on Antarctic Research - Научный комитет по антарктическим исследованиям): «Антарктика и Южный океан в Глобальной системе» в Бремене, Германия, в

2004 г.;

«Полярные исследования - арктические и антарктические перспективы в Международный полярный год» в Санкт-Петербурге, Россия, в 2008 г.;

• на конференции Европейского научного фонда «Полярные регионы и четвертичный климат» в Аквафредда-ди-Маратеа, Италия, в 2005 г.;

• на Всероссийском симпозиуме по полярной гляциологии в Сочи, Россия, в

2005 г.;

• на Открытой научной конференции EPICA2008 «Четвертичный климат - от полюса до полюса» (в 2008 г. в Венеции, Италия);

• на рабочих совещаниях по программе EPICA (в 2004 г в Париже, Франция, в 2005 г. в Аквафредда-ди-Маратеа, Италия, в 2006 г. - в Иль-Чиокко, Италия);

• на рабочих совещаниях по программе NEEM (в 2009 и 2010 гг. в Копенгагене, Дания);

• в 2002-2011 гг. - на семинарах и рабочих встречах в группах гляциологии и палеоклиматологии Института им. Альфреда Вегенера, Бремерхафен, Германия, на географическом факультете университета Брюсселя (Vrije Universiteit Brüssel), в Лаборатории гляциологии и геофизики окружающей среды (LGGE) Универстета им. Жака Фурье, Гренобль, Франция.

Область исследования. Содержание диссертации соответствует паспорту специальности 05.13.18 «Математическое моделирование, численные методы и комплексы программ» (физико-математические науки) по следующим областям исследований: п. 1 «Разработка новых математических методов моделирования объектов и явлений», п. 3 «Разработка, обоснование и тестирование эффективных численых методов с применением современных компьютерных технологий», п. 4 «Реализация эффективных численных методов и алгоритмов в виде комплексов проблемно-ориентированных программ для проведения вычислительного эксперимента», п. 5 «Комплексные исследования научных и технических проблем с применением современной технологии математического моделирования и вычислительного эксперимента» и п. 7 «Разработка новых математических методов и алгоритмов интерпретации натурного эксперимента на основе его математической модели».

Публикации. Основные результаты работы опубликованы в 2003-2010 гг. в двух монографиях, в 20 статьях в отечественных и международных реферируемых журналах, а также в статьях в нереферируемых журналах и в материалах конференций.

Структура диссертации. Работа состоит из двух томов. Первый том содержит введение, пять глав, заключенные, список литературы из 221 наименования и приложения. Объем первого тома составляет 244 страницы, в том числе 60 рисунков и 7 таблиц. Во второй том (128 страниц) вошли тексты программ для ЭВМ, которые использовались для расчетов в ходе написания работы.

СОДЕРЖАНИЕ РАБОТЫ

Во Введении обосновывается актуальность темы диссертационного исследования, формулируются цели и задачи, кратко излагается ее содержание, перечисляются результаты, которые выносятся на защиту. Задачи описания динамики твердых тел, которые ведут себя подобно медленнотекущим жидкостям, в общем случае нелинейны, и за исключением специальных случаев не существует аналитических методов их решения. Это обуславливает бурное развитие методов математического моделирования динамики льдов, вулканической лавы и других природных объектов, которые можно условно отнести к медленнотекущим объектам. Отмечено также, что в связи с тем, что крупномасштабные процессы в атмосфере, океане и криосфере Земли нельзя повторить в естественных условиях, математическое моделирование неизбежно становится единственной возможностью проверить те или иные гипотезы или оценить надежность оценок, относящихся к прошлому или будущему.

В первых двух главах работы приведены основные сведения о физико-механических свойствах льда и о математическом аппарате их описания, а также об основах построения математических моделей течения льда. Целесообразность разделения вводного материала между двумя главами логически обусловлена его спецификой. В Главе 1 кратко рассмотрены физические и механические характеристики природного льда, математический аппарат описания механических свойств льда и его течения. При написании

главы мы опирались на классические монографии У.С.Б. Патерсона (W.S.B. Paterson), P.C.B. Ван дер Веена (R.S.V. van der Veen), П.А. Шуйского, И.А. Зотикова, Б. Ле Меоте (В. Le Mehaute), JI. Ллибутри (L. Lliboutry), К. Хуттера (К. Hutter), Р. Греве (R. Greve) и X. Блаттера (Н. Blatter), в которых, в которых изложены основы физики, механики, и термодинамики ледников, гидродинамики, медленного течения твердых тел, дано строгое математическое изложение динамики ледниковых покровов.

Лед в горных ледниках и континентальных щитах Гренландии и Антарктиды находится в состоянии постоянного движения. Скорость движения во внутренних областях континентов имеет порядок от нескольких сантиметров до нескольких метров в год, а на окраинах может достигать десятков и сотен метров в год. Относительно высокие скорости движения отличают лед от других природных твердых тел. В общем случае, математические модели строятся в предположении о том, что лед представляет собой вязкую, неньютоновскую, теплопроводящую, несжимаемую жидкость.

В трехмерном пространстве три ортогональных вектора напряжений (поверхностных сил, действующих на элементарный объем льда) в каждой точке можно разложить на три независимых компоненты: одну нормальную и две тангенциальных, ау, где i,j = x,y,z. Тензор напряжений

/ °ХХ °ух о Л

°УУ °гу

К°хг °ух

(1.1)

симметричен, и, следовательно, только шесть его компонент независимы. Тензор (c'y) можно представить в виде суммы шарового тензора {<Уkk) и девиатора напряжений (т^

= ХЦ + О-2)

где ôy - символ Кронекера. Деформации льда практически не зависят от гидростатического давления, и соотношение, связывающее скорость деформаций льда с приложенными напряжениями, также не зависит него. Считается, что плотность льда постоянна, const, и динамика его описывается уравнением Навье-Стокса для несжимаемой жидкости

p^-^-Vp + nbv + F (1.3)

где р .-. давление, F - внешние силы (сила тяжести), вязкость г] нелинейно зависит от температуры. Поскольку, как правило, поле температуры не известно a priori, то говорят о сопряженной термомеханической задаче. Величину, обратную г), называют текучестью:

где /(]т|) - функция ползучести, г - девиатор тензора напряжений. Функция А = А(Т,р) задается в форме закона Аррениуса

Л = £аехр(-б//?7'*), (1.5)

где - активационная энергия ползучести, Г*- температура с учетом поправки на давление. Константы, входящие в (1.5) и в последующие выражения, собраны в таблице I. Qи а зависят от пороговой величины Т : Т* < 263,15*" а = 1,14 х 10~5 (3 = 60,

(1-6)

Т г 263,15А: а = 5,47х10_:> 6 = 139, а масштабный множитель Е - определяется особенностями применения (1.5).

Таблица 1 - Список физических констант, используемых уравнениях в главах 1-3.

Символ Константа Значение Единица измерения

а Коэффициент для расчета добавочной температуры льда 8,7 х 10 4 Па"3 год"1

Я Удельная теплота плавления 3,35 х Ю5 Дж кг"'

А Плотность льда 910 кг м"

Ра Плотность астеносферы 3300 кг м'

Рг Плотность подстилающих пород 3300 кг м"

р* Плотность морской воды 1028 кг м"

Ао Постоянная в законе течения льда 9,302 х 107 КПа'3 год"1

Ам Коэффициент трения 1,8 х Ю-'0 Н"3 год"' м8

С Постоянная в законе течения льда 0,16612 К"

О Коэффициент жесткости 10" Нм

А, Коэффициент диффузии 0,5 х 108 Зависит от м2 год"'

Е Коэффициент усиления условий эксперимента -

ср Теплоемкость льда 2009 Дж кг"1 К'1

с, Теплоемкость подстилающих пород 1000 Дж кг"1 К"1

& Ускорение свободного падения 9,81 м с"2

к Показатель степени в законе течения льда 1,17

кг Коэффициент теплопроводности подстилающих пород 1,04 х 108 Дж м"' К"' год"1

1 Удельная теплота плавления 3,35 х 105 Дж кг"'

0 Активационная энергия ползучести 7,88х 104 Дж моль"'

Л Универсальная газовая постоянная 8,314 Дж моль"' К"'

То Температура плавления льда 273,15 К

Тс Критическая температура в законе течения льда 273,39 К

При соблюдении условия изотропности и несжимаемости льда, главные оси тензора скоростей деформации параллельны осям тензора напряжений, а

компоненты обоих тензоров пропорциональны друг другу. Эффективные скорости деформации связаны с эффективными напряжениями Законом Глена:

¿,-АтГхтц. (1.7)

Параметр п в (1.7) определяет характер зависимости скоростей деформаций от напряжений. Значение п=3 соблюдается для широкой полосы напряжений в природных процессах.

В Главе 2 изложены основные закономерности макродинамики

ледниковых щитов как комплексных природных систем. Уравнения,

описывающие состояние ледника, формулируются в виде законов сохранения

массы, импульса и энергии. В декартовой системе координат (х,у,г), в которой

ось г направлена вертикально вверх, а оси х и у параллельны поверхности

геоида, уравнения выглядят следующим образом :

V • V = 0 (2.1)

¿у _ д. =Уст +

Л

Л}

(2.2) (2.3)

где V - вектор скорости течения льда с компонентами и, V и и/, соответственно, по осям х, у и г, 7) - температура льда, Ф - поток тепла за счет внутреннего трения. Считая, что по абсолютной величине градиенты вертикальной скорости много меньше градиентов горизонтальной скорости, компоненты тензора скоростей деформации можно выразить через градиенты скорости течения

(2.4)

При описании течения льда, локальными и конвективными производными в (2.2) пренебрегают. Считая, что да ^¡дх = 0 и да1у/ду = 0, получим да^/дг = р?, после чего можем переписать уравнение (2.2) в скалярной форме:

/ • £хх еху . \

с с сух суу *Л II

\£1Х £гу

ди 1 ди лЛ 1 ди

дх г\ду + дх) 2дг

1(ди лЛ дv 1 ды

2\ду + дх) ду 2 дг

\ди \дv дю

2дг 2 дг ~дг

дх до

да г

да.

У"

ду да

уу

дг да

= 0

уг-

дх дг

ду

дг

= 0

= А£

Исключая гидростатическое давление выражение для девиатора напряжений г,

из

ч

тензора

напряжении, где

кк .

(2.5)

получим

- символ

Кронекера. Подставив его в первые два уравнения (2.5), получим после некоторых алгебраических преобразований:

д_ дх д_ ду

+ Гуу |+

гУ

хуу + х х

, ^.су = 0

дх ду дг

= 0

ду ду дг

(2.6)

Пренебрегая атмосферным давлением в левой части третьего уравнении из (2.5) и интегрируя его от поверхности щита г=з(х,у) до подстилающей поверхности

о

г=0, получим сгц= / Р^Лг = -р$а(х,у). Считая, что на поверхности

-у)

напряжение сггг|г,5(х;),) =0, и подставляя о1г=-р^{х,у) в первые два

уравнения (2.5), получим уравнения для девиаторов напряжений в приближении неполного второго порядка:

; (2гх* + гуу )+ —^ + —5 = А* ■-

дх д_ ду

(2г^+тх

>

ду дх

дх д$

ух ^у;

ду дг ду

(2.7)

Физический смысл (2.7) заключается в пренебрежении компонентами вариационного напряжения. Подставив в (2.5) компоненты тензора скоростей деформации из (2.4), получаем выражение для эффективной вязкости:

диск 1(ди <ЛЛ2 \ldli + дхду + ~Ь\ду + дх) + 4\<%

Н—

+ 4\о£

(1-»)/2л

• (2.8)

Подставляя закон Глена (1.7) в (2.7), и принимая во внимание (2.8), получаем систему нелинейных эллиптических уравнений для компонент горизонтальной скорости течения льда в приближении неполного второго порядка:

д_ дх

д_. ду

„ , _ ди <Л>

Ъх)

„ ¿V ди 2—+ — ду дх

д + —

ду

дх

(ди ди ~ду + ~дх

ди <Л> ^ ду + дх

д + —

дг

д + —

дг

ди

¿V"

л

дх

-■Я 8

дх

(2.9)

Краевое условие на нижней границе - скорость базального скольжения. На верхней границе вводится условие равенства нулю суммы всех напряжений

(2тгг+2г„„^+г —-аг, =0

(2ту

■2ТАх + Т

'ду

4 ду д..с

(2.10)

Л-"--0

Выразив в (2.10) напряжения через градиенты скоростей, получим

'4 — + 2 —'

дх ду,

ds (ди + ЛЛ efe дх дх)ду

ди

~dz

.dv -du\ds 4 —+ 2— — + ду dxjdy

ди dv — + —

кду dxi

ás dv _ q дх dz

(2.11)

Пренебрегая первым и вторым членами в обоих уравнениях системы (2.7), получаем приближение мелкого льда или SIA (Shallow Ice Approximation):

дх

xz

dz дх

= P¡g~

ds

дх

v ds m PiS~~

(2.12)

дг ду

В разделе 23 разобрано строение стандартной комплексной модели ледникового щита, которая используется в дальнейших расчетах для генерирования граничных условий для региональной модели. Для современных моделей ледниковых щитов характерна блочная архитектура. Условно модель можно разбить на 7 блоков (рис. 1). Это два служебных блока: для организации ввода начальных данных, различных параметров и управляющих рядов (блок входных данных) и для вывода результатов расчетов (блок вывода).

Рисунок I - Основные блоки в модели континентального ледникового щита и организация взаимодействия между ними

Климатический блок (блок климатического форсинга) служит для расчета приземной температуры воздуха, составляющих баланса массы ледника, скорости таяния шельфового льда. Три блока необходимы для расчетов динамики соответственно ледникового щита, шельфовых ледников и подстилающих пород. Отдельный блок введен для расчета положения линии налегания - границы между ледниковым щитом и шельфовым ледником.

Уравнение сохраниения массы (2.1) преобразуется в уравнение, описывающее изменение локальной толщины льда:

У-(у „Н)+МЛ-М„, (2.13)

где \!г - вектор осредненной по вертикали горизонтальной скорости течения

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

Через верхнюю границу ледникового щита осуществляется его связь с атмосферой. Температура поверхности льда принимается равной приземной температуре воздуха. У верхней границы ледникового щита задается приземная температура воздуха Г5, которая считается равной температуре льда на поверхности Т(\шг =Т$ и служит граничным условием для (2.3). В качестве нижнего граничного условия для (2.3) задается градиент температуры. Тепло здесь поступает за счет трения при базальном скольжении, если последнее имеет место, и за счет притока тепла из подстилающих пород.

В отличие от ледникового щита, шельфовый ледник находится в гидростатическом равновесии с окружающим океаном. Если уровень моря падает ниже равновесной отметки, то шельфовый ледник переходит в категорию континентального и наоборот - происходит миграция линии налегания. Таким образом, через изменения уровня Мирового океана реализуется один из двух механизмов привязки динамики ледникового щита к динамике палеоклимата (либо к проекциям изменений климата в будущем). В противоположность континентальному льду, шельфовый ледник испытывает пренебрежимо малое трение, так что горизонтальными тангенциальными напряжениями можно пренебречь. Скорости деформации (и компоненты скорости течения) шельфового ледника можно считать не зависящими от глубины, поскольку углы наклона поверхности крайне малы. Таким образом, движущей силой здесь будут горизонтальные градиенты нормальных и тангенциальных напряжений в горизонтальной плоскости. Поскольку механика деформаций на шельфового и покровного ледников определяется различными напряжениями, между ними вводится переходная зона, которая расположена в окрестностях линии налегания.

В разделе 2.4 приведен краткий исторический обзор математических моделей ледниковых щитов и дана оценка перспективных направлений их развития.

В Главе 3 рассмотрена разработанная автором математическая модель течения льда в ограниченной пространственной области. Иногда в русскоязычной литературе для обозначений моделей подобного типа применяется не вполне удачный термин «угнезденные». Мы же остановились на названии «региональная». Применение модели выбранного типа в настоящем исследовании обусловлено тем, что ограничения, накладываемые приближением мелкого льда (2.12) на пространственное разрешение модели, не позволяют адекватно описывать процессы на окраинах ледников, вблизи ледоразделов и в районах со сложным рельефом подстилающей поверхности. Помимо этого, характерное пространственное разрешение в рамках этого

приближения ограничено ~20 км, что исключает детальное описание полей скорости и деформаций - льда в непосредственной окрестности скважин и ассимиляцию в модель полей топографии поверхности и ложа с большим пространственным разрешением. Альтернативой для региональных исследований является аппроксимация неполного второго порядка (2.8)-(2.9). Поскольку детальные расчеты необходимы в относительно небольшой по площади области, то для оптимизации использования вычислительных ресурсов целесообразно разбить задачу на две часта:

1. Имитация динамики всего ледникового щита на комплексной модели с относительно грубым пространственным разрешением;

2. Решение уравнений течения более высокого порядка аппроксимации в ограниченной области (встроенном домене) с большим пространственным разрешением.

Итак, оптимизация расчетов при большом пространственном разрешении достигается разделением решения задачи между двумя субмоделями (рис. 2):

1. Моделирование эволюции динамики ледникового щита на комплексной Б1А-модели с относительно грубым пространственным разрешением («большая модель», БМ).

2. Решение уравнений стационарного течения льда более высокого порядка аппроксимации в ограниченной области (встроенном домене) с большим пространственным разрешением для задаваемых («малая модель», ММ).

Роль БМ заключается в генерировании граничных условий для ММ. В БМ производится расчет эволюции топографии поверхности и толщины льда, полей температуры и скорости, включая скорость базального скольжения, скорости базального таяния, реакции подстилающей поверхности на меняющийся вес ледникового щита, климатических условий у поверхности, баланса массы. Расчеты охватывают область, заведомо превышающую размеры континента и шельфовых ледников. Чтобы избежать дополнительных сложностей, связанных с описанием процесса откалывания айсбергов, считается, что вся область БМ представляет собой лед. Такое упрощение вполне оправдано, поскольку в цели моделирования не входит реконструкция пространственной протяженности морских льдов, и динамическое взаимодействие с океаном не описывается. В рамках БМ решаются уравнения сохранения массы, количества движения и энергии. Домен БМ (рис. 2) имеет квадратную форму 281x281 узлов (5600x5600 км) с вертикальным разрешением 30 слоев, которые расположены неравномерно, уплотняясь по направлению к основанию щита.

В соответствии с изложенным выше принципом построения региональной модели или малой модели (ММ), нами выделены области (домены), в которых решаются уравнения течения льда с большим пространственным разрешением (рис. 3). Размеры домена 1 600x400 км, пространственное разрешение 2,5 км. Размеры домена 2 960x640 км, пространственное разрешение 5 км. Вертикальное разрешение составляет 101 равномерно расположенный слой.

6М:

• Ледниковый щит/шельф/литосфера ■ Свободная линия налегания

* Термомеханическое сопряжение

а Разрешение: 20x20 км по горизонтали, 30 слоев по вертикали

ММ:

• Новейшие данные по топографии ложа и скорости аккумуляции

• Разрешение: 2.5x2.5 к и 5x5 км по горизонтали, 101 слой по вертикали

■ Размер доменов 600 х 400 км и 960 х 640 км

О 140 160 180 200 220 240 260 260

Стационарное поле скорости ОБ

(¡Г + —— + Яу Л: А

<4: А -Ж— ду

ду А

Два домена

БМ: Большая модель (внешний домен)

ММ: Малая модель (встроенный 40 домен)

000 2000

Рисунок 2 - Схема, поясняющая строение региональной модели

Рисунок 3 - Положение встроенных доменов региональной модели (выделены серым цветом) на карте Антарктиды: (I) домен 600x400 км с пространственным разрешением 2,5 км, (2) домен 960x640 км с разрешением 5 км. Во втором домене показан профиль радиозондирования между станциями Конен и Купол Фуджи, совпадающий линией леОораздела

кт

О 1000 2000

Конен

Купол Фувкн

Восгок

Дискретизация уравнений ММ и процедура численного решения изложена в разделе 3.2. Раскрыв скобки в уравнении (2.9) и упорядочив результат, получим следующую систему нелинейных эллиптических уравнений течения льда

¿и

йс\ дх)

+ 2—1 —1 + —| — + ду) ду

д( ду\ д( ди\

*(ч*Ыч*гм

4—1 —1 + 2—1 + + —) + — { —) = дуУдуУ дуУ дх) дху ду) дхУ дх) дгУ дг) Р

А дх

А дх

(3.1)

Введем систему координат (*',)>',£), где вертикальная координата нормирована на толщину льда, £ = (5 - £)/Н:

д( Л<\1 / д ( ди) . \д( ¿»«V

д

д_( ди

М ¿У

ди

+ а.

д ди

+ а,

ди

+^ + Фу ) = ^ - ЪГ)С

дК

-2а

дх

% К

+ а,

+ 4 а2г

д

д I ди

«га

ди

я2

д_(

д_1 ¿V

ьуэу.

-2а,

-п-КУду.

-Ъахау

д( <М1 \д( ЛЛ

-а,

ЛЛ] д ( лЛ

.(3.2)

Второе уравнение из (3.1) преобразуется аналогичным образом. Уравнение для эффективной вязкости (2.8) в безразмерных координатах можно записать как

1 а/Т* Г"" Г7 = -Л(Г ) х

.(ди\2 ди

У

л+(4«?+Дг+1

дх

ду дv ( . 2 1 +4 а„ —— + |4 а;, + —- + 1

+4 а

ду ди

Н'

Н<

ди

К)

+ 5| —

№

дидУ

ди дУ

— +4--+ 4 аг--+

дх ду

д^ду

ди ду . диди „ ду ду .виду

+ Аа.а,. —— + 2 а„—— + 2а, —— + 2--+

* у д£ду д^дх дудх

ди ду „ дУ ди _ ди дУ

+2 а„--+ 2аг--+ 2а„а..---

уд£дх х д^ду у

(3.3)

Легко заметить, что (3.2) состоит в основном из операторов вида

где (р = иу, Я = х,у£, I * у. Произведя соответствующую замену, получим уравнение для и

4£2(кд) + 4ахЧ>{и,х£) + 4ахУ{й£,х) + ^4 а2х + а) + +

+£2{«,у) + ауЧ>{и,у£) + ^(«¿.у) + + 4^ -

(3.5)

Л>

и аналогичное выражение для V

(3.6)

Щу,у)+4ауЧ>(у,у£)+4ауЧ>{ч£,у) + [4а2у + а] ++ +ПМ + + ауФ(у,^) + (Ьх + =

= - 2Ч»(и,у,лс) - 2а ^(и.у.С) - 2в,>Р(и^) - ЗахауП{и£) --4>(и,х,у) - - ахЧ>{и£,х) - 3 с^ц^-

Приведя уравнения (3.1) к такому виду, мы упрощаем процедуру численного решения и программирования алгоритма решения с использованием смещенной пространственной сетки. Система нелинейных уравнений (3.5)-(3.6) решается методом последовательных подстановок. Обозначим номера итераций через /, /+1, 1+2, . Каждая из итераций распадается на три шага:

1. Расчет эффективной вязкости »у в (3.3) с использованием и и V, найденных на предыдущей итерации. На первой итерации используются локальные решения для и и V в соответствии с аппроксимацией мелкого льда, но, в отличие от БМ, рассчитанные на пространственной сетке с высоким разрешением.

2. Расчет и. При этом V и г) в первом из уравнений (3.1) считаются независимыми от и, и, таким, образом уравнение считается линейным. Первые и вторые производные рассчитываются как центральные разности внутри области, и как односторонние разности на границе. Уравнения и граничные условия для и и аналогичное для V в координатной системе (х',у',£), будучи переписанными в конечно-разностной форме, преобразуются в систему линейных алгебраических уравнений. В матричной форме их можно записать как Ах = Ь, где х = {х,,...,х„} - неизвестная переменная, А^Л^хЛ^х/У. -число переменных, Ь = {Ьх,...,Ьп,...,Ьц} - матрица коэффициентов, и А -квадратная матрица элементов а,у. В случае, если ¡А|*0, единственным решением будет х = А"'Ь. Обращение матрицы А требует около Л'3 операций. Однако, поскольку лишь незначительное число элементов ац отлично от нуля,

то матрица А является разреженной. Для решения системы применяется метод

сопряженных градиентов.

3. Расчет v проводится аналогично расчету и.

Сходимость решения оценивается, исходя из значения параметра остановки Q, который должен быть меньше заданной величины допуска toi: тах|Еттх|/max|x| < toi.

Найдя на втором и третьем шагах итерации / оценки и иг, которые обозначим как и* и v*, произведем обновление rj. Эти итерации будем повторять до тех пор, пока не будет найдено стационарное решение. Для ускорения процесса сходимости применим релаксационный алгоритм. На каждой итерации /, /+1,1+2 нелинейной системы генерируются решения и', и'+|, им . Допустим, что эти решения можно представить в виде и'+'= и' +с', где с - корректирующий вектор. Неизвестная ошибка (расстояние) между точным решением и решением и' в текущей итерации до коррекции обозначим как е. Допустим также, что последовательность ошибок е', e/+I, е'+2 сходится на прямой линии в корректирующем пространстве, и что jjc'j =|е' — ei—1|. Это

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

Решения и и V считаются промежуточными между итерациями / и /+1. Во время выполнения итерации /+1 применяется модифицирующий корректирующий вектор ci-сl/a: u'+1= u' +с'/а, где с' вычисляется явным образом как u'-u'. Угол со между двумя последовательными корректирующими векторами рассчитывается как

Г . . Т '//.. -мы - . .л'

(3.8)

3s[(c'f-c'-

fti = arccos|(c*| -с" •/111с'|МсМ

Этот угол должен быть близок к 0 или к л. Тестовые расчеты показали,что наилучшая сходимость получается, если корректирующий множитель равен 0,5 в случае, если ш > Зл/4, и 2, если со < л/4.

Операторы й и *Р из (3.4) содержит эффективную вязкость г] и пространственную производную компоненты скорости, которые целесообразно рассчитывать на смещенной сетке (рис. 3-4). Операторы ЙиЧ'из (3.5) и (3.6), записанные в конечных разностях, выглядят как

1 . 1 +,7. 1 . |

Q(u,x)н-ir,-) = 2 : щ ,, v ' дх\ дх)

2 Ах

»

»7i i + Ч i 1 + Ч 1 . 1 + V. 1 . | П. 1 . 1 + П. 1 . I

-гП ц "ïJ~i ц

^ u'-i+ "'M

(3.9)

2Асг w 2Дх

дх\ ду

4ДхДу

1 1 . I +»7 1 I -V | . 1 ~»7 I I

■ - ,—ц.

2 2 \ /

(+-,/+-- I—,/— )—,;+— * ./—

\ 2 2 2 2 2 2 2 2/

1. 1 . 1 -Ч. 1 . 1

\ 2У 2 2У 2/

»7. I 1 ~г\. 1 . 1

2 2 2/

(-1,7

1.г *?. 1 . I

I—,1— 1+-Л— \ 2 2 2 2/

п. 1 . 1 -»?. 1 . I

/—— I—,;+—

\ 2 2 2 2/

+,7. I . 1 "/+1,у>1 + V. I . 1 "м,;-1 ~ 1 1 МЫ,/-1 1 . 1 "|-1,у+1

l*-J*■ 2 2

, (3.10)

Ыу* и+1

Ыцк+1

и«*

и» I,/+!,*+!

Рисунок 4 - Трехмерная ячейка смещенной сетки, применяемой для решения системы уравнений 3.5-3.6

Чм,¡*1,к

Дх

Рисунок 5 - Проекция на плоскость (х,у) смещенной сетки, примененная для решения системы уравнений 3.5-3.6

Применение смещенной сетки позволило подавлять случайные возмущения, возникающие в процессе численного решения. Расчет эффективной вязкости на смещенной сетке приводит к значительной экономии машинного времени: в диапазоне значений допуска 10~5-10"7, которые использовались на практике для достижения стационарного решения, время расчета с использованием смещенной сетки сокращалось на 20-50% по сравнению с расчетами на обычной сетке. В сериях тестовых экспериментов на моделях с упрощенной геометрий ледниковых объектов экономия времени достигала семикратного значения. Процесс сходимости решения показан на рис. 6.

700

'I 600

со

§"500 Ь 400 Е 300

I 200

С

5 100 о

2 о*

,3 2 Рисунок 6 - Сходимость решения системы (3.5)-(3.6)

л

10"4

Номер ''нелинейной" итерации

Для встраивания ММ в БМ (раздел 3.3) решаются две группы задач.

1. Объединение БМ и ММ в одном модельном коде и обеспечение потока данных в одном направлении - от БМ к ММ. В этих целях на боковых и нижней границе домена ММ производится интерполяция компонент горизонтальной скорости и передача их в качестве краевых условий от БМ к ММ. Температура в ММ не вычисляется, и ее передача производится не только на границе, но и внутри домена ММ в совпадающих узлах с последующей интерполяцией на всю сетку ММ. Поскольку поле температуры транслируется из БМ, то очевидно, что и скорость базального скольжения целесообразно рассчитывать в БМ и передавать таким же образом, как и температуру, в ММ.

2. Обеспечение плавного перехода от БМ к ММ. В ММ ассимилируются данные полевых исследований: толщина льда и скорость аккумуляции (которая используется для расчета вертикальной скорости). Поле топографии поверхности также не обязательно совпадает в обеих моделях. Их пространственное разрешение (2,5 км) в 8 раз превышает пространственное разрешение БМ (20 км). По сути дела, мы имеем дело с двумя наборами или базами данных, которые, вообще говоря, не совпадают ни в общих узлах сетки, ни на границе. Для того чтобы плавно вписать область ММ в область БМ, вводится переходная зона, в которой поля из БМ постепенно трансформируются в поля ММ с помощью бикубических сплайнов.

В численных экспериментах синхронизация БМ и ММ производится через фиксированные промежутки времени (20-100 лет). Для совмещения полей во внешнем и во встроенном домене используется одна из двух схем

объединения процессов в модельных доменах. Запишем эволюцию модельных полей переменных п, А(п,г)$т как

А(п,0р = А(п,1)о5Ь! + А'{п,1), (3.11)

где Л(я,Г)^ - поле наблюденных значений переменной и во встроенном домене (в ММ), А'(п,1) - поправочный член: разница между значениями переменной п в БМ и ММ, также зависящая от времени.

Схема 1 (ОС1). Предполагается, что матрица коррекции (или матрица аномалий) вычисляется на каждом шаге, и представляет собой разницу между значениями поля БМ в момент времени ( и начальным значением поля в этом же узле в момент времени /0, которые по определению совпадают с наблюденными полями в базе данных БМ:

А'(п,1) = А(л,0Г -А°вЫ(п,10). (3.12)

Далее, эта матрица коррекции накладывается на наблюденные значения соответствующих полей. В результате изменения полей ММ такие же, как и полей БМ, хотя по абсолютной величине сами поля различны. Поскольку в БМ и ММ используется одно и то же поле топографии поверхности, подстановка и=5иг сводится к интерполированию модельного поля s(x,y) из сетки 31x21 узел в сетку 241x161 узел.

Схема 2 (ОС2). Матрица коррекции (или матрица аномалий) вычисляется на каждом таге, и представляет собой разницу между значением поля в модельный момент времени / и момент, соответствующий окончанию численного эксперимента в этом же узле 1еги1

А'(п,1) = Л(«,/)Г - АГЫет1). (3.13)

Использование ОС2 предполагает, что сначала производится прогон только БМ, в котором генерируются поля Адт(п,1ет1), после чего в совместном прогоне совмещенной модели эти поля подставляются в (3.12). Модельные поля п в ММ совпадают с наблюденными в начальный момент интегрирования в случае использования ОС1 и в конечный момент интегрирования в случае использования ОС2. ОС1 целесообразно использовать в случае, если поставлена задача описания эволюции поля скоростей в соответствии с эволюцией модельной топографии в БМ. Однако, если необходима точная привязка поля топографии поверхности к современным условиям, предпочтительнее ОС2.

Рассчитанные в соответствии с 0С2 поля скоростей имеют минимумы вдоль гребня, соединяющего станции Конен и купол Фуджи (рис. 7). Окрестности станций Конен и станции Купол Фуджи различаются с точки зрения топографии поверхности: станция Конен расположена вблизи седловины в конце пологого гребня, а станция Купол Фуджи очень близко от вершины топографического купола. Окрестности станции Конен имеют также более сложную топографию ложа. Непосредственно в седловине амплитуда снижается до величины менее 10 см/год, а приблизительно в 1200 км вверх по течению в районе купола Фуджи, - менее 1 см/год. Верификация модели была произведена сравнением с данными СРБ-наблюдений. Амплитуда и

направление модельной скорости составили соответственно 73 см/год и 279.4°., наблюденной 75,6 см/год и 273.4" .

о м 100

Рисунок 7 - Направления векторов скорости течения льда на поверхности ледникового щита в доменах 1 (вверху) и 2 (внизу). Длина стрелок относится к амплитуде скорости как 1:500. Амплитуды скорости течения льда на поверхности ледникового щита в доменах I (вверху) и 2 (внизу). Шкалы справа даны вм/год

Главы 4 и 5 посвящены практическим приложениям математической модели. Была проведена серия численных экспериментов по расчету динамики льда в области на Земле королевы Мод при различных граничных условиях. Для интерпретации результатов численных экспериментов был разработан метод обратного отслеживания, с помощью которого был решен ряд задач, представляющих практический интерес.

Через верхнюю границу системы ледниковый щит-шельфовый лед осуществляется связь системы с атмосферой. Поле приземной температуры воздуха разбивается на две компоненты: меняющуюся со временем -климатический сигнал и стационарную часть, соответствующую современному климату. Стационарная часть определяется с учетом изменений температуры с широтой и абсолютной высотой места. Климатический сигнал рассчитывается по косвенным данным - длинным рядам <$£> из ледяного керна, исправленным на изменение среднего изотопного состава Мирового океана, и с использованием переводного коэффициента с)Тг/ддО = 0,166°С/%о. Величина климатического сигнала считается постоянной в один и тот же момент времени над всей областью интегрирования. Считается, что уровень Мирового океана менялся пропорционально вариациям б18<9 в океанических донных осадках. Мы использовали длинный изотопный ряд и масштабировали его, исходя из того, что во время последнего ледникового максимума уровень падал до отметки -130 ми превышал современный на 6 м во время предыдущего межледниковья около 122 тыс. лет назад (рис. 8).

Применение модели для палеореконструкций требует расчета возраста льда. В основе эффективного метода решения поставленной задачи лежит представление о том, что возраст льда является консервативной примесью, переносимой полем скорости течения, которое, в свою очередь, генерируется математической моделью. Возраст льда рассчитывался методом обратного отслеживания, который представляет собой комбинацию методов Эйлера и Лагранжа. Метод предполагает нахождение траектории частицы льда от ее положения в виртуальном ледяном керне (конечная точка на временной отметке 1=0 - рис. 9) до места ее происхождения на палео-поверхности ледникового щита.

В соответствии с требованиями метода обратного отслеживания, численные эксперименты проводились в два этапа. На первом этапе БМ прогонялась в течение заданного модельного времени «вперед», от временной отметки 1=970 тыс. лет назад до временной отметки £=0. По мере выполнения форвардного эксперимента, данные по скорости течения льда, топографии щита и други характеристикам во встроенном домене сохранялись с заданной периодичностью. На втором этапе модельное время запускалось «назад» - от *=0 до (=-(*, где -г* - временная отметка в модельном «прошлом», до которой производилось восстановление различных характеристик течения льда и состояния ледникового щита вдоль траекторий виртуальных частиц льда (по сохраненным на первом этапе данным) Процесс отслеживания трассера т заканчивался в момент пересечения им поверхности ледникового щита в момент времени 1т появления трассера на поверхности с выпавшими осадками в «точке происхождения», то есть в точке пространства с координатами (Хт,Ут,0) (рис. 10). Фиксируя момент времени, когда трассер пересекает поверхность щита, мы получаем возраст льда на той глубине в керне, с которой стартовал трассер (рис. 11). Так, самый «глубокий» трассер, расположенный на 89% глубины, был изначально отложен на поверхности щита в точке, которая расположена на расстоянии 184 км по прямой от скважины, и его возраст составил 169,9 тыс. лет. Обе кривые очень близки, и разница в оценке возраста не превышает нескольких лет на всем протяжении керна, за исключением двух сегментов, 1700-1850 м (60-75 тыс. лет.) и 2300-2470 м (120-167 тыс. лет). Фактически, модельная кривая возраста, будучи очень близкой к шкале ЕОМЬ1

Рисунок 8 - Климатический форсинг. Горизонтальная ось -время (тыс. лет); левая вертикальная ось - вариации температуры воздуха ('С), правая вертикальная ось -вариации уровня моря (м).

-700 -600 -500 -400 -300 -200 -100 0

за исключением двух сегментов, свидетельствует о внутренней сбалансированности процессов деформации льда, локальной температурной истории и амплитуды гляциально-межгляциальных вариаций скорости аккумуляции.

Функция динамического сжатия - один из важнейших параметров, характеризующая особенности растекания льда. Он рассчитывается как отношение толщины годового слоя льда в керне к количеству выпавших осадков. Возраст льда A(z) на глубине z, скорость аккумуляции P(z) и функция динамического сжатия R(z) связаны соотношением

где R(z) = А(г)/А0 и A<f=PA(t) - начальная толщина годового слоя на поверхности щита.

Рисунок 9 - Схема, иллюстрирующая метод обратного отслеживания: траектории трассеров в вертикальной плоскости от начального положения в виртуальном керне до поверхности щита в прокручиваемом назад во времени потоке льда

140 1«! 142 143 144 145 148 147 148 U9 150 151

Рисунок 10 - Координаты происхождения трассеров на поверхности щита, изначально расположенных в виртуальном керне каждые 1% относительной глубины на фоне современного поля амплитуд (м/год) и направлений вертикально-осредненной скорости потока. Маркеры показывают время выхода трассеров на поверхность (возраст льда в керне)

Для расчета функции динамического сжатия в ходе пост-экспериментальной обработки данных скорость аккумуляции (модельная) ставится в соответствие толщине годового слоя в керне (также модельному) в соответствии с выражением = (//(г))/{^')> ГДе ' средняя скорость палео-

аккумуляции от момента времени до момента времени Ь, (Н(г)) - средняя толщина слоя годового слоя в виртуальном керне. Функция динамического

Расстояние вдоль оси потока

яектории трассеров е вертикальной

движения

сжатия была использована для расчетов скоростей палеоаккумуляции по данным радизондирования (раздел 4.5).

о

500

2000

1000

1500

Рисунок 11 - Модельный возраст льда в виртуальном керне. Для сравнения показан возраст согласно шкале ЕВ МЫ [203]. По горизонтальной оси отложено время, годы, по вертикальной -абсолютная глубина скважины на станции Конен

2500

0 40

ео

1го 160x1а1

Проникающее радиозондирование позволило выделить несколько внутренних слоев, которые удалось проследить от начала до конца профиля за исключением отдельных участков, где радиоэхограмму не удалось расшифровать. Глубина залегания внутреннего слоя несет информацию о возрасте льда в слое и о скорости аккумуляции в течение отрезка времени, начиная от времени выпадения вулканической пыли на поверхность щита до сегодняшнего дня. Датирование слоев на обоих концах профиля осуществлялось в соответствии с официальной хронологией на ст. Конен и Купол Фуджи. Функции динамического сжатия рассчитывались каждые 25 км в сегменте профиля 0-350 км и каждые 50 км в сегменте 350-1300 км. Численно интегрируя выражение для скорости аккумуляции

в интервалах глубин между внутренними слоями, мы получили средние скорости аккумуляци (рис. 12). Как и ожидалось, скорость аккумуляции постепенно понижается по мере движения от станции Конен к станции Купол Фуджи. Скорости аккумуляции в голоцене (первые пять кривых на рис. 12) в 1.5-2 раза выше, чем течение гляциальной эпохи (последние пять слоев).

Восстановленные скорости аккумуляции позволили скорректироать уравнение для параметризации скорости аккумуляции. Гляциально-межгляциальный контраст в скоростях аккумуляции определяется нетермодинамическими факторами, и зависит от географического положения области, особенностей синоптических процессов в регионе и других условий в прошлом, о которых мало известно. Чтобы учесть их, в уравнение для расчета скоростей аккумуляции в прошлом был введен был введен параметр /3:

¿г

(4.2)

(4.3)

змаЦ Фиг

22.47

Т, СОД г/ (01

(4.4)

В связи с тем, что два временных промежутка имеют разрыв на участке 560-800 км, /3 оценивалось отдельно для сегментов 10-560 км (14,7-72,7 тыс. лет назад) и 800-1280 км (15,8-72,7 тыс. лет назад). Средние значения на указанных сегментах составили 0,0457 и 0,0522, среднеквадратические отклонения соответственно 0,0236 и 0,0101 (рис. 13).

8.01-'

о |

=1 (и 7.0е-

О 6.0 Ёу

2 О, Ъ; 5.0р

сс; С 4.0»

з.ор

«I 1

У 2.0 Ь

а. о ^ О 0

- Ч4ЙМ

Возраст {тыс. пег)

0 5.0

.............. 5.0 ■ 79

.......... 7.9-10.3

-•-.......- 10.3-11.2 •

--------------11.2-13.1

-------13.1 14 7

----13.1 -15.84

------14 7-25.5 ]

-----15.8 25.5Н

---- 25.5 -38,0 3

'-•' ( ->'4*Д» ллА/Ч

200 400 600 800 1000 1200 Расстояние вдеть профиле (км)

Рисунок 12 - Средние скорости аккумуляции, рассчитанные методом перерасчета функций динамического сжатия вдоль профиля радиозондирования

0.04 0.02

\..... 114.7 • 72.7 тыс. лет назад, 10-560 ш: — 20-«и схолъзящее соеднее. <0> = 0.0457, о=0.0236 |— <{}>. — <3>*1о

; :

\ л

N л /! ' |

[ 1 V ]

"15.5 72.7 тыс. тт нгзгя. 300-1280 км: - 20-*и ежог.ьэящее средне«. 41> = 0.0522. с=0.ОЮ1!= <[Ы:1а

100 200 300 400 Расстояние вдоль профиля (км)

900 1000 1100 1200 Расстояние вдоль профиля (км}

Рисунок 13 - Оценка /3 для сегментов 10-560 км (слева) и 800-1280 км (справа)

В разделах 5.1 и 5.2 приведены сведения о методах восстановления температуры воздуха и уровня Мирового океана в прошлом по данным глубокого бурения ледниковых щитов и морского дна, а также рассмотрены спектральные характеристики этих рядов. В разделе 5.3 разбирается метод исключения из изотопно-температурных рядов топографического смещения, рассчитываемого на основе результатов моделирования. Топографический эффект (смещение) имеет две составляющие - локальную и адвективную. Локальная составляющая обусловлена вариациями абсолютной высоты поверхности ледникового щита во времени. Адвективная составляющая обусловлена движением льда вдоль осевой линии потока вниз по течению относительно места, где снег первоначальной выпал на поверхность щита. Чтобы вычесть локальную и адвективную топографические составляющие из вариаций изотопного ряда, нужно, фактически, рассчитать путь, пройденный частицей льда, взятой на определенной глубине в керне, до точки

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

(5.1)

Локальная составляющая равна разнице современной абсолютной высоты Бху(0) и палео-высоты в точке происхождения частицы льда 5^/).

Сумма адвективной и локальной составляющих дает общее топографическое смещение

Л.?;о,(0 = Л^(0+д5/Ос(0- (5-3)

Следующий шаг - пересчет вариаций абсолютной высоты в вариации изотопного состава и температуры. На практике, палео-градиенты Л518о/<?5 и дд^о/дБ неизвестны, и для пересчета были использованы современные эмпирические градиенты. С учетом изложенного выше формулы для расчета суммарных поправок к рядам температуры и изотопного составы выглядят следующим образом:

л^до^+ДМО^- (5.4)

А^со^+дмо^, (5.5)

где оТ5/д8= -0,0117! °С/м и д81*о/д8= -0,00957 %о/м. На рис. 14 показаны величины локальной и адвективной составляющих вариаций абсолютной высоты, и их сумма в зависимости от времени. Для построения использована модельная хронологическая шкала.

» 140 120 100 ВО 60 40 20 О Время (тыс. пет назад) В|»ия «т назад)

Рисунок 14 - Слева: Локальная, адвективная и суммарная составляющие топографического эффекта. Справа: суммарный топографический эффект, пересчитанный в вариации температуры воздуха и изотопного состав льда в соответствии с эмпирическгши «пространственными» градиентами. По вертикальной оси отложены отклонения температуры воздуха (ось слева) и отклонения относительной концентрации (ось справа)

Умножение полученных эмпирических градиентов на вариации абсолютной высоты Д5,„,(г) дает величины топографического эффекта, который надо исключить из изотопного или из температурного ряда (рис. 14).

В разделе 5.4 рассморено влияние потока геотермического тепла на распределение возраста льда и величину топографических поправок в придонной части ледяного керна. Модельные расчеты, в которых С составляла среднюю для Антарктиды величину 54,6 мВт/м2, указывали на отсутствие базального таяния, однако, модельная температура, исправленная на давление, была чрезвычайно близка к температуре таяния. Станция Конен расположена в области, где базальное таяние отсутствует. Пограничное положение станции вблизи границы зоны таяния льда не исключает того обстоятельства, что в прошлом зона таяния захватывала место расположения скважины. Для выяснения вопроса о чувствительности возраста льда к вариациям С были выполнены две серии численных экспериментов, в которых С последовательно увеличивался на 2-30% по сравнению с принятой величиной. С помощью метода обратного отслеживания в ходе пост-экспериментальной обработки рассчитывались хронологические шкалы для льда в придонном слое при различных С. Было установлено, что когда О достигает максимального значения в серии Сге/+ 30%, область таяния распространяется на весь всю исследуемую область за исключением северо-западной части, где толщина льда существенно ниже. Интегрирование скорости базального таяния дает суммарный расход льда в придонной части (рис. 15а).

Рисунок ¡5 - Интегральный расход льда за счет таяния на нижней границе ледникового щита на станции Конен (м) в течение 740 тыс. модельных лет (а); зависимость относительного расхода льда от относительного потока геотермического тепла (б)

Суммарный расход растет практически линейно с ростом О (рис. 156). Результаты численных экспериментов позволили построить хронологии придонного льда и рассчитать топографические поправки для диапазона ... +30% (рис. 14). На глубине 99,9% относительной глубины (2.75 м от ложа ледника) в случае равенства ПГТ среднему по континенту возраст льда в керне станции Конен превышал бы 740 тыс. лет. При наиболее вероятных значениях ПГТ Сге/г5-\0% возраст льда на этой глубине составляет 385,6-410,1 тыс. лет (рис. 16 слева). Потери движущейся колонны льда за счет базального таяния в этих условия составляют 60-125 м.

Топографическая поправка к значению температуры воздуха во время последнего ледникового цикла составляет ок. -1,2 "С, а во время межледниковий -125, -240, -320, -405 тыс. лет назад она достигает половины всей восстановленной аномалии температура воздуха. Наличие базального таяния позволяет скорректировать величину поправки в сторону ее уменьшения на 10-20% (рис. 16 справа). Введение топографических поправок дает возможность гораздо более точно выделить климатический сигнал из общего ряда аномалий температуры воздуха, восстановленной по изотопному составу льда, и, таким образом, определить реальную амплитуду гляциально-межгшщиальных вариаций температуры воздуха, что является ключевым моментом для понимания механизмов, управляющих эволюцией климата.

24001

200 250 300 Возраст льда (тью. лет)

Исток геотермического тепла

Стандарт

— .5% -------.20%

.... .10%-----.25%

-■- .15%--30%

-300 -250 Время (тыс. лет)

Рисунок 16 - Возраст придонного льда (слева). Суммарный топографический эффект, пересчитанный в вариации изотопного состав льда в соответствии с эмпирическими «пространственными» градиентами. По вертикальной оси отложены отклонения от относительной концентрации (справа)

В Заключении подводятся общие итоги работы. На основании выполненного исследования разработаны теоретические положения и получены практические результаты, которые можно квалифицировать, как развитие нового направления в математическом моделировании природных ледниковых объектов.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

1. Решена проблема построения трехмерной модели медленнотекущей жидкости. Разработана оригинальная схема дискретизации нелинейных эллиптических уравнений течения льда. Численная схема решения на смещенной сетке обладает большей устойчивостью и намного эффективнее предложенных ранее. Разработанная модель была применена, в частности, для воспроизведения динамики ледникового щита в ограниченной пространственной области на Земле Королевы Мод в Антарктиде для интерпретации данных глубокого бурения в рамках программы EPICA (.European Project for ¡ce Coring in Antartica) и для расчета скоростей течения льда в северо-западной части Гренландского ледникового щита в рамках программы NEEM {North Eemian).

2. Предложен и реализован метод разделения пространственных масштабов моделируемых процессов, для чего были разработаны алгоритмы объединения суб-моделей, отвечающих за описание динамики льда с разным пространственным масштабом. Показано, что современные методы математического моделирования в приложении к задачам геофизики, гляциологии и палеоклиматологии, являются эффективным инструментом исследования. Предложенный подход может быть реализован для модельного исследования в иных географических регионах, а также в целях прогноза поля течения льда, что является необходимым для долгосрочного планирования научных исследований, в том числе глубокого бурения ледниковых щитов. После незначительных модификаций предложенная модель может быть использована также для расчетов поля скоростей других медленнотекущих жидкостей, в частности вулканической лавы.

3. Разработана методика проведения численных экспериментов длительностью до нескольких сот тысяч модельных лет, и таким образом, охватить несколько ледниковых циклов. Методика позволяет включать в численные эксперименты новейшие данные полевых исследований, данные глубокого бурения на антарктических станциях и данные, полученные из кернов морских донных осадков.

4. Разработана методика пост-экспериментальной обработки данных, основанная на обратном отслеживании виртуальных частиц в модельном потоке льда, а также соответствующие ей алгоритмы и компьютерные программы. Методика позволяет определять время, и координаты происхождения частиц льда на поверхности ледникового щита и отслеживать изменение модельных характеристик льда по мере движения частицы от места происхождения до ее конечного положения в виртуальном ледяном керне. Методика была применена для расчета траекторий движения частиц льда, для анализа динамической истории потока льда, для расчета динамического сжатия, то есть тех характеристик течения льда, которые в дальнейшем были использованы для практического использования при интерпретации физико-химического анализа ледяного керна, полученного на станции Конен (Земля Королевы Мод, Антарктида) в рамках программы EPICA. Предложенная методика является универсальной, и может применяться в иных модельных исследованиях.

5. Разработаны методика, алгоритмы и компьютерные программы восстановления скоростей аккумуляции, основанные на анализе данных радиозондирования с помощью модельных функций динамического сжатия. Методика позволила уточнить соотношение между температурой воздуха и скоростями аккумуляции на Земле Королевы Мод на протяжении последних нескольких десятков тысяч лет. Методика может применяться для обработки данных радиозондирования в областях с горизонтальной адвекцией льда.

6. Разработана методика, алгоритмы и компьютерные программы расчета неклиматических (топографических) поправок к рядам изотопного состава льда и реконструируемым по ним рядам температуры воздуха. Введение неклиматических поправок в изотопно-температурный ряд, полученный из

ледяного керна станции Конен в рамках программы EPICA, позволило уточнить реальное изменение температуры воздуха в прошлом, и, таким образом, точнее оценить масштабы климатических вариаций на протяжении сотен тысяч лет.

7. Предложена методика и проведен комплекс модельных исследований возраста льда в нижней части Антарктического ледникового щита в условия неопределенности величины потока геотермического тепла. Численные эксперименты позволили установить наиболее вероятный возраст придонного льда на станции Конен и диапазон неклиматических (топографических) поправок к изотопно-температурному ряду в нижней части ледяного керна. Одновременно определена область залегания наиболее старого антарктического льда в районе Купола Аргус, где, согласно модельным оценкам, возраст льда в придонной части превышает 1,2 млн лет.

ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ

Монографии

1. Рыбак О.О. Математическое моделирование динамики ледникового щита Антарктиды: теория, эксперименты и приложения в палеореконструкциях - М: Физматлит, 2007,223 с.

2. Рыбак О.О. Математическое моделирование эволюции ледниковых щитов -М: Физматлит, 2011, 193 с.

Журналы, рекомендованные ВАК для публикации основных научных результатов диссертации на соискание ученой степени доктора наук

3. Rvbak О.. Huybrechts P. A comparison of Eulerian and Lagrangian methods for dating in numerical ice-sheet models // Annals of Glaciology. 2003. V. 37. P. 150158.

4. Рыбак P.O.. Рыбак E.A. Спектральная структура изменений климата в позднем плейстоцене // Вестник Южного научного центра РАН. 2005. т.1. № 3. С. 51-65.

5. Рыбак О.О. Выделение климатического сигнала из изотопного ряда в антарктическом ледовом керне // Известия ВУЗов. Северо-Кавказский регион. Естественные науки. 2006. №4. С. 104-108.

6. EPICA, community members, Fischer H., Freitag J., Frenzel A., Fritzsche D., Fundel F., Gersonde R., Hamann I., Huybrechts P., Kipfstuhl S., Lambrecht A., Meyer H., Miller H., Oerter H., Ruth U., Rvbak P., Schmitt J., Valero-Delgado F., Wegner A., Wilhelms F. Pne-to-one interhemispheric coupling of polar climate variability during the last glacial // Nature. 2006. V. 444. P. 195-198

7. Рыбак O.O., Рыбак E.A. Спектральный анализ изотопных рядов в ледяных кернах и морских донных осадках // Материалы гляциологических исследований. 2006. вып. 101. С. 17-23.

8. Рыбак О.. Хёбрехтс Ф., Паттэн Ф., Штайнхаге Д. Региональная модель динамики льда. Часть 1. Описание модели, постановка численного эксперимента и современная динамика потока в окрестностях станции Конен // Материалы гляциологических исследований. 2007. Вып. 102. С. 3-11.

9. Рыбак О.. Хёбрехтс Ф.. Паттэн Ф., Штайнхаге Д. Региональная модель динамики льда. Часть 2. Пост-экспериментальная обработка данных // Материалы гляциологических исследований. 2007. Вып. 103. С. 3-10.

10. Рыбак P.O.. Хёбрехтс Ф. Лагранжев метод расчета возраста и изотопного состава льда в трехмерной модели антарктического ледникового щита // КриосфераЗемли. 2007. Т. 11. №4. С. 70-79.

11. Рыбак О.О. Спектральная структура изотопного состава антарктических ледовых кернов: сравнение оценок Ломба-Скаргла и максимальной энтропии // Известия ВУЗов. Северо-Кавказский регион. Естественные науки. 2007. №6, С. 50-56.

12. Ruth U., Barnola J.-M., Beer J., Bigler M., Blunier Т., Castellano E., Fischer H„ Fundel F., Huybrechts P., Kauimann P., Kipfstuhl S., Lambrecht A., Morganti A., Oerter H„ Parrenin F., Rvbak P.. Seven M., Udisti R„ Wilhelms F., Wolff E. "EDML1": A chronology for the EDML deep ice core from Dronning Maud Land, Antarctica, over the last 150000 years // Climate of the Past. 2007. V. 3. P. 475-484.

13. Huybrechts P., Rvbak P.. Pattyn F., Ruth U., Steinhage D. Ice thinning, upstream advection, and non-climatic biases of the upper 89% of the EDML ice core from a nested model of the Antarctic ice sheet // Climate of the Past. 2007. V. 3. P. 577-589.

14. Рыбак О.О. Математические модели континентальных ледниковых щитов. Часть 1. Архитектура моделей // Криосфера Земли. 2008. т. 12. №1. С. 12-23.

15. Рыбак О.О. Математические модели континентальных ледниковых щитов. Часть 2, Сравнительная характеристика // Криосфера Земли. 2008. т. 12. №3. С. 12-21.

16. Rvbak О.. Huybrechts P. Sensitivity of the EDML ice core chronology to geothermal heat flux // Материалы гляциологических исследований. 2008. Вып. 105. С. 35-40.

17. Рыбак О.О. Эйлеров и лагранжев методы восстановления возраста льда в простой численной модели // Вестник Южного научного центра РАН. 2008. №4. С. 38-45.

18. Nemec J., Huybrechts P., Rvbak P.. Oerlemans J. Reconstruction of the surface mass balance of Morteratschgletscher since 1865 // Annals of Glaciology. 2009. V. 50. P. 126-134.

19. Huybrechts P., Rvbak P.. Steinhage D., Pattyn F. Past and present accumulation rate reconstruction in the Eastern Dronning Maud Land, Antarctica // Annals of Glaciology. 2009. V. 51. P. 112-120.

20. Рыбак P.P. Моделирование миграции границы морского ледникового шита в простой численной модели // Известия ВУЗов. Северо-Кавказский регион. Естественные науки. 2010. №4. С. 126-130.

21. Рыбак P.O.. Рыбак Е.А. Алгоритм решения системы уравнений течения льда в трехмерной математической модели // Известия ВУЗов. СевероКавказский регион. Естественные науки. 2010. №6. С. 117-122.

22. Рыбак P.O. Миграция границы морского ледникового щита в численной модели // Вестник Южного научного центра РАН. 2010. №4. С. 50-56.

Нереферируемые журналы, сборники и материалы конференций

23. Rybak О. Spectra of Isotopic Records: Marine Sediments vs Ice Core // ASOF Newsletter. Issue 4. November 2005. P. 20-23.

24. Рыбак P.O.. Рыбак E.A. Стабильные изотопы в ледовых кернах и вариации температуры воздуха в прошлом / Материалы конф. "Проблемы устойчивого развития регионов Юга России". Сочи: СНИЦ РАН, 2006, с. 153-176.

25. Fürst J.J., Rybak P.. Goelzer H., De Smedt В., de Groen P., Huybrechts P. Improved convergence and stability properties in a three-dimensional higher-order ice sheet model // Geoscientific Model Development Discussions. 2011. V. 4. P. 1569-1610.

Свидетельства о государственной регистрации программ для ЭВМ Свидетельство №2011612257 от 17.03.2011 о государственной регистрации программы для ЭВМ homr.f. Правообладатель и автор Рыбак О.О. Свидетельство №2011612285 от 17.03.2011 о государственной регистрации программы для ЭВМ backtracinglsmr.f. Правообладатель и автор Рыбак P.P. Свидетельство №2011615780 от 22.07.2011 о государственной регистрации программы для ЭВМ gl2dr.f. Правообладатель и автор Рыбак P.P.

Формат 60*90/16 Бумага офсетная Гарнитура шрифта Тайме Усл. печ. л. 2,0. Тираж 100 экз. Заказ № 2944 Подписано в печать: 02,09.2011г.

Типография ИП Кривлякин С П. г. Сочи, ул. Новоселов, 5-44