Электромагнитные УНЧ поля на земной поверхности и в ионосфере от подземного сейсмического источника

Полный текст

1. ВВЕДЕНИЕ: УНЧ ПОЛЯ СЕЙСМИЧЕСКИХ ТОКОВЫХ СИСТЕМ

Одно из актуальных направлений в современной геофизике – разработка физических основ оперативного (за дни–недели) прогноза землетрясений по аномальным возмущениям электромагнитных полей и ионосферы, дополняющего стандартные сейсмические методы [Бучаченко и др., 1996]. Ключевую роль в сейсмоэлектромагнитных методах играют излучения разных частотных диапазонов. В настоящее время перспективным представляется мониторинг излучений ультранизкочастотного (УНЧ) диапазона $({10}^{-2}-10$ Гц), источником которых могут быть динамические процессы в области очагов землетрясений. Такими процессами могут быть разделение зарядов при трещинообразовании, иррегулярное течение подземных флюидов через поры в горной породе, акустические импульсы раскрывающихся трещин [Surkov, Hayakawa, 2014]. Однако оценки электромагнитных УНЧ полей, создаваемых системой хаотично ориентированных трещин, показали, что статистический эффект таких микроизлучателей на земной поверхности оказывается слишком малым [Molchanov, Hayakawa, 1995; Surkov, Hayakawa, 2006]. В то же время возможно образование крупномасштабных токовых систем, обусловленных движением горной среды вдоль разломов при активизации сейсмической активности [Lockner et al., 1983; Гохберг и др., 1985; Гульельми, Левшенко, 1997]. Эти крупномасштабные токи могут иметь электрокинетическую или трибоэлектрическую природу либо могут быть связаны с движением заряженных дислокаций. Попытки мониторинга таких движений ведутся с помощью регистрации УНЧ импульсов и шумов на земной поверхности [Hattori, 2004; Freund et al., 2021].

Более того, предпринимаются попытки обнаружить сейсмогенные УНЧ возмущения на низкоорбитальных спутниках [Kodama et al., 2000]. Постоянно растет число сообщений об “аномальных” электромагнитных возмущениях в УНЧ диапазоне, обнаруженных спутниками в верхней ионосфере [Picozza et al., 2021]. Во всех работах спутниковые данные анализировались только для пролетов над изучаемым регионом в ночные часы. Электрические возмущения на частотах <10 Гц, возможно, связанные с предстоящим землетрясением, были обнаружены еще на ранних спутниках ИКБ-1300 и ОРЕОЛ-3. Эти обнадеживающие результаты стимулировали разработку специализированных спутниковых миссий для обнаружения сейсмоэлектромагнитных УНЧ излучений: DEMETER (высота орбиты ∼660 км) [Parrot, 1990; 2018]; CSES (высота орбиты ∼500 км) [Zhima et al., 2022]; ESPERIA [Sgrigna et al., 2008]; TwinSat [Chmyrev et al., 2013].

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

2. ТЕОРЕТИЧЕСКИЕ ОСНОВЫ МОДЕЛИ

Предположим, что механоэлектрические преобразования в земной коре создают электродвижущую силу, генерирующую ток $J\left(t\right),$ который затем замыкается токами проводимости в земной коре. Задача заключается в том, чтобы найти создаваемые этой токовой системой электромагнитные поля как в атмосфере, так и в ионосфере на спутниковых высотах. Рассматриваемая задача не сводится к классическим задачам об электромагнитном излучении токового источника, заглубленного в проводящее полупространство [Baños, 1966; King et al., 1981]. В рассмотренной модели система осциллирующих токов в земной коре самосогласованно связана с возбуждаемыми ею электромагнитными полями $B(x,y,z,t)$ и $e(x,y,z,t)\equiv c^{-1}E$ (где с – скорость света) в атмосфере и ионосфере.

Геометрия модели показана на рис. 1: ось z декартовой системы координат направлена вертикально вверх (z = 0 на поверхности Земли), ось x – на восток, и ось y – на север. Для простоты принимается, что геомагнитное поле вертикально, т.е. его наклонение I = 90°. Атмосферная проводимость увеличивается экспоненциально до высоты 80 км по закону $\sigma \left(z\right)={\sigma }_a\exp (z/z_a),$ при этом величина z_a подбирается такой, чтобы σ(z) сшивалась с проводимостью, даваемой моделью IRI. С высоты 80 км начинается ионосфера, вертикальный профиль параметров которой рассчитан с помощью моделей IRI и MSIS. Такая модель со слоистой горизонтально однородной ионосферой с реалистическим вертикальным профилем аналогична модели, использованной в работах [Fedorov et al., 2020; 2021].

Пусть нестационарный прямолинейный ток $J\left(t\right)$ конечной длины L расположен на некоторой глубине $z=h<0$ под поверхностью Земли и направлен вдоль оси x. Рассмотрим отдельную спектральную гармонику $J\left(t\right)=J_0\exp (-i\omega t).$ Начальный и наиболее важный этап решения задачи – нахождение решения уравнений Максвелла с источником $\mathit{j}$ в виде точечного горизонтального токового диполя:

$\nabla \times B=-i{k}_0\widehat{\epsilon \epsilon }e+{\mu }_{0}j,\nabla \times e=i{k}_{0}B$ (1)

