The Effect of Heavy Inclusions in the Upper Crust on Thermoconvective Instability of the Lithosphere

Full Text

ВВЕДЕНИЕ

Литосфера, определяемая как пограничный слой, формируемый крупномасштабной мантийной конвекцией, значительно толще, чем холодный верхний слой коры, обладающий очень высокой эффективной вязкостью. Литосфера, определяемая таким образом, включает в себя не только верхние вязкие слои Земли, которые при не слишком медленных процессах ведут себя как упругие, но и астеносферу (слой с пониженной эффективной вязкостью, которая соответствует данным о послеледниковых поднятиях земной поверхности). Толщина конвективной литосферы порядка 100 км, а в некоторых регионах достигает 200 км. Несогласованности, возникающие из-за различных определений понятия “литосфера”, подробно обсуждаются в работах [Anderson, 1995; Karato, 2008; Биргер, 2016].

В тектонике плит обычно предполагается, что нижняя граница литосферной плиты определяется изотермой $(\approx 1300°\text{\hspace{0.17em}К}).$ Вертикальный градиент температуры в литосфере оценивается как $1.3°\text{\hspace{0.17em}К/км,}$ а ниже этого пограничного слоя градиент температуры почти нулевой. Актуальное для литосферы число Рэлея близко к критическому, т.е., к такому числу Рэлея, при котором происходит переход от конвективной устойчивости к конвективной неустойчивости. Те участки литосферы (континентальные кратоны), где ее толщина повышена, а градиент температуры понижен, наиболее неустойчивы, и в них возникает вторичная мелкомасштабная конвекция [Birger, 2013]. Критическое число Рэлея зависит от периодичности по горизонтали возникающего течения. Эта периодичность описывается волновым числом k, которое дает горизонтальный размер π/k конвективной ячейки. Волновое число определяется начальным возмущением равновесного состояния. При анализе влияния на конвективную устойчивость тяжелых включений в верхней коре предполагается, что утяжеленный участок верхней коры имеет горизонтальный размер L. Этот размер определяет волновое число $k=\pi /L.$ Когда L заметно меньше толщины литосферы (например, при  литосферу можно моделировать как полупространство, подстилающее верхнюю кору.

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

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

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

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

$\frac{d{p}_0}{dz}=-{\rho }_0\left(z\right)g,$ (1)

где g – ускорение силы тяжести, а плотность определена как:

${\rho }_0\left(z\right)=\rho +\Delta \rho$ при $0<z<d,$

${\rho }_0\left(z\right)=\rho$ при $-\infty <z<0.$

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

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

$-\frac{\partial p}{\partial z}+\frac{\partial {\sigma }_{xz}}{\partial x}+\frac{\partial {\sigma }_{zz}}{\partial z}+\rho g\alpha \theta =\rho \frac{\partial }{\partial t}{v}_z,$ (3)

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

$\frac{\partial \theta }{\partial t}-\frac{dT}{dz}{v}_z-\kappa \left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial z^2}\right)\theta =\mathrm{0,}$ (5)

где: p – возмущение давления; ${\sigma }_{xx},$ ${\sigma }_{xz}$ и ${\sigma }_{zz}$ – компоненты девиаторного тензора напряжений; $v_x$ и $v_z$ – скорости; $\theta$ – возмущение температуры; $\kappa$ – температуропроводность; $\frac{dT}{dz}$ – вертикальный градиент невозмущенной температуры, который предполагается однородным по глубине литосферы; $\rho$ – плотность; $\alpha$ – коэффициент теплового расширения; g – гравитационное ускорение. Скорости, напряжения, возмущения давления и температуры являются функциями вертикальной пространственной координаты z, горизонтальной координаты x и времени t. Уравнения (2) и (3) описывают двумерное движение среды (движение происходит в плоскости xz) с учетом силы Архимеда. Уравнение (4) представляет собой условие несжимаемости среды, а уравнение (5) является уравнением теплового баланса. К уравнениям (1)–(5) добавляются реологические уравнения:

${\sigma }_{xx}=2\eta {\dot{\epsilon }}_{xx},$ ${\sigma }_{zz}=2\eta {\dot{\epsilon }}_{zz},$ ${\sigma }_{xz}=2\eta {\dot{\epsilon }}_{xz},$  (6)

где $\eta$ – коэффициент вязкости. Уравнения (6) связывают девиаторные напряжения со скоростями деформаций, которые определены как:

${\dot{\epsilon }}_{xx}=\frac{\partial v_x}{\partial x},$ ${\dot{\epsilon }}_{zz}=\frac{\partial v_z}{\partial z},$ ${\dot{\epsilon }}_{xz}=\frac{1}{2}\left(\frac{\partial v_x}{\partial z}+\frac{\partial v_z}{\partial x}\right).$  (7)

Уравнения (1)–(7) справедливы и для верхнего, и для нижнего слоя, моделированного как полупространство. Однако плотности и реологические свойства слоев отличаются. Чтобы учесть скачок плотности на границе между слоями, необходимо ввести граничные условия.

