A Spherical Block Model of Lithosphere Dynamics and Seismicity: Current State and Development Prospects

Abstract

A description of the evolution of a spherical block model of the dynamics and seismicity of the lithosphere is given. The main focus is on the current version and the introduction of a constructive automatic calibration (parameter selection) procedure to obtain the best approximation of key properties of regional and/or global seismicity. The paper presents some results of computational experiments.

Full Text

ВВЕДЕНИЕ: АКТУАЛЬНОСТЬ МОДЕЛИРОВАНИЯ

Статистический анализ сейсмичности как пространственно-временной последовательности землетрясений на основе каталогов зарегистрированных событий крайне затруднен ввиду короткой истории надежных инструментальных наблюдений. Явления, обнаруженные в реальных каталогах землетрясений, могут быть единичными и не повторяться в будущем. Синтетические же каталоги, полученные путем численного моделирования, могут покрывать сколь угодно длительные интервалы времени, что позволяет анализировать значимость исследуемых свойств сейсмического потока, в частности, выявлять/подтверждать закономерности, предшествующие сильным толчкам, что может быть востребовано в экспертных системах мониторинга регионального и глобального сейсмического риска [Gabrielov, Newman, 1997; Keilis-Borok, Soloviev, 2003; Ismail-Zadeh et al., 2018; Ismail-Zadeh, Soloviev, 2022]. Основным результатом моделирования сейсмичности литосферы является синтетический каталог землетрясений, в котором каждое событие характеризуется моментом времени, координатами эпицентра, глубиной и магнитудой. Моделирование динамики земной коры предполагает получение поля скоростей движения точек на разных глубинах, действующих сил, обусловленных ими смещений, а также характера взаимодействия структурных элементов.

Общепринято среди различных подходов к моделированию литосферных процессов (см., например, работу [Gabrielov, Newman, 1997] и библиографию к ней) выделять два основных направления. Первое опирается на детальное исследование одного специфического тектонического разлома или, нередко, одного конкретного сильного землетрясения с целью воспроизведения определенных пре- и/или постсейсмических явлений, характерных для данного разлома или события. Модели второго направления трактуют сейсмотектонический процесс гораздо более абстрактно; основной задачей моделирования является получение универсальных свойств сейсмичности, обнаруженных эмпирическим путем (например, степенного закона распределения “размера” событий (закона повторяемости Гутенберга–Рихтера), кластеризации, миграции событий, сейсмического цикла и т.д.). Представляется, однако, что адекватная модель должна не только отражать некоторые общие свойства нелинейных динамических систем, но и учитывать геометрию взаимодействующих тектонических разломов. Блоковые модели динамики и сейсмичности литосферы [Keilis-Borok, Soloviev, 2003; Ismail-Zadeh et al., 2018; Ismail-Zadeh, Soloviev, 2022] разрабатывались с учетом обоих требований. В них рассматривается система абсолютно жестких блоков, находящаяся в состоянии квазистатического равновесия; при этом модельное событие представляет собой резкий сброс напряжений, возникающих на разломах, разделяющих блоки, под действием внешних сил. Два главных механизма, включенных в сейсмотектонический процесс, тектоническое нагружение с характеристической скоростью в несколько см/год и перераспределение упругого напряжения с характеристической скоростью в несколько км/с, трактуются в модели в стандартной временной шкале как, соответственно, равномерное движение и мгновенный сброс напряжения. Волновые процессы остаются вне рамок существующих блоковых моделей. Плоская модель, в которой структура ограничена двумя горизонтальными плоскостями, является наиболее изученной; на ее основе построены аппроксимации реальных сейсмических регионов [Keilis-Borok, Soloviev, 2003; Panza et al., 1997; Peresan et al., 2007; Ismail-Zadeh et al., 2007; Соловьев, Горшков, 2017; 2021; Vorobieva et al., 2017; 2019]. Однако при попытке моделирования динамики глобальных тектонических плит обнаружены существенные неточности, для преодоления которых введена сферическая геометрия [Мельникова и др., 2000; Rozenberg et al., 2005].

