Stability of heavy inclusions in the Earth’s crust

Full Text

Введение

Линейная теория устойчивости прослеживает эволюцию начальных возмущений той или иной пространственной формы. Согласно линейной теории, система устойчива, если она устойчива к малым возмущениям любого типа, и неустойчива, если она неустойчива хотя бы к одному типу малых возмущений. В геодинамике при рассмотрении медленных процессов, связанных с гравитационной неустойчивостью, обычно рассматриваются только малые начальные возмущения, а развитие возмущения представляется в виде $a\cdot e^{\lambda t},$ где a – начальное возмущение, а λ – комплексный инкремент. Если действительная часть λ отрицательна или равна нулю, имеет место устойчивость, а если положительна – неустойчивость. Когда начальное возмущение мало, а $\mathrm{Re}\lambda \le 0$ (устойчивость) линейная теория полностью описывает эволюцию начального малого возмущения. Когда $\mathrm{Re}\lambda >0$ (неустойчивость), линейная теория описывает эволюцию возмущения только на не слишком больших временах, при которых возмущение остается малым.

Реологию коры и мантии можно описать вязкоупругой моделью Максвелла, хотя эта модель и дает упрощенное представление о ползучести геоматериала [Биргер, 2016; Karato, 2008]. Согласно модели Максвелла, среда ведет себя как упругая, когда $\frac{\eta }{\mu }\left|\lambda \right|\gg \mathrm{1,}$ где µ – упругий модуль сдвига, а η – коэффициент вязкости, и как вязкая среда, когда $\frac{\eta }{\mu }\left|\lambda \right|\ll 1.$ Вязкость геоматериала сильно зависит от температуры, которая в Земле быстро растет с глубиной. Поэтому верхняя кора обладает очень высокой вязкостью и ведет себя как упругая среда даже при медленных процессах, а более глубокие слои Земли, где вязкость значительно ниже, ведут себя при медленных процессах как вязкая среда. В работе [Биргер, 2023] была исследована устойчивость тяжелого включения в верхней упругой коре Земли, подстилаемой вязким слоем, который представляет нижнюю кору и мантийную литосферу. Такая модель справедлива только в случае медленных безынерционных течений, при которых модуль комплексного инкремента $\left|\lambda \right|$ мал. Анализ устойчивости, при котором найдены только малые значения $\left|\lambda \right|,$ не является полным. В настоящей работе ищутся и большие значения $\left|\lambda \right|,$ при поиске которых надо учитывать инерционность и упругость среды, но можно пренебречь ее вязкостью.

Анализ устойчивости упругих слоев

Рассматривается двухслойная модель. Нижний слой представлен как полупространство. Начало координат помещено на нижней поверхности верхнего слоя, а ось z направлена вертикально вверх. Слой (0 < z < d) моделирует верхнюю упругую кору, а полупространство (z < 0) – подстилающую литосферу и мантию. Уравнения, описывающие возмущения литостатического равновесия несжимаемой упругой среды, записываются в виде:

$-\frac{\partial p}{\partial x}+\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{xz}}{\partial z}=\rho \frac{{\partial }^2}{\partial t^2}{u}_x,$     (1)

$-\frac{\partial p}{\partial z}+\frac{\partial {\sigma }_{xz}}{\partial x}+\frac{\partial {\sigma }_{zz}}{\partial z}=\rho \frac{{\partial }^2}{\partial t^2}{u}_z,$     (2)

$\frac{\partial u_x}{\partial x}+\frac{\partial u_z}{\partial z}=\mathrm{0,}$     (3)

где: p – возмущение давления; ${\sigma }_{xz}$ и ${\sigma }_{zz}$ – компоненты девиаторного тензора напряжений; $\rho$ – плотность; $u_x$ и $u_z$ – смещения; t – время. Уравнения (1) и (2) описывают двумерное движение, а уравнение (3) представляет собой условие несжимаемости. К этим уравнениям добавляется уравнение состояния для упругой среды:

${\sigma }_{ij}=2\mu {\epsilon }_{ij},$     (4)

связывающее девиаторные напряжения с деформациями ${\epsilon }_{ij},$ которые определены как

${\epsilon }_{xx}=\frac{\partial u_x}{\partial x},{\epsilon }_{zz}=\frac{\partial u_z}{\partial z},{\epsilon }_{xz}=\frac{1}{2}\left(\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x}\right).$     (5)

Сжимаемостью упругой среды можно пренебречь при условии

$\frac{\mu }{K}\ll \mathrm{1,}$ $\nu =\frac{3K-2\mu }{2\left(3K+\mu \right)}=1/\mathrm{2,}$     (6)

где K – модуль всестороннего сжатия; ν – коэффициент Пуассона. Для геоматериала ν ≈ 0.3, $\frac{\mu }{K}\approx 0.5$ и условие несжимаемости (6), строго говоря, не выполняется. Однако учет сжимаемости, как показали проведенные автором расчеты, качественно не изменяет результаты, полученные для несжимаемой упругой среды и представленные в настоящей статье. В уравнения, учитывающие сжимаемость среды, модуль всестороннего сжатия входит в виде безразмерного параметра $\frac{\mu }{3K},$ который можно считать малым.