Условие несжимаемости для вязкой среды, при выполнении которого справедливо уравнение (4), записываются в виде:

$\tau \gg \frac{\eta }{K},$ (8)

где τ – характерное время рассматриваемого течения, а K – модуль всестороннего сжатия. При выполнении условия

$\tau \gg \frac{\rho d^2}{\eta }$ (9)

в уравнениях (2) и (3) могут быть отброшены инерционные члены, что законно при достаточно медленном течении с большим характерным временем $\tau .$ Такие течения называются ползущими, а конвективные течения в Земле являются типичным примером ползущих течений. Кроме того, использовано приближение Буссинеска, в рамках которого можно пренебречь механической сжимаемостью среды, а термическую сжимаемость учитывать только в уравнениях движения. Литосфера, проявляющая и упругие, и вязкие свойства, описывается вязкоупругой реологической моделью Максвелла. Упругостью можно пренебречь при достаточно медленном течении, когда

$\tau \gg \eta /\mu ,$ (10)

где $\mu$ – упругий модуль сдвига, а $\eta /\mu$ называют временем Максвелла.

Для перехода к безразмерным переменным, введем следующие масштабы. Масштаб длины – толщина слоя d; масштаб времени – $d^2\text{}/\kappa ,$ где $\kappa$ – температуропроводность; масштаб скорости $\kappa /d;$ –  масштаб давления (и напряжений) – $\kappa \eta /d^2.$ В качестве масштаба температуры принят перепад температуры $\Delta T=\frac{dT}{dz}d$ между горячей нижней и холодной верхней поверхностями слоя.

Литосфера характеризуется следующими значениями физических параметров [Теркот, Шуберт, 1985]:

$\frac{dT}{dz}=1.3\cdot {10}^{-2}°К\text{}/\text{м,}$ $\alpha =3\cdot {10}^{-5}°{\text{\hspace{0.05em}К}}^{-1},$

$\rho =3\cdot {10}^3\text{\hspace{0.33em}кг\hspace{0.05em}}/{\text{\hspace{0.05em}м}}^3,$ $g=10{\text{\hspace{0.33em}мс}}^{-2},$ (11)

$\kappa ={10}^{-5}{\text{\hspace{0.33em}м}}^2{\text{с}}^{-1},$ $\mu =6\cdot {10}^{10}\text{\hspace{0.33em}Па,}$ $K=3\cdot {10}^{11}\text{\hspace{0.33em}Па.}$  

Эффективная вязкость холодной верхней коры на несколько порядков выше, чем эффективная вязкость подстилающих ее горячих слоев литосферы. Толщина верхней коры, которая принимается в качестве масштаба длины, оценивается как $d={10}^4\text{\hspace{0.33em}м.}$

Исследуя устойчивость тяжелой верхней коры методом линейной теории устойчивости, вертикальную скорость представим в виде:

$v_z=V_z\exp \left(\lambda t\right)\text{exp}\left(ikx\right),$ $k>\mathrm{0,}$ (12)

где λ – комплексный инкремент, k – действительное волновое число. В аналогичном виде представим и все остальные физические переменные. Такое представление позволяет свести систему уравнений в частных производных (2)–(7) к системе обыкновенных дифференциальных уравнений, в которой все переменные, характеризующие скорости деформации, напряжения и давление, зависят только от вертикальной координаты z. Характерное время τ для течения с инкрементом  представимо как $\tau =\frac{1}{\left|\lambda \right|}.$

Перейдя к безразмерным переменным в уравнениях (2)–(7) и подставляя соотношение (12) и аналогичные соотношения для других физических переменных в уравнения (2) и (4)–(7), получим соотношения, связывающие амплитуды горизонтальной скорости, давления и компонент тензора девиаторных напряжений с амплитудой вертикальной скорости:

$V_x=\frac{i}{k}D{V}_z,$ (13)

$P=\frac{1}{k^2}\left(D^3-k^{2}D\right)V_z,$ (14)

${\Sigma }_{zz}=-{\Sigma }_{xx}=2D{V}_z,$ (15)

${\Sigma }_{xz}=\frac{i}{k}\left(D^2+k^2\right)V_z,$ (16)

$\Theta =\frac{{\left(D^2-k^2\right)}^2}{\text{Ra\hspace{0.17em}}{k}^2}.$ (17)

В уравнениях (13)–(17) введен дифференциальный оператор $D=\frac{d}{dz}$ и число Рэлея:

$\text{Ra}=\rho \alpha g\frac{dT}{dz}{d}^4/\eta \kappa .$ (18)

Подставляя соотношения (13)–(18) в уравнение (3), приходим к уравнению для вертикальной скорости в нижнем слое, моделируемом как полупространство:

$-\infty <z<\mathrm{0,}$ $\left[{\left(D^2-k^2\right)}^3-\lambda {\left(D^2-k^2\right)}^2+\text{Ra\hspace{0.17em}}{k}^2\right]V_z=0.$ (19)

Более вязкий верхний слой описывается уравнением:

$0<z<\mathrm{1,}$ $\left[{\left(D^2-k^2\right)}^3-\lambda {\left(D^2-k^2\right)}^2+\frac{\text{Ra}}{m}{k}^2\right]V_z=\mathrm{0,}$ (20)

где $m\approx {10}^3$ – отношение вязкости коры к вязкости подстилающей литосферы.

Граничные условия на верхней деформируемой поверхности слоя $z=1$ имеют вид:

$-P+{\Sigma }_{zz}+\phi U_z=\mathrm{0,}$ (21)

${\Sigma }_{xz}=\mathrm{0,}$ (22)

$\Theta =U_z,$ (23)

где $U_z$ – вертикальное смещение верхней границы слоя, т.е. отклонение граничной поверхности от плоскости $z=1.$ Уравнения (21) и (22) следуют из условия обращения в нуль силы, действующей на единицу площади возмущенной поверхности слоя, а уравнение (23) – из условия обращения в нуль возмущения температуры на этой поверхности. В уравнения (21) и (23) входит смещение границы $U_z$ в силу того, что в состоянии равновесия в слое имеется вертикальный градиент давления и температуры. Условие (23) предполагает, что движение границы определяется движением материальных точек, расположенных на этой границе. В уравнении (21) введен безразмерный параметр:

$\phi =\rho g{d}^3\text{}/m\eta \kappa ,$ (24)

который вместе с числом Рэлея характеризует задачу о конвективной неустойчивости слоя с деформируемой свободной поверхностью. Используя параметр $\phi ,$ число Рэлея можно записать как $\text{Ra}=m\phi \alpha \frac{dT}{dz}d.$ Параметр описывает подвижность границы и представляет собой отношение дополнительного гидростатического давления, вызванного смещением границы, к характерному вязкому напряжению в слое. Чем меньше параметр ϕ, тем подвижнее граница. При очень большом значении ϕ верхняя поверхность слоя ведет себя как неподвижная граница.

В рамках линейной теории устойчивости амплитуда смещения $U_z$ связана с амплитудой скорости как:

$V_z=\lambda U_z.$ (25)

В критическом режиме конвективной неустойчивости, когда  граничные условия (21)–(23) принимают вид:

$V_z={\Sigma }_{xz}=\Theta =0.$ (26)

Уравнения (26) остаются справедливыми в режиме, близком к критическому, когда инкремент $\lambda$ очень мал. Точнее говоря, верхняя граница слоя неподвижна и граничные условия на ней принимают вид (26), когда $\frac{\phi }{\lambda }\gg \mathrm{1,}$ т.е. безразмерный инкремент λ должен удовлетворять условию:

$\lambda \ll \phi .$ (27)

Уравнения (26) описывают рэлеевские граничные условия, при которых граница слоя предполагается неподвижной, изотермической (возмущения температуры на границе обращаются в нуль) и свободной от касательных напряжений. Такую границу обычно называют “свободной” несмотря на то, что нормальные напряжения не обращаются в нуль на этой границе. Как следует из (16) и (17), граничные условия (26) можно записать как

$V_z=D^2{V}_z=D^4{V}_z=0.$ (28)

На границе $z=0$ между верхним слоем, которому приписан индекс 2, и полупространством, которому приписан индекс 1, непрерывны скорости

$V_x^{\left(2\right)}=V_x^{\left(1\right)},$ $V_z^{\left(2\right)}=V_z^{\left(1\right)}$ (29)

и сила, действующая на единицу площади поверхности, откуда следует:

${\Sigma }_{xz}^{\left(2\right)}={\Sigma }_{xz}^{\left(2\right)},$ $-P^{\left(2\right)}+{\Sigma }_{zz}^{\left(2\right)}+R_{\rho }{U}_z=-P^{\left(1\right)}+{\Sigma }_{xz}^{\left(1\right)},$ (30)

где безразмерный параметр $R_{\rho },$ который является аналогом числа Рэлея при гравитационной, а не конвективной неустойчивости, определен как:

$R_{\rho }=\frac{g{d}^3\Delta \rho }{\eta \kappa }=\frac{\Delta \rho }{\rho }m\phi .$ (31)

Кроме того, на границе $z=0$ непрерывны возмущение температуры и теплового потока

${\Theta }^{\left(2\right)}={\Theta }^{\left(1\right)},$ $D{\Theta }^{\left(2\right)}=D{\Theta }^{\left(1\right)}.$ (32)

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

$V_z^{\left(1\right)}=D{V}_z^{\left(1\right)}={\Theta }^{\left(1\right)}=0.$ (33)

Таким образом, нижняя граница литосферы, рассматриваемой как пограничный слой крупномасштабной конвекции, является изотермической (возмущение температуры равно нулю на этой границе) и “твердой” (нулевые компоненты скорости на этой границе). Граничные условия (33) не справедливы в узких областях литосферы, расположенных у границ литосферных плит, под которыми концентрируются вызванные крупномасштабной мантийной конвекцией потоки горячего или холодного материала.

РЕШЕНИЕ ДЛЯ НЬЮТОНОВСКОЙ РЕОЛОГИЧЕСКОЙ МОДЕЛИ

