Трехмерное моделирование конвективных процессов в мантии земли Червов Виктор Васильевич

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

Введение

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

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

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

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

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

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

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

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

2.1 Сетки 46

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

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

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

2-5 Основные тесты 58

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

2-7 Преобразования координат 31

Глава 3 Трехмерное моделирование тепловой конвекции в верхней мантии под литосферой континентов 87

3.1. Основные параметры моделей 88

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

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

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

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

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

Заключение. 176

Список литературы 177

• Основные уравнения Движения несжимаемой жидкости
• Алгоритм решения задач конвекции в верхней мантии
• Тепловая конвекция под литосферной плитой с утолщенной полосой в ее центральной части
• Тепловая конвекция под литосферной плитой с утонением в ее центральной части

Введение к работе

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

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

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

В книге Е.Л. Тарунина [43] также рассматриваются решения уравнений вязкой несжимаемой жидкости в переменных у/?СО\

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

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

естественной конвекции; показана зависимость числа Нуссельта от числа Рэлея и используемых разностных схем,

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

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

Конвективному тепломассообмену в двумерных задачах геодинамики посвящены работы Туркотта, Мак-Кензи, Хауземана, Робертса, Вейса, Кристенсена, Флейто, Йена, Олсона [76, 85, 86, 97, 98, 99, 100, 101, 102, 109, ПО, 111,112, 113, 122, 127, 128, 129, 131, 132, 136] и многих других.

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

В работах Houseman G.A., McKenzie D.P., Molnar P., 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 Ля , температура вычислялась по формуле верхней релаксации: 74" =(1—Я) 7,Я+Л Г 1-1, где параметр релаксации Я менялся от 0.5 до 0.9 в зависимости от числа Рэлея.

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

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

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

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

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

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

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

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

И, наконец, в классических задачах гидродинамики, применяется трехмерный аналог ( У)-подхода; при этом» в роли скалярного поля завихренности со выступает трехмерный вектор завихренности со= iuj + }й)у -i-ko)z t5=V-U. Аналогом двумерной функции тока у/ является векторный потенциал ір = іу/х+ jy/y +ky/z, Условие несжимаемости выполняется тождественно; давление исключается. Задача сводится (в однородной жидкости) к последовательному интегрированию трех диффузионных уравнений для компонент вихря и трех уравнений Пуассона для компонент векторного потенциала.

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

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

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

В работе Blankenbach В,, Busse F. et а], [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) Парментер и Сотин применили конечно - разностный метод для полоидального потенциала.

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

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

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

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

