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

Баранов, Алексей Андреевич

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

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

российская академия наук

институт физики земли им. О.Ю. Шмидта

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

БАРАНОВ Алексей Андреевич

ВЛИЯНИЕ РАДИАЛЬНЫХ И ЛАТЕРАЛЬНЫХ

ВАРИАЦИЙ ВЯЗКОСТИ НА СТРУКТУРУ ТЕПЛОВОЙ КОНВЕКЦИИ В МАНТИИ ЗЕМЛИ

Специальность 25 00 10 Геофизика, геофизические методы поиска полезных ископаемых

Автореферат иисиь4577

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

Москва -2007

003064577

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

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

член-корреспондент РАН, д ф м.н. Валерий Петрович Трубицын

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

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

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

Ведущая организация:

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

Защита состоится 23 октября 2007 г в 12 ч на заседании диссертационного совета К 002 001 02 Института физики Земли Российской Академии Наук (ИФЗ РАН) по адресу: 123995, г Москва, ул Большая Грузинская, 10

С диссертацией можно ознакомиться в библиотеке ИФЗ РАН

Автореферат разослан «10» сентября 2007 г.

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

Э А.Боярский

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

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

Понимание геодинамики базируется на знаниях структуры мантии и ядра Земли и процессов, происходящих в них

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

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

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

Цели работы определили постановку следующих задач:

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

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

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

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

Основные результаты работы, выносимые на защиту:

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

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

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

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

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

Основным методом исследований являлся численный эксперимент Метод основан на использовании конечных элементов

Разработка программного средства велась на языке программирования Фортран Компилятором и средой разработки был выбран пакет Compaq Visual Fortran 6 6

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

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

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

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

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

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

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

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

Обсуждение работы

Результаты исследований докладывались и обсуждались на научных семинарах Института физики Земли РАН и Института Экспериментальной Минералогии РАН, были представлены на международных научных конференциях, в том числе на генеральной Ассамблее Европейского Союза Наук о Земле (Вена, Австрия — 2006,2007)

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

Диссертация состоит из введения, 4 глав, заключения, списка литературы и приложения Общий объем диссертации составляет 99 страниц машинописного текста, в том числе 20 рисунков и список литературы из 79 наименований

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

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

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

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

2 Во второй работе Автором на основе идеи Трубицына В П найдено точное аналитическое решение уравнения Стокса для двумерной модели конвекции с произвольными вариациями вязкости

3 В третьей работе Автором была усовершенствована известная программа расчета двумерной конвекции Китком и протестирована оригинальная программа расчета мантийной конвекции По совету Трубицына В П, исследованы низковязкие области в мантии и получены критерии возникновения автономной конвекции в астеносфере

4 В четвертой работе Автором добавлены активные маркеры в оригинальную программу расчета конвекции Это позволило рассчитывать не только тепловую но и термохимическую конвекцию По совету научного руководителя программа была применена для исследования субдукции океанической и континентальной коры

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

Основное содержание работы Введение.

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

Глава 1, Конвекция в мантии Земли.

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

Агрегатное состояние, в котором находится вещество Земли, зависит от температуры и давления Температура растет с глубиной сначала быстро с градиентом около 30 градусов на километр, затем более медленно с градиентом в среднем около 1 градуса на километр В центре Земли температура порядка 5 тысяч градусов, а давление вышележащих слоев достигает 3 5 млн бар На границе мантии с ядром температура около 4 тысяч градусов и давление порядка 1 5 млн бар Также в полутвердом состоянии вблизи температуры плавления находятся и тугоплавкие силикаты в большей части мантии Поэтому они способны течь подобно высоковязкой жидкости

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

Глава 2. Численное моделирование двумерной мантийной конвекции.

Постановка задачи. Точное аналитическое решение уравнения Стокса

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

температурой 0°С В результате решения взаимосвязанной системы уравнений переноса массы, энергии и импульса для мантии находим распределения температуры, скоростей течений и давления в мантии Все эти величины являются не только функциями координат, но и времени

Уравнения мантийной конвекции. Граничные условия.

Система уравнений тепловой конвекции включает в себя уравнение Навье-Стокса, уравнение переноса тепла и уравнение неразрывности рО\У<й= -Эр/йх.н- дх1}1дк}+ р(Т)§ 5,3, (1)

Ш7сИ= а(кдТ/дх, )Эх, + Н, (2)

др№ + Э(У, р)/3х, = 0, 1=1,2, 3, (3)

где р(Т) — плотность мантии, g, — ускорение силы тяжести, Т — температура, отсчитываемая от адиабатического распределения, к — коэффициент теплопроводности, Н — термометрическая плотность тепловых источников, 8 — символ Кронекера, равный 1 при ¡=] и равный 0 при Щ Девиаторный тензор вязких напряжений т ч

= (4)