Общее решение уравнений (19), (20) записывается как:

$V_z=e^{q_{1}z}{C}_1+e^{q_{2}z}{C}_2+e^{q_{3}z}{C}_3+e^{-q_{1}z}{C}_4+e^{-q_{2}z}{C}_5+e^{-q_{3}z}{C}_6,$ $0\le z\le 1.$ (34)

$V_z=e^{{\tilde{q}}_{1}z}{B}_1+e^{{\tilde{q}}_{2}z}{B}_2+e^{{\tilde{q}}_{3}z}{B}_3,$ $-\infty \le z\le \mathrm{0,}$ (35)

где $C_i$ и $B_i$ – произвольные константы, а $q_i$ и ${\tilde{q}}_i$ являются корнями алгебраических уравнений:

${\left({q_i}^2-k^2\right)}^3-\lambda {\left({q_i}^2-k^2\right)}^2+\frac{\text{Ra}}{m}{k}^2=\mathrm{0,}$ (36)

${\left({{\tilde{q}}_i}^2-k^2\right)}^3-\lambda {\left({{\tilde{q}}_i}^2-k^2\right)}^2+\text{Ra\hspace{0.17em}}{k}^2=0.$ (37)

Раскладывая аналитические выражения для корней уравнений (36) и (37) в ряд по малому инкременту  и сохраняя только два первых члена разложения, получаем

$q_1=\sqrt{k^2-{\left(\frac{\text{Ra}}{m}{k}^2\right)}^{1/3}}+\frac{\lambda }{6\sqrt{k^2-{\left(\frac{\text{Ra}}{m}{k}^2\right)}^{1/3}}},$  

$\begin{array}{c}{q}_{\mathrm{2,}3}=\sqrt{k^2+\frac{1\pm i\sqrt{3}}{2}{\left(\frac{\text{Ra}}{m}{k}^2\right)}^{1/3}}+\\ +\frac{\lambda }{6\sqrt{k^2+\frac{1\pm i\sqrt{3}}{2}{\left(\frac{\text{Ra}}{m}{k}^2\right)}^{1/3}}},\end{array}$ (38)

${\tilde{q}}_1=\sqrt{k^2-{\left(\text{Ra\hspace{0.17em}}{k}^2\right)}^{1/3}}+\frac{\lambda }{6\sqrt{k^2-{\left(\text{Ra\hspace{0.05em}}{k}^2\right)}^{1/3}}},$

$\begin{array}{c}{\tilde{q}}_{\mathrm{2,3}}=\sqrt{k^2+\frac{1\pm i\sqrt{3}}{2}{\left(\text{Ra\hspace{0.05em}}{k}^2\right)}^{1/3}}+\\ +\frac{\lambda }{6\sqrt{k^2+\frac{1\pm i\sqrt{3}}{2}{\left(\text{Ra\hspace{0.05em}}{k}^2\right)}^{1/3}}}.\end{array}$

Уравнение (35) следует как из условий (33), так и из требования ограниченности скоростей в подстилающем слой полупространстве. Согласно (35), (38), при достаточно малом числе Рэлея компоненты скорости и возмущения температуры в полупространстве $z<0$ убывают с глубиной как $e^{kz}.$ Следовательно, смещения проникают в моделирующее литосферу полупространство на глубину, определяемую волновым числом k. Зафиксируем $k=\mathrm{1,}$ что соответствует горизонтальному размеру тяжелого включения в верхней коре $\frac{\pi }{k}d,$ равному 30 км. При значительно меньших значениях k, например, при $k=\mathrm{0.1,}$ конвективное течение проникало бы в слои, подстилающие литосферу, и нельзя было бы моделировать литосферу как полупространство.

Граничные условия (26)–(33) приводят к системе линейных алгебраических уравнений для произвольных констант в уравнениях (34) и (35):

$e^{q_1}{C}_1+e^{q_2}{C}_2+e^{q_3}{C}_3+e^{-q_1}{C}_4+e^{-q_2}{C}_5+e^{-q_3}{C}_6=0$

$\begin{array}{c}{q}_1^2{e}^{q_1}{C}_1+q_2^2{e}^{q_2}{C}_2+q_3^2{e}^{q_3}{C}_3+q_1^2{e}^{-q_1}{C}_4+\\ +q_2^2{e}^{-q_2}{C}_5+q_3^2{e}^{-q_3}{C}_6=0\end{array}$

$\begin{array}{c}{q}_1^4{e}^{q_1}{C}_1+q_2^4{e}^{q_2}{C}_2+q_3^4{e}^{q_3}{C}_3+q_1^4{e}^{-q_1}{C}_4+\\ +q_2^4{e}^{-q_2}{C}_5+q_3^4{e}^{-q_3}{C}_6=0\end{array}$

$C_1+C_2+C_3+C_4+C_5+C_6-B_1-B_2-B_3=0$

