Физика Земли, 2023, № 4, стр. 16-28

ВВЕДЕНИЕ

Вопросы краткосрочного прогноза сильных землетрясений активно разрабатываются с использованием различных подходов, включая технологии дистанционного зондирования [Акопян и др., 2017; Бондур и др., 2007]. Важным звеном в понимании динамики напряженно-деформированного состояния среды перед катастрофическими событиями является разработка моделей разрушения земной коры, направленная на изучение таких особенностей развития данного процесса, которые могли бы рассматриваться как краткосрочные предвестники.

Согласно недавнему обзору [Kato, Ben-Zion, 2021] для объяснения процесса подготовки крупных землетрясений используются три основные модели: каскадно-восходящая (cascade-up), предсейсмического скольжения (pre-slip) и прогрессирующей локализации (progressive localization). Каскадно-восходящая модель предполагает, что подготовка землетрясения связана с прогрессирующим разупрочнением горных пород, вызванным слабой сейсмичностью на протяжении длительного интервала времени в пределах системы ранее существовавших разломов, подверженных фоновым тектоническим напряжениям [Yoon et al., 2019]. В рамках модели предсейсмического скольжения основное внимание уделяется медленным асейсмическим движениям, происходящим в конкретном разломе перед крупным сейсмическим событием [Dieterich, 1992; Кочарян, 2014], а концепция прогрессирующей локализации описывает динамику распределенных напряжений и деформаций, развивающихся таким образом, что перед сейсмическим вспарыванием возникает аномалия деформации сдвига в зоне, не связанной с ранее существовавшими разломами [Rice, 2017; Ben-Zion, Zaliapin, 2020].

Механическая устойчивость сейсмически активных зон с уже имеющимися разрывными нарушениями рассматривается, например, в моделях лавинно-неустойчивого трещинообразования (ЛНТ) и диффузионно-дилатантной концепции (ДД) [Мячкин и др., 1975]. При этом, процесс сводится к разрушению так называемой перемычки, ограниченной с двух сторон начавшимся образованием трещин, что находит отражение в лабораторных экспериментах [Соболев, Пономарев, 2003].

В работе [Гарагаш, 2000] показано, что разрушение перемычки перед Кроноцким землетрясением на Камчатке связано с появлением аномалии гантелевидной формы в касательных напряжениях. Эффективность применения геомеханического моделирования на примере Южной Калифорнии показана в работах [Bondur et al., 2020; Бондур и др., 2010; 2016; 2021], где выявлены краткосрочные предвестники в динамике НДС сильных землетрясений 2010 и 2019 гг. с М = 7.1.

Другое направление в изучении природы сейсмичности связано с конструированием детерминированных механических моделей, используемых для описания процессов в разломно-блоковых структурах земной коры [Burridge, Knopoff, 1967]. Эти модели представляют собой набор отдельных элементов (блоков), лежащих на движущейся плите и соединенных между собой пружинными жесткими связями. Преодоление сил трения на контакте с плитой способствует проскальзыванию отдельных блоков. Такое прерывистое скольжение (stick-slip) и порождает сейсмические события. Д. Туркотт [Turcotte, 1994] провел обзор большого числа детерминированных механических моделей, эволюция которых во времени использовалась широким кругом исследователей для описания такого в целом хаотического природного процесса, как сейсмичность [Huang, Turcotte, 1990; Narkounskaia, Turcotte, 1992; Brown et al., 1991; Carlson et al., 1993].

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

Основой создания модели “Качели” являются результаты, полученные по динамике напряженно-деформированного состояния (НДС) земной коры Южной Калифорнии перед землетрясением Риджкрест 6.07.2019 г. с М = 7.1 и координатами 35.766° с.ш. и 117.605° з.д.

Анализ проводился в рамках пятислойной геомеханической модели, учитывающей разломную тектонику, реологические свойства слоев земной коры, рельеф и внешние тектонические воздействия, что подробно описано в работах [Bondur et al., 2020; Бондур и др., 2010; 2016; 2021]. В плане область моделирования имеет размер 645 × 560 км (31–36° с.ш. и 114–121.2° з.д.) и дискретизируется прямоугольной сеткой с элементами (ячейками) 5 × 5 км. Модель охватывает интервал глубин от дневной поверхности (с учетом рельефа) до 35 км; разбиение в вертикальном направлении увязано с геометрией границ слоев, задаваемых на основе имеющихся геологических данных. Каждому слою присвоены определенные значения геомеханических параметров в рамках упругопластической модели Кулона−Мора (модуль объемного сжатия K, модуль сдвига G, сцепление c и угол внутреннего трения ${{\varphi }}$). Разломно-блоковая тектоника региона задается с использованием данных линеаментного анализа [Bondur et al., 2010]. Распределения некоторых параметров показаны на рис. 1.