На первом этапе исследований эффекта литосферы на эволюцию тепловой конвекции использовались исключительно 2D математические модели конвекции- Из выполненных ранее работ в этом направлении необходимо, прежде всего, отметить цикл работ В.П.Трубицына с коллегами, который первый сформулировал данную проблему и показал ее принципиальную важность для динамики континентальной мантии [9Э39,46-50], Основные результаты выполненных группой работ кратко сводятся к тому, что под утолщенной континентальной плитой из-за теплоэкранирующего эффекта формируется устойчивый восходящий поток верхнемантийного конвективного течения. Структура конвекции также ме няется, вытягиваясь по горизонтали до аспектного отношения более 3 (аспектное отношение есть отношение горизонтального размера ячейки или расчетной области к вертикальному). Время формирования подобной структуры tj зависит от размеров континентальной плиты. Так, для плиты толщиной в L=70 км t] порядка 0,015 в безразмерных единицах или 225 млн. лет в абсолютном исчислении- С повышением толщины L, время ti уменьшается и составляет всего 0,01 или 150 млн. лет для L=210 км. Время формирования восходящего потока зависит также и от горизонтальных размеров плиты, I, уменьшаясь на несколько десятков процентов при увеличении горизонтального размера в два раза. Таким образом, под толстой длинной плитой континента перестройка конвекции идет существенно быстрее, чем под тонкой и короткой. Для организации устойчивого восходящего потока под континентом, его размер должен превышать 1 1 или 700 км в размерных единицах. Так для плиты с 1=914 км и L=91 км, время формирования восходящего потока под плитой составляет 11 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 млрд. лет по изотопно-геохимическим данным, В работе Lenardic А. [101], автор, используя в целом технологию моделирования конвекции в океанических областях, когда вязкость литосферы и конвектирующей мантии задается единой функцией от Р-Т параметров, ввел в модель континентальной литосферы химические слои, отличающиеся между собой только по плотности- Верхний слой со-ответствовал коре, ниже располагался более плотный слой лито-сферной мантии, который подстилался еще более плотной верхней мантией. Использование автором единой зависимости вязкости от температуры для литосферы и конвектирующей мантии, несмотря на наличие химических слоев, привело к тому, что утолщенная литосфера континентов оказалась над нисходящими конвективными потоками, как и в предыдущих работах Schmeling Н,, Marquart G. [128,129]., где применялась подобная реология, В работе не было получено соответствия между наблюдаемыми и рассчитанными латеральными вариациями теплового потока (другие характеристики - рельеф поверхности, гравитационное поле - в работе не вычислялись), что заставило автора предположить существование значительных латеральных вариаций концентрации радиогенных изотопов в коре континентов. В данной работе, кроме того, была представлена модель с континентами, вязкость литосферы которых на три порядка превосходила вязкость конвектирующей мантии, что близко к рассматриваемой в нашей работе постановке задачи. Вместе с тем следует отметить, что автор использовал нере альные размеры континентов по горизонтали и вертикали. Так толщина континентальной литосферы была выбрана им в 60 км, что почти вдвое меньше толщины зрелой океанической литосферы, а длина континента не превышала 150 км. В результате эта модель не показала реальных вариаций теплового потока, что дало основание автору отвергнуть подобную реологию континентальных плит как нереальную. Вообще необходимо отметить, что автор в своей статье довольно резко выступил против введения в модель континентальной высоковязкой литосферы a priori. Это, считает автор, не соответствует принципу саморегулирующихся моделей тепловой конвекции, когда верхний кондуктивный слой -литосфера - получается автоматически из расчетов. Данная точка зрения или подход к созданию моделей динамики недр континентальных областей является весьма распространенным (Fleitout L., Yuen D.A,, Schmeling H., Marquart G., Tackley P.J. [85, 128, 129, 135]). Следует отметить, что все модели динамики континентальной мантии, базирующиеся на подобном подходе, получили взаимоотношение расчетных геофизических характеристик (рельеф поверхности, тепловой поток и геоид), не соответствующее наблюдениям как было показано в [55] (Тычков С.А., Рычкова Е.В., Василевский А.Н., Червов В.В.). Так, например утоненной континентальной литосфере осадочных бассейнов в этих моделях соответствует положительный рельеф поверхности более 1 км и отрицательные значения высот геоида, хотя, как известно, осадочные бассейны (Западная Сибирь, Прикаспийская низменность и т.п.) представляют собой равнины с максимальным рельефом, не превышающим 0,2 км, а высоты геоида здесь положительны. Если принять подобный подход для континентальных областей, то такая модель, кроме тепловой конвекции, необходимо должна описывать также экстракцию континентальной коры из первичной мантии и формирование более легкой, обедненной главными элементами и летучими компонентами мантийной части литосферы, как это следует из наблюдений. Геологические данные подтверждают, что древние платформы или кратоны существуют уже более 2 млрд. лет. Столь длительное существование литосферы этих областей возможно лишь при условии, что вещество литосферы отлично от конвектирующей мантии по химическому составу. Вместе с тем, хотя механизм формирования континентальной литосферы еще до конца не выяснен [55], сводить природу литосферы континентов к чисто тепловой не вполне корректно. При таком допущении она может быть диссипирована всего за 200-400 млн. лет восходящими потоками тепловой конвекции или плюмами.

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

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