$\begin{array}{c}{q}_1{C}_1+q_2{C}_2+q_3{C}_3-q_1{C}_4-q_2{C}_5-\\ -q_3{C}_6-{\tilde{q}}_1{B}_1-{\tilde{q}}_2{B}_2-{\tilde{q}}_3{B}_3=0\end{array}$

$\begin{array}{c}m\left(q_1^2+k^2\right)C_1+m\left(q_2^2+k^2\right)C_2+m\left(q_3^2+k^2\right)C_3+\\ +m\left(q_1^2+k^2\right)C_4+m\left(q_2^2+k^2\right)C_5+m\left(q_3^2+k^2\right)C_6-\\ -\left({\tilde{q}}_1^2+k^2\right)B_1-\left({\tilde{q}}_2^2+k^2\right)B_2-\left({\tilde{q}}_3^2+k^2\right)B_3=0\end{array}$ (39)

$\begin{array}{c}\left[m\lambda \left(-q_1^3+3q_1\right)+R_{\rho }\right]C_1+\left[m\lambda \left(-q_2^3+3q_2\right)+R_{\rho }\right]C_2+\\ +\left[m\lambda \left(-q_3^3+3q_3\right)+R_{\rho }\right]C_3-\left[m\lambda \left(-q_1^3+3q_1\right)+R_{\rho }\right]C_4-\\ -\left[m\lambda \left(-q_2^3+3q_2\right)+R_{\rho }\right]C_5-\left[m\lambda \left(-q_3^3+3q_3\right)+R_{\rho }\right]C_6-\\ -\lambda \left(-{\tilde{q}}_1^3+3{\tilde{q}}_1\right)B_1-\lambda \left(-{\tilde{q}}_2^3+3{\tilde{q}}_2\right)B_2-\lambda \left(-{\tilde{q}}_3^3+3{\tilde{q}}_3\right)B_3=0\end{array}$

$\begin{array}{c}m{\left(q_1^2-k^2\right)}^2{C}_1+m{\left(q_2^2-k^2\right)}^2{C}_2+m{\left(q_3^2-k^2\right)}^2{C}_3+\\ +m{\left(q_1^2-k^2\right)}^2{C}_4+m{\left(q_2^2-k^2\right)}^2{C}_5+m{\left(q_3^2-k^2\right)}^2{C}_6-\\ -{\left({\tilde{q}}_1^2-k^2\right)}^2{B}_1-{\left({\tilde{q}}_2^2-k^2\right)}^2{B}_2-{\left({\tilde{q}}_3^2-k^2\right)}^2{B}_3=0\end{array}$

$\begin{array}{c}m{q}_1{\left(q_1^2-k^2\right)}^2{C}_1+m{q}_2{\left(q_2^2-k^2\right)}^2{C}_2+\\ +m{q}_3{\left(q_3^2-k^2\right)}^2{C}_3-m{q}_1{\left(q_1^2-k^2\right)}^2{C}_4-\\ -m{q}_2{\left(q_2^2-k^2\right)}^2{C}_5-m{q}_3{\left(q_3^2-k^2\right)}^2{C}_6-\\ -{\tilde{q}}_1{\left({\tilde{q}}_1^2-k^2\right)}^2{B}_1-{\tilde{q}}_2{\left({\tilde{q}}_2^2-k^2\right)}^2\times \\ \times B_2-{\tilde{q}}_3{\left({\tilde{q}}_3^2-k^2\right)}^2{B}_3=0.\end{array}$

Литосфера, рассматриваемая как пограничный слой крупномасштабной мантийной конвекции, находится в режиме, близком к порогу конвективной неустойчивости. В этом пороговом режиме $\lambda =\mathrm{0,}$ а число Рэлея, при котором реализуется этот режим, называется критическим. Подставив $m={10}^3,$ $k=\mathrm{1,}$ $R_{\rho }=0$ и $\lambda =0$ в систему линейных уравнений (39) и положив определитель этой системы равным нулю (только в этом случае система уравнений имеет ненулевое решение для $C_i$ и $B_i),$ получаем алгебраическое уравнение с одним неизвестным, которым является критическое число Рэлея, соответствующее волновому числу $k=1.$ Решая это алгебраическое уравнение, находим

$\text{Ra}={\text{Ra}}_{cr}\left(1\right)=4\cdot {10}^{-2}.$ (40)

Как следует из (18), такому значению Ra соответствует эффективная вязкость $\eta =3\cdot {10}^{20}\text{\hspace{0.33em}Па}\cdot с\text{,}$ при которой безразмерные параметры $\phi$ и $R_{\rho }$ оцениваются как:

$\phi \approx {10}^{-2},$ $R_{\rho }=\frac{g{d}^3\Delta \rho }{\eta \kappa }=\frac{\Delta \rho }{\rho }m\phi \approx 10\frac{\Delta \rho }{\rho }.$

При таких значениях параметров Ra и $R_{\rho },$ приравнивая нулю определитель системы уравнений (39), в которых $m={10}^3,$ $k=\mathrm{1,}$ получаем алгебраическое уравнение для неизвестного инкремента $\lambda .$ Решая это уравнение, находим в первом приближении по инкременту ƛ, который априори предполагается малым,

