Бесплатный автореферат и диссертация по наукам о земле на тему
Численные модели мантийной конвекции с переменной вязкостью и фазовыми переходами
ВАК РФ 25.00.10, Геофизика, геофизические методы поисков полезных ископаемых

Автореферат диссертации по теме "Численные модели мантийной конвекции с переменной вязкостью и фазовыми переходами"

РОССИЙСКАЯ АКАДЕШМ НАУК ИНСТИТУТ ФИЗИКИ ЗЕМЛИ им. О.Ю. Шмидта

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

ЕВСЕЕВ Александр Николаевич

ЧИСЛЕННЫЕ МОДЕЛИ МАНТИЙНОЙ КОНВЕКЦИИ С ПЕРЕМЕННОЙ ВЯЗКОСТЬЮ И ФАЗОВЫМИ ПЕРЕХОДАМИ

Специальность 25.00.10 - Геофизика, ■ геофизические методы поиска полезных ископаемых

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

2 з О ИТ 2008

Москва-2008

Работа выполнена в институте физики Землй Российской Академии наук, г. МоСКВа. = . :

Научный руководитель; член-корреспочаент Р АЦ

доктор физико-математических наук, Валерий Петрович Трубицын, ' -1 Институт физики Земли Р АН

Офшдиальныеоппонентыгчлен-керреепондентРАН,"

доктор химических наук, Олег Львович Кускоц

Институт геохимии и аналитической химии Р АН

доктор физико-математических наук, Алексей Леонидович Собисевич, Институт физики Земли РАН

Ведущая организация: Международный институт

теории прогноза землетрясений и математической геофизики РАН

Защита состоится 13 ноября 2008 г. в 14 ч, 00 на заседании диссертационного совета Д 002.001.01 Института физики Земли Российской Академии Наук (ИФЗ РАН) по адресу: 123995, г. Москва, ул. Большая Грузинская 10

С диссертацией можно познакомиться в библиотеке ИФЗ РАН Автореферат разослан «¿3» октября 2008 г.

Ученый секретарь

диссертационного совета /

кандидат физико-математических наук о.В: Пилипенко

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

Актуальность работы

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

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

Цели и задачи работы

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

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

Основные положения работы, выносимые на защиту

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

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

3. Найдено количественно, в какой мере влияют на конвекцию все основные известные в мантии фазовые переходы на глубинах 410 км, 550 км, 660 км, 1500 км, 2700 км. Фазовый переход на границе верхней и нижней мантии (660 км) тормозит конвекцию, уменьшая тепловой поток из мантии на 4-5%, однако

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

Методика исследований

1. Основным методом исследований являлось численное моделирование.

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

3. Разработана программа на языке программирования Си в среде Cygwin, позволяющей проводить расчеты на ПК в среде Windows по программам для Unix.

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

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

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

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

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

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

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

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

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

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

Теоретическая и практическая значимость

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

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

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

Результаты исследований докладывались и обсуждались на научных семинарах Института физики Земли РАЦ Международной конференции EGU General Assembly,Vienna, 13-18 April 2008 и на Девятом международном совещании "Физико-химические и петрофизические исследования в науках о Земле", г. Москва(ИФЗ РАН ГЕОХИРАН), 8 октября 2008 г.

Структура и объем диссертации

Диссертация состоит из введения, 6 глав, заключения и списка литературы. Общий объем диссертации составляет 102 страницы машинописного текста, содержит 40 рисунков и список литературы из 100 наименований.

Выполнение работы

Результаты, изложенные в диссертации, получены автором во время учебы в аспирантуре и в ходе работы в должности научного сотрудника Института физики Земли РАН в период 2006-2008 гг. Диссертация выполнялась при поддержке грантов РФФИ 05-05-64029, 07-05-01002, гранта Президента РФ государственной поддержки научных школ НШ-1901.2003.5 и грантов "Лучшие аспиранты РАН' Фонда содействия отечественной науки.

Личный вклад автора

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

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

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

Автор глубоко признателен своему научному руководителю члену-корреспонденту РАН В.П. Трубицыну за постоянное внимание, помощь при выполнении работы и плодотворные обсуждения по теме исследований. Автор также благодарен всем сотрудникам ИФЗ РАН, особенно Баранову A.A., Боброву A.M., Трубицыну А.П. и Харыбину Е.В. за дружеское внимание, поддержку и полезные обсуждения во время работы над диссертацией.

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

Введение

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

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

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

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

Данная гипотеза держалась в течение двух десятилетий. Анализу этой гипотезы и ее геофизических, геохимических и геологических следствии были посвящены сотни работ. Однако новые данные лабораторных измерений параметров фазового перехода (2003-2005 гг.) поставили под сомнение гипотезу о возможности расслоения конвекции в мантии Земли. В диссертации развита методика количественной оценки влияния фазовых переходов на конвекцию. Были рассчитаны модели конвекции для всего возможного диапазона их значений и найдены области однослойного, переходного и двухслойного режимов конвекции. Данные новых лабораторных измерений вещества мантии на рассчитанной диаграмме попадают в область общемантийной конвекции без расслоения и без аваланчей.

В пиролитовой (ПИРоксен-ОЛивин) модели Рингвуда верхняя мантия состоит из перидотита (90%) и эклогита (10%). Перидотит является ансамблем минералов: 60% оливина, 25% ортопироксена, а также клинопироксена и граната. Эклогит содержит клинопироксен (60%), гранат (30%) и другие малые добавки. Нижняя мантия состоит на 75% из силикатного перовскита (М§0 9ре<н)8Юз, на 15% магнеовюстита (М£о8Реог)0 и 10% кальциевого перовскита (Са81)03.