Настоящая работа, являясь продолжением исследований [Мельникова и др., 2000; Rozenberg et al., 2005; 2020; Мельникова, Розенберг, 2007; 2015; Digas et al., 2010; Melnikova et al., 2017; Розенберг, 2023], преследует две цели. Во-первых, она содержит ретроспективный обзор различных модификаций сферической блоковой модели и краткое описание текущей версии, а во-вторых, в ней обсуждаются возможности получения оптимального набора модельных параметров, обеспечивающего наилучшую аппроксимацию ключевых свойств региональной и/или глобальной сейсмичности, в результате внедрения усовершенствованной процедуры калибровки модели.

РАЗЛИЧНЫЕ МОДИФИКАЦИИ СФЕРИЧЕСКОЙ БЛОКОВОЙ МОДЕЛИ

Подход к моделированию опирается на представление тектонических плит в виде системы абсолютно жестких блоков на сфере. Блоковая структура является ограниченной и односвязной частью шарового слоя глубиной H, заключенного между двумя концентрическими сферами, одна из которых (внешняя) интерпретируется как поверхность Земли, другая (внутренняя) – как нижняя граница упругой литосферы. Разделение структуры на блоки определяется пересекающими этот слой бесконечно тонкими разломами, каждый из которых представляет собой коническую поверхность, наклоненную под определенным углом к внешней сфере. Общие точки двух разломов на внешней и внутренней сферах называются вершинами. Участки разломов, ограниченные соответствующими парами соседних вершин, называются сегментами. Пересечения блока с ограничивающими сферами представляют собой сферические многоугольники, при этом пересечение с нижней (для блока) сферой называется подошвой. Предполагается, что вне блоковой структуры могут находиться граничные блоки, примыкающие к внешним сегментам. Другая возможность состоит в рассмотрении блоковой структуры, замкнутой на сфере. Блоки считаются абсолютно жесткими, все их смещения – бесконечно малыми по сравнению с линейными размерами, поэтому геометрия блоковой структуры не меняется в процессе моделирования, и структура не движется как единое целое. Гравитационными силами можно пренебречь, так как они слабо зависят от смещений блоков и блоковая структура в начальный момент времени находится в состоянии квазистатического равновесия. Блоки (в том числе и граничные) имеют шесть степеней свободы. Смещение каждого блока состоит из поступательной и вращательной компонент. Предполагается, что законы движения граничных блоков и подстилающей среды известны, при этом движение описывается как вращение на сфере, т.е. задаются положение оси вращения и угловая скорость. Во все моменты времени система находится в состоянии квазистатического равновесия; при этом модельным событием является резкий сброс тектонических напряжений, возникающих на разломах, разделяющих блоки, под действием внешних сил. В случае замкнутой структуры единственным источником модельных трансформаций является движение подстилающей среды, которое определяется как вращение на сфере согласно модели HS3-NUVEL1 [Gripp, Gordon, 2002].

Подробное описание процесса усовершенствования сферической блоковой модели динамики и сейсмичности литосферы можно найти в работах [Мельникова, Розенберг, 2007; Мельникова, Розенберг, 2015]; в данной работе ограничимся краткой характеристикой основных этапов эволюции модели.

Разработано несколько модификаций модели, зависящих от способа трактовки глубины сферического слоя. В первой модификации модели (без глубины) литосфера рассматривалась как тонкий слой, все характеристики точек структуры определялись только их координатами и не зависели от глубины, поскольку толщина литосферы значительно меньше линейных размеров блоков [Мельникова и др., 2000; Мельникова, Розенберг, 2007]. Главное преимущество модификации состоит в значительной экономии времени счета при моделировании, что может быть существенно при большом количестве запусков в эксперименте по вариации того или иного параметра; основной недостаток – невозможность учета углов наклона разломов, фактически определяющих характер сейсмичности, и других свойств неоднородной литосферы. Появившаяся позже модификация с постоянной глубиной [Rozenberg et al., 2005] использовала предположение об однородности литосферы по глубине (все блоки имели одну и ту же глубину H, а свойства всех частей блока (разлома) были одинаковыми). Модификация, используемая в настоящее время, предусматривает возможность задания различных глубин (в пределах H) для разных блоков и учета зависимости вязко-упругих свойств разлома от его глубины [Мельникова, Розенберг, 2007; 2015]. Отметим, что, по существу, эта опция явилась первой попыткой учета неоднородности литосферы (например, различий в строении континентальной и океанической коры и уменьшения вязкости коры с глубиной) в блоковых моделях.

