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

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

Червов Виктор Васильевич

ТРЕХМЕРНОЕ МОДЕЛИРОВАНИЕ КОНВЕКТИВНЫХ ПРОЦЕССОВ В МАНТИИ ЗЕМЛИ

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

АВТОРЕФЕРАТ

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

Новосибирск 2004

Работа выполнена в Институте Геологии СО РАН и в Институте Вычислительных Технологий СО РАН (г. Новосибирск)

Научные руководители:

доктор геолого-минералогических наук

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

С.А. Тычков Г.Г. Черных

Официальные оппоненты:

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

А.Ф. Воеводин

кандидат геолого-минералогических наук, старший научный сотрудник

О.П. Полянский

Ведущая организация - Институт Теплофизики им. С.С. Кутателадзе СО РАН (г. Новосибирск)

Защита диссертации состоится 03 марта 2004 года в 14 часов 00 минут на заседании диссертационного совета Д 003.061.02 при Институте Вычислительной Математики и Математической Геофизики СО РАН по адресу: 630090, г. Новосибирск, проспект академика Лаврентьева, 6.

С диссертацией можно ознакомиться в читальном зале библиотеки ИВМ и МГ СО РАН, проспект академика Лаврентьева, 6.

Автореферат разослан 31 января 2004 г.

Ученый секретарь диссертационного совета -

д.ф.-м.н С.Б. Сорокин

aooi-b

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

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

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

При решении трехмерных задач гидродинамики достаточно хорошо зарекомендовали себя переменные «завихренность - векторный потенциал»4,10. Однако при численном моделировании конвективных процессов в мантии Земли указанным переменным уделено недостаточное внимание. Так, например, в статье5, выполнено решение только для постоянной вязкости. В настоящей работе построена численная модель конвективных процессов в мантии Земли, основанная на вышеупомянутых переменных и методе дробных шагов6.

Актуальность темы.

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

1 McKcnzic D.P., Roberts J.M., Weiss N.O, Convection in Earth's mantle: towards a numerical simulation // J. Fluid Mech.; 1974; v.62, part 3, p.465-538.

2 Рыков B.B., Трубицын В.П. Численное моделирование трехмерной мантийной конвекции и тектоника литосферных плит//Вычисл. сейсмология. 1994. Вып. 26. С.94-102.

3 Busse F.H., Christensen U., Clever R., Csercpes L., Gable C., Giannandrea E., Guillou L., Houseman G., Nataf H.-C., Ogawa M., Parmentier M., Sotin C., Travis B. 3D Convection at infinite Prandl number in cartesian geometry - a benchmark comparison.// Geophiys. Astro-phys. Fluid Dynamics. 1993. Vol.75. P.39-59.

4 Роуч X. Вычислительная гидродинамика. M.: Мир, 1980.

5 Sparks D. W., Parmentier E.M, Morgan J.P. Three-dimensional mantle convection beneath spreading centers spreading centers Implications for along- axis variations in crustal thickness and gravity // Journal of Geophysical research; 1993; v.98, No.B12, p.21,977-21,995, December 10.

6 Япенко H.H. Метод дробных шагов решения многомерных задач математической физики. «Наука», СО, Новосибирск, 1967.

ным представлениям, Земля является «тепловой машиной», в которой радиогенное тепло превращается в механические движения вещества ядра и мантии. На поверхности планеты эти движения выражаются в создании горных хребтов, горизонтальных перемещениях материков на тысячи километров и формировании глубоких котловин. Изучение структуры недр возможно только дистанционными методами с помощью геофизических наблюдений. Основным же инструментом в создании динамических эволюционных моделей недр планеты являются математические модели. Первые современные двухмерные численные модели динамики мантии появились тридцать лет назад. В настоящее время, с развитием вычислительных технологий и вычислительной техники, появилась возможность создания трехмерных моделей движений в мантии. Таким образом, в настоящее время весьма актуальной проблемой является создание трехмерных моделей динамики недр конкретных геологических структур с реальной трехмерной геометрией и более сложной реологией недр, что обеспечивает переход моделей на качественно новый уровень достоверности. Цель работы состоит в разработке трехмерной численной модели тепловой конвекции в верхней мантии Земли, основанной на переменных у/,а>, и исследовании особенностей динамики мантии континентальных областей Земли. Научная новизна изложенных в диссертации результатов заключается в следующем:

> построение и детальное тестирование трехмерной численной модели конвекции в верхней мантии Земли;

> создание трехмерных геодинамических моделей тепловой конвекции в верхней мантии Земли под континентальной литосферой переменной толщины;

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

> трехмерная численная модель конвекции в верхней мантии Земли основанная на уравнениях Навье -Стокса в приближении Обербека -Буссинеска и геодинамическом приближении;

> результаты численных экспериментов по созданию трехмерной модели динамики верхней мантии платформенных областей Центральной Азии;

Апробация работы. Основные результаты диссертации докладывались на международных конференциях: совещании по геодинамической эволюции палеоазиатского океана (Beijing, 1991), мемориальной конференции памяти Л.П. Зо-неншайна по Тектонике Плит. (Москва, 1993) и международной конференции «Потоки и структуры в жидкостях» (С.- Петербург, 2003). Полное содержание диссертации докладывалось на семинарах Института Геологии СО РАН (рук.

д.г.-м.н. А.Ю. Казанский), Института Гидродинамики им. М.А. Лаврентьева СО РАН (рук. проф. А.Ф. Воеводин), Института Вычислительных Технологий СО РАН (рук. проф. В.М. Ковеня), Института Вычислительной математики и математической геофизики СО РАН (рук. академик A.C. Алексеев; чл.-корр., РАН Б.Г. Михайленко).

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

Представленные в диссертации результаты получены при проведении исследований в рамках Приоритетного направления СО РАН «Геодинамическая и геохимическая эволюция литосферы и мантии Земли: тектоника, магматизм, флюидный режим и металлогения», по Программе СО РАН №26.2. «Геодинамическая эволюция литосферы Азии: тектоника, магматизм, метаморфизм, геохимия и металлогения основных ее этапов»

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

Краткое содержание диссертации.

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

1

В Главе 1 приводится постановка задачи. Для описания течения в верхней мантии Земли привлекается хорошо известная математическая модель, включающая в себя обезразмеренные уравнения7:

7 Добрецов Н.Л., Кирдяшкнн А.Г., Кнрдяшкин А.А. Глубинная геодинамика, 2-е издание, Новосибирск, Изд-во СО РАН, 409 е., 2001.

(1.2)

(1.1)

Эу Эх

г да Эу> ,Эу Эх,

+ 2—77— +—?7

ду ду дх

'Эу | Эw Эу /

(1.3)

Эр д (да Э\у^ Э ^ Эу Эи^

— = —77 — + — +—77 — + —

Эг Эх ^Эг Эху Эу ^Эг Эу

ЭТ ЭТ ЭТ ЭТ _2_ _

— + и— + V— + \у— = У2Т + р. д1 Эх Эу дг

+ (1.4)

дъ дг

(1.5)

Здесь и, V, - компоненты вектора скорости, р - давление, Т - температура, О - источник тепловыделения, 11а - число Рэлея, 77 - динамическая вязкость.

Система уравнений (1.1)-(1.5) устроена так, что в начальный момент времени г = ¡0 задаются начальные условия лишь для температуры:

Т(х^0)= Т0(х,у,г) •

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

Для построения численной модели осуществляется переход к новым независимым переменным: вектору завихренности

сд= 1й)х + )ах' +к«иг = УхУ (1.6)