Упругий модуль сдвига и плотность недр можно оценить как $\mu \approx 6\cdot {10}^{10}Па,$ $\rho \approx 3\cdot {10}^3кг/{м}^3.$ Все дальнейшие соотношения будут выписаны для безразмерных физических переменных, для которых сохраняем те же обозначения, что и для размерных. В качестве масштаба длины использована толщина верхнего упругого слоя коры, которая оценивается как d = 10 км, масштаб напряжения – упругий модуль сдвига µ, масштаб скорости – $\sqrt{\frac{\mu }{\rho }},$ масштаб времени – $d\sqrt{\frac{\rho }{\mu }}.$ Согласно приведенным оценкам, получаем масштаб времени $d\sqrt{\frac{\rho }{\mu }}\approx 2с.$ Нет смысла рассматривать возмущения, длина волны которых $\frac{2\pi }{k}$ превышает 1000 км. Такое ограничение, наложенное на длину волны, приводит к ограничению $k>6\cdot {10}^{-2},$ наложенному на безразмерное волновое число.

Вертикальное смещение представим в виде:

$u_z=U_z\left(z\right)\exp \left(i\omega t\right)\cos \left(kx\right),$ $k>\mathrm{0,}$     (7)

где  – комплексная частота, k – действительное волновое число. Комплексный инкремент λ, используемый в линейной теории устойчивости, записан в виде $\lambda =i\omega .$ В аналогичном виде представим и все остальные физические переменные, причем в выражения для $u_x,$ ${\epsilon }_{xz}$ и ${\sigma }_{xz}$ входит не $\cos \left(kx\right),$ а $\sin \left(kx\right).$ Такое представление позволяет свести систему уравнений в частных производных (1)–(5) к системе обыкновенных дифференциальных уравнений, в которой все переменные, характеризующие смещения, деформации, напряжения и давление, зависят только от вертикальной координаты z. Рассматриваемая механическая система неустойчива, если действительная часть инкремента положительна, т.е. если $\mathrm{Im}\omega <0.$

Волновое число определяется горизонтальным размером L возмущения плотности (тяжелого включения) в верхней коре. Вертикальный размер d этой аномальной области выбран в качестве масштаба длины и, следовательно, равен 1. Такое возмущение плотности создает возмущение гравитационной силы, которое вызывает течение. Возможность перехода от рассмотрения возмущения с горизонтальным размером L к периодическому по горизонтали возмущению связана с тем, что первая гармоника $\cos \left(\frac{\pi }{L}x\right)$ дает основной вклад в разложение Фурье возмущения с горизонтальным размером L [Cathles, 1975; Биргер, 2016; 2018]. Таким образом, безразмерное волновое число k, характеризующее периодичность по горизонтали, связано с горизонтальным размером аномальной области как $k=\frac{\pi }{L}d$ (полудлина волны равна L). При периодическом движении среда, неограниченная по горизонтали, разбивается на бесконечный набор вертикальных столбов. Движение в каждом из них моделирует движение, которое возникает в случае, когда возмущение плотности происходит в области с горизонтальном размером L.

Подставляя (7) в уравнения (1)–(5), выражаем все физические переменные через вертикальное смещение:

$u_x=\frac{1}{k}D{U}_z\exp \left(i\omega t\right)\sin \left(kx\right),$    (8)

${\sigma }_{zz}=$$-{\sigma }_{xx}=2D{U}_z\exp \left(i\omega t\right)\cos \left(kx\right),$    (9)

${\sigma }_{xz}=\frac{1}{k}\left(D^2+k^2\right)U_z \exp \left(i\omega t\right)\sin \left(kx\right),$     (10)

$p=\frac{D}{k^2}\left(L^2+{\omega }^2\right)U_z\exp \left(i\omega t\right)\cos \left(kx\right),$     (11)

а для амплитуды вертикального смещения в верхнем упругом слое получаем обыкновенное дифференциальное уравнение

$\left(L^4+{\omega }^2{L}^2\right)U_z=0.$     (12)

В уравнениях (8)–(12) введены дифференциальные операторы $D=\frac{d}{dz}$ и $L^2=D^2-k^2.$ Вертикальное смещение в упругом полупространстве описывается уравнением:

$\left(L^4+{\omega }^2\left(1-\alpha \right)L^2\right)U_z=\mathrm{0,}$     (13)

где введен безразмерный параметр $\alpha =\frac{\Delta \rho }{\rho },$ в котором $\Delta \rho$ – скачок плотности на границе между верхним тяжелым слоем и полупространством. Соотношения (8)–(10) справедливы и для полупространства, а соотношение (11) принимает вид:

$p=\frac{D}{k^2}\left(L^2+\left(1-\alpha \right){\omega }^2\right)U_z\exp \left(i\omega t\right)\cos \left(kx\right),$     (14)