Здесь $k_0=\omega /c$ – волновое число в вакууме, $\hat{\epsilon }$ – тензор относительной комплексной диэлектрической проницаемости среды. В координатной системе с осью z вдоль геомагнитного поля $B_0$ тензор$\hat{\epsilon }$  имеет ненулевые элементы ${\epsilon }_{xx}={\epsilon }_{yy}={\epsilon }_{\perp },{\epsilon }_{zz}={\epsilon }_{\parallel },{\epsilon }_{xy}=-{\epsilon }_{yx}=ig.$ Стандартные формулы для элементов тензора $\hat{\epsilon }$ можно найти, например, в книге [Гинзбург, 1967]. Даваемые моделью IRI профили параметров ионосферы дают возможность рассчитать высотные профили элементов тензора диэлектрической проницаемости. В приземной атмосфере элементы тензора ${\epsilon }_{\perp }={\epsilon }_{\parallel }\to 1$ и g → 0. Ниже поверхности Земли ${\epsilon }_{\perp }={\epsilon }_{\parallel }={\epsilon }_g=\mathrm{Re}{\epsilon }_g+i{\sigma }_g/\left(\omega {\epsilon }_0\right)$ и g = 0, где ${\sigma }_g$ и Reε_g – проводимость и диэлектрическая проницаемость Земли.

 

Рис. 1. Схематическое изображение геометрии задачи о подземном токовом источнике c током ${\mathit{J}}_{\mathbf{0}}$ на глубине h при наличии атмосферы и ионосферы.

 

При решении поставленной задачи использовано представление электромагнитного поля через потенциалы. В этом формализме наиболее эффективно удается разбить поле на потенциальную и вихревую составляющие. Это важно при решении данной задачи, которая, очевидно, лишена осевой симметрии; но оказывается, что потенциальная и вихревая составляющие в отдельности такой симметрией обладают. Благодаря этому становится возможным, разделяя переменные с помощью преобразования Ганкеля, прийти к краевой задаче для системы обыкновенных дифференциальных уравнений. Более подробно математический формализм изложен в работе [Федоров и др., 2023]. Итак, электрическое и магнитное поля стандартно выражаются через скалярный и векторный потенциалы $\Phi$ и $A:$

$B=\nabla \times A,e=-\nabla \Phi +i{k}_{0}A.$

Если воспользоваться условием калибровки ${\nabla }_{\perp }\cdot A_{\perp }=\mathrm{0,}$ то поперечная компонента векторного потенциала может быть представлена в виде $A_{\perp }={\left(i{k}_0\right)}^{-1}\nabla \times \Psi \hat{z},$ где $\Psi$ – скалярный магнитный потенциал. В результате имеем представление электромагнитного поля через потенциалы $A=A_z,$ $\Phi$ и $\Psi :$

$\begin{array}{c}B=\nabla A\times \hat{z}+{\left(i{k}_0\right)}^{-1}{\nabla }_{\perp }{\partial }_z\Psi -{\left(i{k}_0\right)}^{-1}{\nabla }_{\perp }^2\Psi \hat{z},\\ e=-\nabla \Phi +i{k}_{0}A\hat{z}+\nabla \Psi \times \hat{z}.\end{array}$ (2)

Использование этого представления приводит систему уравнений Максвелла (1) к следующей системе уравнений, определяющей эти потенциалы:

$\begin{array}{c}{\nabla }_{\perp }^2{\partial }_{z}A=i{k}_0{\epsilon }_{\perp }{\nabla }_{\perp }^2\Phi +k_{0}g{\nabla }_{\perp }^2\Psi +{\mu }_0{\text{Div j}}_{\perp },\\ {\left(i{k}_0\right)}^{-1}{\nabla }_{\perp }^2{\nabla }^2\Psi =-k_{0}g{\nabla }_{\perp }^2\Phi +i{k}_0{\epsilon }_{\perp }{\nabla }_{\perp }^2\Psi +{\mu }_0{\text{Rot\hspace{0.17em}j}}_{\perp },\\ {\nabla }_{\perp }^{2}A=i{k}_0{\epsilon }_{\parallel }(-{\partial }_z\Phi +i{k}_{0}A)-{\mu }_0{j}_{\parallel }.\end{array}$ (3)

Здесь использованы обозначения 2D-операторов $Diva={\nabla }_{\perp }\cdot a={\partial }_x{a}_x+{\partial }_y{a}_y$ и $Rota={({\nabla }_{\perp }\times a)}_z=$${\partial }_x{a}_y-{\partial }_y{a}_x.$

Рассмотрим элемент горизонтального тока – точечный токовый диполь, расположенный в точке с координатами $x=\mathrm{0,}y=\mathrm{0,}z=h.$ Тогда в уравнениях (3) имеем $j_{\left|\right|}=\mathrm{0,}$ а выражение для поперечной плотности тока в декартовых координатах имеет вид $j_{\perp }=j_x\hat{x}=M_0\delta \left(x\right)\delta \left(y\right)\delta (z-h)\hat{x},$ где $M_0=J_{0}L$ – токовый момент. Для дальнейшего анализа перейдем к цилиндрическим координатам $\rho ,\phi ,z,$ в которых выражения для входящих в уравнения (3) источников имеют вид $Div{j}_{\perp }=q\left(\rho \right)\delta (z-h)\cos \phi$ и $Rot{j}_{\perp }=-q\left(\rho \right)\delta (z-h)\sin \phi ,$ где $q\left(\rho \right)=M_0{\partial }_{\rho }\left[{\left(2\pi \rho \right)}^{-1}\delta \right(\rho \left)\right].$ Слагаемое $Div{j}_{\perp }$возбуждает потенциальную часть поля, тогда как слагаемое $Rot{j}_{\perp }$ создает его вихревую часть.

