CC BY
12
GEODYNAMICS & TECTONOPHYSICS
Published by the Institute of the Earth's Crust, Siberian Branch, Russian Academy of Sciences
TECTONOPHYSICS
2021 VOLUME 12 ISSUE 3 PAGES 455-470
ISSN 2078-502X
DOI: 10.5800/GT-2021-12-3-0533
UPPER MANTLE CONVECTION RELATED TO SUBDUCTION ZONE AND APPLICATION OF THE MODEL TO INVESTIGATE THE CRETACEOUS-CENOZOIC GEODYNAMICS
OF CENTRAL EAST ASIA AND THE ARCTIC
L.I. Lobkovsky© 12, M.M. Ramazanov »3^ V.D. Kotelkin 14
1 Shirshov Institute of Oceanology, Russian Academy of Sciences, 36 Nahimovskiy Ave, Moscow 117997, Russia
2 Moscow Institute of Physics and Technology, 9 Institutskiy Ln, Dolgoprudny 141701, Russia
3 Institute for Geothermal Research and Renewable Energy, Branch of Joint Institute for High Temperatures of the Russian Academy of Sciences, 39а Shamil Ave, Makhachkala 367030, Republic of Dagestan, Russia
4 Lomonosov Moscow State University, 1 Leninskie Gory, Moscow 119991, Russia
ABSTRACT. A geodynamic model of upper mantle convection related to the Pacific subduction zone is mathematically substantiated and applied to investigate the Cretaceous-Cenozoic evolution of Central East Asia (CEA) and the Arctic. We present a solution for the two-dimensional stationary problem of thermal convection in the upper mantle layer, considering different Rayleigh numbers and taking into account the influence of the subduction process and lithospheric movements along the upper mantle base. We describe the results of 3D modeling of nonstationary upper mantle convection in a subduction zone. Our data give grounds to propose explanations for the entire spectrum of tectonic-magmatic processes developing within CEA in the Cenozoic and the Arctic in the Upper Cretaceous and Cenozoic. We discuss the reasons why the lithosphere in CEA and the Arctic is generally shifting towards the Pacific subduction zone, considering the presence of separate magmatic provinces and rift zones. In our opinion, this is due to the existence of a large horizontally elongated convective cell, which interior is composed of smaller isometric cells. This long cell creates the effect of conveyor dragging of the lithosphere, while its internal cells produce the effect of upper mantle plumes.
KEYWORDS: lithosphere; subduction zone; upper mantle; thermal convection; 2D convection model; 3D modeling; Cretaceous-Cenozoic evolution; Central East Asia; Arctic; Baikal rift zone
FUNDING: The study was performed under the state assignment of Shirshov Institute of Oceanology (project 01282021-0004) and partially supported by the Russian Foundation for Basic Research (project 18-05-70012 - Arctic).
RESEARCH ARTICLE Received: December 19, 2020
FOR CITATION: Lobkovsky L.I., Ramazanov M.M., Kotelkin V.D., 2021. Convection related to subduction zone and application of the model to investigate the Cretaceous-Cenozoic geodynamics of Central East Asia and the Arctic. Geodynamics & Tectonophysics 12 (3), 455-470. doi:10.5800/GT-2021-12-3-0533
Correspondence: Mukamay M. Ramazanov, mukamay-ipg@mail.ru
Revised: May 3, 2021 Accepted: May 12, 2021
РАЗВИТИЕ МОДЕЛИ ВЕРХНЕМАНТИЙНОЙ КОНВЕКЦИИ, СОПРЯЖЕННОЙ С ЗОНОЙ СУБДУКЦИИ, С ПРИЛОЖЕНИЯМИ К МЕЛ-КАЙНОЗОЙСКОЙ ГЕОДИНАМИКЕ ЦЕНТРАЛЬНО-ВОСТОЧНОЙ АЗИИ И АРКТИКИ
Л.И. Лобковский1,2, М.М. Рамазанов3, В.Д. Котелкин1,4
1 Институт океанологии им. П.П. Ширшова РАН, 117997, Москва, пр-т Нахимовский, 36, Россия
2 Московский физико-технический институт, 141701, Долгопрудный, Институтский пер., 9, Россия
3 Институт проблем геотермии и возобновляемой энергетики, филиал Объединенного института высоких температур РАН, 367030, Махачкала, пр-т И. Шамиля, 39а, Республика Дагестан, Россия
4 Московский государственный университет им. М.В. Ломоносова, 119991, Москва, Ленинские горы, 1, Россия
АННОТАЦИЯ. Рассматривается математическое обоснование геодинамической модели верхнемантийной конвекции, сопряженной с Тихоокеанской зоной субдукции, в приложении к мел-кайнозойской эволюции Центрально-Восточной Азии (ЦВА) и Арктики. Приводится решение двухмерной стационарной задачи термической конвекции в слое верхней мантии при различных числах Рэлея с учетом влияния процесса субдукции и движения литосферного слоя вдоль подошвы верхней мантии. Описываются результаты 3D-моделирования нестационарной конвекции в верхней мантии, сопряженной с зоной субдукции. Полученные результаты позволяют объяснить весь спектр наблюдаемых тектономагматических процессов, развивающихся в пределах ЦВА в кайнозое и Арктики в верхнем мелу и кайнозое, а именно сочетание общего смещения литосферы ЦВА и Арктики в сторону Тихоокеанской зоны субдукции с наличием отдельных магматических провинций и рифтовых зон как следствие существования длинной горизонтально вытянутой конвективной ячейки (создающей эффект конвейерного волочения литосферы), осложненной внутренними изометричными ячейками (создающими эффект верхнемантийных плюмов).
КЛЮЧЕВЫЕ СЛОВА: литосфера; зона субдукции; верхняя мантия; тепловая конвекция; двухмерная модель конвекции; 3D моделирование; мел-кайнозойская эволюция; Центрально-Восточная Азия; Арктика; Байкальская рифтовая зона
ФИНАНСИРОВАНИЕ: Работа выполнена в рамках госзадания Института океанологии им. П.П. Ширшова РАН № 0128-2021-0004 и частично по теме гранта РФФИ «Арктика» № 18-05-70012.
1. ВВЕДЕНИЕ
В работах [Lobkovsky et al., 2011, 2013; Laverov et al., 2013; Lobkovsky, 2016] рассматривалась новая геодинамическая модель эволюции Арктического региона для верхнего мела - кайнозоя и эволюции Центрально-Восточной Азии для кайнозоя на основе представлений о верхнемантийной конвекции, сопряженной с Тихоокеанской зоной субдукции, характеризуемой наличием сильно вытянутой горизонтальной ячейки, в которой субдуцируемая литосфера продолжает горизонтальное движение под континент вдоль подошвы верхней мантии, компенсируемое возвратным подлитосферным потоком астеносферы, направленным от внутрикон-тинентальной области в сторону Тихого океана (рис. 1). Эта модель развивает применительно к Арктике и Центрально-Восточной Азии гидродинамический подход, описанный в работах [Kirdyashkin, Dobretsov, 1991; Юг-dyashkin et al., 2000, 2002] при анализе горизонтально протяженных конвективных ячеек в астеносфере, возникающих при наличии горизонтального градиента температуры.
В рамках гидродинамического подхода в нашей модели естественным образом трактуется выявленный методами сейсмотомографии горизонтальный слой повышенных сейсмических скоростей протяженностью
более тысячи километров, примыкающий к подошве верхней мантии и соединенный с погружающейся литосферой [Lobkovsky, 2016], как слой субдуцирован-ного литосферного вещества, вовлеченный в конвективную циркуляцию и движущийся от зон субдукции Тихоокеанской литосферы в сторону континентальных областей Арктики, Центрально-Восточной Азии (ЦВА). При этом следует отметить, что многие геологи и геофизики придерживаются принципиально иной интерпретации выявленного аномального слоя, исходя из представлений о так называемой «стагнирующей» на подошве верхней мантии литосфере, попавшей туда через зону субдукции и, в силу каких-то причин, потерявшей способность к дальнейшему движению [Zhao, 2009; Zorin, Turutanov, 2005; Zorin et al., 2006].
В развитие нашей геодинамической модели мы поставили перед собой задачу с помощью математического моделирования найти и количественно обосновать те геодинамические условия, при которых имеет место общая циркуляция вещества в горизонтально протяженной верхнемантйной ячейке, сопряженной с зоной субдукции, осложненная возникновением системы изометричных ячеек в пределах конвектирую-щего слоя, гидродинамически играющих роль верхнемантийных плюмов, проявляющихся в существовании
отдельных магматических провинций. Моделирование в основном будет соотноситься с особенностями строения и кайнозойской геодинамики областей ЦВА, характеризуемых формированием в них рифтовых зон и вулканических областей. На рис. 2 показаны смещения поверхности Земли по данным GPS на обширной территории ЦВА от Байкала до Японской островной дуги и места проявления голоценового и четвертичного магматизма. Эти данные прямо указывают на характерную геодинамическую обстановку, в которой сочетается общее смещения литосферы ЦВА в сторону Тихоокеанской зоны субдукции с наличием отдельных магматических провинций и рифтовых зон, что, по нашим представлениям, является следствием существования длинной горизонтально вытянутой конвективной ячейки (создающей эффект конвейерного волочения литосферы), осложненной внутренними изометрич-ными ячейками (создающими эффект верхнемантийных плюмов).
Карта на рис. 2 построена на основе исходных данных ГНСС (глобальные навигационные спутниковые системы), приведенных в работах [Kogan et al., 2003; Steblov et al., 2003; Bürgmann et al., 2005; Apel et al., 2006; Meng et al., 2006; Lobkovsky et al., 2017, 2018]. Результаты ГНСС-измерений представлены в системе координат ITRF - международной земной отсчетной основы [Altamimi et al., 2016].
Несмотря на большой прогресс в изучении глубинного строения литосферы данного региона, связанный,
в частности, с применением методов сейсмической томографии верхней и нижней мантии, а также успехи космических технологий в прямых измерениях поверхностных деформаций и смещений земной коры, у исследователей до сих пор не сформировалась единая точка зрения на геодинамические процессы, определяющие эволюцию ЦВА. Это видно из большого количества публикаций, в которых дается различная геодинамическая интерпретация полученных за последние 50 лет геолого-геофизических данных по ЦВА [Florensov, 1968; Logatchev, 1968; Molnar, Tapponnier, 1975; Zonenshain et al., 1978; Artyushkov et al., 1990; Bott, 1990; Ruzhich, 1997; Ruzhich et al., 2016; Zorin, Turutanov, 2005; Zorin et al., 2006; Zhao et al., 2006; Djadkov et al., 1999; Sankov et al., 2011; Seminsky et al., 2013; Rasskazov, Chuvashova, 2013; Chuvashova et al., 2017; Yarmolyuk et al., 1995, 2013; Mats, 2015]. Такая концептуальная неоднозначность вызвана объективными обстоятельствами и, в первую очередь, связана с многофакторным характером геологической эволюции региона ЦВА и его достаточно сложным глубинным строением.
В последние годы появились работы, в которых предлагаются возможные сценарии эволюции региона ЦВА, исходя из современных представлений тектоники плит с привлечением последних данных по сейсмотомогра-фии коры и верхней мантии, геохронологии магматизма, деформациям коры и осадочной толщи, сейсмологическим измерениям, космическим наблюдениям поверхностных смещений коры и т.д. [Yarmolyuk et al., 2013;
Датун
400
1С
Чанбай ▲
Япония А
800
1200 100°
108°
116°
124°
132°
140°
148° в.д.
Рис. 1. Один из вертикальных разрезов мантии в виде сейсмотомограммы [Zhao, 2009; Zhao et al., 2009, 2010] (а) и модель мантийной конвекции, сопряженной с зоной субдукции [Lobkovsky, 2016; Laverov et al., 2013; Lobkovsky et al., 2011, 2013] (6). Fig. 1. Seismogram showing a vertical section of the mantle [Zhao, 2009; Zhao et al., 2009, 2010] (a); model of mantle convection related to a subduction zone [Lobkovsky, 2016; Laverov et al., 2013; Lobkovsky et al., 2011, 2013] (6).
0
Я^ЫсИ et а1., 2016; СИиуаБИоуа et а1., 2017]. В этих работах описываются возможные глубинные течения вещества в верхней и нижней мантии, объясняющие установленные закономерности распределения магматизма, рифтовых зон, складчатых областей и других тектонических обстановок в пространстве и времени. Однако следует отметить, что указанные работы носят преимущественно качественный описательный характер и это не дает возможности объективно оценить степень адекватности предложенных в них моделей движения вещества в мантии. Вместе с тем они позволяют обосновать соответствующие постановки задач математического моделирования изучаемых процессов, решение которых может подкрепить или существенно скорректировать эти модели.
Целью нашей работы является математическое моделирование эволюции геодинамических процессов, происходящих главным образом в ЦВА в позднемезо-зойско-кайнозойский период, для которого удается достаточно обоснованно сформулировать постановку математической задачи моделирования течений в верхней мантии, исходя из имеющейся совокупности геолого-геофизических данных и геодинамических представлений, развитых в отмеченных выше работах, с акцентом на объяснение происхождения и эволюции Байкальской рифтовой зоны (БРЗ), являющейся в известной степени ключевым элементом для понимания кайнозойской геодинамики ЦВА. Полученные результаты моделирования также применимы к анализу эволюции литоферы Арктики для верхнего мела - кайнозоя.
60° с.ш.
50°
40°
60°
50°
40°
110°
120°
130°
140°
150°
20 мм/год (скорости смещения пунктов GPS) Границы литосферных плит
I
Голоценовый и/или активный вулканизм А Четвертичный вулканизм
Рис. 2. Поле горизонтальных компонент векторов GPS-скоростей пунктов наблюдения в Центрально-Восточной Азии. Скорости приведены в системе координат ITRF. EUR - Евразийская плита, PAC - Тихоокеанская плита, NAM - Североамериканская плита. Обозначения вулканических полей и регионов: Ch - Чангбайшан, D-A - Дарибанга-Абага, H - Хангай, SWB -юго-западная часть Байкальскго рифта, D - Даурия, V - Витим, U - Удокан. Fig. 2. Horizontal components of GPS velocity vectors in Central East Asia.
Velocity values in the ITRF coordinate system. Plates: EUR - Eurasian, PAC - Pacific, NAM - North-American. Volcanic fields and regions: Ch - Changbaishan, D-A - Daribanga-Abaga, H - Hangai, SWB - southwestern part of the Baikal rift, D - Dauria, V - Vitim, U - Udokan.
2. ГЕОЛОГИЧЕСКИЕ И ГЕОДИНАМИЧЕСКИЕ ПРЕДПОСЫЛКИ ДЛЯ ПОСТАНОВКИ МАТЕМАТИЧЕСКОЙ ЗАДАЧИ МОДЕЛИРОВАНИЯ ВЕРХНЕМАНТИЙНОЙ КОНВЕКЦИИ ПОД ЦЕНТРАЛЬНОЙ И ВОСТОЧНОЙ АЗИЕЙ
В работах [Molnar, Tapponnier, 1975; Zonenshain et al., 1978; Yarmolyuk et al., 2013; Ruzhich et al., 2016; Chuva-shova et al., 2017; и др.] анализируются различные геодинамические аспекты влияния межплитных процессов, происходящих в области Индо-Азиатской коллизии и зоне Тихоокеанской субдукции, на внутриплитные процессы в регионе ЦВА. При рассмотрении геодинамического влияния Индо-Азиатской коллизии на регион ЦВА обычно ссылаются на широко известную модель [Molnar, Tapponnier, 1975] о жестком «инденторе» в виде фронтальной части Индийского континента, который вдавливается в Азиатский континент в процессе коллизии, приводя к деформациям и смещениям блоков литосферы ЦВА. При этом деформации и смещения блоков происходят не только в окрестности движущегося фронта индентора, но и вдали от него на расстояниях нескольких тысяч километров вплоть до БРЗ, возникшей, как предполагается, в результате такого процесса вдавливания. Следует, однако, отметить, что в последнее время появляется все больше работ, в которых на основе детального анализа геолого-геофизических данных ставится под сомнение модель образования БРЗ под действием коллизии Индийского континента с Евразией [Mats, 2015; Ruzhich et al., 2016].
Обратимся к другой межплитной границе геодинамического влияния на регион ЦВА - зоне субдукции Тихоокеанской литосферы. Многие исследователи, изучающие эволюцию ЦВА, придают большое значение влиянию Тихоокеанской конвергентной границы на эволюционные процессы в этом регионе [Yarmolyuk et al., 2013; Ruzhich et al., 2016; Chuvashova et al., 2017], предлагая при этом разные геодинамические модели такого влияния. В этой связи следует особо отметить получившую в последние годы популярность модель o стагнирующей литосфере на подошве верхней мантии, попавшей туда через зону субдукции [Zhao, 2009; Zhao et al., 2006]. Эта модель возникла в результате анализа данных по сейсмической томографии верхней мантии в переходной зоне от Тихого океана к Азиатскому континенту, выявивших сейсмическую картину о высокоскоростном горизонтальном участке слоя, расположенном на подошве верхней мантии, который является непрерывным продолжением погружающегося в зоне субдукции слоя литосферы [Zhao, 2009]. Эта гипотеза получила широкое признание среди геологов и геофизиков, изучающих регион ЦВА, несмотря на ее достаточно умозрительный и в известной степени противоречивый характер. Дело в том, что мгновенная «статическая» картина сейсмической томографии мантии при ее динамической интерпретации, включающей анализ движения вещества, должна учитывать законы гидродинамики вязкой несжимаемой жидкости, описывающей мантийные течения. С этой точки зрения
длительная стагнация литосферного участка на подошве верхней мантии, являющегося непрерывным продолжением погружающейся Тихоокеанской литосферы в зоне субдукции, означала бы последовательное накопление литосферного материала и разбухание зоны стагнации в виде громадной скученной литосферной массы, толщина которой в зоне скопления должна была бы значительно превосходить смежные участки слоя. Однако сейсмотомографическая картина показывает наличие однородного горизонтального слоя литосфер-ного материала постоянной толщины [Zhao, 2009]. Такая однородная геометрия слоя скорее свидетельствует в пользу его равномерного распространения в сторону континента со скоростью, сопоставимой с погружением литосферы в зоне субдукции. В этом случае, исходя из гидродинамики несжимаемой жидкости, естественно предположить, что движущееся в сторону континента вдоль подошвы верхней мантии бывшее литосферное вещество вовлечено в циркуляционную конвективную ячейку, у которой существует верхняя возвратная ветвь движения мантийного вещества, текущего в сторону Тихоокеанской зоны субдукции (см. рис. 1).
Такая модель конвективной верхнемантийной ячейки, сопряженной с зоной субдукции, была предложена Л.И. Лобковским с соавторами при геодинамическом анализе эволюции литосферы Арктики для верхнего мела и кайнозоя [Lobkovsky et al., 2011, 2013; Laverov et al., 2013] и использовалась для объяснения различных особенностей строения и геодинамики Центрально-Восточной Азии [Lobkovsky, 2016].
Необходимо отметить, что влияние процесса суб-дукции на структуру конвективных течений в астеносфере под континентом ранее исследовалось в ряде работ аналитически [Kirdyashkin et al., 2000], экспериментально [Kirdyashkin et al., 2002] и численно [Lobkovsky et al., 2014; Kotelkin, Lobkovsky, 2019].
Так, в работах [Kirdyashkin, Dobretsov, 1991; Kirdyashkin et al., 2000] было показано, что субдуцируемая плита за счет охлаждающего эффекта создает горизонтальный градиент температуры в прилегающем слое астеносферы, который приводит к образованию горизонтально вытянутой конвективной ячейки с возвратным верхним потоком, текущим в сторону зоны субдук-ции. Эта модель была применена к случаю Андийской континентальной окраины, причем горизонтальная астеносферная ячейка протягивалась под Атлантический океан вплоть до зоны спрединга. В работе [Lobkov-sky et al., 2014] численно исследовалась нестационарная картина субдукции с меняющимся углом наклона плиты и соответствующими конвективными течениями в смежной области верхней мантии.
Настоящая работа является развитием указанных выше исследований в приложении к региональной геодинамической ситуации в ЦВА, включая анализ различных режимов конвекции, объясняющих наблюдаемые особенности рифтогенеза, магматизма и деформаций литосферы в районе БРЗ и на сопредельных территориях ЦВА.
Совокупность имеющихся на сегодня сейсмотомо-графических данных о строении мантии в зоне перехода от Тихого океана к Восточной и Северо-Восточной Азии, геологических представлений об эволюции ЦВА, а также теоретических моделей развития горизонтальных конвективных ячеек в верхней мантии под окраин-но-континентальной областью литосферы, примыкающей к зоне субдукции, позволяет нам сформулировать постановку математической задачи о мантийной конвекции под регионом ЦВА, сопряженной с Тихоокеанской субдукцией.
При формулировке математической задачи в первую очередь необходимо определиться с ее начальными и граничными условиями, выбор которых может быть сделан, исходя из существующих геологических представлений об эволюции региона ЦВА. Поскольку в данном случае нас будет интересовать в основном кайнозойская эволюция БРЗ, обратимся к геологическим моделям эволюции региона ЦВА в наиболее близкий к кайнозою позднемезозойский период времени, описанным в работах [Уагшо1уик et а1., 2013; ЯаБ5ка7оу СИиуаБИоуа, 2013; СИиуаБИоуа et а1., 2017].
Согласно [Уагшо1уик et а1., 2013], в начале позднего мезозоя (160-150 млн лет назад) начался новый этап тектономагматической активности в ЦВА, который был отмечен возникновением целого ряда рифтовых зон: Западно-Забайкальской, Южно-Хангайской, Восточно-Монгольской. Наибольшая вулканическая активность в ЦВА имела место с раннего кайнозоя, что выразилось в появлении новых вулканических областей, в частности Южно-Байкальской, Удоканской, Дариганга и др. С этими представлениями хорошо согласуются геофизические данные по глубинному строению региона ЦВА, особенностью которого является общее поднятие кровли астеносферы до глубины менее 100 км ро-пп, Тиг^апоу, 2005; 7опп et а1., 2006]. Это региональное поднятие астеносферы осложняется локальными астеносферными выступами, приуроченными к областям новейшего вулканизма и достигающими глубины менее 50 км от поверхности Земли ропп, ТигШапоу 2005; Мо^утоуа et а1., 2008]. В работе [Уагшо1уик et а1., 1995] региональное поднятие астеносферы рассматривается как Центрально-Азиатское горячее поле мантии, связанное с областью концентрации плюмов, являющейся континентальным продолжением Тихоокеанского суперплюма.
Другая трактовка тектономагматических событий в ЦВА дается в работе [СИиуаБИоуа et а1., 2017], в которой важное значение придается так называемым первичным расплавным аномалиям переходного слоя мантии: Гобийской, Западно-Забайкальской и СевероЗабайкальской. Авторы исходят из предположения, что переходный мантийный слой был нарушен в начале новейшего геодинамического этапа около 90 млн лет назад вследствие развития явления мантийного «ава-ланша» - лавинного погружения в нижнюю мантию тяжелой субдуцированной литосферной массы, которая достаточно долго накапливалась и стагнировала
в низах верхней мантии при закрытии Урало-Монгольского и Солонкерского палеоокеанов, а также Монголо-Охотского залива Палеопацифики. Согласно [СИиуа-БИоуа et а1., 2017], после завершения процесса лавинного погружения литосферных масс в нижнюю мантию в верхней мантии под литосферой Азии возникло обратное течение в сторону Тихого океана в пределах Японско-Байкальского коридора, которое привело к растяжению литосферы в районе расположения первичной Западно-Забайкальской расплавной зоны и возникновению Байкальского рифта. Дальнейшее развитие процесса привело к образованию вторичных расплавных областей и соответствующих им магматических провинций: Саянской, Южно-Байкальской, Витимской и других [СИиуаБИоуа et а1., 2017]. Указанная картина возвратного верхнемантийного подлитосферного течения, объясняющая возникновение БРЗ, в общих чертах близка к описанной геодинамической модели образования БРЗ, изложенной в работе [ЬоЬкоуБку, 2016].
При постановке начальных и граничных условий задачи математического моделирования кайнозойской эволюции региона ЦВА и возникновения БРЗ главными являются следующие геодинамические условия, вытекающие из геологических представлений о докай-нозойской эволюции региона ЦВА [СИиуаБИоуа et а1., 2017; Уагшо1уик et а1., 2013]: 1) процесс кайнозойской эволюции региона можно рассматривать, ограничиваясь движениями только в верхней мантии, так как взаимодействие с нижней мантией в виде идущих из нее плюмов [Уагшо1уик et а1., 2013] или, наоборот, лавинного погружения в нее литосферных масс [СИиуаБИоуа et а1., 2017] произошло в докайнозойский период; 2) на восточной границе области моделирования в течение кайнозоя имел место устойчивый процесс субдукции Тихоокеанской литосферы; 3) на западной границе нашей области моделирования вблизи края Сибирского кратона к началу кайнозоя существовал «первичный» верхнемантийный плюм, отвечающий Западно-Забайкальской вулканической провинции [СИиуаБИоуа et а1., 2017; Уагшо1уик et а1., 2013].
Таким образом, мы приходим к постановке математической задачи о развитии конвекции в протяженном горизонтальном слое верхней мантии под регионом ЦВА, расположенным между зоной Тихоокеанской субдукции и краем Сибирского кратона, находящимся в условиях горизонтального градиента температуры, обусловленного, с одной стороны, пониженной температурой восточного края области из-за субдукции холодной Тихоокеанской литосферы, а с другой - повышенной температурой западного края области в силу наличия там Западно-Забайкальского верхнемантийного плюма. Фактически мы будем рассматривать двумерную модель мантийных течений в широтном сегменте, ограниченном с севера примерно 57-м градусом северной широты, с юга - 30-м градусом северной широты, с востока - Японской и Курило-Камчатской островными дугами и с запада - примерно 97-м градусом восточной долготы (рис. 2).
Специального рассмотрения требует вопрос о граничных условиях на подошве области моделирования. На сейсмотомографических разрезах верхней мантии (см. рис. 1, а) видно, что горизонтальный холодный слой, включающий литосферное вещество, имеет значительную толщину - порядка 200 км, примерно вдвое превосходящую мощность погружающегося в зоне суб-дукции слоя литосферы толщиной около 100 км. Эту разницу в толщине наклонного и горизонтального ли-тосферных слоев можно объяснить тем, что наклонный субдуцирующий слой движется вниз как «слэб» с однородной эпюрой скоростей, в то время как горизонтальный слой прилипает снизу к фазовой границе 670 км между верхней и нижней мантией, скорость на которой равна нулю. Отсюда следует, что если скорость на верхней границе горизонтального слоя примерно совпадает со скоростью погружающейся плиты (что вполне естественно допустить), то в условиях несжимаемости и неразрывности среды толщина горизонтального слоя должна быть вдвое больше мощности субдуцируемого слэба, что и наблюдается на упомянутых выше сейсмотомографических разрезах верхней мантии (см. рис. 1, а). Из этих разрезов также видно, что горизонтальный литосферный слой занимает пространство между границами двух фазовых переходов в верхней мантии - на глубине 410 и 670 км. Получается, что этот слой распространяется в переходной части верхней мантии и отделяется от вышележащей мантийной области фазовой границей на глубине 410 км, где происходит скачкообразное увеличение плотности и вязкости среды [Kirdyashkm, Dobretsov, 1991], поэтому, если ограничить область моделирования конвекции слоем верхней мантии выше фазовой границы 410 км, то в качестве граничного условия на нижней границе этой области следует задавать горизонтальную скорость движения вещества, которая равняется скорости погружения литосферы в зоне субдукции. В такой постановке задачи могут возникать различные ситуации в процессе взаимодействия горизонтальной конвективной ячейки с движущимся под ней литосфер-ным слоем, который может играть либо активную, либо пассивную роль в развитии конвективного процесса. В первом случае скорость движения подстилающего ли-тосферного слоя больше характерной скорости в конвективной ячейке, вызванной горизонтальным градиентом температуры, и литосфера снизу действует как
активный фактор, усиливающий естественную циркуляцию вещества. Во втором случае, наоборот, более интенсивная конвективная циркуляция вещества как бы затягивает медленно движущуюся литосферу под континент. Эти и другие возможные режимы верхнемантийной конвекции, сопряженной с процессом субдук-ции, будут рассмотрены ниже.
3. ПОСТАНОВКА МАТЕМАТИЧЕСКОЙ ЗАДАЧИ О ВЕРХНЕМАНТИЙНОЙ КОНВЕКЦИИ, СОПРЯЖЕННОЙ С ЗОНОЙ СУБДУКЦИИ
Для простоты область, показанную на рис.1, б, и моделирующую континентальную астеносферу с кровлей на глубине 100 км и подошвой на глубине 410 км, будем считать прямоугольной, т.е. примем, что океаническая плита (холодный поток жидкости) в зоне субдук-ции погружается вертикально вниз и затем движется влево (рис. 3). Будем рассматривать квазистационарную двухмерную модель, считая, что внешние параметры задачи меняются во времени медленно по сравнению со скоростями конвективных движений. К указанным внешним параметрам относится, например, положение зоны субдукции. Будем рассматривать для астеносферы реологию вязкой жидкости. Горизонтальные границы рассматриваемой области считаются твердыми, и на них задаем условие прилипания. При этом на нижней фазовой границе 410 км задаем условие прилипания к движущейся со скоростью v1 влево литосфере, хотя, как следует из рис. 1, такое условие следовало бы задавать лишь на части основания, занимаемой в данный момент времени литосферным слоем. Однако для оценки максимального эффекта движения плиты будем это условие задавать на всем основании. Граничное условие для скорости на правой границе (рис. 3) отражает условие скольжения, хотя более адекватным было бы использовать условие прилипания к погружающейся плите. Однако можно показать [Lobkovsky Ramazanov, 2021], что картина конвекции от такой подмены слегка меняется лишь в окрестности плиты, а в основной области ее характеристики от этого не зависят.
Рассмотрим условия для температуры. В результате вертикального погружения холодной плиты можно считать, что температура на правой границе слабо меняется, примем ее равной 1100 °С. На нижней границе температура растет от 1100 °С до невозмущенного значения, равного 1500 °С. Для простоты будем считать
v,=0;
v. =0; y, = 0; Т = -x/4L
= 0;
—L
v2 =0; vx = -v,; Т = -x/L
О х
Рис. 3. Постановочный рисунок к задаче с указанием граничных условий, L=10. Fig. 3. Sketch for the problem and boundary conditions, L=10.
указанный рост линейным. Полагаем, что на левой вертикальной границе температура линейно убывает от 1500 до 1200 °С. Наконец, примем, что на верхней границе температура при движении к зоне субдукции линейно убывает от 1200 до 1100 °С. Все эти граничные условия в безразмерном исчислении показаны на рис. 3. При этом температура на правой границе принята за условный ноль, т.е. температура отсчитывается от температуры на правой границе, а за масштаб температуры взята величина 400 °С.
Таким образом, рассматриваемая конвекция в общем случае обусловлена тремя факторами: горизонтальным перепадом температуры, связанным с погружением холодной плиты в зоне субдукции, вертикальным перепадом температуры и движением плиты вдоль основания. Рассматривается и частный случай отсутствия третьего фактора. Все эти факторы учитываются граничными условиями.
Рассматриваемая квазистационарная задача в приближениях Стокса и Буссинеска подчиняется следующей системе уравнений:
-Vp + mAv + ро [1 - в(T - T0) g = 0,
Vv = 0, р0 = const, (1)
vVT = xAT, x = .
P0Cp
Здесь: v - поле скоростей; T - температура; ^ - динамическая вязкость; в - коэффициент теплового расширения; g - ускорение свободного падения; Я - коэффициент теплопроводности; Cp - удельная изобарическая теплоемкость; x - температуропроводность.
Введем функцию тока ф и систему координат, как показано на рис. 3. Используя следующие характерные масштабы величин: h - толщина верхней мантии; р0 -плотность среды; Я/pjCh - скорость; рдСЪ.2/Я - время; AT - перепад температур на левой границе слоя, запишем систему (1) в безразмерном виде:
dT dT dT 32T д2Т
I-, vx--+ vz-=--1--,
дх х дх z dz дх2 dz2 дф дф
dz дх
Д2ф
(2)
pPCp gßATh Ra =-—-; Ra ■
Хц
Граничные условия
число Рэлея.
д2ф
x = 0: ф = 0, —2 = 0, T = 0;
дх2
Я 2
x = -L : ф = 0, = 0, T = 1 ;
дх
z
0: ф = z = 1: ф
0, д-Ф =
dz
=0, дФ=
dz
T = -- ;
0, T =--
3z
4
x
L '
x 4L
(3)
4. РЕШЕНИЕ МАТЕМАТИЧЕСКОМ ЗАДАЧИ
Полное решение данной математической задачи на основе аналитических и численных методов изложено
в статье [ЬоЬкоуБку, Яагс^апоу, 2021]. Следуя этой работе, представим решение в виде конечных сумм Фурье по координате х следующим образом:
(3z - 4)х
Т =---+ 4 sinnzsinakx,
i=i k=i
ф=èФк(z)sinakx' ak = Пк.
k=1 L
(4)
Здесь числа слагаемых в рядах Ых, произвольны, но конечны и фиксированы. Они выбираются исходя из необходимой точности.
В (4) неизвестными являются функции фк(1), к=1,2,... и постоянные числа ¡¡кк, /=1,2,..., к=1,2,... . При этом температура удовлетворяет всем граничным условиям, а функция тока пока только условиям на вертикальных границах. Остальным условиям необходимо удовлетворить с помощью функций фк(Х), к=1,2,.... Подставляя (4) в первое уравнение в системах (2), (3) получим, уравнения относительно этих функций. Пропуская промежуточные выкладки, которые приведены в [ЬоЬкоуБку, Яагс^апоу, 2021], приведем конечный результат:
фк = Ra
ф1к + 22 фк1ашс km ^in
+ v1 ф2к, k = 1,2,...N (5)
ф1к = [ Ak (z -1) + A
chakz + chak (z-1)
1 + chak
+C1k (z -1)shakz + C2kzshak (z -1),
kshak A
k^k + AkShak,
1k (sh2ak -ak2)(1 + chak)
C = akshak { Akshak + A2kak )
2k (sh2ak - ak2 )(l + chak )
_ nk _ 4(1 - cosakL)
ak _~T ' Ak _ 573 '
L akL
ak (l — cosakL)
Ak = , k =1,2,...
Ф
5 r3 ak L
sin mz
( 2 • 2 . (7Г / +(
- + Bkizshak (z -1) + Cki(z - 1)shakz,
nicha,
(\ 2
n2i2 + ak2 ) (shakchak — ak ) 1 „ nk
-—r— Bki, ak =—, k = 1,2,... cha. L
„ 1 — cosnkr„ , , 1
ф2к = 2-Пк-\-C3kZShak (z — 1) + C4k (z — 1)shakzJ '
cha,
shakchak — ak
shakchak — ak
f, k = 1,2,..
В решении (4), (5) неизвестными остались лишь постоянные числа ^ , /=1,2,.Nк=1,2,...Ых. Для их
1=1 m=1
1
нахождения необходимо использовать уравнение переноса тепла.
Рассмотрим это уравнение:
оф от _дф дТ=АТ (6)
дz дх дх дz
Данное уравнение решалось численно методом прогонки. Таким образом, суть примененного в работе метода заключается в следующем: на текущем временном шаге при заданной матрице поле скоростей определяется аналитическими выражениями (5). Затем при заданном поле скоростей нестационарное уравнение переноса тепла на данном временном шаге решается численно методом прогонки. Далее из поля температуры Т , используя (4) и обратное преобразование Фурье, находим новые значения (кГ После определенного числа итерационных шагов переходим на новый временной шаг. Процесс продолжается до выхода на стационарное решение.
5. ОБСУЖДЕНИЕ ПОЛУЧЕННЫХ РЕЗУЛЬТАТОВ
Приведем прежде всего некоторые оценки и соображения, позволяющие обосновать отдельные моменты в постановке задачи. Выбранное распределение температуры на правой и нижней границах основано, в частности, на сейсмотомографическом разрезе верхней мантии (см. рис. 1, а). Из этого разреза видно, как погружающееся в зону субдукции холодное вещество, достигая переходной зоны между верхней и нижней мантией, поворачивает и далее трансформируется в протяженный горизонтальный слой. Сделаем оценку, позволяющую показать, что при погружении холодной плиты она, достигнув границы между верхней и нижней мантией, не успевает существенно нагреться. Для этого приведем сравнение характерного времени погружения холодной плиты и ее разогрева до температуры окружающей среды:
погр. ' нагр. у ^ -1
'погр. Л
Здесь t , t - характерное время погружения ве-
^ погр/ нагр. 11 1 1 ^
щества холодной плиты на глубину верхней мантии и прогрева толщины плиты соответственно; Н , Н -
1 1 ' в.м. т.п.
толщина верхней мантии и толщина погружающейся плиты соответственно; V - скорость погружения пли-
погр. 1 1 ^
ты; х - температуропроводность горных пород.
Примем следующие характерные значения указанных величин:
Н =700км=7-105 м,
в.м.
Н =50-100 км~75 км=0.75-105 м,
т.п.
vпоr =5 см/год~1.6-10~9 м/с, Х=10-6-10-7 м2/с~0.5-10~6 м2/с. (8)
Подставляя величины (8) в формулы (7), получим:
0.039 ~ 4• 10-2.
t
нагр.
Отсюда следует, что для прогрева плиты до температуры окружающей среды необходимо более чем на
порядок больше времени и, следовательно, плита не успевает существенно нагреться, поэтому, благодаря быстрому погружению плиты, температуру на правой границе можно считать равной условному нулю.
Вопрос о том, как влияет субдукция на структуру конвекции в континентальной астеносфере, ранее рассматривался в работах [КМуаБИкт et а1., 2000, 2002; ЬоЬкоУБку et а1., 2014]. Один из выводов который следует из экспериментальной работы [КМуаБИкт et а1., 2002] и предыдущих исследований авторов этой работы, заключается в том, что в условиях горизонтального градиента температуры на подошве астеносферного слоя в нем могут образоваться протяженные «плоские» ячейки. В аналитических исследованиях горизонтально-протяженных конвективных ячеек в астеносфере число Рэлея имело значение 6.3-103. В работе [Кк-dyashkin et а1., 2002] приведены результаты экспериментального моделирования влияния субдукции на пространственную структуру конвективных течений в астеносфере под континентом при числах Рэлея в диапазоне Да=4.8-105-2.4-106, что может отвечать трехмерному режиму течений.
Рассмотрим несколько подробнее этот вопрос. Эксперименты показали, что в подогреваемом снизу горизонтальном слое при больших числах Прандтля область стационарной валиковой конвекции простирается до значений чисел Рэлея Да~2-104. Однако, как отмечается в работе [Kirdyashkin et а1., 2002], близкая к двумерной слабонестационарная конвекция может сохраняться до значений чисел Рэлея, в несколько раз превышающих данное значение Да~2-104. Таким образом, для наших целей, по крайней мере до чисел Рэлея Да~5-104, можно рассматривать двумерную задачу. С другой стороны, в нашей модели, когда учитывается горизонтальный градиент температуры и, в общем случае, горизонтальное движение основания слоя, можно предположить, что порог устойчивости двумерного течения существенно вырастет, поскольку указанные факторы есть элементы вынужденной конвекции. Строго говоря, это утверждение требует отдельного исследования, либо путем анализа сложной задачи устойчивости полученного здесь решения, либо проведением соответствующих экспериментов. Оценки числа Рэлея для астеносферы ввиду того, что параметры среды (особенно вязкость) известны весьма приближенно, могут разниться на два-три порядка. Если, например, взять следующие значения параметров:
Н =3-105 м, ДТ =300 °С,
в.м.
Х=5-10~7 м2/с, Д=10 5 1/°С,
v=5•(1015-1016) м2/с, то для числа Рэлея получим оценку Да~5-104-5-105. В настоящем исследовании за основное взято число Рэлея Да~5-104, хотя затрагиваются и другие его значения. При этом число Рэлея определено по толщине слоя и в качестве перепада температуры взят ее горизонтальный перепад.
Переходя к результатам, сначала рассмотрим случай, когда нижняя граница слоя находится в состоянии
покоя ^=0). На рис. 4 показаны результаты расчета линий тока при различных числах Рэлея. Из рис. 4, а, видно, что при Яа=103 внешние линии тока охватывают всю рассматриваемую область и соответствуют од-ноячеистой структуре. Однако внутри области имеются внутренние ячейки, вращающиеся в ту же сторону, что и внешняя ячейка. Из рис. 4, б, следует, что при йо=5-103 внешние линии тока по-прежнему соответствуют одноячеистой структуре. Это согласуется с экспериментом [Kirdyashkin et а1., 2002], где при йо=6.3-103 в астеносфере образуется всего одна ячейка. Однако в данном случае в центральной части области образуется внутренняя структура течений в виде ряда ячеек, вращающихся в ту же сторону, что и внешняя ячейка. Эта внутренняя структура сложнее, чем при меньшем числе Рэлея (рис. 4, а). Как показывает рис. 4, в, при большем числе Релея (йо=104) течение разбивается на ряд ячеек. При этом образуется одна протяженная ячейка с аспектным отношением, равным 6, и со сложным внутренним течением, а также несколько простых (без внутренней структуры) ячеек. При дальнейшем росте числа Рэлея число внутренних изометричных ячеек растет, а их аспектное отношение уменьшается.
Рассмотрим случай, когда учитывается движение нижней границы слоя, что соответствует горизонтальному движению литосферной плиты со скоростью v1 вдоль основания. При этом движение плиты нелинейным образом влияет на конвекцию, вызванную неоднородностью температуры. На рис. 5 показаны структуры течений при числе Рэлея йо=5-104 и различных скоростях движения плиты по основанию области. Исходный рис. 5, а, соответствует неподвижному основанию ^=0). В этом случае ячейка с максимальным ас-пектным отношением, примерно равным 2, образуется со стороны области субдукции. Неравномерность ячеек вызвана одновременным влиянием горизонтального и вертикального градиентов температуры. Рис. 5, б-д, соответствует случаям движения плиты вдоль основания с различной скоростью v1=(1.56-10.37) см/год. Из этого рисунка видно, что движение плиты переводит внутрь области течения большие ячейки, вращающиеся по часовой стрелке, т.е. по направлению циркуляции жидкости, вызванной плитой, и вытесняет меньшие ячейки, которые при отсутствии плиты вращались в противоположную сторону. На рис. 5, г, две из таких малых ячеек еще сохраняются, а при больших
(о)
(б)
1
-10 z
-10
(в)
z 'rv-2—~ 1 - 1 1Г ■
(У 0 © J 0 1С ^- 1-r-1- x 1 *
-10 z
-8
-6
-2
Рис. 4. Структура конвекции в области верхней мантии без учета движения плиты (v1=0) при числах Рэлея: (а) - Ra=103; (б) - Ra=5-103; (е) - Ra=104, (г) - Ra=5-104, (д) - Ra=5-105.
Fig. 4. Structure of upper mantle convection at different Rayleigh numbers: (а) - Ra=103; (б) - Ra=5-103; (в) - Ra=104, (г) - Ra=5-104, (д) - Ra=5-105. Plate movements are not taken into account, v1=0.
0
0
0
0
(а) 1
^ Л
О
0
a
-10
T
-8
Г
-6
-2
Рис. 5. Структура течений при йа=5-104 с учетом движения плиты вдоль основания со скоростями: (а) - v1=0; (б) - v1=300 (1.56 см/год); (е) - v1=500 (2.59 см/год); (г) - v1=1000 (5.18 см/год); (5) - v1=2000 (10.36 см/год).
Fig. 5. Structure of convection flows atfta=5-104 and different velocities of plate movements: (а) - v1=0; (б) - v1=300 (1.56 cm/yr); (е) - v1=500 (2.59 cm/yr); (г) - v1=1000 (5.18 cm/yr); (5) - v1=2000 (10.36 cm/yr).
x
0
0
скоростях плиты и они разрушаются и исчезают, как это видно из рис. 5, д.
Как показывает рис. 5, д, при больших скоростях движения плиты внешние линии тока охватывают всю рассматриваемую область, однако внутри нее образуется сложная структура течений, включающая ряд ячеек. Для больших чисел Рэлея, например Яа=105, картина качественно не отличается от представленной здесь. Однако, как отмечалось выше, при преодолении числом Рэлея некоторого порога, который в рассматриваемом случае наличия горизонтального градиента температуры и движения границы пока неизвестен, течение станет трехмерным и реальная картина течений может измениться.
Отметим существенное влияние эффекта движения плиты при характерных для верхней мантии числах Рэлея и скоростях субдукции литосферы, т.е. 5-10 см/год. Как показывает рис. 5, б-д, при таких скоростях и характерном для верхней мантии числе Рэлея во внешней области имеет место одноячеистая структура течений в согласии с представлениями, изложенными в работах [ЬоЪкоузку et а1., 2011, 2013; Laverov et а1., 2013; ЬоЪкоузку, 2016].
Существенные отличия полученных нами данных о конвекции в области верхней мантии, примыкающей к зоне субдукции, от результатов исследований аналогичной задачи, выполненных в работах [ЮМуа-зЬкт et а1., 2000, 2002], можно сформулировать следующим образом. Во-первых, в данной работе, наряду с влиянием горизонтального градиента температуры на развитие длинной горизонтальной ячейки, показана важная роль движущегося горизонтально лито-сферного слоя, подстилающего конвективную область выше фазовой границы 410 км, существование которого установлено по данным сейсмотомографии мантии в переходной зоне от Тихого океана к Восточной и Северо-Восточной Азии и Арктике. Во-вторых, в нашей более общей постановке задачи двухмерной конвекции получены режимы не только одноячеистой структуры течений в длинном горизонтальном слое, но и более сложные формы конвекции, сочетающие области длинных горизонтальных ячеек с изометричными зонами конвекции, что определяется небольшими изменениями значений числа Рэлея. В действительности, геолого-геофизические данные о распределении рифтовых зон и магматических провинций во времени
и пространстве, а также рассеянных деформациях и смещениях поверхности литосферы свидетельствуют о реализации обоих режимов конвекции в процессе эволюции Восточной и Северо-Восточной Азии и Арктики за последние 100 млн лет. На рис. 2 показана картина современных смещений поверхности литосферы Северо-Восточной Азии по данным GPS, очень четко демонстрирующая движение литосферы от БРЗ в сторону Японской островной дуги, которое естественным образом объясняется возвратным подлитосферным течением в длинной конвективной ячейке. Из рис. 2 также видно, что на фоне общего смещения литосферы в сторону Тихого океана в разных местах литосферы существуют зоны магматической активности различного возраста, являющиеся результатом развития
изометричных конвективных ячеек, которые можно трактовать как верхнемантийные плюмы.
6. 3D-МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНОЙ ВЕРХНЕМАНТИЙНОЙ КОНВЕКЦИИ, СОПРЯЖЕННОЙ С ЗОНОЙ СУБДУКЦИИ
В заключительном разделе представлены некоторые результаты 3D-моделирования верхнемантийной конвекции, инициируемой процессом субдукции, которое проводилось по описанной выше модели термической конвекции однородной вязкой жидкости в приближении Буссинеска (1). Моделирование было нацелено на изучение нестационарного трехмерного режима конвекции, возникающего при увеличении значения числа Рэлея, которое было взято равным 500000. Расчетная
Глубина 620 км (е) Глубина 40 км
Рис. 6. Распределение температур и продольных скоростей в разные моменты времени.
(а) - расстояния в плане и карты цветов; (б) - распределение T и Vx в 40 км от нижней границы; (е) - распределение T и Vx на глубине 40 км.
Fig. 6. Distribution of temperatures (T) and longitudinal velocities (Vx) at different time points.
(а) - distances in plan and colour codes; (б) - distribution T and Vx 40 km below the lower boundary; (е) - distribution Tand Vx at a depth of 40 km.
область имела форму прямоугольного бокса, который охватывал всю толщину верхней мантии. В данной постановке проводилось численное 3D-моделирование развития конвекции, включая горизонтальное движение литосферного слоя, поэтому область моделирования включала всю верхнюю мантию в отличие от вышеописанного 2D-моделирования, когда движение литосферы задавалось как граничное условие для расположенного выше слоя астеносферы.
Все границы бокса считались непроницаемыми, и задавалось отсутствие касательных напряжений на всех границах, кроме нижней, где учитывалось трение со стороны нижней мантии. Граничные условия для температуры выставлялись в традиционном для термической конвекции виде: постоянные температуры на горизонтальных границах - горячая (Г =1) снизу и холодная (Готн=0) сверху - и отсутствие тепловых потоков на всех вертикальных границах.
Принципиальным вопросом моделирования является выбор начального состояния региона. Для задания реалистичных начальных условий было проведено предварительное моделирование термической конвекции без процесса субдукции и полученное распределение температур и скоростей, соответствующее хаотичной квазистационарной ячейковой термической конвекции, было затем использовано в качестве начального ^0) состояния верхней мантии, при котором «включается» процесс субдукции. Так как модель термической конвекции (1) допускает множество решений, моделирование нуждается в задании дополнительного условия регуляризации решения [Kote1kin, Lobkovsky, 2014]. В качестве такого условия регуляризации решения задачи дополнительно задавалась и поддерживалась в течение всего времени моделирования низкая температура (Готн=0) в правом верхнем углу расчетной области размером около 200 км. Это локальное условие обеспечивало устойчивую последующую субдукцию холодного материала вдоль всей правой границы расчетной области. Результаты моделирования представлены на рис. 6.
На верхних панелях рис. 6, соответствующих начальному моменту времени, видна ячейковая конвекция, соответствующая подогреваемому снизу верхнемантийному слою. Видно также, что направления скоростей: синие участки (движение влево) и красные участки (движение вправо) - распределены случайно и равномерно, но при этом снизу скорости меньше, чем сверху, в силу наличия трения со стороны находящейся в покое нижней мантии.
Отметим, что каждое из четырех следующих по времени состояний, приводимых на рис. 6, было получено после 3000 шагов по времени, которые, в силу условия устойчивости, уменьшаются, когда скорость увеличивается. Вторая строка панелей отвечает времени 12.74 млн лет, показывая, что за это время холодное вещество успело опуститься сверху вниз и пройти такое же расстояние влево, о чем свидетельствуют отрицательные скорости (синий цвет) в нижней части
расчетной области. Одновременно в верхней части расчетной области (рис. 6, в) вещество движется вправо - положительные скорости красного цвета. Таким образом, в численном моделировании образуется примыкающая к субдукционному потоку горизонтально вытянутая конвективная ячейка.
Три последующих строки рис. 6, б, показывают, что холодный фронт на подошве области неуклонно движется влево и проходит более 2000 км. Распределение скоростей показывает, что синхронно с продвижением нижнего холодного фронта влево под континент увеличивается и длина конвективной ячейки, охватывая все большие области возвратного поверхностного течения вправо к зоне субдукции (красная область), при этом скорость фронта внизу постепенно уменьшается, о чем свидетельствует более бледный оттенок синего цвета.
Полученная нестационарная трехмерная картина конвекции показывает сосуществование двух видов течений: горизонтально-вытянутых и противоположно направленных потоков вещества снизу (синий цвет) и сверху (красный цвет), организующих общую циркуляцию, и многочисленных изометрических ячеек, которые можно трактовать как верхнеманийные плюмы.
Еще одно интересное явление, наблюдаемое в этом численном эксперименте на рис. 6, б, заключается в том, что в примыкающей к зоне субдукции конвективной ячейке развиваются вторичные валиковые течения.
7. ЗАКЛЮЧЕНИЕ
Получено решение задачи термической конвекции в верхней мантии в области, сопряженной с зоной суб-дукции. Из решения соответствующих уравнений следует, что при относительно малых числах Рэлея (Яа= =103-5-103) внешнее течение охватывает всю область, даже без учета движения плиты вдоль подошвы верхней мантии. При больших числах Рэлея, характерных для верхней мантии (Да=5-104), возникает ряд неравномерных ячеек. Однако горизонтальное движение погружающейся литосферной плиты вдоль основания с характерной для зон субдукции скоростью 7-10 см/год приводит к возникновению во внешней области линий тока, охватывающих всю область, т.е. к одноячеи-стой структуре внешних (периферийных) линий тока. При этом внутри области сохраняются ячейки, вращающиеся в том же направлении, что и внешняя ячейка, а ячейки, которые при отсутствии плиты вращались в противоположном направлении, вытесняются вверх и разрушаются. Проведено 3D-моделирование нестационарной конвекции при значении числа Рэ-лея = 500000. Полученные результаты позволяют объяснить весь спектр наблюдаемых тектономагматиче-ских процессов, развивающихся в пределах ЦВА в кайнозое и в Арктике в верхнем мелу и кайнозое, а именно сочетание общего смещения литосферы ЦВА и Арктики в сторону Тихоокеанской зоны субдукции с наличием отдельных магматических провинций и рифтовых зон как следствие существования длинной горизонтально-вытянутой конвективной ячейки (создающей
эффект конвейерного волочения литосферы), осложненной внутренними изометричными ячейками (создающими эффект верхнемантийных плюмов).
Таким образом, в рамках классической модели термической конвекции получено дополнительное математическое обоснование предложенной в работах [Lob-kovsky et al., 2011, 2013; Laverov et al., 2013; Lobkovsky, 2016] геодинамической модели мел-кайнозойской эволюции Арктики и Центрально-Восточной Азии, естественным образом объясняющей основные известные тектономагматические процессы, характерные для литосферы рассматриваемых областей.
8. БЛАГОДАРНОСТИ
Авторы выражают благодарность Ю.В. Габсатарову и И.С. Владимировой за помощь в построении карты смещений поверхности Земли под данным ГНСС. Авторы признательны О.П. Полянскому за ценные замечания, которые улучшили изложение ряда аспектов работы.
9. ЛИТЕРАТУРА / REFERENCES
Altamimi Z., Rebischung P., Metivier L., Collilieux X., 2016. ITRF2014: A New Release of the International Terrestrial Reference Frame Modeling Nonlinear Station Motions. Journal of Geophysical Research: Solid Earth 121 (8), 6109-6131. http://doi.org/10.1002/2016JB013098.
Apel E.V., Burgmann R., Steblov G., Vasilenko N., King R., Prytkov A., 2006. Independent Active Microplate Tectonics of Northeast Asia from Gps Velocities and Block Modeling. Geophysical Research Letters 33 (11). https://doi.org/10. 1029/2006GL026077.
Artyushkov E.V., Letnikov F.A., Ruzhich VV., 1990. On the Development of a New Mechanism for the Formation of the Baikal Basin. In: N.A. Logachev (Ed.), Geodynamics of Intra-continental Mountainous Regions. Nauka, Novosibirsk, p. 367-378 (in Russian) [Артюшков Е.В., Летников Ф.А., Ружич В.В. О разработке нового механизма формирования Байкальской впадины // Геодинамика внутри-континентальных горных областей / Ред. Н.А. Логачев. Новосибирск: Наука, 1990. С. 367-378].
Bott M., 1990. Geodynamic Processes in Continental Rift Systems as Applied to the Baikal Rift. In: N.A. Logachev (Ed.), Geodynamics of Intracontinental Mountainous Regions. Nauka, Novosibirsk, p. 317-323 (in Russian) [Ботт М. Геодинамические процессы в континентальных риф-товых системах в приложении к Байкальскому рифту // Геодинамика внутриконтинентальных горных областей / Ред. Н.А. Логачев. Новосибирск: Наука, 1990. С. 317-323].
Burgmann R., Kogan M.G., Steblov G.M., Hilley G., Levin VE., Apel E., 2005. Interseismic Coupling and Asperity Distribution along Kamchatka Subduction Zone. Journal of Geophysical Research: Solid Earth 110 (B7). http://dx.doi.org/ 10.1029/2005JB003648.
Chuvashova I.S., Rasskazov S.V., Yi-min Sun, 2017. The Latest Geodynamics in Central Asia: Primary and Secondary Mantle Melting Anomalies in the Context of Orogenesis,
Rifting, and Lithospheric Plate Motions and Interactions. Geodynamics & Tectonophysics 8 (1), 45-80 (in Russian) [Чувашова И.С., Рассказов С.В., Йи-минь Сунь. Новейшая геодинамика Центральной Азии: первичные и вторичные мантийные расплавные аномалии в контексте орогенеза, рифтогенеза и движения-взаимодействия литосферных плит // Геодинамика и тектонофизика. 2017. Т. 8. № 1. С. 45-80]. https://doi.org/10.5800/GT-2017-8-1-0232.
Djadkov P.G., Mel'nikova V.I., Nazarov L.A., Nazarova L.A., Sankov V.A., 1999. Increase of Seismotectonic Activity in the Baikal Region in 1989-1995: Results of Experimental Observations and Numerical Modeling of Changes in the Stress-Strained State. Russian Geology and Geophysics 40 (3), 373-386.
Florensov N.A., 1968. The Baikal Rift Zone and Some Problems of Its Study. In: Baikal Rift. Nauka, Moscow, p. 4056 (in Russian) [Флоренсов Н.А. Байкальская рифтовая зона и некоторые задачи ее изучения // Байкальский рифт. М.: Наука, 1968. С. 40-56].
Kirdyashkin A.A., Dobretsov N.L., Kirdyashkin A.G., 2002. Experimental Modeling of the Effect of Subduction on the Spatial Structure of Convective Flows in the Asthenosphere under the Continent. Doklady Earth Sciences 384 (5), 682686 (in Russian) [Кирдяшкин А.А., Добрецов Н.Л., Кирдя-шкин А.Г. Экспериментальное моделирование влияния субдукции на пространственную структуру конвективных течений в астеносфере под континентом // Доклады Академии наук. 2002. Т. 384. № 5. С. 682-686].
Kirdyashkin A.A., Kirdyashkin A.G., Dobretsov N.L., 2000. Influence of Subduction on the Structure of Thermal Gravitational Flows in the Asthenosphere under the Continent. Russian Geology and Geophysics 41 (2), 207-219 (in Russian) [Кирдяшкин А.А., Кирдяшкин А.Г., Добрецов Н.Л. Влияние субдукции на структуру тепловых гравитационных течений в астеносфере под континентом // Геология и геофизика. 2000. Т. 41. № 2. С. 207-219].
Kirdyashkin A.G., Dobretsov N.L., 1991. Modelling of Two-Layer Mantle Convection. Doklady of the USSR Academy of Sciences 318 (4), 946-949 (in Russian) [Кирдяшкин А.Г., Добрецов Н.Л. Моделирование двухслойной мантийной конвекции // Доклады АН СССР. 1991. T. 318. № 4. C. 946-949].
Kogan M.G., Burgmann R., Vasilenko N.F., Scholz C.H., King R.W, Ivashchenko A.I., Frolov D.I., Steblov G.M., Kim Ch.U., Egorov S.G., 2003. The 2000 Mw 6.8 Uglegorsk Earthquake and Regional Plate Boundary Deformation of Sakhalin from Geodetic Data. Geophysical Research Letters 30 (3), 1-4. http://doi.org/10.1029/2002GL016399.
Kotelkin V.D., Lobkovsky L.I., 2014. Regularization of Geodynamic Problems Using Geological Data. Fluid Dynamics 49, 310-319. http://doi.org/10.1134/S00154628 14030028.
Kotelkin V.D., Lobkovsky L.I., 2019. Substantiation of a Geodynamic Model of the Evolution of the Arctic Region. In: Materials of the XII All-Russia Congress on Fundamental Problems of Theoretical and Applied Mechanics (August 1924, 2019). Vol. 4. Bashkir State University, Ufa, p. 63-65 (in
Russian) [Котелкин В.Д., Лобковский Л.И. Обоснование геодинамической модели эволюции Арктического региона // Материалы XII Всероссийского съезда по фундаментальным проблемам теоретической и прикладной механики (19-24 августа 2019 г.). Уфа: РИЦ БашГУ 2019. Т. 4. С. 63-65].
Laverov N.P., Lobkovsky L.I., Kononov M.V, Dobretsov N.L., Vernikovsky V.A., Sokolov S.D., Shipilov E.V., 2013. A Geody-namic Model of the Evolution of the Arctic Basin an Adjacent Territories in the Mesozoic and Cenozoic and the Outer Limit of the Russian Continental Shelf. Geotectonics 1, 1-30. http://doi.org/10.1134/S0016852113010044.
Lobkovsky L.I., 2016. Tectonics of Deformable Lithosphere Plates and Regional Geodynamic Model of the Arctic and Northeastern Asia. Russian Geology and Geophysics 57 (3), 371-386. http://doi.org/10.1016/j.rgg.2016.03.002.
Lobkovsky L.I., Inyukhin A.V., Kotelkin V.D., 2014. Subduction and Cyclicity of the Upper Mantle Processes. Doklady Earth Sciences 459, 1348-1352. https://doi.org/10.1134/ S1028334X14110233.
Lobkovsky L.I., Ramazanov M.M., 2021. Investigation of Upper Mantle Convection Thermomechanically Related to a Subduction Zone, and Its Geodynamic Applications for the Arctic and Northeast Asia. Fluid Dynamics 3, 139150 (in Russian) [Лобковский Л.И., Рамазанов М.М. Исследование конвекции в верхней мантии, термомеха-нически связанной с зоной субдукции, и ее геодинамические приложения для Арктики и Северо-Восточной Азии // Известия РАН: Механика жидкости и газа. 2021. № 3. С. 139-150]. https://doi.org/10.31857/S056852812 1030063.
Lobkovsky L.I., Shipilov E.V, Kononov M.V, 2013. Geodynamic Model of Upper Mantle Convection and Transformations of the Arctic Lithosphere in the Mesozoic and Cenozoic. Izvestiya of of the Physics of the Solid Earth 49, 767-785. http://doi.org/10.1134/S1069351313060104.
Lobkovsky L.I., Verzhbitsky V.E., Kononov M.V, Shreider AA., Garagash I.A., Sokolov S.D., Tuchkova M.I., Kotelkin VD., Vernikovsky V.A., 2011. Geodynamic Model of the Evolution of the Arctic Region in the Late Mesozoic - Cenozoic and the Problem of the Outer Boundary of the Continental Shelf of Russia. Arctic: Ecology and Economics 1 (1), 104-115 (in Russian) [Лобковский Л.И., Вержбицкий В.Е., Кононов М.В., Шрейдер А.А., Гарагаш И.А., Соколов С.Д., Тучкова М.И., Котелкин В.Д., Верниковский В.А. Геодинамическая модель эволюции Арктического региона в позднем мезозое - кайнозое и проблема внешней границы континентального шельфа России // Арктика: экология и экономика. 2011. № 1 (1). С. 104-115].
Lobkovsky L.I., Vladimirova I.S., Gabsatarov V.V., Garagash I.A., Baranov B.V, Steblov G.M., 2017. Post-Seismic Motions after the 2006-2007 Simushir Earthquakes at Different Stages of the Seismic Cycle. Doklady Earth Sciences 473, 375-379. http://doi.org/10.1134/S1028334X17 030266.
Lobkovsky L.I., Vladimirova I.S., Gabsatarov Y.V., Steblov G.M., 2018. Seismotectonic Deformations Related to the 2011 Tohoku Earthquake at Different Stages of the Seismic
Cycle on the Basis of Satellite Geodetic Observations. Doklady Earth Sciences 481, 1060-1065. http://doi.org/10. 1134/S1028334X18080159.
Logatchev N.A., 1968. Sedimentary and Volcanogenic Formations of the Baikal Rift Zone. In: N.A. Florensov (Ed.), Baikal Rift. Nauka, Novosibirsk, p. 72-101 (in Russian) [Логачев Н.А. Осадочные и вулканогенные формации Байкальской рифтовой зоны // Байкальский рифт / Ред. Н.А. Флоренсов. Новосибирск: Наука, 1968. С. 72-101].
Mats V.D., 2015. The Baikal Rift: Pliocene (Miocene) -Quarternary Episode or Product of Extended Development since the Late Cretaceous under Various Tectonic Factors. A Review. Geodynamics & Tectonophysics 6 (4), 467-490 (in Russian) [Мац В.Д. Байкальский рифт: плиоцен (миоцен) - четвертичный эпизод или продукт длительного развития с позднего мела под воздействием различных тектонических факторов. Обзор представлений // Геодинамика и тектонофизика. 2015. Т. 6. № 4. С. 467490]. https://doi.org/10.5800/GT-2015-6-4-0190.
Meng G., Shen X., Wu J., Rogozhin E.A., 2006. Present-Day Crustal Motion in Northeast China Determined from GPS Measurements. Earth, Planets and Space 58, 1441-1445. https://doi.org/10.1186/BF03352642.
Molnar P., Tapponnier P., 1975. Cenozoic Tectonics of Asia: Effects of a Continental Collision: Features of Recent Continental Tectonics in Asia Can Be Interpreted as Results of the India-Eurasia Collision. Science 189 (4201), 419426. https://doi.org/10.1126/science.189.4201.419.
Mordvinova V.V., Treusov A.V., Sharova E.V., Grebenshchi-kova VI., 2008. Results of Teleseismic Two-Dimensional P-To-mography: Evidence of the Mantle Plume under Khangai. In: Geodynamic Evolution of the Lithosphere of the Central Asian Mobile Belt (from Ocean to Continent). Proceedings of Scientific Meeting (October 14-18, 2008). Iss. 6. Vol. 2. IEC SB RAS, Irkutsk, p. 41-43 (in Russian) [Мордвинова В.В., Треусов А.В., Шарова Е.В., Гребенщикова В.И. Результаты телесейсмической двумерной Р-томографии: свидетельство мантийного плюма под Хангаем // Геодинамическая эволюция литосферы Центрально-Азиатского подвижного пояса (от океана к континенту): Материалы научного совещания по Программе фундаментальных исследований ОНЗ РАН (14-18 октября 2008 г.). Иркутск: ИЗК СО РАН, 2008. Вып. 6. Т. 2. С. 41-43].
Rasskazov S.V., Chuvashova I.S., 2013. The Most Recent Mantle Geodynamics of Central Asia. ISU Publishing House, Irkutsk, 308 p. (in Russian) [Рассказов С.В., Чувашова И.С. Новейшая мантийная геодинамика Центральной Азии. Иркутск: Изд-во ИГУ, 2013. 308 с.].
Ruzhich V.V., 1997. Seismotectonic Destruction of the Earth's Crust in the Baikal Rift Zone. Publishing House of SB RAS, Novosibirsk, 144 p. (in Russian) [Ружич В.В. Сейсмотектоническая деструкция в земной коре Байкальской рифтовой зоны. Новосибирск: Изд-во СО РАН, 1997. 144 с.].
Ruzhich V.V., Kocharyan G.G., Levina Е.А., 2016. Estimated Geodynamic Impact from Zones of Collision and Subduction on the Seismotectonic Regime in the Baikal Rift. Geodynamics & Tectonophysics 7 (3), 383-406 (in Russian)
[Ружич В.В., Кочарян Г.Г., Левина Е.А. Оценка геодинамического влияния зон коллизии и субдукции на сейсмотектонический режим Байкальского рифта // Геодинамика и тектонофизика. 2016. Т. 7. № 3. С. 383-406]. http://doi.org/10.5800/GT-2016-7-3-0214.
Sankov V.A., Parfeevets A.V., Lukhnev A.V., Miroshni-chenko A.I., Ashurkov S.V., 2011. Late Cenozoic Geodynamics and Mechanical Coupling of Crustal and Upper Mantle Deformations in the Mongolia-Siberia Mobile Area. Geo-tectonics 45, 378-393. http://dx.doi.org/10.1134/S0016 852111050049.
Seminsky K.Z., Kozhevnikov N.O., Cheremnykh A.V., Po-speeva E.V., Bobrov A.A., Olenchenko V.V., Tugarina M.A., Potapov V.V., Zaripov R.M., Cheremnykh A.S., 2013. Interblock Zones in the Crust of the Southern Regions of East Siberia: Tectonophysical Interpretation of Geological and Geophysical Data. Geodynamics & Tectonophysics 4 (3), 203-278 (in Russian) [Семинский К.Ж., Кожевников Н.О., Черемных А.В., Поспеева Е.В., Бобров А.А., Оленчен-ко В.В., Тугарина М.А., Потапов В.В., Зарипов Р.М., Черемных А.С. Межблоковые зоны в земной коре юга Восточной Сибири: тектонофизическая интерпретация геолого-геофизических данных // Геодинамика и тек-тонофизика. 2013. Т. 4. № 3. С. 203-278]. https://doi. org/10.5800/GT-2013-4-3-0099.
Steblov G.M., Kogan M.G., King R.W., Scholz C.H., Burgmann R., Frolov D.I., 2003. Imprint of the North American Plate in Siberia Revealed by GPS. Geophysical Research Letters 30 (18). http://doi.org/10.1029/2003GL017805.
Yarmolyuk V.V., Kovalenko I.I., Ivanov V.G., 1995. Intraplate Late Mesozoic-Cenozoic Volcanic Province of Central-East Asia - Projection of the Hot Mantle Field. Russian Geotec-tonics 5, 41-67 (in Russian) [Ярмолюк В.В., Коваленко И.И., Иванов В.Г. Внутриплитная позднемезозойская - кайнозойская вулканическая провинция Центрально-Восточной Азии - проекция горячего поля мантии // Геотектоника. 1995. № 5. С. 41-67].
Yarmolyuk V.V., Kuzmin M.I., Vorontsov A.A., 2013. West Pacific-Type Convergent Boundaries and Their Role in the Formation of the Central Asian Fold Belt. Russian Geology and Geophysics 54 (12), 1427-1441. https://doi.org/10. 1016/j.rgg.2013.10.012.
Zhao D., 2009. Multiscale Seismic Tomography and Mantle Dynamics. Gondwana Research 15 (3-4), 297-323. https:// doi.org/10.1016/j.gr.2008.07.003.
Zhao D., Lei J., Inoue T., Yamada A., Gao S.S., 2006. Deep Structure and Origin of the Baikal Rift Zone. Earth and Planetary Letters 243 (3-4), 681-691. https://doi.org/10.1016/j. epsl.2006.01.033.
Zhao D., Liu L., Pirajno F., Dobretsov N.L., 2010. Mantle Structure and Dynamics under East Russia and Adjacent Regions. Russian Geology and Geophysics 51 (9), 925-938. https://doi.org/10.1016/j.rgg.2010.08.003.
Zhao D., Tian Y., Ley J., Liu L., Zheng S., 2009. Seismic Image and Origin of the Changbai Intraplate Volcano in East Asia: Role of Big Mantle Wedge above the Stagnant Pacific Slab. Physics of the Earth and Planetary Interiors 173 (3-4), 197-206. https://doi.org/10.1016/j.pepi.2008.11.009.
Zonenshain L.P., Savostin L.A., Misharina L.A., Solonen-ko N.V., 1978. Plate Tectonics of the Baikal Mountain Region and the Stanovoy Ridge. Doklady of the USSR Academy of Sciences 240 (3), 669-672 (in Russian) [Зоненшайн Л.П., Савостин Л.А., Мишарина Л.А., Солоненко Н.В. Тектоника плит Байкальской горной области и Станового хребта // Доклады АН СССР. 1978. Т 240. № 3. С. 669-672].
Zorin Y.A., Turutanov E.X., 2005. Plumes and Geodynamics of the Baikal Rift Zone. Russian Geology and Geophysics 46 (7), 685-699 (in Russian) [Зорин Ю.А., Турутанов Е.Х. Плюмы и геодинамика Байкальской рифтовой зоны // Геология и геофизика. 2005. Т. 46 № 7. С. 685-699].
Zorin YA., Turutanov E.K., Kozhevnikov VM., Rasskazov S.V, Ivanov A.V., 2006. The Nature of Cenozoic Upper Mantle Plumes in East Siberia (Russia) and Central Mongolia. Russian Geology and Geophysics 47 (10), 1056-1070.