Уравнение (12) имеет общее решение:

$U_z=e^{kz}{C}_1+e^{-kz}{C}_2+e^{z\sqrt{k^2-{\omega }^2}}{C}_3+e^{-z\sqrt{k^2-{\omega }^2}}{C}_4,$

$0\le z\le 1.$     (15)

Общее решение уравнения (13) записывается как

$U_z=e^{kz}{B}_1+e^{z\sqrt{k^2-\left(1-\alpha \right){\omega }^2}}{B}_2,$ $-\infty \le z\le \mathrm{0,}$     (16)

где $C_i$ и $B_i$ (i = 1, 2, 3, 4) – произвольные константы. Уравнение (16) следует из требования ограниченности смещений в подстилающем слой полупространстве. Согласно (16), смещения в полупространстве убывают с глубиной как $e^{kz}.$ Полупространство моделирует нижний слой, в который смещения проникают на глубину, зависящую от волнового числа k.

На верхней деформированной поверхности обращается в нуль сила, действующая на единицу площади поверхности. Это условие можно перенести на исходную плоскую поверхность $(z=1)$ как [Биргер, 2023]:

${\sigma }_{xz}=\mathrm{0,}$     (17)

$-p+{\sigma }_{zz}+\phi u_z=0.$     (18)

В уравнении (18) введен безразмерный параметр $\phi =\frac{\rho gd}{\mu }\approx \mathrm{0.005,}$ где g – ускорение силы тяжести, ρ – плотность верхнего слоя. Учитывая соотношения (8)–(11), уравнения (17) и (18) перепишем в виде:

$\left(D^2+k^2\right)U_z=\mathrm{0,}$     (19)

$\left(3D-\frac{D^3}{k^2}+\phi \right)U_z=0.$     (20)

На деформируемой границе $(z=0)$ между верхним слоем и полупространством непрерывны смещения (вертикальные и горизонтальные) и касательные напряжения, а нормальные напряжения связаны соотношением [Биргер, 2023]:

$-p^{\left(2\right)}+{\sigma }_{zz}^{\left(2\right)}+\alpha \phi u_z=-p^{\left(1\right)}+{\sigma }_{zz}^{\left(1\right)},$     (21)

где индекс 1 относится к полупространству, а индекс 2 – к верхнему слою. Если считать, что плотность верхнего слоя выше, чем плотность подстилающего полупространства, перепад плотности $\Delta \rho$ и параметр $\alpha$ положительны. Поскольку в недрах Земли перепады плотности малы, можно считать, что $\alpha \ll 1.$ С учетом уравнений (9), (11) и (14) условия на границе (z = 0), разделяющей слои, записывается в виде:

$U_z^{\left(2\right)}=U_z^{\left(1\right)},$     (22)

$D{U}_z^{\left(2\right)}=D{U}_z^{\left(1\right)},$     (23)

$D^2{U}_z^{\left(2\right)}=D^2{U}_z^{\left(1\right)},$     (24)

$\left(-\frac{D^3}{k^2}+\alpha \phi \right)U_z^{\left(2\right)}=\left(-\frac{D^3}{k^2}-\alpha {\omega }^2\frac{D}{k^2}\right)U_z^{\left(1\right)}.$     (25)

Граничные условия (19), (20), (22)–(24) приводят к системе из шести однородных уравнений для шести произвольных констант, введенных в уравнениях (15) и (16),

$\begin{array}{c}2k^2{e}^k{C}_1+2k^2{e}^{-k}{C}_2+(2k^2-{\omega }^2)e^{\sqrt{k^2-{\omega }^2}}{C}_3+\\ +(2k^2-{\omega }^2)e^{-\sqrt{k^2-{\omega }^2}}{C}_4=\mathrm{0,}\end{array}$     (26)

$\begin{array}{c}\left(\phi +2k-\frac{{\omega }^2}{k}\right)e^k{C}_1+\left(\phi -2k+\frac{{\omega }^2}{k}\right)e^{-k}{C}_2+\\ +\left(\phi +2\sqrt{k^2-{\omega }^2}\right)e^{\sqrt{k^2-{\omega }^2}}{C}_3+\\ +\left(\phi -2\sqrt{k^2-{\omega }^2}\right)e^{-\sqrt{k^2-{\omega }^2}}{C}_4 =\mathrm{0,}\end{array}$     (27)

$C_1+C_2+C_3+C_4-B_1-B_2=\mathrm{0,}$     (28)

$\begin{array}{c}k{C}_1-k{C}_2+\sqrt{k^2-{\omega }^2}{C}_3-\sqrt{k^2-{\omega }^2}{C}_4-\\ -k{B}_1-\sqrt{k^2-\left(1-\alpha \right){\omega }^2}{B}_2=\mathrm{0,}\end{array}$     (29)