Поскольку неоднородность в системе уравнений (3) представляет собой линейную комбинацию векторов с коэффициентами $\cos \phi$ и $\sin \phi ,$ то и решение следует искать в виде таких же комбинаций $F(\rho ,\phi ,z)=F_c(\rho ,z)\cos \phi +F_s(\rho ,z)\sin \phi ,$ где $F=A,\Phi ,\Psi .$ В результате подстановки этих комбинаций в (3) и группировки членов получается система уравнений:

$\begin{array}{c}{\partial }_z{A}_{c,s}=i{k}_0{\epsilon }_{\perp }{\Phi }_{c,s}+k_{0}g{\Psi }_{c,s}+a_{c,s},\\ {\partial }_z{B}_{c,s}=-k_{0}g{\Phi }_{c,s}-{\left(i{k}_0\right)}^{-1}(k_0^2{\epsilon }_{\perp }+R){\Psi }_{c,s}-b_{c,s},\\ {\partial }_z{\Phi }_{c,s}=i{k}_0[1+{\left(k_0^2{\epsilon }_{\parallel }\right)}^{-1}R]A_{c,s},\\ {\partial }_z{\Psi }_{c,s}=i{k}_0{B}_{c,s},\end{array}$ (4)

где: оператор $R={\rho }^{-1}{\partial }_{\rho }\rho {\partial }_{\rho }-{\rho }^{-2};$ неоднородные члены $a_c=b_s={\mu }_0{R}^{-1}q\left(\rho \right)\delta (z-h)$ и $a_s=b_c=0;$ для симметрии введена дополнительная искомая функция $B_{c,s}={\left(i{k}_0\right)}^{-1}{\partial }_z{\Psi }_{c,s}.$ Принципиально важным является то обстоятельство, что система уравнений (4) уже обладает цилиндрической симметрией, а несимметричность исходной задачи преодолена благодаря разделению на две составляющие, потенциальную (с индексом c) и вихревую (c индексом s), имеющие разную угловую зависимость $(\cos \phi$ и $\sin \phi ).$ Это позволяет разделить переменные и прийти к одномерной краевой задаче.

Система уравнений (4) с частными производными приводится к системе обыкновенных дифференциальных уравнений по переменной z с помощью преобразования Ганкеля первого порядка  которое переводит дифференциальный оператор R в умножение на $-k^2.$ После применения преобразования ${\mathbb{K}}_1$ перейдем к переменным ${\tilde{A}}_{c,s}=k{\mathbb{K}}_1\left[A_{c,s}\right],$ ${\tilde{B}}_{c,s}=k{\mathbb{K}}_1\left[B_{c,s}\right],$ ${\tilde{\Phi }}_{c,s}=k{\mathbb{K}}_1\left[{\Phi }_{c,s}\right]$ и ${\tilde{\Psi }}_{c,s}=k{\mathbb{K}}_1\left[{\Psi }_{c,s}\right].$ В результате получается система обыкновенных дифференциальных уравнений:

$\begin{array}{c}{\partial }_z{\tilde{A}}_{c,s}=\alpha {\tilde{\Phi }}_{c,s}+\beta {\tilde{\Psi }}_{c,s}+{\tilde{a}}_{c,s},{\partial }_z{\tilde{\Phi }}_{c,s}=\lambda {\tilde{A}}_{c,s},\\ {\partial }_z{\tilde{B}}_{c,s}=\gamma {\tilde{\Phi }}_{c,s}+{\delta }_k{\tilde{\Psi }}_{c,s}-{\tilde{b}}_{c,s},{\partial }_z{\tilde{\Psi }}_{c,s}=i{k}_0{\tilde{B}}_{c,s},\end{array}$ (5)

где: $\alpha =i{k}_0{\epsilon }_{\perp };$ $\beta =k_{0}g;$ $\gamma =-k_{0}g;$ ${\delta }_k=i{k}_0({\epsilon }_{\perp }-k^2\text{\hspace{0.05em}}/k_0^2);$ $\lambda =i{k}_0[1-k^2/(k_0^2{\epsilon }_{\parallel }\left)\right];$ ${\tilde{a}}_c={\tilde{b}}_s=S_0\delta (z-h);$ ${\tilde{a}}_s={\tilde{b}}_c=0.$      

Ввиду того, что в неоднородности системы (5) присутствует множитель $\delta (z-h),$ ее решение сводится к решению соответствующей однородной системы с условиями сшивки при $z=h.$ Эту однородную систему можно представить наглядно в виде

${\partial }_{z}u=Sv,{\partial }_{z}v=Tu,$ (6)

используя двумерные векторы и матрицы $2\times 2:$