СОВРЕМЕННАЯ ВЕРСИЯ СФЕРИЧЕСКОЙ БЛОКОВОЙ МОДЕЛИ

Приведем краткое описание текущей версии сферической блоковой модели, разработанной с возможностью учета случайных факторов. Основные ее принципы изложены в работе [Мельникова, Розенберг, 2015], версия окончательно сформирована в работе [Rozenberg, 2020].

Поскольку блоки являются абсолютно жесткими, все деформации имеют место на разломах и подошвах блоков, где силы возникают из-за смещений соседних элементов структуры. Для подсчета различных криволинейных интегралов выполняется дискретизация (разбиение на ячейки) сферической поверхности подошв блоков и сегментов разломов, при этом предполагается, что значения сил и смещений совпадают для всех точек ячейки. Приведем формулы для определения упругой силы $(f_t,f_l,f_n),$ действующей на единицу площади разлома:

$\begin{array}{c}{f}_t\left(\tau \right)=K_t\left({\Delta }_t\right(\tau )-{\delta }_t(\tau \left)\right),\\ f_l\left(\tau \right)=K_l\left({\Delta }_l\right(\tau )-{\delta }_l(\tau \left)\right),\\ f_n\left(\tau \right)=K_n\left({\Delta }_n\right(\tau )-{\delta }_n(\tau \left)\right).\end{array}$ (1)

Здесь $\tau$ – время; $(t,l,n)$ – система координат, связанная с точкой приложения силы (оси t, l лежат в плоскости, касательной к поверхности разлома, ось n ей перпендикулярна); ${\Delta }_t,$ ${\Delta }_l,$ ${\Delta }_n$ – компоненты относительного смещения в системе $(t,l,n)$ соседних элементов структуры; ${\delta }_t,$ ${\delta }_l,$ ${\delta }_n$ – соответствующие неупругие смещения, зависимость от времени которых описывается квазилинейными стохастическими дифференциальными уравнениями вида:

$d{\delta }_t\left(\tau \right)=W_t{K}_t\left({\Delta }_t\right(\tau )-{\delta }_t(\tau \left)\right)d\tau +{\lambda }_t{\delta }_t\left(\tau \right)d{\xi }_t\left(\tau \right),$

$d{\delta }_l\left(\tau \right)=W_l{K}_l\left({\Delta }_l\right(\tau )-{\delta }_l(\tau \left)\right)d\tau +{\lambda }_l{\delta }_l\left(\tau \right)d{\xi }_l\left(\tau \right),$

$d{\delta }_n\left(\tau \right)=W_n{K}_n\left({\Delta }_n\right(\tau )-{\delta }_n(\tau \left)\right)d\tau +{\lambda }_n{\delta }_n\left(\tau \right)d{\xi }_n\left(\tau \right),$ (2)

записанными для всех ячеек сегментов. Полагаем, что стандартные независимые винеровские процессы ${\xi }_t,$ ${\xi }_l$ и ${\xi }_n$ (т.е. выходящие из нуля, имеющие нулевое математическое ожидание и дисперсию, равную $\tau )$ являются скалярными, а коэффициенты ${\lambda }_t,$ ${\lambda }_l,$ $r_s$ регулируют амплитуду случайных воздействий на всем сегменте, обеспечивая ее малость по сравнению с величиной неупругих смещений на всем промежутке моделирования. Единственным решением каждого из уравнений (2), понимаемых в смысле Ито, является нормальный марковский случайный процесс с непрерывными реализациями [Oksendal, 2003]. Коэффициенты $K_t,$ $K_l,$ $K_n$ (1)–(2), характеризующие упругие свойства разлома, и коэффициенты $W_t,$ $W_l,$ $W_n$ (2), характеризующие вязкие свойства разлома, могут быть различными для разных разломов и, кроме того, могут изменяться в зависимости от глубины. В начальный момент времени величины, вычисляемые согласно (1) и (2), имеют нулевые значения. Аналогично выглядят формулы для вычисления сил и неупругих смещений на подошвах блоков.

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

$Aw=b.$ (3)

