Доклады Российской академии наук. Науки о Земле, 2022, T. 506, № 2, стр. 252-258

Сопутствующие землетрясениям эффекты не ограничиваются сильными движениями твердой среды, изменением ее механических свойств в очаговой зоне события и формированием сейсмического сигнала. Землетрясения, особенно сильные, вызывают волновые движения в атмосфере и заметные вариации геофизических полей – электрического и магнитного, причем на значительных эпицентральных расстояниях [1, 2]. Изучение геомагнитных вариаций при сейсмических событиях представляет особый интерес, во-первых, с точки зрения их всестороннего описания и, во-вторых, – определения характеристик межгеосферных взаимодействий в системе земная кора–атмосфера–ионосфера [3].

В качестве основного механизма возбуждения геомагнитных вариаций при землетрясениях рассматривается воздействие вызванных в результате вертикальных движений земной коры атмосферных волн¹¹ на ионосферу в эпицентральной области сейсмического события [4–6].

В настоящем сообщении предлагается численная модель развития движений воздушных масс при землетрясении с последующим формированием источника геомагнитных возмущений на ионосферных высотах. В качестве примера приведены результаты расчетов, выполненных для сильного землетрясения 11.03.2011 г. (М = 9.1, Япония) [7].

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

Для описания параметров атмосферных волн решалась система дифференциальных уравнений газовой динамики, записанных в переменных Эйлера [8]:

где $\rho $ – плотность воздуха, ${\vec {v}}$ – вектор скорости, t – время, p – давление, $\vec {g}$ – ускорение свободного падения, $\varepsilon $ – внутренняя энергия, T – температура.

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

Вертикальная u и радиальная v компоненты вектора скорости ${\vec {v}}$ представлялись в виде:

где ${{z}_{0}}$ и ${{z}_{f}}$ – соответственно начальное и конечное значение вертикальной координаты сетки, ${{r}_{0}}$ и ${{r}_{f}}$ – начальное и конечное значение радиальной координаты, А – амплитуда волны, ${{t}_{o}}$ – период волны, ${{r}_{s}}$ – радиус области генерации волны.

Газодинамические параметры атмосферы задавались в соответствии с ATMCIRA²².

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

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

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

В качестве примера на рис. 1 и 2 приведены результаты расчета относительного давления р/р₀ и массовой скорости частиц среды u в атмосфере в зависимости от высоты z для момента времени t = = 450 с (расчеты выполнены для значений параметров: А = 10 м/с, t₀ = 1 c, r_s = 10 км).

Воздействуя на ионосферу, атмосферное возмущение вызывает систему электрических токов в динамо-области и как результат – вариации магнитного поля [4, 10, 11].

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

где$\vec {j}$ – плотность тока, $\vec {H}$и $\vec {D}$ – соответственно напряженность магнитного поля и индукция электрического поля, ρ – плотность заряда, ${{\vec {j}}_{E}}$ – плотность сторонних токов.

В проводящей ионосфере:

Согласно закону Ома для плазмы [12]:

Здесь $\vec {B}$ – вектор индукции магнитного поля, $\vec {b}$ – единичный вектор в направлении магнитного поля, Е и E_d – соответственно напряженность электрического и динамо-поля, ${{\sigma }_{0}}$ – удельная проводимость вдоль магнитного поля, ${{\sigma }_{p}}$ и ${{\sigma }_{H}}$ – соответственно удельные проводимости Педерсена и Холла поперек магнитного поля, ${v}$ – вектор скорости смещения частиц среды, который определялся в результате решения газодинамической задачи.

Рассматривая для конкретности Е-слой ионосферы, в котором продольная проводимость превышает поперечные, плотности токов $\vec {j}$ и проводимости $\sigma $ интегрируются по высоте динамо-слоя.

Поскольку

где U – потенциал электрического поля, то уравнение (3) для полного тока сводится к уравнению второго порядка в частных производных для U.

