Трехмерное моделирование конвекции под кратонами центральной Азии* Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

Вычислительные технологии

Том 12, Специальный выпуск 4, 2007

ТРЕХМЕРНОЕ МОДЕЛИРОВАНИЕ КОНВЕКЦИИ ПОД КРАТОНАМИ ЦЕНТРАЛЬНОЙ АЗИИ*

, В. В. Червов Институт нефтегазовой геологии и геофизики СО РАН,

Новосибирск, Россия e-mail: chervov@ict.nsc.ru Г. Г. Черных

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: chernykh@ict.nsc.ru

In this paper we present the results of the 3D modeling of convection under cratons of Central Asia. Numerical model is based on the "vector potential — vorticity" variables and the method of fractional steps. Computational results that demonstrate the structure of convective fluxes are given.

Введение

Данные о современной структуре мантии, полученные методами сейсмической томографии, дают представление о пространственном распределении сейсмических неоднород-ностей во всей мантии Земли. При численном моделировании конвективных процессов оказывается весьма полезным знание реологии литосферных блоков, полученное в результате анализа глобальных сейсмотомографических моделей. Моделирование конвективных процессов в верхней мантии Земли позволяет не только уточнить имеющиеся на настоящий момент знания о состоянии мантии, но и проследить динамику глубинных процессов во времени. Численному моделированию трехмерных конвективных течений в мантии Земли посвящено большое количество работ (см., например, [1-11] и приведенную в них библиографию). В [8-10] построена и детально протестирована численная модель тепловой конвекции в мантии Земли, основанная на переменных завихренность — векторный потенциал, методе дробных шагов, последовательности сеток и экстраполяции Ричардсона. Основанная на неявном методе расщепления по физическим процессам трехмерная численная модель конвекции в верхней мантии Земли предложена в [11].

В настоящей работе выполнено численное моделирование тепловой конвекции под внутриконтинентальной областью Азии, в которую входят Западно-Сибирская плита, Сибирская платформа, Центрально-Азиатский складчатый пояс, Тарим и часть СевероКитайской платформы. Структура континентальной литосферы исследуемой области

* Работа выполнена при финансовой поддержке интеграционного проекта СО РАН № 116.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

С. А. Тычков

весьма неоднородна. Значения толщины литосферы древних платформ, таких как Сибирская, Тарим и Китайская, составляют 200,,, 250 км, в то время как для палеозойской Западно-Сибирской плиты значения толщины литосферы не превышают 120,,, 130 км, В рифтовых долинах имеет место утонение литосферы до 40 км. Подобные вариации мощности литосферы существенно влияют на характер мантийных течений [12-14], и поэтому их учет имеет первостепенное значение при проведении численного моделирования.

1. Постановка задачи

Для описания течения в верхней мантии Земли привлекается хорошо известная математическая модель, включающая в себя обезразмеренные уравнения [13]:

ди ду дю

7Г + 7Г + 7Г= 0; дх ду дх

др ^ д ди ^ д дх дх ^ дх д]Р др д /ди ду ду дх^ \ду дх др д (ди дю дх дх^

ди ду ду дх

д

+ 2 ® дь д ду ^ ду дх^

ди д т

дх дх

ду д т

дх ду

д /ду дю\ ^д дъи ^ ^ дх дх ) ду ^ \ дх ду ) дх ^ дх

дТ дТ дТ дТ „~ — + и— + у— + ю— = У2Т. дЬ дх ду дх

(1) (2)

(3)

(4)

(5)

Здесь и, у, ю — компоненты вектора скорости; р — давление; Т — температура;

Яа

адх рАТй3

ПХ

— число Рэлея; а — температурный коэффициент расширения; дх —

х АТ

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

Система уравнений (1)-(5) устроена так, что в начальный момент времени Ь = ¿0 задаются начальные условия лишь для температуры: Т(х, у, х, ¿0) = Т0(х, у, х).

Задача решалась в параллелепипеде 0 < х < X, 0 < у < У, 0 < х < Z. В качестве краевых условий на боковых границах задаются условия симметрии, а на нижней и верхней — условия прилипания и фиксированные значения температуры. На границах неоднородной литосферной плиты также задаются условия прилипания как в вертикальном, так и в латеральном направлении. На нижней кромке литосферы при постановке начального распределения температуры учитывается ее первоначальное значение: Т = 1200 °С, Температура рассчитывалась во всем параллелепипеде, движение жидкости — вне литосферы.

