Физика Земли, 2021, № 1, стр. 122-132

ВВЕДЕНИЕ

Универсальные методы для выделения областей готовящихся будущих землетрясений еще не найдены. Однако, вероятней всего, эти области окажутся там, где происходят интенсивные изменения деформационных процессов за реальный промежуток времени. Одним из важных параметров при сейсмическом районировании является сейсмотектоническое течение. Рассматриваемый регион, ограниченный географическими координатами 36:46 с.ш., 56:76 в.д., по известной гипотезе находится в тисках платформенных массивов и интерпретируется как результат раздавливания коры этого пояса в поле латерального сжатия, порожденного коллизией Евразийской с одной стороны, Индийской и Аравийской литосферных плит с другой. Для решения поставленных задач используются уравнения равновесия моментной теории. Механизм очага землетрясения моделируется двумя способами: сосредоточенным внутренним моментом и двойным диполем без момента. Трехмерные уравнения сплошных сред усредняются по толщине литосферы с использованием свойств геодинамической постановки задачи. Правые части усредненных уравнений сейсмотектонического течения содержат тектонические напряжения литосферы и рельеф земной коры. В качестве таковых используются напряжения и рельеф, получаемые решением обратной задачи упругости. С этой целью разыскиваются возможные граничные условия для упругой задачи, создающие современный рельеф земной коры. Полученные результаты сравниваются с полями скоростей, построенными на основе данных GPS.

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

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

где: σ_ij, µ_ij – компоненты силового и моментного тензора напряжений; индекс j после запятой означает дифференцирование по декартовым координатам x_j; F(0, 0, ρg) – массовая сила; ε_ijk – тензор Леви-Чивита; M_i – удельный момент, имеющий размерность момент, деленный на объем.

Основная сложность в применении моментной теории заключается в трудностях, возникающих при определении констант, связывающих обобщенные напряжения с кинематическими параметрами. Известно небольшое количество экспериментов, позволяющих идентифицировать шесть упругих констант Коссера лишь для простейших материалов. Учитывая эту неопределенность относительно констант, можно выразить вектор вращения как и в безмоментной теории ω_k = ε_ijku_i,j. Тогда µ_ij = 0 и формула (2) выражает асимметричность тензора напряжения. В случае точечного очага, смоделированного внутренним моментом, с координатами (x_i0i =1, 2, 3) M_i выражается через дельта функций:

Здесь $M_{i}^{0}$ = const – значение удельного момента. Однако в случае моделирования точечного очага двойной диполью без момента, в уравнение (2) в точке x_i0 приходится ввести два момента, равные по величине, но противоположные по направлению, которые в сумме дают нуль, что затрудняет использование уравнения (2). В этом случае применение теории дислокации дает возможность [Ландау и др., 1965; Аки и др., 1983; Касахара, 1985] введения дополнительных членов в объемные силы F_i в уравнениях (1), эквивалентных действиям диполя. Эти силы по-разному выражаются в зависимости от вида подвижки в очаге землетрясения (сдвиги, сбросы и т.д.). В нашу модель вводили соответствующие члены, исходя из конкретного вида землетрясения.

Решение трехмерных уравнений довольно сложно, поэтому предлагаем один способ упрощения системы (1)–(2) с учетом специфики геодинамической постановки задач, а именно: размеры рассматриваемой территории в плане намного больше по сравнению с толщиной литосферы; поверхность земли свободна от напряжений; на контакте плит с астеносферой касательные напряжения отсутствуют. Эти особенности задачи дают возможность усреднять уравнения (1)–(2) по формуле:

где: черта означает усреднение по x₃ (принята географическая система координат с вертикальной осью, направленной вниз); h(x₁, x₂) – глубина подошвы литосферы, H(x₁, x₂) – высота поверхности Земли.

Усредним систему (1)–(2), принимая ${{\sigma }_{{ij}}} = {{\sigma }_{{ji}}} + {{\varepsilon }_{{ijk}}}{{M}_{k}}$ (i, j, k – круговое индексирование) и обезразмеривая напряжение и удельный момент относительно модуля сдвига G, а линейные переменные относительно h:

При выводе уравнений (5) учтены условия на поверхности Земли и на подошве астеносферы:

а также результаты усреднения дельта-функций:

Согласно правилам интегрирования дельта-функций, при i = 1, 2 имеем

Интеграл в (7) обращается в нуль благодаря антисимметричности производной δ-функции при H/h < x₃₀ < 1.