$\begin{array}{c}u=\left(\begin{array}{c}{\tilde{A}}_{c,s}\\ {\tilde{B}}_{c,s}\end{array}\right),v=\left(\begin{array}{c}{\tilde{\Phi }}_{c,s}\\ {\tilde{\Psi }}_{c,s}\end{array}\right),\\ S=\left(\begin{array}{cc}\alpha & \beta \\ \gamma & {\delta }_k\end{array}\right),T=\left(\begin{array}{cc}\lambda & 0\\ 0& i{k}_0\end{array}\right).\end{array}$

Условие сшивки решений однородной системы (6) при $z=h$ получаются путем интегрирования уравнений (5) по малой окрестности источника, т.е. от $z=h-0$ до $z=h+0.$

При решении однородной системы (6) полезной оказывается матрица $Y\left(z\right),$ преобразующая вектор $v\left(z\right)$ в вектор $u\left(z\right):$ $u\left(z\right)=Y\left(z\right)v\left(z\right).$ Матрица $Y\left(z\right)$ аналогична матрице адмиттанса, связывающей электрическое и магнитное поля и облегчающей решение уравнений Максвелла при наличии резко растущих или затухающих волновых мод. Матрица $Y\left(z\right)$ удовлетворяет нелинейному дифференциальному уравнению типа Риккати ${\partial }_{z}Y=S-YTY.$

Краевыми условиями при $z\to \pm \infty$ для систем (5) и (6) являются условия затухания возмущений. При этом предельные значения матрицы $Y\left(z\right)$ при $z\to \pm \infty$ однозначно определяются путем комбинации затухающих решений системы (6). Это позволяет найти $Y\left(z\right)$ сверху и снизу от уровня источника, решая численно соответствующие задачи Коши для уравнения Риккати в направлении из $+\infty$ или $-\infty$ к источнику. При этом на уровне источника $z=h$ получаются различные матрицы $Y(h+0)$ и $Y(h-0).$ Сопоставление скачка матрицы $Y\left(z\right)$ с условием сшивки при $z=h$ позволяет найти векторы $u(h\pm 0)$ и $v(h\pm 0),$ т.е. начальные данные для численного решения задачи Коши для системы (6). Более подробное обсуждение краевых условий и сведения краевой задачи к задаче Коши для аналогичной системы уравнений можно найти в работах [Fedorov et al., 2020; 2021].

Решая численно задачу Коши для системы (6) с краевыми условиями на уровне источника вверх и вниз, можно найти высотные профили модифицированных потенциалов ${\tilde{A}}_{c,s}(z,k),$ ${\tilde{B}}_{c,s}(z,k),$ ${\tilde{\Phi }}_{c,s}(z,k)$ и ${\tilde{\Psi }}_{c,s}(z,k).$ После этого с помощью обратного преобразования Ганкеля находятся горизонтальные пространственные распределения компонент поля. Наконец, использование формул (2), записанных в цилиндрических координатах, дает решение системы (4).

Для нахождения электромагнитного поля тока с конечной длиной L использованы результаты, полученные для токового диполя. С этой целью большое число N диполей с токовым моментом $M_0$ размещаются в серединах разбиения линии тока на N равных частей, и создаваемые ими поля суммируются. Этой системе диполей соответствует величина тока $J_0=N{M}_0{L}^{-1}.$

3. МОДЕЛИРОВАНИЕ ЭЛЕКТРОМАГНИТНОГО ОТКЛИКА НА СЕЙСМОГЕННЫЙ ТОК

Возбуждаемые подземным источником электрические и магнитные компоненты поля рассчитаны для земной поверхности (z = 0) и для фиксированной высоты в верхней ионосфере (z = 500 км), соответствующей типичной высоте низкоорбитальных спутников. Расчеты проведены для источника с частотой f = 1 Гц и с разными масштабами $(L=1-100$ км) на глубине h = 10 км. Проводимость Земли выбрана ${\sigma }_g={10}^{-3}$ См/м, а диэлектрическая проницаемость $\mathrm{Re}{\epsilon }_g=10.$ Рассчитанные поля соответствуют току источника $J_0=1$ A.

Параметры ионосферной модели IRI были выбраны соответствующими зимним ночным условиям (LT = 23, 08.12.2021 г.) на широте Камчатки (Петропавловcк-Камчатский 53.0° с.ш., 158.6° в.д.). При этих условиях foF2 = 2.37 МГц, а высота максимума плотности h_mF2 = 344 км. Атмосферная проводимость у поверхности Земли принимается равной ${\sigma }_a=1.1\cdot {10}^{-14}$ См/м.

3.1. Электромагнитные поля на поверхности земли

На земной поверхности УНЧ поля измеряются преимущественно магнитометрами, поэтому мы приводим результаты расчетов горизонтальных магнитных компонент $B_x$ и $B_y,$ и вертикальной компоненты B_z. Измерения электротеллурического поля (E_x- и E_y-компоненты) также ведутся, хотя и не столь широко. Возможные УНЧ электромагнитные предвестники были обнаружены на сравнительно малых удалениях от эпицентра [Hattori, 2004], поэтому расчеты проводились до расстояний по горизонтали не более 200 км.

3.1.1. Сопоставление с моделью подземного диполя