Компонентами неизвестного вектора $w=(w_1,w_2\mathrm{,...,}{w}_{6n})$ являются смещения и углы поворота внутренних блоков (n – число таких блоков). Элементы матрицы A (размерности $6n\times 6n)$ не зависят от времени и могут быть вычислены один раз в начале процесса. При реалистичных значениях параметров модели матрица A является невырожденной, и система (3) имеет единственное решение.

В каждый момент дискретного времени ${\tau }_i$ при вычислении компонент силы, действующей на ячейку разлома, определяется безразмерная величина $\kappa ,$ трактуемая как модельное напряжение:

$\kappa =\frac{\sqrt{f_t^2+f_l^2}}{P-f_n}.$ (4)

Здесь P – одинаковый для всех разломов параметр, который может интерпретироваться как разность между литостатическим и гидростатическим давлением. Таким образом, фактически величина $\kappa$ является отношением модуля силы, стремящейся сдвинуть блоки вдоль разлома, к модулю силы, прижимающей блоки друг к другу. Для каждого разлома задаются значения трех порогов прочности, вообще говоря, зависящие от времени:

$B>H_f\ge H_s,$

$\begin{array}{c}B=B\left({\tau }_i\right)=B_0+\sigma X\left({\tau }_i\right),\\ H_f=H_f\left({\tau }_i\right)=aB\left({\tau }_i\right),\\ H_s=H_s\left({\tau }_i\right)=bB\left({\tau }_i\right).\end{array}$ (5)

Полагаем, что $0<B_0<\mathrm{1,}$ $0<\sigma <<B_0,$ для всех i $X\left({\tau }_i\right)$ – случайная величина, распределенная по стандартному нормальному закону $N(0;1),$ $0<a<\mathrm{1,}$ $0<b\le a,$ и начальные условия таковы, что неравенство $\kappa <B$ имеет место во всех ячейках структуры.

Взаимодействие между блоками (между блоком и соседней подстилающей средой) полагается вязкоупругим (нормальное состояние) до тех пор, пока величина $\kappa$ (4) на части разлома, разделяющего элементы структуры, не достигает заданного порога. Такая ситуация интерпретируется как землетрясение. В ячейках, попавших в “критическое” состояние, в соответствии с законом сухого трения, происходит резкий сброс напряжения посредством изменения значений неупругих смещений ${\delta }_t,{\delta }_l,{\delta }_n$ таким образом, чтобы после пересчета сил по формулам (1)–(2) величина модельного напряжения (4) была равна $H_f.$ Затем находится правая часть системы (3), определяются векторы сдвига и углы поворота блоков. Если вновь в какой-либо ячейке $\kappa \ge B,$ то вся процедура повторяется. Когда во всех ячейках на разломах $\kappa <B,$ вычисления продолжаются по обычной схеме. Считается, что ячейки, в которых произошли землетрясения, находятся в состоянии крипа. Это означает, что для них в уравнениях (2) для вычисления значений неупругих смещений используются параметры $W_t^s$ $(W_t^s>W_t),$ $W_l^s$ $(W_l^s>W_l)$ и $W_n^s$ $(W_n^s>W_n),$ обеспечивающие значительно более быстрый, по сравнению с нормальным состоянием, рост неупругих смещений и, следовательно, уменьшение значений сил и напряжений. Состояние крипа продолжается до тех пор, пока $\kappa >H_s$ (это процесс “заживления”), после чего ячейки возвращаются в нормальное состояние с использованием при расчетах $W_t,$ $W_l$ и $W_n$ и снова могут накапливать напряжение. Отметим, что для отражения неподдающихся точному аналитическому описанию зависимостей в модели реализовано две схемы введения стохастической компоненты: во-первых, “зашумление” дифференциальных уравнений (1)–(2), описывающих динамику сил и смещений, и во-вторых, использование случайных величин при задании порогов прочности среды тектонических разломов (5). Вычислительные эксперименты показали, что дополнение процедуры определения модельного события элементами случайности позволило заметно улучшить свойства синтетической сейсмичности [Мельникова, Розенберг, 2015; Melnikova et al., 2017; Rozenberg, 2020].

