CC BY
22
УДК 551.24
Вестник СПбГУ. Сер. 7, 2005, вып. 4
Ю. И. Кудрявцев
КРИСТАЛЛИЗАЦИОННАЯ МОДЕЛЬ ОКЕАНСКОЙ ЛИТОСФЕРЫ С ЛОКАЛЬНЫМ КОНВЕКТИВНЫМ ПОТОКОМ
В кристаллизационной модели формирования океанской литосферы используют различные решения двухмерного уравнения теплопроводности в области, на контуре которой заданы определенные температурные условия. Широко применяется решение для четверти пространства, на горизонтальной границе которой задается температура морской воды, а на вертикальной дивергентной границе - температура астеносферы [1-3]. Вытекающие из этого решения не ограниченное возрастание мощности литосферы и снижение при этом почти до нуля плотности теплового потока не согласуются с современными данными по тепловым потокам [4].
Численные методы решения уравнения теплопроводности в ограниченной области (например, в прямоугольнике [5]) мало пригодны для теоретических обобщений. Достаточно перспективным оказалось решение в виде ряда для полубесконечной полосы [6-10]. На начальном возрастном интервале оно преобразуется к аналитическому выражению для четверти пространства [9].
Из-за неопределенности выбора ширины полосы и температурных условий на нижней ступенчатой границе возникает неоднозначность решения, которую можно снизить лишь при учете экспериментальных данных. Использование решения уравнения теплопроводности в полубесконечной полосе в качестве аппроксимирующего соотношения, позволяющего объединить теоретические представления и опытные данные, оказалось весьма продуктивным при изучении океанской литосферы [10]. Ниже аппроксимационный подход распространяется на моделирование локального конвективного потока в пространстве между склонами литосферных плит.
Исходные соотношения. Введем координатную систему с началом на дивергентной границе, горизонтальной осью Ох, перпендикулярной этой границе, и вертикальной осью Ог, направленной вниз (рисунок). Предполагается, что двухмерное поле температуры Г удовлетворяет уравнению теплопроводности
V дТ а дх
дх2 д.21
(1)
где а
■ к/[ас
ир) - температуропроводность; к - коэффициент теплопроводности; ср - удельная теплоемкость; ст -
плотность вещества; V - полускорость раз движения плит.
Иногда точность моделирования можно улучшить, увеличивая количество параметров путем задания различных температур на сторонах полосы, отмеченных на рисунке пунктиром. При этом появляются уже две точки разрыва температуры (л: = 0,7 = 0)и(х = 0,2:=21), что несущественно ввиду условного характера нагретой ступенчатой границы, приближенно воспроизводящей глубинный тепловой поток в зоне спрединга. Температура Т\ принята для нижней полупрямой 7, а 72- для вертикального отрезка 2, определяющего ширину полосы 2\. В этой задаче лишь ось Ох соответствует (без учета рельефа дна) реальной границе литосферная плита - морская вода с температурой 7"о.
Для таких граничных условий получено решение уравнения (1) в виде [10]
Т'Т0={Т2 ~ ^"Х-ехр
Я /7 = 1 п
тс п
1 Л
+— \х 1 2 а
. . п% Бш! -2
яп1— г| + (7;-:г0)—.
© Ю. И. Кудрявцев, 2005
д, мВт/м2 250
200
150
100
50
50 75
млн лет
100
25 50
X, км
75
V = о
2, КМ
Результаты моделирования океанской литосферы. Гг оме трические условия: 1 - нижняя горизонтальная граница с температурой Гь 2 - вертикальная граница длиной г\ с температурой Т2. Расчетные графики: 5 - плотность теплового потока через дно океана, 4 - глубина подошвы литосферы, 5 ~ изотерма Г= Тт, 6-8- соответственно плотность теплового потока, глубина подошвы литосферы и
изотерма Т= Тт при 0.
В случае равенства Т2 =Т\ формула (2) сводится к известному решению [6-8].
При спрединге кондуктивным механизмом теплопереноса в направлении оси Ох обычно можно пренебречь. Формула (2) при этом после введения возраста литосферы / = х/У на расстоянии х от дивергентной границы приводит к соотношению
Т-Т0={т2 - Т0)— I —ехр
71 п=I П
\2 си
■(г.-т-^хЫИ
Л Я = 1 П
ехр
(3)
/ \2 Ш
■(ля)- — 2\ _
Дифференцируя зависимость (3) и полагая 2 = 0, на основании закона Фурье находим плотность теплового потока через дно океана:
Т.-То
д -2к ——— £ ехр
+ гкЬ^Л £(-1)"ехр
п=1
-И24
+ к-
В дальнейшем используем условие Т0= 0°С и наименьшую плотность теплового потока при / —> °о
(4)
(5) 33
Ятт=кТт/2* шах' (6)
здесь 2лТ]1ах - предельная мощность океанской литосферы, Т,п - температура ее подошвы и примыкающей к ней астеносферы. не возмущенной влиянием зоны спрединга. Параметры задачи 2и Тх и Т2 удобно заменить двумя другими характеристиками
г = Р - ДГ/Г], (7)
где АТ= Т2- Т\ На основании формул (5) и (7) преобразуются соотношения (4) и (3):
1} _ 2/3 £ (- 1)пехр| - (пп)211| + 1 [<7™ , (8)
о Г /
- («тс)2 — ¡Бт] | -
л/лт ;
?2 £ (-0я Г / -1 О 2 I
X --— ехр1 -\п%) — ¡БШ Л7С-= ] + —=>.
ъ»=\ п |_ I ыах ) 4ах\
(9)
Параметры задачи 2[г Т}, Т2 непосредственного физического смысла не имеют и являются характеристиками фиктивных источников, которые приближенно воспроизводят теплоту, подводимую к склонам литосферных плит реальным конвективным потоком. Благодаря введению в формулы для д и Г надежно определяемой предельной плотности теплового потока дтш количество такого рода параметров сократилось до двух (ти р).
Для каждого возраста необходимо определять глубину подошвы литосферы г, = на которой достигается температура плавления (солидуса) мантийного вещества Тх. С увеличением глубины и давления она растет примерно по линейному закону [11]
Т$=Т50+<Х1$ , (Ю)
в котором а- вертикальный градиент температура солидуса, приведенная к поверхности 2ч-= 0.
Методика определения границы = состоит в следующем [9,10]. Возраст I изменялся с малым шагом А/. Для каждого фиксированного значения / по формуле (9) вычисляют температуру на трех глубинах г и интерполированием находят глубину на которой достигается Т = ТК. При вычислениях Гл по формуле (10) используется глубина установленная для предыдущего возраста, например г + А/.
Расчеты удобно начинать от мощности литосферы, близкой к наибольшей которой соответствует самая большая температура солидуса Т,„. Величину получаем на основании формул (6) и (10):
^тах = * ^о/(<7тт ~ *а) • (П)
Расчетные значения мощности литосферы нуждаются в корректировке ввиду наличия изостатически уравновешенной структуры срединно-океанских хребтов. Поэтому к ним целесообразно добавить слой переменной мощности [10]. Если превышение дна с возрастом I под наиболее древней пониженной частью литосферы составляет И, а превышение фебня срединно-океанского хребта - /гтах, то перепад уровня рельефа равен Д/г= /гтах - И. График Д/г приближенно повторяет конфигурацию глубин подошвы литосферы
~ П7.7, ~ (12)
с коэффициентом
= -0а)/(°1 -О» (13)
причем сг/, а- средние плотности литосферы, астеносферы и воды соответственно.
Аппроксимация экспериментальных данных. Как и в работе [10], учитывались плотности теплового потока, измеренные в северной части Тихого океана и Северо-Западной Атлантики и усредненные по интервалам в 2 млн лет [4]. Данные по возрастному промежутку до 50 млн лет отсутствуют из-за значительной погрешности определения д по температурному градиенту в областях, примыкающих к рифтовым зонам. Использовались также сейсмические определения мощности литосферы [12], воспроизведенные в монографии [13].
При приближении расчетной кривой плотности теплового потока к экспериментальным точкам можно ограничиться экспериментальным значением дт-т и двумя варьируемыми параметрами г и Д При расчетах температуры и глубин литосферы, помимо этих параметров, приходится подбирать теплофизические коэффициенты к и а из интервалов наиболее вероятных их значений для магматических пород [14-16].
По данным различных авторов, вертикальный градиент температуры солидуса находится на интервале 2-^5 °С/км [17-22]. Приведенную же к поверхности дна океана температуру солидуса мантийного вещества лерцо-литового состава можно принять равной Т& = 1100 °С [11].
Подбор параметров т, ¡3 и дшщ. учитывающих влияние внешних по отношению к литосфере эквивалентных источников теплоты, и теплофизических характеристик к, а и а по табличным данным проводился в интерактивном режиме путем приближения расчетных графиков д и 2Л к экспериментальным точкам. Качественный уровень при-
ближения вызван значительной дисперсией экспериментальных данных, которым соответствуют большие стандартные отклонения.
В расчетную формулу (8) для плотности теплового потока входят лишь подбираемые параметры г, /? и причем <7ШШ = 49 мВт/м2. Формула (9) для температуры, помимо их. содержит также теплофизические характеристики. Коэффициенты к и а выбирались несколько меньшими, чем для перидотита, из-за наличия у лерцолитов межзернового пространства, заполненного легкоплавкой компонентой.
Подбором выбраны значения: т = 380-Ю6 лет, р = -0.03, к = 3,2 Вт/(м-°С), а = 9,9-Ю"7 м2/с = 31 м2/год и а = 3,5 °С/км. Таким образом, ввиду малости значения ¡5 другие параметры океанской литосферы практически совпали с установленными ранее [10].
На рисунке представлены расчетные графики 3 и 4 плотности теплового потока и глубины подошвы литосферы. Кривая 4 в интервале от 10 до 50 млн лет выражает почти линейный рост мощности литосферы с возрастом. С последующим увеличением t график z, выполаживается и приближается к предельной мощности = 93 км, определенной по формуле (11).
В интервале изменения / до 20 млн лет тепловое влияние полупрямой У мало, и для определения глубин можно использовать уравнение для координатной четверти [9]
L^lO+ /9) J
где erf"' - функция, обратная интегралу вероятностей. В пределах до 10 млн лет вкладом слагаемого azs можно
пренебречь, и мощность литосферы возрастает здесь пропорционально jt.
Используя предельную мощность zATDax и значение к, на основании формулы (6) находим не возмущенную влиянием зоны спрединга температуру астеносферы Тт = 1430 °С.
Расчетная кривая z4 по формуле (12) на начальном временном интервале пересчитывается в графики перепадов уровня рельефа для северной части Тихого океана и Северо-Западной Атлантики [4] при использовании коэффициентов т, соответственно равных 0,029 и 0,033. По среднему значению т с помощью формулы (13) при <х/= 3,31 г/см3 и - 1,0 г/см3 определена инверсия плотности Лс = а - <га = 72 кг/м3, что совпало с прежней оценкой [10].
Переходный слой. В таблице для возраста t, изменяющегося с шагом 5 млн лет на интервале 5 -^80 млн лет, указаны глубины литосферы zs и производная dzj dt. Последняя делением на полускорость спрединга преобразуется в производную dzj ох, которая имеет смысл тангенса угла наклона у касательной кривой zs к горизонтальной оси Ох:
tg y = (dzs/dt)/V. (14)
Параметры переходного слоя
106 лет Z„ км bs> км dzj dt, 10"2см/год т;,° с dbjdt, 10~2см/год bsdzjdt, м2/год dbjdx!s\ п;/0
5 17,9 9,8 19,8 1160 6,0 19,4 0,50
10 27,0 12,2 16,5 1190 3,6 20,1 0,30
15 34,4 13,4 14,0 1220 1,7 18,8 0,14
20 41,0 13,9 12,6 1240 0,3 17,5 0,03
25 47,0 13,8 11,7 1260. -0,6 16,1 -0,05
30 52,7 13,3 11,1 1280 -1,4 14.8 -0,12
35 58,1 12,5 10,5 1300 -2,0 13,1 -0,17
40 63,2 11,3 9,8 1320 -2,6 11,1 -0,22
45 67,9 9,9 8,9 1340 -3,0 8,8 -0,25
50 72,1 8,4 7,9 1350 -3,0 6,6 -0,25
55 75,8 6,9 6,8 1360' -2,8 4,7 -0,23
60 78,9 5,5 5,6 1370 -2,4 3,1 -0,20
65 81,4 4,4 4,6 1380 -2,0 2,0 -0,17
70 83,5 3,5 3,7 1390 -1,6 1,3 -0,13
75 85,1 2,8 3,0 1395 -1,3 0,8 -0,11
80 86,4 2,3 2,4 1400 -1,0 0,5 -0,08
Углы падения склона плиты исключительно малы. Например, при V — 5 см/год и t = 5-Ю6 лет tgу = 0,040, т. е. у = 2,3°, а при / - 30-106лет у = 1,3°.
Каждый фрагмент литосферной плиты, ограниченный соседними трансформными разломами, характеризуется своей линейной скоростью движения. Следовательно, примыкающие друг к другу фрагменты несколько отличаются углами падения, что способствует образованию и существованию трансформных разломов.
Температура солидуса мантийного вещества совпадает с температурой Тт верхней невозмущенной процессом спрединга астеносферы лишь в случае предельной мощности zsmax. При zs < zíniax, в силу формулы (10), температура Ts < Тт. Таким образом, график глубин литосферы должен проходить выше кривой z* =z*(í) изотермы Т= Тт. Значения zs изотермы 5 рассчитывались по формуле (9) с использованием той же программы, что и глубины zS9 но при условиях равенства а — 0 и замены Ts0 на Тт в формуле (10). Кривые zs и z* ограничивают переходный слой астеносферы.
Расстояние между кривыми 4 и 5 в направлении оси Oz обозначим bs = z* -zs. Касательная к графику 4 образует угол у с осью Ох; тот же угол составляет с осью Oz нормаль к этой касательной. Проекция отрезка bs на эту нормаль равна bsoosy мощности слоя. Ввиду исключительной малости угла у с высокой точностью cosy-1, а поэтому bs фактически является мощностью переходного слоя. В таблице указаны мощность bs и производная dbs! 3/, которая может быть пересчитана в производную dbs /дх = (dbs /dt)/V . На большом
интервале t и х центральной части склона bs« 12-И Зкм.
Литосферные плиты, раздвигаясь от дивергентной границы на величину V, на столько же сближаются с ней благодаря кристаллизации вещества на их склонах. Профиль сечения литосферных плит при этом примерно сохраняется с течением времени. Скорость приращения мощности литосферы равна проекции V на нормаль к кривой zs. Касательная к этой кривой составляет угол у с осью Ох и У, а перпендикуляр к касательной прямой - угол 90° - у Потому приращение мощности склона плиты за единицу времени равно V cos(90° - у) = V sin у ~ Vtgy (ввиду малости /). Заменяя тангенс по формуле (14), получаем
^cos(90o-y) = dz5/d¿. (15)
Таким образом, кристаллизационное увеличение толщины плиты за единицу времени не зависит от скорости спрединга.
Как видно из формулы (15) и данных таблицы, утолщение склона плиты происходит исключительно медленно: dzjdt = 0,2 и 0,03 см/год соответственно для возраста t = 5-Ю6 и 75-106 лет. Также очень мала плотность <j¡ dzjdt массы, кристаллизирующейся в единицу времени на единице площади склона. Для тех же моментов времени она составляет 0,66 и 0,1 г/(см2-год). Выделяющаяся при кристаллизации теплота весьма мала, тем более, что основная тугоплавкая гранульная составляющая астеносферного вещества не претерпевает фазового перехода.
Формулу (15) умножим на a¡(dx/eos у)- a¡Vdt и проинтегрируем по времени. В итоге получим кристаллизирующееся в единицу времени на полосе склона единичной ширины количество вещества, равное <7^51ШХКи пропорциональное скорости спрединга.
Как следует из таблицы, температура кровли переходного слоя Ts монотонно возрастает с увеличением глубины zs и возраста t. При малых t разность температур ST= Тт- Ts изменяется около 200 °С, а при t > 70 млн лет приближается к нулю. Температура кристаллизации некоторых минералов [23] пластичной составляющей межгранульных промежутков может попасть на температурный интервал Тт - ST, Тт. Это приведет к возрастанию гра-
нульной составляющей и увеличению плотности переходного слоя. Прогнозная оценка температуры плавления породы по минеральному составу весьма затруднена из-за существенной зависимости ее от присутствия даже очень малого количества свободной воды.
Другой механизм уплотнения переходного слоя проявляется и при средних, и при малых перепадах температуры 5Т. В процессе кристаллизации на поверхности плиты происходит «усадка» пластической составляющей, дефицит которой компенсируется за счет ее фильтрации из переходного слоя. Из-за оттока этой составляющей тугоплавкие гранулы под действием литостатического давления несколько сближаются, увеличивая плотность переходного слоя.
Избыточная плотность переходного слоя в его верхней части возникает при формировании слоя из тугоплавкой истощенной (деплетированной) мантии. Последняя образуется в условиях пониженного давления под срединным хребтом вследствие выплавления из асте-носферного вещества легкоплавкой и менее плотной составляющей (базальтовой магмы).
Избыточная плотность появляется также благодаря термическому сжатию вещества в соответствии с формулой аа/ЗЗТ, где /?- температурный коэффициент объемного расширения, примерно равный утроенному температурному коэффициенту линейного расширения, приближается к значению 10~5 на 1 °С. Следовательно, эффект термического сжатия заметно проявляется для мантийного вещества, охлажденного по отношению к температуре Тт на величину 8Т, составляющей первые сотни градусов.
В результате действия рассмотренных факторов переходный слой обладает избыточной плотностью о* по отношению к изотермической астеносфере. Она возрастает от а = 0 на границе с температурой Тт до значения, близкого к величине инверсии плотности Асу. Следовательно, склон литосферы, являясь поверхностью фазового перехода, может и не быть заметно выраженной плотностной границей.
Возникновение конвективного потока. Заштрихованный на рисунке переходный слой астеносферы между подошвой литосферы 4 и изотермой 5 обладает повышенной по сравнению с ай плотностью и поэтому оказывается гравитационно-неустойчивым. Он скользит вдоль склона литосферной плиты, образуя нисходящую ветвь конвекции. В месте погружения нисходящего потока формируется зона повышенного давления. Восходящий поток конвекции вдоль осевой области спрединга создается более нагретым астеносферным веществом и является общим для обеих плит.
Обозначим плотность переходного слоя а = аа + ст. Предположим, что избыточная плотность сг* линейно убывает с глубиной г от наибольшего значения на границе 2 = 25 до нуля при 1 — г\\
В качестве оценочной величины а*^ можно принять инверсию плотности Ас при переходе от плиты к астеносфере.
Ввиду малой скорости спрединга структура переходного слоя воспроизводится при сохранении литостатической уравновешенности. Давление в пределах слоя можно получить вычитанием из давления р{г5Гаю) на глубине компенсации г = 25тах давления столба 25ТШХ- 2:
/ \ / * \ ^ тах
Подставляя формулу (] 6), находим
Р = Р(-7-5 гпах ) - ё°а (г5 шах - -
-г
(17)
Определим вертикальную производную давления (17)
При вычислении горизонтальной производной имеем в виду, что р зависит отх посредством глубин г* и
*-г д
ор
сх
1
ох 2
-г. дх
Исключим 2г =;
Ъ:
др
ох
023 1 ~дх 1
1+-
■~15]дЬ
С учетом формулы (18) имеем
_£_
дх
-.1 -- + — 1 + -
5х 21
дЬ с
Эх
(19)
Под действием избыточной плотности в переходном слое возникает движение вещества с малой скоростью и относительно склона плиты. Если пренебречь силой инерции, то уравнение Навье-Стокса [24, 25] выражает равенство нулю всех сил, действующих на единицу объема:
+ = 0.
Здесь представлены плотности силы тяжести , силы давления —§гас1 р и силы вязкого трения (77- коэффициент вязкости). Подставим в это уравнение соотношение (19)
+ + ех+1/У2и = 0. (20)
[ох 2 ^ Ь5 ) ох
Так как ось О2 направлена по вертикали, а плотность в переходном слое а = сга + сг*, то + = 0. Следовательно, уравнение (20) сводится к виду
2-2 АдЬ.
gG | -- + — 1+--I--
I дх 21 Ь, дх
ех и = 0.
(21)
Введем систему координат хЮ*г\ ось О'х' которой на интервале г от 10 до 40 млн лет приблизим к табличным значениям 25 по способу наименьших квадратов. Уравнение этой оси в исходной системе координат имеет вид 2=16 + 1,2*. Здесь 2 выражено в километрах, г - в миллионах лет, а начало отсчета О' находится на оси Ог на глубине ОО' = 16 км. Коэффициент же при г имеет размерность км/(млн лет).
Заменяя г -х/У9 преобразуем уравнение оси О'х':
г^ОО' + хЪуь, (22)
где тангенс угла наклона этой оси составляет
1ё/о-ОД2/К (23)
Даже при весьма малом значении У= 1 см/год угол у0 = 6,8°, т.е. при всех реальных значениях V МОЖНО ПОЛОЖИТЬ СО$у0 - 1, БШУо —ЩУо*
Переход к штрихованным координатам производился по формулам х' = хсоъуо + (2 - 00*)ьт/о, 2' = (2 - 00*)со$уо - хьту0.
Последнее равенство выражается также соотношением
= (г - 00'-щуо)со$го. (24)
Рассмотрим интервал расстояний между точками с возрастом гх — 10 млн лет и и = 40 млн лет и координатами х{ = 100- V км их2 = 400- К км. На этом интервале определяемые формулой (22) глубины практически совпадают с а поэтому из формулы (24) при
со5у0 - 1 следует г'- z - z5, а максимальная мощность слоя г'п1ах = г* = Ъв. Соотношение же (16) для избыточной мощности преобразуется к виду
а (25)
При замене в уравнении (21) орта ех =ех, соьу0 -е2, ьту0 слагаемыми, содержащими произведения $ту0 д2$ /дх и бш у0дЬ5/дх, можно пренебречь как величинами второго порядка малости
8ЬС
ge
1 , , -■ + — 1+-
дх 2 ^ Ъ
S )
сх
eos у0ех> + r¡V2u - 0 . (26)
Отброшенные слагаемые, однако, характеризуют силу движения, направленную противоположно еу и способствующую поэтому динамической стабильности слоя.
Из таблицы видно, что на интервале (хь х2) мощность слоя очень мало изменяется около среднего bs = 13 км, а абсолютная величина производной obs/3c более чем в 4 раза меньше практически постоянных значений azjax -ígу0. В пределах этой плоской части слоя возникает ламинарное (слоистое) течение вдоль подошвы плиты со скоростью и = иеХ'9 причем
величина и от координаты х' почти не зависит: и = С учетом этих условий и формулы (25) из уравнения (26) следует
ц ТТТ + ^^Х sin 70 (1- z'/bs) = 0 . dz "
После первого и второго интегрирования находим
" = j + Qz'+CV
При z' = 0 скорость w = 0, а поэтому = 0. Скорость на другой границе слоя z' = bs мала; положим ее также равной нулю и тогда
Q = siny06,.
Следовательно,
±kz'--z,2 + -í- — . (27)
1 * .
sm7o
Ч
v3 2 6 bs j
Дифференцируя равенство (27), находим координату экстремума 2:' — (х — 1/л/з), а затем максимальную скорость конвекции
"шах = Л/(А2 • (28)
П
Положим в формуле (28) 7] = 1019 Па -с, £ = 10 м/с2, Ъ% - 13-10^м, <7^х = 72 кг/м\ т. е. приравняем а*тх ~ А о, и заменим на основании формулы (23) 8шу0 ~ЩУъ • Перейдем от единиц метр в секунду к сантиметру в год и получим
"шах = 0,30/ К,
причем численный коэффициент имеет размерность (см/год)2. Чрезмерно малое значение ытах, которое более чем на порядок меньше полускорости У > 3 см/год, обусловлено исключительно пологим падением склона плиты и, возможно, завышенной оценкой коэффициент та 77. При rj = 10|8Па- с получилось бы и^ в 10 раз большим.
Интегрированием (27) с учетом равенства (28) получим конвективный поток вещества через площадку высотой bs и единичной ширины:
a¡ \udz = 0965urn&xbsal. о
Астеносферное вещество, кристаллизирующееся в единицу времени на полосе склона единичной ширины, составляет C/zsmzxV. Это значение многократно превосходит величину конвективного потока
^тах Щ0,65bs Ы1ШХ) = 1 1 V/Umах. Таким образом, воспроизведение профиля океанской литосферы при спрединге происходит
не за счет притока астеносферного вещества, а благодаря смещению температурных границ *
zs и zs (без учета очень малого потока пластической составляющей межгранульных промежутков). Указанное условие является необходимой предпосылкой применения кристаллизационной модели.
Сила волочения. Касательная составляющая напряжения трения ps, приложенная к склону плиты, на основании закона Ньютона составляетps = r¡{duidzr) при z' = 0. Подставляя сюда формулу (27), находим плотность силы волочения
Ps =(l/3)g<CAsiny0- (29)
Силу волочения можно оценить, распространяя формулу (29) на любую точку склона плиты. При этом необходимо заменить угол уо на у. Тогда на элементарную площадку 1 -ds, где ds - элемент длины подошвы литосферы, действует сила трения
Psds = ^ga^xbs sin yds. j
Затем заменим ds = dx!cos/, а затем tgу выразим по формуле (14):
л 1 * ь 1 Л
3 et У
Далее положим х = Vt и проинтегрируем по склону плиты
\psds=\ga^x[bs^-dt. (30)
3 ot
Произведения bsdzjdr указаны в предпоследнем столбце таблицы. Численное интегрирование с использованием элементов этого столбца (при t = 0, bs = 0) приводит к соотношениям
\bs(dzs/dt)dt = 8,М08 м2 и \psds = 2,7-108g^x Н/м. Полученный результат легко понять, если рассмотреть фрагмент переходного слоя единичной ширины (в направлении оси спрединга), средней мощности ^ ис максимальной глу-
биной погружения =80-10 м. Преобразовыванием первого интеграла определяется средняя мощность
\bs{czs/3t)dt = \bsdzs =bsís , т. е. bs = 10-103м. С использованием данных обозначений выразим силу волочения (30) этим фрагментом
2 f * Л \Psds = -g
^max
bszs.
2
Величина ^^^х/^)^ представляет перепад давления в пределах г%; это давление, умноженное на Ь5, образует движущую силу, возникающую благодаря средней избыточной
плотности Две трети этой силы компенсируются трением фрагмента о подошву
литосферы, а одна треть - трением на другой границе слоя {г' = Ь5). Таким образом, численная оценка силы волочения составляет
(31)
где Ьл- длина дивергентной границы вдоль оси спрединга.
Сила торможения: Под подошвой практически горизонтальной части литосферной плиты, движущейся с абсолютной скоростью щ, скорость астеносферы направлена вдоль оси Ох( и = иех) и зависит в основном от координаты 2: и = и(г). Величина скорости удовлетворяет однородному уравнению (21) при а = 0, поэтому и = В\2 + Коэффициенты В\ и В2 определяются из граничных условий непрерывности скорости. Под подошвой океанской части плиты и{г = zsmгx) = щ. На границе Леман на которой вязкость астеносферы резко возрастает, можно принять, что и(г = глем) = 0. Отсюда следуют представления скорости
и- uf
и плотности поверхностной силы торможения плиты
Р,=-П-Г-• (32)
^ яем ^ 5 тпах
Под континентами мощность части литосферы, отсчитываемой от уровня дна океана, составляет 7Х7У, и формула (32) преобразуется к виду
- , (33)
2 — 2 л ем ^ кн
где г/'-~ коэффициент вязкости около подошвы континентальной части плиты.
Обозначим площадь плоской части плиты океанской литосферы (без учета ее склона) а континентальной части - 5ХИ. Сила торможения на основании соотношений (32) и (33) составляет
^тр ~ ~U0
(34)
Если сила торможения превалирует над другими силами, то именно она уравновешивает силу волочения: + Ртр = 0. Это условие после подстановки в него соотношений (31) и (34) можно представить в форме равенства для скорости
2,7 4О8 от*
"о = —т-!-^-;-^, (35)
П ] П
и
S ок 2 лел!
где L0K = 80К/Ьд- приведенная длина плоской части океанской плиты.
При оценках скорости примем 72 кг/м\ g = 10 м/с2, г/ = tj' = 1019 Па-с, zJeM =
400-103 м. Для океанской литосферы (SKH = 0) с длиной L0K = 3500-103 м на основании фор-
мулы (35) получим щ = 17- 1(У10 м/с, т. е. и0 = 5,3 см/год. Если же L0K = 2000-103 м, а соотношение площадей SKH/SOK = 1, то скорость w0 ~ 2,3 см/год.
С другими границами плиты также связаны определенные силы. На трансформной границе действует сила взаимодействия контактирующих плит. На конвергентной границе погружающаяся часть плиты благодаря избыточной плотности увлекает за собой, затягивает всю плиту, что характерно для современных литосферных плит [26].
Таким образом, основная движущая сила плиты образуется суммированием сил волочения и затягивания. Предельный возраст океанской плиты ограничен тем, что с увеличением / возрастают площадь подошвы плиты и сила торможения, а сила затягивания снижается из-за формирования слэбов в зоне субдукции. Полускорость раздвижения контактирующих плит V зависит от абсолютных скоростей вида и0 как одной, так и другой плиты и от взаимной ориентировки этих скоростей.
Конвекция в случае плоской границы. Рассмотрим конвекцию в плоском полупространстве с переходным слоем, одной границей которого является ось О'х', а другая образуется откладыванием от нее в перпендикулярном направлении значений bs. Хотя эта задача имеет лишь косвенное отношение к реальной ситуации, на основе ее решения можно подтвердить и дополнить некоторые из полученных ранее результатов.
Положим в уравнении (26) (dzjdx) = tgy0, примем cosy0 - 1 и используем тождество
rot rot u = grad div u - V2u при условии div u = 0 :
n\
di
-rot,
u--rotU = -ЯО
ду/ - 1 5
sm 7o
_! 2 ^
b]
db.
s J
дх
Отбрасывая в левой части равенства производную по у' ввиду предполагаемой слабой зависимости и от этой координаты и подставляя равенство (25), после интегрирования получим
rot = sin7o
>7
, 1 2'
z---
2 b
1
+ — 2
/
z --
v
dbs 1
дх siny0
+ p
где Р - постоянная, не зависящая от г'. Множитель (дЬ5/дх)/бш/о при замене х = У( и использовании (23) вычисляется по формуле (дЬ5/д()/0,\2. Его значения приведены в последнем столбце таблицы. Как видно, слагаемое с данным множителем много меньше первых двух.
Следовательно, вектор гош имеет лишь эту проекцию: го1и = С,, причем плотность вихревых источников равна
Г- ё п S---аг
smy0
, 1
z---
2 b.
1
+ — 2
í
1 z13
db.
1
дх! smy0
■ + P
(36)
Второе уравнение и = 0 соблюдается, поэтому можно ввести векторный потенциал: и = го1А. Ввиду независимости вектора С, от координаты у' его дивергенция равна нулю, а векторный потенциал определяется интегралом
4л: # Я
причем Г - область, занятая вихревыми источниками. Дифференцируя под знаком интеграла по координатам х0, _у0, г0 точки наблюдения, находим
4тг Я3
где Я = ех. (д:0 - *')+ е . (у0 - у')+ е,- (г0 - г').
После интегрирования по у' в пределах -оо, +оо, имеем
здесь 5 - сечение переходного слоя.
К функции и = и(х0,20) можно добавить решение однородных уравнений гоШ = 0 и сНуц = 0, учитывающее влияние отрицательных, «отраженных» в плоскости г' = 0 источников с координатами элементарной площадки (х'9 -2Г). Заменяя также £ по формуле (36), получим
£ * ■ гг и = - —^ТШХ ЯП у о Н 2щ
, 1
г---+ —
2 Ьг 2
, I/3
дЪ.
1
дУ $ш у0
■ + Р
е., х
-У)2+(20-/)2 (Х0-У)2+(20+/)2 Таким образом, проекции вектора скорости составляют
(38)
"У = 8И1УоЯ
271?/
2
1
3 Ъ]
дЬ.
1
г0-2!
2Л +2
оУ эт у0 М №,
- + />
(39)
2щ
I/2 1
2 Ь5 2
2
12^ 3 Ь2 .
1
5У sin у0 <ЫсЫ.
■ + Р
(40)
>0 -У)2 +(2о ~/)2 (х0 -У)2 + + _
Благодаря добавлению предложенного решения однородных уравнений формула (40) удовлетворяет граничному чг-(г0 = 0) = 0 . Второе граничное условие их-(г0 = о)= 0 следует из равенства (ЗУ), если принять, что коэффициент Р зависит только от координаты х0 и равен
(2 л ( 1 _.3 >
Р = ~\\
5*
¿<ы<ы
,12" 1
2---+ —
.2 Ъ. 2
дЪ„
1
5 )
дх этуо
(х0 - х'У
+ 2'
-с1х (12 IВ ,
(41)
о -х)
Рассмотрим бесконечный слой постоянной мощности Ъ3\ (дЬ5/дх) = 0. В этом случае в формулах (41) можно проинтегрировать по х' в пределах -оо, +оо, а затем по г' и получить
Р = - ЪД. (42)
В случае бесконечного слоя из формулы (40) следует = 0. Интегрирование же в формуле (39) пох' приводит к соотношению
2г\ —
\
20
(х0-У)2+(20-2')2 (Ло-У)2+(20+2')2_
м 20 - ^ ах - -п.-
ч-А
-п .
(43)
Итак, величина (43) равна нулю при г' < ¿0 и -2п при 2' > Умножим на этот коэффициент результат интегрирования (39) по в пределах (20, Ь5) с учетом условия (дЬ5/дх) = 0 и формулы (42). В итоге приходим к соотношению (29)
Иу *Ы70\\20 С
// ^2 6 3
хотя условие = 65) = 0 при выводе не использовалось, а получилось в качестве след-
ствия.
Учитывая исключительно плавное изменение мощности Ь5, в качестве приближенного решения можно использовать проекции скорости (39) и (40) при выборе в качестве коэффициента Р его значение (42) для бесконечного слоя.
К подошве переходного слоя с внешней стороны > Ь5) на интервале, где Ь5 практически постоянно, примыкает застойная зона с нулевой скоростью, вытекающей из формул (39) и (43).
Из соотношений (39), (40) и (42) можно получить предельные формулы для скорости при > 2,5Ь£. В этом случае значения ых» получаются отрицательными; проекция же иг> < 0 на интервале малых х0 и и2> >0 при относительно больших х0. Такие особенности поля скоростей характерны для обратной, замыкающей ветви конвекции.
Затухание спрединга. Теплопереносом за счет относительного движения плит можно пренебречь, если полускорость спрединга У < 0,05 см/год. В этом случае в температурной зависимости (2) можно опустить слагаемое У/2а и заменить параметры Ти Т2 по формулам (5) и (7) и приравнять Т0 - 0 °С:
Т = Ч г
у[ат
|-(1 + £)£-ехр
ИТ /7=1«
/
-пп
2 » (-IV
П п= 1 П
-пп-
ч
у[ат
БШ
пп
у[ат
(44)
пп
4ат) \ у[ат) л[ат \
Из формулы (44) следует представление плотности теплового потока
/ \ V \
2(1 + /?)1ехр
/7 = 1
-пп
ч
у[ат
-2/?1(-1)" ехр
/7=1
— пп
■у/ат
(45)
По формулам (44) и (45) рассчитаны представленные на рисунке кривые 6-8. Точка на склоне кривой 7, заглубленная на половину 25тах, находится от начала координат лишь на расстоянии хо.5= 23 км. Это примерно в 50 раз меньше, чем в случае полускорости У = 5 см/год. Как следствие имеет место крутое падение склона литосферы - примерно под углом 60°.
Из-за чрезмерного сближения переходных слоев вертикальный поток в верхней его части блокируется, и пространство между краями плит заполняется продуктом кристаллизации мантийного вещества. Образование застойной зоны между границами переходных слоев может быть начинается уже при К, составляющей первые десятые доли 1 см/год.
При затухании спрединга самая интенсивная кристаллизация начинается под частью дивергентной границы, наиболее приближенной к эйлерову полюсу вращения, где линейные скорости раздвижения малы. Конвективное движение здесь возможно только вниз при условии, что на некотором расстоянии просвет между плитами заполнен вертикальным потоком поднимающегося снизу нагретого мантийного вещества. Вверху восходящего потока может сформироваться вулканический центр (горячая точка). По другую сторону его возникает вторая нисходящая ветвь потока. Таким образом, появляется продольная конвективная ячейка, ориентированная вдоль дивергентной границы. Характерной особенностью этой 44
ячейки являются относительно высокие скорость и температура восходящего потока из-за малости площади его сечения по сравнению с параметрами нисходящих ветвей. Очевидна также недолговечность ее существования вследствие обуживания нисходящих ветвей и сокращения протяженности дивергентной границы. Появление следующей горячей точки возможно лишь на некотором расстоянии от исчезающего вулканического центра.
Постепенное затухание спрединга от периферийной части срединно-океанского хребта может отмечаться цепочкой вулканических центров, как, например, в случае Гавайского и Императорского хребтов.
Наибольшее количество современных горячих точек приурочено к дивергентным границам с относительно малой полускоростью спрединга и частично к трансформным разломам, а также к местам пересечения этих границ [26]. Вероятно, их появление обусловлено преобразованием фрагментов срединно-океанских хребтов в вулканические центры в процессе
перехода к механизму теплопереноса посредством продольных конвективных ячеек.
* * *
Используемое при моделировании океанской литосферы решение двухмерного уравнения теплопроводности в прямоугольной полубесконечной полосе предусматривает предварительное задание как ширины полосы, так и температурных условий на ее нижней ступенчатой границе. Эти параметры задачи конкретного физического содержания не имеют, а являются характеристиками фиктивных источников, приближенно воспроизводящих тепловой режим в зоне спрединга. Потому решение задачи использовалось лишь для аппроксимации экспериментальных определений плотности теплового потока через дно океана и глубин литосферы, полученных сейсмическим методом.
Аппроксимационный подход обеспечил сглаживание и объединение данных различных геофизических методов, обремененных значительными погрешностями, а также их согласованность с теоретическими представлениями о кристаллизации литосферы из астеносфер-ной мантии. Полученные расчетные графики мощности литосферы и плотности теплового потока достаточно хорошо согласуются с предшествующими оценками.
Разработанный метод моделирования был применен и для обоснования условий возникновения локального конвективного потока.
Вследствие зависимости температуры солидуса от глубины на склоне литосферной плиты формируется переходный астеносферный слой, ограниченный снизу изотермической поверхностью, которая касается подошвы наиболее древней океанской литосферы. Температура на этой поверхности совпадает с температурой верхней астеносферы. Рассмотрены процессы возникновения в переходном слое избыточной плотности, которая может возрастать от нулевого значения на изотермической поверхности почти до величины инверсии плотности вблизи склона литосферы. Для гравитационно-неустойчивого переходного слоя на основании приближенного решения уравнения Навье-Стокса и строгого решения для «плоской» наклонной границы установлены основные особенности возникающего локального конвективного потока. Оценена сила волочения, приложенная к склону плиты, в сопоставлении с силой торможения, действующей на плоских частях ее подошвы. Высокотемпературное состояние астеносферного вещества между краями литосферных плит является необходимым условием, способствующим образованию срединно-океанского хребта и существованию самого спрединга. Ввиду интенсивного теплового потока через склоны хребта поддержание температуры вблизи значения Тт возможно лишь при конвективном теплопереносе. Нисходящие ветви химико-плоскостной конвекции на склонах литосферных плит возникают в основном благодаря избыточной плотности деплетированной мантии. Последняя формируется под срединно-океанским хребтом при декомпрессионном плавлении астеносферного вещества в итоге отделения его легкоплавкой компоненты (ба-
зальта). В средней же части межплитового пространства появляется замыкающий горячий мантийный поток.
При очень малой скорости раздвижения плит их переходные слои чрезмерно сближаются, что приводит к блокированию восходящего потока. Некоторое время режим конвекции сохраняется в виде продольных ячеек, сопряженных с вулканическими центрами.
Summary
Kudrjavtsev Yu. /. Crystallization model of the ocean lithosphere with the local convectional flow. The parameters of a crystallization model of the ocean lithosphere are obtained by approximation of geophysical data using the solution of a two-dimensional equation of heat conductivity in the semi-infinite layer. As follows from this model, a gravitationally unstable transitive layer is formed on the slope of every lithospheric slab and there arises a local convectional flow connected with it. The value of the plate dragging force in comparison with one of braking is estimated.
Литература
1. Сорохтин О. Г. Образование литосферных плит и природа срединночжеанских хребтов // Океанология. Геофизика океана. Т. 2: Геодинамика / Под ред. А. С. Монина, О. Г. Сорохтина М, 1979. 2. Городницкий А. М, Сорохтин О. Г. Карта расчетных значений теплового потока через дно океана // Проблемы теоретической геодинамики и тектоника литосферных плит. М., 1981. 3. ТеркотД., Шуберт Дж. Геодинамика: Геологические приложения физики сплошных сред / Пер. с англ.; Под ред. В. Н. Жаркова. М., 1985. Т. 1.4. Stein С. A., Stein S. Comparison of plate and asthenospheric flow models for the thermal evolution of oceanic lithosphere // Geophys. Res. Lett. 1994. Vol. 21, N 8. 5. Kusznir N. J. Thermal evolution of the oceanic crust; its dependence on spreading rate and effect on crustal structure // Geophys. J. R. astr. Soc. 1980. Vol. 61. 6. Mc Kenzie D. P. Some remarks on heat flow and gravity anomalies // J. Geophys. Res. 1967. Vol. 72. 7. Sclater J. G. New perspectives in terrestrial heat flow // Tectonophysiks. 1972. Vol. 13. 8. Sclater J. G., Anderson R. N., Bell M. L. Elevation of ridges and evolution of the central eastern Pacific // J. Geophys. Res. 1971. Vol. 76. 9. Кудрявцев Ю. Я Связь теплового потока через дно океана с мощностью литосферы // Рос. геофиз. журн. 1998. № 11-12.10. Кудрявцев Ю. И. Моделирование океанской литосферы // Рос. геофиз. журн. 2000. № 17-18. 11. Городницкий А. М. Строение океанской литосферы и формирование подводных гор. М., 1985. 12. Watts А. ВBodine J. Я, Steckler М. S. Observations of flexure and the state of stress in the oceanic lithosphere // J. Geophys. Res. 1980. Vol. 85. 13. Kearey P., Vine F. J. Global tectonics. Tulane, 1996. 14. Моисеенко У. И. Теплофи-зические свойства горных пород и глубинные температуры // Физические процессы горного производства. 1982. № 12. 15. Ыоисеенко У. Я, Смыслов А. А. Температура земных недр. Л., 1986. 16. Моисеенко У. Я., Негров О. Б. Теплофизические параметры горных пород// Петрофизика: Справочник. Кн. 1. Горные породы и полезные ископаемые / Под ред. Н. Б. Дортман. М., 1992. 17. Йодер Г. С. Образование базальтовой магмы / Пер. с англ.; Под ред.
A. А. Кадика. М, 1979. 18. Йодер Г. С., Тилли К Э. Происхождение базальтовых магм / Пер. с англ.; Под ред.
B. С. Соболева. М., 1965. 19. ho КKennedy G. С. Melting and phase relations in a natural peridotite to 40 kilobars // Amer. J. Sci. 1967. Vol. 265. 20. Kushiro /., Syono У., Akimoto S. Melting of a peridotite nodule at high pressures and high water pressures I I J. Geophys. Res. 1968. Vol. 73, N 18. 21. Ringwood A. E. Composition and evolution of the upper mantle // Тис ЕашГь crust and upper mantle. Washington, 1969 (Geophys. Monogr. Amer. Geophys. Union. N 13). 22. Pot-lack Я. /V., Chapman D. 5. On the regional variation of heat flow, geotherms and lithospheric thickness // Tectonophysics. 1977. Vol. 38. 23. Robie R. AHemingway В. Thermodynamic properties of minerals and reladed substances at 298, 15 К and 1 Bar pressure and at higher temperatures // U.S. Geol. Surv. Bull. 1995. N 2131. 24. Лойцянский 77. Г. Механика жидкости и газа. М.? 1987. 25. Нигматулин Р. Я Основы механики гетерогенных сред. М., 1973. 26. Апло-нов С. В. Геодинамика. СПб., 2001.
Статья поступила в редакцию 12 мая 2005 г.