$\lambda \approx \left(2.4\cdot {10}^{-4}+5.9\cdot {10}^{-10}i\right)R_{\rho }.$ (41)

С ростом скачка плотности неустойчивость нарастает. Таким образом, тяжелая верхняя кора повышает, хотя и очень слабо, конвективную неустойчивость литосферы, и приводит к чрезвычайно медленным колебаниям. Учитывая принятый масштаб времени $\frac{d^2}{\kappa }={10}^{13}ñ\approx 3\cdot {10}^5\text{\hspace{0.33em}лет}$ и оценку $R_{\rho }=10\frac{\Delta \rho }{\rho },$ период этих колебаний при максимальном значении скачка плотности $\frac{\Delta \rho }{\rho }={10}^{-2}$ можно оценить как ${10}^{16}лет\text{.}$ Следует отметить, что относительное изменение плотности, вызванное скачком плотности на границе верхняя кора–литосфера, имеет тот же порядок величины, что и относительное изменение плотности, $\frac{\Delta \rho }{\rho }=\alpha \frac{dT}{dz}d\approx {410}^{-3},$ вызванное градиентом температуры.

Физический смысл возникающих колебаний можно объяснить следующим образом. Верхний тяжелый слой $(\Delta \rho >0)$ останавливает конвективное течение, возникающее из-за подогрева снизу, а затем меняет его направление на противоположное (восходящий и нисходящий потоки меняются местами). Это связано с тем, что тяжелому верхнему слою выгоднее погружаться вниз в том месте, где под ним находится более горячая и легкая среда.

В случае, когда плотность нижнего слоя превышает плотность верхнего слоя, т.е. когда $\Delta \rho <\mathrm{0,}$действительная часть комплексного инкремента $\lambda$ становится отрицательной (конвективная устойчивость) [Richter, Johnson, 1974; Биргер, 1988]. Колебания возникают и в этом случае: легкий верхний слой останавливает нисходящий конвективный поток и обращает направление конвективного течения.

Найденное решение (41) соответствует заранее сделанным предположениям: безразмерный инкремент $\lambda$ является достаточно малым и удовлетворяет условию (27), а характерное время процесса, т.е. его период удовлетворяет требованиям (8)–(10). Отметим, что после введения масштаба времени $d^2/\kappa$ условие (9), позволяющее пренебречь инерционностью, принимает вид $\lambda \ll \text{Pr}.$ Число Прандтля определено как $\text{Pr}=\frac{\nu }{\kappa },$ где $\nu =\frac{\eta }{\rho }$ – кинематическая вязкость. При выполнении условия $\lambda \gg \text{Ra}$ можно пренебречь имеющимся в литосфере вертикальным градиентом температуры. Данное условие выполняется при исследовании достаточно быстрых процессов [Биргер, 2023], но конечно не выполняется при анализе конвективной устойчивости.

РЕШЕНИЕ ДЛЯ РЕОЛОГИЧЕСКОЙ МОДЕЛИ АНДРАДЕ

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

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

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

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

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

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

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

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

$K=df/dt.$ (45)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Для среды Андраде введем масштаб вязкости – ${\eta }_A=A{\left(d^2/\kappa \right)}^{2/3},$ и тогда при экспоненциальной зависимости от времени (20), безразмерная эффективная вязкость принимает вид

$F=F\left(\lambda \right)={\lambda }^{-2/3},$ (51)

а число Рэлея для среды Андраде определено как

$\text{Ra}=\frac{\rho \alpha g\frac{dT}{dz}{d}^4}{{\eta }_A\kappa }=\rho g\alpha \frac{dT}{dz}{d}^2{\left(d^2/\kappa \right)}^{1/3}/A.$ (52)

Подставляя оценки (13) в (52), получаем $\text{Ra}\approx 2.5\cdot {10}^{10}/A.$ Если оценивать значение параметра Андраде как $4\cdot {10}^{12}\text{\hspace{0.33em}Па}\cdot {с}^{1/3},$ число Рэлея оценивается как $\text{Ra}\approx 6.3\cdot {10}^{-3}.$ Масштабная вязкость ${\eta }_A,$ введенная для реологии Андраде, почти на порядок выше, чем вязкость $\eta =3\cdot {10}^{20}\text{\hspace{0.33em}Па}\cdot с\text{,}$ введенная в рассмотренном выше случае ньютоновской реологии. Поэтому значение числа Рэлея почти на порядок меньше того, которое использовались при ньютоновской реологии.

Подставляя соотношения (20) в уравнения (15), (16), (17) и (19), получим соотношения, связывающие амплитуды давления, температуры и компонент тензора девиаторных напряжений с амплитудой вертикальной скорости,

$P=\frac{1}{k^2}F\left(\lambda \right)\left(D^3-k^{2}D\right)V_z,$ (53)

${\Sigma }_{zz}=-{\Sigma }_{xx}=2F\left(\lambda \right)D{V}_z,$

${\Sigma }_{zz}=\frac{i}{k}F\left(\lambda \right)\left(D^2+k^2\right)V_z,$