Для оценки гипотетического тока в области гипоцентра землетрясения, необходимого для возбуждения наблюдаемых на земной поверхности электромагнитных сигналов, авторы работы [Bortnik et al., 2010] использовали относительно простую модель подземного источника тока. В качестве источника рассматривался горизонтальный диполь c моментом M₀ = 1 А ⋅ м на глубине 10 км. Для расчета были использованы известные аналитические решения задачи о поле заглубленного горизонтального диполя [King, 1981; Banos, 1966]. Результаты расчетов для проводимостей ${\sigma }_g<{10}^{-3}$ См/м показали, что для наблюдаемых импульсов с амплитудами ∼30 нТл на частоте 1 Гц [Bleir et al., 2009] ожидаемая величина сейсмогенного тока $J_0$ попадает в диапазон от нескольких десятков до нескольких сотен кА при $L=0.1-1$ км. Если установить порог регистрации на уровне 1 пТл, то согласно этим оценкам даже при больших проводимостях грунта магнитные сигналы должны легко обнаруживаться в радиусе до 30 км от эпицентра.

Для сопоставления результатов нашей модели с результатами работы [Bortnik et al., 2010] на рис. 2 приведен расчет пространственного распределения компоненты магнитного поля $B_y\left(x\right)$ при $y=0$ на поверхности Земли для следующих параметров: $h=-10$ км, ${\sigma }_g={10}^{-3}$ См/м, $L=1$ км. Компонента $B_x\left(x\right)$ в этой области мала и не приводится. Для этих параметров глубина скин-слоя составляет ${\delta }_g=\sqrt{2/{\mu }_0\omega {\sigma }_g}$ ≈ 16 км. Выбранный малый масштаб источника должен дать результаты, соответствующие приближению точечного диполя. В реальных условиях длина разлома, а значит и масштаб гипотетического сейсмогенного тока, может значительно превышать 1 км. Модель точечного диполя не позволяет оценить изменение наземного возмущения в зависимости от масштаба источника. Такую оценку можно сделать с помощью разработанной нами модели источника конечной длины.

 

Рис.  2. Пространственная структура амплитуды магнитной компоненты $\mathbf{\left|}{\mathit{B}}_{\mathbf{y}}\mathbf{\right(}\mathit{x}\mathbf{\left)}\mathbf{\right|}$ на земной поверхности в направлении вдоль тока для источников с $\mathit{L}\mathbf{=}\mathbf{1}$ км и ${\mathit{J}}_{\mathbf{0}}\mathbf{=}\mathbf{50}$ A (соответствует точечному диполю) и c $\mathit{L}\mathbf{=}\mathbf{50}$ км и ${\mathit{J}}_{\mathbf{0}}\mathbf{=}\mathbf{1}$A. Частота источников $\mathit{f}\mathbf{=}\mathbf{1}$ Гц, удельная проводимость Земли σ_g = 10^–3 См/м.

 

С этой целью была рассчитана пространственная структура возбуждаемого поля для $L=50$ км и для L = 1 км, но при тех же параметрах среды. Сила тока подбиралась так, чтобы токовый момент $M_0$ обоих источников был одинаковым (рис. 2). Сопоставление показывает, что поле источника конечных размеров на больших удалениях спадает медленнее, чем поле точечного диполя. Непосредственно над источником с масштабом L = 50 км величина магнитного возмущения на земной поверхности достигает ∼10 пТл.

3.2. Особенности поля подземного излучателя

Ключевая проблема при поиске электромагнитных предвестников – возможность дискриминации магнитосферно-ионосферных и сейсмогенных возмущений. Если не применять специальные критерии для разделения магнитосферно-ионосферных и сейсмогенных возмущений (ряд которых приведен ниже), то стандартный подход может привести к тому, что магнитосферные сигналы будут ложно интерпретированы как сейсмогенные излучения [Костерин и др., 2015].

3.2.1. Поляризационная структура

Одно из отличий поля магнитосферно-ионосферного и подземного источников видно на рис. 3, где показано пространственное распределение горизонтальной и вертикальной магнитных компонент в направлении поперек токового источника (по оси y). Горизонтальная компонента $\left|\text{\hspace{0.05em}}{B}_y\right(y\left)\right|$ имеет максимум при $y=\mathrm{0,}$ т.е. непосредственно над источником. Величина же вертикальной компоненты $\left|B_z\right(y\left)\right|$ над источником проходит через нулевое значение и достигает максимума на удалении ∼15 км.

Проведенные расчеты показывают характерную особенность поля подземного источника по сравнению с полем падающих сверху магнитосферно-ионосферных возмущений. При типичных масштабах такого возмущения и проводимости подстилающей земной поверхности вертикальная компонента поля возмущения B_z мала по сравнению с горизонтальной $B_{\perp }:$ $\left|B_z\right|/\left|B_{\perp }\right|~{\delta }_g/{\lambda }_{\perp }<<\mathrm{1,}$ где ${\lambda }_{\perp }$ – горизонтальный масштаб возмущения [Pilipenko et al., 1998]. В то же время структура поля подземного источника такова, что $\left|B_z\right|>\left|B_y\right|$ вблизи источника, а вдали от него $\left|B_z\right|~\left|B_y\right|$ (рис. 3). Эту особенность использовали при выделении сейсмогенных сигналов перед сильным землетрясением на о. Гуам [Hayakawa et al., 1996]. Далее мы укажем на еще одну особенность, которая может быть использована для идентификации сигналов от подземных источников.

 

