CC BY
59
Трехмерная вязкоупругая модель литосферы Центральной Азии: методология построения и численный эксперимент
П.Г. Дядьков, Л.А. Назаров1, Л.А. Назарова1
Институт геофизики СО РАН, Новосибирск, 630090, Россия 1 Институт горного дела СО РАН, Новосибирск, 630091, Россия
В работе рассматривается методология построения трехмерных моделей литосферы для численного расчета полей деформаций и напряжений и их изменений во времени. Разработаны численные схемы, учитывающие сферичность Земли. Создана трехмерная вязкоупругая модель литосферы Центральной Азии, включающая порядка 40 тысяч конечных элементов. Модель учитывает основные тектонические структуры — плиты, микроплиты, области складчатости и сеть главных тектонических разломов. По вертикали модель включает 11 неравных по мощности слоев вплоть до астеносферы. При этом достаточно детально отражена неоднородность реологических свойств верхней и нижней коры, верхней мантии и зон разломов. Модель подготовлена к проведению численных экспериментов. Первые расчеты смещений литосферы Центральной Азии, вызванных движением Индо-Австралийской плиты со скоростью ~50 мм/год в северо-восточном направлении, позволили оценить роль вязкой компоненты в изменении начальных смещений литосферы. Оказалось, что со временем возрастает восточная компонента смещений участков литосферы с пониженными значениями вязкости (Тибет и близлежащие складчатые области). Последнее подтверждается имеющимися экспериментальными данными.
1. Введение
Одним из важных направлений в изучении современных геодинамических процессов является численное моделирование. Численные модели литосферы основных тектонических регионов Земли позволяют не только исследовать современные геодинамические процессы, например, накопление и высвобождение энергии в сейсмически активных зонах, но также лучше понять особенности их развития в прошлом.
Развитие в последнее десятилетие методов космической геодезии позволяет с большой точностью определять скорости перемещений земной поверхности. Число пунктов GPS-наблюдений постоянно растет и получаемая информация совместно с данными мониторинга напряженного состояния (сейсмологические и тектономаг-нитные данные [1]) позволяет осуществлять процесс верификации модели, подбирая реологические свойства среды и граничные условия, которые в наилучшей степени соответствуют реальности.
Центральная Азия и окружающие районы, включая конвергентные границы основных тектонических плит, представляют особый интерес. Во-первых, в этой части Евразии расположена зона Индо-Евразийской коллизии. Результатом столкновения Евразии с Индо-Австралийской плитой и движения последней (продолжающегося и в настоящее время) в северо-восточном направлении явилось создание практически всех основных горных массивов Центральной Азии и, в первую очередь, Гималаев. Во-вторых, определенное влияние на геодинамику Центральной Азии должна оказывать зона Западно-Тихоокеанской субдукции — район Камчатки, Курил, Японии, где осуществляется поддвиг Тихоокеанской плиты под Евразию. Предполагается также наличие ряда внутримантийных источников — плюмов [2], которые могут создавать дополнительные силы у подошвы литосферы.
Рассматривая имеющийся мировой опыт исследования геодинамики отдельных регионов, отметим публи-
© Дядьков П.Г., Назаров Л.А., Назарова Л.А., 2004
кацию [3], в которой для литосферы Ближнего Востока предложена двухмерная модель из жестких блоков с вязкоупругим взаимодействием между ними. Для расчета смещений, вызванных Индо-Евразийской коллизией, в [4] разработана двухмерная модель Азии, однако были введены искусственные ограничения на тип подвижек каждого из разломов. Предлагаемые нами модели свободны от этих недостатков: блоки в них могут быть упругими и вязкоупругими и не ограничивается число степеней свободы движений по разломным зонам.
Ранее нами была разработана и исследована двухмерная квазиупругая модель Центральной Азии, состоящая из четырех тысяч конечных элементов [5]. В модели учтены основные тектонические структуры — плиты, микроплиты и сеть главных тектонических разломов. В зависимости от характеристик разломных зон (меж-плитные, внутриплитные), типа напряженного состояния, в котором находятся эти зоны, степени их современной сейсмической активности, вводилась трехуровневая градация свойств разломов по толщине контакт-элементов, моделирующих разломную зону и их жесткости. Было выполнено моделирование смещений на разломах по нормальному и касательному направлениям, вызванных Индо-Евразийской коллизией со скоростью движения Индостана порядка 50 мм/год.
Для достаточно широкого диапазона граничных условий и свойств среды результаты двухмерного моделирования показали хорошее совпадение с наблюдаемой в реальности картиной распределения областей сжатия и растяжения в Центральной Азии [5]. Поскольку полученные величины относительных смещений берегов разломов в районе Байкальского рифта имели помимо раздвиговой также достаточно большую сдвиговую составляющую, то на основании этого был сделан вывод о том, что только действием сил Индо-Евразийской коллизии трудно объяснить экспериментально наблюдаемое поле деформаций в районе Байкальского рифта. Это, по всей видимости, указывало на наличие источника активного рифтинга в этом районе. Кроме этого, были получены первые оценки изменений деформаций в разломных зонах Центральной Азии, обусловленные влиянием цикличности процесса накопления и реализации упругой энергии при поддвиге субдуктирующей океанической плиты под континентальную в районах Курило-Камчатской и Японской зон субдукции вблизи восточных границ Амурской плиты. В центральной части Байкальского рифта при закрепленных восточных границах Евразийской плиты (что соответствует возникновению «зацепов» между субдуктирующей океанической Тихоокеанской и континентальной Евразийской плитами) возникают условия слабого сжатия (вариант жестких блоков) или чистого сдвига (вариант плит с упругими свойствами). Таким образом, экспериментальные данные об изменениях напряженного состоя-
ния в Байкальском регионе [1, 6], вызываемых сейсмо-геодинамическими процессами в зонах Западно-Тихоокеанской субдукции, подтверждаются результатами численного моделирования. Следствием подобных процессов могут быть изменения механизмов очагов землетрясений и сейсмического режима в рассмотренных внутриконтинентальных районах Евразии в зависимости от стадий цикла сейсмической активности в зоне субдукции.
Рассмотренная двухмерная квазиупругая модель [5] явилась основой для разработки трехмерной вязкоупругой модели Центральной Азии и ее окружения, методология построения которой и результаты первых численных экспериментов рассматриваются в данной работе.
2. Структура и физические свойства трехмерной вязкоупругой модели Центральной Азии
Разбиение по основным структурным тектоническим элементам — плитам, микроплитам — выполнялось аналогично варианту двухмерной модели [5] и незначительно отличается в некоторых деталях от последней (рис. 1). При этом, в качестве первого приближения в основу было положено разбиение, выполненное в работе Зоненшайна и Савостина [6]. При включении в модель основных разломных зон использовались сведения об их расположении и характере смещений бортов разломов, взятые из [7, 8] с привлечением данных о пространственном распределении эпицентров землетрясений.
Для выявления особенностей расположения основных структурных границ земной коры и литосферы по вертикальному разрезу использовались имеющиеся экспериментальные данные и косвенная информация, в основе которой лежали сведения о тепловом потоке (для оценки мощности литосферы) [9, 10]. С этой же целью принимались во внимание данные о пространственном расположении аномалий сейсмических скоростей, полученные на основе сейсмотомографических методов для разных типов волн [11, 12].
Значение коэффициента Пуассона было принято равным 0.25, следуя работе [13].
Для подбора значений вязкости использовались данные, полученные различными методами. Ряд оценок вязкости нижней коры получен в работах, исследующих релаксационный процесс после сильных землетрясений для верхней коры. Эти осредненные значения составляют 1019-5-1019 Па-с [14]. На такие же величины указывают работы, основанные на подборе и исследовании вязкоупругих моделей коры [15-18].
Вязкость верхней мантии оценивается как 5 -1018-
1019 Па-с также из данных о постсейсмической релаксации в районе разлома Сан-Адреас, Калифорния, а для Канады по данным VLBI-измерений ~1021 Па-с [14, 19].
Рис. 1. Схема разбиения модели по основным структурным тектоническим элементам. Номера элементов на схеме в кружках соответствуют номерам этих тектонических элементов в таблице 1. Для зон разломов введена следующая градация: межплитные границы сдвигового типа — тип 1; межплитные границы в зонах сжатия — тип 2; типы 3-5 — это, в основном, разломные зоны внутриплитные или между микроплитами, причем тип 3 — это разломы в областях близгоризонтального растяжения, тип 4 — в областях сдвига и тип 5 — в областях сжатия или относительно слабоактивные
Для астеносферы под континентами оценки вязкости, основанные на данных о скоростях послеледниковых поднятий, составляют 1020—1021 Па-с [20]. Для океанической астеносферы по скорости распространения изменений напряжений, вызванных сильными землетрясениями, получено значение вязкости 5-1017 Па-с [21].
Используя данные расчетов осредненных по вертикали девиаторных напряжений для литосферы Азии и простую связь между скоростью деформаций и эффективной вязкостью, в работе [22] были рассчитаны средние значения вязкости литосферы для районов Азии.
Как все вышеперечисленные данные, так и сведения о распределении прочности по разрезу литосферы, полученные по данным о тепловом потоке для регионов Азии [10] были использованы нами для задания в трехмерной модели литосферы Азии значений вязкости (таблица 1).
Весьма важными параметрами для модели являются реологические свойства разломных зон. Для оценки их вязкости в верхней и средней коре были использованы экспериментальные данные о характерных временах накопления упругой энергии при действии разного типа нагрузок [23, 24]. По своим реологическим свойствам разломы были поделены на 5 типов (рис. 1). Первые два типа относятся к границам основных тектонических плит. В модели принято, что такие разломные зоны характеризуются ослабленной реологией, пронизывающей всю литосферу вплоть до астеносферы. При этом разломные зоны (межплитные границы) сдвигового типа (тип 1 на рис. 1) имеют более низкую вязкость (1019 Па-с) по сравнению с теми, которые находятся в зонах коллизии или субдукции, где преобладает близгоризонтальное сжатие (тип 2). Для последних принята вязкость 1020 Па-с в коре и 1021 Па-с от границы Мохо до астеносферы. Ширина зоны для первых двух
Таблица 1
Свойства среды для объемной модели Центральной Азии и ее окружения
№ слоя 1 2 3 4 5 6 7 8 9 10 11
Н 10 10 10 5 5 5 35 20 50 100 50
ЕН 10 20 30 35 40 45 80 100 150 250 300
Е 40 50 60 70 70 70 70 70 70 70 70
Евразия 1 Р 2 800 2 800 2900 2 900 2900 3 300 3 300 3 300 3 300 3 300 3 300
^п 24 25 24 22 21 23 23 23 23 22 20
Северная Е 40 50 70 70 70 70 70 70 70 70 70
Америка Р 2 800 2 800 2900 2 900 3 300 3 300 3 300 3 300 3 300 3 300 3 300
2 ^п 24 25 23 21 23 23 23 23 22 20 20
Охотоморская Е 40 50 70 70 70 70 70 70 70 70 70
плита Р 2 800 2 800 2900 3 300 3 300 3 300 3 300 3 300 3 300 3 300 3 300
3 ^п 24 24 21 23 23 23 22 19 19 19 19
Тихоокеанская Е 35 70 70 70 70 70 70 70 70 70 70
плита Р 2 900 3 400 3 400 3 400 3 400 3 400 3 400 3 400 3 400 3 400 3 400
4 ^п 25 24 24 24 24 23 22 18 18 18 18
Филиппинская Е 35 70 70 70 70 70 70 70 70 70 70
плита Р 2 900 3 400 3 400 3 400 3 400 3 400 3 400 3 400 3 400 3 400 3 400
5 ^п 25 24 24 24 24 23 22 18 18 18 18
Китайская Е 40 50 50 50 50 70 70 70 70 70 70
плита Р 2 800 2 800 2900 2 900 3 300 3 300 3 300 3 300 3 300 3 300 3 300
6 ^п 24 24 22 20 23 23 23 22 19 19 19
Индокитайская Е 30 50 50 50 50 70 70 70 70 70 70
плита Р 2 800 2 800 2900 2 900 2900 3 300 3 300 3 300 3 300 3 300 3 300
7 ^п 23 23 22 21 19 23 23 22 19 19 19
Индо-Австра- Е 40 50 60 70 70 100 100 100 80 70 70
лийская плита Р 2 800 2 800 2900 2 900 2900 3 300 3 300 3 300 3 300 3 300 3 300
8 ^п 24 25 24 23 21 24 24 23 23 19 19
Аравийская Е 40 50 50 50 70 70 70 70 70 70 70
плита Р 2 800 2 800 2 800 2 900 3 300 3 300 3 300 3 300 3 300 3 300 3 300
9 ^п 22 23 22 21 20 23 22 21 19 19 19
Афганская Е 25 30 40 50 50 50 70 70 70 70 70
плита Р 2 800 2 800 2900 2 900 2900 2900 3 300 3 300 3 300 3 300 3 300
10 ^п 23 23 22 21 20 19 23 19 19 19 19
Е 40 60 70 70 70 70 70 70 70 70
Тарим 11 Р 2 800 2 800 2 800 2 900 2900 2900 3 300 3 300 3 300 3 300 3 300
lgп 23 24 24 22 21 20 24 23 19 19 19
Е 25 40 50 50 50 50 70 70 10 70 70
Тянь-Шань 12 Р 2 800 2 800 2 800 2 900 2900 2900 3 300 3 300 3 300 3 300 3 300
lgп 22 23 22 21 20 19 22 21 20 19 19
Е 20 20 40 50 50 50 60 70 70 70 70
Тибет 13 Р 2 800 2 800 2 800 2 800 2 800 2 800 3 300 3 300 3 300 3 300 3 300
lgп 22 23 22 21 20 19 21 20 19 19 19
Продолжение таблицы 1
Северный Е 20 30 50 50 50 50 60 70 70 70 70
Тибет Р 2 800 2 800 2 800 2 800 2 800 2900 3 300 3 300 3 300 3 300 3 300
14 23 24 22 21 19 18 22 21 19 19 19
Е 30 50 50 50 50 50 70 70 70 70 70
Гоби 15 Р 2 800 2 800 2 800 2 800 2900 2900 3 300 3 300 3 300 3 300 3 300
23 23 21 20 19 19 23 22 19 19 19
Е 30 50 50 50 50 50 70 70 70 70 70
Джунгария 16 Р 2 800 2 800 2 800 2 900 2900 2900 3 300 3 300 3 300 3 300 3 300
23 24 22 21 20 19 23 22 21 19 19
Е 30 50 50 70 70 100 100 100 100 70 70
Ордос 17 Р 2 800 2 800 2 800 2 900 2900 3 300 3 300 3 300 3 300 3 300
23 24 22 21 21 24 23 23 22 19 19
Монгольская Е 20 40 50 50 50 70 70 70 70 70 70
плита Р 2 800 2 800 2 800 2 900 2900 3 300 3 300 3 300 3 300 3 300 3 300
18 22 24 22 20 19 24 23 23 22 19 19
Е 30 50 50 50 50 70 70 70 70 70 70
Байкал 19 Р 2 800 2 800 2900 2 900 2900 3 300 3 300 3 300 3 300 3 300 3 300
23 24 22 19 18 22 22 19 19 19 19
Амурская плита Е 40 50 70 70 70 100 100 100 70 70 70
Запад Р 2 800 2 800 2 800 2 900 2900 3 300 3 300 3 300 3 300 3 300 3 300
20 23 25 23 22 21 24 23 22 19 19 19
Амурская плита Е 30 50 60 70 100 100 100 100 70 70 70
Восток Р 2 800 2 800 2 800 2 900 3 300 3 300 3 300 3 300 3 300 3 300 3 300
21 23 24 22 20 24 24 23 22 19 19 19
Обозначения: Н — мощность слоя в км; ХН — глубина до подошвы слоя в км; Е — модуль упругости, ГПа; р — плотность, кг/м3; ^ п — логарифм вязкости [п, Па-с]
типов составляет в модели 30 км. Три следующих типа разломов относятся, в основном, к внутриплитным, при этом часть из них — это границы между микроплитами, блоками. Их ширина принята равной 10 км. Наиболее ослабленные из них те, которые находятся в условиях растяжения (тип 3 на рис. 1). Их вязкость составляет
1019 Па-с вплоть до астеносферы. Разломы, преимущественное движение по которым сдвиговое, относятся к типу 4, при этом их вязкость составляет 1019 Па-с в коре и далее до астеносферы 1020 Па-с. Тип 5 — это внутриплитные разломы, находящиеся в условиях сжатия или относительно слабоактивные. Их вязкость равна
1020 Па-с в коре и 1021 Па-с от Мохо до астеносферы.
3. Постановка задачи построения численной трехмерной модели литосферы в сферической системе координат
Если линейные размеры исследуемой области сравнимы с радиусом Земли (Л0 = 6 370 км), то использо-
вание плоских моделей для расчета полей напряжений и деформаций, а также для анализа их эволюции во времени возможно лишь с существенными оговорками. К недостаткам разработанной нами ранее двумерной модели [5] можно отнести:
- невозможность получить информацию о распределении напряжений и деформаций и их изменений во времени по вертикальному разрезу литосферы;
- неучет сил гравитации;
- невозможность описать влияние рельефа на напряженно-деформированное состояние;
- искажение геометрии структурных элементов объекта вследствие представления сферической поверхности Земли на плоскости (картографическая проекция).
Эти обстоятельства показали необходимость создания объемной модели.
В сферической системе координат (Л, 0, ф) (Л — радиус, 0 — широта, ф — долгота) с началом в центре
Рис. 2. Схема расчетной области в сферической системе координат
Земли рассмотрим область G = {Л1 < Л < Л0, 01 < 0 < <02, ф1 <ф<ф2} (рис. 2). Значения 01 = 10° СШ и 02 = 100° (10° ЮШ), ф1 = 30° ВД и ф1 = 180° ВД выбраны так, чтобы в G находился исследуемый объект — Центральная Азия.
Примем, как и для плоской задачи, что среда состоит из упругих блоков, а необратимые деформации могут иметь место только на межблочных границах (контактах), мощность которых много меньше характерных размеров слагающих их блоков. Приведем полную систему уравнений, описывающих деформирование такой среды.
Уравнения равновесия
Эа
RR
2аRR - а99 - афф + аЯфС^Ф
эя
1 Эа
Я9
я 1 Эа
Яф
Я sin ф Э9 Я Эф
+ № = 0,
да Яф 3а Яф + (а99 афф )С^ ф
дя +
Я
1 да9ф
------------+------------------
Я дф Я sin ф д9
1 да
фф
= 0,
(1)
да Я9 + 3а Я9 + 2а9фС^ф + 1 да9ф +
дЯ Я
1 да
Я дф
+
99
Я sin ф д9
= 0.
Уравнения состояния для упругого материала блоков (закон Гука)
аЯЯ = Х8 + 2^8ЯЯ, а99 = Х8 + 2М^899,
афф = Х8 + 2М8фф, а Я9 = 2Ц8 Я9,
а Яф = 2Ц8 Яф, а9ф = 2М89ф.
Связь деформаций со смещениями
(2)
8
диЯ
ЯЯ~ 1т 9
99
фф
иЯ
ди9 иф
+ -^ ^ф,
Я Я sin ф д9 Я
иЯ
Я
8 Я 9 =
8 Я ф =
Эиф
Я дф ’
1 диЯ ди9 и9
(3)
Я sin ф д9 дЯ Я
1 диЯ ди
- + -
ф иф
Я дф дЯ Я
у
1 диф 1 ди9 и9
------------—^ +--------—9--------9 ^ф
Я sin ф д9 Я дф Я
В (1)-(3) агу, 8л — компоненты тензоров напряжений и деформаций (/, ] = Л, 0, ф), и = (ил, иф, и0) — вектор смещений, 8 = 8 ЛЛ + 8фф + 800 — объемная деформация, Л, ц —параметры Ламе, которые выражаются через модуль Юнга Е и коэффициент Пуассона V по формулам
vE
(1 + у )(1 - 2у)
2(1 + v)
(4)
|Ч > Ґ \ w1 ' Г К11 К12 К13 Ї
Т 2 = К w2 , К = К 21 2 2 К 23
,Т3) w3 3 1К31 К32 К33 1
На стадии допредельного деформирования определяющие уравнения для межблочного контакта 3 имеют вид [25]:
(5)
где т к — компоненты вектора напряжений на У; w1 = = (и - и ) • п — сближение берегов контакта (п — вектор нормали к У); w2 и w3 — проскальзывания в касательной к У плоскости в двух взаимно перпендикулярных направлениях; Кт — жесткости (к, т = 1, 2, 3). Эксперименты [26] показали, что внедиагональные элементы К, отвечающие за дилатансию, по крайней мере, на порядок меньше Ктт, поэтому из (5) следует:
т1 = К1^1, т 2 = К 22 т3 = К33 w3•
+
+
+
Для численных оценок использовались предложенные в [27] соотношения
Кп = , К22 = К33 = Д/^
где Е и ц — средние значения соответствующих параметров блоков, слагающих нарушение.
4. Вязкоупругая модель
Для описания эволюции напряженно-деформированного состояния во времени используются реологические модели [28], которые обычно «конструируются» из упругих, вязких и пластических элементов, а также элементов с сухим трением. Среди наиболее известных отметим среду Максвелла, Кельвина-Фойгта, Пойнтин-га-Томпсона, Шведова, Бингама, Сен-Венана [28, 29]. Числовые значения параметров таких моделей определяются из прямых экспериментов на ползучесть (мягкое нагружение) или релаксацию (жесткое нагружение).
Для исследования вязкого деформирования крупномасштабных объектов в [30] предложено использовать модель Максвелла. Здесь выбрана та же модель, поскольку именно по ней выполнены числовые оценки вязкости (таблица) структурных элементов — плит и микроплит.
Основываясь на одномерной модели, выведем определяющие уравнения для вязкоупругой среды Максвелла в общем случае, воспользовавшись методом замороженных коэффициентов [31], а также экспериментально установленным фактом о том, что объем деформируется упруго — первые инварианты тензоров напряжений и деформаций связаны линейной зависимостью, а вязкие эффекты описываются специфической зависимостью между девиаторами этих тензоров.
В одномерном случае деформирование среды Максвелла описывается уравнением
а- = 2Ц8й ,
(11)
£ £ е = —\—.
Ц П
(7)
где £, е — напряжение и деформация; ц — модуль упругости; п — вязкость. Если в момент времени г = 0 к системе внезапно приложена деформация е = е0, решение (7) имеет вид
) = Ц1(г)е0>
) = це~Ц‘1п.
(8)
(9)
Откуда следует, что £(г) ^ 0 при г ^ 0 — происходит релаксация напряжений.
Из закона Гука (2) следует, что
(10)
где а = (аЛЛ + а00 + афф)/3 — среднее напряжение. Вычитая (10) из первых трех соотношений (2), получим
где а- = а- - а и 8- = 8- -8Й — компоненты девиато-ров напряжений и деформаций. Теперь, формально заменяя в (11) ц на ц1(г) и учитывая (10), найдем аналог закона Гука для вязкоупругой среды:
а- (г) = ^(г )е8- + 2ц1(г)8г-/,
Л1(г)=Л+3(ц1(г) -ц)>
(12)
(13)
8- — символ Кронекера. Вид (11) свидетельствует о том, что для решения вязкоупругих задач можно применять методы классической теории упругости, заменяя в них параметры Ламе на обобщенные параметры — функции времени (9) и (13).
5. Граничные условия и метод решения
Граничные условия. Для всех моделируемых процессов на сферических границах задавались одни и те же условия. Верхняя граница — свободная:
алл (Л, 0, ф) = ал0 (Л, 0, ф) = алф (Л, 0, ф) = 0,
на нижней — отсутствуют касательные напряжения и вертикальное смещение (подробное обоснование этого предположения выполнено в [32]):
ал0 (Лх, 0, ф) = алф (Лъ 0, ф) = 0,
иЛ (Л1’ 0> ф) = 0
(14)
На рис. 3 жирными линиями показана разломно-бло-ковая структура на поверхности Земли Л = Л0. Обозначим А1, к, А10 характерные точки на этой поверхности, а под В1, к, В10 будем подразумевать соответствующие им точки на нижней границе области G (на рис. 3 последние не видны). Условия на границах формулировались в зависимости от типа моделируемого процесса (коллизия, субдукция и т.д.). В частности, при исследовании кинематики плит, обусловленной движением Индо-Австралийской плиты:
на Л2В2 В3 А3 и А^^а А10
аЛ0 (Л, 02, ф) = а00 (Л, 02, ф) = = а0ф (Л, 02, ф) = 0,
(15)
на Л3В3В4 Л4
аЛф (Л 0 ф1) =афф (Л 0 ф1) =
= а0ф(Л 0 ф1) = °> (16)
граница A4A5...Aх0 В10В9...В4 неподвижна
иЛ = и0 = иф= °> (17)
Рис. 3. Схема области моделирования, поясняющая задаваемые граничные условия
перемещение границы A1B1B2A2 плиты происходит только вдоль сферической поверхности (рис. 3)
UR (R, 02’ ф) = 0
и0(R, 02, ф) = -U0 cos а, (18)
Мф(R, 02’ Ф) = UoSinа,
где U0 и а — амплитуда подвижки и угол, характеризующий ее направление.
Метод решения. Решение краевой задачи (1)-(3), (14)—(18) осуществлялось методом конечных элементов. Не вдаваясь в детали построения алгоритма реализации метода конечных элементов, которые можно найти в [33], отметим лишь некоторые особенности.
Дискретизация расчетной области (тонкие линии на рис. 3) проводилась в три этапа:
- в сечении R = R0 строилась плоская сетка, ассоциированная с тектонической структурой;
- на разломах число узлов удваивалось (для описания взаимодействия берегов разломов);
- из плоских собиралась трехмерная сетка гекса-эдрических элементов.
Для каждого элемента матрица жесткости находилась численным интегрированием (использовалась десятиточечная схема Гаусса). Нумерация узлов выполнялась так, чтобы минимизировать ширину ленты В глобальной матрицы жесткости. Для регулярной сетки ^, М, N — число узлов по координатам) количество неизвестных составляет 3LMN, а В = 3^М + М + 2). Тогда при выполнении неравенств N > L > М размерность разрешающей системы линейных уравнений (и, следовательно, время расчета) будет минимальной.
Информация о глубинном строении исследуемого объекта наименее детальна, поэтому М — число узлов по Л, L — по 0 и N — по ф. В расчетах принято L = 57, М = 12, N = 67, таким образом, число конечных элементов — 40040, число неизвестных — 135432.
6. Некоторые результаты расчетов
С использованием разработанной модели выполнен ряд первых расчетов смещений в районах Центральной
Рис. 4. Результат модельного расчета векторов смещений поверхности литосферы, вызванного продвижением Индо-Австралийской плиты в северо-восточном направлении на 50 мм (годовое смещение). Черные стрелки — смещение через год, серые — через 200 лет
Азии, вызванных движением Индо-Австралийской плиты. Поля горизонтальных смещений были рассчитаны для разных глубин. На рис. 4 показаны результаты расчета поля смещений, возникших за счет продвижения на северо-восток Индо-Австралийской плиты на ~50 мм. Черные стрелки — величина первоначального упругого смещения литосферы на поверхности Земли. Серыми стрелками показаны смещения через 200 лет. Этот расчет хорошо иллюстрирует роль вязкости в изменении со временем величины и направления смещения. На больших временах вклад вязкой деформации увеличивается. Если не принимать во внимание области, близкие к границам модели (краевой эффект), то из рис. 4 хорошо видно, что наибольшее отклонение на восток серых стрелок (результат смещения через 200 лет) от черных наблюдается в областях Тибета, Гималаев и некоторых соседних с ними складчатых областях. Наиболее вероятная физическая интерпретация этого эффекта — выдавливание в восточном направлении материала с пониженной вязкостью, который зажат между сближающимися жесткими плитами — Индией и жесткой северной частью Евразийской плиты. Этот результат подтверждается имеющимися экспериментальными данными [34]. Полученные расчетные данные свидетельствуют о весьма удовлетворительных первых результатах апробации модели.
7. Заключение
Для численного моделирования полей деформаций и напряжений и их изменений во времени разработана трехмерная вязкоупругая модель литосферы Центральной Азии, включающая порядка 40 тысяч конечных элементов. Разработаны численные схемы, учитывающие сферичность Земли. Эта трехмерная вязкоупругая модель подготовлена к проведению численных экспериментов с использованием метода конечных элементов. По вертикали модель включает 11 неравных по мощности слоев вплоть до астеносферы. При этом достаточно детально отражена неоднородность реологических свойств верхней и нижней коры, верхней мантии и зон разломов. Выполнены первые расчеты смещений на разных глубинах литосферы, вызванные движением Индо-Австралийской плиты со скоростью ~50 мм/год в северо-восточном направлении. Оценена роль вязкой компоненты в изменении со временем начальных смещений литосферы: со временем возрастает восточная компонента смещений участков литосферы с пониженными значениями вязкости (Тибет и близлежащие складчатые области). Последнее связывается с «выдавливанием» материала с пониженной вязкостью, находящегося между сближающимися жесткими плитами — Индией и жесткой северной частью Евразийской плиты.
Работа поддержана проектом РФФИ 03-05-65418, Программой Отделения наук о Земле РАН 6.5 (проект
6.5.1), Интеграционным проектом СО РАН № 101, Программой Президиума РАН № 13.
Литература
1. ДядъковП.Г., Мелъникова В.И., НазаровЛ.А., Назарова Л.А., Санъ-
ков В.А. Сейсмотектоническая активизация Байкальского региона в 1989-1995 годах: результаты экспериментальных наблюдений и численное моделирование изменений напряженно-деформированного состояния // Геология и геофизика. - 1999. - Т. 40. - №2 3. -С. 373-386.
2. Добрецов Н.Л., Кирдяшкин А.Г. Глубинная геодинамика // Труды
ОИГГМ СО РАН. - Новосибирск: Изд. СО РАН, 2002. - Вып. 830. - 300 с.
3. Соболев П.О., Соловъев А.А., Ротвайн И.М. Моделирование динамики литосферы и сейсмичности для региона Ближнего Востока // Вычислительная сейсмология. - 1996. - Вып. 28. - С. 131137.
4. Peltzer G., Saucier F. Present-day kinematics of Asia derived from geologic fault rates // JGR. - V. 101. - No. B12. - P. 27943-27956.
5. Назарова Л.А., Назаров Л.А., Дядъков П.Г. Математическое моделирование кинематики плит Центральной Азии // ФТПРПИ. -2002. - № 5. - С. 3-9.
6. Зоненшайн Л.П., Савостин Л.А. Введение в геодинамику. - М.: Недра, 1979. - 276 с.
7. ТрифоновВ.Г. Неотектоника Евразии. - М.: Научный мир, 1999. -252 с.
8. Trifonov V.G. The map of active faults in Eurasia: principles, methods, and results // Journal of Earthquake Prediction Research. - 1996. -V. 5. - P. 326-347.
9. Artemieva I.M. and Mooney W.D. Thermal thickness and evolution of
Precambrian lithosphere: A global study // JGR. - 2001. - V. 106. -No. B8. - P. 16387-16414.
10. Wang Y. Heat flow pattern and lateral variations of lithosphere strength in China mainland: constraints on active deformation // Physics of the Earth and Planetary Interiors. - 2001. - V. 126. - P. 121-146.
11. Yanovskaya T.B., Antonova L.M., Kozhevnikov VM. Lateral variations of the upper mantle structure in Eurasia from group velocities of surface waves // Phys. Earth Planet. Inter. - 2000. - V. 122. - P. 19-32.
12. Kulakov I., Tychkov S., Bushenkova N. et al. Structure and dynamics of the upper mantle beneath the Alpine-Himalayan orogenic belt, from teleseismic tomography // Tectonophysics. - 2002. - V. 358. - P. 7796.
13. Grollimund B., Zoback M.D. Post glacial lithospheric flexure and induced stresses and pore pressure changes in the northern North Sea // Tectonophysics. - 2000. - V. 327. - P. 61-81.
14. PollitzF.F., Burgmann R. Joint estimation of afterslip rate and postseis-mic relaxation following the Loma Prieta earthquake // JGR. - 1998. -V. 103. - No. B11. - P. 26975-26992.
15. Roy M., Royden L.H. Crustal rheology and faulting at strike-slip plate boundaries. 1. An analitic model // JGR. - 2000. - V. 105. - No. B3. -P. 5583-5597.
16. Roy M., Royden L.H. Crustal rheology and faulting at strike-slip plate boundaries. 2. Effects of lower crustal flow // JGR. - 2000. - V. 105. -No. B3. - P. 5599-5613.
17. McKenzie D., Nimmo F., Jackson J.A. Characteristics and consequences of flow in the lower crust // JGR. - 2000. - V. 105. - No. B5. -P. 11029-11046.
18. Дядъков П.Г., Назаров Л.А., Назарова Л.А. Вязкоупругая модель Байкальского рифта для численного моделирования сейсмогеоди-намических процессов // Физические основы прогнозирования разрушения горных пород: Тезисы докл. 1-й Междунар. школы-семинара / Отв. ред. проф. В.А. Мансуров. - Красноярск: САА, 2001. - С. 128.
19. Pollitz F.F., Peltzer G., Burgmann R. Mobility of continental mantle: Evidence from postseismic geodetic observations following the
Landers earthquake // JGR. - 2000. - V. 105. - No. B4. - P. 80358054.
20. Ботт М. Внутреннее строение Земли. - М.: Мир, 1974. - 374 с.
21. PollitzF.F., Burgmann R., Romanowich B. Viscosity of oceanic asthe-nosphere inferred from remote triggering of earthquakes // Science. -1998. - V. 280. - P. 1245-1249.
22. Flesch L.M., Haines A.J., Holt W.E. Dynamics of the India-Eurasia collision zone // JGR. - 2001. - V. 106. - No. B8. - P. 16435-16460.
23. Доброволъский И.П. Теория подготовки тектонического землетрясения. - М.: ИФЗ АН СССР, 1991. - 220 с.
24. Дядъков П.Г. Вызванная сейсмичность в центральной части оз. Байкал и использование ее характеристик для оценки эффективной вязкости разломной зоны // Тезисы докл. 1-й Междунар. школы-семинара / Отв. ред. проф. В.А. Мансуров. - Красноярск: САА, 2001. - 60 с.
25. Назарова Л.А. Моделирование объемных полей напряжений в разломных зонах земной коры // Доклады РАН. - 1995. - Т. 342. -№ 6. - С. 804-808.
26. Jing L., Nordiund E., Stephansson О. An experimental study on the anisotropy and stress-dependency of the strength and deformability of rock joints // Int. J. of Rock Mech. and Min. Sci. & Geomech. Abstr. - 1992. - V. 29. - No. 6. - P. 535-542.
27. Юфин C.A. Механические процессы в природном массиве и их взаимодействие с подземными сооружениями / Автореф. дис. ... докт. техн. наук. - Москва: МГИ, 1991. - 36 с.
28. Bland D.R. The theory of linear viscoelasticity. - Oxford-London-New York: Pergamon Press, 1960. - 125 p.
29. Руппенейт К.В., Либерман Ю.В. Введение в механику горных пород. - М.: Госгортехиздат, 1960. - 356 с.
30. Carey S.W. The rheid concept in geotectonics // J. Geol. Soc. Aust. -1953. - V. 1. - P. 67-117.
31. Работнов Ю.Н. Механика деформируемого твердого тела. - М.: Наука, 1979. - 744 с.
32. Дядъков П.Г., Назаров Л.А., Назарова Л.А. Численное моделро-вание наряженного состояния земной rapbi и условий возникновения динамической неустойчивости сейсмоактивных pазломов ри pифтогенезе // Геология и геофизика. - 1997. - Т. 38. - №2 12. -С. 2001-2010.
33. Зенкевич O. Метод конечных элементов в технике. - М.: Мир, 1975. - 541 с.
34. Molnar P., Gibson J.M. A bound on the rheology of continental lithosphere using very long baseline interpherometry: the velocity of south China with respect to Eurasia // JGR. - 1996. - V. 101. - No. B1. -P. 545-553.
Three-dimensional visco-elastic model of the Central Asia lithosphere: methodology of development and numerical experiment
P.G. Djadkov, L.A. Nazarov1, and L.A. Nazarova1
Institute of Geophysics SB RAS, Novosibirsk, 630090, Russia 1 Institute of Mining SB RAS, Novosibirsk, 630091, Russia
In the paper consideration is given to the methodology of the development of three-dimensional lithosphere models for calculating numerically the strain and stress fields and their changing with time. Numerical schemes allowing for the Earth’s sphericity are developed. A three-dimensional visco-elastic model of the Central Asia lithosphere including about 40 000 finite elements is constructed. The model allows for the main tectonic structures, namely, plates, microplates, folding regions and a net of major tectonic faults. Vertically the model consists of 11 layers varying in thickness up to the asthenosphere. An inhomogeneity of rheological properties of the upper and lower crust, upper mantle and fault zones are represented in detail. The model is prepared to carry out numerical experiments. The first calculations of the lithosphere displacements in Central Asia induced by the Indo-Australian plate movements at a velocity of ~ 50 mm/year in the north-east direction have allowed us to evaluate the role of the viscous component in changing of the initial lithosphere displacements. It is found that the east component of displacements of the lithosphere parts with low viscosity (Tibet and neighboring folding zones) increases with time. The latter is confirmed by the available experimental data.