Основным результатом процесса моделирования является синтетический каталог землетрясений. Принадлежащие одному разлому ячейки, в которых произошло землетрясение в момент времени ${\tau }_i,$ объединяются в одно событие. Географические координаты его эпицентра и глубина вычисляются как взвешенные суммы координат и глубин ячеек (вес ячейки определяется как отношение ее площади к сумме площадей всех ячеек, вовлеченных в землетрясение). Взвешенная сумма добавок к неупругим смещениям ${\delta }_t$ и ${\delta }_l$ при сбросе напряжения аппроксимирует случившуюся подвижку блоков вдоль разлома и позволяет определить механизм модельного события, который информирует о процессе распространения различных сейсмических волн от очага. В зависимости от направления подвижки и угла наклона разлома принято выделять следующие основные механизмы: сдвиг, сброс и взброс [Аки, Ричардс, 1983]. Магнитуда землетрясения вычисляется в зависимости от его механизма с использованием известных в сейсмологии эмпирических формул [Wells, Coppersmith, 1994]

$M=DlgS+E,$ (6)

где: S – сумма площадей ячеек (в км²); $D=\mathrm{1.02,}$ $E=3.98$ для сдвига; $D=\mathrm{1.02,}$ $E=3.93$ для сброса; $D=\mathrm{0.90,}$ $E=4.33$ для взброса.

Дополнительно модель позволяет получить картину мгновенной кинематики блоков и информацию о характере их взаимодействия на разломах и границе литосфера/мантия.

РЕАЛИЗАЦИЯ ПРОЦЕДУРЫ КАЛИБРОВКИ МОДЕЛИ

Компьютерная реализация сферической блоковой модели потребовала значительных затрат памяти и времени работы процессора, что привело к необходимости применения параллельных вычислительных технологий. Современная версия, с возможностью использования в расчетах реальных геофизических и сейсмических данных, реализована в виде пакета программ, ориентированного на многопроцессорную технику в вычислительных процедурах и на персональную в сервисных (см. подробное описание в работах [Digas et al., 2010; Melnikova et al., 2017]). Отметим программу визуализации состояния разлома [Melnikova et al., 2017], которая, обеспечивая возможность анализа распределения модельного напряжения и миграции событий вдоль структуры, служит эффективным инструментом как для калибровки модели, так и для верификации результатов моделирования.

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

$K*=argmin\left\{\sum _{i=1}^{N^{reg}}\text{\hspace{0.17em}}{\alpha }_i^{reg}\text{}\sum _{j=1}^{N_i^{crit}}{\beta }_{ij}^{crit}{\rho }_j(H_{ij}^{real}\text{},H_{ij}^{model}\text{}(K_{ij}\left)\right),K_{ij}\in S_K\right\}.$ (7)

Здесь: $K*$ – оптимальный набор модельных параметров; $S_K$ – некоторое конечное множество допустимых значений модельных параметров $K_{ij};$ $N^{reg}$ – количество рассматриваемых регионов; ${\alpha }_i^{reg}$ – весовые коэффициенты для регионов; ${\sum }_i{\alpha }_i^{reg}=\mathrm{1,}$ $N_i^{crit}$ – количество региональных критериев; ${\beta }_{ij}^{crit}$ – весовые коэффициенты для региональных критериев, ${\sum }_j{\beta }_{ij}^{crit}=\mathrm{1,}$ $H_{ij}^{real}$ и $H_{ij}^{model}\left(K_{ij}\right)$ – соответственно, реальные и модельные характеристики региона i, составляющие критерий j; ${\rho }_j$ – нормированное, т.е. деленное на максимально возможное для этого критерия, отклонение j-й характеристики в подходящей метрике (очевидно, со значениями от 0 до 1). Отметим, что реализация возможности индивидуального для региона набора критериев объясняется различиями в доступности, полноте и надежности региональных данных.

Для получения множества допустимых значений параметров $S_K$ (каждый его элемент – конечномерный вектор, определяющий набор параметров для расчетов и, соответственно, тестируемый вариант) мы используем вариацию коэффициентов, описывающих вязко-упругое взаимодействие по формулам типа (1), (2), вариацию различных случайных факторов из соотношений (4), (5) и некоторые поправки к движениям подстилающей среды в определенных “разумных” диапазонах изменения с некоторыми “разумными” шагами дискретизации. В качестве региональных/глобальных критериев, характеризующих качество модельных данных, используем сравнительный анализ ключевых характеристик динамики и сейсмичности систем тектонических плит; некоторые их них обсуждались в работе [Rozenberg, 2020].

