CC BY
24
Вычислительные технологии
Том 19, № 5, 2014
Численное моделирование трёхмерной конвекции в верхней мантии Земли под литосферой Евразии*
В. В. ЧЕРВОВ1, Г. Г. ЧЕРНЫХ2'3, Н.А. БУШЕНКОВА1, И.Ю. КУЛАКОВ1
Построена трёхмерная численная модель тепловой конвекции в верхней мантии Земли в сферических переменных с использованием метода искусственной сжимаемости. Представлены результаты трёхмерного моделирования конвекции под кра-тонами Евразии. Результаты расчётов демонстрируют структуру конвективных потоков. Построена численная модель погружающейся под континент океанической плиты.
Ключевые слова: трёхмерное численное моделирование, сферические переменные, мантия Земли, тепловая конвекция, кратон, литосфера, субдукция, спрединг.
Введение
Современное численное моделирование конвекции в мантии Земли в сферических координатах представлено в основном работами по изучению процессов, протекающих во всей мантии и во внешнем ядре с учётом фазовых границ между слоями [1-8]. Анализ цитированных публикаций показывает, что результаты численного моделирования конвекции под кратонами Евразии отсутствуют. В настоящей работе, посвящённой конвекции под континентальной литосферой, область исследования ограничена сферической областью верхней мантии (т. е. фазовой границей, находящейся на глубине 670-700 км) от нулевого меридиана до восточных границ Евразии (145° в.д.) и от экватора до северных границ материка (80° с.ш.). В эту область в качестве объектов исследования входят континентальные платформы, неоднородные по возрасту и мощности. В [9] в декартовых координатах исследовалась область Центральной Азии, основным объектом которой был Сибирский кратон. Тарим и Северо-Китайский кратоны замыкали границы области. При переходе к расчётам всего Евразийского континента возникла необходимость в учёте влияния сферичности. В настоящей работе область исследования включает Сибирский кратон и окружающие его Западно-Сибирскую, Русскую и Китайскую плиты, а также Тарим, комплекс Тувинских кратонов, Центрально-Азиатский складчатый пояс и часть Индийской плиты.
Представленное исследование посвящено разработке трёхмерной численной модели тепловой конвекции в верхней мантии Земли под кратонами Евразии и численной модели погружающейся под континент океанической плиты. Для построения модели субдукции оказалось достаточным задать на некотором участке поверхности Земли наблюдаемые скорости расхождения литосферных океанических плит. В ходе исследо-
1 Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, Новосибирск, Россия
2Институт вычислительных технологий СО РАН, Новосибирск, Россия
3Новосибирский государственный университет, Россия Контактный e-mail: [email protected]
* Работа выполнена при поддержке СО РАН (проекты № 20 и 76) и РФФИ (инициативный научный проект № 13-05-00054).
ваний использовалась концепция Флейто — Йена [10, 11] в совокупности с концепцией Трубицына — Тычкова [12, 13]. Суть этих представлений о формировании структурных элементов заключается в том, что в первом случае движения плит, например, как в нашем случае — спрединг и субдукция, — должны возникать в результате конвективных течений в мантии, в которой литосферные плиты, как океанические, так и континентальные, рассматриваются как жидкость с повышенной (по сравнению со всей мантией) вязкостью. Согласно же концепции Трубицына — Тычкова континенты не участвуют в конвекции и ведут себя как твёрдые тела, погружённые в жидкую мантию. В рассматриваемой постановке задачи океанические плиты моделировались сильновязкой жидкостью (по Флейто — Йену), континентальная плита моделировалась твёрдым телом.
Работа является продолжением и развитием [9].
1. Математическая постановка задачи
Конвективные течения в верхней мантии Земли описываются известной математической моделью, включающей в себя обезразмеренные уравнения в сферических координатах [5, 14]:
1 -V * 1 д (sin 9Ve) 1 d(r2Vr) „
+ . „ „„—¿ + —¿ = 0, (1)
r sin 9 др r sin 9 д9 r2 dr
d_ д<р
2r¡r
дУ *
( -р
+ cos 9Vе + sin 9Vr
д_ -r
r¡r3 sin2 9
дV * -r
д_ + д9 1
2 9fдV* vrsin 9\~df
+ — ^ - 1V *
r sin 9 др r
1 дVв \
+ —fí - ctg 9V* sin 9 др J
2 ■ a-P = r sin 9^—,
др
+
д_
др
r¡r
sin 9 \ др
дVе -V*
+ sin 9-^- cos 9V*
д9
д_ + д9
. дVе r¡r sin 9 ( + Vr
д9
+
д_ -г
Пг3 sin 9 +2n r cos 9
дVе 1 дVr 1
+ —
-г r д9
- - V' r
+
-Vr -Vе ЛГГ r—--+ + Vr
-r
-9
2 • a-P
r sin 9W
d_ др
П
-Vr
sin 9 др -r
+ sin 9
-V * r-Z--V *
2rqr2 sin 9
-r -V r -r
+ 2n sin 9
+ д9 -V r
П sin 9
-V r ~39
-r
+r
-p
-V е -r
-Vе
+
r sin 9 I —--RaT
-r
)
-T 1 --t r sin 9 др
(V* T) +
1
- 2T
+
1 д_ r sin 9 д9 1 -
1д
(sb 9V е T) + T2 Q-rir2V r T)
r2 sin2 9 др2 r2 sin9 д9
• 9-t sin 9
1
r2 -r
,-t
-r
Здесь Т — температура, р — дефект давления, V*, V6, Vr — компоненты вектора скорости V, Ra = RaR — число Рэлея (Ra = agzpR3ДТ/п0х), R — радиус Земли, ДТ — разность между температурами на подошве мантии и на поверхности, х — коэффициент температуропроводности, a — коэффициент теплового расширения, p, п — характерные плотность и динамическая вязкость, gz — ускорение силы тяжести.
Тензор вязких напряжений определён формулами [14]
т
**
0 / 1 dV * Vr V6 ctg 6\
= 2п —— — + — +-
\т sin6 r r /
/1 dVr dV6 V6 \ Tre = v\z—г + —----),
r d6 dr Te* = П
« /1 dV6 гее = 2n -
V Л _ dVr r d6 r / rr dr
Tr* п
'dV * 1 dVr V *\ r sin 6 r y '
dr
1 dV6 1 dV * +
r sin 6 r d6
V*ctg 6 \ r
Обезразмеренные переменные вводятся следующим образом:
r = Rr, t = — i, T = ДТТ, п = П0 П, Р = V = X V.
X R2 R
(6)
На рис. 1 схематично изображена расчётная сферическая область. Для вектора скорости на боковых гранях задаются условия проскальзывания, а на нижней и верхней — условия прилипания. Для температуры на боковых гранях ставятся условия теплоизоляции (адиабатическая стенка), т. е. первые нормальные производные на вертикальных стенках равны нулю. На верхней и нижней гранях ставятся условия Дирихле: нулевая температура на верхней и 1800 0 С (в обезразмеренных переменных равная единице) на нижней:
Рис. 1. Схема расчётной сферической области (с кратоном). Представлены граничные условия для компонент скорости и температуры. Мощность литосферной плиты составляет 120 км, литосферы с кратоном — 240 км
на поверхностях p = pi, p = p2 (9i < 9 < 92, ri < r < r2)
V v = dve = ^ = дт = 0,
dp dp dp '
на поверхностях 9 = 9^ 9 = 92 (pl < p < p2, rl < r < r2)
dV! - V-ctg 9 = v »= dV: = dT = 0
d9 V ctg 9 V 89 d9
На подошве верхней мантии ставятся условия прилипания
V! = Vе = Vr = 0, r = ri, pi < p < p2, 9i < 9 < 92.
Если литосфера не выделена, т. е. мантийная жидкость занимает всю расчётную область, граничные условия прилипания ставятся на верхней поверхности:
V! = Vе = Vr = 0, r = r2, pi < p < p2, 9i < 9 < 92,
а в случае присутствия литосферы — на подошве литосферы, при этом скорости внутри литосферной плиты соответствуют скоростям, заданным на границах литосферного блока.
Для температуры на горизонтальных поверхностях расчётной области были поставлены граничные условия Дирихле: нулевые на дневной поверхности и равные единице на подошве верхней мантии:
T = 0, r = r2, pi < p < p2, 9i < 9 < 92,
T =1, r = ri, pi < p < p2, 9i < 9 < 92.
В начальный момент времени t = t0 задаются начальные условия лишь для температуры: T(p, 9, r, to) = To(p, 9, r).
2. Алгоритм расчёта
Для построения численной модели применялся метод установления с использованием неявной схемы искусственной сжимаемости [15-17] и метода дробных шагов [16, 18]. Схематично алгоритм запишем в виде
dpH+1 ~2i-
dt
f
—c2div(Vn ), (7)
dVi d
+ Vipn+1 = 7—rffc+1 + RaT nßi, dtf dxk
dT
+ V (Vn+1T) = V2T, (9)
1 d
где em — компоненты вектора с координатами e = {0, 0,1}; Vm = -— ——, hm — коэф-
hm dxm
фициенты Ламэ: h1 = r sin 9, h2 = r, h3 = 1; t — время, tf — фиктивное время; n — номер шага по времени t; i = 1, 2, 3, m = 1, 2, 3, k = 1, 2, 3.
Численная реализация (7)-(9) включает в себя следующие этапы:
1) в исследуемой области задаётся начальное распределение температуры, удовлетворяющее граничным условиям. Дефект давления р и компоненты вектора скорости V полагаются нулевыми на первом слое по времени, а далее берутся с предыдущего слоя по времени;
2) из векторного уравнения (8) находится поле скорости Vra+1 с привлечением уравнения (7);
3) путём решения уравнения (9) вычисляется поле температуры.
Процесс повторяется до некоторого значения £ = NД£, где N — число слоёв по времени.
Конечно-разностный алгоритм решения задачи основан на применении метода дробных шагов [16, 18]. Для реализации (8), (9) использовалась схема стабилизирующей поправки. Для решения задач вводилась равномерная в каждом направлении пространственная сетка.
На каждом дробном шаге схемы (7)-(9) применялись скалярные трёхточечные прогонки. Метод искусственной сжимаемости применительно к задачам конвекции в верхней мантии Земли в декартовой системе координат реализован в работе [19].
В процессе решения задачи проводилась проверка сходимости применяемого алгоритма. Для этой цели привлекались модельные задачи для конечно-разностных аналогов уравнений (1)-(5), имеющие аналитическое решение. Одним из основных этапов тестирования является интегрирование уравнения теплопроводности (5). Для тестирования разностных аналогов уравнений (2)-(4) применялось интегрирование модельных уравнений Пуассона для компонент скорости и давления.
Итак, тестирование численной модели осуществлялось путём решения модельных задач. Тестирование же численного алгоритма в целом проведено в [19] путём решения задачи конвекции в мантии Земли в декартовой системе координат. Рассмотренная там модельная задача является международным тестом [20].
3. Численное моделирование трёхмерной конвекции под литосферой Евразии
Как отмечалось выше, основанное на методе слабой сжимаемости в декартовой системе координат исследование трёхмерной конвекции в модельной постановке проведено в [19]. В настоящей работе выполнено моделирование тепловой конвекции под внут-риконтинентальной областью Евразии, в которую в виде литосферных блоков входят Русская платформа, Западно-Сибирская плита, Сибирская платформа, Центрально-Азиатский складчатый пояс, Тувинский комплекс микрократонов, Тарим, Китайские платформы, Индийская и Аравийская плиты (рис. 2).
Краевые условия задаются, как и в модельной постановке на рис. 1.
На подошве литосферы при постановке начального распределения температуры учитывается первоначальное значение температуры: Т = 1200 0С. Температура рассчитывалась во всём параллелепипеде: кондуктивно в пределах литосферных блоков и конвективно в остальной области.
Таким образом, движение жидкости, т. е. поле скорости, вычислялось вне литосферы. Число Рэлея, характеризующее режим конвекции, было выбрано как йа — Иая — 152 000 000, что отвечает современным представлениям об условиях в недрах Земли. Число Рэлея, рассчитанное от вертикального размера d конвектирующей области, вы-
Рис. 2. Схема расположения литосферных блоков. Тувинский комплекс микрократонов находится в Центральном складчатом поясе в юго-западном направлении от Сибирской платформы
числяется следующим образом: Rad = Ra (d/R)3 = 201706. Основные параметры задачи в системе СИ, пригодные для верхней мантии, выбирались следующими:
R = 6 371000 м, d = 700 000 м, AT = 1800 °C, х = 10-6 м2/с, а = 2 ■ 10-5 °С, р = 3300 кг/м3, gz = 10 м/с2, По = 2.025 ■ 102i кг/(м ■ с). (10)
Вязкость мантийного вещества задавалась в виде
П (p, 9, r) = exp (br — aT (p, 9, r, t)). (11)
Здесь параметры a = 3.89 и b = 5.84 обеспечивают перепад вязкости Птах/Птт от 20 до 200, что характерно для верхнемантийных характеристик течений.
В уравнении (7) полагалось с2 = 25 000, величина шага по фиктивному времени Atf = min(Ap, A9, Ar)/с.
Моделирование показало, что, как и в случае прямоугольных в плане кратонов [12, 21, 22], реальные кратоны Евразии порождают аналогичные структуры. Наблюдаются устойчивые восходящие потоки в виде плюмов, а также нисходящие потоки и прогретые области по периферии кратонов (рис. 3-5). Перенос мантийного вещества от оснований кратонов к верхним горизонтам (обтекание) проявляются в виде мелкомасштабной моды конвекции вдоль бортов кратонов (см. рис. 3, 4).
Представленная трёхмерная численная модель конвекции обнаруживает повышение средней мантийной температуры под кратоном на 100 °C.
По геолого-геофизическим данным [23] в районе южнее Сибирского кратона и севернее Тарима и Северо-Китайской платформы мощность литосферы составляет от 40 до 95 км. В численной модели толщина литосферного блока в указанном районе принималась равной 80 км.
В результате численного моделирования было показано, что в зоне ловушки (утонение литосферы), как правило, преобладают нисходящие потоки холодного мантийного
Рис. 3. г#-сечение температурного поля плоскостью ^ = 450. Плоскость пересекает Аравийский кратон и Русскую платформу. Наблюдаются устойчивые восходящие потоки (1200-14000С) под Аравией и Русской платформой и прогретые области по периферии кратонов
Рис. 4. г#-сечение температурного поля плоскостью ^ = 300. Плоскость пересекает Русскую платформу. Наблюдается перенос мантийного вещества от основания кратона к верхним горизонтам (обтекание), что проявляется в виде мелкомасштабной моды конвекции вблизи бортов кратона
Рис. 5. Сечение температурного поля на глубине 320 км, £ = 400 млн лет. Кратоны обозначены сплошной линией. В области ловушки (Центральный складчатый пояс) наблюдаются вЁоснов-ном нисходящие потоки (900-10000С). Наиболее прогретые области (1200-14000С) установлены под кратонами
вещества. На глубине 320 км обнаруживается достаточно холодное (800-900 0С) мантийное вещество (рис. 5). Следует заметить, что под территорией Западно-Сибирской плиты, где мощность литосферы составляет 120 км, комплекс нисходящих потоков имеет температуру в среднем на 100 0С выше.
4. Формирование слэба
В рамках задачи моделирования конвекции на основе уравнений Навье — Стокса в приближении Обербека — Буссинеска и геодинамическом приближении в сферических координатах (1)-(5) была построена модель погружающейся под континент плиты. Для построения модели субдукции оказалось достаточным задать на некотором участке поверхности Земли наблюдаемые скорости расхождения литосферных океанических плит. Известно, например, что в районе срединно-атлантического медленно-спредингового хребта [24] наблюдаемые скорости составляют 2-4 см в год.
В задаче скорость спрединга задавалась даже меньше наблюдаемых на планете значений [24]. Однако расчёты показали, что в течение короткого геологического промежутка времени (5-10 млн лет) сформировался мощный восходящий поток, устремлённый к области раздвижения плит. В районе сочленения океанической и континентальной плит конвективные потоки, спровоцированные спредингом, немедленно и естественным образом порождают нисходящий конвективный поток, который можно интерпретировать как погружающуюся под континент плиту (слэб). Угол погружения зависит от угла наклона континентальной плиты, на передней части которой задавались либо естественные условия скольжения, либо скорости, рассчитанные от скорости раздвижения.
В рассматриваемой задаче, как уже отмечалось во введении, использовалась концепция Флейто — Йена [10, 11] в совокупности с концепцией Трубицына — Тычкова [12, 13]. Суть этих взглядов на формирование структурных элементов заключается в том, что в первом случае движения плит, например, как в нашем случае — спрединг и субдук-ция, — должны возникать в результате конвективных течений в мантии, в которой литосферные плиты, и океанические, и континентальные, рассматриваются как жидкость с повышенной (по сравнению со всей мантией) вязкостью. В случае концепции Трубицына — Тычкова континенты не участвуют в конвекции и ведут себя как твёрдые тела, погруженные в жидкую мантию.
В нашей задаче океанические плиты моделировались сильновязкой жидкостью (по Флейто — Йену), континентальная плита моделировалась твёрдым телом.
Постановка задачи. В узкой сферической области мантии (рис. 6) рх < р < р4, Вх < В < В2, тх < т < г2, рх = 0°, р4 = 60°; Вх = -2.5°, В2 = +2.5°; тх = 3571 км, г2 = 6371 км (тх = 5671 км для верхней мантии) моделируется спрединг (расхождение плит). Справа от оси спрединга (р2 = 15°) расположен континент, протяжённость которого составляет Л/2, Л = р4 — рх = 60°. Координаты континента на поверхности
Мощность континента 220 км. Угол наклона левой континентальной грани j задан равным 45°.
Граничные условия ставились следующими.
На вертикальных границах Г^, Г front, Г^аск для вектора скорости поставлены условия проскальзывания, а для температуры — условия симметрии (адиабатическая стенка):
на границе rfe/t : Опр = = 0 ^
r = r2, p3 < p < p4, 91 < 9 < 92, p3 = 30°, p4 = 60°.
V р
8Ve dVr dT
r
0, p = pi, ri < r < r2, 9i < 9 < 92,
dp dp dp
Г2 -1 CM/ГОД <Р2 +1 СМ/ГОД % Г2
0° 5° 10° 15° 20р 25° 3(Г 35° 40° 45а 5(Г 55° 60° Р
Рис. 6. Схема раздвижения литосферных плит со скоростью 2 см/год. Шкала расстояний здесь (а также ниже на рис. 7, 8) представлена отсчётом от дневной поверхности: 0 км соответствует расстоянию г2 = 6371 км от центра Земли, -2800 км — п = 3571 км от центра Земли
на границах Г front и ГЬаск : аге = aev = 0 ^ dVv dVr dT
Vе = -yj- - VVctg в = = = 0, в = в1 , в = 02, П < Г < Г2, P1 < V < V4.
На правой вертикальной границе rright для вектора скорости установлены условия протекания, а для температуры, как и в случае противоположной границы, — условия симметрии:
BVV dT
-К— = Vе = Vr = —- = 0, V = V4, Г1 < Г < Го, 01 < в < 02. dp dp
На всех остальных границах приняты условия прилипания:
а) под океанической литосферой
Vr = Vе = 0, г = г2, V1 < V < V3, в1 < в < в2,
VV = —u0 = —1 см/год, г = г2 , V1 < V < V2, в1 < в < в2, VV = +u0 = +1 см/год, г = г2 , v2 < V < V3, в1 < в < в2;
б) на нижней горизонтальной границе всей расчётной области
Vr = Vе = VV = 0, г = Г1, V1 < V < V4, в1 < в < в2;
в) под континентом
Vr = Vе = VV = 0, г = Го, Vo <V < V4, в1 < в < в2,
где V0 = Vo (180°/п), Vo = v3 + (г2 — r0) ctg (Y), причем переменные v', г' , выраженные в радианах, получены из размерных по формулам v' = (п/180°)v; г' = г/R; V0 — координата вершины тупого угла (v0, в, r0) в фигуре континента, r0 — координата подошвы континентальной плиты, находящейся на глубине 220 км: r0 = (6371 — 220) км = 6151 км (см. рис. 6);
г) на боковой границе континента, расположенной под углом y
Vе = 0, V* = u0 cos y, Vr = u0 sin y; r0 < r < r2, р3 <р < р0, 01 < 0 < 02;
д) на линии срединно-океанического хребта
Vr = Vе = V* = 0, r = Г2, р = Р2, 01 < 0 < 02;
е) на линии субдукции
Vr = Vе = 0, V* = 0.5u0, r = Г2, р = P3, 01 < 0 < 02.
На горизонтальных поверхностях расчётной области были поставлены граничные условия Дирихле для температуры: нулевые — на дневной поверхности и равные единице — на подошве мантии. В случае верхней мантии это 1800 °C, в случае нижней — 3700 °C, следовательно
T = 0, r = r2, р1 < р < р4, 01 < 0 < 02,
T = 1, r = r1, < р < р4, 01 < 0 < 02.
Расчёты температурного поля проводились сквозным счётом через границы между континентальной литосферой и мантийной жидкостью.
Значения параметров задачи для верхней мантии использовались из (10), для всей мантии — принимались следующими:
R = 6 371000 м, d = 2 800 000 м, AT = 3 700 °C, х = 10-6 м2/с, а = 10-5 °C, р = 4000кг/м3, gz = 10 м/с2, П0 = 3.375 ■ 1022 кг/(м ■ с). (12)
Вязкость задавалось с помощью стандартной формулы (11).
На рис. 7 приведены результаты моделирования спрединга и субдукции в верхней мантии Земли. Структура зоны субдукции получена в результате естественного процесса раздвигания плит и влияния геометрии континентальной плиты. Затенённая область (сгущение изотерм в пределах р3 < р < р' = 35° в центре фигуры) соответствует тонущей части океанической литосферной плиты (слэбу). Слэб в нашем случае представляет собой конвективный нисходящий поток мантийного вещества, температура которого находится в пределах 0-1200 °C.
Вначале задача решалась для верхней мантии. Условия, поставленные в данном случае на подошве верхней мантии, не позволяют погружающемуся конвективному потоку
Рис. 7. гр-сечение поля температуры, рассчитанного в верхней мантии. Время — 70 млн лет. На нижней границе (г = 5671 км) температура равна 1800 0С, на верхней (г = 6371 км) — 0 0С. Длина расчётной области составляет 60°. Угол наклона континентальной плиты равен 450. Расчёты выполнены на сетке с числом узлов 664 х 11 х 71
Рис. 8. гр-сечения полей температуры, рассчитанных для всей мантии: а — время 70, б — 240, в — 270 млн лет
Таблица 1. Динамика изоповерхности T = 1200 °С
Глубина 800 км 925 км 1175 км 1300 км 1400 км
Время 70 млн лет 120 млн лет 200 млн лет 240 млн лет 270 млн лет
Скорость 1 см/год 0.25 см/год 0.312 см/год 0.312 см/год 0.33 см/год
Таблица 2. Динамика изоповерхности T = 800 °С
Глубина 700 км 825 км 1025 км 1125 км 1205 км
Время 70 млн лет 120 млн лет 200 млн лет 240 млн лет 270 млн лет
Скорость 1 см/год 0.25 см/год 0.25 см/год 0.25 см/год 0.27 см/год
проникнуть за пределы верхней мантии (см. рис. 7). Далее задача решалась по всей мантии Земли (см. рис. 6). Для расчётов течений во всей мантии таких ограничений, о которых говорилось выше (рис. 8), нет. Слэб достигает границы примерно за 70 млн лет, т.е. средняя скорость равна 10 км за 1 млн лет (что соответствует 1 см/год), заданной для океанической плиты на поверхности. Совпадение скоростей указывает на то, что слэб вместе с плитой на поверхности движутся как единое твёрдое тело (см. рис. 8).
Прослеживая рассчитанную горячую границу слэба (T = 1200 °C, табл. 1), видим, что за первые 70 млн лет изоповерхность 1200 °C сместилась с глубины 100 км до глубины 800 км (средняя скорость движения составляет 10 км/млн лет, рис. 8, а); за следующие 50 млн лет продвинулась в глубь нижней мантии на 125 км (средняя скорость 2.5 км за 1 млн лет), а к отметке 1300 км изоповерхность 1200 °C подошла через 120 млн лет, имея среднюю скорость 3.125 км за 1 млн лет, т.е. слэб за счёт конвекции внутри нижней мантии получает на этом участке некоторое значительное ускорение (рис. 8, б). Скорость нарастает до 3.333 км/млн лет, так как следующие 100 км изоповерхность преодолевает за 30 млн лет (рис. 8, в).
Динамику изоповерхности T = 800 °С можно проследить по табл. 2.
Основные результаты работы сводятся к следующему. Построена трёхмерная численная модель конвекции в мантии Земли. Приведены результаты численного моделирования конвекции под кратонами Евразии в сферической системе координат. Построена численная модель погружающейся под континент океанической плиты. Представлены результаты численного моделирования спрединга и субдукции в узкой сферической области.
Настоящую работу авторы посвящают светлой памяти Сергея Анатольевича Тычкова.
Список литературы
[1] Peter van Keken. Cylindrical scaling for dynamical cooling models of the Earth // Phys. of the Earth and Planetary Interiors. 2001. Vol. 124. P. 119-130.
[2] Kageyama A., Yoshida M. Application of the Yin-Yang grid to a thermal convection of a Boussinesq fluid with infinite Prandtl number in a three-dimensional spherical shell // Geophys. Res. Lett. 2004. Vol. 31. L12609. doi:10.1029/2004GL019970.
[3] Zhong S., Zhang N., Li Z-X., Roberts J.H. Supercontinent cycles, true polar wander, and very long-wavelength mantle convection // Earth and Planetary Sci. Lett. 2007. Vol. 261. P. 551-564.
[4] Трувицин В.П., Евсеев А.Н., Баранов А.А., Трубицын А.П. Структура конвекции при различной ширине зон фазовых переходов // Физика Земли. 2008. № 8. С. 3-14.
[5] Трувицин В.П. Сейсмическая томография и дрейф континентов // Там же. 2008. № 11. С. 3-19.
[6] Трувицин В.П. Уравнения тепловой конвекции для вязкой сжимаемой мантии Земли с фазовыми переходами // Там же. 2008. № 12. С. 83-91.
[7] Lassak T.M., MoNamara A.K., Edward J. et al. Core-mantle boundary topography as a possible constraint on lower mantle chemistry and dynamics // Earth and Planetary Sci. Lett. 2010. Vol. 289. P. 232-241.
[8] Nakagawa T., Taokley P.J. Effects of low-viscosity post-perovskite on thermochemical mantle convection in a 3-D spherical shell // Geophys. Res. Lett. 2011. Vol. 38 (L04309). doi:10.1029/2010GL046494.
[9] Chervov V.V., Chernykh G.G. Numerical modeling of 3d convection in the upper mantle of the Earth beneath lithosphere of Central Asia // J. of Eng. Thermophys. 2012. Vol. 21, No. 1. P. 78-89.
[10] Fleitout L., Yuen D.A. Steady state, secondary convection beneath lithospheric plates with temperature- and pressure-dependent viscosity //J. Geophys. Res. 1984. Vol. 89, No. B11. P. 9227-9244.
[11] Fleitout L., Yuen D.A. Secondary convection and the growth of the oceanic lithosphere // Phys. of the Earth and Planetary Interiors. 1984. Vol. 36. P. 181-212.
[12] Тычков С.А., Червов В.В., Черных Г.Г. Численная модель трёхмерной конвекции в верхней мантии Земли // Физика Земли. 2005. № 5. С. 48-64.
[13] Трубицын В.П. Основы тектоники плавающих континентов // Там же. 2000. № 9. С. 4-40.
[14] Кочин Н.Е. Векторное исчисление и начала тензорного исчисления. М.: Наука, 1965.
[15] Владимирова Н.Н., Кузнецов Б.Г., Яненко Н.Н. Численный расчёт симметричного обтекания пластинки плоским потоком несжимаемой жидкости // Некоторые вопросы прикладной и вычислительной математики. Новосибирск: ВЦ СО АН СССР, 1966. С. 186-192.
[16] Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука. Сиб. отд-ие, 1967. 195 c.
[17] Peyret R., Taylor T.D. Computational Methods for Fluid Flow. New York: SpringerVerlag, 1983.
[18] Самарский А.А. Теория разностных схем. М.: Наука, 1977.
[19] Червов В.В. Моделирование трёхмерной конвекции в мантии Земли с применением неявного метода слабой сжимаемости // Вычисл. технологии. 2009. Т. 14, № 3. С. 86-92.
[20] Busse F.H., Christensen U., Clever R. et al. 3D convection at infinite prandl number in cartesian geometry — a benchmark comparison // Geophys. Astrophys. Fluid Dynamics. 1993. Vol. 75. P. 39-59.
[21] Тычков С.А., Червов В.В., Черных Г.Г. О численном моделировании тепловой конвекции в мантии Земли // Докл. АН. 2005. Т. 402, № 2. С. 248-254.
[22] Тычков С.А., Червов В.В., Черных Г.Г. Трёхмерное моделирование конвекции под кратонами Центральной Азии // Вычисл. технологии. 2007. Т. 12. Спец. выпуск 4: Труды V Совещания российско-казахстанской рабочей группы по вычислительным и информационным технологиям. С. 85-95.
[23] Бушенкова Н.А. Неоднородности верхней мантии и современная структура литосферы центральной Сибири по данным сейсмотомографии на отражённых волнах. Автореф. дисс. ... к.г.-м.н. Новосибирск: Ин-т геологии СО РАН, 2004. 20 c.
[24] Добрецов Н.Л. Основы тектоники и геодинамики. Новосибирск: НГУ, ИГМ CO PAH, 2011. 491 c.
Поступила в 'редакцию 20 января 2014 г., с доработки — 4 апреля 2014 г.
Numerical modelling of 3D convection in the upper Earth mantle under lithosphere of Eurasia
Chervov Viktor V.1, Chernykh Gennadiy G.2'3, Bushenkova Natalia A.1, Koulakov Ivan Yu.1
The present study addresses convection beneath the continental lithosphere. The target area is bounded by a spherical region of the upper mantle (i.e., by the phase boundary, down to the depth of 670-700 km) from the zero meridian to eastern boundaries of Eurasia (145 °E) and from the equator to northern boundaries of the continent (80 °N). This large region includes several continental platforms as investigated objects, which are not homogeneous in their age and thickness. During the calculating for Eurasian continent, it becomes necessary to take into account the influence of sphericity. The area of our investigation includes the Siberian craton and its surrounding territory: West Siberian, East European and Chinese plates, Tarim, a complex of Tuva-Mongolian micro plates, Central Asia folded zone, and a segment of the Indian Plate. We constructed a three-dimensional numerical model of thermal convection based on Navier — Stokes equations in the geodynamic approach, in spherical variables using the artificial compressibility method. We were interested in determination of the influence of the craton geometry on the structure of mantle flows beneath the continental lithosphere. The calculation results illustrate the structure of convective flows. We also have constructed the numerical model for oceanic plate, which immersed under continental plate. The results of computations agree with geophysical understanding on the process of subduction.
Keywords: 3D numerical modelling, spherical coordinates, upper mantle, thermal convection, lithosphere, subduction, spreading.
Received 20 January 2014 Received in revised form 4 April 2014
1 A.A. Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Novosibirsk, Russia
2Institute of Computational Technologies SB RAS, 630090, Novosibirsk, Russia
3 Novosibirsk State University, 630090, Novosibirsk, Russia Corresponding author: Chernykh Gennadiy G. e-mail: [email protected] Aknowlegements: The work was supported by SB RAS (joint integration projects nos. 20, 76) and RFBR (initiative scientific project No. 13-05-00054).
