Physical and mathematical modeling dynamic pressure forces and viscous stresses in lithospheric subduction zones

Full Text

Введение

Целью данной работы является изучение роли сил динамического негидростатического давления и вязких напряжений, действующих в зонах литосферной субдукции посредством физико-математического моделирования. Приведем некоторый обзор опубликованной литературы по данному направлению исследований. Проблеме физико-математического моделирования сил динамического давления и вязких напряжений в зонах литосферной субдукции посвящено достаточно много работ как отечественных [1-5, 8], так и зарубежных [9-14] ученых.

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

Результаты и обсуждение

Рассмотрим модель конвекции в прямоугольной ячейке 0 < x < L, 0 < z < d, с началом координат в основании верхней мантии на глубине d, вертикальной осью z, направленной вверх, и горизонтальной осью x вдоль основания верхней мантии. Ячейка заполнена однородной жидкостью плотностью ρ с коэффициентом вязкости η, горизонтальные границы z = 0 и z = d изотермичны с температурами T(z = 0) = T₀ и T(z = d) = T₁, а вертикальные границы x = 0 и x = L считаются адиабатичными, на которых ∂T/∂x = 0. Безразмерные линеаризованные уравнения, определяющие возмущения термомеханического состояния среды в ячейке при бесконечном числе Прандтля в приближении Буссинеска. имеют вид уравнений (7.3.11) – (7.3.14) в [14], которые в обозначениях настоящей работы могут быть записаны как

$0=-{\partial }_{x}p+{\partial }_x{\text{τ}}_{xx}+{\partial }_z{\text{τ}}_{xz}$ (1)

$0=Ra\times \text{θ}-{\partial }_{z}p+{\partial }_x{\text{τ}}_{xz}+{\partial }_z{\text{τ}}_{zz}$ (2)

$0={\partial }_x{v}_x+{\partial }_z{v}_z$ (3)

${\partial }_t\text{θ}=-v_z{\partial }_{z}T+\text{χ}\Delta \text{θ}$ (4)

где сохранен член ${\partial }_t\text{θ}$, описывающий нестационарную задачу, и знак при $Ra\times \text{θ}$ изменен, так как ось z направлена вверх. Уравнения (1)–(4) есть соответственно $x$ - и $z$-компоненты уравнения движения, уравнение неразрывности и уравнение теплопереноса, в которых $x$ и $z$ – декартовы координаты, p – динамическое (негидростатическое) давление, ${\text{τ}}_{ik}$ – тензор вязких напряжений,  $\text{ρ}$– плотность,  $g$– ускорение силы тяжести,  $c_p$– удельная теплоемкость при постоянном давлении,  $T$– абсолютная температура,  $\text{κ}$– коэффициент теплопроводности,  $\Delta$– оператор Лапласа, а символ $\partial$ с индексом обозначает частную производную по координатам x, z и времени $t$.

В (1)–(4)  = (κ / ρc_p) – коэффициент температуропроводности, и для приведения уравнений к безразмерной форме в качестве новых единиц измерения координат $x$ и $z$ выбраны мощность слоя $d$, скорость – величина $\text{χ}/d$, время – $d^2/\text{χ}$, температура $T$ и ее возмущения $\text{θ}$ – характерный перепад температуры  $\text{δ}T$= (T₀ – T₁) > 0, напряжения и давление – величина ($\text{ηχ}/d^2$ ). В (2) безразмерное число Рэлея есть

$Ra=\frac{\text{ρα}g{d}^3\text{δ}T}{\text{ηχ}}$> 0, (5)

где α – коэффициент теплового расширения.

Рассматривая двумерную конвекцию в плоском горизонтальном слое $0\le z\le d$ первоначально покоящейся жидкости, в которой имеется вертикальный градиент температуры $T_z=(T_1-T_0)/d<0$, с невозмущенным термомеханическим состоянием покоя с постоянным вертикальным градиентом температуры $T_z=(T_1-T_0)\times d^{-1}$ и чисто кондуктивным переносом тепла, можно искать решение уравнений (1)–(4) с экспоненциальной зависимостью от времени по закону exp(γt).

При условии свободных непроницаемых изотермических горизонтальных и адиабатических вертикальных границ ищем решение (1)–(4), при постоянных (безразмерных) ${\partial }_{z}T$ < 0 и $\text{χ}$ в виде:

$v_x=A\sin kx\cos \pi z$, $v_z=B\cos kx\sin \pi z$, $\theta =C\cos kx\sin \pi z$, $p=D\cos kx\cos \pi z$, ${\partial }_{x}p=-Dk\sin kx\cos \pi z$, ${\partial }_{z}p=-D\pi \cos kx\sin \pi z$, ${\partial }_z{v}_x=-A\pi \sin kx\sin \pi z$, ${\partial }_x{v}_x=Ak\cos kx\cos \pi z$, ${\partial }_x{v}_z=-Bk\sin kx\sin \pi z$, ${\partial }_z{v}_z=B\pi \cos kx\cos \pi z$, ${\tau }_{xx}=2\eta Ak\cos kx\cos \pi z$, ${\tau }_{xz}=-\eta (A\pi +Bk)\sin kx\sin \pi z$, ${\tau }_{zz}=2\eta B\pi \cos kx\cos \pi z$, ${\partial }_x{\tau }_{xx}=-2A\eta k^2\sin kx\cos \pi z$, ${\partial }_z{\tau }_{zz}=-2B\eta {\pi }^2\cos kx\sin \pi z$, ${\partial }_x{\tau }_{xz}=-\eta (A\pi +Bk)k\cos kx\sin \pi z$, ${\partial }_z{\tau }_{xz}=-\eta (A\pi +Bk)\pi \sin kx\cos \pi z$,  (6)

где все не зависящие от координат величины A, B, C, D в (6) зависят от времени t по экспоненциальному закону $\exp \left(\text{γ}t\right)$, а k = πd×L^–1 есть (безразмерное) волновое число. Подставляя (6) в (1)–(4), находим для безразмерного инкремента $\text{γ}$

$\text{γ}=-\frac{Ra{k}^2{T}_z}{\text{η}{({\text{π}}^2+k^2)}^2}-\text{χ}({\text{π}}^2+k^2)$ (7)

Условие возникновения конвекции $\text{γ}$= 0 дает $Ra(\text{γ}=0)=-{({\text{π}}^2+k^2)}^3/k^2{T}_z$. Эта величина достигает минимума при $k=\text{π}/\sqrt{2}$, и при $T_z$=1, $R{a}_{\min }=\frac{27}{4}{\pi }^4$ ≈ 658. Если конвекция происходит в горизонтальном слое неограниченной длины, то возникают ячейки с пространственным периодом d $\sqrt{2}$. В случае если $T_z$, $\text{η}$, $\text{χ}$ переменны, можно для оценки инкремента конвективной неустойчивости воспользоваться формулой (7), подставив в нее средние значения $T_z$, $\text{η}$, $\text{χ}$.

Рассмотрим подробнее вывод формулы (7) из уравнений (1)–(4). Пусть начальное возмущение температуры задается в (6) как $\text{θ}=C\cos kx\sin \text{π}z$ с C>0. Это означает, что возмущение температуры в левой части ячейки положительно, а в правой части – отрицательно, т. е. в левой части ячейки вещество всплывает, а в правой – опускается, и, следовательно, конвективное движение жидкости происходит по часовой стрелке. Из (3) следует B = –Ak/π. Так как Δθ = –(π² + k²) C·cos kx·sin πz, ∂_tθ = γ C cos kx sin πz, то из (4)

C =$\frac{(k/\text{π})A{T}_z}{\text{γ}+\text{χ}({\text{π}}^2+k^2)}$ .                                                                

где при C > 0 и T_z < 0 должно быть A < 0. Подставляя выражения (6) в уравнения (1) и (2), вычитая одно из уравнений из другого и сокращая полученный результат на A, приходим к формуле (7). Из (1) находим D=Aη(k² + π²)/k, где k = π/L, т. е. D < 0.

Согласно выражениям в верхней строке (6)

$v_x=A\sin kx\cos \pi z$, $v_z=B\cos kx\sin \pi z$, $\theta =C\cos kx\sin \pi z$, $p=D\cos kx\cos \pi z$,              