$\begin{array}{c}2k^2{C}_1+2k^2{C}_2+\left(2k^2-{\omega }^2\right)C_3+\\ +\left(2k^2-{\omega }^2\right)C_4-2k^2{B}_1-\left(2k^2-\left(1-\alpha \right){\omega }^2\right)B_2=\mathrm{0,}\end{array}$     (30)

$\begin{array}{c}\left(\alpha \phi -k\right)C_1+\left(\alpha \phi +k\right)C_2+\\ +\left(\alpha \phi -\frac{1}{k^2}(k^2-{\omega }^2)\sqrt{k^2-{\omega }^2}\right)C_3+\end{array}$      

$+\left(\alpha \phi +\frac{1}{k^2}(k^2-{\omega }^2)\sqrt{k^2-{\omega }^2}\right)C_4+\left(k-\alpha \frac{{\omega }^2}{k}\right)B_1+$       

$+\frac{1}{k^2}\left(k^2-{\omega }^2\right)\sqrt{k^2-\left(1-\alpha \right){\omega }^2}{B}_2=0.$     (31)

Приравнивая нулю определитель этой системы линейных уравнений (если определитель отличен от нуля, все константы равны нулю), получаем характеристическое уравнение, связывающее $\omega ,$ k и α.

На рис. 1 представлена зависимость $\omega$ от k, полученная при фиксированном значении скачка плотности $\alpha =\frac{\Delta \rho }{\rho }=\mathrm{0.01.}$ При любых k характеристическое уравнение имеет решение $\omega \left(k\right)=0$ и решение $\omega \left(k\right)\approx 0.95931k.$ Поскольку все найденные значения $\omega$ – действительные (инкремент $\lambda =i\omega$ – чисто мнимый), тонкий упругий слой, лежащий на упругом слое с меньшей плотностью, устойчив: малые начальные возмущения вызывают незатухающие колебательные движения, аналогичные стоячим волнам Рэлея. Слабое затухание этих колебаний определяется реологической моделью Ломнитца, которая описывает неустановившуюся ползучесть при высоких частотах [Биргер, 2016; Karato, 2008].

Рис. 1. Зависимость частот $\omega$ от волновых чисел k при скачке плотности $\alpha =0.01$  и пунктирной линией – при $\alpha =0.$

На рис. 2 представлена найденная из характеристического уравнения зависимость $\omega$ от $\alpha$ при фиксированном малом значении волнового числа $k=\mathrm{0.001.}$ Как видно из этого рисунка, при малых значениях $\alpha$частота $\omega$ очень слабо зависит от этого параметра (верхняя прямая имеет малый наклон).

Рис. 2. Зависимость частот $\omega$ от скачка плотности $\alpha$ для волнового числа $k=\mathrm{0.001.}$

Рисунок 3, который соответствует расчету, проведенному при $k=\mathrm{0.1,}$ показывает, что частота $\omega$ возрастает с ростом скачка плотности $\alpha .$

Рис. 3. Зависимость частот $\omega$ от скачка плотности $\alpha$ для волнового числа $k=\mathrm{0.1.}$

Геофизическим приложениям соответствуют очень малые значения безразмерных параметров: $\alpha =\frac{\Delta \rho }{\rho }\approx 0.01$ и $\phi =\frac{\rho gd}{\mu }\approx \mathrm{0.005.}$ Если положить $\alpha =\phi =\mathrm{0,}$ полученное характеристическое уравнение сводится к виду

$4k^3\sqrt{k^2-{\omega }^2}-{\left(2k^2-{\omega }^2\right)}^2=\mathrm{0,}$     (32)

которое имеет решения $\omega \left(k\right)=0$ и $\omega \left(k\right)\approx 0.95529k.$ Учет влияния гравитации  незначительно повышает частоты  Рисунок 1 и рис. 3 показывают поправки к частоте Рэлея, которые дает скачок плотности. При скачке плотности $\alpha =0.01$ частоты повышаются до значений $\omega \left(k\right)\approx 0.95931k.$ Уравнение (32) описывает волну Рэлея в несжимаемой среде. С учетом сжимаемости уравнение (32) принимает вид

$4k^2\sqrt{\left(k^2-{\omega }^2\right)\left(k^2-b{\omega }^2\right)}-{\left(2k^2-{\omega }^2\right)}^2=\mathrm{0,}$     (33)

где $b=\mu /\left(K+\frac{4}{3}\mu \right)\approx \mathrm{0.3.}$ Этот параметр можно считать малым, поскольку он слабо влияет на результат. Частоты $\omega ,$ найденные из уравнений (32) и (33), отличаются приблизительно на 3%.

Анализ устойчивости с учетом ползучести коры

