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

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

Том 14, № 5, 2009

Численное моделирование трехмерной конвекции под кратонами Центральной Азии*

В. В. червов Учреждение Российской академии наук Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, Новосибирск, Россия e-mail: [email protected]

Г. Г. Черных, А. В. Червов Учреждение Российской академии наук Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: [email protected]

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

Ключевые слова: тепловая конвекция в мантии Земли, кратоны Центральной Азии, численное моделирование.

Введение

Данные о структуре недр, т. е. о пространственном положении мантийных неоднородно-стей, являются одним из важнейших источников информации о современных процессах в недрах, определяющих тектонический режим территорий. Центральная Азия включает в себя ряд платформенных областей, среди которых можно выделить Таримскую плиту, Северо-Китайский и Южно-Китайский кратоны. В северной части область включает Западно-Сибирскую палеозойскую платформу, или плиту, и древнюю Сибирскую платформу. Эти области существенно влияют на стиль деформирования и тектонический режим литосферы Центральной Азии. Структура континентальной литосферы исследуемой области весьма неоднородна. Значения толщины литосферы древних платформ, таких как Сибирская, Тарим и Китайская платформы, составляют 200... 250 км, в то время как для палеозойской Западно-Сибирской плиты толщина литосферы не превышает 120... 130 км. В рифтовых долинах имеет место утонение литосферы до 40 км. Подобные вариации мощности литосферы существенно влияют на характер мантийных течений [1-3] и поэтому имеют важное значение при проведении численного моделирования трехмерных конвективных течений в мантии Земли. Численное моделирование, которому посвящены многие работы (см., например, [3-9] и приведенную в

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 08-05-00276-а) и СО РАН (интеграционные проекты № 116 и 44).

© ИВТ СО РАН, 2009.

них библиографию), позволяет построить тепловые и скоростные поля под обозначенными структурами. В [7] построена и детально тестирована численная модель тепловой конвекции в мантии Земли, основанная на переменных завихренность—векторный потенциал, методе дробных шагов, последовательности сеток и экстраполяции Ричардсона. Основанная на неявном методе расщепления по физическим процессам трехмерная численная модель конвекции в верхней мантии Земли предложена в [8]. В настоящей работе выполнено моделирование тепловой конвекции под азиатской внутриконтинен-тальной областью, в которую входят Западно-Сибирская плита, Сибирская платформа, Центрально-Азиатский складчатый пояс, Тарим и часть Северо-Китайской платформы. Работа является продолжением [9].

1. Математическая постановка задачи

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

V- V = 0,

V ■ р = Е + Иа ■ Т ■ е,

дТ

V2T.

Здесь р — Рэлея, V

а ■ р ■ дх ■ С3 ■ АТ

давление, Т — температура, Ь — время, И,а

= (и, V, и>) — вектор скорости, е = (0, 0,1), дг — вертикальный размер конвектирующей области, АТ температуропроводности, а — коэффициент теплового расширения, р, г/ плотность и динамическая вязкость, Е— вектор:

число

По ■ X

ускорение силы тяжести, с — Ттах - Тт1п, X — коэффициент характерные

Fx

я,

2 д ди д /ди дv\ д /ди дт дх"Пдх дуП \ду дх) дг^ \дг дх

д

дv

ди

дх'П \ду дх д (ди дт

дх'П \дг+ дх

2 д дv д / дv дт дуП ду дг'П удг ду

д /дv дт\ 2 д дт дуП \дг ду) дгП дг'

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

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

с применением переменных завихренность—векторный потенциал [7]. Число Рэлея, характеризующее режим конвекции, было выбрано как И,а = 2.72 • 105, что отвечает современным представлениям об условиях в недрах Земли. Основные параметры задачи в системе СИ, пригодные для верхней мантии, выбирались следующими: d = 700 000 м, ДТ = 1800 °С, х = 10-6 м2/с, а = 10-5 °С-1, р = 3300 кг/м3, дх = 10 м/с 2, По = 3 • 1021 кг/(м-с).

Вязкость мантийного вещества задавалась в виде

ф,у,г) = еЬх-аТ (х'у'х).