Перечислим основные критерии.

1. Пространственное распределение эпицентров сильных событий. Для сравнения реального и модельного распределений вводится ограничение снизу по магнитуде для сильных событий, реальное распределение ограничивается окрестностями разломов (поскольку внутри блоков модельных событий нет по определению). После этого вычисляется расстояние по Хаусдорфу (в естественной метрике на поверхности Земли) между двумя множествами эпицентров. Отметим, что два множества близки в метрике Хаусдорфа, если каждая точка первого множества близка к некоторой точке второго и наоборот. Расстоянием по Хаусдорфу называется наибольшее из расстояний от точки одного множества до ближайшей точки другого. Таким образом, для двух непустых множеств эпицентров, реального и модельного $E_m,$ расстояние по Хаусдорфу определяется так:

$\begin{array}{c}{d}_H(E_r,E_m)=\\ =max\left\{\underset{e_r\in E_r}{sup}\underset{e_m\in E_m}{inf}\text{}d\right(e_r,e_m),\underset{e_m\in E_m}{sup}\underset{e_r\in E_r}{inf}\text{}d(e_r,e_m\left)\right\},\end{array}$ (8)

где $d(e_r,e_m)$ – расстояние между точками $e_r$ и $e_m$ на поверхности Земли.

2. Распределение землетрясений по глубине. Для сравнения реального и модельного распределений в этом аспекте мы не проводим тонкий количественный анализ, ограничиваясь определением класса события (например, поверхностное, промежуточное или глубокое). Поэтому, после введения характеристических глубин и магнитудных порогов, мы вычисляем долю событий в каждом интервале по глубине относительно общего количества событий. Тогда отклонение между реальными $Dr$ и модельными $D_m$ данными можно определить так:

$d_D(D_r,D_m)=\sum _{i=1}^n|D_r^i-D_m^i|,$ (9)

где n – выбранное количество интервалов по глубине; $D_r^i$ и $D_m^i$ – доли событий, принадлежащих соответствующему интервалу, для реального и модельного каталогов.

3. Закон повторяемости Гутенберга–Рихтера. Соотношение характеризует степенное распределение землетрясений по магнитуде, при этом внимание уделяется линейности и углу наклона графика зависимости аккумулированного количества событий N от магнитуды M. Все такие зависимости аппроксимируются линейной регрессией $lgN=c-SM,$ построенной по методу наименьших квадратов. Величина S служит оценкой угла наклона, а среднее расстояние между точками графика и построенной прямой трактуется как ошибка аппроксимации A. Для сравнения параметров, реальных $G_r=(S_r,A_r)$ и модельных $G_m=(S_m,A_m),$ используется следующее соотношение:

$d_G(G_r,G_m)={\alpha }_1|S_r-S_m|+{\alpha }_2|A_r-A_m|,$ (10)

где $S_r,$ $S_m$ и $A_r,$ $A_m$ – соответственно, реальные и модельные углы наклона линий регрессии и ошибки аппроксимации, ${\alpha }_1,{\alpha }_2>\mathrm{0,}$ ${\alpha }_1+{\alpha }_2=1$ – весовые коэффициенты.

4. Смещения на границах плит. Фактически речь идет о количественной оценке качества аппроксимации характера межплитовых границ (зоны растяжения, сжатия, сдвига). На знаковых границах плит, в локальных системах координат, описанных в формулах (1)–(2), рассматриваются, как правило, на поверхности Земли, трехмерные векторы относительного смещения $\Delta =({\Delta }_t,{\Delta }_l,{\Delta }_n)$ соседних блоков. Суммарное отклонение модельных величин от реальных (посчитанных в соответствии с некоторой общепринятой моделью кинематики плит, например, HS3-NUVEL1) вычисляется в естественной метрике:

$d_{\Delta }({\Delta }_{},{\Delta }_{})={||{\Delta }_{}-{\Delta }_{}||}_{R^3},$ (11)

где ${\Delta }_r$ и ${\Delta }_m$ – соответственно, реальный и модельный векторы смещений.