при C > 0, A < 0, B > 0, D <0 компоненты скорости v_x и v_z соответствуют движению жидкости по часовой стрелке, т. е. всплыванию жидкости в левой части ячейки и опусканию жидкости в правой части ячейки. На верхней границе ячейки (при z = 1) возмущение динамического давления `p = – D·coskx= – Aη(k² + π²)/k·coskx. Сила давления действующая изнутри жидкости на верхнюю границу положительна в левой части ячейки (т.е. «подпирает» границу снизу) и отрицательна в правой части ячейки (т. е. «засасывает» границу вниз). Сравним силу негидростатического давления на верхней границе ячейки с вертикальной силой вязких напряжений, действующей со стороны жидкости на верхнюю границу ячейки.

Нормальная компонента тензора вязких напряжений

τ_zz = – 2Dk²/(π² + k²)×coskxcosπz = –2ηAk×coskxcosπz,                                   

и на верхней границе z = 1, cos πz = –1, τ_zz = 2ηAk×cos kx, т.е. при A < 0 оказывается, что τ_zz на верхней границе отрицательна в левой части ячейки и положительна в правой части ячейки. Так как сила, действующая со стороны жидкости на единицу обтекаемой поверхности границы с внешней нормалью n_i, равна [5, формула (15.14), в которой изменен знак нормали n_i]:

f_i = – pn_i + τ_ikn_k,                                                                   

то вертикальная сила, соответствующая вязким напряжениям, равна –τ_zz, так как направленная внутрь жидкости нормаль на верхней границе n_z = –1. Следовательно, вязкие напряжения в левой части ячейки действуют на верхнюю границу, как сила «подпора» снизу (в положительном направлении оси z), а в правой части ячейки – как сила «подсоса» вниз. Сравним конвективные силы вязких напряжений и негидростатического давления, действующие на верхнюю границу ячейки. Отношение этих сил

$\frac{f_{visco}}{f_{press}}=\frac{2k^2}{k^2+{\pi }^2}$, (8)

откуда видно, что при k = π (т. е. в изометрической ячейке с отношением сторон 1 : 1, т. е. при L = d) эти силы равны между собой. В случае, например, удлиненной ячейки, для которой k < π, на верхней границе ячейки преобладает сила возмущенного негидростатического давления. В сильно удлиненной ячейке, для которой k << π, действие вязких напряжений на верхней границе пренебрежимо мало по сравнению с действием сил возмущенного динамического давления. Соотношение (8) справедливо не только на поверхности ячейки, но и во всем ее объеме. Следует отметить, что силы негидростатического давления и вязкие напряжения в земных недрах действуют, так сказать, «однонаправленно», т. е. в одну сторону, и этот вывод не связан именно с конвективной природой движения, а приложим к движениям различной природы.

Поскольку основные результаты работы были получены на модельных примерах, то, по-видимому, необходимо привести ссылки на опубликованные нами работы [1-4] и привести графический пример о сравнении результатов наших модельных расчетов с опубликованными результатами практических геолого-геофизических исследований [6].

 

Рисунок 1 – Глубинный геолого-геофизический разрез квазистационарного распределения безразмерной функции тока в мантийном субдукционном клине

 

На рисунке 1 изображен глубинный геолого-геофизический разрез квазистационарного распределения безразмерной функции тока в мантийном субдукционном клине, расположенном на акватории Черноморской литосферной плиты, погружающейся под Русскую континентальную литосферную плиту, с учетом эффектов диссипативного нагрева и конвективной неустойчивости для неньютоновской реологии мантии и концентрации воды C_w = 3×10^-1 весовых %% в горных породах мантии. Вертикальная стрелка (с) показывает восходящий конвективный поток. В верхней части рисунка приведен реальный геолого-геофизический разрез мантии [6] по профилю от полуострова Крым по направлению на северо-запад, к Русской литосферной плите. Условные обозначения: 1 – осадочные породы; 2 – дислоцированные породы молодого фундамента; 3 – поверхность молодого фундамента; 4 – поверхность дорифейского фундамента (V=6.2-6.5 км/с); 5 – «гранитный» слой; 6 – породы основного состава (V=7.0 км/с); 7 – породы коро-мантийного слоя; 8 – поверхность коро-мантийного слоя (V=7.5-7.6 км/с); 9 – граница Мохоровичича; 10 – сейсмические горизонты верхней мантии; 11 – слои с пониженной скоростью в верхней мантии по данным ГСЗ; 12 – поверхность астеносферного слоя по геотермическим данным; 13 – поверхность астеносферного слоя по данным МТЗ; 14 – глубинные тектонические разломы; 15 – очаги землетрясений.

Из рисунка 1 видно, что в расчетной модели субдукционной зоны (нижняя часть рисунка) с погружением Черноморской литосферной плиты под Русскую был правильно рассчитан угол наклона этой плиты (27^о), который соответствует наклону верхней кромки этой литосферной плиты, фиксируемому на реальном геолого-геофизическом разрезе (верхняя часть рисунка). Кроме того, из рис. 1 видно, что точно рассчитано местоположение конвективной вихревой структуры квазистационарного распределения безразмерной функции тока, связанной с формированием восходящего конвективного потока (вертикальная стрелка – с). Над этой конвективной зоной, в верхней части рисунка, на реальном геолого-геофизическом разрезе наблюдается зона плавления горных пород мантии, которую мы ассоциируем с восходящим мантийным термическим диапиром.

В качестве примера рассмотрим тектонически активную окрестность зоны субдукции, и качественно сравним силу динамического (негидростатического) давления и вязкие напряжения, действующие на субдуцирующий литосферный блок и подошву динамической топографии в этой области. Так как крупномасштабные циркуляционные движения под движущейся и субдуцирующей Черноморской субокеанической литосферной плитой и Русской континентальной плитой, с которой сталкивается субокеаническая плита, происходят внутри тех частей верхней мантии, которые сильно удлиненны в горизонтальной направлении, то в рамках рассмотренной конвективной модели эти циркуляционные движения характеризуются условием d << L, или, в безразмерном виде, k << π в формуле (8). Учитывая изложенное выше и ранее опубликованные работы [1-4] это может означать, что в окрестности зон субдукции на субдуцирующие литосферные блоки и «подошву» настилающей литосферы действуют преимущественно силы динамического давления, а вязкие напряжения несущественны. Это условие (хотя и без всякого обоснования), используется в [7] в параграфе § 6.11 об угле субдукции. Можно, также, отметить, что расчет динамической топографии в [10] выполняется путем вычисления упругого изгиба верхней части коры, подпираемой снизу динамическим давлением в вязком течении, происходящим в слое нижней коры, причем этот достаточно тонкий слой очень сильно вытянут в горизонтальном направлении. При этом также не учитываются вязкие напряжения, действием которых авторы пренебрегают, не приводя каких-либо обоснований.

Однако в окрестности зон субдукции микро-плит, например, Черноморской [1], Амурской [2], Адриатической [3] и Палео-Уральской [4], движения в астеносфере оказываются примерно изометричны, и роли сил динамического давления и вязких напряжений оказываются сравнимы между собой. Этим, возможно, объясняется то, что субдукция микро-плит и малых плит происходит под достаточно малыми углами к горизонту, так как субдуцирующий блок поддерживается снизу и «подсасывается» сверху не только силами динамического давления, но и сравнимыми силами вязких напряжений. Напротив, динамическая топография, формирующаяся над восходящими термическими диапирами, относительно узкими в горизонтальном направлении, связана с астеносферными потоками, в которых k>>π в формуле (8). Следовательно, динамическая топография над диапирами оказывается обусловлена преимущественно вязкими напряжениями.

Сравнивая с другими [7-14] (и огромным количеством ссылок по ним) полученные в данной работе результаты, следует отметить, что многочисленные термомеханические модели мантии в зонах субдукции показали, что конвекция в виде поперечных валообразных вихрей и формирующихся термических диапирах никогда ранее не анализировалась, поскольку модели с чрезвычайно малым углом субдукции и достаточно большой скоростью субдукции ранее никем не исследовались.

Заключение и выводы

Показано, что относительная роль сил динамического давления и вязких напряжений, действующих в областях верхней мантии, характеризуемых астеносферными течениями в окрестности активных тектонических зон, зависит от соотношения горизонтального и вертикального масштабов течений в астеносфере. Если горизонтальный масштаб движений значительно превышает их вертикальный масштаб, то роль сил динамического давления существенно преобладает над ролью вязких напряжений, и последними можно пренебречь. Так, в окрестности зон субдукции протяженных плит при вычислении динамической топографии и угла субдукции можно пренебречь вязкими напряжениями и учитывать только силы динамического давления. В зонах субдукции микро-плит следует учитывать как динамическое давление, так и вязкие напряжения, роли которых сравнимы, чем, по-видимому, объясняются малые углы субдукции микро-плит. Динамическая топография над термическими диапирами, напротив, обязана своим происхождением преимущественно вязким напряжениям. Новизна настоящей работы заключается в том, что в окрестности зон субдукции литосферных плит движения в астеносфере оказываются примерно изометричны, и роли сил динамического давления и вязких напряжений оказываются сравнимы между собой. Кроме того, новизна результатов данной работы может быть обусловлена еще и тем, что динамическая топография, формирующаяся над восходящими термическими диапирами, относительно узкими в горизонтальном направлении, связана с астеносферными потоками, в которых k >> π в формуле (8). Полученные в настоящей работе результаты могут быть использованы при решении некоторых задач теории упругости литосферных плит.

About the authors

Sergey V. Gavrilov

Schmidt Institute of Earth Physics of the Russian Academy of Sciences

Author for correspondence.
Email: gavrilov@ifz.ru

Doctor of Physical and Mathematical Sciences, Main Researcher of the Laboratory 102

Russian Federation, Moscow

Andrey L. Kharitonov

Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation of the Russian Academy of Sciences

Email: gavrilov@ifz.ru

Candidate of Physical and Mathematical Sciences, Leading Researcher of the Laboratory of the Main Magnetic Fields of the Earth

Russian Federation, Moscow

References

Supplementary files