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

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

ВИЛЬЧИНСКИЙ Александр Владиславович

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ МОРСКИХ ледников И ИХ ОСОБЫХ ЗОН

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

Автореферат диссертации па соискание ученой степени кандидата физико - математических наук

КАЗАНЬ - 1997

Работа выполнена на кафедре прикладной математики Казанского государственного университета.

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

профессор В.А. Чугунов.

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

доцент А.Б. Мазо, кандидат физико - математических наук, доцент М.Ф. Павлова.

Ведущая организация: Институт физико-технических проблем

Севера СО РАН, г. Якутск.

Защита состоится " г 9\иС2я|рчД 1998 г. в \ ^ часов на заседании диссертационного совета Д053.29.20 при Казанском государственном университете по адресу: 420008, г. Казань, ул. Университетская, 17, НИИММ, ауд. 324.

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

Автореферат разослан " 1997 г.

Ученый секретарь диссертационного совета, кандидат физико-математических наук, доцент —

Е.М. Федотов

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

Одним из направлений в исследовании естественных масс льда является изучение динамики ледников, взаимодействующих с морем. Данный аспект определяется тем, что крупнейшие ледниковые системы, такие как Антарктида, Гренландия и Земля Франца-Иосифа, являются морскими и оказывают значительное влияние на климат как ближайших материков, так и Земли в целом.

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

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

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

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

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

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

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

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

Апробация работы Основные научные результаты

диссертации докладывались на Международной конференции "International Symposium on Ice Sheet Modelling" (Chamonix MontBlanc, France, 18-22 сентября, 1995), Международном гляциологическом симпозиуме (Пущино, 22-26 мая, 1996),

Международной гляциологической конференции "Эволюция и прогноз развития полярного оледенения: новые результаты и модели"(Москва, 19-22 мая, 1997), 2 Международной конференции по математическому моделированию (Якутск, 28 июня - 2 июля, 1997), летней школе 'Modelling of Glaciers and lee Sheets' (Karthaus-Certosa, Italy, 2-13 сентября, 1997), итоговых научных конференциях Казанского государственного университета 1993-1996, семинарах кафедры прикладной математики Казанского государственного университета.

Публикации. По основным результатам диссертации опубликовано 5 работ и одна принята к опубликованию.

Структура и объем работы. Диссертация состоит из введения, четырех' глав, заключения и списка литературы из 103 наименований. Общий объем работы 156 страниц.

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

Во введении обосновывается актуальность темы, дается обзор литературы и современного состояния проблемы. Ставится цель работы и приводятся основные результаты, выносимые на зашиту. Описывается содержание диссертации.

Глава I посвящена построению упрощенной модели плоского изотермического нестационарного течения морского ледника и определению граничных условий на его линии налегания.

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

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

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

Ледниковый шельф характеризуется своими собственными пространственными масштабами. Для исследования его динамики необходимо перейти к безразмерным координатам, нормированным на данные масштабы. Несмотря на то, что задача течения шельфа была решена и ранее при помощи приближения тонкого слоя, масштабы его длины и толщины рассматривались как независимые, что, вообще говоря, неверно. Поэтому данная задача рассматривается в §1.3 с целью построения модели, включающей и определение связи между данными величинами. Масштаб толщины шельфа [г] находится из анализа уравнений и зависит от масштаба его длины [х]:

где [я] - масштаб потока массы, г) - реологический параметр, а -коэффициент ползучести, плотности льда и воды,

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

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

Устанавливается, что характерный наклон верхней поверхности в переходной зоне пропорционален нормированной разности между плотностями воды и льда 5=0.1. Анализ уравнений

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

L(1 - 5) = D + O(S) + О(е), р = ôL(a+2)/a / eQ1/a. Здесь L, D - толщина льда и глубина воды на линии налегания, с - характерный наклон поверхности наземного ледника, р -параметр, определяемый решением внутренней задачи, Q - поток массы через сечение, положение которой совпадает в данный момент с линией налегания.

После разложения функций задачи в асимптотические ряды по малым параметрам s и 5 можно получить более простую (относительно формы области и граничных условий) задачу для главных членов внутреннего разложения функции тока у и функции толщины льда h в переходной зоне:

( -Г>

ду~

1 \

d~\\i д'\\1 \ду2 дЕ,

+ 4

З2

д2ш З2 I 3~ш д~и>

дудЪ, дуд^ 1-а

= О,

ду-

ц =

32у

W

Л

О \|/

/2 1 \1 д у

ду1 д^

О, \|/ = 1 при у = 1;

, - од < Ç < +«, 1 - h < у < 1,

ц/ = о, — = 0 при у = 0. £, < 0 ; (1) ду

д ¿ у дуг 3Ç2

(-h").

-4h

д^ду

,v = О при у = l-h(Ç), Ç > 0;

j "у

р—+ h —= f Bdy при Ç > 0 ; h— 1 при Ç < 0 ; 7 п? J

dï i

В = 4ц

3~1[/ 3 f [ д V|/ 3~\|;

яр i

ô^Sy 3Ç;'Азу2 3£,

l-h 2,

dy,T = {

i ( =2 л

3 у 3 \|/

i-h ^дУ~

dy.

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

Процедура определения равномерно пригодного разложения для функции профиля поверхности ледника описана в §1.5.

Общий алгоритм решения задачи динамики ледника при плоском течении льда описан в §1.6 и представляет собой совместное решение задач переходной зоны, наземной части ледника и шельфа.

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

= 0 при 0 < ^ < +оо, 0 < у < 1, д2Ч> дЧ1

—- = = ^ при у = 1; ¥ = 0,-— = 0 при у = 0; (2)

ду~ 8у

эй* а2 у

- = 0,^-4=0 при 5 = 0.

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

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

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

Анализ задач плоского и осесимметричного течения льда сходен по своей структуре, однако осесимметричное течение

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

В §2.1 описана постановка задачи осесимметричного течения. После обезразмеривания параметров задачи на масштабы, характерные для' наземной части ледника, уравнения, записанные в безразмерном виде, включают малый параметр при старшей производной. Решение, описывающее динамику наземного ледника при осесимметричном течении, было получено ранее другими авторами. Его нахождение методами возмущений описано в §2.2. Исследование динамики шельфа при осесимметричном течении, который ранее не приводился, позволяет найти масштабы шельфа и уравнение для определения толщины льда (§2.3).

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

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

Моделирование стационарного осесимметричного течения льда в области ледораздела представлено в §2.7. Полученные результаты аналогичны найденным при исследовании плоского течения: при моделировании течения льда в области ледораздела уравнения Стокса в цилиндрической системе координат должны

решаться в бесконечном слое с горизонтальными параллельными границами.

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

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

Анализ стационарного течения льда с постоянной вязкостью вблизи линии налегания представлен в §3.1. Рассматриваемая математическая модель включает бигармоническое уравнение для функции тока в бесконечной полосе с неизвестной нижней границей и уравнение для определения данной границы (1). При численном моделировании полоса усекается с обоих концов. Дифференциальная задача четвертого порядка с однородными граничными условиями для функции тока сводится к проблеме минимизации функционала

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

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

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

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

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

В §3.2 моделируется течение льда со степенным реологическим законом в окрестности ледораздела, описываемое уравнениями (2). После усечения исходной области определения функции тока -

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

1+а

Ф(чО = I

1 а

п1 + а

5 VII д у

дг~ дх2

( -,->

+ 4

2а

сШ

^дхдг)

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

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

Глава IV посвящена моделированию стационарного течения ледников Земли Вильчека, которые взаимодействуют с морем и характеризуются наличием точки перегиба функции профиля верхней поверхности, что не характерно для большинства ледников. Цель исследования - установить факторы, определяющие данную особенность течения, положение линии налегания, а также найти неизвестные параметры граничного условия на линии налегания по реальными данными. Анализ распределения функции высоты поверхности льда позволяет установить, что у ледников Земли Вильчека нет шельфов, так как условие гидростатического равновесия льда в воде достигается только на их линиях скола. В §4.1 приводится математическая

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

Метод численного решения задачи представлен в §4.2. Разностная схема строится методом сумматорных тождеств. Полученная система уравнений решается методом Ньютона.

Численные эксперименты, описанные в §4.3, показывают, что в областях, где ложе находится ниже уровня моря, скорость скольжения сильно возрастает. Это может быть вызвано проникновением морской воды между льдом и ложем и образованием там полостей, наполненных водой, что приводит к существенному увеличению значения коэффициента скольжения. Данный факт определяет существование точки перегиба функции профиля поверхности ледников. Идентификация неизвестных параметров граничного условия на линии налегания ледника подтвердила оценки их характерных величин, полученные в главах 1,11.

В заключении подводится итог выполненному исследованию.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ

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

2. Определены граничные условия на линии налегания и влияние на них различных типов течения: плоского, осесимметричного, стационарного и нестационарного. Показано, что при плоском течении с достаточной степенью точности задача для переходной зоны является стационарной. При осесимметричном течении это не всегда верно.

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

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

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

Автор выражает глубокую признательность научному руководителю профессору В.А. Чугунову за поддержку при выполнении данных исследований, профессору Саламатину А.Н. за обсуждение аспектов механики льда, сотрудникам Института Географии РАН Глазовскому А.Ф. и Мачерету Ю.Я. за предоставленные материалы по ледникам Земли Вильчека (Земля Франца-Иосифа), а также всем участникам семинара кафедры

фикладной математики Казанского государственного университета за полезные замечания.

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

1. Chugunov, V.A., and A.V. Wilchinsky, Modelling of a marine glacier and ice-sheet - ice shelf transition zone based on asymptotic analysis // Annals of Glaciology, 1996, Vol. 23, p.59-67.

2. Вильчинский A.B., Чугунов B.A.. Моделирование динамики морского ледника и его особых зон // Материалы Гляциологических Исследований, 1997, Вып.83, с. 92-97.

3. Wilchinsky A.V., Chugunov V.A. Modelling ice-divide dynamics by perturbation methods // Journal of Glaciology, 1997, Vol. 43, №114 , p.352-358.

4. Вильчинский A.B., Чугунов B.A. Математический анализ условий течения ледников Земли Вильчека/ Казан, ун-т. - Казань, 1997. - 13 е.: 3 ил. - Библ. 11 назв. - Рус. - Деп. в ВИНИТИ № 2076-В97.

5. Вильчинский А.В., Чугунов В.А. Анализ динамики ледораздела при осесимметричном течении льда/ Вильчинский А.В., Чугунов В.А.; Казан, ун-т. - Казань, 1997. - 16 с.: ил. - Библ. 9 назв. - Рус. - Деп. в ВИНИТИ №2078-В97.

6. Вильчинский А.В, Чугунов В.А. Построение теоретической модели стационарного морского ледника с законом течения Глена // ПММ (принято к опубликованию 23.10.1996).