Материал земной коры, как и любой поликристаллический материал, обладает не только упругостью, но и ползучестью, которая определяется имеющимися дефектами в идеально правильной кристаллической структуре. Упруго-ползучая реология хорошо описывается моделью Максвелла, которая представляет полную деформацию в виде суммы упругой деформации и вязкой деформации. Ползучесть геоматериала, конечно, не описывается вязкой ньютоновской моделью, однако можно ввести эффективную ньютоновскую вязкость, зависящую от характерной продолжительности рассматриваемого процесса. Материал, описываемый моделью Максвелла, ведет себя как вязкий при медленных процессах, характерное время которых значительно превышает время Максвелла, представляющее собой отношение коэффициента вязкости к упругому модулю сдвига. В задаче об устойчивости характерное время процесса можно определить как $\tau =1/\left|\lambda \right|.$ Поэтому среда ведет себя как вязкая при малых значениях инкремента. В процессе с большим значением инкремента характерное время значительно короче, чем время Максвелла, и среда ведет себя как упругая.

С помощью данных о послеледниковых поднятиях эффективный коэффициент вязкости земных недр оценивается как $\eta \approx 4\cdot {10}^{19}Па\cdot с$ [Cathles, 1975; Биргер, 2018]. Упругий модуль сдвига оценивается как $\mu \approx 6\cdot {10}^{10}Па.$ Согласно этим оценкам, время Максвелла для нижнего менее вязкого слоя составляет $\eta /\mu \approx 6\cdot {10}^8с\approx 20 лет.$ На рис. 1, рис. 2, рис. 3 безразмерные частоты ω лежат в интервале от 0.01 до 1. Таким частотам соответствуют размерные периоды колебаний $\frac{2\pi }{\omega }\cdot 2с,$ которые лежат в интервале от 120 до 12 000 с. Поскольку эти периоды значительно короче, чем время Максвелла, пренебрежение вязкостью среды при таких частотах является вполне оправданным. Таким образом, найденный в предыдущем разделе статьи инкремент $\lambda =i\omega \left(\alpha ,k\right)$ достаточно велик для того, чтобы соответствующее этому инкременту движение описывалось чисто упругой моделью. Полученное характеристическое уравнение при любых k и $\alpha$ имеет еще корень $\lambda =i\omega =\mathrm{0,}$ показанный на рис. 1. Этот корень не удовлетворяет исходному предположению об упругости нижнего слоя (частоте $\omega =0$ соответствует бесконечно большой период, значительно превышающий время Максвелла) и поэтому должен быть отброшен.

Удобно изменить масштаб времени, введенный для упругой среды, и принять в качестве масштаба время Максвелла, определив его как $\eta /\mu \approx 6\cdot {10}^8с\approx 20 лет.$ После введения такого масштаба времени можно сказать, что нижний слой ведет себя как вязкий, если безразмерный инкремент удовлетворяет условию:

$\left|\lambda \right|\ll 1.$     (34)

Если вязкость верхнего слоя в m раз выше, чем вязкость нижнего слоя, то верхний слой ведет себя как вязкий при условии:

$m\left|\lambda \right|\ll 1.$     (35)

Поскольку инкремент λ заранее неизвестен, надо сначала решить задачу устойчивости для случая, когда оба слоя вязкие, и если условие $m\left|\lambda \right|\ll 1$ не выполняется, надо решать задачу для случая, когда верхний слой является упругим. Задача об устойчивости упругого слоя, под которым находится вязкий слой с меньшей плотностью, решена в работе [Биргер, 2023], где показано, что такая система устойчива.

В предыдущем разделе статьи анализ устойчивости проведен для упругой среды, что предполагает достаточно большое (по модулю) значение искомого инкремента. При малом инкременте можно пренебречь инерционностью и упругостью среды, но необходимо учитывать ее ползучесть. Уравнения, используемые для анализа устойчивости в случае, когда искомый инкремент предполагается малым, отличаются от уравнений предыдущего раздела статьи. Соотношения (8)–(10) не изменяются, а в соотношениях (11)–(14) следует положить $\omega =0.$ Уравнение (12), которое при $\omega =0$ совпадает с уравнением (13), имеет общее решение:

$U_z=e^{kz}{C}_1+z{e}^{kz}{C}_2+e^{-kz}{C}_3+z{e}^{-kz}{C}_4,$ $0\le z\le 1.$     (36)

$U_z=e^{kz}{B}_1+z{e}^{kz}{B}_2,$ $-\infty \le z\le 0.$     (37)

Граничные условия (19)–(23) остаются прежними, а условия (24) и (25) принимают вид:

$m{D}^2{U}_z^{\left(2\right)}=D^2{U}_z^{\left(1\right)},$     (38)

$\left(-m\frac{D^3}{k^2}+\alpha \phi \right)U_z^{\left(2\right)}=-\frac{D^3}{k^2}{U}_z^{\left(1\right)},$     (39)

где m – отношение вязкости верхнего слоя к вязкости нижнего слоя (полупространства).

Граничные условия приводят к системе линейных уравнений для произвольных констант, входящих в общее решение (36) и (37):

$k{e}^k{C}_1+\left(1+k\right)e^k{C}_2+k{e}^{-k}{C}_3+\left(-1+k\right)e^{-k}{C}_4=\mathrm{0,}$