$\Theta =\frac{F\left(\lambda \right)}{\text{Ra\hspace{0.05em}}{k}^2}{\left(D^2-k^2\right)}^2$

Соотношения (53) заменяют соотношения (13)–(17), соответствующие ньютоновской реологии. Исключая из уравнений амплитуды все физические переменные, кроме амплитуды вертикальной скорости $V_z,$ приходим к обыкновенным дифференциальным уравнениям

$\left[{\lambda }^{-2/3}{\left(D^2\text{}-k^2\right)}^3-{\lambda }^{1/3}{\left(D^2-k^2\right)}^2+\frac{\text{Ra}}{m}{k}^2\right]V_z=\mathrm{0,}$ (54)

$\left[{\lambda }^{-2/3}{\left(D^2-k^2\right)}^3-{\lambda }^{1/3}{\left(D^2-k^2\right)}^2+\text{Ra\hspace{0.05em}}{k}^2\right]V_z=0.$ (55)

Уравнение (54) относится к верхнему тяжелому слою, а уравнение (55) – к подстилающему полупространству. Эти уравнения сильно отличаются от уравнений (36) и (37), соответствующих вязкой ньютоновской реологии. Заранее предполагая, что искомый инкремент мал $\left(\left|\lambda \right|\ll 1\right),$ в уравнениях (54) и (55) можно пренебречь членами, содержащими ${\lambda }^{1/3}.$ Результат расчетов подтвердил это предположение.

Решения уравнений (54) и (55) по-прежнему ищем в виде (36), (37), но теперь выражения для $q_i$и ${\tilde{q}}_1$ принимают вид

$q_1=k-\frac{1}{2k}{\lambda }^{2/9}{\left(\frac{\text{Ra}}{m}{k}^2\right)}^{1/3},$ 

$q_{\mathrm{2,3}}=k+\frac{1\pm i\sqrt{3}}{4k}{\lambda }^{2/9}{\left(\frac{\text{Ra}}{m}{k}^2\right)}^{1/3},$ (56)

${\tilde{q}}_1=k-\frac{1}{2k}{\lambda }^{2/9}{\left(\text{Ra\hspace{0.05em}}{k}^2\right)}^{1/3},$

${\tilde{q}}_{\mathrm{2,3}}=k+\frac{1\pm i\sqrt{3}}{4k}{\lambda }^{2/9}{\left(\text{Ra\hspace{0.05em}}{k}^2\right)}^{1/3},$

где число Рэлея определено формулой (52). Уравнения (39) сохраняют свой вид и в случае реологии Андраде, но в том уравнении, которое содержит $R_{\rho },$ надо заменить $\lambda$ на ${\lambda }^{1/3}.$ Подставляя в уравнения (39) соотношения (56) и определяя безразмерный параметр $R_{\rho }$ как

$R_{\rho }=\frac{\Delta \rho }{\rho }\frac{\rho gd}{A}{\left(\frac{d^2}{\kappa }\right)}^{1/3}\approx 1.6\frac{\Delta \rho }{\rho },$ (57)

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

Решая характеристическое уравнение (расчеты проведены для фиксированного значения $k=\mathrm{1,}$а отношение коэффициента Андраде верхнего слоя к коэффициенту Андраде подстилающих слоев оценивается $m={10}^3),$ находим в первом приближении по инкременту ƛ, который априори предполагается малым,

$\lambda ={\lambda }_0+R_{\rho }{\lambda }_1,$ (58)

где

${\lambda }_1=\left(4.1+0.21i\right)\cdot {10}^{-14},$ (59)

а ${\lambda }_0$ – значение инкремента, соответствующее конвективной устойчивости литосферы, в которой присутствует вертикальный градиент температуры, характеризуемый числом Ra, но нет скачка плотности. В работах [Birger, 1998; 2013] получены численные оценки этого конвективного инкремента, который в среде Андраде является комплексным числом, и, следовательно, конвективная неустойчивость литосферы имеет колебательный характер. Вертикальный градиент температуры в литосфере достаточно высок, а эффективная вязкость, соответствующая неустановившейся ползучести в литосфере, достаточно мала для того, чтобы вызвать мелкомасштабную конвективную неустойчивость тех регионов континентальной литосферы, где ее толщина велика. Наиболее утолщенные регионы литосферы расположены под континентальными кратонами, где ее толщина порядка 200 км.

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

Оценка инкремента ${\lambda }_0,$ полученная в результате решения характеристического уравнения, в принятом масштабе времени $3\cdot {10}^5$ лет дает

${\lambda }_0\approx 0.15+0.06i.$ (60)

Таким образом, период конвективных колебаний составляет $\frac{2\pi }{\text{Im}{\lambda }_0}\cdot 3\cdot {10}^5\text{\hspace{0.33em}лет}=3\cdot {10}^7$ лет. Когда $\frac{\Delta \rho }{\rho }={10}^{-2}$ и, следовательно, $R_{\rho }={1.610}^{-2},$

$R_{\rho }{\lambda }_1=\left(6.6+0.33i\right)\cdot {10}^{-16}.$ (61)