и векторному потенциалу у/ = \у/* +)у/'/ + к.у/г, У = Ух^.

В результате система уравнений (1.1 )-(1.4) переходит в следующую:

где

V V = ~а>

У2[т7<ах] = тгта>х + (77Х + 77 усоух + 17ж®хж)+ Р1 " КаТу> У2[77й)у] = т]уусоу + {т1хсо;+т1усо1 + 77,©;)+ Р2 + ЫаТх,

У2[т]сог) = + (+ + )+

=2(77^ -77ууУг)+2т7уг(уу -ид-77ху(и2 +^ух)+77х2(ух -иу); =2С7яЦ. -7Л)+2/7хгК -их)-77уг(ух +иу)+?/.,у(\уу -уД Рз=2(/?^х -?/ххиу)+277ху(их -Уу)-7?«^ +У2)+^(иг -wx)..

(1.7)

(1.8)

(1.9)

(1.10) (1.11)

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

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

1) При известном распределении температуры вычисляются завихренность, векторный потенциал и компоненты скорости (до сходимости итерационных процессов);

2) вычисляется поле температуры по схеме стабилизирующей поправки:

т-в+1/3_гг«

+и1^Г+1/3-ЬххГ+,/3 =-уЦТ ч-Г^Г-и^Г +ЬаГ

'I

- + уЬуТ5+2/3 - Ь>уГ+2/3 = VЬуТ5 - ЦГ

рв+г/з

П5+1

+ \\гЬ,Т - Ь„Т5+1 = wL .Г - ЬТ5

где г, - величина шага по времени. В схеме трехточечные разностные операторы ЬХХ,Ь ,ЬЯ и Ьх,Ьу,Ьг аппроксимируют дифференциальные следующим образом:

Зх2 хх Ь2 ах"ЧГ~ 2ЬХ '

££^Т + и

дуг-^- Ь2 ау-ЦГ" 2Ьу '

э2? ~ г ^ _ 4-и - + 4-1 и _ 4+< - 4-1 • = ъ1 = 2ь

2

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

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

В качестве примера тестирования численной модели ниже представлены результаты решения модельной задачи3. Решение отыскивалось в единичном кубе. При этом задавались следующие параметры: масштабный множитель при вязкости:

% = 1.20165 -1024; 7?(т)= ехр[#/(Т + ©)~ 0/(0.5 + 0)]; в = [225/ 1п(г)] - 0.251п(г); ©=15/1п(г) -0.5,

Число Рэлея:

Яа = Ь.'Х Iг)йх = 20 000 •

Вычислялись следующие параметры: (¡) среднеквадратичная скорость

где А - объем параллелепипеда со сторонами X, У и

(и) число Нуссельта (Ыи) по формуле

Ми = -(Х-У)"1 [Г—(Ыу> где ^„р" веРхняя поверхность параллелепипеда;

(¡и) значение вертикальной компоненты скорости и температуры Т в угловых точках среднего сечения конвективного слоя;

(¡V) значение теплового потока 19 = _5Т/ в угловых точках верхней по/ дг

верхности куба;

(V) интегральный параметр, вычисляемый по формуле

Т(Х,2) = (—<1у ¿дг

вдоль линии, параллельной оси У, начинающейся точками (0,1/4), (1/2,1/4), (1,1/4) фронтальной (Хг)-плоскости;

(уО средняя температура у _ вычисляемая на горизонтальных

сечениях области Зг_0 м и 8г_0 50 , на глубинах ъ = 3/4 и ъ = 1/2;

(ун) значение вертикальной компоненты вектора завихренности сог в точке (0.75,0.25,0.75).

Размерные значения (в системе СИ), которые были использованы в настоящей работе, принимались следующими:

й = 2 700 000 м, АТ = 3700 °С, % = Ю^м2 /с = 10 -'"С"1, р = 3300 кг / м3, §2 = 10 м /с2 , 770 = 10 24 Па -с.

Краевые и начальные условия для завихренности и векторного

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

Т(х,у,2,1:0)= Т0(х,у,2) = (1 - г) + 0.2(соз(ях /X) + соз(яу /У))Бт(яг) • Начальные значения векторных функций завихренности и потенциала - нулевые.

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

—; ах'

ди ~ду'

Граничные условия симметрии для функции со,у/

> на поверхностях х = 0; х = Х, 0<у<У, 0йг<1'-

др* V 2 /Ч V л г

-X— = у/* = у/1 = 0 > -= й)У = 0, й) = •

Эх Эх

> на поверхностях у = 0;у = У, О^х^Х, 0<г<1:

ЭуУ X г Л дО)У . Л г

—£— = Щ-Х =у/г =0> --ф = 0, СО ■=-

ду ду

Условия прилипания для функций а>,у/~.

> на поверхностях х = 0; х = X, 0 < у < У,

х ди/у ч ду/г 2 Л л „ ды , ду .

дх дх дх Эх

> на поверхностях у = 0;у = У. 0 < х ^ X, 0^251:

ду/* , у ду/г 2 . дч/ у . , Эи.

ду ду ду ду

> на поверхностях г = 0; г = 1, 0 < х < X, 0 < у У:

х у г л х у да , .

гага гага

Результаты расчетов сопоставляются с данными Кристенсена, как наиболее полными из имеющихся в статье3, и представлены в Табл. 1,2. Ошибка вычислялась по формуле:

Err =

Che - Chr

100%- (L1)

Chr

Достаточно хорошо известным эффективным подходом к решению задач математической физики является метод последовательности сеток8. Результаты его применения в настоящей работе иллюстрируются Табл. 2. При непосредственных вычислениях на сетке 48x48x48 потребовалось 660 минут и 24 секунды на ПК с процессором Athlon 1000. Таким образом, выигрыш во времени счета на последовательности сеток более чем восьмикратный. Существенный выигрыш во времени счета может быть достигнут также с применением экстраполяции по Ричардсону9.

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

8 Фсдорснко Р.П. О скорости сходимости одного итерационного процесса // Журн. Вы-числ. Математики и мат. Физики. 1964. Т.4. № 3. С.559-564.

' Марчук Г.И., Шайдуроц В.В. Повышение точности решения разностных схем. М., «Наука», 320 стр., 1979.

Таблица 1

Наименование Результаты Результаты Относительная

параметров Кристенсена автора (Che) ошибка(Егг)

(СЬг) на сетке на сетке 32x32x32 в%

32x32x64

№ 3,03927 3,040 0,0240

УГГПБ 35,132 35,18 0,1366

-58,23 -58,23 0,0000

Т( 1,1,0.5) 0,23925 0,2379 0,5643

5(1,1) 0,7684 0,7731 0,6117

х(1,1) -0,1388 -0,1356 2,3055

Тт(0.50) 0,58158 0,5799 0,2889

СО1 (0.75,0.25,0.75) -11,125 -11,25 1,1236

Таблица 2

Наименование параметров Результаты Кристенсена (СЬг) на сетке 32x32x64 Результаты автора (СЬе) на сетках

12x12x12 24x24x24 48x48x48

N11 3,03927 3,245 3,0397 3,040

Угшв 35,132 37,97 35,360 35,07

№(1,1,0.5) -58,23 -60,68 -58,262 -58,42

Т( 1,1,0.5) 0,23925 0,2326 0,2377 0,2393

5(1,1) 0,7684 0,80172 0,7804 0,7726

т( 1,0.25) -0,1388 -0,06761 -0,1341 -0,140

Тга(0.50) 0,58158 0,5844 0,5792 0,5815

Сй! (0.75,0.25,0.75) -11,125 -11,03 -11,203 -11,35

Время счета 65 сек. 119 сек. 4811 сек.

Отдельные этапы алгоритмов тестировались на задачах, имеющих точные аналитические решения.

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

10 Андерсон Д., Таннехнлл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. М.: Мир, 1990.

X = X, у = у,

где В = —In 2 Я

z = B + — arsh Я

l + (eA-l)(ze/Z)

-1 -sh(AB)

0<Я<оо.

_1 + (е-*-1)(2с/г)_ Значения параметра растяжения Я влияют на измельчение сетки вблизи г = гс • Растяжение применялось по направлению Ъ для модельной задачи. Для значений гс = 1, Z = 1 и Я = 2 получены результаты, представленные в Табл. 3. Для значений г =3/4, 2 = 1 и Я=4 -в Табл.4.

Таблица 3. Переменная вязкость. Преобразованные уравнения

Наименование Результаты Результаты Относитель-

параметров Кристенсена автора(СЬе) ная

(СЬг) на сетке ошибка

на сетке 32x32x64 32x32x32 (Err) в %

№ 3,03927 3,0890 1,607

Ушв 35,1320 35,800 1,868

М1,1,0.5) -58,230 -60,840 4,285

Т( 1,1,0.5) 0,23925 0,2431 1,596

5(1,1) 0,76840 0,7862 2,271

Тт(0.50) 0,58158 0,5912 1,631

фг (0.75,0.25,0.75) -11,1250 -11,520 3,446

Таблица 4. Переменная вязкость. Преобразованные уравнения

Наименование параметров Результаты Кристенсена (СЬг) на сетке 32x32x64 Результаты автора (Che) на сетке 32x32x64 Относительная ошибка (Err) в %

№ 3,03927 3,0910 1,684

УГШБ 35,1320 35,710 1,631

м1,1,0.5) -58,230 -60,340 3,497

Т(1,1,0.5) 0,23925 0,2415 0,919

>9(1.1) 0,76840 0,7816 1,689

Тт(0.50) 0,58158 0,5890 1,263

йГ (0.75,0.25,0.75) -11,1250 -11,460 2,955

Приведенные выше результаты численных экспериментов свидетельствуют о достаточно высокой эффективности построенной численной модели.

В Главе 3 представлены результаты трехмерного численного модели-

рования процесса тепловой конвекции в верхней мантии Земли под континентальной литосферой переменной толщины. Последовательно рассмотрены следующие пять задач эволюции мантийного вещества в параллелепипеде (4200*4200x700 км3):

1) Тепловая конвекция под литосферной плитой постоянной толщины;

2) Тепловая конвекция под литосферной плитой с утолщенной полосой в ее

центральной части;

3) Тепловая конвекция под литосферной плитой с квадратным, в плане

утолщением в ее центральной части (кратоном);

4) Тепловая конвекция под литосферной плитой, содержащей два утолщен-

ных массива;

5) Тепловая конвекция под литосферной плитой с утонением в ее центральной

части.

Число Рэлея, характеризующее режим конвекции, было выбрано как Ла = 3.04x105, что отвечает современным представлениям об условиях в недрах Земли. Основные параметры задачи в системе СИ, для верхней мантии, принимались следующими:

<1 = 700 000м, ДТ = 1800°С,лг = 10"6м2/с, а = Ю-5 = 3300кг= Юм /с2, ?]0 = 3 -Ю^Па -с.

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

г/(х, у, г, 0 = ехр(Ьг - аТ(х, у, г, 1)). Здесь параметры а = 3.89 и Ь = 5.84 обеспечивают перепад вязкости от 20 до 200, что присуще верхнемантийным характеристикам течений:

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

лои:

Все задачи данной главы были рассчитаны на равномерных сетках 66x66*35 (размер ячейки: 63,6x63,6x20 км3). Задача с одним и с двумя кратона-ми была пересчитана на сетке со сгущением на глубине 300 км с целью выявления более тонких структур под кратоном/кратонами. Полученные решения показали близость результатов на равномерной и неравномерной сетках, что позволило сделать вывод о достаточности разбиения расчетной области.

Для оценки погрешности вычислений была посчитана задача тепловой конвекции под литосферой постоянной толщины на сетке вдвое более мелкой: 132x132x70. Здесь применялись как последовательность сеток, так и экстраполяция по Р.ичардсону.

На Рис. 1, Рис. 2 изображены для качественного сравнения (без указания значений изолиний) горизонтальные срезы на глубинах 220, 260,420 и 460 км |

Рис. 1. Распределение температуры в расчетной области на различных горизонтальных уровнях. Сетка 132x132x70, верхний левый квадрат-срез на глубине 460 км, верхний правый - срез на глубине 420 км, нижний левый - срез на глубине 260 км и нижний правый - срез на глубине 220 км. Изолинии температуры ' изображены с шагом 150 °С.

I

Сравнение срезов Рис. 1 с соответствующими срезами Рис. 2 показывают, что структуры теплового поля на существенно различающихся сетках достаточно близки.

Рис. 2. Распределение температуры в расчетной области на различных горизонтальных уровнях. Сетка 66*66x35, верхний левый квадрат - срез на глубине 460 км, верхний правый - срез на глубине 420 км, нижний левый - срез на глубине 260 км и нижний правый - срез на глубине 220 км. Изолинии температуры изображены с шагом 150 °С.

Кроме того, производился расчет на сетке 132x132x70 для неустановившегося течения. Была применена процедура экстраполяции по Ричардсону на двух сетках для получения решения, соответствующего решению на сетке 264x264x140. На Рис. 3 - Рис. 5 изображены для качественного сравнения (без указания значений изолиний) горизонтальные срезы на глубине 460 км

течения. Время I = 1375 млн лет. Горизонтальный срез на глубине 460, рассчитанный на сетке 66*66x35. Изолинии температуры изображены с шагом 100 °С.

Рис. 4. Распределение температуры в расчетной области для неустановившегося течения. Время I = 1375 млн лет. Горизонтальный срез на глубине 460, рассчитанный на сетке 132x132x70. Изолинии температуры изображены с шагом 100 °С.

Рис. 5. Распределение температуры в расчетной области для неустановившегося течения. Время I = 1375 млн лет. Горизонтальный срез на глубине 460, полученный в результате экстраполяции по Ричардсону на сетках 66x66*35 и 132x132x70. Изолинии температуры изображены с шагом 100 °С.

Сравнение Рис. 4, 5 позволяет сделать вывод о близости полученных решений.

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

платформы, средний возраст которых составляет 2000 млн лет". Типичным примером подобных структур в России является Сибирская платформа, расположенная восточнее р. Енисей и простирающаяся в этом направлении до гор Верхоянья, у подножия которых течет р. Лена. На юге эта платформа ограничена оз. Байкал, а на севере - Енисей-Хатангской низменностью. Характерные размеры платформы составляют около 2000 км по широте и долготе. Именно такой размер кратона был использован в этой модели, причем расположение его было смещено относительно центра верхней плоскости на 700 км по оси у и на 300 км по оси х. Это было сделано для того, чтобы избежать возможного появления ложных структур конвективного течения из-за симметрии.

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

Рис. б. Расположение горизонтальных сечений (ХУ) на глубинах 160-(а), 220-(Ь), 380-(с) и 560 км-(с1), в модели конвекции под кратоном, при I = 0 млн. лет. Мощность литосферы Ь = 120 км, мощность кратона - 200 км.

" Condie K.C. Episodic continental growth models: afterthoughts and extensions // Tectono-phys.; 2000; v.322, p. 153-162.

С глубиной эти изолированные потоки сливались в один с одновременным заметным уменьшением по толщине, что хорошо видно на Рис. 7Ь. Температура восходящих потоков была выше своего среднего значения на данном уровне глубины на 300° С.

а b

Рис. 7. Сечения в плоскости (XY) на глубинах 220 км -(а) и 380 км -(Ь) в модели конвекции под кратоном при t = 0 млн. лет.

Характерной чертой пространственно-временной эволюции конвективных течений под кратоном является формирование выраженного нисходящего потока конвекции под центральной частью кратона, вытянутого в меридиональном направлении с температурой 850° - 1000° С (Рис. 8а). Этот поток простирается до глубин 450-500 км (Рис. 8с). Следует отметить, что температура в нем на 100° С выше, чем в подобных нисходящих потоках вне кратона. Более высокая температура на ту же величину отмечена и в восходящих потоках. В целом особенностью конвекции под кратоном является локализация ячеек по глубине: если вне кратона вертикальный размер ячеек оценивается в 500 км, то под кратоном он менее 400 км.

Стационарное течение достигается в модели по прошествии 2000-2500 млн лет. Из сравнения этой структуры течения с начальной на малых глубинах, непосредственно под кратоном, видно, что теперь существуют два основных восходящих потока, разделенных вытянутым нисходящим в его центральной части (Рис. 7а. и Рис. 8а).

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

Главное отличие состоит в том, что под кратоном сохраняются как восходящие, так и нисходящие потоки, эволюционирующие во времени. В двухмерном же случае под кратоном формировался только восходящий поток.

Столь принципиальное отличие можно объяснить трехмерной природой конвекции при числах Рэлея более 105, что наглядно было продемонстрировано в модели с кратоном при числе Рэлея 300 ООО.

300 1000 1900 2000 2500 3000 3500 4000км

0 500 1000 1500 2000 2500 3000 3500 4000 км

(Ь)

-600 км

0 500 1000 1500 2000 2500 3000 3500 4000 км

-200-

-400-4,1

-600км И—=-r~--1 -г ....

О 500 1000 1500 2000 2500 3000 3500 4000 кы

Рис. 8. Сечения в плоскости (XY) на глубине 220 км (а) и на глубине 380 км (Ь), и сечения в плоскости (XZ): (с) в модели конвекции под литосферой с кратоном при t = 2500 млн. лет. Мощность литосферы L = 120 км, мощность кратона - 200 км. Сечение с(а) взято непосредственно перед кратоном (540 км),

сечение с(с) — за кратоном (2660 км), сечение с(Ь) пересекает кратон в срединной его части (1600 км)

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

'Рис. 9. Горизонтальный срез поля температуры на глубине 175 км, рассчитанный на сетке 72x72x40. Положение кратона выбиралось таким образом, чтобы, по возможности, восходящие потоки находились на периферии кратона. Левый нижний угол кратона имеет координаты (1100,1600) (км), правый верхний - (2900, 3600) (км).

Задача была решена в боксе 4200 х 4200 х 700 км3 на сетке 72x72x40. В начале выполнялся расчет до установления под однородной литосферной плитой, мощность которой составляла 105 км. Кратон, мощностью 192 км, был помещен в область нисходящих потоков. Длина кратона вдоль оси X составила 1800 км, вдоль оси У - 2000 км. Через 2000 млн. лет восходящие потоки на периферии кратона сблизились. К примеру, на глубине 245 км расстояния между потоками сократились на 315 км (верхняя пара, соединенная двойной стрелкой) и на 630 км (нижняя пара). В среднем - 472,5 км.

Рис. 10. Горизонтальные срезы поля температуры на глубине 245 км, рассчитанные на сетке 72x72x40 для времен 5 млн. лет (а), и 2000 млн. лет (Ь). Центры восходящих потоков, соединенных двойными стрелками сблизились за 1995 млн. лет в среднем на 472.5 км

В заключении перечислены основные результаты, полученные в диссертации. Они сводятся к следующему:

1. Разработана численная модель трехмерной конвекции в мантии Земли, основанная на переменных «векторный потенциал -завихренность» и методе дробных шагов. Осуществлено детальное тестирование построенной численной модели.

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

3. Принципиальный результат трехмерного моделирования тепловой конвекции в верхней мантии под кратоном состоит в выявлении мелкомасштабной моды конвекции непосредственно под литосферой на «астеносферном» уровне глубин 200-350 км Данная мода развивается по периферии утолщенного участка литосферы (кратона) и может объяснить особенности режима траппо-вого магматизма древних кратонов по их периферии.

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

Публикации по теме диссертации

Tychkov S.A., Kulakov I.Yu., Chervov V.V., Mantle structure and dynamics in Baikal and surround areas // IGCP Project No.2S3: Geodynamic evolution ofthePaleoasian ocean: Report No.2- Beijing, 1991, p.67. Tychkov S.A., Chervov V.V., Thinning of the lithosphere by upper mantle convection // L.P. Zonenshain memorial conference on plate tectonics: Abstr. - Moscow, 1993, p. 149.

Червов В. В. Сравнение некоторых аппроксимаций конвективных членов в уравнении теплопереноса для задачи конвекции в мантии Земли. // Сборник научных трудов ИВТ СО РАН «Вычислительные технологии»; 1995; т.4, № 13, с.295-305.

Тычков С.А., Рычкова Е.В., Василевский А.Н., Червов В.В. Тепловая конвекция в верхней мантии континентов и ее эффект в геофизических полях // Геология и геофизика; 1999; т.40, №9, с. 1275-1290. Червов В. В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением завихренности и векторного потенциала // Вычислительные технологии; 2002; т.7, №1, с. 114-125. Червов В. В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением последовательности сеток. // Вычислительные технологии; 2002; т.7, № 3. С.85-92.

Тычков С.А., Червов В.В., Черных Г.Г. Численная модель трехмерной конвекции в верхней мантии Земли // Тез. докл. международной 1 конф. «Потоки и структуры в жидкостях». (23-26 июня 2003 г., С-Петербург). - 2003.-Москва: Институт проблем механики РАН.- С.215.

Технический редактор О.М. Вараксина

Подписано к печати 26.01.2004 Формат 60x84/16. Бумага офсетЫ 1. Гарнитура Тайме. Офсетная печать. Печ. л. 1,2. Тираж 100. Заказ 15.

Издательство СО РАН. 630090 Новосибирск, Морской пр. 2 Филиал "Гео". 630090 Новосибирск, пр. Ак. Коптюга. 3

РНБ Русский ф

2007-4

17078

ФЕВ 200^

JZ&JSJZ&ttJ&JE;.

Глава 1 Постановка задачи трехмерного моделирования

• конвекции в мантии Земли.:.

1.1 Основные уравнения.

1.2 Основные уравнения Движения несжимаемой жидкости.

1.3 Обезразмеривание.

1.4 • Граничные условия.

1.5 Начальные данные.

1.6 Описание трехмерных течений через завихренность и векторный потенциал.

Глава 2 Численные модели и их тестирование.

• 2.1 Сетки.

2.2 Метод дробных шагов.

2.3 Алгоритм решения задач конвекции в верхней мантии.

2.4 Модельные уравнения.

2.5 Основные тесты.

2.6 Последовательность сеток.

2.7 Преобразования координат.

Глава 3 Трехмерное моделирование тепловой конвекции в верхней мантии под литосферой континентов. ф 3.1' Основные параметры моделей.

3.2 Тепловая конвекция под литосферной плитой постоянной толщины

3.3 Тепловая конвекция под литосферной плитой с утолщенной полосой в ее центральной части.

3.4 Тепловая конвекция под литосферной плитой с квадратным в плане утолщением в ее центральной части (кратоном).

3.5 Тепловая конвекция под литосферной плитой, содержащей два утолщенных массива.

3.6 Тепловая конвекция под литосферной плитой с утонением в ее центральной части.

Исследование конвекции в недрах Земли является одной из центральных задач геофизики. Работы, в этом направлении, выполненные в последние годы [1,2,22-21,48-52 и др.], значительно расширили наши представления о строении и составе недр. Весьма важная роль в процессе накопления полезной информации по данному вопросу принадлежит численным экспериментам.

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

В работах В.И. Полежаева [19, 33, 36] рассмотрены вопросы двумерного моделирования конвективных процессов, описаны разностные схемы интегрирования уравнений Навье - Стокса в переменных у/,со, подробно рассмотрены различные аппроксимации граничных условий для вихря СО.

В книге Е.Л. Тарунина [43] также рассматриваются решения уравнений вязкой несжимаемой жидкости в переменных у/,60', применяется метод последовательности сеток; излагается метод фиктивных областей для уравнений Навье - Стокса и проанализированы его возможности; приведены примеры решения задач свободной конвекции в замкнутых объемах.

В работе Берковского, Ноготова [7], дается изложение конечно - разностных алгоритмов, предназначенных для исследования задач конвекции. Книга [8] Берковского Б.М., Полевикова В.К. посвящена обсуждению методов решения уравнений гравитационной и термомагнитной конвекции. Приведены примеры задач естественной конвекции; показана зависимость числа Нуссельта от числа Рэлея и используемых разностных схем.

В работе Воеводина А.Ф., Остапенко В.В., Пивоварова Ю.В., Шугрина С.М. [15] для двумерных уравнений вязкой жидкости, записанных в переменных «завихренность - функция тока» рассматривается конвективное течение в замкнутой прямоугольной области. Для завихренности использовалась факторизованная разностная схема; для функции тока - классическая пятиточечная схема. Решение двумерной задачи в нелинейном случае сводится к решению одномерных задач путем применения при аппроксимации конвективных слагаемых подхода Г.И. Марчука [29], основанного на методе слабой аппроксимации. На основе методов дробных шагов и преобразования Фурье рассматриваются прямые и итерационные методы решения системы разностных уравнений. Обоснование сходимости итерационных методов в линейной постановке подробно рассмотрено в статье А.Ф. Воеводина [14]. Граничные условия для вихря вычислялись по формуле Тома [45]. Дальнейшее развитие эти исследования получили в работах А.Ф. Воеводина, Т.В. Юшковой [17] и А.Ф. Воеводина, Т.В. Протопоповой [16].

Пирсоном, в работе [35] рассматривались двумерные нестационарные уравнения (в переменных СО, у/) для описания течения вязкой однородной несжимаемой жидкости, приведено точное аналитическое решение. Интегрирование нестационарного уравнения для вихря осуществлялось неявным конечно - разностным методом Писмана — Рекфорда с центральными разностями для пространственных производных, а уравнение Пуассона для функции тока интегрировалось методом верхней релаксации [12].

Конвективному тепломассообмену в двумерных задачах геодинамики посвящены работы Туркотта, Мак-Кензи, Хауземана, Робертса, Вейса, Кристенсена, Флейто, Йена, Олсона [76, 85, 86,

97, 98, 99, 100, 101, 102, 109, 110, 11 1,112, 113, 122, 127, 128, 129, 131, 132, 136] и многих других.

Torrance & Turcotte [136] использовали двухполевой метод (переменные со,у/) для системы уравнений Навье - Стокса при решении задач конвекции с вязкостью, зависящей от температуры и/или давления. Применялся метод конечных разностей с центральными разностями для производных по пространству и времени. Для аппроксимации конвективных слагаемых в уравнении теп-лопереноса использовались трехточечные несимметричные разностные операторы.

В работах Houseman G.A., McKenzie D.P., Molnar Р., Moore D.R., Weiss N.O., Roberts [97, 109, 127] также применялся двухполевой метод при интегрировании уравнений Навье - Стокса. Moore, Weiss [109], в частности, описывают двумерную конвекцию в жидкости между свободными границами с числами Прандт-ля от 0.01 до бесконечности при постоянной вязкости.

Моделирование двумерной конвекции в мантии Земли при постоянной вязкости выполнено Houseman G.A., McKenzie D.P., Molnar P. [97]. Численное решение нестационарного уравнения теплопереноса осуществлялось по схеме «чехарда» с центральными разностями по времени и пространству и со вторым порядком аппроксимации. Уравнение Пуассона для функции тока интегрировалось методом Фурье и трехдиагональной прогонкой.

Houston, Bremaecker [98] рассматривали двумерные уравнения Навье - Стокса в переменных у/,Т, с уравнением 4-го порядка для функции тока, в процессе моделирования двумерной конвекции с внутренними источниками тепла и с вязкостью, зависящей от глубины. Применялись разнесенные неравномерные сетки. Использовалась схема переменных направлений. Аппроксимация конвективных слагаемых в уравнении теплопереноса осуществлялась направленными разностями. Бигармоническое уравнение для функции тока интегрировалось по схеме стабилизирующей поправки [65,79].

Работа Кристенсена [76] посвящена двумерному моделированию конвекции в мантии Земли на основе уравнений Навье - Сто-кса с вязкостью, зависящей от температуры и давления. Использовалась схема с применением уравнения 4-го порядка для функции тока, которое интегрировалось методом конечных элементов. Функция тока была представлена бикубическими сплайнами, а поле температуры - биквадратными сплайнами. Для чисел Рэлея, превышающих 1000* Racrit, температура вычислялась по формуле верхней релаксации: Тп = (1 — Л)*Тп + X* Т"~[, где параметр релаксации X менялся от 0.5 до 0.9 в зависимости от числа Рэлея.

В работах В.П. Трубицына с соавторами рассмотрены численные модели конвективных процессов, записанных в естественных переменных как для двумерного случая, так и для трехмерного [38,39, 46-51, 138]. Применялся неявный метод расщепления по физическим процессам.

Решая аналогичную задачу Schmeling, Marquart [128, 129] применяли декомпозицию Холесского симметричной матрицы.

Методы конечных элементов применили Hansen, Ebel [93]. Для решения бигармонического уравнения для функции тока использован некомфорный тип конечных элементов и билинейный -для уравнения теплопереноса с коррекцией «по потоку»: Heinrich J.С., Huyakorn P.S., Zienkiewicz О.С., [94].

Fleitout, Yuen [85,86], для численного интегрирования уравнений Навье - Стокса в переменных U,V, р с вязкостью, зависящей от температуры и давления, применили метод, основанный на распределении прямоугольных конечных элементов в Лагранже-вой системе координат.

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

Метод искусственной сжимаемости (Владимирова, Кузнецов, Яненко (1965, [13]); Чорин (1967, [3]) состоит в том, что уравнение несжимаемости заменяется уравнением неразрывности слабо-сжимаемой жидкости. Система уравнений становится системой типа Коши - Ковалевской и для нее может быть применен метод дробных шагов. Однако расчеты ведутся при конечном параметре сжимаемости £ и возникают проблемы, связанные с >0. По -видимому, один из возможных подходов к устранению этого недостатка является метод Ричардсона [30] (экстраполяция по малому параметру с).

Весьма популярным является в настоящее время метод расщепления по физическим процессам (Белоцерковский [6], Пейре, Тейлор [34]). Однако, существенной стороной данного подхода, является отыскание давления (избыточного давления); при этом необходимо решать задачу Неймана для трехмерного уравнения Пуассона.

Для численного интегрирования трехмерных задач гидродинамики несжимаемых жидкостей может быть также применен подход, основанный на переменных О, со (вектор скорости, вектор завихренности) (Роуч, [37]). Однако, как отмечается в [3], при этом могут возникнуть проблемы, связанные с выполнением условия несжимаемости (Азиз [66], Роуч [37], Полежаев [33]).

И, наконец, в классических задачах гидродинамики, применяется трехмерный аналог (^,<г>)-подхода; при этом, в роли скалярного поля завихренности аз выступает трехмерный вектор завихренности &='шх + ]а>у + кй?2 3 = У-0. Аналогом двумерной функции тока у/ является векторный потенциал ip = iy/x+ jу/у +ky/z. Условие несжимаемости выполняется тождественно; давление исключается. Задача сводится (в однородной жидкости) к последовательному интегрированию трех диффузионных уравнений для компонент вихря и трех уравнений Пуассона для компонент векторного потенциала.

Подход с использованием трехмерных векторов завихренности и потенциала успешно применен, в частности, в работе O.A. Бессонова, В.А. Брайловской, В.И. Полежаева, [9], где рассматривалось течение, вызванное градиентами температуры и концентрации в поле силы тяжести в прямоугольном параллелепипеде. Математическое моделирование основывалось на нестационарных уравнениях Навье - Стокса в приближении Обербека - Буссине-ска. Для решения применялись разнесенные сетки, с размещением физических величин в различных местах вычислительной ячейки, что позволило обеспечить консервативность для завихренности на дискретном уровне. Уравнение переноса завихренности интегрировалось с применением неявного метода переменных направлений; уравнения Пуассона для векторного потенциала решалось с помощью быстрого преобразования Фурье по двум пространственным направлениям и прогонкой по третьему.

Тесты двумерной мантийной конвекции анализировались в работах Blankenbach В., Busse F. et al. [70], а также в работах Мошкина Н.П., Рычковой Е.В., Тычкова С.А.,Черных Г.Г. [32,75]. В [32,75] рассматривались два подхода. В первом из них использовались переменные у/, Т. Уравнение переноса тепла интегрировалось с применением метода предиктор-корректор; на этапе предиктора использовались направленные разности. Второй подход основан на явном методе расщепления по физическим процессам.

В работе автора [62] применен подход с использованием разностных схем интегрирования уравнений Навье - Стокса в переменных l//,CD; рассмотрены четыре аппроксимации конвективных слагаемых в уравнении теплопереноса, сделан сравнительный анализ различных аппроксимаций при решении тестовой задачи о конвективном теплопереносе в мантии Земли.

В работе Blankenbach В., Busse F. et al. [70] приведены тесты двумерной мантийной конвекции, в работе Busse F.H., Christensen U. et al. [73] - система тестов трехмерной конвекции Данные, представленные в [73] получены в результате применения следующих методов:

- СВ - (Cleves & Busse) Клевес и Буссэ использовали скалярный потенциал для описания течения; применены спектральный, Галеркина и Ньтона-Рафсона методы решения системы уравнений Навье - Стокса и теплообмена;

- Ch - (Christensen) Кристенсен применил гибридный спектрально - конечно - разностный метод и также как и СВ использовал скалярный потенциал для описания течения;

- Cs - (Cserepes) Церепес получил результаты на основе подхода, аналогичного Ch, но только с постоянной вязкостью;

- Ga - (Gable) Гейбл применил гибридный спектрально -конечно - разностный метод для переменных «скорость -давление»; температура вычислялась явным конечно -разностным методом;

- Но - (Houseman) Хауземан применил гибридный спектрально - конечно - разностный метод для переменных «вихрь - векторный потенциал».

- Og - (Ogava) Огава применил конечно - разностный метод для переменных «скорость - давление»; температура вычислялась неявным конечно - разностным методом.

- PS - (Parmentier & Sotin) Парментер и Сотин применили конечно - разностный метод для полоидального потенциала.

- Тг - (Travis) Трэвис использовал скалярный потенциал для описания течения; применен гибридный спектрально - конечно - разностный метод; температура вычислялась явным конечно - разностным методом.

Результаты трехмерного тестирования можно найти в работах автора [63,64].

По современным представлениям, тепловая конвекция в мантии Земли, эффективно выносящая радиогенное тепло из недр, принимающая непосредственное участие в движении литосферных плит и формировании тектонических режимов территорий, имеет сложную трехмерную структуру, которая определяется строением оболочки, распределением физико-химических характеристик с глубиной и фазовыми переходами вещества. Главной особенностью структуры принято считать слоистость конвекции, т.е. существование замкнутых изолированных ячеек отдельно в верхней и нижней мантии (Добрецов H.JI., Кирдяшкин А.Г., Тычков С.А., Dubuffet F., М. Rabinowicz, М. Monnereau, Hewitt J.С., McKenzie D.P., Weiss N.O. Richter F. [20, 52, 84, 95, 123]). Эта изоляция, однако, может локально нарушаться, что обусловлено периодически возникающей неустойчивостью теплового граничного слоя в районе эндотермического фазового перехода на границе верхней и нижней мантии (Solheim L.P.,Peltier W.R., [131]). Ключевым вопросом при изучении конвекции является выяснение причин и условий, определяющих пространственно-временную эволюцию структуры конвекции, поскольку именно эта характеристика во многом определяет кинематику литосферных плит и геологическую историю развития континентальных областей.

Эволюция конвективных течений в пространстве и времени зависит от целого ряда параметров модели, причем важное место занимают условия на верхней границе конвектирующей области, которые в континентальных областях формирует литосфера. Необходимо отметить, что резкие вариации толщины верхнего кон-дуктивного слоя - литосферы - присущи только континентальным областям, и, как показали последние работы в этом направлении, они существенно влияют на режим верхнемантийной конвекции. Главная задача настоящего исследования состоит в том, чтобы исследовать эффект трехмерности конвекции на ее пространственно-временные характеристики в присутствии неоднородного по толщине верхнего кондуктивного слоя, где утолщенные до 200 км участки литосферы соответствуют древним континентам и/или кратонам, таким как Африка, Восточно-Сибирская платформа, Северная Америка или Индостан.

На первом этапе исследований эффекта литосферы на эволюцию тепловой конвекции использовались исключительно 2и математические модели конвекции. Из выполненных ранее работ в этом направлении необходимо, прежде всего, отметить цикл работ В.П.Трубицына с коллегами, который первый сформулировал данную проблему и показал ее принципиальную важность для динамики континентальной мантии [9,39,46-50]. Основные результаты выполненных группой работ кратко сводятся к тому, что под утолщенной континентальной плитои из-за теплоэкранирующего эффекта формируется устойчивый восходящий поток верхнемантийного конвективного течения. Структура конвекции также меняется, вытягиваясь по горизонтали до аспектного отношения более 3 (аспектное отношение есть отношение горизонтального размера ячейки или расчетной области к вертикальному). Время формирования подобной структуры tj зависит от размеров континентальной плиты. Так, для плиты толщиной в L=70 км ti порядка 0,015 в безразмерных единицах или ~ 225 млн. лет в абсолютном исчислении. С повышением толщины L, время tj уменьшается и составляет всего 0,01 или 150 млн. лет для L=210 км. Время формирования восходящего потока зависит также и от горизонтальных размеров плиты, 1, уменьшаясь на несколько десятков процентов при увеличении горизонтального размера в два раза. Таким образом, под толстой длинной плитой континента перестройка конвекции идет существенно быстрее, чем под тонкой и короткой. Для организации устойчивого восходящего потока под континентом, его размер должен превышать 1>1 или 700 км в размерных единицах. Так для плиты с 1=914 км и L=9 1 км, время формирования восходящего потока под плитой составляет tj~240 млн. лет, правда при этом ширина континентальной плиты увеличивалась со временем, что могло сокращать время перестройки. Было также показано, что при меньших горизонтальных размерах, 1<1, под толстой континентальной плитой формируется устойчивый нисходящий поток конвекции.

Несмотря на то, что 2D модели тепловой мантийной конвекции являются существенным упрощением реальной трехмерной ситуации, практически .до настоящего времени продолжают появляться работы в двухмерном варианте, посвященные изучению взаимоотношения верхнемантийной конвекции под континентами и океанами. В этих работах исследуется, как правило, влияние одной определенной характеристики недр на конвективное перемешивание в мантии. Так, например, расчет общемантийной тепловой конвекции в прямоугольной области с аспектным отношением 8 был выполнен Nakakuki Т., Yuen D.A., Honda S. в [110] с целью изучения эффекта фазовых переходов на 400 км и 660 км и присутствия толстой континентальной литосферной плиты, занимающей половину верхней границы по длине и имеющую вязкость вещества на три порядка выше, чем конвектирующая мантия, на структуру и эволюцию конвективных течений в мантии. Из результатов, полученных авторами, можно выделить формирование восходящего конвективного потока в нижней мантии под континентальной плитой и рассеянного нисходящего под плитой океана. В верхней мантии под континентом также формировался восходящий поток, причем толщина его из-за меньшей вязкости верхней мантии оказывалась существенно меньше, а скорости подъема вещества мантии возрастали (эффект «фокусировки»). После того, как авторы мгновенно убирали континентальную плиту, моделируя тем самым ее интенсивный горизонтальный дрейф, восходящий поток в верхней мантии становился еще более узким, напоминая хвост плюма, и сохранял свое пространственное положение теперь уже в океанической области не менее 1 млрд. лет. Необходимо отметить, что введение в рассмотрение фазовых переходов, и, прежде всего эндотермического перехода на границе 660 км, обусловило режим слоистой конвекции. Doin М.-Р., Fleitout L., Christensen U., в [83], изучали влияние тепловой верхнемантийной конвекции на стабильность или «время жизни» более вязкого, химически более легкого корня континента. Авторы предполагали, задавая соответствующие значения энергии активации и предэкс-поненциального множителя в выражении зависимости вязкости вещества корня континента от Р-Т условий, что вязкость вещества корня в 6-7 раз выше, чем вязкость низлежащей мантии. В результате они получили скорость эрозии континентальной литосферы в

200 км за 800 млн. лет, причем эрозия осуществлялась как снизу, так и сбоку из-за мелкомасштабной конвекции у края континента. Несмотря на малую скорость эрозии, очевидно, что контраст вязкости, а значит и степень деплетирования по летучим, у вещества корней должна быть существенно выше, не менее 2-3 порядков, поскольку «время жизни» корней оценивается в 1-3 млрд. лет по изотопно-геохимическим данным. В работе ЬепагсПс А. [101], автор, используя в целом технологию моделирования конвекции в океанических областях, когда вязкость литосферы и конвектирующей мантии задается единой функцией от Р-Т параметров, ввел в модель континентальной литосферы химические слои, отличающиеся между собой только по плотности. Верхний слой соI ответствовал коре, ниже располагался более плотный слой лито-сферной мантии, который подстилался еще более плотной верхней мантией. Использование автором единой зависимости вязкости от температуры для литосферы и конвектирующей мантии, несмотря на наличие химических слоев, привело к тому, что утолщенная литосфера континентов оказалась над нисходящими конвективными потоками, как и в предыдущих работах БсЬшеПпд Н., МагяиагЧ в. [128,129], где применялась подобная реология. В работе не было получено соответствия между наблюдаемыми и рассчитанными латеральными вариациями теплового потока (другие характеристики - рельеф поверхности, гравитационное поле - в работе не вычислялись), что заставило автора предположить существование значительных латеральных вариаций концентрации радиогенных изотопов в коре континентов. В данной работе, кроме того, была представлена модель с континентами, вязкость литосферы которых на три порядка превосходила вязкость конвектирующей мантии, что близко к рассматриваемой в нашей работе постановке задачи. Вместе с тем следует отметить, что автор использовал нереальные размеры континентов по горизонтали и вертикали. Так толщина континентальной литосферы была выбрана им в 60 км, что почти вдвое меньше толщины зрелой океанической литосферы, а длина континента не превышала 150 км. В результате эта модель не показала реальных вариаций теплового потока, что дало основание автору отвергнуть подобную реологию континентальных плит как нереальную. Вообще необходимо отметить, что автор в своей статье довольно резко выступил против введения в модель континентальной высоковязкой литосферы а priori. Это, считает автор, не соответствует принципу саморегулирующихся моделей тепловой конвекции, когда верхний кондуктивный слой -литосфера - получается автоматически из расчетов. Данная точка зрения или подход к созданию моделей динамики недр континентальных областей является весьма распространенным (Fleitout L., Yuen D.A., Schmeling Н., Marquart G., Tackley P.J. [85, 128, 129, 135]). Следует отметить, что все модели динамики континентальной мантии, базирующиеся на подобном подходе, получили взаимоотношение расчетных геофизических характеристик (рельеф поверхности, тепловой поток и геоид), не соответствующее наблюдениям как было показано в [55] (Тычков С.А., Рычкова Е.В., Василевский А.Н., Червов В.В.). Так, например утоненной континентальной литосфере осадочных бассейнов в этих моделях соответствует положительный рельеф поверхности более 1 км и отрицательные значения высот геоида, хотя, как известно, осадочные бассейны (Западная Сибирь, Прикаспийская низменность и т.п.) представляют собой равнины с максимальным рельефом, не превышающим 0,2 км, а высоты геоида здесь положительны. Если принять подобный подход для континентальных областей, то такая модель, кроме тепловой конвекции, необходимо должна описывать также экстракцию континентальной коры из первичной мантии и формирование более легкой, обедненной главными элементами и летучими компонентами мантийной части литосферы, как это следует из наблюдений. Геологические данные подтверждают, что древние платформы или кратоны существуют уже более 2 млрд. лет. Столь длительное существование литосферы этих областей возможно лишь при условии, что вещество литосферы отлично от конвектирующей мантии по химическому составу. Вместе с тем, хотя механизм формирования континентальной литосферы еще до конца не выяснен [55], сводить природу литосферы континентов к чисто тепловой не вполне корректно. При таком допущении она может быть диссипирована всего за 200-400 млн. лет восходящими потоками тепловой конвекции или плюмами.

Пример лабораторных исследований тепловой конвекции содержится в [90] (виШои Ь., .1аира1Ч С.), где изучался только тепловой эффект континентальной литосферы. Континент был представлен теплоизолирующей пластиной, которая была вставлена в медную плиту. Эта плита поддерживалась при постоянной температуре, моделируя тем самым океаническую литосферу. Подобные граничные условия впервые были использованы в [46] (Трубицын В.П., Белавина Ю.Ф., Рыков В.В.) в численных экспериментах. Результаты лабораторных экспериментов совпали с численными расчетами В.П. Трубицына и подтвердили обнаруженный им эффект формирования устойчивого восходящего конвективного потока под частично теплоизолированной поверхностью, соответствующей континентальной плите.

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

Главный вывод, полученный из 3D моделей, с учетом сжимаемости, состоял в том, что в этом случае нарушалась симметрия между формой восходящих и нисходящих потоков, которая характерна для конвектирующей несжимаемой жидкости. Теперь структура течения представляла собой изолированные узкие восходящие потоки (плюмы), окруженные основной массой медленно погружающейся «холодной» мантии (Bercovici D., G. Shubert, G.A. Glatz-maier, Machetel P., M.Rabinowicz, P. Bernadet, Ratcliff J.Т., P.J. Tackley, G. Schubert, A. Zebib. [69, 88, 103, 119]). На данном этапе удалось исследовать общие закономерности конвектирования жидкости с учетом, помимо сжимаемости, роли фазовых переходов в мантии или ее стратификации по плотности и реологии [80,116] (Cserepes L., М. Rabinowicz, С. Rosemberg-Borot, Ceule-neer G., Monnereau M., Rosemberg С.). Но, вместе с тем, уровень возможностей вычислительной техники того времени, не позволял строить модели с такой разрешающей способностью, чтобы их можно было использовать для изучения эволюции конкретных геологических структур даже в региональном плане, на уровне, например, Альпийско-Гималайского складчатого пояса или Сибирской платформы.

Начальный этап 3D моделирования характеризовался тем, что во всех упомянутых выше моделях вязкость вещества не зависела от температуры и глубины, что было существенным упрощением реальной ситуации [92] (Hager В. Н.). Первые модели трехмерной конвекции с вязкостью вещества, зависящей от температуры и давления, были построены в квазикубических областях (small boxes) и преследовали цель, главным образом, математически корректно промоделировать данные лабораторных экспериментов с температурнозависящей жидкостью (Booker J.R., [71]), чтобы убедиться в надежности разностных схем (Busse F.H., Н. Frick,

Christensen U., H. Hager, Ogawa M., G. Shubert, A. Zebib. [74,77,111]). Главный вывод из моделирования конвекции в квадратных боксах состоял в том, что зависимость вязкости от температуры формирует структуру течения в виде восходящих плюмов, окруженных нисходящими потоками. Структура нисходящих течений в плане была представлена гексагонами и треугольниками.

Следующий этап в 3-D моделировании тепловой конвекции связан с изучением структуры и эволюции течений в больших боксах, горизонтальные размеры которых в 4 и более раз превышают вертикальные. Именно такие расчетные области наиболее удобны или более всего подходят для моделирования конвективных течений в верхней мантии Земли, т.е. в слое до 670 км глубины. Во-первых, увеличение горизонтальных размеров в 4 или 8 раз по сравнению с вертикальными позволяет существенно ослабить влияние боковых границ на структуру и динамику конвективных течений. Во-вторых, при подобном существенном увеличении размеров расчетной области, вычислительные возможности современных персональных компьютеров все же позволяют проводить вычисления одной задачи за относительно короткий интервал времени, 12-24 часов. В-третьих, при данных горизонтальных размерах можно еще пренебрегать эффектом сферичности: так как высота расчетной области верхней мантии составляет всего 10 % от радиуса Земли, то площади верхней и нижней горизонтальных поверхностей разнятся незначительно - не более, чем на 12%. Для примера можно оценить, как соотносятся площади горизонтальных плоскостей при моделировании общемантийной конвекции в слое 3000 км толщины. В этом случае нижняя поверхность из-за сферичности в 4 раза меньше, чем верхняя, поэтому изучать динамику всей мантии в прямоугольных боксах не вполне корректно.

Сейчас существует ряд исследований, посвященных изучению конвективных течений в больших боксах: Dubuffet F., М. Rabinowicz, М. Monnereau. Ratcliff J.T., P.J., Tackley, G. Schubert, A. Zebib - [84,119], среди которых выделяется своей полнотой и последовательностью работа Р. Tackley [134]. Автор начал со случая постоянной вязкости в боксе 8x8x1 при Ra=105. Структура течения характеризуется приблизительной симметрией между восходящими и нисходящими потоками, планформа которых варьирует между квадратами и гексагональной структурами. Когда была введена зависимость вязкости от температуры, нисходящие потоки объединились в один квазицилиндрический, и структура течения трансформировалась в единую ячейку с квадратной планфор-мой, которая заняла, практически, всю расчетную область. С другой стороны, структура тепловой конвекции, вязкость которой зависит только от глубины, представляла собой двухмерные валы, медленно эволюционирующие во времени. И, наконец, был рассмотрен случай одновременной зависимости вязкости и от температуры, и от давления. Структура течения теперь представляла собой суперпозицию случаев зависимости вязкости от температуры и от давления в виде ансамбля ячеек, горизонтальные размеры которых были в два раза больше чем вертикальные и восходящие потоки были окружены нисходящими течениями вещества. Далее автор выполнил тест на влияние аспектного отношения на структуру конвекции. Численные эксперименты на различных размерах расчетной области показали, что структура течений в области 4x4x1 с рефлективными граничными условиями («скользкие» границы) оказалась идентичной структуре в области 8x8x1 с периодическими граничными условиями. Этот тест позволяет сделать вывод о том, что вертикальные границы несущественно влияют на структуру конвекции при ширине области более чем в 4 раза превосходящей ее вертикальный размер. Увеличение же числа Рэлея с 105 до 106 также принципиально не изменило структуру течения, хотя вертикальные потоки становились чуть уже и появлялись модуляции локальных структур конвекции и нестабильности в горизонтальных тепловых слоях. Введение в рассмотрение зависимости вязкости от температуры, при граничном условии на верхней грани в виде жесткой границы, не изменило в значительной мере ни значения числа Нуссельта, ни среднеквадратичную скорость движения вещества в боксе. Таким образом, данная работа была базовой при выборе размеров расчетной области, граничных условий и параметров конвекции в настоящей работе. Более подробно параметры нашей модели будут представлены в главе «Математическая формулировка задачи».

Настоящая работа является продолжением исследований динамики земных недр, выполняющихся совместно Институтом геологии СО РАН и Институтом вычислительных технологий СО РАН в течение последних 8 лет. За это время были получены заметные результаты в изучении особенностей тепловой конвекции в верхней мантии внутриконтинентальных областей на двухмерных численных моделях [24,25,32,40,62]. В частности было показано, что минимальная мощность утолщенной литосферы, начиная с которой латеральные неоднородности литосферы начинают воздействовать на структуру конвективного течения, составляет 11=170 км, что соответствует амплитуде увеличения мощности А 11=50 км. Эффект воздействия утолщенного участка литосферы с 11 > 170 км на структуру течения состоит в формировании устойчивого восходящего потока верхнемантийной конвекции под этим участком. Время выхода системы на устойчивый режим оценивается величиной 2109 лет, что сопоставимо с возрастом древних платформ, хотя формирование восходящего потока под толстой плитой происходит в течение всего 0,2*109 лет. Устойчивый режим структуры конвекции сохраняется длительное время. В рамках экспериментов эта длительность составила 7' 109 лет. Структура конвекции в устойчивом режиме остается практически неизменной при увеличении толщины плиты от 170 до 280 км за исключением формирования локальных конвективных ячеек у краев утолщенного участка при h > 190 км. Размеры этих локальных ячеек определяются, в основном, амплитудой скачка мощности Ah, а скорость перемешивания вещества в них не превосходит среднюю скорость движения в глобальных ячейках, т.е. ~ 1-2 см/год. Участки с аномально тонкой литосферой также влияют на структуру верхнемантийной тепловой конвекции, когда толщина литосферы здесь не более 60 км. При этом образуется нисходящий поток конвекции у одного из краев ловушки из-за интенсивного остывания вещества мантии в ней. Если ловушка расположена у края литосферы докембрийских платформ с литосферой не менее 190 км толщины, то в ней формируются локальные ячейки, что ведет к появлению относительно холодной области ловушки, прилежащей к платформе. Архейские ядра с толщиной литосферы до 350 км отмечены существованием под ними восходящих потоков верхнемантийной конвекции, что не согласуется с некоторыми моделями сейсмотомографии.

Выявленные закономерности имели существенный недостаток в виду того, что динамика мантии моделировалась плоскими структурами течений. Необходимо отметить, что в геологии существуют ситуации, когда двухмерное моделирование оправдано и вполне адекватно описывает реальные ситуации. Это касается, прежде всего, тех геологических структур, размер которых по горизонтали в одном направлении намного превышает размеры в других и также превышает толщину верхней мантии. В качестве примеров можно привести срединно-океанические хребты, протягивающиеся на тысячи километров по дну океанов, или зоны суб-дукции океанической литосферы, имеющие подобные геометрические характеристики. Во внутриконтинентальных же областях, динамику мантии которых мы исследуем, ситуация сложнее. Характерный размер геологических структур, отличающихся по своему глубинному строению и вещественному составу коры, колеблется в пределах от 300 до 1500 км, что сопоставимо с толщиной верхней мантии. Поэтому для корректного сопоставления результатов математического моделирования с наблюдаемыми геологическими и геофизическими характеристиками структур, что является главным критерием адекватности моделей, необходимо учитывать трехмерность объектов исследования. Таким образом, настоящая работа является логическим продолжением исследований мантийной динамики континентальных областей, выполняемых в Институтах Геологии и Вычислительных Технологий СО РАН.

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

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

3D моделирования динамики недр с учетом континентальных плит, следует отметить работу В.П.Трубицына с соавторами [138]. Континенты здесь моделировались жесткими плитами с бесконечной вязкостью, окруженные океанической литосферой, вязкость которой зависела от температуры. Главная цель этой работы состояла в выяснении возможности горизонтального перемещения континентов тепловыми конвективными течениями в мантии Земли. Другая работа (Honda S., М. Yoshida, S. Ootorii, Y. Iwase, [96]) посвящена оценке времени формирования горячего восходящего потока в мантии под жесткой кондуктивной плитой, задававшейся a priori. Данная плита моделировала Пангею - древний суперматерик, включающий в себя большинство современных материков, который распался, как предполагается по геологическим данным около 160 млн. лет. Считается, что один из возможных мантийных механизмов распада Пангеи - расходящиеся конвективные течения возникшего под суперматериком плюма. Авторы проводили эксперименты с 2- и 3-D моделями в прямоугольных областях и со сферическим слоем. Оказалось, что для модели со сферической геометрией и для модели в большом прямоугольном боксе время формирования плюма оценивается величиной в 1-2 млрд. лет, в то время как геологические данные указывают на то, что в течение геологической истории суперматерики образовывались с периодичностью в 0,2-0,4 млрд. лет. Расхождение между экспериментами и фактическими данными можно объяснить тем фактом, что в качестве начального состояния здесь была выбрана модель кон-дуктивного распределения температуры и конвекция возникала и эволюционировала уже в присутствии материков, хотя известно, что средний возраст континентальных образований всего около 2 млрд. лет, а Земли - вдвое больше. То есть конвекция в недрах Земли уже существовала 2 млрд. лет и только потом на ее поверхности образовались заметные континентальные массы, сопоставимые по объему с современными. С учетом этих результатов в нашей модели в качестве начального состояния предполагалось, что конвекция до введения в модель кондуктивной литосферы переменной мощности уже существовала длительное время, и распределение температуры в верхней мантии отвечало квазистационарному режиму конвекции. Идеология наших исследований состоит в последовательном расширении как физической области или размеров моделей, так и во введении в рассмотрение все большего числа эффектов или условий, влияющих на структуру и динамику земных недр в континентальных областях. Выбор континентов в качестве объекта исследований не случаен. Именно в коре и в литосфере в целом континентальных образований содержится информация об эволюции недр в течение практически всей геологической истории Земли, насчитывающей почти 4 млрд. лет. Для примера можно сказать, что средний возраст океанических котловин не превышает 50-70 млн. лет, а самая древняя океаническая литосфера, в районе северо-западной части Тихого океана, имеет возраст не более 200 млн. лет.

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

Анализ цитируемой литературы позволяет сделать вывод о недостаточной полноте численных моделей трехмерной тепловой конвекции в верхней мантии Земли.

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

В двумерных задачах геодинамики этот подход получил широкое распространение [76, 85, 86, 97, 98, 99, 100, 101, 102, 109, 110, 1 11,112, 113, 122, 127, 128, 129, 131, 132, 136].

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

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

Научная новизна изложенных в диссертации результатов заключается в следующем:

• построение и детальное тестирование трехмерной численной модели конвекции в верхней мантии Земли;

• создание геодинамических моделей тепловой конвекции в верхней мантии континентов.

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

- трехмерная численная модель конвекции в верхней мантии Земли;

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

Апробация работы Основные результаты диссертации докладывались на международных конференциях: совещании по геодинамической эволюции палеоозианского океана (Beijing, 1991), мемориальной конференции памяти Зоненшайна по Тектонике Плит. (Москва, 1993) и международной конференции «Потоки и структуры в жидкостях» (С. Петербург, 2003). Полное содержание диссертации докладывалось на семинарах Института Геологии СО РАН (рук. д.г.-м.н. А.Ю. Казанский), Института Гидродинамики им. М.А. Лаврентьева СО РАН (рук. проф. А.Ф. Воеводин), Института Вычислительных Технологий СО РАН (рук. проф. В.М. Ковеня), Института Вычислительной математики и математической геофизики СО РАН (рук. академик A.C. Алексеев). Публикации.

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

Разработанная численная модель может быть использована для моделирования широкого класса задач геодинамики недр планеты. Представленные в диссертации результаты получены при проведении исследований в рамках Приоритетного направления СО РАН «Геодинамическая и геохимическая эволюция литосферы и мантии Земли: тектоника, магматизм, флюидный режим и металлогения», по Программе СО РАН №26.2. «Геодинамическая эволюция литосферы Азии: тектоника, магматизм, метаморфизм, геохимия и металлогения основных ее этапов» Структура и объем дисертации.

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

Основные результаты диссертации сводятся к следующему.

1. Разработана численная модель трехмерной конвекции в мантии Земли, основанная на переменных «векторный потенциал - завихренность» и методе дробных шагов. Осуществлено детальное тестирование построенной численной модели.

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

3. Принципиальный результат трехмерного моделирования тепловой конвекции в верхней мантии под кратоном состоит в выявлении мелкомасштабной моды конвекции непосредственно под литосферой на «астеносферном» уровне глубин 200-350 км. Данная мода развивается по периферии утолщенного участка литосферы (кратона) и может объяснить особенности режима траппового магматизма древних кратонов по их периферии.

Заключение

1. Алексеев A.C., Лаврентьев М.М., Мухометов Р.Г., Нерсесов И.Л., Романов В.Г. Численный метод определения структуры верхней мантии Земли // Математические проблемы геофизики; 1971; Вып. 2. Новосибирск: ВЦ СО АН СССР. с. 143-165.

2. Алексеев A.C., Рябой В.З. Модель строения верхней мантии по объемным сейсмическим волнам // Строение земной коры и верхней мантии по данным сейсмических исследований. Киев: Наукова думка; 1977; с.67-82.

3. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. М., «Мир», 1990.

4. Арфкен Г. Математические методы в физике. М., Атомиздат, 1970.

5. Белоцерковский О.М. Обтекание кругового цилиндра с отошедшей удар-ной.волной // ДАН СССР; 1957; т. 113, №3, с.509-512.

6. Белоцерковский О.М. Численное моделирование в механике сплошных сред. М., «Наука», 1984.

7. Берковский Б.М., Ноготов Е.Ф. Разностные методы исследования задач теплообмена. Минск, Наука и техника, 1976.

8. Берковский Б.М., Полевиков В.К. Вычислительный эксперимент в конвекции. Минск, «Университетское», 1988.

9. Бессонов O.A., Брайловская В.А., Полежаев В.И. Пространственные эффекты конвекции в расплавах: концентрационные неоднородности, возникновение несимметрии и колебания // Механика жидкости и газа; 1997; № 3. с.74-82.

10. Бобров A.M., Трубицын В.П. Времена перестроек структуры мантийных течений под континентами // Физика Земли; 1995; №7, с.5-13.

11. Булгаков В.К., Соловьев C.B. Модели тепловой конвекции в мантии и ядре Земли. М.: Наука, 2001, 239 с.

12. Вазов В., Форсайт Дж. Разностные методы решения дифференциальных уравнений в частных производных. М.: ИЛ, 1963.

13. Владимирова H.H., Кузнецов Б.Г., Яненко H.H. Численный расчет симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости. В кн.: Некоторые вопросы вычислительной и прикладной математики. Новосибирск, «Наука», 1966, с. 186-192.

14. Воеводин А.Ф. Устойчивость и реализация неявных схем для уравнения Стокса // ЖВМ и МФ; 1993; т.33, №1. с.119-130.

15. Воеводин А.Ф., Остапенко В.В., Пивоваров Ю.В., Шугрин С.М. Проблемы вычислительной математики. Новосибирск: Изд-во СО РАН, 1995.

16. Воеводин А.Ф., Протопопова Т.В. Метод расчета вязких течений в замкнутых областях // Сиб. Журнал Индустриальной Математики; 2001; т.4, №1, с.29-37.

17. Воеводин А.Ф., Юшкова Т.В. Численный метод решения начально краевых задач для уравнений Навье — Стокса в замкнутых областях на основе метода расщепления // Сиб. Журнал вычисл. Математики; 1999; т.2, №4, с.321-332.

18. Годунов С.К., Рябенький B.C. Разностные схемы. М., «Наука», 1973.

19. Грязнов В.Л., Полежаев В.И. Исследование некоторых разностных схем и аппроксимаций граничных условий для численного решения уравнений тепловой гравитационной конвекции. М., 1974.

20. Добрецов H.JI., Кирдяшкин А.Г. Глубинная геодинамика. НИЦ ОИГГМ СО РАН, Новосибирск, 1994.

21. Добрецов H.JI., Кирдяшкин А.Г., Кирдяшкин A.A. Глубинная геодинамика, 2-е издание, Новосибирск, Изд-во СО РАН, 409 е., 2001.

22. Добрецов H.JI. Пермо-трасовый магматизм и осадконакопление в Евразии как отражение суперплюма, ДАН, 1997,т.354, с.220-223

23. Коновалов А.Н. Численное решение задачи теории упругости. Лекции для студентов НГУ Новосибирск; 1968; 128 стр.

24. Коробицина Ж.Л., С.А.Тычков. Численное моделирование процессов тепло- и массопереноса с учетом фазового перехода в геодинамике // Журнал вычислительной математики и математической физики; 1997; т.37, №6, с. 733-741.

25. Коробицына Ж.Л., Тычков С.А. Численное моделирование взаимодействия тепловых конвективных течений в мантии Земли и континентальной литосферы с учетом фазового перехода // Вычислительные технологии; 1995; т.4, №13, с.212-223.

26. Ладыженская O.A. Математические вопросы динамики вязкой несжимаемой жидкости. М., «Наука», 288 стр., 1970.

27. Ландау Л.Д., Лифшиц Е.М. Гидродинамика, М.: Наука, 1986.

28. Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1978.

29. Марчук Г.И. Методы расщепления. М.: Наука, 1988.

30. Мярчук Г.И., Шайдуров В.В. Повышение точности решения разностных схем. М., «Наука», 320 стр., 1979.

31. Монин A.C., Яглом A.M. Статистическая гидродинамика. Механика турбулентности. Часть I. М.: Наука, 1965.

32. Мошкин Н.П., Рычкова Е.В., Тычков С.А., Черных Г.Г. Тестирование некоторых численных моделей конвективных течений применительно к задачам геодинамики // Вычислительные технологии; 1995; т.4, №13, с.224-232.

33. Пасконов В.М., Полежаев В.И., Чудов Л.А. Численное моделирование процессов тепло- и массообмена. М., «Наука», 1984.

34. Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. Ленинград, Гидрометеоиздат, (1986).

35. Пирсон К. Е. Численный метод для задач вязкого потока. // Механика, №6, с. 65-77, (1965).

36. Полежаев В.И., Грязнов В.Л. Метод расчета граничных условий для уравнений Навье -Стокса в переменных вихрь, функция тока. //ДАН СССР, т. 219, № 2., (1974).

37. Роуч X. Вычислительная гидродинамика. М., «Мир», (1980).

38. Рыков В.В., Трубицын В.П. Трехмерная модель мантийной конвекции с движущимися континентами // Вычислительная сейсмология; 1994; Вып. 27, с. 21-41.

39. Рыков В.В., Трубицын В.П. Численное моделирование трехмерной мантийной конвекции и тектоника литосферных плит // Вычислительная сейсмология; 1994; Вып. 26, с. 94-102.

40. Рычкова Е.В., Тычков С.А. Численная модель тепловой конвекции в верхней мантии Земли под литосферой континентов // Вычислительные технологии; 1997; т.2, №5, с.66-81.

41. Самарский A.A. Введение в численные методы. М., «Наука», 1982.

42. Самарский A.A. Теория разностных схем. . М., «Наука», 1977.

43. Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. Издательство Иркутского университета, 1990.

44. Тарунин E.JI. Метод последовательности сеток для задач свободной конвекции // Журн. Вычисл. матаматики и мат. физики; 1975; т. 15, №2, стр. 436-445.

45. Том А., Эйплт К. Числовые расчеты полей в технике и физике. М., 1964.

46. Трубицын В.П., Белавина Ю.Ф., Рыков В.В. Тепловое и механическое взаимодействие мантии с континентальной литосферой // Физика Земли; 1993; №11, с. 3-15.

47. Трубицын В.П., Белавина Ю.Ф.,Рыков В.В. Тепловая конвекция в мантии с переменной вязкостью и континентальной плитой конечных размеров // Физика Земли; 1994; №7, с. 5-17.

48. Трубицын В.П., Бобров A.M., Кубышкин В.В. Влияние континентальной литосферы на структуру мантийной тепловой конвекции // Физика земли; 1993; №5, с. 3-11.

49. Трубицын В.П., Рыков В.В. Механизм формирования наклонных зон суб-дукции// Физика Земли; 1997; №6, с. 3-14.

50. Труб ицын В.П., Рыков В.В., Трубицын A.U. Конвекция и распределение вязкости в мантии // Физика Земли; 1997; №3, с.3-10.

51. Трубицын В.П., Фрадков А.С. Конвекция под континентами и океанами // Физика Земли; 1985; №7, с. 3-13.

52. Тычков С.А. Конвекция в мантии и динамика платформенных областей. . Новосибирск, «Наука» СО, 1984.

53. Тычков С.А., Рычкова Е.В., Василевский А.Н. Взаимодействие плюма и тепловой конвекции в верхней мантии под континентом // Геология и геофизика; 1998; т.39, №4, с. 413-425.

54. Тычков С.А., Рычкова Е.В., Василевский А.Н. Тепловая конвекция в верхней мантии Земли под литосферой переменной мощности // Физика Земли; 1999; №9, с. 38-51.

55. Тычков С.А., Рычкова Е.В., Василевский А.Н., Червов В.В. Тепловая конвекция в верхней мантии континентов и ее эффект в геофизических полях // Геология и геофизика; 1999; т.40, №9, с. 1275-1290.

56. Тычков С.А., Василевский А.Н., Рычкова Е.В. Эволюция плюма под континентальной литосферой с резкими вариациями толщины //Геология и геофизика. 1999. Т.40, №8. с. 1 182-1 196.^

57. Федоренко P.П. О скорости сходимости одного итерационного процесса // Журн. Вычисл. Математики и мат. Физики; 1964; т.4, № 3, стр. 559-564.

58. Федорюк М.В. Характеристики течений несжимаемой жидкости в гравитационном поле//Мат. Сборник 1988; 137(179), №4(12).- С. L/?¿

59. Флетчер К. Вычислительные методы в динамике жидкостей; 1991; т. 1,2. М., «Мир».

60. Червов В.В. Сравнение некоторых аппроксимаций конвективных членов в уравнении теплопереноса для задачи конвекции в мантии Земли. // Журн. Вычислительные технологии; 1995; т.4, №13, с. 295-305.

61. Червов В.В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением завихренности и векторного потенциала // Вычислительные технологии; 2002; т.7, №1, с. 114-125.

62. Червов В.В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением последовательности сеток. // Вычислительные технологии; 2002; т.7, № 3, с. 85-92.

63. Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. «Наука», СО, Новосибирск, 1967.

64. Aziz К., Heliums J.D. Numerical solution of the three- dimensional equations of motion for laminar natural convection//Phys. Fluids; 1967; v. 10, N. 2, p. 314-324.

65. Balachadar S., Yuen D.A., Reuteler D.M., Lauer G.S. Viscous dissipation in three- dimensional convection with temperature -dependent viscosity // Science; 1995; v. 267.

66. Basu A.R., Poreda R.J., Renne P.R., Teichman F., Vasiliev Y.R., Sobolev N.V., Turin B.D. High 3He plume origin and temporal-spatial evolution of the Siberian flood basalts // Science; 1995; v.269, p.822-825.

67. Bercovici D., Shubert G., Glatzmaier G.A. Three-dimensional spherical models of convection in the Earth's mantle // Science; 1989; v.244, p.893-1016.

68. Blankenbach В., Busse F. et al. A benchmark comparison for mantle convection codes // Geophys. J.; 1989; Int 98, p. 23-38.

69. Booker J.R. Thermal convection with strongly temperature dependent viscosity // J. Fluid Mech.; 1976; v.76, p.741-754.

70. Busse F.H. High Prandtle number convection. // Physics of the Earth and Planetary Interiors; 1979; v. 19, p. 149.

71. Busse F.H., Frick H. Sguare-pattern convection in a layer heated from below // J. Fluid Mech.; 1985; v.150, p.451-465.

72. Christensen U. Convection with pressure- and temperature-dependent non-Newtonian rheology // Geophys. J. R. astr. Soc.; 1984; v.77, p. 343-384.

73. Christensen U., H. Hager. 3-D convection with variable viscosity // Geophys. J.; 1991; Int., 104, p.213-226.

74. Condie K.C. Episodic continental growth models: afterthoughts and extensions // Tectonophys.; 2000; v.322, p.153-162.

75. Conte S.D., Dames R.J. An alternating directions method for solving the bi-harmonic equation // Math. Tables Aids Comput.; 1958; v.12. p.198-205.

76. Cserepes L., M. Rabinowicz, C., Rosemberg-Borot Three dimensional infinite Prandtle number convection in one and two layers with implications for the Earth gravity field // J. Geophys. Res.; 1988; v.43, p. 12009-12025.

77. Davies G.F., Richards M.A. Mantle convection // J. Geol.; 1992; v.3, p. 151 -206.

78. De Bremaecker J. CI. Convection in the Earth's mantle // Tectonophysics; 1977; v.41, p.195.

79. Doin M.-P., Fleitout L., Christensen U. Mantle convection and stability of depleted and undepleted continental lithosphere // J. Geophys. Res.; 1997; v.102, p.2771-2787.

80. Dubuffet F., Rabinowicz M., Monnereau M. Multiple scales in mantle convection // Earth Planet. Sci. Lett.; 2000; v.178, p.351-366.

81. Fleitout L., Yuen D.A. Steady state, secondary convection beneath litho-spheric plates with temperature- and pressure- dependent viscosity // J. Geo-phys. Res.; 1984; v.89, N Bll, p.9227-9244.

82. Fleitout L., Yuen D.A. Secondary convection and the growth of the oceanic lithosphere // Physics of the Earth and Planetary Interiors; 1984; v.36, p. 181 -212.

83. Gable W., O'Connell J. Convection in three dimensions with surface plates: generation of toroidal flow // J. of geophysical research; 1991; v.96, No. B5, p.8391-8405.

84. Glatzmaier G.A. Numerical simulations of mantle convection: Time-dependent, three-dimensional, compressible, spherical shell // Geophys. Astro-phys. Fluid Dyn.; 1988; v.43, p. 223-264.

85. Glatzmaier G.A., Schubert G. Three dimensional spherical model of layered and whole mantle convection // J. of geophysical research; 1993; v.98, No. B12, p.21969-21976.

86. Guillou L., Jaupart C. On the effect of continents on mantle convection // J. of Geophysical Research; 1995; v. 100, p.24217-24238.

87. Gurnis M. Large-scale mantle convection and the aggregation and dispersal of supercontinents // Nature; 1988; No.332, p.695-699.

88. Hager B. H. Subducted slabs and'the Geoid; constraints on mantle reology and flow // J. Geophys. Res.; 1984; v.89, p.6003-6015.

89. Hansen U., Ebel A. Time-dependent thermal convection possible explanation for a multiscale flow in the Earth's mantle //Geophsical Journal; 1988; v.94, No.2, p.181-191.

90. Heinrich J.C., Huyakorn P.S., Zienkiewicz O.C. An upwind finite element scheme for two-dimensional convective equation // Int. J. Numer. Math. Eng.; 1977; v.l 1, p.131-143.

91. Hewitt J.C., McKenzie D.P., Weiss N.O. Large aspect ratio cells in two-dimensional thermal convection//Earth Planet Sci. Lett.; 1980; v.51, p.370-380.

92. Honda S., Yoshida M., Ootorii S., Iwase Y. The timescales of plume generation caused by continental aggregation // Earth Planet. Sci. Lett.; 2000; v. 176, p.31-43.

93. Houseman G.A., McKenzie D.P., Molnar P. Convective instability of a thickened boundary layer and its relevance for the thermal evolution of continental convergent belts // J. of Geophysical Research; 1981; v.86, No.B7, p.6115-6132.

94. Houston M.H., Jr., De Dremaecker J.CI. ADI solution of free convection in a variable viscosity fluid // J. of Computational Physics; 1974; v.16. p.221-239.

95. Jarvis G.T., Peltier. W.R. Mantle convection as a boundary layer phenomenon // Geophys. J. Roy. Astr. Soc; 1982; v.68, p.389-427.

96. Karato S., Wu P. Rheology of the upper mantle: A synthesis. // Science; 1993; v.260, p.771-778.

97. Lenardic A. On the head flow variation from Archean cratons to Proterozoic mobile belts // J. Geophys. Res.; 1997; v.102, p.709-721.

98. Lowman J.P., Jarvis G.T. Mantle convection flow reversals due to continental collisions // J. Geophys. Res.; 1993; v.20, p.2087-2090.

99. Machetel P., M.Rabinowicz, P. Bernadet. Three-dimensional convection in spherical shells, Geophys. Astrophys. Fluid Dyn.; 1986; v.37, p.57-84.

100. Machetel P., Thoravan C., Brunet D. Spectral and geophysical of 3-D spherical mantle convection with an endothermic phase change at the 670 km discontinuity // Phisics of the Earth and planetary interiors; 1995; v.88, p.43-51.

101. Malevsky A.V. Pattern of convective turbulence: an effect of Prandtl number// Phisics of the Earth and planetary interiors; 1995; v.88, p.31-41.

102. McKenzie D., Weiss N. Speculation on the thermal and tectonic history of the earth // Geophys. J. R. Astr. Soc.; 1975; v.48, p. 131-174.

103. McKenzie D.P., Roberts J.M., Weiss N.O. Convection in Earth's mantle: towards a numerical simulation // J. Fluid Mech.; 1974; v.62, part 3, p.465-538.

104. Megnin C., Bunge H-P., Romanowicz B., Richards M.A. Imaging 3-D spherical convection models: What can seismic tomography tell us about mantle dynamics? // Geophysical research letters; 1997; v.24, No. 11, p. 1299-1302.

105. Moore D.R., Weiss N.O. Two-dimensional Rayleigh-Benard convection // J. Fluid Mech.; 1973; v.58, part 2, p.289-3 12.

106. Nakakuki T., Yuen D.A., Honda S. The interaction of plumes with the transition zone under continents and oceans // Earth Planet. Sci. Lett.; 1997; v. 146, p.379-391.

107. Ogawa M., G. Shubert, A. Zebib. Numerical simulations of three-dimensional thermal convection in a fluid with strongly temperature dependent viscosity // J. Fluid Mech.; 1991; v.233, p.299-328.

108. Olson P, Silver P.G., Carlson R.W. The large-scale structure of convection in the Earth's mantle // Nature; 1990; N.344, p.209-215.

109. Olson P. Mantle convection with spherical effects // J. Geophys. Res.; 1981; v.86, p.4881.

110. Olson P., Schubert G., Anderson C., Goldman P. Plume formation and lithosphere erosion: A comparison of laboratory and numerical experiments // J. Geophys. Res.; 1988; p. 15065-15084.

111. Parmentier E.M., Sotin C., Mravis B.J. Turbulent 3D thermal convection in an infinite Prandtl number, volumetrically heated fluid: implications for mantle dynamics // Geophis. J. Int.; 1994; v. 116, p.241-251.

112. Rabinowicz M., Ceuleneer G., Monnereau M., Rosemberg C. Three dimensional models of mantle flow across a low- viscosity zone: implications for hotspot dynamics // Earth and planetary science letter; 1990; v.99, p. 170-184.

113. Rabinowicz M., Rouzo S., Sempere J.C. Three dimensional mantle flow beneath mid - ocean ridges // J. of geophysical research; 1993; v.98, No.B5, p.7851-7869.

114. Rabinowicz R. A., Lago B. Thermal transfer between the continental asteno-sphere and oceanic subducting litosphere: it's effect on subcontinental convection // J. Geophys. Res.; 1980; v.85, p.1839.

115. Ratcliff J.T., P.J. Tackley, G. Schubert, A. Zebib. Transitions in thermal convection with strongly variable viscosity // Phys. Earth Planet. Inter.; 1997; v. 102, p.201-212.

116. Richter F. Convection and the large-scale circulation of the mantle // J. Geophys. Res.; 1973; v.78, p.8735.

117. Richter F. Experiments on the stability og convection rolls in fluid, whose viscosity depeds upon temperature // J. Fuid Mech.; 1978; v.89, p.533.

118. Richter F. Finite amplitude convection throuth a phase boundary // J. Geoph. Roy. Astr. Soc.; 1973; v.35, p.265.

119. Richter F., Mckenzie D.P. On some consequence and possible causes of layered mantle convection // J. Geophys. Res.; 1981; v.86, p.6133-6147.

120. Richter F., Parson B. On the interaction of two scale of convection in the mantle // J. Geophys. Res.; 1985; v.80, p.2529.

121. Ritsema J., Nyblade A.A., Owens T.J. et al. Upper mantle seismic velocity structure beneath Tanzania, east Africa: implication for the stability of craton lithosphere // J. Geophys. Res.; 1998; v.103, p.21201-21213.

122. Ritzwoller M.H., Levshin AL. Eurazian surface wave tomography: group velocities // J. Geophys. Res.; 1998; v. 103, p.4839-4878.

123. Schmalze J., Hansen U. Mixing properties of mantle convection // Eug. VII, Terra Abst. Suppl. N 1 to Terra Nova. 1993. N 5. p.56

124. Schmeling H., Marquart G. Mantle flow and evolution of the lithosphere // Phys. Earth Planet. Inter.; 1993; v.79, p.241-267.

125. Schmeling H., Marquart G. The influence of second-scale convection on the thickness of continental lithosphere and crust // Tectonophys.; 1991; v.189, p.281-306.

126. Simons F.J., Zeihuis A, van der Hilst R.D. The deep structure of the Australian continent from surface wave tomography // Lithos, 1999; v.48, p. 17-43.

127. Solheim L.P.,Peltier W.R. Mantle phase transitions and layered convection // Can. J. Earth Sci.; 1993; v.30, p.881-892.

128. Spalding D.B. A novel finite difference formulation for differential expressions involving both first and second derivatives // Int. J. Numer. Math. Eng.; 1972; v.4, p.551-559.

129. Tackley P.J. Effects of strongly variable viscosity on three-dimensional compressible convection in planetary-mantles // J. Geophys. Res.; 1996; v. 101, p.331 1-3332.

130. Tackley P.J. Self- consistent generation of tectonic plates in three- dimensional mantle convection // Physics of the Earth Science letters; 1998; v. 157, p.9-22.

131. Torrance K.E., Turcotte D.L. Thermal convection with large viscosity variations // J. fluid Mech.; 1971; v.47, pat 1, p.113-125.

132. Travis B., Olson P., Shubert G. The transition from two-dimensional to three-dimensional planform in infinite Prandtl number thermal convection // J. Fluid Mech.; 1990; v.216, p.71-91.

133. Trubitsyn V.P., Rykov V.V. A 3D numerical model of the Wilson cycle // J. Geodynam.; 1996; v.20, p.63-75.

134. Van der Lee S., Nolet G. Upper mantle S velocity structure of Noth America //J. Geophys. Res., 1997; v.102, p.228 15-22838.

135. White R., McKenzie D. Magmatism at rift zones: the generation of volcanic continental margins and flood basalts // J. Geophys. Res.; 1989; v.94, p.7685-7729.

136. Whitehead J.A.Jr. Convection models: laboratory versus mantle // Tectono-physics; 1976; v.35, p.215.

137. Widmer-Schnidrig R. Free oscillations illuminate the mantle // Nature; 1999; v.398, p.292-293.

138. Zhang K., Busse F.H. Convection in spherical fluid shells with an outer crust of variable thickness // Physics of the Earth and planetary interiors; 1997; v. 104, p.283-294.

139. Zhang S., Yuen D.A. Intense local toroidal motion generated by variable viscosity compressible convection in 3-D spherical -shell // Geophysical research letters; 1996; v.23, No.22, p.3135-3138.

140. Zhong S., Zuber M. Role of temperature- dependent viscosity and surface plates in spherical shell models of mantle convection // J. of geophysical research; 2000; v.105, No.B5, p.11063-1 1082.