$\left(\phi +2\lambda km\right)e^k(C_1+C_2)+\left(\phi -2\lambda km\right)e^{-k}(C_3+C_4)=\mathrm{0,}$

$C_1+C_3-B_1=\mathrm{0,}$     (40)

$k{C}_1+C_2-k{C}_3+C_4-k{B}_1-B_2=\mathrm{0,}$

$k{C}_1+C_2+k{C}_3-C_4-k{B}_1-B_2=\mathrm{0,}$

$\left(\alpha \phi +2\lambda km\right)C_1+\left(\alpha \phi -2\lambda km\right)C_3-2\lambda k{B}_1=0.$

Приравнивая нулю определитель этой системы линейных уравнений, получаем характеристическое уравнение. Положив в этом уравнении $\alpha =\mathrm{0.01,}$ $\phi =0.005$ и $m={10}^3,$ находим показанную на рис. 4 зависимость инкремента  от волнового числа k. Как видно из рис. 4, наиболее неустойчивым является тяжелое включение, которому соответствует безразмерное волновое число $k\approx 1.5$ (такое включение имеет горизонтальный размер $\frac{\pi }{k}d\approx 20 км)$ и безразмерный инкремент $\lambda \approx {10}^{-8},$ который удовлетворяет условиям (34) и (35), когда $m={10}^3.$ Поскольку зависимость смещения от времени имеет вид $\exp \left(\lambda t\right),$ характерное время развития неустойчивости можно определить как 1/λ. Безразмерному инкременту $\lambda \approx {10}^{-8}$ соответствует огромное характерное время развития неустойчивости $\frac{\eta }{\lambda \mu }\approx 6\cdot {10}^{16}с\approx 2\cdot {10}^9лет.$ Если верхний холодный слой коры рассматривается как упругий, ему приписывается бесконечно большая вязкость, при которой тяжелые включения в этом слое устойчивы. Однако вязкость этого слоя хотя и велика, но не бесконечна. Поэтому тяжелые включения в этом верхнем слое коры неустойчивы, но эта неустойчивость развивается так медленно, что тяжелые включения в коре практически не меняют своего положения за период времени, сравнимый с возрастом Земли.

Рис. 4. Зависимость инкремента λ от волнового числа k для вязкой ньютоновской реологической модели.

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

 $2{\epsilon }_{ij}={\sigma }_{ij}f\left(t\right),$     (41)

где $f\left(t\right)$ – функция ползучести, дающая аналитическое описание неустановившейся ползучести; ${\epsilon }_{ij}$ – тензор деформаций, отсчитываемых от состояния в момент приложения напряжения. Для горных пород функция ползучести при высоких температурах хорошо описывается законом Андраде:

 $f\left(t\right)=t^{1/3}/A,$     (42)

где A – реологический параметр Андраде. На малых временах, неустановившаяся ползучесть подчиняется закону Ломнитца, но уже на временах порядка суток становится справедливым закон Андраде [Birger, 1998; Биргер, 2007]. Поэтому именно закон Андраде будет применяться при анализе неустойчивости коры, который предполагает исследование течения на больших временах, прошедших с момента возникновения начального возмущения.

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

$2{\epsilon }_{ij}=\underset{0}{\overset{t}{\int }}K\left(t_1\right){\sigma }_{ij}\left(t-t_1\right)d{t}_1,$     (43)

где t – момент наблюдения; $K\left(t\right)$ – интегральное ядро ползучести, определяемое функцией ползучести:

$K=df/dt.$     (44)

Как следует из (42) и (44), ядро ползучести, соответствующее закону Андраде, имеет вид:

$K\left(t\right)=\frac{1}{3A}{t}^{-2/3}.$     (45)

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

Когда деформации и напряжения зависят от времени как $\exp \left(\lambda t\right),$ правая часть уравнения (43) принимает вид:

$2{\epsilon }_{ij}={\sigma }_{ij}{e}^{\lambda t}\underset{0}{\overset{\infty }{\int }}K\left(t_1\right)e^{-\lambda t_1}d{t}_1,$     (46)

$\underset{0}{\overset{\infty }{\int }}K\left(t_1\right)e^{-\lambda t_1}d{t}_1=K^{*}\left(\lambda \right),$

где звездочка обозначает преобразование Лапласа, которое использовано здесь только для того, чтобы вычислить интеграл в уравнении (46). Преобразование Лапласа ядра ползучести (45) дает:

$K^{*}\left(\lambda \right)=\frac{\Gamma \left(1/3\right)}{3}\frac{{\lambda }^{-1/3}}{A}\approx \frac{{\lambda }^{-1/3}}{A},$     (47)

где гамма-функция $\Gamma \left(1/3\right)\approx 3.$ Линейная теория устойчивости рассматривает поведение механической системы на больших временах, прошедших с момента возникновения малого возмущения. Поэтому в уравнении (46) верхним пределом интегрирования является $t=\infty .$ Таким образом, на больших временах эффективный модуль сдвига среды Андраде имеет вид:

${\mu }_{eff}=A{\lambda }^{1/3}=A{\tau }^{-1/3},$     (48)

а эффективная ньютоновская вязкость записывается как

${\eta }_{eff}=A{\lambda }^{-2/3}=A{\tau }^{2/3},$     (49)

где $\tau =\frac{1}{\lambda }$ – характерное время процесса в случае, когда инкремент λ – положительное действительное число. Эффективная вязкость для среды Андраде зависит от характерного времени рассматриваемого процесса и, следовательно, эффективная вязкость, найденная для послеледниковых течений с характерным временем 1000 лет, не может применяться при исследовании более медленных процессов [Биргер, 2016].

Следует отметить, что соотношения (46)–(49), строго говоря, перестают быть справедливыми, когда λ – отрицательное число, поскольку в этом случае интеграл в правой части уравнения (46) расходится. Неустойчивость в среде Андраде описывается функцией $\exp \left(\lambda t\right),$ где $\lambda >\mathrm{0,}$ а устойчивость, которая имеет место в отсутствие скачка плотности, описывается заранее неизвестной функцией, убывающей с ростом t. Эту неизвестную функцию можно найти, последовательно применяя преобразование Лапласа, как это сделано в работе [Биргер, 2018]. Однако, если заменить $\lambda$ на $\left|\lambda \right|$ в уравнениях (48) и (49), эти уравнения можно применять в случае устойчивости как довольно грубое приближение (см. Приложение). Такое приближение позволяет исследовать устойчивость в среде Андраде, не применяя преобразование Лапласа и, используя обычную для линейной теории устойчивости зависимость от времени $\exp \left(\lambda t\right),$ находить не только положительные, но и отрицательные значения инкремента

Реологический параметр Андраде, сильно зависит от температуры. В нижнем горячем слое коры $A=4\cdot {10}^{12}Па\cdot {с}^{1/3},$ а в верхнем холодном слое этот параметр приблизительно на 3 порядка выше [Birger, 2013; Karato, 2008]. Как следует из (50), на временах $\tau \approx 1000лет\approx 3\cdot {10}^{10}с,$ характерных для мелкомасштабных послеледниковых течений, значению реологического параметра $A=4\cdot {10}^{12}Па\cdot {с}^{1/3}$ соответствует эффективная вязкость ${\eta }_{eff}\approx 4\cdot {10}^{19}Па\cdot с.$ Такая оценка согласуется с оценкой вязкости, полученной при рассмотрении мелкомасштабных послеледниковых течений в рамках реологической модели ньютоновской жидкости [Cathles, 1975].

Переходя в формуле (48) к безразмерному инкременту λ, получаем:

$\frac{{\mu }_{eff}}{\mu }=\frac{A}{\mu }{\left(\frac{\eta }{\mu }\right)}^{-1/3}{\lambda }^{1/3}\approx 0.07{\lambda }^{1/3},$     (50)

где $\eta =4\cdot {10}^{19}Па\cdot с$ и $\mu =6\cdot {10}^{10}Па$ – введенные ранее вязкость и упругий модуль, с помощью которых определен масштаб времени. При эффективном модуле упругости (50) уравнения (40) принимают вид:

$k{e}^k{C}_1+\left(1+k\right)e^k{C}_2+k{e}^{-k}{C}_3+\left(-1+k\right)e^{-k}{C}_4=\mathrm{0,}$

$\begin{array}{c}\left(\phi +0.14m{\lambda }^{1/3}k\right)e^k\left(C_1+C_2\right)+\\ +\left(\phi -0.14m{\lambda }^{1/3}k\right)e^{-k}\left(C_3+C_4\right)=\mathrm{0,}\end{array}$

 $C_1+C_3-B_1=\mathrm{0,}$     (51)

$k{C}_1+C_2-k{C}_3+C_4-k{B}_1-B_2=\mathrm{0,}$

$k{C}_1+C_2+k{C}_3-C_4-k{B}_1-B_2=\mathrm{0,}$

$\begin{array}{c}\left(\alpha \phi +0.14m{\lambda }^{1/3}k\right)C_1+\\ +\left(\alpha \phi -0.14m{\lambda }^{1/3}k\right)C_3-0.14{\lambda }^{1/3}k{B}_1=0.\end{array}$

Приравнивая нулю определитель этой системы линейных уравнений, получаем характеристическое уравнение. Положив в этом уравнении   $\alpha =\mathrm{0.01,}$ $\phi =0.005$ и $m={10}^3,$ находим показанную на рис. 5 зависимость инкремента $\lambda$ от волнового числа k.

Рис. 5. Зависимость инкремента $\lambda$ от волнового числа k для реологической модели Андраде.