Рис.  3. Пространственная структура амплитуды магнитных компонент $\mathbf{\left|}{\mathit{B}}_{\mathbf{y}}\mathbf{\right(}\mathit{y}\mathbf{\left)}\mathbf{\right|}$ и $\mathbf{\left|}{\mathit{B}}_{\mathbf{z}}\mathbf{\right(}\mathit{y}\mathbf{\left)}\mathbf{\right|}$ на земной поверхности в направлении поперек тока для источника c $\mathit{L}\mathbf{=}\mathbf{50}$ км и ${\mathit{J}}_{\mathbf{0}}\mathbf{=}\mathbf{1}$ A на частоте $\mathit{f}\mathbf{=}\mathbf{1}$ Гц.

 

3.2.2. Кажущийся импеданс

Для дискриминации магнитосферных и сейсмогенных возмущений может быть применен кажущийся импеданс, т.е. отношение ${\mu }_0\left|E\right|/\left|B\right|$ [Pilipenko et al., 2005]. Для магнитосферных источников кажущийся импеданс практически для всех реальных условий совпадает с поверхностным импедансом Земли $Z_g\left(\omega \right)=\sqrt{-i\omega {\mu }_0/{\sigma }_g}.$ Кажущийся импеданс сейсмогенных возмущений был рассчитан по модельным значениям магнитного и электро-теллурического полей $Z_{xy}={\mu }_0{E}_x\text{/}{B}_y$ и $Z_{yx}={\mu }_0{E}_y\text{/}{B}_x$ (рис. 4). Расчет показывает, что для подземного источника кажущийся импеданс вблизи источника на порядок может превышать импеданс Земли $Z_g.$ Согласно этим расчетам, синхронные электротеллурические и геомагнитные наблюдения, дающие возможность определять импеданс зарегистрированных возмущений, могут быть применены для выделения возмущений от сейсмогенных источников.

 

Рис.  4. Профиль кажущихся импедансов ${\mathit{Z}}_{\mathbf{x}\mathbf{y}}\mathbf{\left(}\mathit{x}\mathbf{\right)}$ и ${\mathit{Z}}_{\mathbf{x}\mathbf{y}}\mathbf{\left(}\mathit{y}\mathbf{\right)}$ в сравнении с импедансом Земли ${\mathit{Z}}_{\mathbf{g}}$ (штриховая линия) на различных удалениях от источника. Расчет проведен для значений h = –10 км, ${\mathit{\sigma }}_{\mathbf{g}}\mathbf{=}{\mathbf{10}}^{\mathbf{-}\mathbf{3}}$ См/м, L = 50 км, f = 1 Гц.

 

3.3. Электромагнитные поля в верхней ионосфере

Спутниковые измерения УНЧ и КНЧ полей преимущественно ведутся электрическими антеннами, которые более чувствительны к электромагнитным возмущениям, чем бортовые магнитометры. Для интерпретации спутниковых наблюдений сигналов, предположительно вызванных сейсмическими источниками, проведены расчеты амплитуды поперечных электрических составляющих $\left|E_x\right|$ и $\left|E_y\right|$ на высоте $z=500$ км (рис. 5). При $f=1$ Гц максимальное возмущение непосредственно над источником $(x=0)$ с J₀ = 1 A и L = 50 км достигает $\left|E_x\right|~3\cdot {10}^{-3}$ мкВ/м.

 

Рис.  5. Горизонтальная пространственная структура амплитуды горизонтальных электрических компонент $\mathbf{\left|}{\mathit{E}}_{\mathbf{x}}\mathbf{\text{}}\mathbf{\right(}\mathit{x}\mathbf{\left)}\mathbf{\right|}$ и $\mathbf{\left|}{\mathit{E}}_{\mathbf{y}}\mathbf{\text{}}\mathbf{\right(}\mathit{x}\mathbf{\left)}\mathbf{\right|}$ в верхней ионосфере на высоте $\mathit{z}\mathbf{=}\mathbf{500}$км излучения подземного источника c $\mathit{L}\mathbf{=}\mathbf{50}$ км и ${\mathit{J}}_{\mathbf{0}}\mathbf{=}\mathbf{1}$ A на частоте $\mathit{f}\mathbf{=}\mathbf{1}$ Гц. Источник ориентирован либо по оси x (ϕ = 0°, сплошная линия), либо по оси y (ϕ = 90°, штриховая линия).

 

В ранних публикациях сообщалось об обнаружении квазистационарных (∼1 Гц) электромагнитных возмущений с характерными амплитудами порядка первых мВ/м по электрическому полю [Chmyrev et al., 1989]. Из результатов модели следует, что для создания в верхней ночной ионосфере возмущения с амплитудой Е∼3 мВ/м необходим крупномасштабный подземный ток на глубине 10 км с интенсивностью $J_0~{10}^6$ А. При этом на земной поверхности при тех же параметрах источника возникнет возмущение геомагнитного поля  нТ (рис. 2).

4. ОБСУЖДЕНИЕ