Здесь параметры а = 3.89 и Ь = 5.84 обеспечивают перепад вязкости Птах/Птт от 20 до 200, что присуще верхнемантийным характеристикам течений.

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

В настоящей работе моделирование процессов в верхней мантии было ограничено вну-триконтинентальной областью Азии (рис. 1), в которую вошли Сибирская платформа (лежащая восточнее реки Енисей и простирающаяся в этом направлении до гор Верхо-янья, у подножия которых течет река Лена; на юге эта платформа ограничена озером Байкал, а на севере — Енисей-Хатангской низменностью), Западно-Сибирская плита (примыкающая с запада к Сибирской платформе), Центрально-Азиатский складчатый пояс, Тарим и часть Северо-Китайской платформы — они расположены южнее. Между Таримом и Сибирской платформой находится Тувинский комплекс из мелких крато-нов. Вычисления проводились в параллелепипеде П = [0, 4200] х [0, 4200] х [0, 700] км, на последовательности равномерных сеток 67 х 67 х 36, 101 х 101 х 57 и 133 х 133 х 71 ячеек; величина шагов по времени на соответствующих сетках — 5, 2.5 и 1.25 млн лет.

Значения температуры при этом на двух последних сетках различаются не более чем на 6 % в равномерной сеточной норме.

Схема на рис. 2, где изображена область моделирования конвекции в верхней мантии под Центральной Азией, показывает расположение элементов литосферы в расчетной области. Сибирская платформа, как единый кратон мощностью 220 км, включает в себя два архейских кратона, мощность которых задавалась при моделировании

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

Рис. 2. Схема расположения кратонов и ловушки с указанием мощностей литосферных элементов; между Сибирской платформой и южными кратонами — ловушка, мощность литосферы над которой 60 км; мощность литосферы, не занятой кратонами и ловушкой, — 120 км

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

На рис. 3 и 4 даны градусные сетки — параллели и меридианы. Границы расчетной области построены по геодезическим прямым. Прямое наложение на градусную (в, ф) сетку изображено на рис. 5.

Величина и расположение кратонов приближенно соответствуют реальным данным [10]. Моделирование показало, что, как и в случае прямоугольных в плане кратонов [7, 3, 11], реальные кратоны Цетральной Азии порождают аналогичные структуры. Наблюдаются устойчивые восходящие потоки в виде плюмов; нисходящие потоки и прогретые области по периферии кратонов. Перенос мантийного вещества от оснований кратонов к верхним горизонтам (обтекание) проявляется в виде мелкомасштабной моды конвекции около бортов кратонов.

Реологические особенности Центральной Азии по данным геотермии и сейсмотомо-графии [10] достаточно хорошо прослеживаются на рассчитанных глубинных тепловых полях. Например, совпадают положения центрального плюма и восходящего потока в северо-западной части Сибирской платформы под архейским кратоном, ответственные за трапповый магматизм, который имел место в прошлом Сибирской платформы с середины пермского периода и до начала триасового. Предполагается, что в то время кратон проходил над нижнемантийным плюмом. Сложение температур верхне- и нижнемантийного плюмов привело к излиянию платобазальтов на западе Сибирской платформы. Участие нижнемантийного вещества в излияниях сейчас подтверждено геохимическими и изотопными исследованиями лав Сибири [12].

Рис. 3. Сечение температурного поля в плоскости (ХУ) на глубине 350 км в модели конвекции под литосферой Центральной Азии. Число Рэлея Иа = 271 656

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

с.ш

70°

60° 50° 40°

50° 60° 70° 80° 90° 100° 110° 120° 130° 140° 150° и.д.

Рис. 5. Сечение температурного поля в плоскости (0, р) на глубине 150 км в модели конвекции под литосферой Центральной Азии с наложением на физическую карту

На рис. 4 показаны сечения на глубинах 150 км. Видно, что северо-западный плюм под архейским кратоном мощностью 320 км (см. рис. 3) изменил направление, обтекая подошву кратона, породив восходящий поток около северо-западного угла Сибирской платформы. Под другим архейским кратоном, расположенным на юге Сибирской платформы, также наблюдается восходящий поток, который обтекает кратон и порождает разогретую область у борта кратона, в районе озера Байкал, а также между Тувинским комплексом кратонов и Сибирской платформой. Здесь скорее всего происходит сложение восходящих обтекающих потоков. Ареалы гранитного и бимодального магматизма по периферии кратонов напрямую связаны с восходящими потоками в этих областях, что и подтверждает численный эксперимент.

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

