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},$ вызва