где лСрЛО — вязкость. Уравнения содержат неизвестные функции вектор скорости У,(х,Д), температуру Т(х,Д), давление р(х,д) и тензор вязких напряжений тц.

Здесь Ка=ар0§Т0О7кг10— число Рэлея, определенное по значению вязкости % (5)

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

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

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

У2(х,г)= У08т(я г)соз(ях), (6)

Ух(х,2)= Ух(г) зт(ях) (7)

Неизвестную функцию распределения температуры представим в виде Т(х,г)=(1-г)+ ДТ(х,г)

Подставим эти формулы в уравнение несжимаемости и после дифференцирования получим лУх(г)+У0псо8(лг)=0 или Ух(х,г)= -\осо$(пг) 81п(ях)

Подставим формулы (6) и (7) в выражения для тензора напряжений а**—р+2т]ЭУх/ск~р-2г1яУ0со8(я;г)со8(ях), (8)

0x2= 0, (9)

-р+2т!дУ2 р+2г|лу0соз(7с2)со5(я;х) (10)

"Для простоты решение выписано для квадратной области Однако заменой переменных х'= (1/с1)х оно легко обобщается на прямоугольную область с произвольным отношением ширины 1 к толщине <1 Полученное решение также легко обобщается на более сложные течения, представляемые более высокими гармониками типа У2(х,2)=У081п(ш12)со8(шях), а также, в виду линейности уравнения Стокса, любыми их линейными комбинациями

Решение иногда удобнее записать, если амплитуду в распределении температуры обозначить через То, т е если сделать замену 4л(УоЛ1а)= Т0 или Уо= 11аТо/(4я) В результате частное решение уравнения Стокса для распределения скоростей течений, напряжений и температуры примет вид

Уг(х,г)= Ка(То/4я2)8т(и:г)со8(1сх) , (11)

Ух(х,г)= -Ка(То/4я2)со8(ж2;) зт(лх), (12)

ахх=Ка(г2/2-г+1/2+Т0Т101/2я), 0, (13)

Яа{(г2/2-г+1/2)+ (То/л) [г|соз(и:г)со8(лх)+Т1 м /2]}, (14)

р = -Ка{(г2/2-2+1/2) + (То/2я)[г|со8(яг)со8(ях)+г|о1]}, (15)

Т=1-г+ Т0со8(ях) [т|зш(яг)—(Эт(/яЭг)со8(712)] (16) При этом компоненты вязкого тензора напряжений и функция тока будут равны

тн—Ка(То/2я)г|со8(яг)со8(71х), (17)

0, (18)

тш= Ка(То/2п)г] соз(яг)соз(ях), (19)

^(х^)—Ка(То/4я2)81п(яг) вт(лх) (20) Полученное решение справедливо для произвольной вязкости, заданной

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

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

Результаты тестирования программы расчета мантийных течений

Для численного тестирования возьмем их равными Ка=10 и То= 0.5, так что КаТ0=1.

1 0.8

0.6

0.4

02 о

Рис. I Рассчитанные поля вязкое™, температуры Т{х,7.), скоростей течений У(х,г) п напряжений для модели 1 при 41=1 НООхг в безразмерных единицах. Скорости

течений показаны стрелками с максимальным значением безразмерной скорости ^'пш, =0.0253. На рис. 1Л температура изображена изолиниями от 0 до I, а вязкость серыми тонами. На рис. 1 й безразмерное напряжение показано серыми тонами и одновременно изолиниями.

100 wo тш :')<■ -153 .50 Sfl 15а -2iM -tto ü

