Физика Земли, 2023, № 1, стр. 24-38

ВВЕДЕНИЕ

При проектировании социально значимых объектов проводится тщательное изучение сейсмической опасности региона. Местности, находящиеся в различных геодинамических состояниях по-разному реагируют на действие землетрясений. Известно, что при взбросах амплитуда колебаний в ближайшей зоне землетрясений примерно в два раза больше чем при сдвигах. Таким же образом по-разному изменяется статическое напряженное состояние земной коры после землетрясения по направлению простирания плоскости разрыва и по перпендикулярному направлению. Для территории Центральной Азии были установлены выборочные эмпирические оценки величин тектонических напряжений [Гзовский, 1964; Уломов, 1974], которые не потеряли свою актуальность. Существует достаточное число публикаций, где проводится оценка напряженного состояния по механизмам землетрясений и по геологическим данным [Осокина и др., 1979а; Шерман, Днепровский, 1989; Ребецкий и др., 2017]. Обширный обзор напряжений по всему миру приводится в работах [Zoback, 1992; Zoback, Zoback, 2007; Heidbach, 2018]. Такие данные имеются также по отдельным территориям Центральной Азии [Юнга и др., 2010; Rebetsky et al., 2012; Лукк, Юнга, 1988; Осокина и др.,1979б; Сычева, Мансуров, 2016]. Однако все они ограничиваются локальной реконструкцией геодинамического типа деформирования (растяжение–сдвиг–сжатие) и не дают количественную характеристику напряженного состояния земной коры Центральной Азии в целом. Более того, подобные работы зачастую не учитывают основное положение механики - уравнения равновесия и граничных условий, что справедливо критикуется в статье [Мухамедиев, 2019]. Основная трудность решения задачи о напряжениях на основе уравнений сплошных сред в геодинамике состоит в отсутствии достоверных граничных условий. В данной работе напряженное состояние земной коры Центральной Азии реконструируется решением обратной задачи на базе известных эмпирических данных. Для изучения изменения напряжений после землетрясений предлагается модель механизма очага землетрясений на основе моментной теории.

ПОСТАНОВКА ЗАДАЧИ

Тектонические землетрясения обусловлены разрушением сплошности земной коры в виде смещения под действием накопившегося дисбаланса объемных сил или возобновления бывшего разрыва, сцепленного благодаря цементированию поверхностей. Установлено, что наблюдаемые смещения на земной поверхности сильных, но не глубокофокусных землетрясений удачно объясняются с помощью модели диполя. Разновидности моделирования процесса разрыва зависят от постановки вопроса в рамках кинематической или динамической модели. В случае кинематической модели вектор сдвига не связан с порождаемой его силой. Разрыв полностью описывается векторной функцией сдвига бортов Δu(х_i) = u₁ – u₂ плоскости сдвига Σ, положение которой определяется азимутом 0° ≤ φ < 360°, наклоном 0° ≤ δ ≤ 90° и углом скольжения –180° ≤ λ ≤ 180° (рис. 1). Зона разрыва в простейших случаях считается в виде прямоугольника с площадью S = LD, где L – длина; D – ширина прямоугольника или круга с радиусом R; S = πR².

Моделированием сил, приводящих к разрыву земной коры при землетрясениях занимались многие исследователи (обширная библиография об их усилиях приведена, например, в работах [Nikano, 1923; Гоцадзе и др., 1957; Касахара, 1985]). Твердо установлено, что действия ординарного диполя с моментом (тип 1) или пары противоположно направленных диполей без момента (тип 2) в рамках упругости удовлетворительно объясняют квадрантное расположение областей растяжения и сжатия на поверхности Земли, замеченное еще профессором Oshi Shida из университета г. Киото в 1910 г. С помощью этих моделей получено хорошее совпадение между высвобождаемой сейсмической энергией и изменением упругой энергии, а также статическим и динамическим моментами в очаге землетрясения. Результаты расчета перемещений на поверхности сферы, порожденные дислокацией по обеим моделям дают одинаковую векторную картину и трудно различать разницу между ними. Проводились два симпозиума [Hodgson, 1959; 1961], чтобы выбрать, который из двух типов диполей дает более реальную картину сейсмического излучения очага, и специалисты не смогли придти к консенсусу. Наконец, подтверждение о близости угла поляризации поперечной волны по записям сейсмических волн с одноименным параметром модели, решило выбор в пользу 2-го типа.