С ростом давления и температуры вещество мантии испытывает фазовые превращения (см. рис.1). Поскольку вещество мантии многокомпонентное, то даже в оливине переход происходит не скачком, а размазан на некоторый интервал давлений Ар или глубин (1=Ар/(р§). В мантии на глубине 410 км оливин превращается в вадслеит со скачком плотности 5р/р0=0.07 при наклоне кривой фазового равновесия ур=(1р/сГГ= 1.6-3.0-2.3 МПа/К и ширине фазового перехода с!=13 км. На глубине 520-580=550 км вадслеит превращается в рингвудит со скачком 8р/ро=0.03 при ур=4.3 МПа/К и ширине фазового перехода до <1=60 км. На глубине 660 км рингвудит превращается в смесь перовскита и магнеовюстита с большим скачком плотности 8р/ро=0.09 при отрицательном наклоне кривой фазового равновесия (см рис.1). Его величина по последним данным лабораторных измерений оказывается существенно меньшей (Катсура и др. 2003, 2005; Литасов и др. 2005), чем обычно принималось в большинстве работ и составляет лишь ур=с1р/с1Т=-(0.5-2.0)~-1.3 МПа/К при ширине с!~3 км. В интервале глубин 1250-1750 км железо в перовските переходит из состояния с высоким спином электрона в состояние с низким спином, при этом плотность меняется на 5р/ро~0.01. В этом переходе также происходит перераспределение железа и магния между перовскитом и окисью магния. Наклон кривой равновесия оказывается достаточно высоким Ур=(1р/сГГ~16 МПа/К, и переход имеет ширину до (1—500 км. Наконец, на глубине около 2700 км (на верхней границе слоя Б") перовскит переходит в пост-перовскит при 5р/ро=0.02, ур~8 МПа/К и ширине перехода (1-40 км. В Земле еще имеется фазовый переход базальта в эклогит на глубинах 60-100 км. Однако в отличие от рассматриваемых переходов он происходит не во всем веществе мантии, а только в океанических базальтах опускающейся литосферной плиты, т.е. только в одной компоненте. В диссертации приводятся

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

Т°С

2400-

400

600

глубина

(км)

давление

10 15 20 25 30(ГПа)

Рис. 1 Фазовая диаграмм для оливина, пунктирная кривая - распределение средней температуры мантии.

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

Фазовый переход, в отличие от химического перехода, несмотря на скачок плотности, может или не влиять на конвекцию, или тормозить ее, или ускорять. Области устойчивости состояния фаз (см. рис. 1) разделяются кривой фазового равновесия р=р(Т). Поскольку в мантии давление зависит от глубины р=р§Ъ, то слои с разными фазами разделяются границей Ь=Ь(Т). Если наклон фазовой кривой ур=с1р/с1Т в переменных давление-температура или уь=сШ/сГГ=ур/р§ в переменных глубина-температура равен нулю, то граница перехода не зависит от температуры и находится на одной и той же глубине в восходящем и нисходящем мантийном потоке. Пересекая фазовую границу, вещество мгновенно (если не учитывать малое время кинетики процесса) переходит в новую фазу. Поэтому вес столбов нисходящего и восходящего мантийных потоков не изменяется при таком фазовом переходе. Заметим, что в этом случае фазовый переход происходит без выделения и поглощения тепла, т.к. наклон фазовой кривой связан с теплотой выделяющейся при фазовом переходе, соотношением Клапейрона-Клаузиуса q=yp5pт/(plp2), где р) и рг -плотности фаз, и 6р=р2-р1 и ур=<1р/сГГ.

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

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

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

Влияние фазового перехода на мантийную конвекцию изучается уже много лет, начиная с работ Христенсена 1984-1986 гг. Результаты двумерного численного моделирования показали, что при отрицательном наклоне фазовой кривой при значениях модуля |ур|>2.5-3.0 МПа/К торможение вертикальных мантийных потоков приводит к переходному режиму конвекции (при котором конвекция попеременно расслоена во времени и пространстве), а при более высоких значениях |ур| затем происходит глобальное расслоение конвекции. До 2002 г. лабораторные измерения оливина давали значение наклона ур в диапазоне от -2 МПа/К до -2.5 МПа/К. Эти значения по модулю несколько меньше критического, но находятся на грани точности измерений и расчетов. Поэтому гипотеза о расслоении конвекции интенсивно обсуждалась. После 2002 г. была существенно улучшена методика измерений, т.к. фазовый переход стали фиксировать непосредственно во время опыта, а не по остаточным образцам после снятия давления и температуры. Эти измерения дали меньшее по модулю значение наклона ур от -0 Д° -2 МПа/К. Поэтому гипотеза о возможности глобального расслоения конвекции в мантии Земли стала пересматриваться. Но поскольку фазовый переход с отрицательным наклоном все-таки тормозит конвекцию и появились возможности более точного моделирования конвекции, то стала актуальной проблема более детального расчета влияния фазовых переходов, причем как эндотермических, так и экзотермических.

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

Глава I. Уравнения мантийной конвекции с переменной вязкостью и Фазовыми переходами

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

имеют вид:

ЭУх/Эх +дУг/дг =0, (1)

0= до^/дх+ до^/дг, (2)

до 2Х/Зх+ Эо ¡в/Э^-ЛаТ+КрьДГ, (3)

ЭТ/а+УхЭТ/йх+ У2ЭТ/дг= ^Т/д^+д^Т/дг? , (4) где

ст^-р+гг^Ух/дх, о22=-р+2т!аУ2 /дг. (5)

ат= ц(д\х/дг +дУ2/дх), (6)

Ка=аоР(£ТоВ3/(кпо), ЯРь= 5рёВ3/(ктю), КР= Ирь^а = (8р/ро)/(аоТо) (7)

Здесь Ла и Ырь - тепловое и фазовое числа Рэлея, Ир - плотностное отношение, к - коэффициент температуропроводности. За единицы измерения приняты: толщина слоя О - для длины, У0=к/О - для скорости, 1о=Б2/к - для времени, г|о - для вязкости и Т0 - для над адиабатического перепада температуры.

Основные неизвестные функции V - скорость течений, Т - температура, р - давление находятся из решения трех уравнений. Компоненты тензора напряжений а,] находятся дифференцированием скоростей вязких течений. Вязкость г] считается заданной функцией координат, давления, температуры и напряжений. Ось г считается направленной вверх. Функция ДГ=Т(х,г)-Го отлична от нуля только в области между фазовой границей и уровнем Ь=Ь0. При этом она равна ДГ=1 в области, где фазовая граница поднята и равна ДГ=-1 в области, где фазовая граница опущена. Если фазовых переходов несколько, то в правой части уравнения (3) стоит суммаX ЯрыАГ,.

При заданной фиксированной глубине Ьо фазового перехода система уравнений конвекции (1-7) содержит три безразмерных параметра: Ыа, Ярь и у, входящий в ДГ. Тепловое число Рэлея 11а, характеризующее интенсивность тепловой конвекции, и фазовое число Рэлея характеризующее влияние скачка плотности при фазовом переходе, являются комбинацией параметров: коэффициента теплового расширения осо, толщины слоя Б, температуропроводности к и средней вязкости Т10.

Безразмерный наклон кривой фазового равновесия у, входящий в выражение для функции распределения фаз Г, характеризует величину смещения фазовой границы. Он связан с размерным наклоном в переменных глубина-температура и в переменных давление-температура соотношениями у=у„То/0=урТо/р8В (8)

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

латеральных вариаций температуры Т(х,г*). Поэтому приращение ДГ можно записать как сложную функцию ДГ=(дГ/с^)(Э^'ЭТ)ДТ. Учитывая, что в уравнениях конвекции через Т фактически обозначена поправка к температуре за счет конвекции ДТ=Т, можно переписать правую часть уравнения Стокса (7) в виде:

-КаТ+Ярь-ДГ = -К.аТ+КрЬ(ЗШОдеТ)Т = -Ка(1+Рдт0-Т = -Яае.Т (9) где Яае{ = Яа(1+РЗГ/о^) - эффективное число Рэлея с учетом фазового перехода и Р - фазовый параметр, характеризующий влияние фазового перехода на конвекцию. С учетом (4) и (11) его можно записать в виде: Р = -Крь(9С/ЭТ)Л1а = -Кр(5С/ЭТ) = 1^/ (10)

Таким образом, уравнение Стокса (3), входящее в систему уравнений конвекции (1-7) с учетом фазового перехода в безразмерных переменных может быть записано как в форме (3), выраженной через параметры и у и функцию Г, так и в форме, выраженной через параметр Р и производную дГ/дС,: до2Х/Зх+ ао22/& = -Ка(1+РЭШ0 (11)

Безразмерный параметр фазового перехода Р,являющийся произведением скачка плотности фазового перехода 8р/ро и наклона кривой равновесия у, характеризует избыточный вес столба нисходящего или восходящего потока, влияющий на скорости течений.

Поскольку вещество мантии является смесью (твердым раствором) компонент, то фазовые переходы в мантии осуществляются не скачком, а размазаны на интервал давлений или глубин 81г = с1. В этом интервале плавно уменьшается содержание одной фазы и увеличивается содержание другой. Для аналитического описания функцию распределения фаз Г(^(х,г)) можно представить в виде Г(х,г) = 1(^/(1). Здесь (1 - ширина зоны фазового перехода и - какая-либо ступенчатая функция, например, функция гиперболического тангенса = 0.5(1+1Ц) или экспоненциальная функция вида: Щ = 1/(1+ехр(-44)).

Функция распределения фаз Г равна 0.5 в середине зоны фазового перехода, а на краю зоны перехода (при значении ^=Ь-Ь(х)=(1/2) уменьшается до значения Г(0.5)~0.1, что соответствует содержанию второй фазы 10%. Производная дГ/д£„ входящая в уравнение Стокса (11), будетравна

= (12) где Г©= 4ехр(-4^)/[1+ехр(-4^)]2=4[ехр(24)+ехр(-2^)]"2 (13)

Функция дГ/дС, принимает максимальное значение, равное 1/бЬ, при §=0, т.е. в точках в середине зоны фазового перехода. При удалении от этой границы она экспоненциально стремится к нулю как ехр(-2^с1).

В общем случае при произвольных вариациях вязкости и функциях тепловых источников эта система уравнений решается только численно. Для анализа решений и тестирования численных решений были использованы точные аналитические частные решения уравнения Стокса, которые можно найти, например, следующим простым способом. Сначала подберем какие-либо скорости V, удовлетворяющие граничным условиям, с неизвестными коэффициентами. Затем дифференцированием найдем вязкие напряжения. Неизвестные коэффициенты найдем из уравнений (1) и (2). Затем из (3) найдем

температуру. В результате все функции в (1-3), (5,6) станут известными, а исходные скорости можно рассматривать как точное специальное решение уравнения Стокса при некоторой специально заданной температуре при произвольном распределении вязкости.

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

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

Глава II. Разработка программы численного моделирования двумерной мантийной конвекции методом конечных элементов

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

Л [01 <14)

к в вт о

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

Вязкость на каждом элементе задается в каждой точке интегрирования интерполированием поля температур внутри элемента. Значения элементов матрицы К находились численным интегрированием по Гауссовой схеме по сетке 2x2 на каждом элементе.

Для решения матричного уравнения (14) был применен алгоритм У завы. Комбинируя уравнения и исключая V, получим уравнение для давления Р в виде (0т-К"1-0)-Р=0т-К"1-Р. Матрица У=От-К"'-С симметрична и положительно определена. На практике давление Р не получают напрямую из уравнения из-за сложностей связанных с получением К"1, поэтому для получения коррекции давления был использован алгоритм метода сопряженных градиентов который не требует вычисления матрицы У"1.

Зададимся малой сжимаемостью г. На нулевом шаге выберем давление Ро = 0. Скорости на первом шаге получим как решение уравнения КУ0 + ОРо = КУ0=Р, т.е. УрЬТ^. Начальная невязка г0 будет равна Г0=О1У0=ОтК"1Р. Начальное направление поиска в) выберем равным 81=Го-

Давление на к-ом шаге уточняется по формуле Рк =Рк.1+акзк, скорость - по формуле Ук=Л/1<.1-акК"'05к=Уы-акик. Согласно методу сопряженных градиентов, параметры ак и зк на к-ом шаге выбираются как:

8к=Гк-1+гтк-1Гк-1/(гтк-2Гк-2)зк-1 и ак=гтыгк.1/(08к)тК"|08к, где гы - невязка с предыдущего шага. Для упрощения вычислений параметра ак сначала решаем

уравнение Ки^ОЭк относительно ик. Это дает возможность быстро вычислять

Т Т 1 Т 1 т

произведения: в ^кНОвк) К" 08к=(05к) ик, акК" 05к=акик и акО К' 1С8к=а1сОтиь Процесс останавливается, когда текущее значение нормы невязки гк становится меньше параметра е.

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

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

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

Глава Ш. Влияние эндотермического фазового перехода на структуру тепловой конвекции и массобмен между слоями.

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

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

На рис.3 представлена рассчитанная структура конвективных переходов при Ка=105 при двух значениях фазового параметра (10), равных Р=-0.3 и Р=-0.44. Вверху даны распределения температуры Т(х,г) и поля скоростей У(х.г). Внизу - зависимости от глубины (\\=1-г) для температуры Т(г) и модуля вертикальной скорости у=|\г2(г)|, усредненные по латерали. Поток массы вещества через границу 2=г* в каждой ее точке равен рУ2(х,г*). В рассматриваемом приближении Буссинеска не учитывается сжимаемость вещества. Поэтому массобмен характеризуется самим значением вертикальной скорости У2(х,г*). При этом поток вещества через границу является функцией

х. Вертикальная скорость течений У2(х,2) равна нулю на верхней и нижней границе и без фазового перехода достигает экстремума на полувысоте. Поскольку скорость У2(х,г) положительна в восходящем мантийном потоке и отрицательна в нисходящем мантийном потоке, то ее среднее значение равно нулю. Поэтому скорость перемешивания мантии через рассматриваемую границу можно характеризовать усредненным по латерали модулем вертикальной скорости у=|У2(г)|. Как видно из рис.3, при значениях фазового параметра Р = -уЯр, по модулю меньших 0.3, конвекция остается однослойной, т.е. охватывает всю область, и течения свободно проходят через фазовую границу. Эффект фазового перехода определяется всего одним параметром Р, зависящим от значения произведения убр. Поэтому, несмотря даже на большое значение наклона кривой фазового равновесия у, фазовая граница может не влиять на структуру конвективных течений, если мал скачок плотности. При значении фазового параметра |Р|>0.44 конвекция становится расслоенной, и конвективные течения в верхнем и нижнем слое замыкаются внутри своей области практически без обмена веществом.

Рис. 3 Рассчитанные распределения температуры и конвективных скоростей для 1*а=105 при однослойной и расслоенной конвекции при значениях фазового параметра, соответственно равных Р=-0.3 (слева) и при Р=-0.44 (справа). Фазовая граница показана черной линией. Безразмерная температура показана серыми тонами, безразмерные скорости течений - стрелками с максимальным значением Ут=340. Внизу приведены усредненные по горизонтали распределения температуры Т(х) и модуля вертикальной скорости |Уг(х)|.

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

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

На рис, 4 приведена структура конвективных течений также при Ка=10 дня промежуточного значения фазового параметра Р=-0.4 для двух моментов безразмерного времени 11=1.00 и 12=1.02.

Рис. 4 Рассчитанные распределения температуры и конвективных скоростей для К;1"" 1 !Г , аналогичные рис. 3, но для переходного режима конвекции при значении фазового параметра Р=-0.4 для двух моментов времени.

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

Как видно из рис.4, в момент ^=1.00 степень расслоенности невелика, при этом часть конвективных течений замыкается в нижнем слое, а часть проникает в верхний слой. На кривой зависимости среднего модуля вертикальной скорости на фазовой границе у=|У2(г)| от глубины при г~т* имеется небольшой провал. В момент 12=1.02 конвекция становится более расслоенной.

На рис.5 приведены кривые зависимости от значения фазового параметра |Р| для у^^У^г*^/^^*,?^)!, характеризующего относительную степень расслоения конвективных течений, а также для амплитуды 8у1* и периода колебаний т при переходном режиме конвекции. В отношение V]* входит величина У2(г*,Р=0), равная вертикальной скорости течений через фазовую границу при нулевом наклоне фазовой кривой, т.е. когда фазовый переход не оказывает влияние на конвекцию.

Как ввдно из рис.5, в интервале значений параметра |Р| от 0.34 до 0.44 массообмен через фазовую границу и степень расслоенности У1*(Х) во времени меняется. Амплитуда значений У1*(Т) в левой части рис.5 показана вертикальными отрезками. В средней части рис.5 приведена рассчитанная амплитуда колебаний 8у[* в зависимости от фазового параметра |Р|. В правой части рис.5 приведена вычисленная зависимость периода колебаний VI* при переходном режиме конвекции. Проведенные расчеты показывают, что при рассматриваемом числе Рэлея 11а=105 период колебаний достаточно плавно возрастает с ростом значений фазового параметра |Р|.

0.024

0.016-

О 0.2 0.4 0.6 0.8

0.34

0.38

л 0.008-

0.42

0.34

0.38

т*

0.42

1*1 |Р| |Р|

Рис. 5 Рассчитанный относительный средний модуль вертикальной скорости на фазовой границе VI*, а также амплитуда и период его колебаний во времени в зависимости от значения модуля фазового параметра |Р|.

V?

Уг

Р=-0-34

—I-1-г-1-г

О Й,1 0,2

t

Р=-0.36

I

V*

р—0.344

—I-!

а2

I

Р-0.44

—г чг

{

Рис. 6 Рассчитанные значения относительного среднего модуля вертикальной скорости на фазовой границе V!* как функции времени для двух значений фазового параметра Р=-0.34 и Р=-0.44, соответствующих крайним значениям области переходного режима конвекции и значениям Р=-0.344 и Р=-0.36 для внутренней области.

На рис.6 представлены рассчитанные значения относительного среднего модуля вертикальной скорости У[*(1) в зависимости от времени при различных

значениях фазового параметра |Р|. При изучении влияния фазового перехода на конвекцию сначала была рассчитана структура течений, устанавливающаяся при рассматриваемом числе Рэлея Яа=103 без фазового перехода Затем включался фазовый переход с конкретным значением фазового параметра Р. При резком изменении параметров расчитываемой конвекции всегда возникают колебания структуры конвективных течений. Однако такие колебания не связаны с собственными колебаниями системы и затухают со временем. Как видно из рис.6, при значении параметра |Р|<0.34 и |Р|>0.44 рассчитанные колебания среднего модуля вертикальной скорости на фазовой границе затухают со временем, а в интервале 0.34<|Р|<0.44 колебания продолжаются бесконечно долго, что соответствует переходному режиму конвекции. Амплитуды этих колебаний и их периоды, приведенные на рис.5, были построены на основании зависимостей У1*(Т), рассчитанных для большого диапазона значений фазового параметра Р.

11 Т IVZJ " " Т ™ JVI! ° "т Ж IVit

Рис. 7 Рассчитанные распределения температуры и конвективных скоростей для Ra=lü6 и значений фазового параметра, равных Р=-0.15, -0.2 и -0.4.

На рис.7 приведены аналогичные распределения, но при большем значении числа Рэлея Ra=106 и меньших значениях фазового параметра, равных Р=-0.15, -0.2 и -0.4. Эти значения параметра Р выбраны, чтобы показать различные режимы конвекции. При |Р|<0.15 фазовый переход почти не влияет на конвекцию. При Р=-0.2 структура конвективных течений также почти не меняется, но начинает уменьшаться интенсивность конвекции. Как видно из рис.7, при Р=-0.15 максимальная вертикальная скорость (в безразмерных единицах) равна 640, а при Р=-0.2 она уменьшается до значения 580. При Р=-0.4 возникает режим расслоенной конвекции, при котором в нижней мантии устанавливается одна деформированная конвективная ячейка, а в верхней мантии — две ячейки в соответствии с условием их энергетически более выгодной квазиквадратной формы. При этом течения в левой ячейке верхней мантии более интенсивны, т.к. их направление против часовой стрелки скоординировано с течениями по часовой стрелке в нижней мантии. Течения по часовой стрелке в правой ячейке верхней мантии более слабые. Благодаря большому трению между встречными скоростями на фазовой границе

тормозятся течения в правой части рассматриваемой квадратной области, как в верхней, так и в нижней мантии. При расслоенной конвекции при Р=-0.4 распределение температуры имеет характерный скачок в окрестности фазовой границы, который представляет собой сдвоенный тепловой (кондуктивный) погранслой конвекции нижней и верхней мантии. На кривой распределения модуля вертикальной скорости |У2(г)| в области фазовой границы возникает провал, соответствующий уменьшению массообмена между нижней и верхней мантией.

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

Рис. 8 Эволюция ноля скоростей и температуры в переходном режиме конвекции при Иа^Ю6 и Р=-0.22.

На рис.8 приведены результаты расчета эволюции мантийных течений при Яа=106 и значении фазового параметра Р=-0.22, соответствующему переходному режиму конвекции. Распределения температуры и скоростей течений приведены для четырех значений безразмерного времени ^=0, \.2~1'1(У4, 1.6-10"3 и 14=2.5• 10~3.

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

является неустойчивым и в момент 12=7-10" вещество верхней мантии прорывается в нижнюю. Возникает аваланч, который увлекает за собой значительное количество вещества верхней мантии. Скорости течений и в верхней и в нижней мантии увеличиваются в полтора-два раза. При этом аваланч создает в нижней мантии две ячейки конвекции. Благодаря обмену веществом между горячей нижней мантией и холодной верхней мантией температура верхней мантии повышается, ее плотность уменьшается, и аваланч прекращается. Как видно из рис.8, средняя температура верхней мантии после аваланча может повышаться на 10%, что в размерных единицах соответствует 250 К. К моменту 13=1.6-10"' опять наступает относительно спокойная фаза расслоенной конвекции. Однако в нижней мантии остаются две конвективные ячейки. В верхней мантии после аваланча течения остаются настолько интенсивными, что воздействуют на нижнюю мантию, усиливая течения в ее правой ячейке. После момента времени 14=2.5-10"3 формируются условия для возникновения второго аваланча

1=0

1=1.2-10

Рис. 9 Эволюция поля скоростей и температуры в переходном режиме конвекции при Ка~Ш и Р=-0.18

На рис.9, аналогично рис.8, приведены результаты расчета эволюции мантийных течений, но при еще более высоком значении числа Рэлея Яа=107 и при значении фазового параметра Р=-0.18, также соответствующему переходному режиму конвекции. Распределения температуры и скоростей течений приведены для четырех значений безразмерного времени ^=0, 12=1.2-10ч, 1з=8 .0-10"4 и 14=2.2-10"3.

Как видно из рис.9, при значении числа Рэлея 11а=107 прорывы вещества из верхней мантии в нижнюю становятся более мощными, и соответственно во

время аваланча резко увеличивается перемешивание всей мантии. Повышение температуры верхней мантии после аваланча может достигать 20% или 500 К.

Для количественной характеристики степени расслоения мантийных течений, вызванного эндотермическим фазовым переходом, использались две величины. Первая величина была введена ранее и равна v*i=|V2(z*)|/|Vz(z*,P=0)|. Отношение v*i показывает, в какой степени фазовый переход изменяет массообмен между верхней и нижней мантией. Вторая величина v*2, равная отношению потока массы через фазовую границу к усредненному по вертикали потоку массы: v*2= |Vz(z*)]/<]V2(z)|>, где z*-координата фазовой границы.

Величина v*2 использовалась ранее в работе (Solheim and Peltier, 1994), но она является просто нормированной скоростью. В то время как введенная в диссертации величина v*i характеризует изменение степени массобмена, вызванное фазовым переходом

X1

02 04 Об

Рис. 10 Степень расслоения конвективных течений у*2 (сплошная линия) и относительный массообмен У, (пунктирная линия) в зависимости от фазового параметра при различной интенсивности тепловой конвекции. Точками показаны рассчитанные значения.

На рис.10 приведены рассчитанный относительный массообмен у*1 и степень расслоения у*2 при конвекции с эндотермическим фазовым переходом в зависимости от фазового параметра Р для значений числа Рэлея 11а от 104 до 107. Как видно из рис.10, обе кривые имеют сходное поведение, поскольку расслоение и массообмен определяются величиной усредненного по горизонтали модуля вертикальной скорости течений. По определению относительный массообмен у*1 нормирован на единицу при Р=0. По кривым зависимостей у*|(Р) и у*2(Р) от фазового параметра Р можно выделить три режима конвекции: общемантийный режим конвекции с торможением течений в верхней мантии, переходный режим конвекции и двухслойный режим конвекции. При этом с ростом числа Рэлея границы режимов становятся резче

и ширина зоны переходного режима конвекции уменьшается. Зависимость у*1(Р) показывает, что уже до расслоения при общемантийной конвекции фазовый переход уменьшает массобмен между верхней и нижней мантией. Это обусловлено тем, что при эндотермическом фазовом переходе течения из нижней мантии меньше проникают в верхнюю мантию. Причину ослабления общей интенсивности конвекции от фазового перехода можно понять также из выражения для эффективного числа Рэлея Ка,,р11а( 1+РЗГУЗО, которое уменьшается при Р<0.

Расслоение конвективных течений и смена режима конвекции проявляются также в значениях относительного теплового потока, выходящего через верхнюю границу, (числа Нуссельта N11) и в значении среднеквадратичной скорости мантийных течений УГШ8. Рассчитанные зависимости №(Р) и Утк(Р) от величины фазового параметра приведены на рис.11.

Из рис. 10 и 11 также видно, что с ростом числа Рэлея сужается область переходного режима. Для Ыа^Ю4 при изменении фазового параметра Р от 0 до -1 число Нуссельта № плавно меняется от 5 до 1.6. Для числа Рэлея 11а=105 в диапазоне изменений Р от -0.3 до -0.45 появляется область переходного режима конвекции, в которой число Нуссельта № меняется от 10 до 5. Для числа Рэлея К.а=106 переходная область резко сужается и находится в диапазоне изменений Р от -0.2 до -0.22, в которой число Нуссельта N11 меняется от 20 до 10. Для числа Рэлея Яа=107 переходная область находится в диапазоне изменений Р от -0.159 до -0.16, в которой число Нуссельта Ми меняется от 40 до 20.

Рис. 11 Рассчитанные зависимости от фазового параметра для относительного теплового потока №(Р) и среднеквадратичной скорости конвекции Уг[т(Р) при значениях числа Рэлея от 104 до 107.

В этом переходном режиме конвекции происходят колебания структуры конвективных течений и всех связанных с ними характеристик. На рис.11 приведены рассчитанные числа Нуссельта № (относительного поверхностного теплового потока) как функции времени для различных режимов конвекции при значении числа Рэлея от 104 до 107. При всех числах Рэлея с увеличением модуля фазового параметра |Р| уменьшается интенсивность конвекции и соответственно уменьшается число Нуссельта

При низком значении числа Рэлея 11а=104 не только расслоение, но и переходный режим конвекции практически не возникают. Поэтому число Нуссельта практически не меняется со временем. При Ка=105 возникают все

три режима. Как указывалось при обсуждении рис.3, переходный режим конвекции при Яа=105 ограничен значениями параметра Р от -0.3 до -0.45. Как видно из рис.12, число Нуссельта почти не меняется со временем при общемантийной конвекции и расслоенной конвекции, а в переходном режиме конвекции возникают колебания числа Нуссельта, спектр которых может иметь одну или две моды.

Ва=10"

Ра=105

-р=оо

- Р= -0.2 ""

- Р= -0 5

- Р= -1 0

0.04

Ва=106

\TVWVYVWVVV р= '°-42

-Р= -0.44

01 02

- Р=0 0

г Р= -0 202 —,длл р= -О 5

(Ча=107

002 0 04 » 006

Рис. 12 Эволюция теплового потока из мантии при различных значениях фазового параметра и чисел Рэлея

При Яа= 106 переходный режим конвекции находится в более узком диапазоне изменений параметра Р от -0.2 до -0.22. В этом режиме спектр колебаний, связанный с эндотермическим фазовым переходом, имеет много мод с различными периодами. При этом колебания также имеют место и в режиме двухслойной конвекции (например, при Р =-0.5), но они более регулярные, чем в переходном режиме конвекции. При еще большем числе Рэлея 11а=107 спектр колебаний в переходном режиме конвекции становится очень нерегулярным и появляются моды с острыми пиками, соответствующие лавинообразным прорывам вещества (аваланчам) из верхней мантии в нижнюю и обратно. При этих аваланчах мгновенные значения теплового потока из мантии, несмотря на тормозящий эффект эндотермического фазового перехода, могут быть даже большими, чем в отсутствие фазового перехода

Как отмечалось выше, с ростом числа Рэлея расслоение конвективных течений начинается при меньших значениях модуля безразмерного фазового параметра Р. Для количественной характеристики этого эффекта по рассчитанным моделям была построена диаграмма степени расслоения у*2 в зависимости от фазового параметра Р и числа Рэлея Яа.

Приведенные на рис.13 кривые соответствуют условию у*2(Ка, Р)=соп81 Сплошные кривые проведены по результатам рассчитанных моделей.

Пунктирная кривая предложена в работе Христенсена 1985 г. как кривая критических значений фазового параметра Рсг, разделяющих режимы общемантийной и расслоенной конвекции. Она аппроксимируется формулой Pcr=-4.4-exp(-0.2Ra). Как видно из рис.13, эта кривая близка к рассчитанной кривой v*2(Ra,P)=0.5, соответствующей степени расслоения 0.5 для рассматриваемых в настоящей работе двумерных моделей в квадратной области с постоянной вязкостью. При больших числах Рэлея расслоение конвективных течений наступает при меньших значениях модуля фазового параметра Р, т.е. при меньшем скачке плотности 5р/ро фазового перехода или при меньших значениях модуля наклона кривой фазового равновесия у. При этом область значений фазового параметра, соответствующая переходному режиму конвекции, сужается с ростом числа Рэлея.

1дЯа

Рис. 13 Диаграмма степени расслоения \*г конвективных течений в зависимости от значений числа Рэлея Иа и безразмерного фазового параметра Р.

Выводы к гл.111

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

Глава IV. Тепловая конвекция при различной ширине фазовых переходов

Расчеты структуры конвективных течений проводились на модели нагреваемой снизу вязкой жидкости в вытянутой по горизонтали области при интенсивности конвекции, характеризуемой значениями числа Рэлея, равными 105 и 106. Их численное решение позволяет найти распределение температуры и скорости при конвекции. По этим распределениям можно затем рассчитать число Нуссельта, характеризующее эффективность выноса тепла в жидкости по отношению теплового потока на поверхности при конвекции и без конвекции Ш=(-ЭТ/&+У2Т)/(ДТ/Е)), где Б - толщина слоя жидкости и ДТ - перепад

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

у=—С. 13 -=-0.25

Рис. 14 Рассчитанная структура конвективных течений с фазовым переходом при К а 1(Г и К 2 при различных значениях ширины перехода (.1 = 0.025, 0.125, 0.25. Левая колонка для у=-0.18, правая для у=-0.25. Относительная температура Т показана серыми тонами. Безразмерные скорости течений показаны стрелками при максимальном значении У=310для т=-0.18 и У=240 для -у=-0.25

На рис. 14-16 представлены результаты расчета структуры мантийных течений с фазовыми переходами при интенсивности конвекции, характеризуемой безразмерным числом Рэлея Яа=10:>. В проведенных расчетах плотностное отношение принималось равным Яр=2, что примерно соответствует фазовому переходу на глубине 660 км.

На рис.14 представлены рассчитанные поля температуры и скоростей при конвекции для двух значений безразмерного наклона фазовой кривой у=-0.18 (левая колонка) и у=-0.25 (правая колонка) при различных безразмерных ширинах эндотермического фазового перехода, равных <1=0.025, 0.125 и 0.25. При значении наклона у=-0.18 узкий фазовый переход (<1=0.025) частично тормозит конвекцию. При увеличении ширины фазового перехода торможение уменьшается, а очень размазанный переход при <1=0.25 почти не влияет на конвекцию. При значении наклона у=-0.25 узкий фазовый переход вызывает значительное расслоение конвективных течений. Однако при увеличении ширины фазового перехода и в этом случае его влияние на конвекцию резко ослабевает.

На рис.15 представлены рассчитанные поля температуры и скоростей при конвекции для случая без фазового перехода (слева) и при экзотермическом фазовом переходе со значением наклона у=0.18 (справа) при небольшой

ширине перехода (1=0.025. Как видно из рис.15, при экзотермическом фазовом переходе граница в нисходящем потоке поднимается, что ведет к его утяжелению. При этом конвекция ускоряется, и максимальная скорость течений возрастает с У0=36О до значения Уг=470.

Рис. 15 Рассчитанная структура конвективных течений при Ка=10э без фазового перехода (слева) и с экзотермическим фазовым переходом Кр=2 при у=0.18 и (1=0.025 (справа). Относительная температура Т показана серыми тонами. Безразмерные скорости течений показаны стрелками при максимальном значении У=470.

м ы

Рис. 16 Рассчитанные числа Нуссельта N11 и массопереноеа через фазовую границу \'*2 при К л 111 в зависимости от значения безразмерного наклона фазовой кривой у при различной ширине перехода с). Сплошные кривые для с!=0.025, штриховые-для <1=0.125 и штрих-пунктирные для (1=0.25

На рис. 16 представлены интегральные характеристики эффекта фазового перехода по выносу тепла (число Нуссельта Ни) и массопереносу через фазовую границу у*2. При экзотермическом переходе (у>0) рост наклона кривой фазового равновесия увеличивает вынос тепла и массоперенос через фазовую границу. При этом увеличение размазанности фазового перехода несколько усиливает этот эффект. При эндотермическом фазовом переходе рост абсолютного значения наклона фазовой кривой |у| замедляет конвекцию. Область резкого замедления соответствует переходному режиму конвекции с эпизодическими аваланчами - прорывами вещества из слоя в слой. При больших абсолютных значениях наклона |у| устанавливается режим частично расслоенной конвекции. При этом, как видно из рис.16, число Нуссельта и массоперенос через фазовую границу уменьшаются. При очень больших значениях |у| массоперенос через границу стремится к нулю, а число Нуссельта стремится к конечному числу, соответствующему кондуктивному переносу тепла в зоне, где конвекция запирается. Влияние ширины фазового перехода для эндотермических фазовых переходов оказывается более заметным.

0Д5 7—0.35

Рисунок 17 То же, что и на рис. 14 при К но для числа Рэлея Ка=106. Левая колонка для ;--').!5, правая для т=-0.25. Безразмерные скорости течении показаны стрелками при максимальном значении У=1400 для -(=-0.15 и \ = 120(1 для -р-0.25.

На рис.17 и 18 представлены результаты расчета структуры мантийных течений с фазовыми переходами, подобные рис. 14-15 при том же значении плотностного отношения Р=Кру=2у, но при большей интенсивности конвекции, характеризуемой безразмерным числом Рэлея Яа=10б. Как видно из сравнения рис.14 и 15, при большей интенсивности конвекции фазовый переход сильнее влияет на ее структуру. При увеличении числа Рэлея в 10 раз конвекция при тех же параметрах становится более расслоенной. Поэтому теплоперенос через фазовую границу уменьшается, что ведет к перегреву нижней мантии и охлаждению верхней мантии. Влияние же ширины фазового перехода имеет похожие тенденции. Так при очень широкой зоне фазового перехода с1=0.25 даже при большом значении наклона |у| =0.25 фазовый переход практически

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

а а

Рис. 19 Рассчитанная зависимость массопереноса через фазову» границу \'*г в зависимости от ширины перехода (1 при различных значениях наклона фазовой кривой [у| при интенсивности конвекции 1?а=Н)5, 106.

Влияние ширины фазового перехода на структуру конвекции более отчетливо можно проследить на рис.19, на котором зависимость У*((1) представлена в явном ввде. Кривые У*(<1) приведены для различных значений наклона фазовой кривой |у|. Как видно из рис.19, при стремлении ширины фазового перехода к нулю массоперенос через фазовую гранипу плавно стремится к некоторому конечному значению. Выводы к ui.IV

Впервые рассчитанные зависимости теплового потока и степени массобмена от ширины ширины зоны фазового перехода во всем диапазоне значений для мантии Земли.

Глава V.Модели конвекции при высоких числах Рэлея в вытянутой по горизонтали области

Поскольку расчеты эволюции конвекции требуют большого компьютерного времени, а моделей потребовалось рассчитать несколько сот, то результаты предыдущих глав относились к моделям конвекции в квадратной области. Более приближены к реальной Земле модели в области, вытянутой по горизонтали в 5 и более раз. Кроме того, число Рэлея, характеризующее интенсивность конвекции, дая мантии Земли составляет 10 -108. Для этих случаев был просчитан ряд эволюционных моделей конвекции. На рис.20 приведены результаты расчета конвекции при Ка=107 в области с аспектным отношением х:г=5:1 для некоторых характерных моментов времени. Наклон фазовой кривой был взят у=-1.5 МПа/К и -3 МПа/К. Как видно из рис.20 при наклоне, по модулю меньшем критического -3.5, и при высоком числе Рэлея в вытянутой области нет ни глобального расслоения, ни аваланчей, но возникает небольшое торможение конвекции на границе 660 км, которое проявляется в структуре мантийных течений, и зависимости среднего модуля вертикальной скорости и распределения температуры по глубине.

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

йа=1 е7, д=-1.5 МПа/К

0 1 2 3 4 5

Рис. 20 Рассчитанные распределения температуры (серые гона) и конвективных скоростей (стрелки с максимальным значением 3 см/год) для постоянной вязкости, числе Рэлея 11а=107 при значениях фазового параметра у=-1.5 МПа/К (сверху) и у=-3 МПа/К (внизу). Также приведены графики зависимости от глубины для средней нададиаба тической температуры Т(2) и модуля вертикальной скорости |У,(г)(.

0 1 2 3 4 5

Рис. 21 Рассчитанные распределения температуры и конвекгнвных скоростей для переменной вязкости со скачком в 100 раз на ранице 660 км.

Выводы к гл-V

Показано, что в вытянутой по горизонтали области расслоение несколько затрудняется, а при увеличении интенсивности конвекции облегчается в соответствии ранее рассчитанной в гл.Ш диаграммой (рис.13). Скачок вязкости в 100 раз вызывает небольшое (несколько %) торможение массобмена через границу.

Глава VI. Фазовые переходы в мантии Земли и их влияние на конвекцию

Как указывалось в главе 1, структура мантийных течений зависит в основном не порознь от у и 8р, а от их произведения. Поэтому, наряду с безразмерным параметром P=Rpy, впервые введен несколько отличный размерный фазовый параметр Рр=(5р/ро)-ур, который выражается непосредственно через данные измерений относительного скачка плотности 8р/ро и наклона кривой равновесия фазового перехода yp=dp/dT.

Все результаты расчетов в предыдущих главах были приведены в безразмерных единицах. Для получения размерных величин их нужно умножить на соответствующие единицы измерения. Для мантии Земли D=2.9-106m, g=l0 м/с2, ро=4000 кг/м3, Т0 =2500 К, Оо=2-10"5 1/К, Ло=5'Ю21 Пах. При этих значениях число Рэлея равно Ra=1.5-107, единица измерения скорости Vo=k/D=1.1-10"' см/год единица измерения времени t0=D2/K= 2.6- 10й лет. При этом размерный ур и безразмерный у парамегры наклона кривой равновесия и фазовые параметры Рр и Р связаны соотношениями ур = (p0gD/T0)-y ~ (48 МПа/К)-у, Рр =p0gDa -Р =(2.4 МПа/К) -Р.

Как указывалось во введении, в мантии известно 5 основных фазовых переходов, на глубинах 410, 550, 660, 1500 и 2700 км. Значения параметров этих переходов (и наклоны кривых равновесия ур и скачки плотности 8р) сильно различаются между собой . Однако их размерные фазовые параметры, от которых зависит их влияние на конвекцию, оказываются почти равными для всех экзотермических переходов. Для фазового перехода в мантии на глубине 410 км размерный фазовый параметр равен Рр=0.07-2.3 МПа/К=0.16 МПа/К, на глубине 550 км он равен Рр=0.03-4.3 МПа/К=0.13 МПа/К, на глубине 660 км

(при максимально допустимом значении наклона ) он равен Рр=0.09-(-2.0) МПа/К=-0.18 МПа/К, на глубине 1500 км он равен Рр =0.01-16 МПаЛС=0.16 МПа/К, наконец, на глубине 2700 км он равен Рр=0.02-8 МПа/К=0.16 МПа/К.

На рис. 18 была приведена рассчитанная диаграмма зависимости теплового потока от безразмерного фазового параметра Nu(P). Для размерного параметра Рр она приведена на рис.22. Как видно из рис. 22, значения фазового параметра для всех экзотермических переходов попадают на начало кривой Nu(Pp). Они примерно на 2% увеличивают конвективный вынос тепла из мантии, ускоряя конвекцию. Значение параметра Рр для эндотермического перехода на глубине 660 км (даже при максимально возможном значении) также попадает на начало кривой, но при Рр<0. Причем со значительным запасом попадает в область однослойной конвекции, не достигая области переходного режима конвекции.

Рр= ТрАр/р

Рис. 22 Рассчитанные зависимости числа Нуссельта № от модуля размерного фазового параметра Рр=(6р/р0)7р при различной ширине перехода А

Как видно из рис.22 переходный и двухслойный режимы конвекции могли бы возникнуть только при достаточно больших значениях фазового параметра Рр. Например, при значении наклона кривой равновесия ур=-4 МПа/К и при скачке плотности (5р/р0)=Ю% размерный фазовый параметр будет равен Рр=-0.4 МПа/К. На диаграмме рис.22 это значение Рр для узкого фазового перехода

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

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

Результате; приведенные в диссертации, были получены для двумерных декартовых моделей. Поскольку эффект фазового перехода для мантии оказывается малым и проявляется только в некоторых локальных областях, то очевидно, что учет сферичности мантии не может существенно изменить полученные результаты Эти эффекты должны проявиться только при высоких значениях параметров, когда действительно возникает не локальное, а глобальное расслоение конвекции. При этом в трехмерных моделях расслоение еще более отодвигается в область больших значений модуля фазового параметра, т.к. слои разделяются не линией, а плоскостью, и имеется больше возможности проникнуть жидкости из одного слоя в другой. Этот аргумент согласуется с результатом высокоточного расчета трехмерной сферической модели Бунге, 1997 (потребовавшего недельной работы на одном из лучших суперкомпьютеров США). Даже дам значения наклона фазовой кривой ур=-4 МПа/К (в несколько раз превышенного по модулю измеренное значение) отсутствуют какие-либо проявления не только расслоения конвекции, но и переходного режима конвекции.

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

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

Заключение

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

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

позволяющий определять .. точность численного решения путем непосредственной подстановки его в уравнение Стокса

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

4. В применении к мантии Земли показано, в какой степени влияют на конвекцию фазовые переходы на глубине 410 км (оливин-вадслеит), 550км (вадслеит-рингвудит), 660 км (рингвудит-перовскит), 1500км (спиновый переход в железе), 2700 км (перовскит-постперовскит). Получен вывод о том, что эндотермический фазовый переход на глубине 660 км тормозит конвекцию, но несильно. Он не может приводить ни к расслоению конвекции, ни к глобальным аваланчам. На этой границе возможна лишь небольшая задержка и деформация нисходящих мантийных потоков (опускающейся океанической литосферной плиты) в некоторых регионах в соответствии с данными сейсмических измерений. При этом полный тепловой поток уменьшается на 4-5%. Каждый из всех остальных фазовых переходов немного ускоряет конвекцию, увеличивает перемешивание вещества мантии и увеличивает вынос тепла из мантии на ~ 2%.

Список опубликованных работ до теме диссертации

1) Трубицын В.П., Баранов А.А., Евсеев А.Н., Трубицын А.П.. Точные аналитические решения уравнения Стокса для тестирования уравнений мантийной конвекции с переменной вязкостью //Физика Земли.2006. № 7. С. 311.

2) Трубицын В.П., Баранов А.А., Евсеев А.Н., Трубицын А.П., Харыбин Е.В. Влияние низковязкой астеносферы на мантийные течения //Физика Земли.2006. №.12. С. 3-13.

3) Трубицын В.П., Евсеев А.Н., Баранов А.А., Трубицын А.П. Мантийная конвекция с эндотермическим фазовым переходом//Физика Земли.2007, №.12. С. 3-13.

4) Трубицын В.П., Евсеев А.Н., Баранов А.А., Трубицын А.П. Влияние эндотермического фазового перехода на массообмен между верхней и нижней мантией//ФизикаЗемли.2008, №6, с. 3-16

5) Трубицын В.П., Евсеев А.Н., Баранов А.А., Трубицын А.П. Структура конвекции при различной ширине зон фазовых переходов // Физика Земли.2008, №8, С.3-14

Тезисы докладов

1) V. Trubitsyn, A. Baranov,A. Evseev, A. Trubitsyn. Influence of endothermic phase transition on mantle convection. //EGU General Assembly,Vienna, 13-18 April 2008, Geophysical Research Abstracts, Vol. 10, EGU2008-A-04601,2008

2) Evseev A.N., Trubitsyn V.P., Baranov A.A., Trubitsyn A.P. (2008) - Numerical models of subduction of the oceanic and continental crust // EGU General Assembly, Vienna, 13-18 April 2008, Geophysical Research Abstracts, Vol. 10, EGU2008-A-08001, 2008

3) Евсеев A.H. "Мантийная конвекция и фазовые переходы" // Девятое международное совещание "Физико-химические и петрофизические исследования в науках о Земле", г. Москва (ИФЗ РАН, ГЕОХИ Р АН), п. Борок, Ярославская область (ГО "Борок" ИФЗ РАН). 7-11 октября 2008 г.

Содержание диссертации, кандидата физико-математических наук, Евсеев, Александр Николаевич

Введение

Глава I. Уравнения мантийной конвекции с переменной вязкостью и фазовыми переходами.

1.01 Уравнения мантийной конвекции с переменной вязкостью и фазовыми переходами.

1.02 Точное аналитическое решение уравнений мантийной конвекции

1.03 Влияние низковязкой астеносферы на структуру конвекции.

Глава II. Разработка программы численного моделирования двумерной мантийной конвекции методом конечных элементов.

Глава III. Влияние эндотермического фазового перехода на структуру мантийной конвекции и массобмен между нижней и верхней мантией.

3.01 Свойства перемежающейся конвекции: зависимость амплитуды и периода колебаний от числа Рэлея и наклона фазовой кривой.

3.02 Диаграмма критических значений наклона фазовой кривой, числа Рэлея и степени массообмена.

Г лава IV. Тепловая конвекция при различной ширине фазовых переходов

Г лава V. Модели конвекции при высоких числах Рэлея в вытянутой по горизонтали области.

Глава VI. Фазовые переходы в мантии Земли и их влияние на конвекцию.

Введение Диссертация по наукам о земле, на тему "Численные модели мантийной конвекции с переменной вязкостью и фазовыми переходами"

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

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

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

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

Данная гипотеза держалась в течение двух десятилетий. Анализу этой гипотезы и ее геофизических, геохимических и геологических следствии были посвящены сотни работ. Однако новые данные лабораторных измерений параметров фазового перехода (2003-2005 гг.) поставили под сомнение гипотезу о возможности расслоения конвекции в мантии Земли. В диссертации развита методика количественной оценки влияния фазовых переходов на конвекцию. Были рассчитаны модели конвекции для всего возможного диапазона их значений и найдены области однослойного, переходного и двухслойного режимов конвекции. Данные новых лабораторных измерений вещества мантии на рассчитанной диаграмме попадают в область общемантийной конвекции без расслоения и без аваланчей.

В пиролитовой (ПИРоксен-ОЛивин) модели Рингвуда верхняя мантия состоит из перидотита (90%) и эклогита (10%). Перидотит является ансамблем минералов: 60% оливина, 25% ортопироксена, а также клинопироксена и граната. Эклогит содержит клинопироксен (60%), гранат (30%) и другие малые добавки. Нижняя мантия состоит на 75% из силикатного перовскита (Mg0.9Fe0.i)SiO3, на 15% магнеовюстита (Mg0.8Fe0.2)O и 10% кальциевого перовскита (CaSi)C>3.

С ростом давления и температуры вещество мантии испытывает фазовые превращения (см. рис.1). Поскольку вещество мантии многокомпонентное, то даже в оливине переход происходит не скачком, а размазан на некоторый интервал давлений Ар или глубин d=Ap/(pg). В мантии на глубине 410 км оливин превращается в вадслеит со скачком плотности 5р/р0-0.07 при наклоне кривой фазового равновесия yp=dp/dT=1.6-3.0-2.3 МПа/К и ширине фазового перехода d=13 км. На глубине 520-580=550 км вадслеит превращается в рингвудит со скачком 5р/ро=0.03 при ур=4.3 МПа/К и ширине фазового перехода до d=60 км. На глубине 660 км рингвудит превращается в смесь перовскита и магнеовюстита с большим скачком плотности 8р/ро=0.09 при отрицательном наклоне кривой фазового равновесия (см рис.1). Его величина по последним данным лабораторных измерений оказывается существенно меньшей (Катсура и др. 2003, 2005; Литасов и др. 2005), чем обычно принималось в большинстве работ и составляет лишь yp=dp/dT=-(0.5-2.0)~-1.3 МПа/К при ширине d—З км. В интервале глубин 1250-1750 км железо в перовските переходит из состояния с высоким спином электрона в состояние с низким спином, при этом плотность меняется на Sp/po^O.Ol. В этом переходе также происходит перераспределение железа и магния между перовскитом и окисью магния. Наклон кривой равновесия оказывается достаточно высоким yp=dp/dT~ 16 МПа/К, и переход имеет ширину до d=500 км. Наконец, на глубине около 2700 км (на верхней границе слоя D") перовскит переходит в пост-перовскит при 5р/р0=0.02, Ур~8 МПа/К и ширине перехода d^=40 км. В Земле еще имеется фазовый переход базальта в эклогит на глубинах 60-100 км. Однако в отличие от рассматриваемых переходов он происходит не во всем веществе мантии, а только в океанических базальтах опускающейся литосферной плиты, т.е. только в одной компоненте. В диссертации приводятся результаты только для однокомпонентной тепловой конвекции с фазовыми переходами. Поэтому переход базальт-эклогит не рассматривается.

Рис. 1 Фазовая диаграмм для оливина, пунктирная кривая — распределение средней температуры в мантии.

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

Фазовый переход, в отличие от химического перехода, несмотря на скачок плотности, может или не влиять на конвекцию, или тормозить ее, или ускорять. Области устойчивости состояния фаз (см. рис.1) разделяются кривой фазового равновесия р=р(Т). Поскольку в мантии давление зависит от глубины p=pgh, то слои с разными фазами разделяются границей h=h(T). Если наклон фазовой кривой yp=dp/dT в переменных давление-температура или yh=dh/dT=Yf/pg в переменных глубина-температура равен нулю, то граница перехода не зависит от температуры и находится на одной и той же глубине в восходящем и нисходящем мантийном потоке. Пересекая фазовую границу, вещество мгновенно (если не учитывать малое время кинетики процесса) переходит в новую фазу. Поэтому вес столбов нисходящего и восходящего мантийных потоков не изменяется при таком фазовом переходе. Заметим, что в этом случае фазовый переход происходит без выделения и поглощения тепла, т.к. наклон фазовой кривой связан с теплотой q, выделяющейся при фазовом переходе, соотношением Клапейрона-Клаузиуса q=Yp5pT/(p,p2), где pi и р2-плотности фаз, и 6p=p2~pi и yp=dp/dT.

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

Если наклон фазовой кривой уь отличен от нуля, то граница фазового перехода в горячем восходящем и холодном нисходящем мантийных потоках будет смещаться, при этом различно. На рис.2 схематически изображена конвективная ячейка и смещения фазовой границы на глубине 660 км при отрицательном значении yh. В области горячего восходящего мантийного потока фазовая граница приподнята. Поэтому восходящий мантийный поток утяжеляется. Аналогично облегчается холодный нисходящий мантийный поток. В результате конвекция замедляется. Аналогичные рассуждения приводят к тому, что при положительном наклоне фазовой кривой конвекция должна ускоряться.

Влияние фазового перехода на мантийную конвекцию изучается уже много лет, начиная с работ Христенсена 1984-1986 гг. Результаты двумерного численного моделирования показали, что при отрицательном наклоне фазовой кривой при значениях модуля |ур|>2.5-3.0 МПа/К торможение вертикальных мантийных потоков приводит к переходному режиму конвекции (при котором конвекция попеременно расслоена во времени и пространстве), а при более высоких значениях |ур| затем происходит глобальное расслоение конвекции. До 2002 г. лабораторные измерения оливина давали значение наклона ур в диапазоне от -2 МПа/К до -2.5 МПа/К. Эти значения по модулю несколько меньше критического, но находятся на грани точности измерений и расчетов. Поэтому гипотеза о расслоении конвекции интенсивно обсуждалась. После 2002 г. была существенно улучшена методика измерений, т.к. фазовый переход стали фиксировать непосредственно во время опыта, а не по остаточным образцам после снятия давления и температуры. Эти измерения дали меньшее по модулю значение наклона ур от -0 до -2 МПа/К. Поэтому гипотеза о возможности глобального расслоения конвекции в мантии Земли стала пересматриваться. Но поскольку фазовый переход с отрицательным наклоном все-таки тормозит конвекцию и появились возможности более точного моделирования конвекции, то стала актуальной проблема более детального расчета влияния фазовых переходов, причем как эндотермических, так и экзотермических.

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

Цели и задачи работы

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

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

Актуальность работы

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

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

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

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

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

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

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

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

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

Теоретическая и практическая значимость

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

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

Личный вклад автора

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

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

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

B.П. Трубицына была использована форма записи уравнений конвекции через новый параметр, что позволило количественно представить эффекты всех основных известных в мантии Земли фазовых переходов в единой форме.

Список опубликованныхработ по теме диссертации

1) Трубицын В.П., Баранов А.А., Евсеев А.Н., Трубицын А.П. Точные аналитические решения уравнения Стокса для тестирования уравнений мантийной конвекции с переменной вязкостью //Физика Земли, 2006, № 7,

C. 3-11.

2) Трубицын В.П., Баранов А.А., Евсеев А.Н., Трубицын А.П., Харыбин Е.В. Влияние низковязкой астеносферы на мантийные течения //Физика Земли.2006, №.12, С. 3-13.

3) Трубицын В.П., Евсеев А.Н., Баранов А.А., Трубицын А.П. Мантийная конвекция с эндотермическим фазовым переходом //Физика Земли.2007, №.12, С. 3-13.

4) Трубицын В.П., Евсеев А.Н., Баранов А.А., Трубицын А.П. Влияние эндотермического фазового перехода на массообмен между верхней и нижней мантией // Физика Земли.2008, №6, с. 3-16

5) Трубицын В.П., Евсеев А.Н., Баранов А.А., Трубицын А.П. Структура конвекции при различной ширине зон фазовых переходов // Физика Земли.2008, №8, С.3-14

Тезисы докладов

1) V. Trubitsyn, А. ВагапоуД. Evseev, A. Trubitsyn. Influence of endothermic phase transition on mantle convection. //EGU General Assembly,Vienna, 13-18 April 2008, Geophysical Research Abstracts, Vol. 10, EGU2008-A-04601, 2008

2) Evseev A.N., Trubitsyn V .P., Baranov A.A., Trubitsyn A.P. (2008) - Numerical models of subduction of the oceanic and continental crust // EGU General Assembly, Vienna, 13-18 April 2008, Geophysical Research Abstracts, Vol. 10, EGU2008-A-08001, 2008

3) Евсеев A.H. "Мантийная конвекция и фазовые переходы" // Девятое международное совещание "Физико-химические и петрофизические исследования в науках о Земле", г. Москва (ИФЗ Р АН, ГЕОХИРАН), п. Борок, Ярославская область (ГО "Борок" ИФЗ РАН). 7-11 октября 2008 г.

Структура и объем диссертации

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

Заключение Диссертация по теме "Геофизика, геофизические методы поисков полезных ископаемых", Евсеев, Александр Николаевич

Заключение

В диссертационной работе были получены следующие наиболее важные результаты.

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

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

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

4. В применении к мантии Земли показано, в какой степени влияют на конвекцию фазовые переходы на глубине 410 км (оливин-вадслеит), 550км (вадслеит-рингвудит), 660 км (рингвудпт-перовскит), 1500км (спиновый переходе железе), 2700 км (перовскит-постперовскит). Получен вывод о том, что эндотермический фазовый переход на глубине 660 км тормозит конвекцию, но несильно. Он не может приводить ни к расслоению конвекции, ни к глобальным аваланчам. На этой границе возможна лишь небольшая задержка и деформация нисходящих мантийных потоков (опускающейся океанической литосферной плиты) в некоторых регионах в соответствии с данными сейсмических измерений. При этом полный тепловой поток уменьшается на 4-5%. Каждый из всех остальных фазовых переходов немного ускоряет конвекцию, увеличивает перемешивание вещества мантии и увеличивает вынос тепла из мантии на ~ 2%.

Библиография Диссертация по наукам о земле, кандидата физико-математических наук, Евсеев, Александр Николаевич, Москва

1., Система Mg0-Fe0-Si02 при высоких давлениях и температурах-фазовые равновесия и упругие свойства, в сб. Верхняя мантия, ред. А. Ритсема, 1975, «Мир» Москва* с. 60-80.

2. Браун Д., Массет А., Недоступная Земля, Москва: Мир, 1984.

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

4. Жарков В.Н. Внутреннее строение Земли и планет. М.Наука,1978

5. Каракин А.В. Магницкий В.А., Калашникова И.В. О горизонтальных и вертикальных перемещениях литосферы // Изв. АНСССР, Физика Земли, 1974, № 9.

6. Каракин А.В. К вопросу о природе длинноволновых вертикальных и горизонтальных смещений земной поверхности // Геология и геофизика, 1976, № 6, с. 75-81.

7. Каракин А.В., Магницкий В. А., Калашникова И.В. Об эффектах "смазочного слоя", возникающих при горизонтальных перемещениях литосферы // Докл. АН СССР, 1974, т. 214, № 3, с. 561-564.

8. Котелкин В.Д., Лобковский Л.И. Общая теория Мясникова эволюции планет и современная термохимическая модель эволюции Земли// Физика Земли. 2007. № 1, с. 26-44.

9. Ландау JI.Д. и Е.М. Лифшиц, Статистическая физика, 1964, Наука Москва. С. 567

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

11. Марчук Г.И., Методы вычислительной математики, М.Наука, 1989

12. Самарский А.А. Гулин А.В. Численные методы. М.Наука, 1986, 430с.

13. Сорохтин О.Г.,УшаковС.А., Глобальная эволюция Земли. М.:МГУ,1991.

14. Трубицын В.П. , Баранов А.А., Евсеев А.Н., Трубицын А.П. Точные аналитические решения уравнения Стокса для тестирования уравнений мантийной конвекции с переменной вязкостью//Физика Земли. 20066. № 7. С. 3-11.

15. Трубицын В.П. ,Симакин А.Г., Баранов А.А. Влияние пространственных вариаций вязкости на структуру мантийных течений //Физика Земли. 2006а. № 1. С. 3-15.

16. Трубицын В.П. Геодинамическая модель эволюции Тихого океана //Физика Земли. 2006. №2. С. 3-25.

17. Трубицын В.П. Тектоника плавающих континентов//Вестник РАН 2005, № 1. с. 10-21.

18. Трубицын В.П., Баранов А.А., Евсеев А.Н., Трубицын А.П., Харыбин Е.В. Влияние низковязкой астеносферы на мантийные течения. //Физика Земли. 2006. № 12. С. 11-19.

19. Трубицын В.П., Баранов А.А., Харыбин Е.В. Численные модели субдукции океанической коры с базальтовыми плато. //Физика Земли. 2007. №7. С. 3-10.

20. Трубицын В.П., Евсеев А.Н., Баранов А.А., Трубицын А.П. Мантийная конвекция с эндотермическим фазовым переходом //Физика Земли. 2007. № 12. С. 3-10.

21. Трубицын В.П., Рыков В.В., численные модели эволюции мантийной конвекции//сб. Глобальные изменения природной среды, ред. H.JI. Добрецов. Наука, Новосибирск, 2002, 42-56.

22. Трубицын В.П., Фазовые переходы, изотермическая сжимаемость и тепловое расширение Земли, Физика Земли, 1979, № 1, с. 21-27.

23. Трубицын В.П., Фазовые переходы, сжимаемость, тепловое расширение Земли, теплоемкость и адиабатическая температура в мантии, Физика Земли, 2000, № 2, с. 3-16.

24. Трубицын В.П., Харыбин Е.В. Конвекция в вязкой жидкости с оседающими частицами // Физика Земли, 2005, №12. с. 3-11

25. Anderson D. L., Theory of the Earth, Blackwell Scientific Publ.,Boston, Oxford, London, Edinburgh, Melourne, 1989, P. 366.

26. Baumgardner J. R., "Application of supercomputers to 3-D mantle convection," in The Physics of the Planets, S. K. Runcorn, ed., John Wiley and Sons, pp. 199-231, 1988.

27. Bercovici D., 2003. The generation of plate tectonics from mantle convection, Earth Planet. Sci. Lett., 205, 107-121.

28. Blankenbach В., Busse F., U. Christensen, L. Czerepes, D. Gunkel, U. Hansen, G. Jarvis, M. Koch, G. Marquart. D. Moore, P. Olson, H. Schmeling, and T. Schnaubelt, 1989. A benchmark comparison for convection codes// Geophys. J. Int. 98, P. 23-38.

29. Boyd, F.R., 1989. Compositional distinction between oceanic and cratonic lithosphere. Earth planet. Sci. Lett., 96, 15-26.

30. Brunet D. and Ph. Machtel, Large-scale tectonic features induced by mantle avalanches with phase, temperature, and pressure lateral variations of viscosity. J. Geph.Res., 1998, V. 103, pp. 4920-4945.

31. Brunet D., and P. Machete1, 1998. Large-scale tectonic features induced by mantle avalanches with phase, temperature and pressure lateralvariations of viscosity, J. Geophys. Res., 103, 4929-4945.

32. Bunge H.P., and Baumgardner J.R., Mantle Convection Modeling on Parallel Virtual Computers., Computers in Physics, v.9, pp. 207-215, April 1995.

33. Bunge H.P., Richards M.A., and Baumgardner J.R.,A Sensitivity Study of 3-D Spherical Mantle Convection at 10A8 Rayleigh Number, Journal of Geophysical Research, in press, 1997.

34. Bunge H.P., Richards M.A., and Baumgardner J.R.,The Effect of Depth- Dependent Viscosity On The Planform Of Mantle Convection, Nature, v. 379, pp. 436-438, February 1, 1996.

35. Bunge H.P., Richards M.A.,"The Origin of Large Scale Structure in Mantle Convection," Geophysical Research Letters, v. 23, pp. 2987-2990, October 15, 1996.

36. Christensen U., and D. A.Yuen, Layered convection induced by phase transition // J. Geophys. Res., 1985, vol. 90, pp. 10291-10300.

37. Christensen U., Effects of phase transitions on mantle convection// Annu. Rev. Earth Planet. Sci., 1995, vol. 23, pp. 65-87.

38. Christensen, U.R., and D. A. Yuen, 1984. The interaction of a subducting lithospheric slab with a chemical or phase boundary, J. Geophys. Res. 89, 4389-4402.

39. Christensen, U.R., and D. A. Yuen, 1985. Layered convection induced by phase transition, J. Geophys. Res. 90, 10291-10300.

40. Christensen, U.R., Heat transport by variable viscosity convection for the Earth's thermal evolution// Phys. Earth and Planet. Interiors. 1984. V. 35. P. 264-282.

41. Davies G.F. Dynamic Earth, Cambridge University Press, Cambridge, UK, 1999.

42. Davies G.F., 1988.Role of the lithosphere in mantle J. Geophys. Res., 93, 10451-10466.

43. Davies G.F., 1999. Dynamic Earth, Plates, Plumes and Mantle Convection, Cambridge University Press, p. 458

44. Fortin M.and Fortin A.A. Generalization of Uzawa's algorithm for the solution of the Navier-Stokes equations//Comm. Appl. Numer. Methods. 1985.V.l.P .205-210.

45. Gurnis M. and B.H. Hager, 1988. Control of the structure of subducted slab, Nature 335,317-321.

46. Gurnis, M, J. Ritsema, H. van Heijst, and S. Zhong, Tonga slab deformation: The influence of a lower mantle upwelling on a slab in a young subduction zone, Geophys. Res. Lett., Vol 27, p 2373, 2000.

47. Gurnis, M. and S. Zhong, Generation of long-wave length wavelength heterogeneity in the mantle by the dynamic interaction between plates andconvection, Geophys. Res. Lett., 18, 581-584, 1991

48. Hansen, G. Jarvis, M. Koch, G. Marquart. D. Moore, P. Olson,H. Schmeling, and T. Schnaubelt, 1989. A benchmark comparison for convection codes, Geophys. J. Int., 98, p, 23-38.

49. Hawkesworth С .J., Kempton P .D., Rogers N.W., Ellam R.M. and van Calsteren P.W., 1990. Continental mantle lithosphere, and shallow level enrichment processes in the Earth's mantle. Earth plan. Sci. Lett., 96, 256268.

50. Huang, J., S. Zhong, and J. van Hunen, Controls on sub-lithospheric small-scale convection, JGR, 108, 2405, doi: 10.1029/2003JB002456,2003.

51. Hughes T.J.R., Liu W.K.Brooks A.Finite Element Analysis of Incompressible Viscous Flows by the Penalty Function Formulation// J.Comput.Phys. 1979. V.30.P.1-60.

52. Ita J, Stixtrude L., Petrology, elasticity and composition of the mantle transition zone//J. Geophys. Res. 1992, V. 97, pp. 6849-6866.

53. Khodakovskii G., Rabinowicz M., Ceuleneer G., Trubitsyn V .P., Melt percolation in a partially molten mush: effect of a variable viscosity //Earth and Planetary Sci. Lett.1995,, V. 134, P. 267-282.

54. Korenaga J. Firm mantle plumes and the nature of core-mantle boundary region //Earth Planet. Sci. Lett. 2005. V. 232. P. 29-37.

55. Lenardic A. and W.M. Kaula, 1995. Mantle dynamics and the heat flow into the Earth's continents, Nature, 378, 709-711.

56. Marsh В., "Magma Chambers," Ann. Rev. Earth Planet. Sci. 17, 439-474(1989).

57. McNamara A. K., and Zhong S., Degree-1 mantle convection: Dependence on internal heating and temperature-dependent rheology, GRL, 32, L01301, 10.1029/2004GL021082, 2005

58. McNamara A. K., and Zhong S., The influence of thermochemical convection on the fixity of mantle plumes, EPSL, vol 222, 485-500, 2005.

59. Moresi L., Zhong S., and Gurnis M., The accuracy of finite element solutions of Stokes' flow with strongly varying viscosity //Phys. Earth Planet. Inter., 97, 83-94, 1996

60. Moresi L.N. and Solomatov S.S, Numerical investigation of 2-D convection with extremely large viscosity variations //Phys. Fliid. 1995. V. 7, №9, P. 2154-2162.

61. Moresi, L., M. Gurnis, and S. Zhong, Plate tectonics and convection in the Earth's mantle: Toward a numerical simulation, Computing in Sci. Eng., 2, no. 3, 22-33, 2000.

62. Parsons D., Daly S. The relationship between surface topography, gravity anomalies and temperature structures of convection//J. Geophys. Res. 1983. V. 88. P. 1129-1144.

63. Pelletier D.Fortin A.Camarero R. Are FEM solutions of incompressible flows really incompressible? /ЯЛ. forNumer. Methods in Fluids. 1989.V.9.P.99-112.

64. Peltier W.R., Solheim L.P. Mantle phase transitions and layered chaotic convection // Geophysical research letters, Vol.19, No.3, P.321-324, 1992

65. Poirier W.R. Introduction to the Physics of the Earth's interior. Cambridge Univ. Press, NewYork, 1991, P. 410.

66. Ringwood A.E., Composition and Petrology of the Earth's mantle, McGraw-Hill, New York, 1975, P. 618.

67. Ritzwoller M.H., N.M. Shapiro, and S. Zhong, Cooling history of the Pacific lithosphere, EPSL, 226, 69-84, 2004.

68. Rudman M., "Two-Phase Natural Convection for Crystal Settling in Magma Chambers," Phys. Earth Planet. Inter. 72, 153-172 (1992).

69. Schott В., Yuen D.A., and H. Schmeling. 2000. The diversity of tectonics from fluid-dynamical modeling of the lithosphere-mantle system, Tectonophysics, 322(1-2), 35-51.

70. Schubert G., D.L. Turcotte and P. Olson. Mantle Convection in the Earth and Planets. Cambridge Univ. Press. 2001. p. 940.

71. Simakin A., Schmeling H., TrubitsynV. Convection in meltsduetosedimentative crystal flux from above // Phys. Eart Planet. Inter., 1997.1. V.102, P. 185-200.

72. Solheim L.P., Peltier W .R. Avalanche effects in phase transition modulated thermal convection: A model of Earth's mantle // Journal of geophysical research, vol.99, no.B4, April 10, 1994

73. Solheim L.P. and W.R. Peltier, Phase boundary deflections at 660-km depth and episodically layered isochemical convection in the mantle, J. Geophys. Res., 1994, vol. 99, pp. 15861-15875.

74. Solomatov V.S. Grain size-dependent viscosity convection and the thermal evolution of the Earth//Earth Planet. Sci. Lett. 2002. V. 191. P. 203212.

75. Solomatov V.S., El-Khozondar R., Tikare V. Grain size in the lower mantle: constraints from numerical modelling of grain growth in two-phase system //Phys. Earth. Planet. Inter. 2002. V. 129. pp. 265-282.

76. Stixrude L., Lithgow-Bertelloni C. Influence of phase transformations on lateral heterogeneity and dynamics in Earth's mantle // Earth and Planetary Science Letters 263 (2007), pp. 45-55

77. Taira A., Mann P., Rahardiavan R. Incipent subduction of the Ontong Java plateau along the North Colomon trench//Tectonophysics. 2004. V. 389. P. 247-266/

78. Tan E., P. Thoutireddy, E. Choi, M. Gurnis, and M. Aivazis. 2000. GeoFramework Parti: Coupling models of mantle convection with Pythonframework. Geochemistry, Geophysics, Geosystems.

79. Trubitsyn V. P. and E. V. Kharybin, "Convection in Magma Chambers Due to an Inversion of the Depth Distribution of Settling Crystals," Fiz. Zemli, No. 5, 47-52 (1997).

80. Trubitsyn V. P. and E. V. Kharybin, "Convective Instability of the Sedimentation Regime in the Mantle," Izv. Akad. Nauk SSSR, Fiz. Zemli, No. 7, 21-30 (1987).

81. Trubitsyn V. P. and E. V. Kharybin, "Thermosedimentary Convective Instability of a Two-Component Viscous Liquid," Fiz. Zemli, No. 2, 3-17 (1991).

82. Trubitsyn V. P., A. M. Bobrov, and E. V. Kharybin, "Convection in Magmas Due to a Horizontal Gradient of Density," Izv. Akad. Nauk SSSR, Fiz. Zemli, No. 6, 3-15 (1989a).

83. Trubitsyn V. P., A. M. Bobrov, and V. V. Kubyshkin "Thermal Convection in the Mantle Due to Horizontal and Vertical Gradients of Temperature," Fiz. Zemli, No. 5, 12-23 (1991).

84. Trubitsyn V .P., Kaban M., Mooney W., Reigber Ch., Schwintzer P. Simulation of active tectonic processes for a convective mantle with moving continents//Geophys. J. Intern. 2006. V. 164. P. 611-623.

85. Trubitsyn, V.P., 2000. Phase Transitions, Compressibility, Thermal Expansion, and Adiabatic Temperature in the Mantle, Izvestiya, Physics ofthe Solid Earth, V. 36, P. 101-113.

86. Watts А. В., and S. Zhong, Observations of flexure and the rheology of the oceanic lithosphere, Geophys. J. Int., 142, 855-875, 2000.

87. Weidber D.J., Mantle models based on measured physical properties of minerals, in Chemistry and physics of terrestrial planets, edd. By S.K. Saxena, Springer-Verlag, New York, 1996 , p. 1-30.

88. Zhang S. and Christensen U. Some effects of lateral viscosity variations on geoid and surface velocities induced by density anomalies in the mantle//Geophys .J. Intern. 1993.V. 170. P. 531-547.

89. Zhong S. and M. Gurnis, 1995. Mantle convection with plates and mobile, faulted plate margins, Science, 267, 838-843.

90. Zhong, S. and M. Gurnis, 1994. Role of plates and temperature-dependent viscosity in phase change dynamics, J. Geophys. Res. 99, 15,903

91. Zhong, S. and M. Gurnis, Viscous flow model of a subduction zone with a faulted lithosphere: long and short wavelength topography, gravity and geoid, Geophys. Res. Lett., 19, 1891-1894, 1992

92. Zhong, S., and M. Gurnis, Dynamic feedback between an non-subducting raft and thermal convection, J. Geophys. Res., 98,12219-12232, 1993

93. Zhong, S., and M. Gurnis, Dynamic interaction between tectonic plates, subducting slabs, and the mantle, Earth Interactions, 1997.

94. Zhong, S., M. T. Zuber, L. Moresi, and M. Gurnis, Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection, J. Geophys. Res., 105, 11063-11082, 2000. 2000.

95. Zhong, S., M.T. Zuber, L. Moresi and M. Gurnis, 2000. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection, J. Geophys. Res. 105, 11,063 -11,082.

96. Zhong S. Analitic solutions for Stokes' flow with lateral variations in viscosity //Geophys. J. Intern. 1996,V. 124.№1, P. 18-28.