Рис.2. Рассчитанные ноля вязкости. температуры Т(х,7_), скоростей течений У(>:,г), напряжений лДх,г) и давления р для модели 1 при г|1=1000-999/[1+схр(-100(х-0,5)) в безразмерных единицах. Скорости течений показаны стрелками с максимальным значением безразмерной скорости \Гп,я>.=;0.0253. Температура на рис. 2Л указана изолиниями от 0 до 1. вязкость изображена серыми тонами. На рис. 2Б и 2В напряжение т„ и давление р изображены серыми тонами и одновременно изолиниям^-

Рис-3. Го же. что и на рис. 2, но для распределения вязкости г|3 --1000-999/[ 1 +ехр(-100(7 -0.5})J.

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

002*3

1 1 " 1 - 1 4 1 1 I - 1

- 1 . 1 1 1 • ч - V,..

............... Г" 1...........Г".......1 ........

120 №

Рис 4 Зависимость рассчитанной максимальной скорости течений от шага расчетной сетки N4 Кривая А — для модели 2, Б — для модели 3

Результаты тестирования можно проиллюстрировать кривой зависимости рассчитанной максимальной скорости течений от сетки Согласно аналитическому решению для всех моделей вязкости эта скорость в безразмерных единицах равна Утах=0 025327 На рис 4 приведены значения полученные численным решением уравнения Стокса при различных расчетных сетках соответственно для модели вязкости 2 и 3 Как видно из рис 4, для обеих моделей со скачками вязкости в три порядка приемлемая точность достигается уже на сетке 1ЧХ N¿=50x50

Выводы

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

Создана, отлажена и протестирована программа, которая позволяет рассчитывать двумерные модели конвекции в мантии с большими скачками вязкости (до 6 порядков) На сетке 10x10, 50х50, 100x100 достигается точность совпадения вычисленной величины с аналитической соответственно 5, 2 и 0 5 процента

Глава 3. Влияние горизонтальных и латеральных вариаций вязкости на

структуру конвекции

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

На рис 5 представлены рассчитанные скорости конвективных течений в квадратной области с высоковязким включением типа плиты, занимающей область г>0 8 и 0 2<х<0 8, а также типа слоя г>0 8 Безразмерная вязкость в основной области равна г)=1, вязкость включения г(=100 В модели А плита фиксирована в пространстве (с условием прилипания к верхней границе области) В модели Б вся верхняя граница скользкая, и плита может перемещаться в пространстве вместе с конвективными течениями

Из рис 5 видно, что скорости течений в высоковязкой плите в моделях А и Б кардинально различны В первом случае закрепленной плиты течения слабо проникают в высоковязкую область и концентрируются в основной области, обтекая высоковязкое включение Во втором случае свободная плита под влиянием течений в основной области начинает двигаться со скоростью, примерно равной средней скорости течений, которые были в отсутствие плиты Однако благодаря повышенной вязкости скорости течений внутри плиты сравниваются между собой

Среднюю скорость течений, как обычно, будем характеризовать среднг " ,ю У=Ути

Повышенная вязкость.

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

6;:

} о.в

-о«,, . . ■ ■ -

Г! ¡1 . !;02 •:

.... 1

о- - Г • • г.Т:..... .£, - - .--г............

0-2 0.4 О-в 1 0.2 0.4 0 6 0 £ 1 6 2 А< 6 в В> !

д 6 в

Рис,5. Рассчитанные мгновенные скорости конвективных течений в неравномерно нагретой жидкости с высоковязкими включениями (выделенными темным тоном). Безразмерная вязкость жидкости 11=1, вязкость включений 11= 100. В модели А плита фиксирована а пространстве, в мидели В плита может свободно перемещаться, в модели В высоковязкое включение занимает нссь верхний слой расчетной области.

В литературе распространено мнение, что скорости течений в высоко вязкой области всегда примерно обратно пропорциональны вязкости. [Могек!, Зо1отаюу, 1995]

Поэтому, при решении уравнений конвекции рекомендуется вводить новую функцию, равную произведению скорости течений на вязкость и=Уг[, которая должна быгь более гладкой по сравнению со скоростью V. Однако, как следует из рис, 5, функция 1Дх,г)=Уг5 на границе с высоковязким подвижным включением не только не становится более гладкой по сравнению с обычной скоростью У(х,г), но, наоборот, имеет больший скачок. Модель с вы со ко вязки м слоем была просчитана для различных значений вязкости в диапазоне от 1 до Ю5. Скорости течений внутри слоя убывают обратно пропорционально его вязкости. Скорости течений вне слоя сначала также несколько убывают, пока его вязкость остается меньшей т| -510". Однако при дальнейшем увеличении вязкости скорости течений вне слоя практически перестают уменьшаться. Таким образом, высоко вязкий слой ведет себя по отношению к основной жидкости как квази-твердое включение, в которое течения основной жидкости уже практически не проникают. При таких больших виз костях: высоковязкий слой по отношению к основной жидкости играет роль непроницаемой границы с прилипанием (с нулевыми значениями нормальных и тангенциальных компонент скоростей). Как видно из рис,6, в модели В скорости течений основной области на границе с высоковязким слоем практически равны нулю. Па

нижней границе со скольжением, наоборот, горизонтальные скорости максимальн ы.

Э_2 ' :4 4 it ' 'if " 't " ¿2 o'l ^ C6 M 1 ° 0? <M a it ii I

Рис Ь. Рассчитанные мгновенные скорости конвективных течений в жидкости с высоковятким слоем (выделенным темным тоном). Безразмерная вязкость жидкости г|=1. В модели А вязкость слоя т|=10, модели Б вязкость г|=Н) . в модели В вязкость i]=I0 . Но сравнению с моделью А масштаб скоростей в модели Б и В увеличен в 2 раза. Скорости течений внутри слои убывают обратно пропорционально его вязкости. Скорости течений вне слоя сначала также несколько убывают, пока его вязкость остается меньшей i] - 5' 10'

