CC BY
72
Трехмерная модель литосферы Земли: методология построения и первые численные эксперименты
П.Г. Дядьков, Л.А. Назаров1, Л.А. Назарова1
Институт геофизики СО РАН, Новосибирск, 630090, Россия 1 Институт горного дела СО РАН, Новосибирск, 630091, Россия
Сформулированы основные положения при построении трехмерной вязкоупругой модели литосферы Земли для численных расчетов смещений, дефомаций и напряжений, вызванных как глобальными, так и региональными сейсмогеодинамическими процессами. С учетом разработанной методологии выполнено построение этой модели — разбиение на сеть конечных элементов, которая учитывает геометрическую структуру основных тектонических плит и разломных зон между ними, а также в первом приближении отражает реологию основных тектонических элементов. Первые тестовые расчеты показали работоспособность этой модели и ее перспективность для использования в моделировании современных сейсмогеодинамических процессов.
3D-model of the Earth’s lithosphere: development methodology and first numerical experiments
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
The paper formulates basic development principles of a 3D viscoelastic model of the Earth’s lithosphere for numerical calculations of displacements, strains, and stresses induced by both global and regional seismic-geodynamic processes. With consideration for the developed methodology we have constructed the model trough the formation a finite element mesh that accounts for the geometric structure of main tectonic plates and fault zones between them and reflecting rheology of main tectonic elements to a first approximation. First test calculations show that the model is efficient and promising for simulation of modern seismic-geodynamic processes.
1. Введение
Разработка геомеханических моделей литосферы Земли имеет важное фундаментальное и прикладное значение. Лучшей иллюстрацией необходимости разработки таких моделей являются происходящие сейсмо-геодинамические процессы, влияние которых явно превосходит региональный уровень. Это, во-первых, наличие в настоящее время глобальной сейсмической активизации, начавшейся, вероятнее всего, в конце 80-х годов прошлого столетия [1]. Наличие этой активизации подтвердилось рядом сильнейших землетрясений последнего времени, в том числе катастрофическим землетрясением у берегов острова Суматры в конце 2004 года. Моделирование подобных процессов позволит понять их природу и в дальнейшем поможет прогнозировать
их влияние на окружающие регионы и особенности глобального состояния литосферы Земли. Основой для создания обсуждаемой здесь глобальной модели литосферы явился ряд разработанных ранее нами моделей для отдельных регионов [2-4].
Экспериментальной основой проверки правильности задания структурных и реологических свойств модели являются, в первую очередь, данные спутниковой геодезии. В последнее десятилетие банк данных GPS-наблюдений постоянно пополняется и позволяет с большой точностью определять скорости перемещений земной поверхности. Получаемая информация о смещениях и деформациях земной поверхности совместно с сейсмологическими данными о действующих в литосфере напряжениях (механизмы очагов землетрясений)
© Дядьков П.Г., Назаров Л.А., Назарова Л.А., 2006
позволяет осуществлять процесс верификации модели, подбирая реологические свойства среды и граничные условия, которые в наилучшей степени соответствуют реальности.
2. Обоснование модели: структура и физические свойства
Отметим основные факторы, которые следует учитывать при построении трехмерной модели литосферы Земли.
2.1. Структура
Во-первых, алгоритм расчета должен быть адаптирован для решения геомеханических задач с учетом сферичности Земли. Во-вторых, структура модели должна учитывать основные структурные элементы литосферы — тектонические плиты и разломные зоны на их границах. Следует включить в модель и некоторые основные внутриплитные неоднородности — жесткие блоки и микроплиты, а также зоны повышенной деформируемости. В-третьих, необходимо, чтобы модель учитывала вертикальную расслоенность литосферы, хотя бы тех слоев, которые наиболее сильно отличаются друг от друга по своим реологическим свойствам — это верхняя и нижняя кора и верхняя мантия. Необходимо также учесть в основании модели наличие астеносферы.
Зоны субдукции, океанического спрединга и коллизии являются важнейшими тектоническими элементами, структурные особенности которых должны найти отображение в модели. Если границы плит в зонах спре-динга и коллизии можно в грубом приближении моделировать близвертикальными контактами, то зона суб-дукции требует учета хотя бы среднего наклона зоны Беньофа в литосфере висячего литосферного крыла.
При проведении численного моделирования может возникнуть необходимость изучения влияния процессов, имеющих глобальный или мегарегиональный уровень, на геодинамическую ситуацию в отдельном ре-
гионе. В этом случае изучаемый локальный район может быть представлен на модели более детальной сеткой разбиения на конечные элементы.
2.2. Реология
На данном этапе, также как и в [4], выбрана вязкоупругая модель литосферы (модель Максвелла), которая позволяет выполнять моделирование как быстрых, фактически упругих, процессов (акты землетрясений), так и более медленных, например релаксационных процессов, при которых требуется учитывать неупругие деформации среды. Ранее опробование подобной модели для региона Центральной Азии показало вполне удовлетворительные результаты [4].
Для подбора значений вязкости коры и мантии использовались данные из литературных источников, полученные различными методами: на основе изучения релаксационного процесса после сильных землетрясений, подбора и исследования вязкоупругих моделей коры, данных о скоростях послеледниковых поднятий [5-11].
Весьма важными параметрами являются реологические свойства разломных зон (границ плит и микроплит). Для оценки их вязкости в верхней и средней коре были использованы экспериментальные данные о характерных временах накопления упругой энергии при действии разного типа нагрузок.
2.3. Граничные условия
Для моделируемых процессов на сферических границах задаются следующие условия. Верхняя граница — свободная, а на нижней отсутствуют касательные напряжения и вертикальное смещение (подробное обоснование этого предположения выполнено в [2]). Условия на границах плит должны задаваться в зависимости от типа моделируемого процесса в данной части модели (спрединг, коллизия, субдукция и т.д.). Другой (или параллельный) вариант задания действующих сил — приложение их к основанию литосферы. Таким способом
Рис. 1. Тектонические плиты Земли
Рис. 2. Схема разбиения земной поверхности на сеть конечных элементов с учетом границ тектонических плит
возможен учет сил волочения плит за счет вязкого трения у основания литосферы при конвективных течениях вещества мантии.
3. Модель и результаты первых численных экспериментов
Основные тектонические плиты и схема разбиения земной поверхности на сеть конечных элементов показаны на рис. 1, 2.
Распределение упругих и вязких свойств задавалось по глубине аналогично модели [4], но разбиение по вертикали сделано на 6 основных слоев (вместо 11 в модели [4]). Свойства разломных зон также задавались подобные тем, которые использовались для межплит-ных и межблочных контактов в модели литосферы
90
50 О 50 100 150
Рис. 3. Смещения земной поверхности, вызванные частичным сцеплением (увеличением жесткости в межплитном контакте в 2 раза) Евразийской и Тихоокеанской плит в районе Японской зоны суб-дукции
Центральной Азии. На данном этапе первых тестовых расчетов внутренние свойства плит были заданы однородными. Упругие и вязкие свойства разломных зон отличались в зависимости от их принадлежности зонам коллизии или субдукции, сдвиговым зонам, например, трансформным разломам, либо зонам дивергенции меж-плитных контактов. Ширина зон контактов варьировала от 30 до 40 км. Модуль упругости имел для первых двух типов значения 40 ГПа и для зон дивергенции — 20 ГПа. Вязкость в зонах конвергенции составляла 1020Па- с, а для зон сдвига и дивергенции — 1018^1019Па- с и 1017Па- с соответственно.
В основу первого численного расчета были положены экспериментальные данные о векторах скоростей движения земных плит. Эти данные были взяты из модели REVEL [12]. Было принято условие однородного движения в пределах каждой из плит. После задания годовой подвижки были рассчитаны нормальные и касательные деформации в межплитных контактах в районах основных зон коллизии и субдукции. Значения нормальных составляющих деформации оказались следующими: -1.5 -п10-6 — Японская зона субдукции; -1.0- 10-6 — Южно-Американские зоны субдукции; до -1.0 - 10-6 — Идзу-Бонинская зона субдукции; до -0.9 - 10-6 — Курило-Камчатская зона субдукции; -0.6 -10-6 — Зондская зона субдукции; -0.5 - 10-6 — Индо-Евразийская зона коллизии. Значения касательных к плоскости контакта деформаций варьировали от 0.7 -10-7 в Японской зоне субдукции (минимальное значение) до 3^5 - 10-7 в Идзу-Бонинской, Курило-Камчатской и Южно-Американских зонах суб-дукции. Полученные значения нормальных деформаций зон межплитных контактов отображают потенциальную возможность каждой из зон к скорости накопления упругих деформаций. Различаясь не более чем в 3 раза по
максимальным значениям, полученные данные указывают на главенство по этому параметру Западно-Тихоокеанских и Южно-Американских зон субдукции.
Другой численный эксперимент связан с оценкой эффекта увеличения внутриплитных смещений среды в пределах Евразийской плиты за счет повышения жесткости в межплитном контакте на границе Евразийской и Тихоокеанской плит (Японская зона субдукции). Таким увеличением жесткости (в 2 раза) моделируется частичное сцепление между плитами. Подобные явления в зонах субдукции имеют место, что показано в настоящее время на основе анализа данных GPS-наблюдений [13]. На первом этапе были рассчитаны смещения, вызванные годовой подвижкой Тихоокеанской плиты на северо-запад при обычных значениях жесткости разломных зон. На втором этапе тот же самый расчет был сделан для условий увеличения жесткости в 2 раза на отмеченной выше границе плит в районе Японии. Далее, из величин смещений, вычисленных на втором этапе, были вычтены те значения, которые рассчитывались на первом этапе. Результат, который демонстрирует распределение векторов дополнительного смещения, вызванного частичным сцеплением плит в зоне субдук-ции, представлен на рис. 3. Обращает на себя внимание относительно слабое затухание эффекта с удалением от зоны субдукции на запад. В районе Байкала величина приращения смещений снижается не более чем в 5^10 раз относительно максимальных смещений вблизи границы плит.
4. Заключение
Сформулированы основные положения для построения трехмерной вязкоупругой модели литосферы Земли для численных расчетов смещений, деформаций и напряжений, вызванных как глобальными, так и региональными сейсмогеодинамическими процессами. Разработаны численные схемы, учитывающие сферичность Земли. С учетом разработанной методологии выполнено построение этой модели — разбиение на сеть конечных элементов, которая учитывает геометрическую структуру основных тектонических плит и разлом-ных зон между ними. Модель в первом приближении отражает упругие свойства и реологию основных тектонических элементов, а также изменение этих свойств по вертикали.
Первые численные эксперименты, основанные на данных GPS-наблюдений, во-первых, позволили сделать количественные оценки нормальных и касательных деформаций в межплитных контактах и выявили их отличия в скоростях накопления упругих деформаций
в основных зонах субдукции и коллизии. Во-вторых, получены данные о величине дополнительных смещений в пределах Евразийской плиты, обусловленных частичным сцеплением этой и Тихоокеанской плиты в районе Японской зоны субдукции.
Работа выполнена при поддержке проектов РФФИ №№ 03-05-65418, 06-05-64921 и программ Президиума РАН 13 и 6.5.1.
Литература
1. Дядьков П.Г. Закономерности в изменениях напряженного состояния земной коры и сейсмического режима в Байкальском регионе в условиях влияния процессов на ближайших конвергентных границах плит // Современная геодинамика и опасные природные процессы в Центральной Азии. Иркутск, 2004. - Вып. 1. -С. 174-180.
2. ДядьковП.Г., НазаровЛ.А., НазароваЛ.А. Детальное моделирование поля напряжений земной коры и динамической неустойчивости сейсмоактивных разломов при рифтогенезе // Геология и геофизика. - 1997. - № 12. - С. 2001-2010.
3. Дядьков П.Г., Мельникова В.И., Назаров Л.А., Назарова Л.А., Саньков В.А. Сейсмотектоническая активизация Байкальского региона в 1989-1995 годах: результаты экспериментальных наблюдений и численное моделирование изменений напряженно-деформированного состояния // Геология и геофизика. - 1999. - Т. 40. -№ 3. - С. 373-386.
4. Дядьков П.Г., Назаров Л.А., Назарова Л.А. Трехмерная вязкоупругая модель литосферы Центральной Азии: методология построения и численный эксперимент // Физ. мезомех. - 2004. -Т.7. - № 1. - С. 91-101.
5. Grollimund B., Zoback M.D. Post glacial lithospheric flexure and induced stresses and pore pressure changes in the nothern North Sea // Tectonophysics. - 2000. - V. 327. - P. 61-81.
6. Pollitz F.F., Burgmann R. Joint estimation of afterslip rate and post-seismic relaxation following the Loma Prieta earthquake // JGR. -1998. - V. 103. - No. B11. - Р. 26975-26992.
7. Roy M., Royden L.H. Crustal rheology and faulting at strike-slip plate boundaries. 1. An analitic model // JGR. - 2000. - V. 105. - No. B3.- Р. 5583-5597.
8. 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.
9. McKenzie et al. Characteristics and consequences of flow in the lower crust // JGR. - 2000. - V. 105.- No. B5. - Р. 11029-11046.
10. 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. - Р. 8035-8054.
11. Pollitz F.F., Burgmann R., Romanowich B. Viscosity of oceanic asthenosphere inferred from remote triggering of earthquakes // Science. - 1998. - V. 280. - Р. 1245-1249.
12. Sella G.F., Dixon T.H., Mao A. REVEL: A model for recent plate velocities from space geodesy // J. Geophys. Res. - 2002. - V. 107. -No. B4. - 10.1029/2000JB000033.
13. Mazzotti S., Le Pishon X., Pierre H., Miyazaki S.-I. Full interseismic locking of the Nankai and Japan-West Kurile subduction zones: An analysis of uniform elastic strain accumulation in Japan constrained by permanent GPS // J. Geophys. Res. - 2000. - V. 105. - No. B6. -P. 13159-13177.