МЕТОДЫ РЕШЕНИЯ

Задача о напряженном состоянии бесконечного упругого пространства под действием диполя, находящегося в конкретной точке в рамках безмоментной теории упругости, поставила вопрос об отыскания объемных сил, выражающие действия момента. Решение этого вопроса осуществлено впервые [Maruyama,1963; Ландау, Лифшиц, 1965] на основе уравнения равновесия сплошной среды для бесконечного пространства.

Другой подход получения эквивалентных сил основан на представлении перемещения бесконечно упругого тела через функцию Грина [Введенская, 1969; Костров, 1970; Райс, 1980; Аки, Ричардс, 1983]. Возникновение очага трактуется как результат взаимодействия разрыва и окружающей его деформации среды. При этом область землетрясения выделяется как предел сингулярной внутренней области, на границе Σ которой имеются разрывы перемещений, равные дислокациям. На границе Σ смещения терпят разрыв непрерывности, а напряжения остаются непрерывными в силу условия равновесия. Чтобы удерживать Σ неизменной, к ее точкам должны быть приложены эквивалентные массовые силы.

Мы предлагаем способ определения эквивалентных сил, опираясь на моментную теорию. Уравнение равновесия в моментной теории имеет следующий вид:

Здесь: σ_ij – компоненты тензора напряжений; индекс j после запятой означает дифференцирование по декартовым координатам x_j; F_i – массовая сила; µ_ij – компоненты моментного тензора напряжений; ε_ijk – тензор Леви–Чивита; M_i – составляющие удельного момента, имеющие размерность момента. По повторяющимся индексам проводится суммирование.

При учете механизма очага землетрясения полное напряжение σ_ij состоит из упругих τ_ij- и неупругих m_ij-составляющих, т.е.

В случае сосредоточенного двойного диполя M_i = 0. Уравнение (2) при расположении плоскости разрыва, например, показанный на рис. 1 (δ = 90°, φ = 0°, λ = 0°) имеет следующий вид:

Для упругой составляющей напряжения примем, как прежде, τ₁₂= τ₂₁, а моментный тензор выразим по его определению для сосредоточенного момента M:

Приравняв выражения, относящиеся к одноименным площадкам, получим:

Подставляя (7), (8) в (1) и переведя производные m_ij в правую часть, получаем эквивалентные силы, которые дают равные по величине, но противоположные по направлению моменты относительно оси перпендикулярной к плоскости:

Таким образом при учете механизма очага землетрясения к массовым силам F_i в системе (1) добавляются f_i.

Для расчета напряженного состояния земной коры для конкретного региона необходимо решать уравнения (1) в ограниченной области, привлекая соотношения между тензором напряжений и тензором деформации, выраженные через перемещения. Притом решение уравнений в частных производных (1) всецело зависит от граничных условий. Установления этих граничных условий являются камнем преткновения решения подобных задач. Нахождение напряжений по (1) даже без механизма очага землетрясения, само по себе является трудной задачей. Современный рельеф Центральной Азии сформировался в течении орогенического периода, который длится более 26 млн лет. Земная кора ведет себя как упругое тело по отношению к кратковременно действующим силам, но в течение геологического времени она деформируется пластически и проявляет свойства ползучести. Анализируя ориентацию и величину напряжений в различных поясах мира, авторы работы [Zoback, Zoback, 2007] пришли к выводу, что поле напряжений в литосфере является результатом современных сил, а не остаточных напряжений от прошлой тектонической активности. Поэтому в работе [Atabekov, 2020], напряженное состояние литосферы Центральной Азии реконструировано, исходя из сегодняшней картины. Граничные напряжения построены решением обратной задачи упругости на основе известной гипотезы тектоники плит. Согласно этой гипотезе, геодинамическая обстановка земной коры рассматриваемой территории определяется действием с одной стороны Евразийской, и с другой Индийской и Арабской плит. При этом решение задачи свелось к определению напряженного состояния призматического тела с толщиной литосферы, верхняя поверхность которого соответствует рельефу Центральной Азии (рис. 2). Это тело подвергается латеральным сжатиям. Верификация задачи проводилась сравнением решения с известными эмпирическими значениями напряжений.