Автор статьи [Молчанов, 1991] оценивал просачивание поля подземного источника в верхнюю ионосферу с помощью многослойной модели, включающей низкопроводящий слой в области очага землетрясения, высокопроводящий поверхностный слой земной коры, изотропно проводящую атмосферу, анизотропно проводящую ионосферу и магнитосферу. На границе каждого слоя распространяющиеся волны с частотами 0.01–100 Гц испытывали частичное отражение. Чтобы свести задачу к одномерной, геометрия предполагалась азимутально-симметричной. Также численные модели проникновения УНЧ полей от подземного источника в ионосферу предлагались в работах [Tian, Hata, 1996; Dong et al., 2005; Wang et al., 2021]. Однако все эти модели имеют те или иные существенные ограничения и не могут дать исчерпывающее решение проблемы. Провести прямое сопоставление их расчетов с нашей моделью не представляется возможным из-за существенной разницы в геометрии источников.

В представленной работе кратко описан теоретический формализм для расчета электромагнитных полей в системе Земля–атмосфера–ионосфера, создаваемых подземным токовым источником конечной длины. На основании этой теории построена численная модель, дающая возможность оценить проникновение в реалистичную ионосферу УНЧ излучения подземного источника, который моделируется как горизонтальный ток конечного масштаба. Важное достоинство модели состоит в возможности согласованной оценки возмущений в приземной атмосфере и в верхней ионосфере при заданном источнике. В представленной модели геомагнитное поле считалось вертикальным. Однако учет наклона геомагнитного поля существенно не изменяет результаты [Fedorov et al., 2021]. Геомагнитное поле может частично канализировать электромагнитное УНЧ–КНЧ излучение, поэтому на средних широтах должно наблюдаться смещение к экватору максимального отклика в ионосфере на сейсмогенный ток.

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

Из полученных на космических аппаратах результатов по регистрации предвестников землетрясений наиболее цитируемыми являются данные спутников ОРЕОЛ-3 и ИКБ-1300 [Чмырев и др., 1986; Chmyrev et al., 1989]. На обоих спутниках эффекты были выделены в ночное время, когда геомагнитные условия были умеренно возмущенными. При пролете ИКБ-1300 на высоте 800 км над очагом землетрясения за 15 мин до главного удара были зарегистрированы вариации магнитного и электрического поля в УНЧ диапазоне 0.1–8 Гц с амплитудами порядка первых единиц нТл и мВ/м. Спутник ОРЕОЛ-3 зарегистрировал подобные УНЧ возмущения над тем же эпицентром землетрясения на высоте 1970 км за ∼5 ч до главного удара. Наличие последовательных измерений с двух спутников позволило авторам заключить, что сейсмоэлектромагнитные шумы длительно присутствовали в области над эпицентром готовящегося землетрясения.

Много интересных результатов дал специализированный спутниковый проект DEMETER [Parrot, 2018]. Чувствительность электромагнитного комплекса аппаратуры была настолько велика, что на авроральных широтах датчики переходили в насыщение. Перед сильным землетрясением с М = 7.9 на глубине 10 км было обнаружено усиление электрической компоненты шумов в окрестности 5.8 Гц [Walker et al., 2013]. Измерения электрической компоненты показали увеличение в среднем энергии УНЧ (<20 Гц) флуктуаций до E²∼1–20 (мВ/м)² над районами с высокой сейсмичностью (M > 5.0), особенно над границами тектонических плит [Athanasiou et al., 2014]. Статистика за несколько лет регистрации электрической компоненты показала небольшое, но статистически значимое, увеличение числа УНЧ сигналов в частотном диапазоне 0.05–1 Гц над эпицентрами будущих землетрясений с М > 5.0 и H < 70 км в среднем за неделю до толчка [Ouyang et al., 2020].

Интригующие результаты поступают со специализированного прогностического спутникового проекта CSES [Zhu et al., 2021]. Измерения электромагнитных шумов в полосе частот 75–90 Гц показали усиление мощности шумов на 10–30% за несколько дней до землетрясений с М = 6.4 и М = 7.4 [Wang et al., 2022]. Было обнаружено увеличение УНЧ флуктуаций электрической компоненты над уровнем дисперсии 0.2 мВ/м перед землетрясением с М = 7.7 и Н = 10 км [Li et al., 2022].

Авторы всех приведенных результатов неявно предполагали, что зарегистрированные на спутниках в верхней ионосфере УНЧ возмущения с типичными амплитудами порядка единиц мВ/м генерируются крупномасштабными сейсмогенными источниками. Разработанная численная модель предсказывает, что электрическая составляющая УНЧ излучения с частотой 1 Гц может достигать таких амплитуд в верхней ночной ионосфере над горизонтальным источником с силой тока ${10}^6$ А и длиной 50 км на глубине 10 км в низкоомной земной коре $(\rho ={10}^3$ Ом ⋅ м). Результат моделирования соответствует наиболее благоприятным условиям, когда спутник находится точно над источником на небольшой глубине. При этом ионосферный профиль соответствует ночным условиям, когда поглощение УНЧ излучения в ионосферной плазме мало по сравнению с дневными условиями [Fedorov et al., 2020].