Качественно поведение высоковязких включений представляется очевидным. Однако количественный критерий, а именно значение вязкости, при котором включение начинает вести себя по отношению к основной жидкости как квази-твердое тело, впервые было найден в работе [Moresi, Solomatov. 1995]. Этот результат в настоящее время широко отмечается в литературе и включается в учебники, т.к. имеет непосредственное отношение к вопросу о влиянии литосферных плит на мантийную конвекцию.

Пониженная вязкость

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

И';', - , ,",.,- - -

. , ... .. - ЖЖ... Ж И

; Ц-. : . ■■

-■ Е-.-:-. : , :, Ш ;.".•.'.■

0.1 0 5 0.5 0-7 0.3 А

Рис,7. Мгновенные скорости конвективных течений в жидкости с низковязким слоем (выделенным светлым тоном). Безразмерная вязкость основной жидкости 1] -1. в модели А вязкость слоя т)-0.1, в модели Б вязкость слоя т|" 10 в модели В — г|=10 '. По сравнению с моделью А масштаб скоростей н модели Б уменьшен в 3 раза, в модели В — и 10 раз. Скорости в основной области перестают изменяться при убывании вязкости слоя ниже ^ = 510"4.

Скорости течений внутри слоя возрастают обратно пропорционально его вязкости. Скорости течений вне слоя сначала также несколько возрастают, пока его вязкость остается большей Т] ~ 5*10 . Однако при дальнейшем уменьшении вязкости этот слой ведет себя по отношению к основной вязкости как квази-идеальная жидкость, которая уже практически не влияет на течения основной жидкости, играя роль свободной проницаемой границы.

На рис. 8 приведены поля конвективных скоростей при значениях вязкости слоя, близких к критическим значениям г|=5-103 (модель А) и г|=5-10_1 (модель Б).

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

Рис,8, Структура конвективных течений в жидкости с высоко- и низковязким слоем (выделенных соответственно темным и светлым тонами). Безразмерная вязкость основной жидкости п = !- в модели А вязкость слоя t]=5 !03. в модели Б вязкость слоя г|=5-Ю'4. Масштаб скоростей в высоковязком слое по сравнению с основной жидкости в модели А увеличен к 2000 pat. В модели Б масштаб в высоковязкой основной жидкости по сравнению с низковязким слоем увеличен в 2 раза. R модели А высоковязкий смой играет роль непроницаемой стенки с прилипанием, в модели Б низконизким слой играет роль свободной границы. При тшм в самом низковязком слое возникает автономная конвекция.

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

Их совместное влияние можно проиллюстрировать на трехслойной модели.

На рис. 9 представлены результаты расчета конвективных течений в вязкой жидкости с заданным распределением температуры) при Ra=1000 и Т0 =0.001. Заметим, что если распределение температуры фиксировано, то и структура конвективных течений не меняется во времени. Как видно из рис.8, в модели А конвективные течения охватывают основную мантию и астеносферу, практически не входя в литосферу. Хотя вязкость верхнего слоя невелика и всего в раз 10 превосходит вязкость основной жидкости, но. поскольку ее скачок на границе астеносфера-литосфера составляет 4 порядка, течения в литосферу не приникают. Она ведет себя как квазн твердая крышка. Перепад вязкости на границе мантия-астеносфера составляет три порядка, что несколько ниже критического. Поэтому конвективные течения астеносферы и нижележащей мантии взаимосвязаны.

Рис.9. Структура течений в трехслойных моделях. Толщины верхнего и среднего слоя равны 0.15. Вязкость основной жидкости равна т|=[. а верхнего лнтосферного слоя 1]=10. Вязкость среднего астеносферного слоя а модели Л Г1=10~3, а в модели в г\=НГЛ, что соответственно выше и ниже критического .значения для возникновения автономной конвекции в астеносфере. По сравнению с астеносферным слоем масштаб скоростей в основной мантии увеличен в 2 раза, а и литосфере в 100 раз.

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

Модели с высоковязкими и низковязкими слоями были просчитаны для различных значений вязкости в диапазоне от 10~8 до 10':. На рис.10 представлены зависимости среднеквадратичных скоростей течений во внешней области V, (рис. 10/1} и внутри высоковязкого слоя V-. (рис. /О/Г) в зависимости от значения вязкости слоя.

¡004

'°9 V.

|°д "Л

« 1од Т|

Рис 10 Среднеквадратичные скорости течений в области в основной жидкости V] вне аномально вязкого слоя (рис 10А) и внутри его У2 (рис 105) в зависимости от вязкости слоя Внутри слоя скорости конвективных течений изменяются примерно обратно пропорционально вязкости Во внешней области скорости течений перестают зависеть от вязкости аномального вязкого слоя при критических значениях его вязкости т)= 104 и Г| =

ю-4