Как уже отмечено выше, в работе расчеты, которые приводят к формулам (59)–(61), проведены для фиксированного значения $k=1$ (горизонтальный размер тяжелого включения в верхней коре $\frac{\pi }{k}d$ составляет 30 км), а отношение коэффициента Андраде верхнего слоя к коэффициенту Андраде подстилающих слоев оценивалось как $m={10}^3.$

Таким образом, скачок плотности  увеличивает действительную и мнимую части инкремента $\lambda ,$повышая неустойчивость и увеличивая частоту колебаний. Как следует из (61), поправка $R_{\rho }{\lambda }_1$ к значению инкремента ${\lambda }_0,$ определяемому из анализа конвективной устойчивости, настолько мала, что ею можно пренебречь. Тяжелые включения в верхней коре практически не меняют состояние колебательной неустойчивости, в котором находится подогреваемая снизу литосфера. Можно рассматривать полученную поправку к значению ${\lambda }_0$ как инкремент $\lambda =R_{\rho }{\lambda }_1,$ определяющий гравитационную неустойчивость, которая вызвана повышенной плотностью  верхней коры. Формула (61), дающая значение этого инкремента, получена с учетом вертикального градиента температуры. Без учета градиента температуры этот инкремент оценивается как $\lambda \approx 3\cdot {10}^{-18}.$ Сравнение этого значения $\lambda$ с формулой (61) показывает, что учет вертикального градиента температуры увеличивает действительную часть инкремента и приводит к появлению мнимой части инкремента, т.е., к колебательной неустойчивости. Инкремент (58), описывающий неустойчивость литосферы, представляет собой сумму инкрементов, первый из которых описывает конвективную устойчивость, а второй – гравитационную неустойчивость. Причем вклад гравитационного инкремента очень мал.

Результаты, полученные с учетом неустановившейся ползучести, описанной с помощью реологической модели Андраде, сильно отличаются от результатов, полученных в рамках ньютоновской реологической модели. В ньютоновской модели вязкость не зависит от длительности рассматриваемого процесса и, следовательно, оценку вязкости, найденную при рассмотрении послеледниковых течений, можно использовать при изучении значительно более медленного процесса развития конвективной неустойчивости. В модели Андраде это не так: эффективная вязкость при конвективном движении значительно выше, чем при послеледниковых поднятиях [Биргер, 2016; 2018]. Поэтому инкремент ƛ, найденный в модели Андраде, значительно ниже, чем тот, который получен в вязкой ньютоновской модели. Но главное отличие в другом. В вязкой ньютоновской жидкости инкремент ƛ, описывающий конвективную устойчивость и неустойчивость, является действительным числом, и $\lambda =0$ в критическом режиме перехода от устойчивости к неустойчивости [Гершуни, Жуховицкий, 1972]. В среде Андраде этот инкремент является комплексным числом, а в критическом режиме $\lambda =i\omega ,$ где $\omega$ – действительная частота [Birger, 1998; 2013]. Поэтому, если в рамках ньютоновской модели колебания возникают только при тяжелых включениях в верхнем слое, то в рамках модели Андраде конвективные колебания имеют место и без тяжелых включений, появление которых дает только небольшую поправку к частоте конвективных колебаний. Если под каким-то регионом земной поверхности литосфера находится в режиме пороговой устойчивости (действительная часть конвективного инкремента равна нулю), тяжелые включения приводят к очень слабой неустойчивости, которая развивается так медленно, что верхняя кора практически не меняют своего положения за период времени, сравнимый с возрастом Земли.

ЗАКЛЮЧЕНИЕ

Методом линейной теории для малых возмущений решена задача об устойчивости тяжелых включений, которые расположены в верхнем слое коры и имеют небольшой горизонтальный размер. Тяжелые включения вызывают быстрые и медленные процессы, характеризуемые большими и малыми значениями инкрементов. Если при исследовании быстрых процессов можно не учитывать тепловые эффекты, то при исследовании медленного процесса, исследуемого в этой работе, можно не учитывать упругость и инерционность среды, но необходимо учитывать ползучесть среды и вертикальный градиент температуры, присутствующий в литосфере. При малых деформациях, рассматриваемых в рамках линейной теории устойчивости, имеет место неустановившаяся ползучесть, которая описывается реологической моделью Андраде. При такой реологии конвективная неустойчивость утолщенных участков литосферы имеет колебательный характер: инкремент является комплексным числом, действительная часть которого близка к нулю (пороговая неустойчивость). Тяжелые включения в верхней коре дают малую комплексную поправку к инкременту, определяющему конвективную неустойчивость литосферы, чуть увеличивая его действительную и мнимую части, повышая неустойчивость и частоту колебаний. Таким образом, тяжелые включения в верхней коре, которые являются устойчивыми при быстрых процессах, связанных с этими включениями [Биргер, 2023], дают свой вклад в медленно развиваюшуюся конвективную неустойчивость в утолщенных частях литосферы.

About the authors

B.  I. Birger

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

References

Supplementary files