По представлениям, развиваемым в работах [Финкель, 1977; Гохберг и др., 1985], во время магистрального разрушения может возникать упорядоченное движение заряженных дислокаций. Плотность дислокационного тока может варьироваться в пределах $j~{10}^{-3}\text{}-{10}^3$ А/м². Поперечный размер крупномасштабных разрывов составляет $10-{10}^3$ м. Если принять $j~1-10$ А/м², то полный ток может составлять $J_0~{10}^2-{10}^7$ А. Если такой всплеск тока величиной $~{10}^6$ А и с временным масштабом 1 с в очаге готовящегося землетрясения возможен, то в принципе он может вызвать возмущение Е-поля в верхней ионосфере с амплитудой около единиц мВ/м, доступное для измерения антеннами низкоорбитальных спутников. Однако согласно численной модели при этом на земной поверхности возникнет возмущение геомагнитного поля $B~{10}^4$ нТ! Такие магнитные возмущения, сопоставимые с возмущениями при сильных суб-бурях на авроральных широтах, были бы обнаружены существующей достаточно плотной сетью магнитометров. Поэтому зарегистрированные ионосферные УНЧ возмущения перед землетрясениями в спутниковых экспериментах вряд ли можно связать с прямым излучением источников сейсмической природы.

Современные низкоорбитальные космические аппараты позволяют при низкой фоновой активности регистрировать на частотах порядка первых единиц – десятков Гц слабые излучения с амплитудами >1 мкВ/м [Simões et al., 2012], поэтому спутниковая регистрация сейсмогенных УНЧ–КНЧ излучений возможна с аппаратурой на 3 порядка лучшей чувствительностью. Поскольку интенсивность сейсмического источника и величина затухания при распространении через ионосферу известны плохо, то для тестирования модели и проверки возможности обнаружения сейсмических излучений на околоземной орбите можно использовать скоординированные спутниковые и наземные наблюдения 50/60 Гц излучений линий электропередач [Пилипенко и др., 2021].

Помимо количественных расчетов распространения возмущений в системе Земля–атмосфера–ионосфера, крайне важен вопрос о возможном физическом механизме подземных источников этих электромагнитных полей. Попытки построить более реалистичную модель физических процессов, протекающих в массивах горных пород при их деформировании, предпринимались в ряде работ. Генерация наземного электромагнитного отклика на линейно нарастающий во времени подземный ток вдоль сейсмического разлома была численно смоделирована в работе [Losseva, Nemchinov, 2005]. Генерируемый импульс распространяется в проводящей земле диффузионным образом $~\sqrt{Dt},$ где $D={\left({\mu }_0{\sigma }_g\right)}^{-1}$ – коэффициент магнитной диффузии. В работе [Лосева и др., 2012] построена численная модель генерации электрических импульсов системой диполей с переменным моментом, распределенных по поверхности релаксирующего структурного блока при его стесненном повороте. Результаты регистрации электрического поля в приповерхностном слое земной коры в зоне влияния тектонически активной структуры согласуются с моделью генерации электромагнитных импульсов при поляризации горных пород в результате дифференциальной подвижки структурного блока при его релаксации [Лосева и др., 2010]. Рассмотрение импульсных возмущений нами планируется провести в последующих работах.

5. ЗАКЛЮЧЕНИЕ

Разработан теоретический формализм для расчета электромагнитных полей в связанной системе Земля–атмосфера–ионосфера, создаваемых подземным токовым источником конечной длины. На основании этой теории создана численная модель с реалистичным профилем ионосферы, находящейся в вертикальном геомагнитном поле. С помощью этой численной модели рассчитана пространственная структура поля на земной поверхности в окрестности источника. Показано, что характерные особенности такого поля подземного источника (превышение на порядок кажущегося импеданса над импедансом Земли и конечное значение отношения $\left|B_z\right|/\left|B_{\perp }\right|)$ могут быть использованы для дискриминации возмущений от сейсмогенных источников. Разработанная модель позволяет рассчитать ожидаемую амплитуду УНЧ излучений на спутниковых высотах для гипотетического сейсмогенного тока. Результаты численного моделирования позволяют связать возмущения приземного магнитного поля и электрического поля в ионосфере, создаваемые крупномасштабным подземным источником. Хотя сильные токи в области очага на заключительной фазе подготовки землетрясения представляются возможными, модельные расчеты показывают, что найденные в спутниковых данных возмущения электрического УНЧ поля перед землетрясениями нельзя связать с прямым излучением сейсмогенных источников, т.к. они должны сопровождаться нереально интенсивными возмущениями геомагнитного поля на поверхности Земли.

БЛАГОДАРНОСТИ

В работе использована модель IRI, распространяемая через сайт НАСА (https://irimodel.org). Выражаем признательность рецензенту за конструктивные замечания.

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа поддержана грантом РНФ 22-17-00125.

Об авторах

Н. Г. Мазур

Институт физики Земли им. О.Ю. Шмидта РАН

Автор, ответственный за переписку.
Email: ngmazur@mail.ru
Россия, Москва

Е. Н. Федоров

Институт физики Земли им. О.Ю. Шмидта РАН

Email: enfedorov1@yandex.ru
Россия, Москва

В. А. Пилипенко

Институт физики Земли им. О.Ю. Шмидта РАН

Email: space.soliton@gmail.com
Россия, Москва

К. Е. Боровлева

Институт физики Земли им. О.Ю. Шмидта РАН

Email: k.borovlyova@gmail.com
Россия, Москва

Список литературы

Дополнительные файлы