Доклады Российской академии наук. Науки о Земле, 2021, T. 501, № 2, стр. 204-209

С использованием методики моделирования динамики блоковых структур и сейсмичности ([1–5] и др.) выполнено численное моделирование для региона Алтай–Саяны–Прибайкалье. Блоковая структура представляется в виде слоя толщины H между двумя параллельными горизонтальными плоскостями. Ее разделение на блоки построено на основе приведенной на рис. 1 схемы морфоструктурного районирования региона [6]. Структура состоит из 14 блоков, которые ограничены линеаментами только первого и второго рангов (рис. 2). Подробное описание методики моделирования динамики блоковых структур приведено в [5]. Одним из ее основных принципов является предположение о том, что блоки разделены узкими ослабленными и менее консолидированными зонами разломов, по сравнению с которыми на характерных временах сейсмотектонического процесса (порядка нескольких тысяч лет) блоки могут считаться абсолютно жесткими. Землетрясения в модели возникают в зонах разломов. Целью моделирования являлось получение путем подбора параметров модели пространственного распределения эпицентров модельных землетрясений, которое согласуется с положениями эпицентров землетрясений, произошедших в регионе. В этом случае модельная сейсмичность может дать дополнительные аргументы в пользу достоверности результата распознавания мест возможного возникновения сильных землетрясений в регионе [6], а каталоги землетрясений, полученных в модели, могут быть использованы для оценки сейсмической опасности [7].

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

где f – вектор плотности упругой силы в рассматриваемой точке, принадлежащей подошве блока или плоскости разлома, Δr – вектор смещения блока относительно подстилающей среды в нижней плоскости или вектор относительного смещения соседних блоков в плоскости разлома, d – вектор неупругого смещения. В случае плоскости разлома коэффициент K может быть интерпретирован как модуль сдвига зоны разлома, деленный на ее ширину. Коэффициент W определяет скорость роста неупругих смещений, и при его увеличении большая часть упругой энергии высвобождается через неупругие смещения. Наряду с упругой силой, лежащей в плоскости разлома, действует сила реакции, которая перпендикулярна плоскости разлома.

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

Для выполнения численного моделирования поверхности, на которых действуют силы (плоскости разломов и подошвы блоков), разбиваются на ячейки с линейными размерами, не превосходящими величину ε, значение которой задается. Предполагается, что вектор неупругого смещения d, координаты X, Y и векторы плотности сил одинаковы во всех точках одной ячейки. Состояние блоковой структуры рассматривается в дискретные моменты времени t_i = t₀ + iΔt (i = 1, 2, …), где t₀ – начальный момент, а Δt – постоянный шаг.

Момент возникновения землетрясения определяется значением отношения

где f и p – величины действующих на единицу площади плоскости разлома упругой силы, вектор которой лежит в плоскости разлома, и силы реакции, нормальной к плоскости разлома, а постоянная P имеет одно и то же значение для всех разломов и интерпретируется как разность между литостатическим и гидростатическим давлением. В случае растяжения в направлении, перпендикулярном плоскости разлома, в знаменателе (2) стоит “–”, в случае сжатия – “+”. Для каждого разлома задаются значения трех порогов B > H_f> H_s. Начальные значения смещений и поворотов внутренних блоков и неупругих смещений на разломах задаются таким образом, чтобы везде было выполнено соотношение κ < B. Если в процессе моделирования величина κ на каком-либо разломе (в одной или нескольких ячейках) достигает порога B соответствующего разлома, то возникает разрыв (“землетрясение”). Разрыв означает проскальзывание, в результате которого неупругие смещения в ячейке резко изменяются так, что величина κ уменьшается до значения порога H_f, т.е. землетрясение происходит согласно модели сухого трения, а коэффициент B играет роль порога прочности. После пересчета неупругих смещений определяются новые положения блоков из условия квазистатического равновесия. В случае, когда после этого значение κ превосходит порог B для какой-либо ячейки или ячеек, процесс пересчета неупругих смещений и определения новых положений блоков повторяется. Ячейки, в которых значение κ превосходило порог B, образуют землетрясение. В качестве координат эпицентра и глубины очага землетрясения вычисляются средние взвешенные значения координат и глубин таких ячеек с весами, пропорциональными их площадям. Магнитуда землетрясения M вычисляется по формуле где S – сумма площадей (в км²) ячеек, образующих землетрясение, а значения коэффициентов выбраны в соответствии с [8].