В идеальном варианте все характеристики типа (8)–(11) для всех регионов следует использовать, после нормирования, в качестве отклонений ${\rho }_j$ для решения задачи (7), но в действительности надежных данных наблюдений не хватает (степень достоверности можно регулировать значениями весовых коэффициентов ${\alpha }_i^{reg}$ и ${\beta }_{ij}^{crit}.$ Далее, в соответствии с (7), следует выбрать вариант с наименьшей агрегированной ошибкой. Поскольку ${\rho }_j$ – нормированные характеристики со значениями от 0 до 1, то наименьшее возможное значение этой ошибки равно 0 (для реальных данных), а наибольшее – 1 (для варианта, наихудшего по всем критериям). Очевидно, существуют и другие пути конструирования агрегированного критерия качества аппроксимации и процедуры ранжирования вариантов. Дополнительно отметим, что отсутствие анализа временных характеристик синтетической сейсмичности объясняется затруднениями в сравнении их с реальными данными ввиду невозможности надежно определить соответствие начального момента моделирования реальной дате во всех модификациях блоковых моделей.

НЕКОТОРЫЕ РЕЗУЛЬТАТЫ ЧИСЛЕННОГО ЭКСПЕРИМЕНТА

Вычислительные эксперименты по тестированию процедуры калибровки модели (7)–(11) требуют большого количества запусков на достаточно мелкой сетке по множеству параметров. Единственным способом точного решения задачи (7) является полный перебор допустимых значений входных параметров, что объясняется отсутствием явных аналитических зависимостей между входом и выходом блоковых моделей. Естественные ограничения на объем вычислений приводят к необходимости оптимизации перебора, трудновыполнимой без участия эксперта. Отметим, что в процессе калибровки одновременно применяются две схемы распараллеливания: во-первых, одновременный запуск вариантов на различных вычислительных узлах и, во-вторых, внутренняя параллелизация модели на каждом узле.

На самом первом этапе приложения сферической блоковой модели было проведено моделирование динамики и сейсмичности небольшой подсистемы плит, в которую входили Южноамериканская, Карибская, Кокос и Наска, при этом другие плиты, находящиеся рядом, представляли собой граничные блоки, законы движения которых известны. Выбор данного региона объясняется тем, что он включает в себя различные типы границ плит с достаточно контрастными движениями и высокой сейсмической активностью. Полученные результаты в целом соответствовали реальным данным [Мельникова и др., 2000]. Значительная часть последующих исследований связана с исследованием динамики и сейсмичности глобальной системы тектонических плит, покрывающей всю поверхность Земли. Геометрия рассматриваемой структуры (и фактически список регионов для функционала из задачи (7)) полагается неизменной: 15 плит (Наска, Южноамериканская, Кокос, Карибская, Североамериканская, Тихоокеанская, Африканская, Антарктическая, Евразийская, Аравийская, Индийская, Сомалийская, Филиппинская, Австралийская и Хуан де Фука), 186 вершин и 199 разломов. Граничные блоки отсутствуют, движение подстилающей среды определяется как вращение на сфере согласно модели HS3-NUVEL1 [Gripp, Gordon, 2002].

Кратко представим результаты численного эксперимента с глобальной системой тектонических плит, в котором варьировались два множества параметров, характеризующих взаимодействие блоков вдоль разломов. Для оптимизации перебора использовались общепринятые в блоковых моделях соображения, обусловленные наблюденной сейсмичностью [Keilis-Borok, Soloviev, 2003; Ismail-Zadeh et al., 2018; Ismail-Zadeh, Soloviev, 2022], именно, коэффициенты $K_t,$ $K_l$ и $K_n$ увеличиваются, а коэффициенты $W_t,$ $W_l$ и $W_n$ уменьшаются для разломов с высоким уровнем сейсмической активности и наоборот. Что касается случайных факторов, полагаем, что критическое значение  (4) должно быть приближенно равно 0.1 и амплитуда шума в квазилинейном стохастическом уравнении должна быть по крайней мере на порядок меньше величин смещений [Rozenberg, 2020].

Приведем “разумные” (из прошлого опыта вычислений) диапазоны и шаги изменения упомянутых параметров разломов: для коэффициентов, $K_t,K_l,K_n\in \left[\mathrm{1,}10\right],$ ${\Delta }_K=\mathrm{1,}$ $W_t,W_l,W_n\in \left[\mathrm{0.01,0.1}\right]$ ${\Delta }_W=\mathrm{0.01,}$ и $W_t^s,W_l^s,W_n^s\in \left[\mathrm{1,5}\right],$ ${\Delta }_s=1;$ для случайных факторов, ${\lambda }_t,{\lambda }_l,{\lambda }_n\in \left[\mathrm{0,}0.01\right],$ ${\Delta }_{\lambda }=0.005$ и $B_0\in \left[\mathrm{0.1,0.3}\right],$ ${\Delta }_B=\mathrm{0.1,}$ $\sigma \in \left[\mathrm{0,0.1}\right],$ ${\Delta }_{\sigma }=\mathrm{0.05.}$ Таким образом, формально мы имеем сетку размерности $10\times 10\times 5\times 3\times 3\times 3=13500.$ Это означает, что для точного решения задачи (7) на указанном множестве параметров следует перебрать 13 500¹⁹⁹ вариантов (199 – количество разломов), что, конечно, совершенно невозможно. На текущем этапе значительное уменьшение количества тестируемых вариантов в такой ситуации (назовем это сужением множества $S_K)$ достигается, главным образом, вручную посредством объединения разломов в группы с одинаковыми значениями параметров и использования соображений, приведенных выше.

Вычисления выполнялись в Институте математики и механики УрО РАН им. Н.Н. Красовского на гибридном вычислителе кластерного типа “Уран”, который имеет более 14 Тб оперативной памяти и пиковую производительность порядка 250 Тфлопс. Одновременно запускались 10 вариантов, 50 ядер на вариант для внутренней параллелизации модели. Высокая эффективность и хорошая масштабируемость задачи [Digas et al., 2010; Melnikova et al., 2017] позволили просчитывать вариант за приемлемое время: 100 “модельных” лет примерно за 80 минут (вместо 64 часов на одном ядре). Таким образом, 1000 вариантов можно просчитать за 135 часов.

В результате расчетов был выбран оптимальный вариант, решение задачи (7) на сужении множества $S_K$ его основные параметры: $K_t=K_l=K_n$ от 3 до 7, $W_t=W_l=W_n$ от 0.01 до 0.03, $W_t^s=W_l^s=W_n^s$ от 3 до 5, для различных групп разломов, ${\lambda }_t={\lambda }_l={\lambda }_n=\mathrm{0.01,}$ $B_0=\mathrm{0.1,}$ $\sigma =0.$ В качестве примера приведем результаты моделирования пространственного распределения сильных событий (рис. 1, рис. 2; при их подготовке использовалась программа Seismic Eruption). Отметим, что в синтетическом каталоге очевидным образом можно идентифицировать ряд черт, присущих реальной сейсмичности, а именно, два основных сейсмических пояса, Тихоокеанский и Средиземноморско-Трансазиатский, где происходит большая часть сильных событий; протяженную, но менее выраженную, сейсмичность срединно-океанических хребтов; увеличение сейсмической активности вблизи точек, где сходятся три и более плит. Как следствие, имеет место соответствие многих сейсмически активных и “спокойных” регионов реальным. Рассмотрено по 20 сильнейших по магнитуде событий реального (оказалось, что $M\ge 8.6)$ и модельного $(8.7\le M\le 8.9)$ каталогов. Большинство таких событий в модели произошли примерно в тех же местах, что и реальные (например, чилийские землетрясения 1960 г. и 2010 г. и Суматра-Андаманское событие 2004 г. хорошо отражены в модельной сейсмичности). В то же время, даже с учетом отсутствия модельных событий внутри блоков (а это специфика всех блоковых моделей), различия также довольно значительны.

Рис. 1. Зарегистрированная сейсмичность (NEIC, 01.01.1900–31.12.2022 [NEIC…]): эпицентры сильных землетрясений с магнитудой $M\ge 6.0$ показаны кругами, цвет которых зависит от глубины (шкала в правом нижнем углу), 20 сильнейших по магнитуде событий каталога отмечены звездочками.

Рис. 2. Синтетическая сейсмичность (100 единиц модельного времени): эпицентры сильных землетрясений с магнитудой $M\ge 6.0$ (см. (6)) показаны кругами, цвет которых зависит от глубины, 20 сильнейших по магнитуде событий каталога отмечены звездочками.

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

ЗАКЛЮЧЕНИЕ

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