Степень горизонтальной неоднородности литосферы учитывалась выбором различных упругих параметров модели. С этой целью использовалось проведенное ранее другими авторами деление литосферы Центральной Азии на 19 блоков, границами которых служат тектонические разломы (Приложение 1).

Решение полных систем трехмерных уравнений (1) даже численными методами довольно проблематично. С другой стороны, наиболее нерешенная проблема в геодинамике — это величина горизонтальных сил, ответственных за напряженное состояние литосферы, перемещения литосферных плит и взаимодействие блоков разного ранга. Согласно уравнению, обобщившему данные, имеющиеся в литературе о напряжениях в складчатом и кристаллическом фундаментах [Булин, 1971], наибольшие горизонтальные силы в глубинах 15 км имеют значения 461.9 МПа. Вертикальная компонента напряжений, взятая как вес вышележащей толщи на глубинах h в км и выраженная по формуле 0.27h, равна соответственно 405–621 МПа. При этом значения максимальных касательных напряжений на порядок ниже и находятся в пределах 0.175–120 МПа [Гзовский, 1964] соответственно для наиболее подвижных районов Средней Азии и Тянь-Шаня, Памира. В работе [Atabekov, 2020] предложен один способ усреднения, позволяющий свести уравнения (1) к средним горизонтальным напряжениям с возможностью сохранения трехмерной специфики решений. Методы решения (1) приводятся в Приложении 2. В результате получаем тензор средних горизонтальных напряжений σ₁₁, σ₁₂, σ₂₂, средние перемещения u₁, u₂ и вертикальное перемещение на свободной поверхности.

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

Отметим, что реконструированные модельные значения напряжений верны только для выбранных глубин. Для других глубин (например, 3–4 км) при проведении численного эксперимента нужно ориентироваться на известные другие значения напряжений. Естественно ожидать, что хотя картина распределения изолиний при этом качественно сильно не изменится, но в количественном отношении будут совершенно другие значения. В результате численного решения получилось, что верхняя поверхность, построенная по значениям вертикальных перемещений несильно уклоняется от реального рельефа, и находится в пределах допустимой погрешности (рис. 3).

На основе полученных решений для средних горизонтальных перемещений Центральной Азии проведем геодинамический анализ состояния земной коры.

Напряженное состояние сплошной среды в общем случае определяется тремя главными векторами напряжения (условно σ₁ > σ₂ > σ₃). Их значения и направления определяются из уравнения:

Для симметричного тензора σ_ijуравнение (10) всегда имеет три действительных корня следующего уравнения [Ильюшин, 2004]:

В нашем случае, из решения (11) мы получаем горизонтальные средние напряжения σ₁₁, σ₁₂, σ_22. Для геодинамического анализа к ним добавим вертикальное напряжение σ₃₃, принятое как вес покрывающей породы. Величины σ₃₁ и σ₃₂ выразим как силы трения через σ₃₃. Согласно “закону Байерли” [Byerlee, 1978] коэффициент фрикционного скольжения в земной коре выбирали в диапазоне 0.6–1.0. В преимущественных работах геодинамическое состояние земной коры оценивается по значениям коэффициента Лоде–Надаи