Главный вывод, полученный из 3D моделей, с учетом сжимаемости, состоял в том, что в этом случае нарушалась симметрия между формой восходящих и нисходящих потоков, которая характерна для конвектирующей несжимаемой жидкости. Теперь структура течения представляла собой изолированные узкие восходящие потоки (плюмы), окруженные основной массой медленно погружающейся «холодной» мантии (Bercovici D., G. Shubert, G.A- Glatz-maier, Machetel P., M.Rabinowicz, P. Bernadet, Ratcliff J,Т., PJ. Tackfey, 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, Chnstensen IL, 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., М. Rabinowlcz, М. Monnereau. Ratcliff J .Т., P.J., Tackley, G. Schubert, A- Zebib - [84,119], среди которых выделяется своей полнотой и последовательностью работа P. Tackley [134]. Автор начал со случая постоянной вязкости в боксе 8x8x1 при Ra=l05. Структура течения характеризуется приблизительной симметрией между восходящими и нисходящими потоками, планформа которых варьирует между квадратами и гексагональной структурами. Когда была введена зависимость вязкости от температуры, нисходящие потоки объединились в один квазицилиндрический, и структура течения трансформировалась в единую ячейку с квадратной планфор-мой, которая заняла, практически, всю расчетную область, С другой стороны, структура тепловой конвекции, вязкость которой зависит только от глубины, представляла собой двухмерные валы, медленно эволюционирующие во времени- И, наконец, был рассмотрен случай одновременной зависимости вязкости и от температуры, и от давления. Структура течения теперь представляла собой суперпозицию случаев зависимости вязкости от температуры и от давления в виде ансамбля ячеек, горизонтальные размеры которых были в два раза больше чем вертикальные и восходящие потоки были окружены нисходящими течениями вещества. Далее автор выполнил тест на влияние аспектного отношения на структуру конвекции. Численные эксперименты на различных размерах расчетной области показали, что структура течений в области 4x4x1 с рефлективными граничными условиями («скользкие» границы) оказалась идентичной структуре в области 8x8x1 с периодическими граничными условиями. Этот тест позволяет сделать вывод о том, что вертикальные границы несущественно влияют на структуру конвекции при ширине области более чем в 4 раза пре восходящей ее вертикальный размер. Увеличение же числа Рэлея с 105 до 106 также принципиально не изменило структуру течения, хотя вертикальные потоки становились чуть уже и появлялись модуляции локальных структур конвекции и нестабильности в горизонтальных тепловых слоях- Введение в рассмотрение зависимости вязкости от температуры, при граничном условии на верхней грани в виде жесткой границы, не изменило в значительной мере ни значения числа Нуссельта, ни среднеквадратичную скорость движения вещества в боксе. Таким образом, данная работа была базовой при выборе размеров расчетной области, граничных условий и параметров конвекции в настоящей работе. Более подробно параметры нашей модели будут представлены в главе «Математическая формулировка задачи».

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

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

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

Необходимо отметить, что большая часть современных работ по 3D конвекции, была посвящена изучению самого процесса конвекции при физических параметрах, отвечающих состоянию недр. Из тех немногих публикаций, где представлены результаты 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, ПО, 111,112, 113, 122, 127, 128, 129, 131, 132, 136].

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

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

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

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

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

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

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

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

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

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

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

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

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

Глава 2 посвящена тестированию численной модели- Описана .геометрия тестовых задач, методы решений, сеточные представления. Представлены результаты расчетов тестовых задач высоковязких течений в прямоугольном параллелепипеде при постоянной и переменной (зависящей от температуры и давления) вязкости. Получено хорошее согласие с результатами теста [73]. Приведены результаты расчетов с применением метода последовательности сеток для течений в кубе при переменной вязкости и для модельной задачи решения уравнения Пуассона, иллюстрирующие его достаточно высокую эффективность. Продемонстрировано повышение точности экстраполяцией по Ричардсону [30],

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

Последовательно рассмотрены следующие пять задач эволюции з мантийного вещества в параллелепипеде (4200x4200x700 км ):

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Необходимо отметить, что большая часть современных работ по 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, ПО, 111,112, 113, 122, 127, 128, 129, 131, 132, 136].

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

Действительно, при существующей концентрации дол гожи-вущих радиоактивных изотопов в недрах, обеспечивающих тепловой режим планеты, без эффективного конвективного выноса радиогенного тепла из области глубиной в 3000 км, температура в ней поднялась бы выше 4000 С за первые сотни тысяч лет и началось бы повсеместное плавление [52],

Конвективное перемешивание вещества, существующее столь длительное время, привело к тому, что в большей части мантии, за исключением граничных тепловых слоев, вертикальный градиент температуры близок к адиабате. Поэтому для изучения влияния неоднородностей литосферы континентов, средний возраст которой составляет 2000 млн. лет, на режим конвекции, в качестве начального распределения температуры и скоростей в расчетной области необходимо использовать распределение этих параметров для развитой субстационарной конвекции. Термином «субстационарный» условимся называть режим конвекции, при котором количество ячеек остается постоянным, а меняются только их размеры. По классификации, представленной в [137] этот режим назван осциллирующим. Наиболее подходящей моделью такой конвекции на наш взгляд является конвекция под плитой постоянной мощности, рассчитанная для «больших» времен, в которой влияние неодно- родностей литосферы по толщине отсутствует, так как отсутствуют сами эти неоднородности- Как показали 2D модели конвекции [56] для выхода на субстационарное состояние конвекции в верхней мантии под кондуктивной литосферой толщины 120 км необходимо время порядка 6000 млн. лет. Структура течения при этом представляет собой вытянутую по горизонтали ячейку с аспектным отношением, доходящим до 10.