Как видно из рис. 5, наиболее неустойчивым является тяжелое включение, которому соответствует безразмерное волновое число $k\approx 1.5$ и безразмерный инкремент $\lambda \approx 4\cdot {10}^{-21}.$ Такому крайне малому значению безразмерного инкремента соответствует сверхмедленное развитие неустойчивости: характерное время развития неустойчивости $\frac{\eta }{\lambda \mu }\approx 5\cdot {10}^{21}лет$ на много порядков превышает возраст Земли. В этом процессе эффективная ньютоновская вязкость нижнего горячего слоя коры чрезвычайно высока: как следует из (49), она приблизительно на 12 порядков выше эффективной вязкости, соответствующей процессу послеледниковых течений.

Заключение

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

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

В линейной теории устойчивости зависимость возмущения от времени ищется в виде $e^{\lambda t},$ где λ – комплексный инкремент. Исследуемой системе обычно соответствует целый набор инкрементов λ. Если действительные части всех инкрементов λ отрицательны или равны нулю, система устойчива, а если действительная часть хотя бы одного из инкрементов положительна, система неустойчива. Каждому значению инкремента λ соответствует процесс, длительность которого определяется временем $1/\left|\lambda \right|.$ В настоящем исследовании устойчивости тяжелых включений в верхней коре найдены три значения инкремента: $\left|{\lambda }_1\right|\gg \left|{\lambda }_2\right|\gg \left|{\lambda }_3\right|.$ Чисто мнимый инкремент ${\lambda }_1=i\omega$ соответствует упругому поведению верхней коры и подстилающего слоя. Отрицательный инкремент ${\lambda }_2,$ найденный в работе [Биргер, 2023], соответствует процессу, при котором верхняя кора ведет себя как упругая, а подстилающий слой демонстрирует вязкость. Положительный инкремент ${\lambda }_3$ характеризует очень медленный процесс, при котором оба слоя ведут себя как вязкие. Поскольку инкремент ${\lambda }_3$ положителен, тяжелые включения в верхнем слое коры неустойчивы. Однако эта неустойчивость развивается так медленно (инкремент ${\lambda }_3$ очень мал), что тяжелые включения практически не меняют своего положения за период времени, равный возрасту Земли.

Приложение

Рассмотрим случай, когда нет ни скачка плотности $\left(\alpha =\frac{\Delta \rho }{\rho }=0\right),$ ни скачка реологического параметра $(m=0).$ В рамках эффективной ньютоновской реологии зависимость от времени вертикального смещения в точке на верхней поверхности слоя записывается в виде:

$\frac{U_z\left(t\right)}{U_z\left(0\right)}=\exp \left(\lambda t\right),\lambda =-1/\tau ,$     (I)

 $\tau =\frac{2k{\eta }_{eff}}{\rho g},$     (II)

где λ – декремент, а $\tau$ – характерное время процесса, за которое смещение уменьшается в e раз. Применение реологической модели Андраде приводит к асимптотической зависимости вертикального смещения от времени:

$\frac{U_z\left(t\right)}{U_z\left(0\right)}=\frac{2kA{t}^{-1/3}}{\Gamma \left(2/3\right)\rho g},$ $\Gamma \left(2/3\right)\approx \mathrm{1.35.}$     (III)

Формула (III) получена в работе [Биргер, 2018] с помощью преобразования Лапласа и справедлива на больших временах, прошедших с момента начального возмущения. Хотя зависимость от времени, даваемая (III), не экспоненциальная, можно ввести время  за которое смещение уменьшается в e раз. Как следует из (III),

$\frac{1}{e}\approx \frac{2kA{\tau }^{-1/3}}{\Gamma \left(2/3\right)\rho g}.$     (IV)

Подразумевая под эффективной вязкостью вязкость такой ньютоновской среды, которая характеризуется таким же временем  как среда Андраде, и подставляя в (IV) соотношение:

$\frac{2k}{\rho g}=\frac{\tau }{{\eta }_{eff}},$

которое следует из (II), получим:

${\eta }_{eff}\approx \frac{e}{\Gamma \left(2/3\right)}A{\tau }^{2/3},$ $\frac{e}{\Gamma \left(2/3\right)}\approx 2.$     (V)

Поскольку $\tau =1/\left|\lambda \right|,$ эффективную вязкость и эффективный модуль сдвига для среды Андраде можно описать приближенными формулами

${\eta }_{eff}\approx A{\left|\lambda \right|}^{-2/3},$ ${\mu }_{eff}\approx A{\left|\lambda \right|}^{1/3}.$     (VI)

About the authors

B. I.  Birger

Author for correspondence.
Email: birgerbor@mail.ru
Russian Federation, Moscow

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Dependence of frequencies on wave numbers k at a density jump and with a dotted line at

Download (100KB)

Indexing metadata

3. Fig. 2. Dependence of frequencies on density jump for wave number

Download (93KB)

Indexing metadata

4. Fig. 3. Dependence of frequencies on density jump for wave number

Download (102KB)

Indexing metadata

5. Fig. 4. Dependence of the increment λ on the wave number k for the viscous Newtonian rheological model.

Download (90KB)

Indexing metadata

6. Fig. 5. Dependence of the growth rate on the wave number k for the Andrade rheological model.

Download (106KB)

Indexing metadata