с трактовкой, что сжатию соответствует µ_σ = +1, растяжению µ_σ = –1 и сдвигу µ_σ = 0. Покажем, что данное утверждение в общем случае неверно. Этот коэффициент изначально предложен в качестве индикатора процесса простого нагружения при проведении экспериментов на пластичность тонкостенных труб на одновременное растяжение с внутренним давлением и кручением. В этих экспериментах установлено, что в противовес выше указанной трактовке равномерному растяжению в двух направлениях соответствует µ_σ = 1, одноосному растяжению µ_σ = –1. В случае сдвига µ_σ = 0. В случае осевого сжатия µ_σ = 1 и в случае равномерного сжатия µ_σ = –1 [Ильюшин, 2004, с. 63]. Избегая дискуссионный вопрос, сформулированный в работах [Введенская, 1961; 1969; Костров, 1970; Мухамедиев, 2019] – можно ли использовать ориентацию снятых напряжений в очагах сильных землетрясений в качестве оценки тектонических напряжений, сделаем несколько замечаний по поводу указанных значений µ_σ на основе формулы (12). Значения µ_σ = +1, как видно по формуле (12), достигаются во всех случаях, когда эллипсоид напряжений является эллипсоидом вращения σ₁ = σ₂, не важно, что при этом эти напряжения сжатия или растяжения. Аналогично µ_σ = –1 при σ₂ = σ₃ независимо от вида (растяжения или сжатия). Таким образом, утверждение о геодинамическом состоянии земной коры на основе коэффициента µ_σ (сжатия µ_σ = +1, растяжения µ_σ = –1) пригодно только при одноосном сжатии или растяжении (эллипсоид превращается в прямую линии), и в общем случае надо быть осторожным. Более того, по иронии природы, как было отмечено еще в работе [Введенская, 1969], в вышеуказанных случаях эллипсоида, однозначное определение полного напряженного состояния на основе ориентации главных осей напряжений, действующих в очагах, выполнено быть не может.

Мы решили привести изолинии этого коэффициента на рис. 4 нашей задачи для того, чтобы с одной стороны продемонстрировать, что крайние значения µ_σ на нашем примере совершенно не соответствуют правилу (сжатия µ_σ = +1, растяжения µ_σ = –1), а, скорее, подтверждают факт двухстороннего сжатия µ_σ = –1 и двухстороннего растяжения µ_σ = +1). И с другой стороны, с целью сопоставления результатов геодинамического состояния земной коры, полученных другим методом. А именно по тому методу, по которому сравниваются отношение вертикального напряжения σ_ver к горизонтальным наибольшим σ₁ и наименьшим σ₂ напряжениям, предложенным еще в работе [Anderson, 1951] и использованным при построении карты напряжений по всему миру [Zoback, 1992; Zoback, Zoback, 2007; Heidbach, 2018]. В этих работах выделяются три различных случая: 1) преобладает вертикальное напряжение σ₁ = σ_ver, когда гравитация вызывает нормальные разломы, создавая горизонтальную деформацию растяжения; 2) случай σ₂ = σ_ver, при котором разница между двумя горизонтальными напряжениями создает сдвиговую деформацию; 3) оба горизонтальные напряжения превышают вертикальное напряжение σ₃ = σ_ver, при котором деформация сжатия компенсируется за счет надвига или взброса. В отличие от работ [Zoback, 1992; Zoback, Zoback, 2007; Heidbach, 2018] мы используем правило знаков напряжений, принятое в классической механике, где считается, что сжимающие напряжения отрицательны, а растягивающие – положительны. В силу этого главные напряжения σ₁ являются минимальными сжимающими, а σ₃ – максимальными сжимающими. Вследствие этого случай 1) и случай 3) в нашем примере поменяются местами.

Значения главных горизонтальных напряжений σ₁, σ₂ и максимальное касательное напряжение σ_max вычисляется по нашим решениям σ₁₁, σ₁₂, σ₂₂, по формулам:

Значения максимальных тангенциальных напряжений в 6–7 раз меньше, чем у наибольших напряжений σ₁. Ориентации σ₁ показаны направленными отрезками, которые удовлетворительно согласуются с направлениями максимального горизонтального напряжения сжатия по данным сгенерированными на сайте WSMP 2016 г. для областей, ограниченных географическими координатами (36:39 с.ш., 56:60 в.д.) и (36:44 с.ш., 66:76 в.д.) [Мухамедиев, Галыбин, 2021], а также направлениями главных современных сжимающих напряжений, построенных по записям сильных землетрясений на территории СССР, приведенных в работе [Гзовский, 1964, с. 418]. При этом возможностям реализации вариантов 1)–3) при землетрясениях соответствуют взброс (R), сдвиг (G) и сброс (B). По результатам реконструкции проведенной численной моделью, преобладающим геодинамическим типом напряженного состояния для глубин земной коры 15 км является режим горизонтального двухстороннего сжатия. В зависимости от глубин направления максимального напряжения немного изменяются. Для глубин ниже 40 км значения σ₁, σ₂, и σ_ver становятся почти одного порядка, и возможным механизмом осуществления землетрясений, с небольшим отклонением, становится сдвиг.

