CC BY
11
GEODYNAMICS & TECTONOPHYSICS
Published by the Institute of the Earth's Crust, Siberian Branch, Russian Academy of Sciences
TECTONOPHYSICS
2021 VOLUME 12 ISSUE 1 PAGES 84-99 ISSN 2078-502X
DOI: 10.5800/GT-2021-12-1-0514
TECTONIC DEPRESSIONS ON THE EAST-EUROPEAN AND SIBERIAN PLATFORMS: NUMERICAL MODELING OF CONVECTION BENEATH THE EURASIAN CONTINENT
V.V. Chervov© \ N.A. Bushenkova© 13G.G. Chernykh © 23
1 Trofimuk Institute of Petroleum Geology and Geophysics, Siberian Branch of the Russian Academy of Sciences, 3 Academician Koptyug Ave, Novosibirsk 630090, Russia
2 Federal Research Center for Information and Computational Technologies, 6 Academician Lavrentiev Ave, Novosibirsk 630090, Russia
3 Novosibirsk State University, 1 Pirogov St, Novosibirsk 630090, Russia
ABSTRACT. In modern concepts, the upper mantle of the Earth is a highly viscous incompressible liquid, and its flow is described using the Navier - Stokes equations in the Oberbeck - Boussinesq and geodynamic approximations. Convec-tive flows in the upper mantle play a decisive role in the kinematics of lithospheric plates and the geological history of continental regions. Mathematical modeling is a basic method for studying convective processes in the mantle. Our paper presents a numerical model of convection, which is based on the implicit artificial compressibility method. This model is tested in detail by comparing our calculation results with the results of a well-known international test. It is demonstrated that the Fedorenko grids sequence method is highly efficient and reduces the computing time almost by a factor of eight. The numerical model is generalized in order to state the problem in a spherical system of coordinates. It is used to analyse the distribution of convective flows in the upper mantle underneath the Eurasian continent. The analysis shows that the thickness and geometrical parameters of the lithospheric blocks are the factors of significant influence on the distribution of convective flows in the upper mantle. The resulting structure of convective flows is manifested in the surface topography of large platform areas wherein the lithosphere thickness is increased. Thus, the locations of extended downward convection flows under the East European and Siberian platforms are clearly comparable to syneclises observed in the study area.
KEYWORDS: thermal convection; mathematical modeling; mantle; lithosphere thickness; craton; depression; syneclise; platform
FUNDING: The study received a partial financial support under Interdisciplinary Integration Project of SB RAS 44 and the State Assignment Project 0331-2019-0010.
RESEARCH ARTICLE Received: March 18, 2020
Revised: December 23, 2020
Correspondence: Natalia A. Bushenkova, [email protected] Accepted: January 11, 2021
FOR CITATION: Chervov V.V., Bushenkova N.A., Chernykh G.G., 2021. Tectonic depressions on the East-European and Siberian platforms: numerical modeling of convection beneath the Eurasian continent. Geodynamics & Tectonophysics 12 (1), 84-99. doi:10.5800/GT-2021-12-1-0514
ТЕКТОНИЧЕСКИЕ ПРОГИБЫ НА ВОСТОЧНО-ЕВРОПЕЙСКОЙ И СИБИРСКОЙ ПЛАТФОРМАХ: ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНВЕКЦИИ ПОД ЕВРАЗИЙСКИМ КОНТИНЕНТОМ
В.В. Червов1, Н.А. Бушенкова13, Г.Г. Черных2,3
1 Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, 630090 Новосибирск, пр-т Академика Коптюга, 3, Россия
2 Федеральный исследовательский центр информационных и вычислительных технологий, 630090 Новосибирск, пр-т Академика Лаврентьева, 6, Россия
3 Новосибирский государственный университет, 630090 Новосибирск, ул. Пирогова, 2, Россия
АННОТАЦИЯ. В соответствии с современными представлениями верхняя мантия Земли рассматривается как высоковязкая несжимаемая жидкость, для описания течения которой привлекаются уравнения Навье - Стокса в приближении Обербека - Буссинеска и геодинамическом приближении. Конвективные течения в верхней мантии Земли играют определяющую роль в кинематике литосферных плит и геологической истории развития континентальных областей. Основным методом исследования конвективных процессов в мантии Земли является математическое моделирование. В настоящей работе представлена численная модель конвекции, базирующаяся на неявной реализации метода искусственной сжимаемости. Приведены результаты детального тестирования модели путем сопоставления результатов расчетов с результатами известного международного теста. Продемонстрирована высокая эффективность метода последовательности сеток Р.П. Федоренко, позволившего сократить примерно в восемь раз время компьютерного счета. Представлено обобщение численной модели по постановке задачи в сферической системе координат. На основе построенной численной модели проанализировано распределение конвективных течений в верхней мантии Земли под Евразией. Показано, что мощность и геометрия блоков литосферы оказывают заметное влияние на распределение конвективных течений в верхней мантии Земли. Установившаяся структура этих течений проявляется в рельефе дневной поверхности обширных платформенных областей с увеличенной мощностью литосферы. Так, расположение протяженных нисходящих потоков конвекции под ВосточноЕвропейской и Сибирской платформами в плане явно сопоставимо с наблюдаемыми в рельефе синеклизами.
КЛЮЧЕВЫЕ СЛОВА: тепловая конвекция; математическое моделирование; мантия; мощность литосферы; кратон; прогиб; синеклиза; платформа
ФИНАНСИРОВАНИЕ: Работа выполнена при частичной финансовой поддержке Междисциплинарного интеграционного проекта СО РАН № 44 и Госзадания № 0331-2019-0010.
1. ВВЕДЕНИЕ
Исследование структуры и динамики мантии Земли является одной из основных задач геофизики. Важнейшая роль в исследовании верхнемантийной конвекции отводится математическому моделированию с применением различных аппроксимаций и получением численных моделей разной степени достоверности для различных масштабов и регионов [Blankenbach et al., 1989; Busse et al., 1993; Trubitsyn et al., 1994; Bobrov, Trubitsyn, 1995; Trubitsyn, Rykov, 1995; Chervov, 2002a, 2002b, 2009; Tychkov et al., 2005a, 2005b, 2005c, 2007; Chervov, Chernykh, 2014; Chervov et al., 2014; Heister et al., 2017].
Областью изучения в настоящем исследовании является обширная территория Евразии, включающая Восточно-Европейскую и Сибирскую платформы (рис. 1), в рельефе которых отчетливо выражены тектонические прогибы.
Наши предыдущие исследования влияния неоднородной структуры литосферы, а именно геометрии границ блоков с увеличенной мощностью, на структуру конвективных потоков в верхней мантии продемонстрировали заметную взаимосвязь этих структур и их
влияние на геодинамические характеристики сейсмоактивных районов с неутолщенной литосферой [Cher-vov, Chemykh, 2014; Chervov et al., 2014; Bushenkova et al.,
2016, 2018]. Логично предположить, что структура конвекции в верхней мантии оказывает определенное влияние и на геодинамические характеристики утолщенных блоков литосферы, если такой блок простирается над несколькими конвективными ячейками одновременно, как большие по площади Восточно-Европейская и Сибирская платформы. Представляется интересной природа возникновения Московского и Балтийского прогибов, прогибов платформенной части Северного Каспия и Волжско-Камского речного бассейна (см., например [Lobkovsky et al., 2004] и ссылки в ней), а также наиболее значительно проявленного на Сибирской платформе Вилюйского прогиба, ширина которого достигает 600 км. В работах О.П. Полянского с соавторами предлагается рифтовый механизм становления Вилюйского прогиба. Приводятся убедительные доводы, подкрепленные численными экспериментами, указывающие на природу возникновения депрессии в результате растяжения земной коры [Polyansky et al., 2013,
2017, 2018].
В работе Е.В. Артюшкова объяснения появления прогибов сводятся в основном к истории осадконакоп-ления в них [Artyushkov, 1993], но глубинные процессы и условия, способствующие их появлению, до сих пор не изучены.
В настоящей работе авторами представлена численная модель, разработанная с применением неявной реализации метода искусственной сжимаемости [Yanenko, 1971; Peyret, Taylor, 1983], приведены результаты ее детального тестирования путем сопоставления с данными международного теста [Busse et al., 1993]. В качестве приложения численной модели рассмотрена задача о конвекции под Евразийским континентом.
В результате численного моделирования удается выявить вероятную причину возникновения прогибов, наблюдаемых на Сибирской и Восточно-Европейской платформах. Данная работа продолжает проведенные ранее исследования [Chervov, 2009; Chervov, Chernykh, 2014; Chervov et al., 2014].
2. МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ В ДЕКАРТОВЫХ КООРДИНАТАХ
Для математического моделирования конвекции в верхней мантии Земли привлекается известная математическая модель, включающая в себя обезразмеренные уравнения Навье - Стокса в приближении Обербека -
Буссинеска и геодинамическом приближении (см. например [Dobretsov et al., 2001]):
du + dv + dw = 0, (1)
dx dy dz
Vp = V-s + (RaT )e, (2)
dT 2
-+ V- (VT) = V2T, (3)
dt
где V=(u, v, w) - вектор скорости, u, v, w- его компоненты, p - дефект давления [Landau, Lifshits, 1986], T - температура, о - тензор вязких напряжений, e=(0,0,1) - безразмерный единичный вектор, t - время, Ra - число Рэлея:
Ra = а p0d3AT / Vox. (4)
В (4) величина d - расстояние от поверхности Земли до подошвы верхней мантии, c - коэффициент температуропроводности, п0 и р0 - характерные значения вязкости и плотности соответственно, AT - разность между температурой на подошве верхней мантии и температурой на поверхности, a - температурный коэффициент объемного расширения мантийной жидкости, gz - z-компонента вектора силы тяжести.
Используется стандартная декартова система координат, в которой ось z направлена вертикально вверх, против силы тяжести (рис. 2). Компоненты тензора
75°
65°
55°
45°
35°
25°
м 8000
6000
3000
2000
1000
500
300
100
50
0
50
100 -300 -500 -1000
2000 -4000 -6000 -8000
130° в.д.
Рис. 1. Рельеф Евразии. Сплошными линиями выделены границы блоков с существенно различной мощностью литосферы. Цифрами обозначены: 1 - Тувино-Монгольский микроконтинент, 2 и 3 - Северо- и Южно-Китайский кратоны соответственно, 4 - Тарим, 5 - Индийский кратон, 6 - Аравийская платформа.
Fig. 1. Relief of Eurasia. Solid lines - boundaries of blocks that significantly differ in the lithospheric thickness.
Numbers: 1 - Tuva-Mongolia microcontinent, 2 and 3 - North and South China cratons, respectively, 4 - Tarim, 5 - Indian craton, 6 -
Arabian platform.
v _ du _ dw _ dT = 0 dy dy dy
и _—_ — _ dT _ о
dx dx dx
Рис. 2. Условия на гранях параллелепипеда для температуры и вектора скорости в декартовых координатах. Fig. 2. Conditions on the faces of the parallelepiped for temperatures and velocity vectors in the Cartesian coordinates.
вязких напряжении о в декартовых координатах имеют следующий вид [Landau, Lifshits, 1986]:
du dx
du + dv dy dx ^
V
du dy
dv dx
V
dv
2 dy; a y
О = n
zx /
du dz
dw
dx,
dv dz
dw dy
du + dw dz dx
dv + dw dz dy ^
~ dw
; 0zz = .
dz
Здесь вязкость жидкости п - нелинейная функция температуры, которая будет задана ниже.
Уравнения (1) - (3) были получены в результате применения обезразмеривающих соотношений (5):
(x,y,z) = d•(x',y',z'), t = —t', T = AT-T',
n = ПоП > p =
n* p, v = XV'. d2 d
(5)
Вопрос о применимости приближения Обербека -Буссинеска возникает в связи с глобальными масштабами течения, существенным варьированием температуры и огромными значениями вязкости. Теоретические оценки применимости приближения даны В.В. Пухначе-вым [Andreev et al., 1998]. Согласно этим оценкам, приближение применимо, если параметр N = gz d3 / п0Х> 1. Простые оценки показывают, что для задач конвекции в верхней мантии Земли N ^100.
Первоначально для нашего исследования рассматривалась тестовая модельная задача [Busse et al., 1993], в которой решение определялось в декартовой прямоугольной системе координат в прямоугольном параллелепипеде (рис. 2): П=[0Д]х[0,У]х[0,2Г].
Уравнение движения (2) в декартовых координатах может быть записано следующим образом:
да„ da„„ ^
(6)
dx dy
dz
dp dx
9(Jyx , dSyy , dSyz
dx
dy dz dy
^ - Ra T.
dx dy dz dz Ставились следующие граничные условия: на поверхностях x=x^0, x=x2=X; (y<y<y2, z1<z<z2):
dv dw
u = — =-= 0,
dx dx
на поверхностях y=y1=0, y=y2=^; (x1<x<x2, z1<z<z2):
du dw
v = — =-= 0,
dy dy
(7)
(8)
(9)
(10)
на поверхностях z=z1=0; z=z2=Z; (x1<x<x2, y <y<y^J:
u = v = w = 0. (11)
Для температуры, как и в [Busse et al., 1993], на боковых гранях ставятся условия теплоизоляции (адиабатическая стенка), т.е. первые производные по нормалям на вертикальных стенках равны нулю. На верхней и нижней гранях ставятся условия Дирихле: нулевая температура на верхней и некоторая фиксированная температура (в обезразмеренной постановке равная единице) на нижней. Система дифференциальных уравнений (1) - (3) требует задания начального распределения температуры T( x, y, z ,0) = T (x, y, z ) Компоненты скорости и давление при t=0 не задаются ввиду стационарности уравнений (1), (2).
3. ЧИСЛЕННАЯ МОДЕЛЬ ТРЕХМЕРНЫХ КОНВЕКТИВНЫХ ТЕЧЕНИЙ В ВЕРХНЕЙ МАНТИИ ЗЕМЛИ
Одним из известных подходов к решению трехмерных задач гидродинамики несжимаемых жидкостей в естественных переменных является метод искусственной сжимаемости [Vladimirova et al., 1966; Yanenko, 1971; Peyret, Taylor, 1983; Rizzi, Eriksson, 1985]. Численное моделирование трехмерных конвективных течений в верхней
мантии Земли, основанное на решении системы дифференциальных уравнений (1) - (3) с применением неявной реализации метода искусственной сжимаемости [СИегуоу, 2009] и метода дробных шагов [Уапепко, 1971], на каждом слое по времени осуществлялось согласно уравнениям:
Р
m,n+1
-pmn + Дт • c2Vh • V"
(12)
Дт
m,n+1
+
Ox,
dv. dv, —L +---
dx. Ox,
At
3xx ^
m,n+1
+
Ra •Tm • e,, (13)
-VlhTn
(14)
где т - номер шага по времени % п - номер шага по фиктивному времени т; с2 - релаксационный параметр;
- разностный оператор градиента; Дт, Дt - шаги сетки по фиктивному (т.е. использовались итерации) и физическому времени; ¡, к=1, 2, 3; по повторяющимся индексам производится суммирование.
Дискретные уравнения (12) - (14) записаны в символическом виде. Для решения задачи вводилась равномерная в каждом направлении сетка с шагами по пространственным переменным:
К = Х+1 - х-' К = У+1 - Уj > К = - гк. (15)
Ограничения на с и Дт
-тп следующие [Taru-
nin, 1990]:
Ат < min(hx,hy,hz)/ c.
Численный алгоритм сводится к выполнению следующих этапов.
1. В исследуемой области задается начальное распределение температуры (т=1, t=0), удовлетворяющее граничным условиям. Исходя из полученного начального распределения температуры вычисляется поле вязкости (формула приведена ниже). Дефект давления и компоненты вектора скорости при t=0 полагаются нулевыми.
2. Из векторного разностного уравнения (13) находится сеточное поле скорости Vm,n+1 с привлечением уравнения (12), из которого дефект давления на верхнем слое по фиктивному времени выражается через дивергенцию поля скорости и подставляется в (13); поле температуры на временном слое т вычисляется путем решения (14).
3. Итерационный процесс повторяется до некоторого значения tm =(т — 1) Д t.
Решения получающихся многомерных разностных уравнений (13) находятся с применением итерационной схемы стабилизирующей поправки [Уапепко, 1971] с итерационным параметром Дт = тп+1 — тп. При решении многомерного разностного уравнения (14), аппроксимирующего дифференциальное уравнение температуропроводности (3), в качестве схемы интегрирования применялась схема стабилизирующей поправки.
Как уже отмечалось выше, тестирование численной модели осуществлялось путем решения модельной задачи [Busse et al., 1993] в единичном кубе. Табл. 1 показывает величины, используемые в [Busse et al., 1993].
Интегралы в табл. 1 вычислялись с применением квадратурной формулы трапеций. Размерные значения (в системе СИ), которые были использованы в [Busse et al., 1993] и в настоящей работе:
d=2700000 м, Д7=3700 °C, х=10-6 м2/с, а =105 °C-1, р0=3300 кг/м3,g=10 м/с2, ^°=1.2065-1024 кг/(м-с), Ra=20000. В качестве начального выбиралось следующее распределение температуры [Busse et al., 1993]:
T(x,y,z,tQ J = T0[x, y,z) =
= (1 - z) + 0.2(cos(^x / X) + cos(^y / Y ))sin(^z).
Вязкость мантийного вещества задавалась по формуле:
П (T ) = exp( / (T + 6)- в / (0.50 + 6)),
в = 225/1n(q)-0.251n(q), (16)
0 = 15/ln(q)-0.50, q = T=0/v[
20.
Результаты расчетов авторов (Che) сопоставляются в табл. 2 с расчетными данными Кристенсена (Chr), как наиболее полными из имеющихся в статье [Busse et al., 1993]. Кристенсен (Christensen U.) при моделировании применял гибридный спектрально-конечно-разностный метод и скалярный потенциал для описания полученных течений. Данные, приведенные во втором столбце табл. 2, были получены им при помощи экстраполяции и перенесены на сетку с числом узлов (32x32x64).
Ошибка в настоящих расчетах определялась по формуле:
Che - Chr
Err =
Chr
•100 %.
Результаты расчетов в настоящей работе получены с применением хорошо известного метода последовательности сеток [Fedorenko, 1964] (успешное применение этого метода можно найти в работах [Tarunin, 1975; Tychkov et al., 2005a, 2005b, 2005c; Chervov, 2002a, 2002b, 2009; Chervov et al., 2014], и др.) и на сетке с числом ячеек 128x128x128 практически совпадают с численным решением [Busse et al., 1993]. Вычисления проводились с шагом по времени At, зависящим от пространственных характеристик численного алгоритма, а именно, исходя из соотношения At/h2=const. Здесь h=hx=hy=hz - обезразмеренный шаг по пространству. На сетке (32x32x32) величина обезразмеренного шага по времени At=0.2, и до получения стационарного решения понадобилось 240 слоев по времени.
Решение на сетке (32x32x32) было проинтерполи-ровано на сетку (64x64x64), и понадобилось 400 слоев с временным шагом At=0.2/4=0.05 для получения решения, близкого к стационарному. Далее осуществлялся переход на сетку с числом ячеек (128x128x128);
д
m
П
величина шага по времени полагалась равной At=0.05/4= =0.0125; выполнялось 780 слоев по времени.
Время расчетов с применением метода последовательности сеток: 16 с, 113 с, 1627 с. При прямом вычислении на сетке (128x128x128) потребовалось 3 часа 54 минуты на ПК с процессором AMD Ryzen 5 1500Х, т.е. примерно в восемь раз больше.
Табл. 2 иллюстрирует сопоставление полученных в настоящей работе результатов расчетов на последовательности сеток с результатами Кристенсена. Результаты в табл. 2 свидетельствуют о сходимости последовательности сеточных решений. Можно видеть также, что результаты наших тестовых расчетов достаточно близки к результатам [Busse et al., 1993], что
Таблица 1. Перечень величин, по которым осуществлялось сопоставление Table 1. Comparison list
I число Нуссельта по формуле [Blankenbach et al., 1989] где Nu = — (X Y )—1 ff T dxdy, T = dL. - верхняя поверхность параллелепипеда V dz
Л среднеквадратичная скорость Vrms = ^ где A - замкнутая область - параллеле х.1.zШи + v2 + w2]dxdydz A пипед со сторонами X=1, Y=1 и Z= 1
III значение вертикальной компоненты скорости w и температуры T в угловых точках среднего сечения конвективного слоя
IV dT значение теплового потока Q =--в dz угловых точках верхней поверхности куба
V г дТ интегральный параметр, вычисляемый по формуле X(х,у)= 1 —ёу вдоль линии, параллельной оси У, ■{ д начинающейся точками (0, (%, (1, фронтальной плоскости Х)
VI средняя температура Тт = Tdxdy, вычисляемая на горизонтальных сечениях области 5 у и 5 у , на глубинах и х=У
VII значение компоненты вектора завихренности ы" в точке (%, %)
Таблица 2. Сопоставление результатов расчета авторов (Che) на последовательности сеток: n32~(32x32x32); n64~(64x64x64); n128~(128x128x128) с результатами [Busse et al., 1993] (Chr - результаты, полученные Кристенсеном на сетке (32x32x64) Table 2. Comparison of the authors' calculations (Che) in grid sequences n32~(32x32x32), n64~(64x64x64), and n128~(128x128x128) with the data from [Busse et al., 1993] (Chr - Christensen' results in grid 32x32x64)
Имена Chr Che Относительная ошибка в процентах
n32 n64 n128 n32 n64 n128
Nu 3.0393 3.0473 3.0431 3.0393 0.2632 0.1250 0.0000
Vrms 35.1320 35.2439 35.1546 35.1280 0.3185 0.0643 0.0114
w(0,0,%) 165.9100 166.2674 166.0415 165.8993 0.2154 0.0793 0.0064
w(0,Y,y2) -26.7200 -27.1400 -26.8274 -26.7199 1.5719 0.4019 0.0004
w(X,Y,i/2) -58.2300 -59.0922 -58.4458 -58.2301 1.4807 0.3706 0.0002
T(0,0,i/2) 0.9053 0.9073 0.9058 0.9053 0.2209 0.0552 0.0000
T(0,Y,/) 0.4957 0.5007 0.4969 0.4956 1.0087 0.2421 0.0202
T(X,Y,/) 0.2393 0.2438 0.2404 0.2392 1.8805 0.4597 0.0418
Q[0,0) 5.8339 5.8395 5.8356 5.8339 0.0960 0.0291 0.0000
Q[0,Y) 1.7136 1.7313 1.7179 1.7136 1.0329 0.2509 0.0000
Q[X,Y) 0.7684 0.7752 0.7700 0.7684 0.8850 0.2082 0.0000
l(0,/) -0.5059 -0.4954 -0.5032 -0.5056 2.0755 0.5337 0.0593
l(X/2,/4) -0.1921 -0.1986 -0.1937 -0.1922 3.3837 0.8329 0.0521
l[X,/) -0.1388 -0.1400 -0.1391 -0.1388 0.8646 0.2161 0.0000
wz(%,%,%) -11.1250 -10.9746 -11.0903 -11.1251 1.3519 0.3119 0.0009
Время счета 16 с 113 с 1627 с
свидетельствует о высокой эффективности численной модели. Достаточно хорошее совпадение получено уже на сетке (32x32x32).
Следует отметить также, что вычисления для наших тестовых исследований проводились на ПК и использовался существенно отличающийся от примененного в работе [Busse et al., 1993] алгоритм, где вычисления проводились на суперкомпьютере. Численные эксперименты по оптимизации шага по времени At и величины шага по фиктивному времени Ат в настоящей работе не выполнялись. Существенное повышение эффективности численной модели может быть получено с применением экстраполяции по Ричардсону [Marchuk, Shaidurov, 1979].
Авторы разработали три подхода к решению задач конвекции в верхней мантии Земли, основанные на применении неявных методов расщепления по пространственным переменным. Это метод с применением векторного потенциала и завихренности [Chervov, 2002а]; метод, основанный на неявных методах расщепления по пространственным переменным с коррекцией давления [Chervov, 2006], и, наконец, изложенный выше подход, основанный на неявной реализации метода искусственной сжимаемости. По мнению авторов, метод искусственной сжимаемости в рассмотренных задачах наиболее прост в реализации по сравнению с рассмотренными. Возможности метода иллюстрируются сопоставлениями рассчитанных на его основе параметров решения модельной задачи с результатами международного теста. Как уже отмечалось выше, результаты расчетов, приведенные в табл. 2, демонстрируют достаточно высокую эффективность разработанной численной модели.
Изложенный выше алгоритм численного моделирования и его модификация в сферических переменных были успешно применены при расчетах конвективных
Рис. 3. Схематическое изображение связи сферических и декартовых прямоугольных координат. Fig. 3. Schematic relationship between the spherical and Cartesian rectangular coordinates.
процессов в верхней мантии Земли под континентальной литосферой. Вначале это была значительная часть Азии, включающая в себя Сибирский кратон и ЗападноСибирскую платформу, Тарим и Северо-Китайский кратон. Впоследствии область была расширена и стала вмещать практически всю Евразию [СИегуоу СИегпукИ, 2014]. Аналогичный подход применялся также при решении задачи конвекции в зонах спрединга и субдукции [СИег-уоу, СИегпукИ, 2019].
4. МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ В СФЕРИЧЕСКИХ КООРДИНАТАХ
Известно, что декартовы прямоугольные и сферические координаты связаны соотношениями: x=r■smв■шsф, у=г^пв^тф, z=r■cosв. Здесь полярный угол в отсчиты-вается от положительного направления оси 2, а азимутальный угол ф в плоскости ху - от положительного направления оси х (рис. 3).
Из системы уравнений (1) - (3) в результате привлечения сферических переменных могут быть получены известные уравнения (17) - (21):
Г2 sin в
dV f d(sin evв) d(r V r--+ r—-- + sin в-
dp
др
2i]r
dV p
dp
дв dr
+ cos 6V6 + sin в Vr
= 0; (17)
+ дв
idV p i
nr sin e de sin e
+— dr
Пr3 sin2(
dVp +1 dr r
or
дф
1 dVr
- cose ■Vf
sin в dp
-V p
2 ■ a dp
= r Sinfl
+ (18)
dp
dp'
dp
П r
sin в
dVв . ndVp
--+ sin в--cos eVp
dp дв
д +—
2 n r sin в
dVв
+ Vr
+
+— dr
+2 n r cos в
■3 ■ л dvв 1 i— -v в \
Пr sin в —
dr r . дв
(19)
dVe d(rV
+
дв dr
d_ dp
sin в
dVr dp
+ sin в
dV p
r —--V p
dr
д + de
dVr \r »Vi-V e A
n sin e --h
de dr )
+
-¡dp
+
+
(20)
+— dr
dT — +
2 n r2 sin в 1
dVr
dr
dVr 2
+ 2n r sin в-= r sin в
dr
- Ra T
dr
dt r2 sin в 1
r—(vVT) + r—(sin eVeT) + sin в—(r2VrT)
дру ' двк ' dr '
r sin в
1 d2T д
sin в др дв
■ адТ sin в—
дв
д
+ Sin в— dr
dT_ dr
, (21)
где Vр, V , Vr- компоненты вектора скорости V.
П
Математическая модель обезразмерена следующим образом:
R
r = R■ r', t =—t', T = AT-T', n =
X
p = p, V = XV'. R2 R
Число Рэлея в приведенных уравнениях:
Ra = " gr Por3at / ПоХ,
(22)
(23)
где й - радиус Земли, gr - г-компонента вектора силы тяжести g, которая ориентирована вдоль оси г (рис. 3) в направлении к центру планеты.
Основные вычислительные эксперименты проводились в сферической расчетной области
П = [^1^2 ]х[01, в2 ]х[г1,г2 ]
(^=5700, г2=й=6370 км; <^=0°, ф2=145° в.д. и ^=20°, 02=8О° с.ш.), охватывающей большую часть верхней мантии под Евразийским континентом (рис. 4).
Для вектора скорости на боковых границах задавались условия проскальзывания, на нижней и верхней гранях - условия прилипания. Для температуры на боковых гранях ставились условия теплоизоляции (адиабатическая стенка). На верхней и нижней гранях - условия Дирихле: нулевая температура на верхней и 1800 °С (в обезразмеренных переменных равная единице) на нижней:
на поверхностях ф=ф1, Ф=Ф2; (в1<в<в2, г1<г<г2):
V *=dVL=dXL=дТ=0 д* д* д*
на поверхностях в=в1, в=в2, (ф1<ф<ф2, r1<r<r2):
dV*
de
V * ctge = V9 = dV! = dL=0,
g de de
(24)
(25)
на поверхности r=r1, (/1<j'</2, q1<q<q2):
V r = Vû = Vr = 0, T = 1, на поверхности r=r2, (/1</</2, q1<q<q2):
Vv = Vв = Vr = T = 0.
(26)
В начальный момент времени при t=t0 для температуры задавались начальные условия:
Т(<р,9, г, С0) = Тв(ф, 9, г). (27)
Компоненты вектора скорости полагались нулевыми. При выборе основных параметров задачи в системе СИ, пригодных для верхней мантии (23), число Рэлея равно 4х105:
й=6370000 м, ¿=670000 м, Д7=1800 °С, X =10-6 м2/с, а =2-10-5 °С, р0=3300 кг/м3, (28)
gz=10 м/с2, п0=1.01871-1021 кг/(м-с). Обезразмеренная вязкость мантийного вещества задавалась, как и в работе [Tychkov et а1., 2005а]:
П (р, в ,г ,t) = ехр(г - аТ (,в, г ,ь )), (29)
где Ь=3.47; а=5.34.
Аналогично модельной задаче в декартовых координатах при численном решении за основу выбран метод искусственной сжимаемости.
Уравнения (17) - (20) с граничными условиями (23) -(25) после дискретизации проинтегрированы на каждом временном слое при помощи метода стабилизирующей поправки, где на дробных шагах применены скалярные трехточечные прогонки. В вычислительной схеме полагалось с2=2500, Дт = тт(Д<^ ,Дв ,Аг)/ с.
Тестирование численной модели в сферической системе координат осуществлялось путем ее сопоставления с результатами моделирования в декартовых координатах для близких по геометрическим размерам областей и идентичных начально-краевых условий.
Ф.А/,
v = V = V = Т = 0
V „ = дТ = дт = дт_ = 0 5ф 5ф 5ф
дТ_ V, 0 = V о = дТ = дТ =0
50 V g 0 V д0 д0 0
Рис. 4. Граничные условия в сферической области: на вертикальных гранях - условия проскальзывания, на горизонтальных плоскостях - условия прилипания.
Fig. 4. Boundary conditions in a spherical area. Vertical - slip conditions; horizontal - stick conditions.
r
r
Численная модель в сферических координатах реализована посредством комплекса программ. Расчеты в настоящей работе выполнены на основе этого комплекса программ, содержащего код для вычислений в МР1-технологии (на многопроцессорной ЭВМ) с применением метода последовательности сеток [СИегуоу 2018].
5. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТРЕХМЕРНЫХ КОНВЕКТИВНЫХ ТЕЧЕНИЙ В ВЕРХНЕЙ МАНТИИ ЗЕМЛИ ПОД ЕВРАЗИЕЙ
В литосферную плиту континента, модельная мощность которой задана значением 120 км, включены древние блоки с увеличенной мощностью: Восточно-Европейская (на западе) и Сибирская (на востоке) платформы, Тувино-Монгольский комплекс (микроконтинент), находящийся к юго-западу от Сибирского кратона, Тарим, Северо- и Южно-Китайский кратоны в юго-восточной области Евразии, Индийская плита на юге, Аравийская - в юго-западной части материка. Мощность утолщенных элементов литосферы в расчетах составляет 220 км, за исключением Аравийской плиты (в модели -180 км). В центральной части Евразии расположен Центрально-Азиатский складчатый пояс, мощность литосферы которого в нашей модели принималась равной 80 км, так как согласно геолого-геофизическим наблюдениям и оценкам [ВшИепкоуа, 2004; ВшИепкоуа et а1., 2008] мощность литосферы в его пределах составляет 40-95 км.
Расположение и конфигурация модельных элементов литосферы континента приближенно соответствуют общепринятым геолого-геофизическим данным,
плоская проекция представлена на рис. 1. Погруженная часть каждого заданного кратона представляет собой перевернутую усеченную пирамиду, верхнее основание которой и изображено на схеме (см. рис. 1).
К востоку от Сибирского кратона (Предверхоянье), согласно геолого-геофизическим данным, мощность литосферы снижается до 80 км (ловушка). Это изменение мощности обусловлено тем, что под складчатыми поясами (Верхоянский хребет) мощность литосферы значительно ниже средней мощности литосферы для континентов (~120 км).
Мощность литосферы складчатых поясов может варьироваться в пределах от 40 до 100 км. Кроме того, на основе анализа доступных для изучаемой территории геолого-геофизических данных, в том числе сей-смотомографических (см. []акоу1еу et а1., 2012; и др.]), можно заключить, что мощность литосферы восточнее Сибирского кратона в районе сочленения Евразийской и Северо-Американской плит (вдоль Верхоянского хребта) составляет не более 80 км. Аналогично так называемыми «ловушками» [Tychkov et а1., 2005с] заданы мощности литосферы других складчатых областей в расчетной области.
В нашей постановке глубина расчетной области составляет 670 км, что соответствует общепринятому значению глубины залегания границы верхней и нижней мантии.
Численные эксперименты проводились в сферической области 0°<ф<145°, 20°<0<80°, 5700 км<г<6370 км, на сетке (806x334x135) точек. Шаг по времени варьировался от 1 до 2.5 млн лет. Результаты расчетов (рис. 5,
S 0 -200
I -400
о
Е? -600
« -200 J -400
-600
5 0 «- -200 J -400 -600
20°
24°
28°
32°
36°
40°
44°
48°
52°
56°
60°
64°
68°
72°
76°
80° с.ш.
0
300
600
900
1200
1500
1800°C
Рис. 5. Вертикальные сечения поля температуры через Аравийскую и Восточно-Европейскую платформы по долготе 40°(а), 45°(б) и 50° (е); красной линией показан схематичный рельеф подошвы литосферы; шкала температуры актуальна и для рис. 6, 7, 8, 9.
Fig. 5. Vertical sections of the temperature field across the Arabian and the East European platforms along 40° (a), 45° (б) and 50° (е) longitudes. Red line - schematic relief of the lower lithospheric boundary. The same temperature scale in Figs. 6, 7, 8, 9.
0
6, 7, 8, 9) соответствуют значению модельного времени t=400 млн лет.
На верхней границе расчетной области была задана неоднородная мощность литосферы по блокам, соответствующим геологическим структурам разного возраста и генеза, наибольшие значения мощности были заданы для кратонов (до 220 км).
Уточнение конфигурации и неоднородности мощности кратонов по всем доступным геолого-геофизическим данным, включая сейсмотомографические, позволило получить при числе Рэлея, равном 4х105, следующие распределения температурного поля (рис. 5, 6, 7, 8, 9). Расчеты выполнены на достаточно мелкой сетке -20x20x5 км (20 км вдоль экватора и меридианов, 5 км по глубине, вдоль радиуса Земли, т.е. 0.18°х0.18°х5 км). Относительные перемещения структурных блоков, составляющих Евразийский континент, в модели не учитывались.
При сопоставлении полученной в настоящей работе численной трехмерной модели тепловой конвекции
в верхней мантии под Евразией с осредненным рельефом дневной поверхности Восточно-Европейской и Сибирской платформ выявлены следующие характерные особенности.
1. Прогиб на Восточно-Европейской платформе, который пересекает платформу с юга на север (Прикаспийская низменность, долины рек Волги и Камы), в плане явно согласуется с расположением нисходящего потока тепловой конвекции в верхней мантии.
2. Пересечение нисходящих конвективных потоков (в точке 43° в.д, 60° с.ш. на рис. 1 и 7) попадает в центр Московской синеклизы (см. рис. 5, б; 10).
3. Многие возвышенности Восточно-Европейской платформы в плане совпадают с полученными в расчетах восходящими потоками конвекции.
При сопоставлении депрессии, протяженной от 35° в.д. на юге Восточно-Европейской платформы в направлении на северо-северо-запад (см. рис. 1; рис. 10), и нисходящего теплового потока в рассчитанной структуре верхнемантийной конвекции (см. рис. 5, 7, 8, 9)
^-400 Я -600
20°
24°
28°
32°
36°
40°
44°
48°
52°
56°
60°
64°
68°
72°
76°
80° с.ш.
0
300
600
900
1200
1500
1800 °C
Рис. 6. Вертикальные сечения поля температуры через Северо-Китайский кратон и восточную часть Сибирской платформы по долготе 120° (а), 122° (б), 125° (е) 127° (г) и 130° (д).
Fig. 6. Vertical sections of the temperature field across the North Chinese craton and the eastern Siberian platform along 120° (a), 122° (б), 125° (е), 127° (г), 130° (д) longitudes.
850 950 1050 1150 1250 1350 1450 1550 °С
Рис. 7. Сечение температурного поля на глубине 220 км (красными линиями показаны контуры кратонов, голубыми -ловушек).
Fig. 7. Horizontal section of the temperature field at a depth of 220 km. Lines: red - contours of cratons with the thick lithosphere; blue - areas with the thin lithosphere).
75°-с.ш.
65°-'
55°-
45°-
35°
25°-
1 10
800 900 1000 1100 1200 1300 1400 1500 1600 °С Рис. 8. Сечение температурного поля на глубине 270 км. Fig. 8. Horizontal section of the temperature field at a depth of 270 km.
75°-с.ш.
65°-55°-45°-35°-25°-
Рис. 9. Сечение температурного поля на глубине 320 км.
Fig. 9. Horizontal section of the temperature field at a depth of 320 km.
800 900 1000 1100 1200 1300 1400 1500 1600 °С
прослеживается их явное пространственное соответствие.
Так, границы нисходящих потоков очерчивают и границы областей повышенного рельефа на ВосточноЕвропейской плите, а восходящие потоки конвекции достаточно точно «попадают» в центральные части видимых на карте возвышений рельефа (см. рис. 1, 6, 7, 8, 9), таких как Белорусская, Воронежская и Волго-Уральская антеклизы. В непосредственной близости от нисходящего потока находится и Балтийская синеклиза.
На рис. 5 изображены три сечения, пересекающие Восточно-Европейской плиту: (а) западнее долины реки Волги, (б) восточнее реки и (в) приблизительно вдоль самой реки. Под Восточно-Европейским кратоном на сечении (б) видно, что интенсивность восходящего потока здесь в среднем заметно ниже, чем западнее (а) и восточнее (в). Вероятно, именно это и определяет прогиб в данном районе платформы.
На рис. 6 вертикальные сечения пересекают Севе-ро-Китайский и Сибирский кратоны. Анализируя нисходящие потоки в широтах от 66.7° (рис. 6, а) до 63.2° (рис. 6, д), можно проследить депрессию, отображающую Вилюйскую синеклизу (рис. 11). Нисходящий поток, расположенный северо-восточнее (вдоль р. Лены), сопоставим с Предверхоянским (Приверхоянским) передовым прогибом.
Прогиб отделяет Сибирскую платформу от Верхо-яно-Колымской области (часть Верхояно-Чукотской складчатой области), мощность литосферы которой в
нашей модели задана ловушкой мощностью 60-80 км (на рис. 7, 8, 9 выделена голубым контуром).
Логично предположить, что именно благодаря существующему перепаду мощностей соседствующих ли-тосферных блоков и их геометрии, также и в конвекции под Сибирской платформой образовалась структура, схожая в плане с Вилюйской синеклизой.
Существующие незначительные несоответствия рельефа (см. рис. 1, 10, 11) и сечений теплового поля (см. рис. 5, 6, 7, 8, 9) можно объяснить как недостаточно полными данными о глубинах подошвы литосферы платформы и окружающих территорий, так и факторами геологической природы.
На юго-востоке от Вилюйской синеклизы (59.8°) восходящий поток определяет поднятие, которое хорошо прослеживается на горизонтальных сечениях (см. рис. 7, 8, 9), выделяется на карте рельефа (см. рис. 1) и соответствует Алданскому щиту на схеме кратона (рис. 11). На севере, около 70-й широты (см. рис. 6, 7), восходящий конвективный поток может быть взаимосвязан с Анабарским выступом (антеклиза).
Значительная часть поднятий и прогибов Сибирской платформы (Тунгусская синеклиза, Байкитская антеклиза и др.) перекрыта мощными траппами, возникшими 230-250 млн лет назад. Образованию этих магматических излияний на поверхность плиты способствовало обилие под платформой интенсивных восходящих потоков мантийной конвекции в совокупности с вероятным прохождением плиты в ту эпоху над
ВОСТОЧНО-ЕВРОПЕЙСКАЯ ПЛАТФОРМА ТЕКТОНИЧЕСКАЯ СХЕМА
Выступы (выходы) архейского фундамента (>2500 млн пет)
Выступы (выходы) карельского фундамента (>1600 млн лет)
Эпнкарепьский чехол
[C-CC^j Складчатое сооружение Донбасса
A v \ Синеклизы (БС-Балтийская,, МС— у а у —Московская, ПС-Прикаспийская)
Антеклиэы (БА-Белорусская, ВА -Воронежская. ВУА —Волго-Уральская)
Авлакогены (ВК-ВерхнекамскНй, ДД-Днепровско -Донецкий, КС— Каэанско-Сергиевский, М—Московский, Пч-Пачелмский, СР-Среднерусский)
Границы платформы
Области с байкальским складчатым фундаментом
Границы территорий с байкальским складчатым фундаментом
Рис. 10. Тектоническая схема Восточно-Европейской платформы [Shatsky, 1946]. Fig. 10. Tectonic scheme of the East European platform [Shatsky, 1946].
Рис. 11. Схема районирования Сибирской платформы [Bogoyavlenskaya et al., 1991].
1 - щиты и поднятия: А - Анабарский, А1 - Оленекское, Б - Алданский, Б1 - Становая область; 2 - погружения склонов щитов и древних складчатых областей: В - Среднеленская мо-ноклиза, В1 - Юдомо-Майская впадина, Г - Прибайкальская и Д - Приенисейская моноклизы; 3 - синеклизы: Е - Ангаро-Тасеевская, Ж - Вилюйская и З - Тунгусская; 4 - прогибы: И -Приверхоянский, К - Лено-Анабарский, Л - Хатанго-Пясин-ский (Енисей-Хатангский), М - Ангаро-Вилюйский, М1 - При-саянский; 5 - авлакогены (цифры в кружках): 1 - Уджинский,
2 - Вилюйский, 3 - Иркинеевский; 6 - впадины: а - Чульман-ская, б - Токинская, в - Рыбинская, г - Иркутская; 7 - Попи-гайская астроблема; С.О. - складчатая область.
Fig. 11. Zoning scheme of the Siberian platform [Bogoyavlenskaya et al., 1991].
1 - shields and uplifts: A - Anabar, A1 - Olenek, Б - Aldan, Б1 -Stanovoy area; 2 - slopes of shields and ancient orogens: В - Middle Lena monoclise, B1 - Udoma-Maysk depression, Г - Pre-Baikal monoclise, Д - Pre-Yenisei monoclise; 3 - syneclises: E - Angara-Taseeva, Ж - Vilyui, З - Tunguska; 4 - depressions: И - Pre-Verkhoyansk, K - Lena-Anabar, Л - Khatanga-Pasin (Yenisei-Khatanga), M - Angara-Vilyui, M1 - Pre-Sayan; 5 - aulacogens (numbers in circles): 1 - Udzha, 2 - Vilyui, 3 - Irkinei; 6 - basins: a - Chulman, б - Toka, в - Rybin, г - Irkutsk; 7 - Popigai astro-bleme; C.O. - orogen.
горячей точкой планеты, плюмом [Zonenshain, Kuz'min, 1976, 1983].
В целом же следует отметить, что результаты численного эксперимента и характеристики построенной модели конвекции в верхней мантии Евразии могут быть привлечены в дальнейшем для согласованных комплексных геодинамических интерпретаций результатов различных геолого-геофизических исследований, а сам алгоритм и комплекс программ - для расчетов разнообразных геодинамических обстановок на нашей планете.
6. ЗАКЛЮЧЕНИЕ
При решении трехмерных задач классической гидродинамики хорошо зарекомендовали себя численные алгоритмы с применением неявных методов расщепления по пространственным переменным и метода искусственной сжимаемости [Yanenko, 1971; Peyret, Taylor, 1983; Tarunin, 1990]. Уравнения конвективных процессов в верхней мантии Земли существенно отличаются от уравнений классической гидродинамики прежде всего бесконечными значениями числа Прандтля, значительным варьированием и нелинейной зависимостью вязкости мантийного вещества от литостатиче-ского давления и температуры. Представленная в настоящей работе численная модель конвекции в верхней мантии Земли под литосферой Евразийского континента построена с применением неявных методов расщепления по пространственным переменным и метода искусственной сжимаемости.
На основе детального анализа полученных результатов численных экспериментов можно сделать вывод о том, что мощность и геометрия блоков литосферы оказывают заметное влияние на распределение конвективных течений в верхней мантии, а установившаяся структура течений может явно прослеживаться в осредненном рельефе дневной поверхности платформенных областей, то есть больших по площади территорий с литосферой увеличенной мощности, таких как Восточно-Европейская и Сибирская платформы.
В предыдущих работах авторов на примере простейших геометрических объектов, прямоугольных в плане кратонов и ловушек [Tychkov et al., 2005a, 2005b, 2005c; и др.] было показано, что вариации мощностей и конфигурация объектов, внедренных в мантию, однозначно определяют структуру мантийных течений под ними и не зависят от начального распределения температуры. В случае численного моделирования конвекции под ровной литосферой и последующего внедрения в расчетную область кратонов и/или ловушек конвективное мантийное течение подстраивается под геометрические особенности объектов, и поэтому при правильном подборе геометрических характеристик, найденных на основе известных геолого-геофизических моделей, можно получить достаточно правдоподобные результаты.
Таким образом, детальный анализ результатов настоящей работы, а также результатов предыдущих
исследований авторов (например [Chervov et al., 2014; Bushenkova et al., 2018]) демонстрирует высокую эффективность разработанного подхода в решении задач взаимного влияния элементов литосферы и конвекти-рующей мантии, пригодного для расчетов разнообразных обстановок на нашей планете.
Комплексы программ расчетов были протестированы и запускались на ЭВМ Сибирского суперкомпьютерного центра ИВМиМГ СО РАН под управлением кластера НКС-1П и компилятора-Intel Parallel Studio XE 2017 Update 4, а также под управлением гибридного кластера НКС-30Т+GPU.
7. ЛИТЕРАТУРА/REFERENCES
Andreev VK, Kaptsov O.V., Pukhnachev V.V, Rodionov AA., 1998. Application of Group-Theoretical Methods in Hydrodynamics. Springer, Dordrecht, 397 p. https://doi.org/10. 1007/978-94-017-0745-9.
Artyushkov E.V, 1993. Physical Tectonics. Nauka, Moscow, 456 p. (in Russian) [Артюшков Е.В. Физическая тектоника. М.: Наука, 1993. 456 с.].
Blankenbach B., Busse F., Christensen U., Cserepes L., Gun-kel D., Hansen U., Harder H., Jarvis G., Koch M., Marquart G., Moore D., Olson P., Schmeling H., Schnaubelt T, 1989. A Benchmark Comparison for Mantle Convection Codes. Geophysical Journal International 98 (1), 23-38. https://doi.org/10.11 11/j.1365-246X.1989.tb05511.x.
Bobrov A.M., Trubitsyn V.P., 1995. Times of Mantle Flow Restructuring underneath Continents. Physics of the Earth 7, 5-13 (in Russian) [Бобров А.М., Трубицын В.П. Времена перестроек структуры мантийных течений под континентами // Физика Земли. 1995. № 7. С. 5-13].
Bogoyavlenskaya O.V., Puchkov V.N., Fedorov M.V, 1991. Geology of the USSR. Nedra, Moscow, 240 p. (in Russian) [Богоявленская О.В., Пучков В.Н., Федоров М.В. Геология СССР. М.: Недра, 1991. 240 с.].
Bushenkova N.A., 2004. Inhomogeneities of the Upper Mantle and the Modern Structure of the Lithosphere in Central Siberia from the Data of Seismic Tomography in Reflected Waves. BriefPhD Thesis (Candidate of Geology and Mineralogy). Novosibirsk, 20 p. (in Russian) [Бушенкова Н.А. Неоднородности верхней мантии и современная структура литосферы Центральной Сибири по данным сейсмото-мографии на отраженных волнах: Автореф. дис. ... канд. геол.-мин. наук. Новосибирск, 2004. 20 с.].
Bushenkova N.A., Deev E.V., Dyagilev G.S., Gibsher A.A., 2008. The Upper Mantle Structure and Cenozoic Volcanism of Central Mongolia. Doklady Earth Sciences 418 (1), 128131. https://doi.org/10.1007/s11471-008-1028-5.
Bushenkova N.A., Kuchai O.A., Chervov V.V., 2016. The Role of the Inhomogeneous Thickness of the Lithosphere in Subduction Processes: Comparison of Seismotomographic and Thermogravitational Models of the Upper Mantle with Seismicity and Seismotectonic Deformation on the Example of the Kamchatka Region and Japan. In: Tectonophysics and Topical Problems in Geosciences. Proceedings of 4th Tec-tonophysical Conference at Institute of Physics of the Earth (October 03-08, 2016). Vol. 1. IPE RAS Publishing House,
Moscow, p. 369-374 (in Russian) [Бушенкова НА., Кучай ОА., Червов В.В. Роль неоднородной мощности литосферы в процессах субдукции: сопоставление сейсмотомографи-ческой и термогравитационной моделей верхней мантии с характером сейсмичности и сейсмотектоническими деформациями на примере Камчатского региона и Японии / / Тектонофизика и актуальные вопросы наук о Земле: Материалы Четвертой тектонофизической конференции в ИФЗ РАН (03-08 октября 2016 г.). М.: Изд-во ИФЗ РАН, 2016. Т. 1. С. 369-374].
Bushenkova N.A., Kuchay O.A., Chervov V.V, 2018. Sub-meridional Boundary Zone in Asia: Seismicity, Lithosphere Structure, and the Distribution of Convective Flows in the Upper Mantle. Geodynamics & Tectonophysics 9 (3), 10071023 (in Russian) [Бушенкова Н.А., Кучай О.А., Червов В.В. Субмеридиональная пограничная зона в Азии: сейсмичность, структура литосферы и распределение конвективных потоков в верхней мантии // Геодинамика и тектонофизика. 2018. Т. 9. № 3. С. 1007-1023]. https://doi. org/10.5800/GT-2018-9-3-0381.
Busse F.H., Christensen U., Clever R., Cserepes L., Gable C., Giannandrea E., Guillou L., Houseman G., Nataf H.-C., Ogawa M., Parmentier M., Sotin C., Travis B., 1993. 3D Convection at Infinite Prandtl Number in Cartesian Geometry - a Benchmark Comparison. Geophysical & Astrophysical Fluid Dynamics 75 (1), 39-59. https://doi.org/10.1080/03091929 408203646.
Chervov VV, 2002a. Numerical Modeling of Three-Dimensional Convection Problems in the Earth's Mantle Using Vor-ticity and Vector Potential. Computational Technologies 7 (1), 114-125 (in Russian) [Червов В.В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением завихренности и векторного потенциала // Вычислительные технологии. 2002. Т. 7. № 1. С. 114-125].
Chervov VV, 2002b. Numerical Modeling of Problems of Three-Dimensional Convection in the Earth's Mantle Using Grid Sequences. Computational Technologies 7 (3), 85-92 (in Russian) [Червов В.В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением последовательности сеток // Вычислительные технологии. 2002. T. 7. № 3. С. 85-92].
Chervov VV., 2006. Modeling of Three-Dimensional Convection in the Earth's Mantle Using the Implicit Method of Splitting by Physical Processes. Computational Technologies 11 (4), 73-86 (in Russian) [Червов В.В. Моделирование трехмерной конвекции в мантии Земли с применением неявного метода расщепления по физическим процессам // Вычислительные технологии. 2006. Т. 11. № 4. С. 73-86].
Chervov VV., 2009. Modeling of Three-Dimensional Convection in the Earth's Mantle Using the Implicit Method of Weak Compressibility. Computational Technologies 14 (3), 86-92 (in Russian) [Червов В.В. Моделирование трехмерной конвекции в мантии Земли с применением неявного метода слабой сжимаемости // Вычислительные технологии. 2009. Т. 14. № 3. С. 86-92].
Chervov V.V., 2018. Software for Calculating Three-Dimensional Convection Underneath Continental Plates of the
Earth in Spherical Coordinates, Navie_Spherical_Coords/ 2017. Computer Software State Registration Certificate № 2018616280 Dated May 28, 2018. ROSPATENT, Moscow (in Russian) [Червов В.В. Программа расчета трехмерной конвекции под континентальными плитами Земли в сферических координатах Navie_Spherical_Coords/ 2017: Свидетельство о государственной регистрации программ для ЭВМ № 2018616280 от 28.05.2018. М.: РОСПАТЕНТ, 2018].
Chervov V.V., Chernykh G.G., 2014. Numerical Modeling of Three-Dimensional Convection in the Upper Mantle of the Earth beneath Eurasia Lithosphere. Journal of Engineering Thermophysics 23 (2), 105-111. https://doi.org/10. 1134/S1810232814020039.
Chervov V.V., Chernykh G.G., 2019. Numerical Modeling of Convection in the Zone of Spreading and Subduction. Journal of Engineering Thermophysics 28, 14-25. https://doi. org/10.1134/S1810232819010028.
Chervov VV, Chernykh G.G., Bushenkova N.A., Kulakov I.Yu., 2014. Numerical Modeling of Three-Dimensional Convection in the Upper Mantle of the Earth underneath the Lithosphere of Eurasia. Computational Technologies 19 (5), 101114 (in Russian) [Червов В.В., Черных Г.Г., Бушенкова Н.А., Кулаков И.Ю. Численное моделирование трехмерной конвекции в верхней мантии Земли под литосферой Евразии // Вычислительные технологии. 2014. Т. 19. № 5. С. 101-114].
Dobretsov N.L., Kirdyashkin A.G., Kirdyashkin A.A., 2001. Depth Geodynamics. GEO, Novosibirsk, 409 p. (in Russian) [Добрецов Н.Л., Кирдяшкин А.Г., Кирдяшкин А.А Глубинная геодинамика. Новосибирск: ГЕО, 2001. 409 c.].
Fedorenko R.P., 1964. The Speed of a Convergence of One Iterative Process. USSR Computational Mathematics and Mathematical Physics 4 (3), 227-235. https://doi.org/10. 1016/0041-5553(64)90253-8.
Heister T., Dannberg J., Gassmoller R., Bangerth W., 2017. High Accuracy Mantle Convection through Modern Numerical Methods - II: Realistic Models and Problems. Geophysical Journal International 210 (2), 833-851, https://doi.org/ 10.1093/gji/ggx195.
Jakovlev A.V, Bushenkova N.A., Koulakov I., Dobretsov N.L., 2012. Structure of the Upper Mantle in the Circum-Arctic Region from Regional Seismic Tomography. Russian Geology and Geophysics 53 (10), 963-971, https://doi.org/10. 1016/j.rgg.2012.08.001.
Landau L.D., Lifshits E.M., 1986. Hydrodynamics. Nedra, Moscow, 736 p. (in Russian) [Ландау Л.Д., Лифшиц Е.М. Гидродинамика. М.: Наука, 1986. 736 с.].
Lobkovsky L.I., Nikishin A.M., Khain V.E., 2004. Current Problems of Geotectonics and Geodynamics. Nauchny Mir, Moscow, 612 p. (in Russian) [Лобковский Л.И., Никишин А.М., Хаин В.Е. Современные проблемы геотектоники и геодинамики. М.: Научный мир, 2004. 612 с.].
Marchuk G.I., Shaidurov V.V., 1979. Improving the Accuracy of Solving Differentiation Schemes. Nauka, Moscow, 320 p. (in Russian) [Марчук Г.И., Шайдуров В.В. Повышение точности решения разностных схем. М.: Наука, 1979. 320 с.].
Peyret R., Taylor T.D., 1983. Computational Methods for Fluid Flow. Springer-Verlag, Berlin and Heidelberg, 358 p. https://doi.org/10.1007/978-3-642-85952-6.
Polyansky O.P., Prokop'ev A.V., Babichev A.V., Korobey-nikov S.N., Reverdatto V.V, 2013. The Rift Origin of the Vilyui Basin (East Siberia), from Reconstructions of Sedimentation and Mechanical Mathematical Modeling. Russian Geology and Geophysics 54 (2), 121-137. https://doi.org/10. 1016/j.rgg.2013.01.001.
Polyansky O.P., Prokopiev AV, Koroleva O.V, Tomshin M.D., Reverdatto V.V., Babichev A.V., Sverdlova V.G., Vasiliev D.A., 2018. The Nature of the Heat Source of Mafic Magmatism during the Formation of the Vilyui Rift Based on the Ages of Dike Swarms and Results of Numerical Modeling. Russian Geology and Geophysics 59 (1), 1217-1236. https://doi. org/10.1016/j.rgg.2018.09.003.
Polyansky O.P., Prokopiev AV, Koroleva O.V, Tomshin M.D., Reverdatto V.V., Selyatitsky A.Y, Travin A.V., Vasiliev D.A., 2017. Temporal Correlation between Dyke Swarms and Crustal Extension in the Middle Palaeozoic Vilyui Rift Basin, Siberian Platform. Litos 282-283, 45-64. https://doi.org/10.1016/ j.lithos.2017.02.020.
Rizzi A.W., Eriksson L.E., 1985. Computation of Inviscid Incompressible Flow with Rotation. Journal of Fluid Mechanics 153 (12), 275-312. https://doi.org/10.1017/S00 22112085001264.
Shatsky N.S., 1946. The Main Features of the Structure and Development of the East European Platform. Comparative Tectonics of Ancient Platforms. Bulletin of the USSR Academy of Sciences. Geological Series 1, 5-62 (in Russian) [Шатский Н.С. Основные черты строения и развития Восточно-Европейской платформы. Сравнительная тектоника древних платформ // Известия АН СССР. Серия геологическая. 1946. № 1. С. 5-62].
Tarunin E.L., 1975. The Method of Sequence of Nets for Free Convection Problems. USSR Computational Mathematics and Mathematical Physics 15 (2), 148-156. https://doi. org/10.1016/0041-5553(75)90049-X.
Tarunin E.L., 1990. Computational Experiment in Free Convection Problems. ISU Publishing House, Irkutsk, 228 p. (in Russian) [Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск: Изд-во ИГУ, 1990. 228 с.].
Trubitsyn V.P., Belavina Yu.F., Rykov V.V., 1994. Thermal Convection in the Mantle with Variable Viscosity and a Continental Plate of Finite Dimensions. Izvestiya, Physics of the Solid Earth 7, 5-17 (in Russian) [Трубицын В.П., Белави-на Ю.Ф., Рыков В.В. Тепловая конвекция в мантии с переменной вязкостью и континентальной плитой конечных размеров // Физика Земли. 1994. № 7. С. 5-17].
Trubitsyn V.P., Rykov VV., 1995. A 3D Numerical Model of the Wilson Cycle. Journal of Geodynamics 20 (1), 63-75. https://doi.org/10.1016/0264-3707(94)00029-U.
Tychkov S.A., Chervov VV, Chernykh G.G., 2005а. Numerical Modeling of 3D Convection in the Earth Mantle. Russian Journal of Numerical Analysis and Mathematical Modelling 20 (5), 483-500. https://doi.org/10.1163/156939805775186677.
Tychkov S.A., Chervov V.V, Chernykh G.G., 2005b. Numerical Modeling of Thermal Convection in the Earth's Mantle. Doklady Earth Sciences 402 (4), 596-601.
Tychkov S.A., Chervov V.V., Chernykh G.G., 2005c. A Numerical Model of Three-Dimensional Convection. Izvestiya, Physics of the Solid Earth 41 (5), 383-398.
Tychkov S.A., Chervov V.V., Chernykh G.G., 2007. Three-Dimensional Modeling of Convection underneath Cratons of Central Asia. In: Computational Technologies. Proceedings of the V Meeting of the Russia - Kazakhstan Working Group for Computational and Information Technologies (February 06-08, 2007). Vol. 12. (Spec. Iss. 4). Novosibirsk, p. 85-95 (in Russian) [Тычков С.А., Червов В.В., Черных Г.Г. Трехмерное моделирование конвекции под кратонами Центральной Азии // Вычислительные технологии: Труды V швещания российско-казахстанской рабочей группы по вычислительным и информационным технологиям (06-08 февраля 2007 г.). Новосибирск, 2007. Т. 12 (Спецвыпуск 4). С. 85-95].
Vladimirova N.N., Kuznetsov B.G., Yanenko N.N., 1966. Numerical Calculation of a Symmetric Flow of a Viscous Incompressible Fluid around a Plate. In: G.I. Marchuk (Ed.), Some Problems of Computational and Applied Mathematics. Nauka, Novosibirsk, p. 186-192 (in Russian) [Владимирова Н.Н., Кузнецов Б.Г., Яненко Н.Н. Численный расчет симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости // Некоторые вопросы вычислительной и прикладной математики / Ред. Г.И. Мар-чук. Новосибирск: Наука, 1966. С. 186-192].
Yanenko N. N., 1971. The Method of Fractional Steps: The Solution of Problems of Mathematical Physics in Several Variables. Springer-Verlag, Berlin and Heidelberg, 160 p. https://doi.org/10.1007/978-3-642-65108-3.
Zonenshain L.P., Kuz'min M.I., 1976. Global Tectonics, Magmatism and Metallogeny. Nedra, Moscow, 231 p. (in Russian) [Зоненшайн Л.П., Кузьмин М.И. Глобальная тектоника, магматизм и металлогения. М.: Недра, 1976. 231 c.].
Zonenshain L.P., Kuz'min M.I., 1983. Intraplate Volcanism and Its Significance for Understanding the Processes in the Earth's Mantle. Geotectonics 1, 28-45 (in Russian) [Зоненшайн Л.П., Кузьмин М.И. Внутриплитовый вулканизм и его значение для понимания процессов в мантии Земли // Геотектоника. 1983. № 1. С. 28-45].
