Физика Земли, 2021, № 3, стр. 19-29

ВВЕДЕНИЕ

Одной из наиболее актуальных проблем современной геофизики является предупреждение и снижение последствий таких опасных природных катастроф, как землетрясения, которые вызывают большое число человеческих жертв и приводят к значительному материальному ущербу. Для решения этой проблемы необходимо изучать закономерности сейсмичности в различных сейсмоопасных регионах и выявлять предвестники землетрясений, проявляющиеся в аномалиях различных физических полей [Соболев, Пономарев, 2003], например, в вариациях параметров ионосферы [Бондур, Смирнов, 2005; Бондур и др., 2007], в изменениях полей температуры [Бондур, Воронова, 2012], деформациях линеаментов, регистрируемых по космическим изображениям [Бондур, Зверев, 2005; 2007]. При этом необходимо использовать различные подходы к мониторингу сейсмоопасных территорий [Соболев, Пономарев, 2003; Акопян и др., 2017; Бондур, Смирнов, 2005; Бондур, Зверев, 2005; 2007; Бондур и др., 2007].

Изучению закономерностей сейсмичности посвящено большое количество исследований, в том числе опирающихся на анализ тех или иных геомеханических моделей. В частности, для сильнейших землетрясений, происходящих в субдукционных зонах, одной из наиболее распространенных моделей является модель “асперити” [Lay, Kanamori, 1981]. Эта модель предлагает объяснение различия характерных размеров очагов и периодов повторяемости катастрофических землетрясений, приуроченных к разным зонам субдукции, однако в ней не рассматриваются длительные переходные процессы, которые могут оказывать серьезное влияние на подготовку последующих сейсмических событий. Ограничения этой модели были частично преодолены в клавишно-блоковой модели сейсмотектонических циклов [Lobkovsky et al., 1991], позволяющей описать геодинамическую эволюцию системы нависающих литосферных блоков, установить параметры цикличности сильнейших субдукционных землетрясений и проследить движение блоков на постсейсмической стадии. В последние годы данная модель нашла свое подтверждение прямыми наблюдениями c использованием глобальных навигационных спутниковых систем, в том числе и для Курило-Камчатской зоны субдукции [Лобковский и др., 2017]. Следует отметить, что названные выше модели описывают общий механизм сейсмотектонических процессов без привязки к особенностям конкретного региона.

В то же время для лучшего понимания геодинамического режима и сейсмической опасности той или иной территории могут оказаться полезны детальные геомеханические модели, учитывающие разломно-блоковое строение геологической среды и основные протекающие в ней деформационные процессы. Сейсмическая активность связана с медленными тектоническими движениями земной коры, формирующими геологические структуры и приводящими к накоплению в них значительной упругой энергии, которая разряжается в окружающее пространство в результате разрушения материала земной коры в тех ее местах, где тектонические напряжения достигают предела прочности естественным образом или в результате вмешательства человека [Триггерные эффекты…, 2010]. Поэтому в основе всех видов геомониторинга должен лежать анализ напряженно-деформированного состояния земной коры исследуемого района, которое складывается под действием многих факторов. Среди основных можно назвать горизонтальные и вертикальные тектонические движения, нескомпенсированный вес горного рельефа, вариации плотности, неравномерный нагрев. Перечисленные виды нагрузок взаимодействуют со сложной геофизической средой, главным признаком которой является разломно-блоковое строение, вертикальная реологическая расслоенность, наличие флюидов и горного рельефа.

Изучение сейсмотектонических процессов на Камчатке имеет важное практическое значение, связанное с оценкой сейсмической опасности и прогнозом землетрясений. Хотя сейсмичность этого региона преимущественно связана с зоной субдукционного слэба, а основная масса гипоцентров событий в пределах полуострова и у его Тихоокеанского побережья локализуется на достаточно большой глубине (свыше 40–50 км), имеются также проявления мелкофокусной сейсмичности, ассоциируемой с разломно-блоковой тектоникой верхних горизонтов земной коры [Гордеев и др., 2006]. В этой связи представляется, что для мониторинга и прогноза опасных геодинамических процессов должна использоваться подробная геомеханическая модель, отражающая разломно-блоковое строение и основные закономерности современной тектоники региона, обеспечивая возможность изучения процессов перераспределения напряженного состояния и выделения зон аккумуляции напряжений.