Представленная здесь трехмерная численная модель конвекции также обнаруживает повышение средней мантийной температуры под кратоном на 100 °C, но вместе с тем показывает более сложные формы рельефа кратона, обусловленные динамическим воздействием конвекции.

Для сравнения рельефа Сибирской платформы (см. рис. 1) с результатами вычислений воспользуемся картой распределения температуры в литосфере кратона на глубине 150 км (см. рис. 4). Как было показано в [13], более высокая температура в литосфере соответствует приподнятым участкам литосферы, а пониженная температура — относительно опущенным участкам поверхности. Из сопоставления рельефа платформы с полем температур видно что на платформе, как и в модели, существуют два региональных поднятия — плато Путорана на северо-западе и Патомское нагорье, Алданский щит на юго-востоке Сибирской платформы. Эти два поднятия разделены вытянутой в центральной части низменностью Вилюйской синеклизы (вдоль реки Вилюй) и низменностью вдоль рек Нижняя Тунгуска и Подкаменная Тунгуска. На рис. 4 как следствие

численного эксперимента можно видеть протяженную зону нисходящего потока: темная холодная область (примерно 650 °C) в центральной части Сибирской платформы с направлением от северо-восточной части к юго-западной, между архейским кратонам под плато Путорана и южным архейским кратоном. Таким образом, обнаруживается соответствие существующего рельефа Сибирской платформы результатам трехмерного моделирования конвекции под Сибирским кратоном.

По геолого-геофизическим данным [10], в районе южнее Сибирского кратона и севернее Тарима и Северо-Китайской платформы мощность литосферы составляет от 40 до 75 км. В численной модели толщина литосферного блока в указанном районе принималась равной 60 км. В результате численного моделирования было показано, что в зоне ловушки, как правило, преобладают нисходящие потоки холодного мантийного материала. И в конкретной геологической обстановке, а именно, в случае взаимодействия четырех кратонов, в самой ловушке также наблюдаются цепи классических нисходящих потоков. На глубине 350 км обнаруживается достаточно холодное (650... 750 °C) мантийное вещество (см. рис. 3). Следует заметить, что под территорией ЗападноСибирской плиты, где мощность литосферы составляет 120 км, комплекс нисходящих потоков в среднем на 100 °C выше. В районе озера Байкал, в области ловушки (мощность литосферы 60 км), в непосредственной близости от Сибирского кратона (мощность литосферы 320 км), наблюдается тепловая аномалия в виде мелкомасштабной конвективной ячейки, которая имеет вытянутую форму и может в какой-то степени обьяснить повышенный тепловой поток в Байкальском регионе (рис. 3).

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

Авторы благодарят Н.А. Бушенкову за полезные обсуждения.

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

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

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

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

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

[5] 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. Р. 39-59.

[6] 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.

[7] Tyohkov S.A., Chervov V.V., Chernykh G.G. Numerical modeling of 3D convection in the Earth mantle // Russ. J. Numer. Anal. Math. Modelling. 2005. Vol. 20, N 5. P. 483-500.

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

[9] Тычков С.А., Черных Г.Г., Червов В.В. Трехмерное моделирование конвекции под кратонами Центральной Азии // Вычисл. технологии. 2007. Т. 12. Спец. выпуск 4: Труды V Совещания российско-казахстанской рабочей группы по вычислительным и информационным технологиям, Новосибирск, 6-8 февраля 2007 г. С. 85-95.

[10] Бушенкова Н.А. Неоднородности верхней мантии и современная структура литосферы центральной Сибири по данным сейсмотомографии на отраженных волнах: Автореферат дис. ... к.г.-м.н. Новосибирск, 2004. 20 с.

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

[12] Basü A.R., Poreda R.J., Renne P.R. et al. High 3He plume origin and temporal-spatial evolution of the Siberian flood basalts // Science. 1995. Vol. 269. P. 822-825.

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

Поступила в редакцию 11 января 2009 г., в переработанном виде — 30 марта 2009 г.