При решении полученного уравнения в конечных разностях применялась схема Кранка-Николсона [8] при заданных нулевых значениях электрического потенциала на границах. Численные расчеты на этом этапе выполнялись в декартовых координатах x, y, z (координата z направлена вверх). Корректность работы расчетной программы проверялась решением уравнения Пуассона для заряженного шара.

Напряженность магнитного поля определяется из уравнения Максвелла (4) и условия (6) с применением операции rot:

Последнее выражение можно представить в виде:

Компоненты напряженности магнитного поля H_x, H_y и H_z, определялись путем решения системы уравнений:

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

где $\mu {}_{0}$ – магнитная постоянная, $\mu $ – относительная магнитная проницаемость среды.

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

Основной интерес представляют вызванные геомагнитные возмущения в приповерхностной атмосфере (вблизи земной поверхности), где в подавляющем большинстве случаев осуществляются инструментальные наблюдения. В качестве примера на рис. 3 приведены результаты расчета компонент индукции возмущений магнитного поля на земной поверхности для одного из моментов времени (t = 360 с), полученные с использованием результатов решения газодинамической задачи при заданных выше значениях параметров А, t₀ и r_s. Из рис. 3 следует, что сформированный в эпицентральной области ионосферы источник магнитных возмущений вызвал заметные геомагнитные вариации вблизи земной поверхности амплитудой ~0.35–0.5 нТл.

С целью верификации полученных результатов расчета проводилось их сравнение с данными инструментальных наблюдений, выполненных при Великом Японском землетрясении 11.03.2011 г. Основной толчок, определивший мегакласс события, произошел в 05:46 UTC. Размеры области выделения энергии (очага землетрясения), оцененные по облаку афтершоков первых суток, оцениваются как ~200 × 500 км [12]. Землетрясение вызвало сильные деформации среды, мощный сейсмический сигнал и цунами. Одновременно многочисленными магнитными станциями сети INTERMAGNET в период землетрясения были зарегистрированы вариации магнитного поля [10]. В частности на 4 ближайших к эпицентру землетрясения станциях MMB, КАК, KNY и CYG (рис. 4). В качестве примера на рис. 5 представлены записи вертикальной B_z и основной, наиболее чувствительной к внешним возмущениям горизонтальной компоненты индукции магнитного поля B_x (направление юг–север), полученные в период землетрясения на станциях КАК и ММВ, расположенных соответственно на эпицентральных расстояниях ~300 и ~640 км.

Из рис. 5 следует, что сильное сейсмическое событие вызвало через ~ 10 мин изменение характера вариаций магнитного поля: появились квазипериодические колебания амплитуды с периодом ~40 мин. С точки зрения причинно-следственной связи между наблюдаемыми вариациями магнитного поля и рассматриваемым событием важно отметить, что, во-первых, наблюдаемый период вызванных геомагнитных вариаций существенно превосходит период геомагнитных пульсаций Рс5 космического происхождения (период 2.5–10 мин), и, во-вторых, колебания с периодом около 40 мин не зарегистрированы ни до, ни после рассматриваемого сейсмического события. Амплитуда первой фазы вызванных вариаций магнитного поля, которую разумно рассматривать как отклик на воздействие атмосферного сигнала, составляет для компоненты B_x около 4–5 нТл относительно тренда. Максимальная амплитуда вариаций B_x относительно тренда достигает 7–8 нТл. Вариации вертикальной компоненты B_z оцениваются в ~3 нТл.

Оценки, выполненные на основе предложенной расчетной модели с использованием данных о параметрах движения среды и характерных размерах активной области земной коры, определяющей основной вклад в вертикальные движения земной поверхности (А = 10 м/с, t_w = 10 c, r_w = 20 км) [7, 13–15], дают значения амплитуды геомагнитных вариаций вблизи земной поверхности в диапазоне 2–4 нТл, что вполне согласуется с приведенными выше результатами инструментальных наблюдений.

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