Настоящая работа посвящена созданию трехмерной (3D) геомеханической модели земной коры Камчатки. Построение модели выполнялось на основе данных о структуре коры и геодинамических процессах, разломной тектоники, неоднородного строения фундамента и осадочной толщи. Такая модель дает возможность рассчитать напряженно-деформированное состояние региона. Знание напряженного состояния земной коры позволяет количественно оценить такие важные параметры, как сдвиговые напряжения и энергонасыщенность, отвечающие за сейсмическую активность [Гарагаш, 1991]. Полученные результаты в той или иной степени позволяют выделить области, в которых могут произойти сильные землетрясения, аналогично тому, как это было сделано при анализе напряженно-деформированного состояния Прикаспия [Лобковский и др., 2013а]. Также они могут послужить основой для организации системы наземно-космического мониторинга напряженного состояния участков земной коры на сейсмоопасных территориях.

ГЕОМЕХАНИЧЕСКАЯ МОДЕЛЬ ЗЕМНОЙ КОРЫ КАМЧАТКИ

Авторами ранее был проведен цикл исследований сейсмичности Южной Калифорнии на основе ежемесячного расчета динамики параметров напряженного состояния в рамках детальной геомеханической модели Южной Калифорнии в работах [Бондур и др., 2010; 2016а; 2016б; 2020; Bondur et al., 2017]. Входными данными для модели является текущая сейсмичность (магнитуды событий) в движущемся с заданным шагом по времени трехмесячном окне. Каждое землетрясение рассматривается как новый дефект земной коры, приводящий к изменению эффективных механических свойств среды, что вызывает изменение напряженного состояния, перераспределение деформаций и плотности упругой энергии. Для оценки динамики этого процесса предложен параметр, являющийся мерой близости горных пород к пределу прочности. Анализ аномалий данного параметра позволил выделить предвестники ряда крупных сейсмических событий (M > 7), произошедших на данной территории в 2010 и 2019 гг. Вышеописанная технология и результаты ее внедрения, начиная с 2009 года, подробно обсуждаются в работах [Бондур и др., 2010; 2016а; 2016б; 2020; Bondur et al., 2017].

При построении геомеханической модели земной коры Камчатки исследовалась область, включающая юго-восточную часть Камчатского полуострова (Охотоморская литосферная плита), северо-западный край Тихоокеанской литосферной плиты, зону сочленения Курило-Камчатской и Алеутской зон субдукции. В полярной стереографической проекции область представляется прямоугольником с характерными размерами 700 на 850 км (рис. 1). Трансформация координат структурных элементов модели в рамках используемой проекции осуществлялась со следующими параметрами: осевой меридиан – 150° E; базовая широта – 70° N; масштабный фактор на базовой широте – 1; единицы измерения – метры; эллипсоид – WGS84. Прямоугольная сетка модели состоит из порядка 100 × 100 ячеек по горизонтали и 9 слоев по вертикали. Горизонтальные размеры ячеек составляют 5 × 5 км.

По глубине область моделирования разделялась на три слоя: кора, надастеносферная мантия и астеносфера. Верхней граничной поверхностью является дневная поверхность (рис. 2а). Граница Мохо построена по материалам ФГУП “ВНИИОкеангеология им. И.С. Грамберга” (рис. 2б).