Поле температур под кондуктивной литосферой толщины 120 км. Временной срез при t=6000 млн. лет. Максимальное значение области по вертикали -670 км. В трехмерном варианте процесс конвективного перемешивания идет более интенсивно [108]. Проведенные нами эксперименты показали, что по прошествии 2500 млн. лет структура течения устанавливается, количество ячеек при этом остается постоянным. Эти результаты дали основание выбрать в качестве субстационарного состояния структуру течения при t=2500 млн, лет. Таким образом, следует отметить, что переход от 2D к 3D моделированию показал, что учет трехмерности в моделях кардинально меняет не только структуру течения, но и ее эволюционные характеристики, сокращая, в частности, время выхода системы на устойчивый режим более чем в два раза. Прежде всего, кратко опишем эволюцию структуры конвективного течения под литосферной плитой постоянной толщины во времени. Начальный этап развития конвективной неустойчивости соответствует временам порядка первых сотен миллионов лет. На Рис. 19 для данного этапа дано расположение горизонтальных плоскостей XY на различных глубинах, для которых показана геометрия структуры течения в плане.

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

Восходящий поток со всех сторон окружен холодными нисходящими потоками (более темные области) с температурой не более 1200С. Структура ячеек в плане представляет собой пятиугольники, хотя встречаются четырех- и шестиугольники. Количество ячеек при эволюции структуры течения медленно уменьшается и находится для данного этапа в интервале значений 30 - 33. " Характерный горизонтальный размер ячеек колеблется в пределах 500-800 км, таким образом можно говорить об их изомет-ричности, принимая во внимание вертикальный размер конвектирующей области в 575 км. Следующее сечение выполнено на глубине 220 км (Рис 20Ь). При сравнении с предыдущим сечением видно, что структура течения в плане практически не изменилась. Это говорит о субвертикальном положении восходящих и нисходящих потоках в данном диапазоне глубин. Следует отметить, что с увеличением глубины обнаруживается некоторая фокусировка (сужение по толщине) восходящих потоков. Структура течения начинает существенно меняться на глубине 380 км, в центральной части конвектирующей области (Рис. 20с)- О ячеистой структуре течения теперь можно говорить с некоторой долей условности, поскольку некоторые из многоугольников нисходящих потоков оказались разрушенными и, вместе с тем появилась тенденция к объединению ранее изолированных восходящих потоков в цепочки. Эта тенденция еще более выражена в нижней части конвектирующей области, на глубине 560 км (Рис. 20d). Такая эволюция структуры течения в плане говорит о том, что, начиная с центральной части конвектирующей области и глубже, начинает ощущаться влияние нижнего горячего теплового слоя, где распределение температуры конвекции под ршшой литосферой, L 120 км. при с-:- 0 млн лет.

На Рис. 21 дано расположение ляхи вертикальных разрезов Б плоскости YZ7 для которых показана структура течения по вертикали. Представленные разрезы показывают, что вдоль горизонтальных осей расчетной области располагаются 6-7 ячеек расстояние между восходящими и нисходящими вертикальными новиками в которых составляет 400-500 к,м Толщина нижнего теплового граничного слоя не превышает 100 км (Рис. 22).

Результаты 2D моделирования конвекции в верхней мантии, перекрытой сверху жестким кондуктивньш слоем, показали, что е течением Бремени эволюция структуры конвекции идет в направлении сокращения количества ячеек с одновременным удлинением их по горизонтали \5Ц. Та же тенденция наблюдается н в 3D моделях. Рассмотрим структуру течения в исследуемой области для времени t=2500 млн- лет. Как следует из сопоставления Рис, 20 и Рис, 26, где представлены структуры течений в плане для начального этапа конвекти-рования и для времени 2500 млн. лет, количество ячеек со временем сократилось вдвое и теперь вместо 30 ячеек в расчетной области существуют только 16. Увеличился средний размер ячеек, который составляет теперь 1000 - 1500 км. Вместе с тем, характер течения с глубиной практически не изменился; течение так же фокусируются с глубиной и в нижней части области сохранилась тенденция к объединению восходящих потоков в непрерывные цепочки, но теперь они не такие длинные как на начальном этапе конвекции. На вертикальных сечениях ясно видно удлинение горизонтального масштаба ячеек (Рис. 24,25).

Как показали результаты моделирования, начиная с этого времени, t=2500 млн. лет, Рис. 26, структура конвекции практически не меняется. Для иллюстрации этого утверждения рассмотрим результаты моделирования трехмерной конвекции при t=5000 млн. лет, Рис. 27. Как видно из сопоставления Рис. 23?26. и Рис, 27,28, течение в общих чертах сохранило свою форму и количество ячеек.