Как видно из рис 10, при изменении вязкости верхнего слоя можно выделить три режима конвекции В области значений относительной вязкости слоя 104>т1>10"4 конвективные течения внутри и вне слоя существенно влияют друг на друга При очень высокой вязкости или очень низкой вязкости слоя конвекции внутри слоя и вне слоя становятся как бы независимыми Скорости внутри слоя изменяются примерно обратно пропорционально его вязкости, а скорости вне слоя перестают изменяться при дальнейшем изменении вязкости слоя При этом аномально вязкий слой при относительной вязкости г)>104 играет роль непроницаемой стенки с прилипанием, а при т| < 10"4 играет роль свободно проницаемой границы

Выводы

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

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

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

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

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

Глава 4. Затягивание океанической коры в зонах субдукции

В главе речь идет о затягивании коры и океанической литосферы в мантию

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

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

В современной Земле нормальная океаническая кора толщиной 7 км свободно затягивается в мантию При внедрении мантийных плюмов в движущуюся плиту на литосферной плите образуются базальтовые плато, т е участки утолщенной океанической коры Наибольшую толщину имеют базальтовые плато Онтонг Джава и Карибское плато с толщинами до 30 км Эти плато находятся на океанических плитах и в настоящее время, несмотря на большую толщину, все-таки медленно погружаются в мантию В то же время при толщине коры всего 20 км островные дуги плавают, не погружаясь в мантию

Погружение базальтовых плато и континентальной коры в зонах

субдукции.

t=0 110

t=0115

10 20 30 40 50 10 20 30 40 50 10 20 30 40 50 10 20 30 40 50 X

10 20 30 40 50 10 20 30 40 50 10 20 30 40 50 10 20 30 40 50

50 40 30 20 10

10 20 30 40 50 10 20 30 40 50 ю 20 30 40 50 10 20 30 40 50