Поверхность астеносферы для изучаемого региона рассчитывалась на основании данных о тепловом потоке с учетом внутренней структуры коры и тепловых свойств образующих пород, в предположении, что вертикальный тепловой поток через земную кору постоянен. Для этого использовалась база данных теплового потока Геологического института РАН, дополненная данными из базы Global Heat Flow Database (https:// ihfc-iugg.org/products/global-heat-flow-database/data) (рис. 3а). Предполагалось, что температура на границе астеносферы равна температуре солидуса 1200°С. В роли второго граничного условия выступает дневная поверхность с фиксированной температурой на ней 0°С. В результате расчета определялась мощность литосферы в регионе, положение поверхности астеносферы (рис. 3б) и распределение температуры. Следует отметить, что построенная таким образом карта поверхности астеносферы является не слишком подробной и содержит погрешности, обусловленные недостаточной детальностью глобальных моделей теплового потока, но в целом отражает закономерный переход от мощной литосферы со стороны Евразии к более тонкой со стороны Тихого океана.

На рис. 4 показана окончательная расчетная геомеханическая модель Камчатского региона с включением границы астеносферы, созданная в программном коде FLAC^3D [Itasca, 2006].

Считается, что в общем случае горные породы ведут себя как упругопластические среды с предельным условием Друкера–Прагера:

где: $\sigma $ и $\tau $ – cоответствующие инварианты тензора напряжений и его девиаторной частиж; ${{q}_{\varphi }}$ – коэффициент трения; ${{k}_{\varphi }}$ – сцепление; условие (1) описывает поверхность в пространстве напряжений.

Выбранные для расчетов значения геомеханических параметров определяются инвариантными характеристиками напряженно-деформированного состояния и их комбинациями, без которых невозможно объяснить особенности рельефа, разломной тектоники и особенности возможного разрушения среды [Качанов, 1969]:

Кора: плотность $\rho = 2875\,\,{{{\text{кг}}} \mathord{\left/ {\vphantom {{{\text{кг}}} {{{{\text{м}}}^{{\text{3}}}}}}} \right. \kern-0em} {{{{\text{м}}}^{{\text{3}}}}}}$, модуль сдвига G = 2.085 × 10¹⁰ Па и объемный модуль K = 6.254 × 10¹⁰ Па, коэффициент трения ${{g}_{\phi }} = 0.64$, сцепление ${{k}_{\phi }} = 6 \times {{10}^{7}}$ Па и прочность на растяжение ${{\sigma }^{t}} = 3 \times {{10}^{7}}$ Па.

Литосфера: плотность $\rho = 3225\,\,{{{\text{кг}}} \mathord{\left/ {\vphantom {{{\text{кг}}} {{{{\text{м}}}^{{\text{3}}}}}}} \right. \kern-0em} {{{{\text{м}}}^{{\text{3}}}}}}$, модуль сдвига G = 3.42 × 10¹⁰ Па и объемный модуль K = = 1.025 × 10¹¹ Па, коэффициент трения ${{g}_{\phi }} = 0.1$, сцепление ${{k}_{\phi }} = 4 \times {{10}^{8}}$ Па и прочность на растяжение ${{\sigma }^{t}} = 2 \times {{10}^{8}}$ Па.

В силу того, что материал астеносферы частично расплавлен и на больших временах испытывает вязкие деформации, астеносфера моделируется вязкоупругой средой Максвелла со следующими параметрами: плотность $\rho = 3225\,\,{{{\text{кг}}} \mathord{\left/ {\vphantom {{{\text{кг}}} {{{{\text{м}}}^{{\text{3}}}}}}} \right. \kern-0em} {{{{\text{м}}}^{{\text{3}}}}}}$, модуль сдвига G = 3.42 × 10¹⁰ Па и объемный модуль K = = 1.025 × 10¹¹ Па, вязкость η = 1.5 × 10²² Па ⋅ с.

Для решения задачи о напряженно-деформированном состоянии земной коры важно учитывать ее разломно-блоковое строение. При введении в модель разломов мы исходим из того, что область пониженной прочности земной коры не ограничивается поверхностью того или иного разрывного нарушения, а распространяется по всем трем измерениям в некоторой его окрестности. Поврежденность характеризуется функцией неоднородности $g({{x}_{s}})$, меняющейся от единицы до нуля, что позволяет задать все механические параметры в виде:

где ${{\Pi }^{0}}$ – однородное начальное значение параметра для ненарушенной среды, κ ≤ 1 – параметр малости.

Расчет функции неоднородности предполагает проекцию карты разломов на сетку модели в плане. При этом максимальные значения поврежденности присваивались тем ячейкам области моделирования, которые непосредственно пересекаются разломами. На следующем шаге выполнялось пространственное сглаживание полученного распределения поврежденности, радиус которого возрастает для более глубоких слоев модели. Для определения функции неоднородности использовались схемы разломной тектоники, приведенные в статьях [Kozhurin et al., 2006; Avdeiko et al., 2007]. Соответствующее нормированное распределение поврежденности среды $g({{x}_{s}})$ показано на рис. 5.

Для изучаемого района характерна сильная сейсмичность. На карту разломов, приведенную на рис. 5, по данным каталога USGS ComCat (https:// earthquake.usgs.gov/earthquakes/search/#earthquakes) вынесены эпицентры землетрясений с магнитудой $M \geqslant 6.0$, известные для этого региона с 1904 г. Следует отметить, что большая часть гипоцентров сильных землетрясений локализуется в пределах зоны Беньофа и связана с процессом субдукции Тихоокеанской литосферной плиты. Таким образом, основная сейсмическая опасность в регионе связана с сильными субдукционными землетрясениями, происходящими на расстояниях около 100 км от Тихоокеанского побережья Камчатки, способными вызывать цунами. В то же время, на полуострове также имеются проявления мелкофокусной сейсмичности, связанной с тектоническими нарушениями в верхней части земной коры [Гордеев и др., 2006].

Регион находится под действием горизонтальных тектонических усилий, порождаемых движениями тектонических плит. Для моделирования движения плит используется рекомендованная Международной службой вращения Земли (IERS) модель литосферной кинематики NNR-NUVEL-1A. При переводе скоростей движения плит из глобальной (сферической) модели в полярную стереографическую проекцию (на плоскость) учитывалось изменение соответствующих азимутальных углов. На рис. 5 показано окончательное распределение скоростей, принятое для расчетов, полученное путем интерполяции значений глобальной модели в пределах области моделирования.

Основные закономерности современных тектонических движений в рассматриваемом регионе определяются процессом субдукции Тихоокеанской литосферной плиты в северо-западном направлении с характерными скоростями около 5 см/год. В зоне сочленения Курило-Камчатской и Алеутской островных дуг движение Тихоокеанской плиты происходит практически по касательной по отношению к Североамериканской плите, движущейся в этой зоне в южном направлении, что выражается в существовании зоны сдвига в районе Командорских островов.

Необходимо отметить, что обсуждаемая в настоящей работе модель не содержит явного структурного элемента, отражающего Курило-Камчатскую зону субдукции (поверхности слэба), однако учитывает ее наличие опосредованно путем задания рельефа границ Мохо и литосфера–астеносфера, а также современных скоростей движения литосферы. При анализе напряженно-деформированного состояния мы в основном ограничивались рассмотрением литосферной части модели, не принимая во внимание конвекционные движения в верхней мантии [Лобковский и др., 2013б]. Тем не менее, включение в модель горизонтальных скоростей литосферных плит в некотором приближении отражает циркуляцию вязкоупругой астеносферы под Охотоморской плитой, что выражается в направленном движении последней на юго-восток, в сторону Курило-Камчатской зоны субдукции со скоростями порядка 1–2 см/год вследствие действия возвратного потока верхнемантийной ячейки, питаемой субдуцирующей литосферой Тихоокеанской плиты.

Согласно имеющимся геологическим и геофизическим данным [Баранов и др., 2010; Kozhurin et al., 2006], большинство разломов продольного направления (субпараллельных оси полуострова) являются сбросами, что отражает поперечное растяжение в данной области, по крайней мере в верхней части земной коры. Для поперечных разломов характерны сдвиговые смещения, что, видимо, связано с клавишно-блоковой структурой внешнего края островодужной гряды.

РАСЧЕТ НАПРЯЖЕННО-ДЕФОРМИРОВАННОГО СОСТОЯНИЯ

Для расчетов использовался программный код FLAC^3D [Itasca, 2006], реализующий явную конечноразностную схему решения трехмерных задач механики сплошных сред.

Напряженно-деформированное состояние земной коры складывается под действием многих факторов. В расчете учитывались основные: некомпенсированный вес горного рельефа, горизонтальные и вертикальные тектонические движения, вариации плотности, неравномерный нагрев.

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

2. Рассчитывались напряжения в коре с учетом распределения температуры. Так как свойства твердой среды в значительной степени определяются температурой, на втором этапе учитываются напряжения, вызванные температурным расширением.

3. Выполнялась нагрузка горизонтальными тектоническими усилиями, порождаемыми движениями тектонических плит. Для этого использовалась стандартная модель распределения скоростей движения плит NNR-NUVEL-1A.

Все расчеты выполнялись с учетом ползучести астеносферы.

После того, как были произведены расчеты, проводился анализ ряда параметров напряженно-деформированного состояния. Прежде всего, это инварианты тензоров напряжений – среднее напряжение

и интенсивность сдвиговых напряжений где ${{s}_{{ij}}} = {{\sigma }_{{ij}}} - {{\delta }_{{ij}}}\sigma $ – компоненты девиатора напряжений.

Заметим, что второй инвариант напряжений практически совпадает по величине с максимальным касательным напряжением.

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

Накопленную объемную энергию можно вычислить по формуле:

а упругую энергию сдвига по формуле:

Здесь: G – модуль сдвига; $K = \frac{{2G(1 + \nu )}}{{3(1 - 2\nu )}}$ – модуль объемного сжатия; ν – коэффициент Пуассона.

Еще один полезный коэффициент – это параметр Лоде–Надаи, который позволяет оценить к какому типу напряженного состояния ближе всего полученное распределение напряжений. Параметр Лоде–Надаи вычисляется по формуле:

и связан со вторым и третьим инвариантами девиатора напряжений.

Параметр μ_σ меняется в пределах от –1 до +1. Так, для чистого растяжения (${{\sigma }_{1}} > 0$, ${{\sigma }_{2}} = {{\sigma }_{3}} = 0$) ${{\mu }_{\sigma }} = 1$, для чистого сжатия (${{\sigma }_{1}} = {{\sigma }_{2}} = 0$, ${{\sigma }_{3}} < 0$) ${{\mu }_{\sigma }} = - 1$, для чистого сдвига (${{\sigma }_{1}} > 0$, ${{\sigma }_{2}} = 0$, ${{\sigma }_{3}} = - {{\sigma }_{1}}$) ${{\mu }_{\sigma }} = 0$.

Пластическое течение происходит при выполнении предельного условия Друкера–Прагера (1).

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Полученные в результате моделирования параметры напряженно-деформированного состояния представлены на серии иллюстраций (рис. 6).

Согласно полученному распределению параметра Лоде–Надаи, в верхних горизонтах земной коры имеет место сложная картина напряжений, отражающая взаимодействие блоков и обнаруживающая ряд областей растяжения. В более глубоких горизонтах земной коры характер этой картины меняется, отражая разнообразное сочетание зон растяжения, сжатия и сдвига (рис. 6б).

На рис. 6в приводится распределение плотности потенциальной энергии сдвига, а на рис. 6г – плотности потенциальной энергии изменения объема в земной коре на глубине 25 км.

В диапазоне глубин 25 км основная часть территории полуострова находится в условиях сжатия, что, видимо, обусловлено конвергентным характером взаимодействия Охотоморской и Тихоокеанской плит. При этом наибольшие значения плотности объемной энергии сжатия локализуются в центральной и восточной частях полуострова в области сочленения ряда глубинных разломов. Максимумы плотности энергии сдвига концентрируются в центральной части полуострова (Центральная Камчатская депрессия). Можно отметить, что энергия сдвига (рис. 6в) имеет наибольшие значения в центральной части полуострова; в то же время, согласно распределению параметра Лоде–Надаи (рис. 6б), эта область соответствует условиям сжатия. В пределах полуострова выделяется ряд зон сдвига: в его юго-восточной части, в окрестности разломов поперечного простирания, а также в зоне продольных разломов, ограничивающих с востока Центральную Камчатскую депрессию. В зоне глубоководного желоба на рассматриваемой глубине также прослеживаются условия сдвига. Область наиболее интенсивного растяжения отмечается к северу от места сочленения Алеутской и Курило-Камчатской зон субдукции.

На рис. 7 показано распределение плотности потенциальной энергии сдвига и параметра Лоде–Надаи в литосфере, с нанесенными на них землетрясениями с магнитудой $M \geqslant 6.0$, известными с 1904 г. Можно отметить, что области максимальной энергонасыщенности литосферы локализуются в окрестности зоны субдукции как со стороны Охотоморской, так и со стороны Тихоокеанской плиты. Повышенные значения энергии наблюдаются непосредственно под границей Мохо в южной части п-ва Камчатка, а также вблизи западной окраины Алеутской островной дуги в пределах Тихоокеанской литосферной плиты. Анализ показывает, что землетрясения тяготеют к зонам с повышенной энергонасыщенностью (рис.7а) и преобладанием напряжений сдвига и сдвига со сжатием (рис.7б). Эти признаки выделяют области, в которых могут произойти сильные землетрясения.

ЗАКЛЮЧЕНИЕ

В рамках настоящей работы построена трехмерная геомеханическая модель литосферы Камчатки и выполнен анализ напряженного состояния земной коры, вызванного действием сил собственного веса, температурой и движением тектонических плит. Установлено, что картины энергонасыщенности и параметра Лоде–Надаи коррелируют с распределением сильной сейсмичности в регионе. Так, из распределения параметров напряженного состояния в литосфере видно, что зоны ее максимальной энергонасыщенности локализуются как со стороны Охотоморской, так и со стороны Тихоокеанской плиты. Повышенные значения энергии наблюдаются непосредственно под границей Мохо в южной части п-ва Камчатка, а также вблизи западной окраины Алеутской островной дуги в пределах Тихоокеанской литосферной плиты. Можно сделать вывод, что землетрясения тяготеют к зонам с повышенной энергонасыщенностью и преобладанием напряжений сдвига и сдвига со сжатием. Эти признаки выделяют области, в которых могут произойти сильные землетрясения.

Построенная геомеханическая модель сможет быть использована для оценки изменений напряженного состояния земной коры, связанных с развитием сейсмического процесса. В рамках подобного исследования модель будет последовательно корректироваться на основе данных сейсмического каталога, при этом каждое землетрясение рассматривается как новый дефект земной коры, приводящий к изменению механических параметров среды и, как следствие, к изменению напряженно-деформируемого состояния, определяемому итерационно по приращениям деформаций и напряжений в течение заданного промежутка времени. Такой подход позволит сделать вывод о том, на каких участках земная кора приближается к пределу прочности и следует ожидать сильного землетрясения [Бондур и др., 2016; Bondur et al., 2017; 2020].