После землетрясения соответствующие ячейки разломов находятся в состоянии крипа. Это означает, что для них скорость роста неупругих смещений больше, чем в обычном состоянии при той же величине плотности упругой силы, т.е. в (1) вместо постоянной W используется постоянная W_s (W_s> W). Состояние крипа прекращается, когда наступает момент времени, в который κ < H_s.

Заданы следующие значения параметров блоковой структуры региона Алтай–Саяны–Прибайкалье: расстояние между верхней и нижней плоскостями, ограничивающими блоковую структуру, т.е. толщина всех блоков структуры, H = 30 км, что соответствует толщине земной коры и глубинам землетрясений, зарегистрированных в регионе; параметр P в выражении (2) равен 2 Кбар; дискретизация времени и пространства определяется значениями Δt = 0.001 и ε = 1 км; значения параметров в выражениях (1) и порогов для величины κ, определенной выражением (2), одни и те же для всех разломов: K = 1 бар/см, W = 0.005 см/бар, W_s = 0.5 см/бар, B = 0.1, H_f = 0.085, H_s = 0.07; значения параметров в выражениях (1) одни и те же для всех блоков: K = 1 бар/см, W = 0.01 см/бар. Отметим, что поскольку безразмерное время используется при моделировании, то перечисленные значения параметров соответствуют одной единице безразмерного времени.

Численное моделирование выполнялось с нулевыми начальными условиями (нулевые смещения границ структуры и подстилающей среды под блоками, а также нулевые неупругие смещения во всех ячейках) в течение 500 единиц безразмерного времени. Несколько сотен численных экспериментов выполнено с моделью. Установлено, что распределение эпицентров модельных землетрясений в наибольшей степени зависит от задаваемых движений границ структуры и подстилающей среды под блоками. Соответственно эти движения варьировались с целью получить эпицентры модельных землетрясений в местах, где известны землетрясения, произошедшие в регионе. Использованный каталог землетрясений региона получен на основе регионального каталога [9] и данных Национального центра информации о землетрясениях (NEIC) США за период 1900–2013 гг. С целью исключить влияние нестационарности начального периода моделирования искусственные каталоги рассматривались на отрезке в 400 единиц безразмерного времени со 100 до 500.

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

На рис. 3а показаны положения эпицентров землетрясений с M ≥ 6.0, полученных в численном эксперименте, в котором они присутствуют в окрестностях проекций на ближайшие разломы эпицентров всех реальных землетрясений с M ≥ 6.0 (рис. 3б). Аналогичное сравнение эпицентров землетрясений с M ≥ 5.0 показано на рис. 4. В этом случае проекции большинства эпицентров реальных землетрясений на ближайшие разломы имеют в окрестностях эпицентры модельных землетрясений. Максимальная магнитуда (7.88) землетрясений, полученных в модели, меньше, чем максимальная магнитуда (8.2) землетрясений, известных в регионе, но эпицентры сильнейших (с M ≥ 7.85) модельных землетрясений расположены вблизи эпицентра сильнейшего землетрясения региона (рис. 3).

Из рис. 3 и 4 следует, что распределение эпицентров модельных землетрясений достаточно хорошо воспроизводит положения эпицентров реальных землетрясений. В окрестностях эпицентров всех известных землетрясений с M ≥ 6.0 и практически всех с M ≥ 5.0 имеются эпицентры соответствующих модельных землетрясений.

Сравнение положений эпицентров модельных землетрясений с M ≥ 6.0 (рис. 3а) с полученными в [6] результатами распознавания мест возможного возникновения таких землетрясений (рис. 1) показывает их согласованность. Особо отметим, что модельные землетрясения с M ≥ 6.0 возникли на северной оконечности оз. Байкал, где такие землетрясения не известны, но по результатам распознавания [6] возможны. В целом результаты моделирования дают дополнительные аргументы в пользу достоверности результатов распознавания.

На рис. 5 представлены кумулятивные графики Гутенберга–Рихтера для модельной и реальной сейсмичности, наклоны которых близки. Из сравнения кривых, представленных на рис. 5, можно оценить продолжительность периода времени, покрытого искусственным каталогом землетрясений. Рассмотрим, например, события с M ≥ 6.0. Число таких землетрясений, произошедших в регионе, равно 15, а в модели получено 2045 таких событий. Если предположить, что интенсивности потоков реальных и модельных землетрясений близки, то получаем, что период времени, покрытый искусственным каталогом, в 136 раз длиннее продолжительности реального каталога, которая составляет порядка 110 лет. Следовательно, искусственный каталог соответствует примерно 14 000 годам, а одна единица безразмерного времени – 35 годам.