На первом этапе геомеханического моделирования выполнялся расчет начального напряженного состояния модели, отражающего ее начальную поврежденность, связанную с разломно-блоковой тектоникой, учитывающего собственный вес и действие региональных тектонических сил. Вычисления проводились в квазистатическом приближении с использованием программного пакета FLAC3D [Itasca…, 2006].

Последующее моделирование предполагает учет текущей сейсмической активности и реализуется в виде итерационной процедуры пересчета текущего состояния (НДС) модели, осуществляемой путем последовательной коррекции значений геомеханических свойств (модулей объемного сжатия, сдвига, сцепления и угла внутреннего трения) тех элементов (ячеек) модели, которые оказались затронуты сейсмическими событиями с магнитудами M > 1 (по данным каталога USGS ComCat, https://earthquake.usgs.gov/data/comcat/), в трехмесячный временной интервал, предшествующий моменту расчета НДС. Корректировка осуществляется на основе оценок выделившейся энергии землетрясений и максимального касательного напряжения. Скорректированная модель используется для расчета нового НДС. Затем описанный цикл повторяется для нового трехмесячного интервала, сдвинутого на 0.5 мес. относительно предыдущего (рис. 2).

Таким образом, “работа” модели в режиме мониторинга сводится к последовательному расчету НДС с полумесячным интервалом по времени. Полученные массивы параметров НДС используются для расчета и визуализации ключевых кинематических характеристик – приращения вектора смещения $u = \left\{ {{{u}_{x}},{{u}_{y}},{{u}_{z}}} \right\}$, рассчитываемого во всех узлах сетки, а также интенсивности деформации сдвига (сдвиговой деформации – СД) ${{\gamma }}$:

где ${{\varepsilon }_{{ij}}}$ – компоненты тензора деформации.

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

Динамика НДС исследовалась при введении в модель текущей сейсмичности, когда каждое землетрясение, начиная с М > 1 по каталогу USGS рассматривалось как новый дефект земной коры соответствующих размеров, приводящий к перераспределению деформаций [Бондур и др., 2016].

Используемые для построения модели результаты относятся к пространственно-временным распределениям приращений деформации сдвига (СД) и векторов смещения на глубине 3–10 км второго слоя верхней земной коры, где локализуется основная текущая сейсмичность Южной Калифорнии.

На рис. 3. представлено распределение максимальных значений приращения СД и смещений с двух сторон от будущего разрыва землетрясения Риджкрест 06.07.2019 г. за 1.5 мес. и 5 дней, соответственно. Размеры областей составляют порядка 20 км при детальности по размерам расчетных ячеек в модели 5 × 5 км, а величины приращения смещений за 0.5 мес. достигают порядка десятков сантиметров.

Полученное распределение в окрестности реального очага находится в согласии с теоретической моделью развития трещины из работы [Мячкин и др., 1975] (рис. 4).

В работе [Бондур и др., 2022] приводятся результаты по динамике НДС в эпицентральной зоне (ЭЗ) на 2018–2019 гг. примерно за 1.5 года до землетрясения, когда началось интенсивное формирование очага (рис. 3). Как видно, аномальное развитие приращения смещений в СД с двух сторон от формирующегося очага принимает гантелевидную форму готовой к разрушению перемычки.

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

В соответствии с рис. 5 и рис. 6 приведено последовательное пространственное распределение характерных моментов максимальных смещений, начиная с марта 2018 г., когда, как показано в работах [Бондур и др., 2021; 2022], проявляется интенсивное формирование очага. Попеременное появление максимальных смещений на разных оконечностях будущего разрыва имитирует процесс раскачивания эпицентральной зоны (ЭЗ) вплоть до момента события.

Детальное распределение всех величин приращения смещений за данный промежуток времени приведено на рис. 7 [Бондур и др., 2022].

Отрицательные величины, показанные красным цветом, соответствуют местоположению аномальных смещений с ЮВ стороны продолжения простирания очага, а положительные (синего цвета) с СЗ стороны в соответствии с рис. 5.

Таким образом, приведенные результаты указывают на наличие процесса раскачивания области, готовящейся к разрушению перемычки. При этом, период раскачивания таких “качелей” составляет порядка одного–нескольких месяцев, протяженность порядка 100 км, размеры аномальных областей – в пределах 20 км. Величины приращения смещений, моделируемых на глубинах 3–10 км второго слоя верхней земной коры в геомеханической модели, достигают десятков см.

МОДЕЛЬ “КАЧЕЛИ”

Согласно современным представлениям, земная кора представляет собой иерархическую дискретную среду [Садовский, Писаренко, 1991]. Такое строение неоднородной среды определяет развитие сейсмотектонического процесса. Подтверждением этого вывода является модель разлома [Burridge, Knopoff, 1967], объясняющая эволюцию сейсмического процесса (рис. 8).

Модель представляет собой набор отдельных элементов (блоков) с массой $m$, лежащих на неподвижной горизонтальной плите и соединенных между собой пружинными жесткими связями ${{k}_{c}}$. С подвижной горизонтальной плитой, движущейся со скоростью $V$, элементы соединены пластинчатыми элементами с жесткостью ${{k}_{l}}$. Преодоление силы трения ${{F}_{i}}$ между нижней плитой и блоками приводит время от времени к проскальзыванию отдельных блоков. Проскальзывать может как отдельный блок, так и несколько блоков. Такое прерывистое скольжение (stick-slip) и возникающие при этом фрикционные автоколебания связывают с возникновением землетрясений [Короновский и др., 2015]. Примером эффективности такого подхода может служить клавишная модель сильных землетрясений в островных дугах и активных континентальных окраинах [Керчман, Лобковский, 1988; Лобковский и др., 2017].

Оставаясь в изложенной парадигме, рассмотрим следующую модель движения блоков, объясняющую раскачивание перемычки перед разрушением. Рассмотрим движение трех блоков, размещенных на движущейся со скоростью $V$ плите (рис. 9а).

Средний блок представляет собой перемычку. Блоки скреплены между собой упругими связями с коэффициентами ${{k}_{2}}$, ${{k}_{3}}$ и вязкими элементами с коэффициентом вязкости $\mu $. Крайние блоки с помощью упругих элементов с жесткостью ${{k}_{1}}$, ${{k}_{4}}$ взаимодействуют с неподвижными преградами. Блоки прижаты к плите усилиями $P$ и взаимодействуют с ней по закону сухого трения (рис. 7б), согласно которому при достижении усилием $F$ на контакте критического значения ${{F}_{{up}}}$ сила трения уменьшается до значения ${{F}_{{dw}}}$, и под действием сброса нагрузки

начинается проскальзывание.

Предположим, что сначала приходит в движение 1-й блок (рис. 10а), который проскальзывает на величину ${{u}_{1}}$ под действием сброса нагрузки $\Delta {{F}^{{(1)}}} = F_{{up}}^{{(1)}} - F_{{dw}}^{{(1)}}$.

По мере того, как нагрузка $\Delta {{F}^{{(1)}}}$ передается на вязкоупругую связь между блоками, начинает возрастать нагрузка на средний 2-й блок (перемычку). Поскольку $F_{{up}}^{{(2)}} > F_{{up}}^{{(1)}}$ и $F_{{up}}^{{(2)}} > F_{{up}}^{{(3)}}$, это не приведет к активизации перемычки. Поэтому через некоторое время растяжение упругой связи ${{k}_{4}}$ приведет в движение 3-й блок, и он сместится на величину ${{u}_{3}}$. В результате к нагрузке на перемычку со стороны 1-го блока добавится давление 3-го блока. Если этого будет достаточно для преодоления силы трения $F_{{up}}^{{(2)}}$, то 2-й блок стронется с места. Если нет, то опять повторится нагружение 1-го блока и дальнейшее возрастание нагрузки на перемычку. Описанный выше процесс будет продолжаться до тех пор, пока не будет преодолено сопротивление 2-го блока.

Выпишем уравнение движения 1-го блока:

Его решение, удовлетворяющее условию ${{u}_{1}} = 0$ при $t = 0$, имеет вид:

Движение 1-го блока остановится, когда нагрузка слева на 1-й блок уменьшится:

Согласно (4) время остановки $t_{1}^{{st}}$ (stick) найдем из условия

Движение 3-го блока удовлетворяет уравнению

аналогичному уравнению (2).

Используя решения уравнений (2) и (6), из условия

найдем время $t_{3}^{{sl}}$ (slip), при котором начнется движение 2-го блока.

Зададимся следующими параметрами модели. Модуль Юнга блоков – ${{E}_{1}} = {{E}_{4}} = 2 \times {{10}^{{10}}}$ Па, ${{E}_{3}} = 1.5 \times {{10}^{{10}}}$ Па, коэффициент Пуассона $\nu = 0.25$, вязкость – ${{\eta }_{1}} = {{\eta }_{2}} = {{\eta }_{3}} = 6 \times {{10}^{{17}}}$ Па ⋅ с, среднее напряжение – ${{\sigma }_{{33}}} = 1.25 \times {{10}^{8}}$ Па, ${{\alpha }_{{up}}} = 0.5$ и ${{\alpha }_{{dw}}} = 0.3$ – коэффициенты трения. Длина и ширина блоков $L = 1 \times {{10}^{5}}$ м и $h = 5000$ м соответственно. Скорость движения плиты $V = 0.5\,\,{{\text{м}} \mathord{\left/ {\vphantom {{\text{м}} {{\text{год}}}}} \right. \kern-0em} {{\text{год}}}}$.

В этом случае время активации (проскальзывания) 1-го блока $t_{{sl}}^{{(1)}} = 0.358$год. Затем происходит активация 3-го блока, которая длится $t_{{sl}}^{{(3)}} = 0.339$ год. Нагрузка, возникающая со стороны 1-го и 3-го блока, приводит в движение 2-й блок (перемычку). Активация 2-го блока длится $t_{{sl}}^{{(2)}} = 0.373$ год. Таким образом, использование характерных значений механических и геометрических параметров системы коровых блоков приводит к обоснованным значениям промежутков времени их раскачивания.

В рамках приведенной модели, используя наблюдаемые (по геодезическим данным) и/или смоделированные (с учетом сейсмичности) смещения (или их скорости) для ряда выделенных блоков, путем решения обратной задачи можно уточнить значения параметров, включающих жесткости k и коэффициенты вязкости $\mu $ для системы блоков таким образом, чтобы ее наблюдаемое поведение соответствовало экспериментальным параметрам, включая период раскачивания и расстояние между фокусами раскачки. На основании полученных численных оценок установлено, что по одним лишь величинам периода раскачивания достаточно устойчиво определяются только межблоковые вязкости $\mu $, в то время как для других геомеханических параметров системы блоков характерна неоднозначность в достаточно широком интервале значений.

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

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

Подход, использованный при построении исходной геомеханической модели Южной Калифорнии [Bondur et al., 2020; 2022], частично включает как элементы каскадно-восходящей модели [Yoon et al., 2019], так и концепции прогрессирующей локализации [Rice, 2017; Ben-Zion, Zaliapin, 2020]. В данном случае поле фоновых напряжений и деформаций, определяемое региональными тектоническими силами и ранее существовавшим распределением разломов, дополняется локальным ослаблением (поврежденностью) и последующим заживлением земной коры, связанными с текущей слабой сейсмичностью, что позволяет моделировать динамические поля напряжений и деформаций для всей области Южной Калифорнии. Фоновая сейсмичность, отражающая непрерывное разрушение земной коры под действием региональных тектонических сил, в основном, связана с разломами Сан-Андреас, Гарлок и Маунт. Однако, поскольку рассматриваемая модель включает данные по сейсмичности малой магнитуды, преобразованные в коэффициент поврежденности горных пород, а фактическая (рассеянная) сейсмичность не ограничивается конкретными разломами (хотя они играют важную роль в распределении стационарных напряжений), расчет эволюции распределенных напряжений и деформаций путем итеративного обновления геомеханических свойств (модулей упругости и т. д.) обеспечивает изменение НДС в достаточно широкой области вокруг будущего разрыва.

Анализ в рамках модели прогрессирующей локализации для Южной Калифорнии, приведенный в работах [Ben-Zion, Zaliapin, 2019; 2020; Kato, Ben-Zion, 2021], свидетельствует о прогрессирующем ослаблении коры вокруг возможных зон вспарывания перед крупными (M > 7) событиями, произошедшими в этом регионе за последние десятилетия, с четкой локализацией в течение последних 2–3 лет перед ними. Вышеупомянутые исследования в основном сосредоточены на прямом преобразовании наблюдаемой сейсмичности в объем и локализацию поврежденности горных пород с применением геостатистических оценок [Molchan, Keilis-Borok, 2008] без использования полного геомеханического моделирования для получения распределений НДС.

Описываемая в настоящей работе модель опирается на предложенную нами ранее концепцию, сочетающую принципы каскадно-восходящей модели и прогрессирующей локализации с геомеханическим моделированием параметров НДС с учетом фоновой сейсмичности, позволившую выявить циклический процесс миграции аномальных деформаций на флангах очага готовящегося сильного землетрясения Риджкрест. С другой стороны, модель “Качели” содержит элементы, близкие к моделям медленного скольжения [Кочарян, 2014].

Выполненные количественные оценки показывают, что обнаруженные периоды процесса “раскачивания” эпицентральной зоны совпадают с периодичностью срыва блоков в рассмотренной детерминированной механической модели. Данный процесс т.н. “раскачивания” является признаком заключительного этапа формирования очага перед началом разрушения перемычки, определяющей само землетрясение. При этом, размеры области “раскачивания” (длина разрыва L) связаны с магнитудой M будущего события, которая может быть оценена, например, в соответствии с формулой Касахары [Касахара, 1985]:

где L задается в см. В случае землетрясения Риджкрест, когда наблюдаемый размер перемычки по данным геомеханического мониторинга составляет порядка 100 км, можно получить оценку максимальной ожидаемой магнитуды около 7.6, что несколько выше фактической. В то же время, вслед за авторами исследования по анализу поверхностных смещений [Chen et al., 2020], можно предположить, что размеры очага и магнитуда события могли оказаться больше фактических, если бы развитие разрыва не было остановлено региональными тектоническими особенностями рассматриваемой территории: разломом Гарлок с одной стороны и зоной вулканов с другой. Таким образом, установление характерных пространственных размеров зон “раскачивания” может дать важнейшую информацию о потенциальной магнитуде готовящегося землетрясения.

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

Следует отметить, что используемые данные при получении результатов по процессу раскачивания, лежащие в основе создания модели “Качели”, были рассчитаны в рамках геомеханической модели, куда вводилась текущая сейсмичность. В процессе развития самой текущей сейсмичности такие особенности динамики НДС выделить весьма затруднительно.

Действительно, слабая сейсмичность, которая связана с локальным разрушением пород, вызывает перераспределение деформаций в сложно построенной земной коре с ее тектоническими особенностями – разломно-блоковой, достаточно мозаичной структурой. Именно в геомеханической модели, в соответствии с развитием сейсмического процесса и особенностями локальной тектоники, выделяются области, где возникают максимумы упругой энергии, касательных напряжений и сдвиговых деформаций. Эти области зачастую не совпадают с местами максимального проявления текущей сейсмичности, а обуславливаются так называемой “чувствительностью” среды, что демонстрируется на рис. 11.

Примеры такой “чувствительности” среды отмечаются в работах [Гохберг и др., 2004; Бондур и др., 2021], когда в модели воздействие в виде локальной поверхностной нагрузки отражается на концентрации процесса приближения пород к пределу прочности на разных глубинах, в достаточно удаленных областях согласно сложному тектоническому строению.

Из приведенного рисунка видно, что для одного достаточно характерного момента времени (01.08.2018 г., примерно за год до землетрясения), пространственное распределение зон с максимальными значениями приращений СД не совпадает с распределением текущей сейсмичности. Надо отметить, что на всем исследуемом промежутке времени 2018–2019 гг. кластер локальной слабой сейсмичности, в основном, приурочен к северо-западному крылу очага. При этом, эффект раскачивания, связанный с последовательной миграцией аномальных зон СД, в определяемой таким образом и подготовленной к разрушению перемычке выявлен при геомеханическом моделировании.

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

ЗАКЛЮЧЕНИЕ

Представленная в работе концепция позволяет глубже понять процессы, происходящие в земной коре перед сильными землетрясениями и связанные с динамикой напряженно-деформированного состояния, мониторинг которого может использоваться для краткосрочного прогнозирования. Важнейшей особенностью выделенного предвестника землетрясения M = 7.1 Риджкрест 2019 г. является т.н. “раскачивание” на флангах будущего разрыва, ассоциированного с землетрясением, выражающееся в циклическом появлении аномалий деформации сдвига на его северо-западном и юго-восточном продолжениях. Представляется, что подобные закономерности изменения напряженно-деформированного состояния перед землетрясением Риджкрест, выявленные на основе геомеханической модели сплошной среды с динамическим учетом текущей сейсмичности, могут быть объяснены в рамках относительно простых блочных моделей. Дальнейшее развитие такого подхода может включать оценку геомеханических параметров системы блоков на основе решения обратных задач по экспериментальным и симуляционным данным, а также последующее моделирование ее дальнейшей эволюции с целью краткосрочного прогноза сильных сейсмических событий.

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