Система в напряжениях (5) является общей для задач упругости и линейно-вязкой жидкости. Для вязкой жидкости система преобразуется в уравнения Стокса, которые относительно усредненных скоростей перемещения ${{\bar {v}}_{i}}({{x}_{1}},{{x}_{2}})$ имеют следующий вид:

Усредненное уравнение несжимаемости в предположении v₃(x₁, x₂, h) = 0 принимает вид:

и выражает вертикальную скорость поверхности Земли.

Правые части уравнений (8)–(9) содержат усредненные тектонические напряжения ${{\bar {\sigma }}_{{ij}}}$ и рельеф поверхности H(x₁, x₂). Вообще-то говоря, для их определения нужно решать задачу пластичности с учетом особенностей всего тектонического периода, что невозможно из-за сложности процесса и неизвестности физических параметров. Вместо этого в качестве тектонических напряжений мы использовали обычное решение упругой литосферы, деформация которой создает вертикальные перемещения, соответствующие современному рельефу.

В упругой постановке система (1)–(2) с M_i = 0 преобразуется в уравнения Ламе для усредненных перемещений:

где правая часть имеет следующие компоненты:

Здесь: Δ – оператор Лапласа; ν – коэффициент Пуассона; U – двумерный вектор с компонентами усредненных горизонтальных перемещений; $u_{3}^{h}$, $u_{3}^{H}$ – вертикальные перемещения на подошве астеносферы и на поверхности земной коры.

При выводе уравнений (10) использованы следующие соотношения для i = 1, 2:

Система уравнений (10) решается методом итерации. За первое приближение средние напряжения и разность $u_{3}^{h} - u_{3}^{H}$ приняты нулевыми. В последующих итерациях эта разность находится по следующей формуле, полученной усреднением напряжения σ₃₃ через перемещения (λ, μ – упругие константы):

Отметим некоторые особенности литосферы исследуемого региона. В работе [Абидов, 1997], где обобщены труды коллектива Института геологии и разведки нефтегазовых месторождений (г. Ташкент), в том числе и автора данной статьи, приведены параметры литосферы Центральной Азии. При анализе имеющейся геолого-геофизической и сейсмологической информации для данного региона установлено, в частности, что литосфера состоит из 19-ти блоков первого ранга (рис. 1). Рассчитаны скорости Индийской и Арабской плит относительно Евразийской, которые оказались равными соответственно 46 и 25 км в 1 млн лет. Поэтому для решения рассматриваемых задач в качестве литосферы принято призматическое тело, состоящее из 19-ти зонально однородных частей. Скорости Индийской и Аравийской плит для первого приближения выбирались согласно вычисленным нами значениям, а направления действующих сил p_j = σ_ijn_i (n_i – компоненты нормали к границе) такими, как указаны стрелками на рис. 1 (последующими приближениями мы их корректировали). В каждой итерации уравнения (10) решались методом граничных интегральных уравнений относительно компонентов перемещения и вектора напряжений.

Соотношения метода граничных интегральных уравнений имеют стандартный вид [Бреббия и др., 1987]:

где: b_j – компонента вектора, представляющая собой правую часть (10); $p_{{ij}}^{*}$, $u_{j}^{*}$ есть фундаментальное решение уравнения Ламе (черта над усредненными переменными опущена):

Здесь: $S = \bigcup\nolimits_{i = 1}^n {{{S}_{i}}} $;$\Omega = \bigcup\nolimits_{i = 1}^n {{{\Omega }_{i}}} $; S_i – граница двумерной области Ω_i; r_j = x_j – ξ_j.

S_i разбиты на линейные элементы, Ω_i на треугольники, вершинами которых являются внутренние и граничные точки. Дискретизация уравнения (14) с учетом специфики рассматриваемой задачи приводит к линейной алгебраической системе относительно неизвестных значений u_j и p_j. В объемных интегралах с сингулярной особенностью интегрирование проводится полуаналитически, с переходом на полярную систему координат с началом в особой точке. Решение обратной задачи, как правило, некорректно и соответствующая матрица плохо обусловлена. Она регуляризирована известным методом А.Н. Тихонова [Тихонов и др., 1979]. Исходную систему с матрицей H и правой частью F можно переписать в эквивалентной форме:

Регуляризованная задача формулируется следующим образом:

где: α – малый положительный параметр. Варьируя α, получаем следующее уравнение: где Е – единичная матрица, H^H – эрмитова сопряженная матрица. Решая его методом исключения Гаусса, найдем регуляризованное решение Y_α, зависящее от параметра α. Если α = 0, то система переходит в плохо обусловленную систему. Если же велико, то регуляризованная система будет хорошо обусловленной благодаря присутствию в левой части хорошо обусловленной матицы αЕ; но преобразованная система при большом α сильно отличается от исходной, и регуляризованное решение Y_α не будет близким к искомому решению. Оптимальным будет наименьшее значение α, при котором обусловленность системы еще удовлетворительна.

Для нахождения оптимального значения α вычисляли невязку r_α= HY_α– F и сравнивали ее по норме с известной погрешностью правых частей δF и с влиянием погрешности коэффициентов матрицы δH ⋅ Y. Оптимальным считали тот, в котором выполнялось условие ||r_α|| ≈ ||δF|| + ||δH ⋅ Y||.

По значениям u_j вычисляются σ_ij, и по ним наибольшие и наименьшие напряжения во внутренних точках. На границах блоков касательное σ_τ и нормальное σ_n напряжения находятся по формулам σ_n = p₁n₁+ p₂n₂, σ_τ = –p₁n₂+ p₂n_1. Численный эксперимент заключался в следующем. Граничные перемещения, соответствующие Евразийской плите на рис. 1, принимались нулевыми (Т-образный символ). В первом приближении компоненты вектора напряжений на нижних углах четырехугольника выбирались пропорционально скоростям Индийской и Аравийской плит, а между ними равномерно распределялись с левого до правого угла. На боковых границах компоненты плавно уменьшались с приближением к верхней границе. Вектор напряжений p на верхней границе и перемещения u на остальных границах считались искомыми.

После нахождения недостающих граничных переменных, перемещения u во внутренних точках определяются опять по уравнениям (11). Для следующего шага вычислительного эксперимента, предполагая u₃(x₁, x₂, h) = 0 находится u₃(x₁, x₂, 0) = H(x₁, x₂) из (10) и сравнивается со значением рельефа физической карты. Далее интуитивно варьируются напряжения на границах и упругие константы блоков до достижения приемлемой формы реальных высот рельефа. Найденные окончательные значения ${{\bar {\sigma }}_{{ij}}}$ и H(x₁, x₂) использовались для решения сейсмотектонического (M_i = 0) и сейсмического (M_i ≠ 0) течения. Уравнения (8) также решаются вышеописанным способом. В данном случае фундаментальные решения задачи Стокса (относительно усредненных v_j и f _j = σ_ijn_i) имеют следующий вид [Ладыженская, 1970]:

После нахождения усредненных скоростей перемещения усредненные напряжения строятся по формулам:

Для проведения расчетов составлена авторская программа в среде Delphi. Изолинии построены с помощью пакета Surfer.

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

На рис. 2 приводятся окончательный вариант рельефа (рис. 2а) из множества решений, полученных в ходе вычислительного эксперимента упругой задачи. Рисунок 2б соответствует современному рельефу, построенному по физической карте, которая служила для верификации решения математической задачи. Для получения окончательного варианта задачи, которая наиболее близко соответствует верифицируемой, приходилось последовательно корректировать значения граничных напряжений и упругих констант каждого из 19-ти блоков со средними значениями G = 0.31 × 10⁵ МПа и ν = 0.25. Максимальная высота, равная 7500–8000 м, получилась при значениях безразмерных компонент вектора напряжения в правом нижнем углу, равных 0.0005G, 0.035G и угол α = arctg(f₂/f₁)= 89.1°, отсчитываемый от параллели против часовой стрелки, который вполне соответствует направлению движения Индийской плиты, установленному по данным GPS в работе [Зубович и др., 2007].

В работах [Гзовский, 1964; Николаев и др., 1967] проведены оценки максимальных касательных напряжений на глубине 15–20 км для Центральной Азии 400–1500 бар. Учитывая [Ильюшин, 1948] известные соотношения 0.816 < σ_int/σ_max < 0.941 между максимальными касательными напряжениями σ_max и интенсивностью напряжений σ_int можно утверждать, что приведенные на рис. 3 величины интенсивности напряжений, рассчитанные по нашей модели, неплохо коррелируются с их оценками.

Для анализа результатов тектонического (M_i = 0) и сейсмотектонического (M_i ≠ 0) течения горных масс значение вязкости варьировалось около 10²² Н с/м². Масштаб времени выбирался из равенства безразмерных касательных напряжений в уравнениях Ламе и Стокса:

Полученное значение, равное 10⁴ лет, вполне согласуется с тем фактом, что максимальное время релаксации напряжений горных пород измеряется отрезком не более 10⁴ лет [Тёркот, 1985].

Скорости перемещений при M_i = 0, рассчитанные по нашей модели (рис. 4–рис. 5, в см/год), по направлениям и по значениям удовлетворительно совпадают с результатами других авторов.

Здесь надо оговорить, что наши расчеты отражают скорости по отношению 46° с.ш. (на этой широте скорости модели принимались нулевыми), когда в большинстве построений по данным GPS на этой параллели имеются еще скорости, направленные на север. Только, в работе [Современная…, 2005] представлен случай, где данные для Ферганской впадины пересчитаны относительно центра впадины. Поэтому наши картины неплохо соответствуют построениям этой работы. При учете реальных скоростей, направленных на север на параллели 46° с. ш., векторы скоростей также стали приблизительно совпадать с результатами остальных работ [Зубович и др., 2007; Кузиков и др., 2010; Костюк и др., 2010].

Интенсивность тектонического течения и скорости вращения определялись следующими формулами:

Интенсивность (22) выбрана аналогично интенсивности сдвига в теории упругости, которая равна с точностью до множителя второму инварианту тензора деформации. Такой выбор коэффициента в упругом случае обусловлен тем, что интенсивность деформации должна равняться удвоенному углу сдвига площадки результирующих сдвигов. Такая площадка наклонена равным образом к главным осям деформации. В нашем случае под e_ij подразумеваются средние скорости деформации, определенные по средним скоростям (рис. 6), при этом учитываются формулы (12)–(13).

Наибольшие значения интенсивности (~1.6 × × 10^–7 1/год) отчетливо наблюдаются (рис. 6) в Ферганской впадине вдоль Восточно-Ферганской и Южно-Ферганской зон тектонических нарушений (границ блока № 12), которые отличаются высокой сейсмичностью. Наименьшие значения ~10^–8 1/год свойственны Центрально-Казахстанскому щиту. В Центральных Кызылкумах величина интенсивности скоростей деформаций возрастает до 2 × 10^–8 1/год.

Изолинии вращения (10^–7 радиан/год) на рис. 7 показывают, что территория Ферганской впадины вращается против часовой стрелки вокруг центра с географическими координатами 40.6° с.ш., 69.9° в.д. со скоростью угла вращения, равном 4 нанорадиан/год. Незначительное вращение наблюдается и на части Таджикской депрессии, только в обратном направлении. Полученные нами числовые значения угловой скорости сравнивались с аналогичными значениями, построенными по данным GPS для Ферганской впадины в работе [Зубович и др., 2010]. В этой работе угловая скорость получена около 12–13 нанорадиан/год. Такая заметная разница при сравнении результатов нашей модели и фактических данных GPS заставила нас скорректировать нашу модель. При уменьшении вязкости для Ферганского блока (№ 12 в рис. 1) на один порядок получились вращения, вполне сопоставимые с данными упомянутой работы.

Вертикальные скорости на поверхности Земли в масштабе мм/год, построенные по формуле (9), приведены на рис. 8. Они почти везде восходящие. Незначительные нисходящие движения (0.1–0.5 мм/год) имеются в пределах Таджикской депрессии и на северной территории, начиная от Чуйской долины. Области резкого изменения градиента скорости могут дать ценную информацию для выявления опасных зон при составлении карт районирования. Для этой цели можно использовать горизонтальную компоненту средней скорости по формулам (14) для задачи Стокса и вертикальную скорость на поверхности Земли по формуле (9). На рис. 9 приводится сеймотектоническое течение, т.е. относительное изменение интенсивности скоростей деформации (в %) вследствие Ташкентского землетрясения (26.04.1966, М = 5.25, h = 8 км, φ = 41.35° N, λ = 69.27°E; нодальные плоскости: I – азимут – Az = 215°, угол падения – e = 28°; II – Az = 40°, e = 61°; объем очага землетрясения V = 62 км³; амплитуда смещения 50 см [Ташкентское…, 2071]). После землетрясения картина тектонического течения немного (менее 2%) изменилась поблизости небольшой территории рядом с очагом землетрясения, почти не меняясь на остальных территориях, что является результатом небольшой глубины гипоцентра.

ЗАКЛЮЧЕНИЕ

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