При проведении численных экспериментов по поиску граничных условий для краевой задачи мы опирались на установленные значения максимальных горизонтальных касательных напряжений на территориях Центральной Азии, приведенных в работах М.В. Гзовского и В.И. Уломова (см. выше). По их оценкам на глубинах 15–20 км значения горизонтальных касательных напряжений на равнинной части Центральной Азии не превышают 175–200 МПа, а на Тянь-Шане могут достигать 1200 МПа. При решении задачи мы искали граничные условия, которые дают такие напряжения. Поэтому наша модель соответствует значениям напряжений в указанных работах и верна настолько, насколько верны их расчеты. Наша модель также учитывает конкретные разломные зоны в качестве границ блоков, на которые разделена наша расчетная область. На контактах задавались непрерывность нормальных напряжений и касательные напряжения в виде трения. Надо отметить, что в наших расчетах получаются разные результаты, когда вся территория рассматривается как одна зона или делится на несколько зон. Только благодаря подбору различных упругих констант для различных блоков удалось получить те значения напряжений, которые соответствуют работам [Гзовский, 1964; Уломов, 1974]. Наши результаты вертикальных перемещений в Копетдаге соответствуют средним значениям за рассматриваемый период 10⁴ год, принятый при выборе коэффициента вязкости в модели. Поэтому их можно использовать при рассмотрении относительно больших интервалов времени. Расхождения модельных и геодезических измерений за короткий период могут вполне не сойтись один к одному, как и нет соответствия между GPS и геодезическими измерениями [Кузьмин, 2017] Однако, используя полученные напряжения, мы можем построить локальные модели и сделать более конкретные расчеты для коротких периодов. В данной статье мы используем их в качестве статических напряжений для решения задачи о вариации напряжений при землетрясениях с предлагаемыми нами эквивалентными силами, полученными на основе моментной теории.

В географической системе координат тензор m_ij преобразуется в $m_{{ij}}^{'}$ через параметры нодальной плоскости к виду:

В случае примера с параметрами (δ = 90°, λ = 0°, φ = 0°) в правую часть усредненных уравнений (1) прибавляются усредненные производные от m_13.2 и m_31.1, имеющие следующий вид:

Здесь: Δσ – сброшенное напряжение; r₀– радиус сферического объема очага землетрясения; n_eff – коэффициент затухания, который определяется эмпирически [Ризниченко, 1974].

В качестве примера применения предлагаемого способа механизма очага землетрясения, введенного по формуле (16), приведен расчет вариации максимальных касательных напряжений после землетрясения в Ташкенте 1966 г. (69°21′ в.д.; 41°26′ с.ш.) для упрощенного варианта расположения нодальной плоскости с параметрами φ = 90°, δ = 90°, λ = 0° (рис. 6). Согласно измеренным данным в работе [Ташкентское …, 1971] объем землетрясения принят в виде сферы с объемом 62 км³ и сброшенное напряжение, равным ∆σ = 6 × 10⁶ Па.

На рис. 7 сопоставлены значения напряжений, получаемые по безмоментной [Atabekov, 2021] и по моментной теории упругости. Их относительные значения по всей территории не превышают 0.5%, не считая территории вокруг эпицентра.

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

ЗАКЛЮЧЕНИЕ

Предложена математическая модель механизма очага землетрясения на основе моментной теории упругости.

При анализе экспериментальных опытов при сложном напряженном состоянии подвергается сомнению метод геодинамической оценки напряжений земной коры (сжатие, растяжение, сдвиг) на основе коэффициента Надаи–Лоде.

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

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