Рис 11 Рассчитанная эволюция движения легкой океанической коры (показанной темными маркерами) относительной толщиной 0 01 при тепловой мантийной конвекции для безразмерного времени to=0 100, t(=0 105, t2=0 110, t3=0 115 при различных значениях плотности, характеризуемых значениями композиционного числа Рэлея Re, равными О, -3 106и -6 106

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

Л = 200ехр[-9.2Т+2 3(1-2)], (22)

где Т - температура и р - давление, в безразмерных единицах зависящее от глубины по закону р=1-г. После установления тепловой конвекции в жидкость была добавлена вторая компонента, описывающая океаническую кору Аналогично конвективное движение второй компоненты, описываемой маркерами, характеризуется композиционным числом Рэлея Яс= Дрд03/к%, где Ар — разность плотностей дополнительной и основной жидкости Толщина океанической коры задана равной (1/0=0 01 от толщины всей области, что для модели верхней мантии соответствует толщине 7 км Вязкость коры задана равной вязкости литосферы Как видно на рис 11 океаническая кора вместе с литосферой затягивается в мантию в зоне субдукции, но благодаря своей плавучести в своем движении несколько отстает от литосферы За то же время она меньше погружается в мантию Более легкая кора при 11с=-6-10б и при рассматриваемой относительной толщине (1/0=0 01 еще более отстает в своем движении от погружающейся литосферы и накапливается над зоной субдукции, оставаясь плавать на поверхности плиты, несмотря на постоянное движение нижележащей плиты В мантию могут затягиваться только тонкие размытые слои.

Влияние толщины океанической коры

Ш 20 ЭО 40

10 30 ЭО «> 50

10 20 30 40 90

Рис 12 Рассчитанная эволюция Движения легкой океанической коры (показанной темными маркерами) при тепловой майтийной конвекции при Яа=б 105 для безразмерного времени ^=0 100,11=0 105, 12=0 110, 13=0 115 при значениях плотности, характеризуемом числом Рэлея Яс=-6 106 для толщин коры 0002 (верхняя строка), 0 01 и 01 (нижняя строка) 22

Для осадков с плотностью 2000 кг/м3 плотностное число Рэлея равно 11с ~ —4 10б Таким образом, в рассмотренной модели осадки, несмотря на высокую плавучесть, все-таки затягиваются в мантию в зонах субдукции Как видно из рис 12, даже при высокой плавучести тонкая кора (с!Л>= 0 002,) легко затягивается в мантию С увеличением толщины легкая кора начинает накапливаться над зоной субдукции и только частично затягивается в мантию При относительной толщине с1/Е)= 0 1 кора остается плавать на мантии. При этом она даже не собирается над зоной субдукции, а расплывается в виде слоя вдоль поверхности мантии В манию затягивается только малая часть коры, та часть, которая размывается с ее нижней границы

Влияние утолщений океанической коры

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

¡=0,100 <=0101 «) 102 ч «=0104-1 1=0105

те5 <1=0 01 м>

«=0104-1 1=010

^ 1

1=0 100

Йс«-3 10 й <1=901 >1=0 07

{=0.102

(=0 104

»«0106

Рис 13 Рассчитанная эволюция движения легкой океанической коры толщиной <1=0 01 Верхняя строка — при тепловой мантийной конвекции, Яа=6 105, Яс—3 10б Средняя строка — при тех же параметрах, но с базальтовым утолщением <1=0 07 Нижняя строка — с таким же утолщением, но для более легкой коры континентальной коры при Яс=-б 106

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

Накопление вещества океанической коры на дне мантии

На глубинах 60-100 км базальт океанической коры переходит в эклогит со скачком плотности 15-20% Поэтому вещество коры становится более тяжелым по сравнению с веществом окружающей верхней мантии Эта отрицательная плавучесть эклогита частично сохраняется и в нижней мантии Как впервые было показано в работе Христенсена и Хофманна, вещество океанической коры может накапливаться на дне мантии В последующем (после нагревания) это вещество поднимается с плюмами горячих точек На рис 14 приведены результаты расчета для модели мантийной конвекции при тепловом числе Рэлея В.а=6105, и композиционном числе Рэлея 11с=-8105 для глубины выше 100км (базальтовая кора) и Кс=106 для глубины более 100 км (после перехода в эклогит) Как видно из рис 13, благодаря плавучести базальтовой коры происходит некоторое ее утолщение в зоне субдукции Благодаря утяжелению коры при переходе базальта в эклогит океаническая кора может накапливаться на дне мантии

О 02 04 06 08

О 0.2 04 06 08

0 02 04 06 06 1

02 04 0.6 08

Рис 14 Рассчитанная эволюция движения океанической коры толщиной <1=0 01 с фазовым переходом базальт-эклогит для безразмерного времени 1о=0 100, ^=0 102,12=0 104,1з=0 110

Выводы

В главе 4 исследуется термохимическая конвекция Результаты проведенного моделирования можно применить к конвекции верхней мантии Земли

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

Как видно из рис. 10-14 несмотря на несовершенство рассмотренной простой модели, результаты численного моделирования не только 24

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

Заключение

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

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

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

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

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

1 Трубицын В П, Симакин А Г, Баранов А А

Влияние пространственных вариаций вязкости на структуру мантийных течений

//Физика Земли 2006 №1 С 3-15.

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

3 Трубицын В П, Баранов А А, Евсеев А Н, Трубицын А П, Харыбин Е В. Влияние низковязкой астеносферы на мантийные течения

//Физика Земли 2006 № 12. С 11-19

4 Трубицын В П, Баранов А А., Харыбин Е В

Численные модели субдукции океанической коры с базальтовыми плато //Физика Земли. 2007 №7. С 3-10

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

6 Трубицын В П, Евсеев А Н, Баранов А А, Трубицын А П Влияние параметров фазового перехода и чисел Релея на мантийную конвекцию

//ФизикаЗемли 2008 №2 С 3-12

7 Rogozhma I, Kaban М К., Trubitsyn V, Baranov А (2006) Importance of lateral viscosity variations for a modelmg of geoid and dynamic topography EGU General Assembly,

Vienna, 2-7 April 2006,Geophysical Research Abstracts, Vol 8, 07351

8 Baranov A .A, Trubitsyn V P, Kaban M.K, Rogozhma I. (2007) Effect of strong viscosity variations on the global mantle flow EGU General Assembly, Vienna, 16-20 April 2007, Geophysical Research Abstracts

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

ВВЕДЕНИЕ

Глава 1. Конвекция в мантии Земли.

1.1 Тепловая конвекция в мантии Земли.

1.2 Свойства тепловой конвекции.

Глава 2. Численное моделирование двумерной мантийной конвекции 14 Постановка задачи

2.1 Уравнения мантийной конвекции. Граничные условия.

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

2.3 Точное частное аналитическое решение уравнения Стокса.

2.4 Тестирование различных программ численного решения. 24 Поля температуры скоростей напряжений.

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

3.1 Повышенная вязкость.

3.2 Пониженная вязкость.

Глава 4. Затягивание океанической коры в зонах субдукции

4.1 Введение.

4.2 Погружение базальтовых плато и коры в зонах субдукции.

4.3 Выводы.

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

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

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

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

Цели работы

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

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

Задачи исследования

Цели работы определили постановку следующих задач:

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

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

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

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

5. Проанализировать скорость затягивания и сравнить ее с реальной.

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

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

2. Численный метод основан на использовании конечных элементов.

3. Разработка программного средства велась на языке программирования Фортран. Компилятором и средой разработки был выбран пакет Compaq Visual Fortran 6.6. С помощью разработанного программного обеспечения были изучены модели конвекции с различными параметрами и исследованы задачи массопереноса вещества.

Основные защищаемые положения диссертации

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

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

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

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

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

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

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

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

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

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

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

Обсуждение работы

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

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

1. Трубицын В.П. ,Симакин А.Г., Баранов А.А.

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

Физика Земли. 2006. № 1. С. 3-15.

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

3. Трубицын В.П., Баранов А.А., Евсеев А.Н., Трубицын А.П., Харыбин Е.В. Влияние низковязкой астеносферы на мантийные течения.

Физика Земли. 2006. № 12. С. 11-19.

4. Трубицын В.П., Баранов А.А., Харыбин Е.В.

Численные модели субдукции океанической коры с базальтовыми плато. //Физика Земли. 2007. № 7. С. 3-10.

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

6. Трубицын В.П., Евсеев А.Н, Баранов А.А., Трубицын А.П. Влияние параметров фазового перехода и чисел Релея на мантийную конвекцию.

Физика Земли. 2008. № 2. С. 3-12.

7. Rogozhina I.; Kaban М.К.; Trubitsyn V.; Baranov A. (2006). Importance of lateral viscosity variations for a modeling of geoid and dynamic topography. EGU General Assembly,

Vienna, 2-7 April 2006,Geophysical Research Abstracts, Vol. 8,07351.

8. Baranov A.A., Trubitsyn V.P., Kaban M.K., Rogozhina I.(2007).Effect of strong viscosity variations on the global mantle flow. EGU General Assembly, Vienna, 16 - 20 April 2007, Geophysical Research Abstracts.

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

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

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

4.3 Выводы

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

1-,

0.80.60.40.20

0.80 60.40.2 0

0.80.6 0.40.2-1 0 ,

0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

0 0.2 04 0.6 0.8 1 химические реакции не рассматриваются, и свойства компонент считаются заданными, такую конвекцию для механической смеси обычно называют термо-композиционной.

Результаты проведенного моделирования можно применить к конвекции верхней мантии Земли. Для верхней мантии D=700km, Ро=3600кг/м3, а=3-10"5К"1, g=10M/c2, k=10"6 м2/с, Т1о=Ю21Пас. Эффективный нададиабатический перепад температуры в верхней мантии, соответствующий наблюдаемому тепловому потоку из мантии, равен ДТ~1500К [Trubitsyn et al., 2006]. При этих значениях параметров единицы измерения будут равны для длины D=700km, для скорости k/D=10" 6/7* 105~1.3 • 10"6м/с=4-10"3см/год, для времени D2/k~1.7-1010 лет. Тепловое число Рэлея Ra=apogToD3/kr|o~6T05. Рассчитанная скорость течений тепловой конвекции, в безразмерных переменных составляющая V=250, в Л размерных переменных равна V=250- 4-10" =1см/год. Характерное безразмерное время t=0.005 соответствует размерному времени t=0.005* 1.7Т0Ш=85 млн. лет.

Плотности океанической и континентальной коры соответственно равны р0=3000 кг/м , ро=2800 кг/м . Эффективные вязкости океанической и континентальной коры в первом приближении можно взять равными вязкости океанической литосферы. Плотность и вязкость осадков

1 | с соответственно равны ро=2000 кг/м , Г)0=Ю Пас. Принимая плотность мантии, находящейся непосредственно под корой, равной 3200 кг/м3 [Schubert et al., 2002], получим для перепада плотности, определяющего плавучесть океанической коры, значение Др=200 кг/м3. При этих значениях композиционное число Рэлея для океанической коры будет равно Rc=ApgD3/kr|o~-7-105. Для континентальной коры с перепадом плотности

1 ft Ч

Ар=400 кг/м число Рэлея равно Rc~ -1.4-10 . Для осадков при Др=1200 кг/м композиционное число Рэлея равно Rc~- 4.2-106.

Средняя толщина нормальной океанической коры равна ~7км, а континентальной коры ~40км. Средняя толщина океанических осадков равна -0.5-1.0 км [Schubert et al., 2002]. В относительных единицах для модели конвекции в верхней мантии толщины океанической, континентальной коры и осадков будут соответственно равны 0.01, 0.06 и 0.001. Средние относительные толщины базальтовые плато Онтонг Джава и Карибского плато равны -0.04.

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

Благодаря переходу базальта в эклогит со скачком плотности 15-20% вещество океанической коры может накапливаться на дне мантии.

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

Заключение

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

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

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

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

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

1. Добрецов Н.Л., Кырдяшкин А.Г., Кирдяшкин А.А.

2. Глубинная геодинамика. Новосибирск. Изд-во СОРАН. 2001. 408 с.

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

4. Каракин А.В., Магницкий В.А., Калашникова КВ.

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

6. Докл. АН СССР, 1974, т. 214, № 3, с. 561-564.

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

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

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

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

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

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

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

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

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

16. Трубицын В.П., Баранов А.А., Харыбин Е.В.

17. Численные модели субдукции океанической коры с базальтовыми плато. //Физика Земли. 2007. № 7. С. 3-10.

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

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

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

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

22. Baranov А.А., У.P. Trubitsyn, М.К. Kaban, I. Rogozhina.(2007).Effect of strong viscosity variations on the global mantle flow. EGU General Assembly,

23. Vienna, 16-20 April 2007, Geophysical Research Abstracts.

24. 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. 199231,1988.

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

26. 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.

27. A benchmark comparison for convection codes// Geophys. J. Int. 98, P. 23-38.

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

29. Brunei D., and P. Machetel, 1998. Large-scale tectonic features induced by mantle avalanches with phase, temperature and pressure lateral variations of viscosity, J. Geophys. Res., 103,4929-4945.

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

31. 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.

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

33. 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.

34. 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.

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

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

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

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

39. 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.1.P.205-210.

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

41. Gurnis, M. and S. Zhang, Generation of long-wavelength wavelength heterogeneity in the mantle by the dynamic interaction between plates and convection, Geophys. Res. Lett., 18, 581-584,1991

42. Gurnis, M, J. Ritsema, H. van Heijst, andS. 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.

43. Hansen, G. Jarvis, M. Koch, G. Marquart. D. Moore, P. Olson, H. Schmeling, andT. Schnaubelt, 1989.

44. A benchmark comparison for convection codes, Geophys. J. Int., 98, p, 23-38.

45. Hawkesworth C.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,256-268.

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

47. Khodakovskii G., Rabinowicz M., Ceuleneer G., Trubitsyn V.P., Melt percolation in a partially molten mush: effect of a variable viscosity,

48. Earth and Planetary Sci. Lett. 1995,, V. 134, P. 267-282.

49. Korenaga J. Firm mantle plumes and the nature of core-mantle boundary region

50. Earth Planet. Sci. Lett. 2005. V. 232. P. 29-37.1.nardic A. and W.M. Kaula, 1995. Mantle dynamics and the heat flow into the Earth's continents, Nature, 378,709-711.

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

52. McNmnaraA. К., and Zhong S., Degree-1 mantle convection: Dependence on internal heating and temperature-dependent rheology, GRL, 32, L01301, 10.1029/2004GL021082,2005

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

54. 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.

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

56. 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.

57. 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.

58. Pelletier D.Fortin A.Camarero R.

59. Are FEM solutions of incompressible flows really incompressible? // IJ. for Numer. Methods in Fluids. 1989.V.9.P.99-112.

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

61. Rogozhina /.; Kaban M.K.; Trubitsyn V.; Baranov A. (2006). Importance of lateral viscosity variations for a modeling of geoid and dynamic topography. EGU General Assembly, Vienna, 2-7 April 2006,Geophysical Research Abstracts, Vol. 8, 07351.

62. Rudman M., "Two-Phase Natural Convection for Crystal

63. Settling in Magma Chambers," Phys. Earth Planet. Inter. 72,153-172 (1992).

64. 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.

65. Schubert G., D.L. Turcotte and P. Olson.

66. Mantle Convection in the Earth and Planets. Cambridge Univ. Press. 2001. p. 940.

67. A.Simakin, H. Schmeling, V. Trubitsyn,

68. Convection in melts duetosedimentative crystal flux from above

69. Phys. Eart Planet. Inter., 1997. V.102, P. 185-200.

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

71. 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. P. 265-282.

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

73. Tan E., P. Thoutireddy, E. Choi, M. Gurnis, and M. Aivazis. 2000.

74. GeoFramework Parti: Coupling models of mantle convection with Pythonframework. Geochemistry, Geophysics, Geosystems.

75. Trubitsyn V. P. andE. V. Kharybin, "Convective Instability of the Sedimentation

76. Regime in the Mantle," Izv.

77. Akad. Nauk SSSR, Fiz. Zemli, No. 7, 21-30 (1987).

78. 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.

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

80. 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, 3151989a).

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 V. V. Kubyshkin, "Thermal Convection in the Mantle Due to Horizontal and Vertical Gradients of Temperature," Fiz. Zemli, No. 5,12-23 (1991).

83. Trubitsyn V. P. and E. V. Kharybin, "Convection in

84. Magma Chambers Due to an Inversion of the Depth Distribution of Settling

85. Crystals," Fiz. Zemli, No. 5,47-52 (1997).

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

87. Uyeda D. and H. Kanamori, 1979. Back-arc opening and the mode of subduction. J. Geophys.1. Res., 84,1049-1061

88. Zhang S. and Christensen U. Some effects of lateral viscosity variations on geoid and surface velocities induced by density anomalies in the mantle//

89. Geophys .J. Intern. 1993.V. 170. P. 531-547.

90. Zhang, 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.

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

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

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

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

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

96. 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

97. 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.

Информация о работе

Баранов, Алексей Андреевич
кандидата физико-математических наук
Москва, 2007
ВАК 25.00.10

Диссертация

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

Автореферат

Похожие работы