2. Результаты расчетов

В качестве моделей неоднородной литосферы в [15] рассматривались к ратины — толстые и прямоугольные в плане участки кондуктивной литосферы с длиной граней от

500 до 5000 км, окруженные литосферой с нормальной толщиной, В настоящей работе учтены результаты моделирования простейших кратонов и приведены результаты моделирования ловушек — утоненных участков литосферы. Основные результаты получены на основе математической модели в декартовых координатах на сетке 67 х 67 х 36 с применением переменных завихренность — векторный потенциал [8-10], Число Рэлея, характеризующее режим конвекции, выбрано как Ra = 3.04 ■ 105, что отвечает современным представлениям об условиях в недрах Земли, Основные параметры задачи в системе СИ, пригодные для верхней мантии, выбирались следующими: d = 700 000 м, AT = 1800 °С, х = 10-6 м2/с, а = 10-5 °С-1, р = 3300 кг/м3, gz = 10 м/с2, По = 3 ■ 1021 кг/(м-с).

Зависимость вязкости от температуры и глубины выражена следующей формулой: n(x,y,z,t) = ebz-aT(x'y'z't\ Здесь параметры а = 3.89 и b = 5.84 обеспечивают перепад вязкости от 20 до 200, что присуще верхнемантийным характеристикам течений,

2.1. Эволюция конвективных течений под кратономи

В качестве первой модельной задачи рассмотрена задача о влиянии утолщенных участков литосферы (кратонов) на движение мантийного вещества [15], Вычисления проводились в параллелепипеде П = [0, 4200 км] х [0, 4200 км] х [0, 700 км], на равномерных сетках 67 х 67 х 36 и 101 х 101 х 54 ячеек. Рассматривались задачи с утолщенной полосой, с кратоном 2000 х 2000 км и с двумя кратонами — 1500 х 1000 км и 1000 х 1000 км. Мощность литосферы составляла 120 км, мощность кратонов — 220 км. Характерной чертой пространственно-временной эволюции конвективных течений под кратонами является

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

°

°

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

Таким образом, к одному из принципиальных результатов трехмерного моделирования тепловой конвекции в верхней мантии под литосферой переменной мощности можно отнести обнаруженную мелкомасштабную моду конвекции непосредственно под литосферой на "астеносферном" уровне глубин 200,,, 350 км. Мелкомасштабная конвекция существует в виде вытянутых ячеек с горизонтальным размером в 500 км между восходящим и нисходящим потоками конвекции. Данная мода развивается по периферии утолщенного участка литосферы и может объяснить особенности режима траппового магматизма древних кратонов и их периферии [15],

2.2. Тепловая конвекция под литоеферной плитой

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

Геофизические наблюдения [13, 18] обнаруживают существенные латеральные вариации толщины континентальной литосферы. Если толщина литосферы докембрийских платформ составляет 200,,, 240 км, то окружающие их более молодые складчатые пояса

характеризуются более значительными вариациями мощности — от 50 до 120 км в зависимости от истории формирования и современного тектонического режима. Участки с утоненной континентальной литосферой часто маркируют древние зоны активных континентальных окраин, как это предполагается, например, для района Байкальского рифта. Другим примером тонкой литосферы в рассматриваемой области служит район Великих Озер в Западной Монголии, где участок с утоненной литосферой шириной около 200 км простирается в субмеридиональном направлении на юг от оз. Убсу-Нур, Выше отмечалось, что утолщение литосферы существенно влияет на структуру и режим тепловой верхнемантийной конвекции под континентами. Поэтому, как представляется, логично исследовать также влияние и утоненных участков литосферы (ловушек) на конвекцию в мантии, С целью анализа конвективного течения иод утоненным участком литосферы выполнены численные эксперименты. Расчеты проводились в параллелепипеде П = [0, 4200 км] х [0, 4200 км] х [0, 700 км], на равномерных сетках 66 х 66 х 35 и 80 х 80 х 50 ячеек. Рассматривалась задача, в которой мощность литосферы составляла 200 км, мощность литосферы в области ловушки 120 км.

Области с тонкой литосферой обнаружены практически на всех континентах и включают крупные осадочные бассейны, протяженные сдвиговые зоны но периферии кра-тонов, иостколлизионные швы. Расположение утоненной полосы и состояние конвергирующей мантии на глубине 160 км показаны на рис, 1 (другое расположение утоненной полосы рассмотрено в |19|), Влияние тонкого участка литосферы, как показывают результаты расчетов, начинает сказываться на структуре течения уже через десятки миллионов лет (характерный шаг сетки но времени принимался равным 6,25 млн лет). Общая тенденция изменения структуры течения состоит в формировании нисходящих

У

400

300

250

200

350

100

150

50

км

X

1800 °С 1700 °С 1600 °С 1500 °С 1400 °С 1300 °С 1200 °С 1100 °С 1000 °С 900 °С5: 800 °С5: 700 °С5: 600 °С5: 500 °С5: 400 °С5: 300

200':°С:. 100':°С:.

0 И

0 500 1000 1500 2000 2500 3000 ~3500 4000 км

Рис. 1. Сечение в плоскости (ХУ) на глубине 160 км в модели конвекции под литосферой с ловушкой при Ь = 2500 млн лет. Мощность литосферы Ь = 200 км, мощность полосы 120 км

Y

t

Ч'800 °с

=¡1700 Щ =11воо °е

=¡1500 *Ц = 14оо = 1эоо щ

— 1200 °С

— 1100 °С = 1000 °с

= 900 °С = SOU °с I 700 Щ

— 500 щ

— 500 °с

— 400 °с

— 300 °с

= 200 °с

— 100 °с

■а сс-

О 500 1000 1500 2000 2500 3000 3500 4000км

Рис. 2. Сечение в плоскости (XY) на глубине 220 км в модели конвекции под литосферой с ловушкой при t = 2500 млн лет. Мощность литосферы L = 200 км, мощность полосы 120 км

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

В результате численных экспериментов показано, что нисходящие потоки но прошествии 625 млн лет постепенно "выдавливают" восходящие в область с толстой литосферой. Структура течения приходит в квазистациопарное состояние через 2500 млн лет. На глубине более 380 км существует охлажденная область с температурой пе более 1050.,, 1100 °С под всей полосой. Выше, на "астеносферных" глубинах (120.,, 350 км), как и в предыдущих моделях, в приграничных областях или областях с высоким градиентом изменения толщины кондуктивной литосферы развивается мелкомасштабная мода конвекции, что хорошо видно па рис. 1 и 2.

°

°

но в плоскости, параллельной границе области тонкой литосферы. Размер этих ячеек составляет 500 км между восходящим и нисходящими потоками. Структура литосферы в данной модели позволяет снова вернуться к обсуждению вопроса о правомерности применения двухмерного моделирования для изучения динамики мантии под вытянутыми геологическими структурами. С одной стороны, под всей полосой в целом фиксируются нисходящие конвективные течения, что было получено и в двухмерных моделях [17|, Но вместе с тем трехмерные модели выявили мелкомасштабную моду конвекции в приграничных районах, развивающуюся вдоль полосы. Это может оказаться принципиальным моментом при сопоставлении результатов моделирования с реальными геологическими структурами (см. также |19|).

3. Результаты моделирования тектонических зон Центральной Азии

В Центральной Азии имеется ряд платформенных областей, среди которых можно выделить Таримскую плиту, Северо-Китайский и Южно-Китайский кратоны. С юга к области Центральной Азии примыкает Индийская платформа, играющая важную роль в формировании современной тектонической обстановки региона. В северной части область включает Западно-Сибирскую палеозойскую плиту и древнюю Сибирскую платформу. Эти области оказывают существенное влияние на стиль деформирования и тектонический режим литосферы Центральной Азии [20].

В настоящей работе моделирование процессов в верхней мантии было ограничено внутриконтинентальной областью Азии (рис. 3), в которую вошли Сибирская платформа, расположенная восточнее р. Енисей и простирающаяся в этом направлении до гор Верхоянья, у подножия которых течет р. Лена. На юге эта платформа ограничена оз. Байкал, а на севере — Енисей-Хатангской низменностью. Западно-Сибирская плита, примыкающаяя с запада к Сибирской платформе, Центрально-Азиатский складчатый пояс, Тарим и часть Северо-Китайской платформы расположены южнее. Вычисления проводились в параллелепипеде П = [0, 4200 км] х [0, 4200 км] х [0, 700 км], на равномерных сетках 66 х 66 х 35 и 100 х 100 х 53 ячеек; величина шага по времени 5 и 2.5 млн лет. Отличие решений при этом не превосходило 12% в равномерной сеточной норме.

Чтобы проиллюстрировать положение кратонов и ловушки в пространстве, приводится рис. 4, на котором можно увидеть конфигурацию ловушки и соотношения мощностей литосферных блоков в исследуемом районе. Схема на рис. 5 показывает расположение элементов литосферы в расчетной области. Кратоны при моделировании выбирались по размерам несколько меньше одноименных платформ. На рис. 6. изображена область моделирования конвекции в верхней мантии под Центральной Азией. Величина и расположение кратонов приближенно соответствуют реальным данным. Моделирование показало, что, как и в случае прямоугольных в плане кратонов [15], Сибирский кратон порождает аналогичные структуры, а именно устойчивые восходящие потоки

Рис. 3. Рельеф внутриконтинентальной области Азии: черная сплошная линия — границы Сибирской платформы, Тарима и Северо-Китайского кратона

Север

запад

Восток

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

Рис. 5. Схема расположения кратонов и ловушки с указанием мощностей литосферных элементов. Вне кратонов и ловушки толщина литосферы 120 км

1800

1700 с

1600

1500 "С

1400'С

1300

1200 '"С

1100 °с

1000

900 'С

800 С

700

600 с

00-

500 С

400

300 С

200 ÙC

100

О 500 1000 1500 2000 2500 3000 3500 4000

Рис. 6. Сечение температурного поля в плоскости (ХУ) на глубине 317 км в модели конвекции под литосферой Центральной Азии с тремя кратонами. Между центральным и южными кратонами ловушка, мощность литосферы над которой 60 км. Мощность литосферы, не занятой кратонами и ловушкой, 120 км

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

Все наблюдаемые на поверхности особенности тепловых полой достаточно хорошо прослеживаются на рассчитанных тепловых полях. Например, совпадают положения центрального плюма и восходящего потока в северо-западной части кратопа, ответственные за трапповый магматизм, который имел место в прошлом Сибирской платформы в интервале 240... 250 млн лет назад, т. е. с середины пермского периода и до начала триасового, когда сам кратон проходил над нижнемантийным шпомом. Сложение температур верхне- и нижнемантийного шпомов привело к излиянию платобазаль-тов па западе Сибирской платформы. Участие пижпемаптийпого вещества в излияниях сейчас подтверждено геохимическими и изотопными исследованиями лав Сибири [16|, На рис. 7 показано сечение температурного поля па глубине 251 км. Видно, что северозападный плюм под кратопом (см. рис. 6) изменил направление, обтекая подошву кратопа, породив два восходящих потока около северо-западного угла Сибирской платформы. Такая же "участь" постигла и другие восходящие непосредственно под кратопом

Рис. 7. Сечение температурного поля в плоскости (ХУ) на глубине 251 км в модели конвекции под литосферой Центральной Азии с тремя кратонами

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

Представленная здесь трехмерная численная модель конвекции также обнаруживает повышение средней мантийной температуры под кратоном на 100 °С, но вместе с тем показывает более сложные формы рельефа кратопа, обусловленные динамическим воздействием конвекции. Для сравнения рельефа древнего кратопа Сибирской платформы (см. рис. 3) с результатами вычислений воспользуемся картой распределения температуры в литосфере кратопа на глубине 251 км (см. рис. 7). Как было показано в [14|, более высокие температуры в литосфере соответствуют приподнятым участкам литосферы, а пониженная температура — относительно онущеппым участкам поверхности. Из сопоставления рельефа платформы с нолем температур видно, что на платформе, как и в модели, есть два региональных поднятия — вдоль западной границы Сибирской платформы с юга на север вплоть до плато Путорана на северо-западе и Патомское нагорье, Алданский щит па юго-востоке и далее к Верхояпыо до устья Лены на северо-востоке Сибирского кратоиа. Эти два поднятия разделены вытянутой в центральной части низ-

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

По геолого-геофизическим данным в районе южнее Сибирского кратона и севернее Тарима и Северо-Китайской платформы мощность литосферы составляет от 40 до 75 км, В численной модели толщина литосферного блока в указанном районе принималась равной 60 км, В результате численного моделирования было показано, что в зоне ловушки преобладают нисходящие потоки холодного мантийного материала, И в конкретной геологической обстановке, а именно в случае взаимодействия трех кратонов, в

самой ловушке также наблюдается цепь классического нисходящего потока. На глубине

°

рис, 6), Наиболее прогретыми оказываются места, расположенные южнее и восточнее Сибирского кратона, и узкая полоса вблизи (севернее) Таримского и Северо-Китайского кратонов, В районе оз. Байкал, в области ловушки (мощность литосферы 60 км), в непосредственной близости от Сибирского кратона (мощность литосферы 250 км), наблюдается тепловая аномалия в виде мелкомасштабной конвективной ячейки, которая имеет вытянутую форму, напоминающую Байкал,

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

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

Список литературы

[1] Рыков В.В., Трубицин В.П. Численное моделирование трехмерной мантийной конвекции и тектоника литосферных плит // Вычисл. сейсмология. 1994. Вып. 26. С. 94-102.

[2] Рыков В.В., Трубицин В.П. Трехмерная модель мантийной конвекции с движущимися континентами // Вычисл. сейсмология. 1994. Вып. 27. С. 21-41.

[3] Gable W., O'Connell J. Convection in three dimensions with surface plates: generation of toroidal flow //J. Geophys. Research. 1991. Vol. 96. N B5. P. 8391-8405.

[4] Busse F.H., Christensen U., Clever R. et al. 3D Convection at infinite Prandtl number in cartesian geometry — a benchmark comparison // Geophys. Astrophys. Fluid Dynamics. 1993. Vol. 75. P. 39-59.

[5] Glatzmaier G.A., Schubert G. Three-dimensional spherical model of layered and whole mantle convection //J. Geophys. Research. 1993. Vol. 98, N B12. P. 21969-21976.

[6] Machetel P., Thoravan C., Brunet D. Spectral and geophysical consequences of 3D spherical mantle convection with an endothermic phase change at the 670 km discontinuity // Phys. Earth and Planetary Interiors. 1995. Vol. 88. P. 43-51.

[7] Zhong S., Zuber M. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection //J. Geophys. Research. 2000. Vol. 105, N B5. P. 11063-11082.

[8] Червов B.B. Численное моделирование трехмерных задач конвекции в мантии Земли с применением завихренности и векторного потенциала // Вычисл. технологии. 2002. Т. 7, № 1. С. 114-125.

[9] Червов В.В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением последовательности сеток // Вычисл. технологии. 2002. Т. 7, № 3. С. 85-92.

[10] Тычков С.А., Червов В.В., Черных Г.Г. Численная модель трехмерной конвекции в мантии Земли // Вычисл. технологии. 2006. Т. 11. Ч. 2. С. 45-52.

[11] Червов В.В. Моделирование трехмерной конвекции в мантии Земли с применением неявного метода расщепления по физическим процессам // Вычисл. технологии. 2006. Т. 11, № 4. С. 73-86.

[12] Добрецов H.J1. Пермотриасовый магматизм и осадконакопление в Евразии как отражение суперплюма // Докл. РАН. 1997. Т. 354. С. 220-223.

[13] Добрецов H.JL, Кирдяшкин А.Г., Кирдяшкин A.A. Глубинная геодинамика. 2-е изд. Новосибирск: Изд-во СО РАН, 2001. 409 с.

[14] Тычков С.А., Рычкова Е.В., Василевский А.Н., Червов В.В. Тепловая конвекция в верхней мантии континентов и ее эффект в геофизических полях // Геология и геофизика. 1999. Т. 40, № 9. С. 1275-1290.

[15] Тычков С.А., Червов В.В., Черных Г.Г. Численная модель трехмерной конвекции в верхней мантии Земли // Физика Земли. 2005. № 5. С. 48-64.

[16] Basij A.R., Poreda R.J., Renne P.R. et al. High ЗНе plume origin and temporal-spatial evolution of the Siberian flood basalts // Science. 1995. Vol. 269. P. 822-825.

[17] Тычков С.А., Рычкова E.B., Василевский A.H. Взаимодействие плюма и тепловой конвекции в верхней мантии под континентом // Геология и геофизика. 1998. Т. 39, № 4. С. 413-425.

[18] Тычков С.А., Владимиров А.Г. Модель отрыва субдуцированной океанической литосферы в зоне Индо-Евразийской коллизии // Докл. РАН. 1997. Т. 354, № 2. С. 238-241.

[19] Тычков С.А., Червов В.В., Черных Г.Г. О численном моделировании тепловой конвекции в мантии Земли // Докл. РАН. 2005. Т. 402, № 2. С. 248-254.

[20] Тычков С.А. Конвекция в мантии и динамика платформенных областей. Новосибирск: Наука, 1984. 96 с.

Поступила в редакцию 18 июня 2007 г.