CC BY
79
Вычислительные технологии
Том 2, № 5, 1997
ЧИСЛЕННАЯ МОДЕЛЬ ТЕПЛОВОЙ КОНВЕКЦИИ В ВЕРХНЕЙ МАНТИИ ЗЕМЛИ ПОД ЛИТОСФЕРОЙ КОНТИНЕНТОВ*
Е. В. РычковА Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: [email protected]
С. А. Тычков Объединенный институт геологии, геофизики и минералогии СО РАН, Новосибирск, Россия e-mail: [email protected]
The spatial-temporal evolution of thermal convection in the continental upper mantle is numerically investigated by the predictor-corrector method with majorant splitting scheme on the predictor (for the equation of temperature transport) and the iterating method of stabilizing corrections (to find the stream function from the equation of the fourth order with variable coefficients). Lithosphere is presented by rigid conductive layer with constant thickness (120 km), the viscosity of the upper mantle depends on the depth. The non-slip conditions are given at the bottom of the lithosphere, and convection is excited by constant heat flux from the lower boundary. At this boundary conditions the structure of the convection in 7000 M years attains to the steady state as one long cell.
1. Введение
По современным представлениям, тепловая конвекция в мантии Земли, эффективно выносящая радиогенное тепло из недр, принимающая непосредственное участие в движении литосферных плит и формировании тектонических режимов территорий, имеет сложную трехмерную структуру, которая определяется строением оболочки, распределением физико-химических характеристик с глубиной и фазовыми переходами вещества.
Главной особенностью структуры принято считать слоистость конвекции, то есть существование замкнутых изолированных ячеек отдельно в верхней и нижней мантии. Эта
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-05-65970.
© Е. В. Рычкова, С. А. Тычков, 1997.
изоляция, однако, может локально нарушаться, что обусловлено периодически возникающей неустойчивостью теплового граничного слоя в районе эндотермического фазового перехода на границе верхней и нижней мантии [1]. Ключевым вопросом при изучении конвекции является выяснение причин и условий, определяющих пространственно-временную эволюцию структуры конвекции, поскольку именно эта характеристика во многом определяет кинематику литосферных плит и геологическую историю развития континентальных областей. Принципиальными понятиями эволюции конвекции принято считать стационарность и степень устойчивости процесса. Под стационарностью подразумевается режим течения, при котором в каждой точке пространства, занятого жидкостью, скорость движения не меняется во времени. Устойчивость процесса характеризует его отношение к малым возмущениям: если, раз возникнув, возмущения затухают со временем, то данный режим конвекции устойчив [2].
Режимы тепловой конвекции в мантии Земли интенсивно изучаются уже более 30 лет. К настоящему времени разработано несколько классификаций режимов конвекции в зависимости от числа Рэлея (И,а) [3-7]. Хотя авторы этих работ при создании классификаций исходили из моделей с различными граничными и начальными условиями, реологией, способами возбуждения течений, геометрией расчетной области, а также привлекали данные лабораторных экспериментов, общим для большинства построений явилось утверждение о том, что переход от стационарных течений к нестационарным осуществляется в пределах 104 < И,а < 105. Необходимо отметить, что авторы классификаций использовали в своих построениях различные характеристики течения, с помощью которых отслеживалась эволюция, и резкое их изменение трактовалось ими как смена режима. Среди рассматриваемых характеристик наиболее часто употребляемыми являются средние по объему температура или скорость движения вещества конвектирующей области [6, 7], зависимость числа Нуссельта от времени и эволюция температуры в средней части восходящего конвективного потока [5], количество и размер ячеек, а также эволюция траекторий состояний в фазовом пространстве [5, 8]. При И,а > 105 и дальнейшем увеличении числа Рэлея структура эволюционирует к хаотическому перемешиванию, проходя через несколько достаточно четко выраженных устойчивых нестационарных состояний или режимов. Физическая природа выхода системы из стационарного состояния и последующая смена режимов остаются не понятыми до конца и определяются, по всей видимости, различными причинами для каждой конкретной модели. В настоящей работе обсуждаются возможные причины смены режимов тепловой конвекции в верхней мантии под литосферой постоянной мощности.
Эволюция конвективных течений в пространстве и времени и сама структура зависят, помимо числа Рэлея, от целого ряда параметров модели. Так, в [9] показано, что механически прочная литосфера подавляет возникновение нестабильности в верхнем тепловом граничном слое, тем самым отодвигая переход к нестационарности в область более высоких значений числа И,а. В качестве граничного условия на подошве литосферы автором [9] было использовано постоянное значение горизонтальной составляющей скорости, что моделировало плитные движения. В нашей модели мы используем условие прилипания, т.е. равенство нулю горизонтальной компоненты скорости, как частный случай плитных движений, когда плита неподвижна. Наличие горизонтальной скорости на подошве литосферы приводило не только к стабилизации структуры течения, но и к вытягиванию горизонтального размера ячейки до аспектного отношения, равного 6. Другим важным параметром, влияющим на устойчивость структуры течения, является размерность модели и аспектное отношение исследуемой области Л, которое есть отношение горизонтального размера области к вертикальному. В работе [10] численно изучались условия перехода
и взаимосвязь существенно двумерных структур течения и трехмерных. Для жидкости при И,а ~ 105, структура течения в виде валов, которая является классическим примером двумерной конвекции, оказалась нестационарной и эволюционировала со временем в стационарную трехмерную гексагональную структуру. Этот результат говорит о том, что изучение конвекции в двумерном варианте для величины порога стационарности дает оценку снизу.
Максимальные значения числа Рэлея для стационарного режима были получены при моделировании конвекции в квадратной области. Так, в [11] было обнаружено, что для определенных начальных условий конвекция остается стационарной при И,а = 108, а в [12] получили подобный режим конвекции даже для более высоких значений И,а. Этот результат не нашел подтверждения в последующих численных экспериментах, поэтому вслед за [8] можно сказать, что результаты изучения структуры и эволюции тепловой конвекции, полученные в квадратной области, имеют довольно ограниченную область применения. Подобную геометрию области можно использовать, например, для тестирования численных схем [13], но вряд ли она пригодна для моделирования процессов в мантии. Чтобы избавиться от выраженного стабилизирующего эффекта боковых границ, наблюдаемого в квадратном боксе, временную зависимость конвекции изучают теперь в областях с ас-пектным отношением от 5 до 10 [14, 15] или в кольце, моделирующем всю мантию, что вообще исключает влияние вертикальных границ [16].
Помимо перечисленных выше свойств моделей, важную роль в эволюции системы при численном моделировании играют начальные условия. В работе [8] предпринято специальное исследование влияния начального распределения температуры на эволюцию конвективного течения. Взяв в первом эксперименте при И,а = 5 • 105 в качестве начального распределения классический установившийся кондуктивный профиль температуры (линейная зависимость температуры с глубиной), авторы получили при £ ^ то вытянутую ячейку в области с Л = 3, которая периодически разбивалась на две. Во втором эксперименте вместо кондуктивного профиля в качестве начального распределения температуры ими была выбрана тепловая структура от промежуточного состояния первого эксперимента, когда в расчетной области существовали три примерно изометрические конвективные ячейки. Полученная в результате структура даже при И,а = 106 оказалась стационарной и состояла из все тех же трех ячеек с аспектным отношением 0.8. Интерпретацию полученных результатов авторы выполнили с помощью траекторий состояния на фазовой плоскости в координатах "производная кинетической энергии системы по времени — кинетическая энергия, ", хотя, по-видимому, разумнее с физической точки зрения фазовую плоскость строить в пространстве Ми— [5]. Траектория состояний для второго эксперимента представляла собой спираль, сходящуюся к точке в центре. Эта точка, называемая предельной точкой [2], отвечает абсолютно стационарному состоянию системы, а область фазового пространства в окрестности этой точки — областью притяжения, или аттрактором. Начинающиеся в пределах этой области траектории в конце концов сходятся к этой точке или к замкнутой линии, отвечающей состоянию периодического, колебательного характера. Ясно, что начальное состояние во втором эксперименте находилось в области притяжения. Это и обусловило выход системы на стационар.
Авторы заключили, что вообще неуместно задавать вопрос о том, является ли конвекция при определенном числе Рэлея зависящей от времени или стационарной. Она может быть всякой в зависимости от выбора аспектного отношения и начального условия. Подобный вывод следует из нелинейности уравнений, описывающих процесс тепловой конвекции. Именно это свойство уравнений и, вероятно, сама природа данного процесса обладают
внутренней неоднозначностью.
На наш взгляд, выбор указанных характеристик (аспектного отношения и начальных условий) диктуется ситуацией в недрах Земли и задачей, которую поставил перед собой исследователь. Аспектное отношение при этом может отличаться на порядки, поскольку механизм конвекции работает, например, в моделях изометрических магматических камер с Л ~ 1 [17] или во всей верхней мантии, где ячейки тепловой конвекции могут увеличиваться до аспектного отношения порядка 10. Начальное распределение температуры, по-видимому, также необходимо выбирать с учетом современных представлений об эволюции теплового состояния, динамики недр и возраста Земли. Учитывая, что ранняя мантия Земли была более горячей с более интенсивным конвективным перемешиванием в ограниченных по толщине приповерхностных слоях и весьма нестационарна, разумно, на наш взгляд, выбрать начальное распределение температуры, которое порождается нестационарными изометрическими ячейками. Наиболее подходящим распределением температуры представляется то, которое было предложено в [18], где на профиль температуры от остывающего полупространства накладывалось малое возмущение в виде суммы десяти синусоид разных периодов. На начальном этапе структура конвекции в вытянутой области состояла из нескольких ячеек с различным аспектным отношением. Таким образом, данное начальное условие предоставляло системе самой выбрать вид течения и его эволюцию в зависимости от других характеристик модели.
2. Математическая модель 2.1. Постановка задачи
Уравнения, описывающие движение конвектирующей жидкости и переноса тепла в ней (приближение Обербека—Буссинеска) для случая бесконечного числа Прандтля, задаются следующим образом в безразмерном виде:
дТ дТ дТ _ _д_ / дТ \ д_ / дТ \ дЬ дх ду дх \ дх) ду \ ду) '
д2 д2
ду2 дх2
V
д2ф_д2ф\\ {д^ \ _ кадТ (2)
ду2 дх2) ) дхду \ дхду / дх '
и _ — V _ - — (3)
ду' дх'
Здесь х, у — оси декартовой системы координат (х _ у _ 0 в левом верхнем углу), ось у направлена вверх; и, V — компоненты вектора скорости; Ь — время; Т — температура; а — коэффициент температуропроводности; А — скорость генерации радиоактивного тепла; ф — функция тока; V(х,у) — кинематическая вязкость; И,а — число Рэлея.
Обезразмеривание величин выполнялось следующим образом (безразмерные величины обозначены штрихом):
(х, у) _ ¿(х, у)', V _ ^ ■ _ ¿¡30Т', А _ А'д/С,
(и, V) _ (и, ь)'ао/д,, ф _ а0 ■ ф', Ь _ Ь'СС2/а0.
Здесь й — характерный размер области, в данном случае он равен мощности верхней мантии минус мощность литосферы, ао = Ам/доср — характерное значение коэффициента температуропроводности (термодиффузии), ср — удельная теплоемкость при постоянном давлении, до — плотность, Ам — коэффициент теплопроводности в мантии (в коре литосферы а равно отношению коэффициентов теплопроводности коры (Ак) и мантии), и0 = до/до — характерное значение вязкости (до — динамическая вязкость), во = —дТ/ду — вертикальный градиент температуры, Q = Амво — тепловой поток на дне верхней мантии.
Значения параметров среды в модели: й = 660 км — 120 км = 540 км, ао = 8 • 10-7 м2/с, Q = 27 • 10-3 Вт/м2, Ам = 4 Вт/м0С, Ак = 2.2 Вт/м0С, во = 6.75 • 10-3 0С/м, а = 3 • 10-51/0С, До = 1021 Па-с, д = 9.8 м/с2, до = 3.3 • 103 кг/м3; Иа = дадовЫ4 / (аоДо) = 7.103 • 105.
Граничные условия и геометрия расчетной области выбирались следующим образом: Ьх = 6000 км — горизонтальный и Ьу = 660 км — вертикальный размеры области, Ьх/Ьу = 9 : 1 (рис. 1); Т = 00С на дневной поверхности (у = 0); дТ/ду = во на дне верхней мантии (у = —660 км); дТ/дх = 0 на вертикальных границах (х = 0, 6000 км).
Для скоростей заданы условия проскальзывания: ф = 0, д2ф/ди2 = 0 на всех границах, кроме подошвы литосферы, где задано условие прилипания ф = 0, дф/ди = 0.
°1
я" 200-
к
й К р 400-
ю >> [-1 600-
т = о
6000 км
120 км Литосфера
- и = V = 0 дТ дх
- Конвектирующая ди
- мантия и = — = 0 ох
дт_
"^М "7Г~ - Ч
ду
ди ду
= V = 0
Рис. 1. Геометрия области и граничные условия.
Величина Q выбиралась таким образом, чтобы при расчетах получить значения теплового потока на поверхности, совпадающие с наблюдениями.
Начальное распределение температуры было выбрано в виде [19] (г = — у)
Тн = ТовгГ
+ (Ть — То)
1 — ег£
гь — г 4л/йоГа
где То = 13500С, Ть = 18000С, гь = 660 км, га = 14.04 млн лет (1 млн лет « 3.15 • 1013 с). Тогда дТ/дг = во при г = 660 км.
После обезразмеривания Тн к нему добавляется малое возмущение в [18]:
1 Л1
пу\ /кпх\ , , ят I — €0^1 —-— . 10 к \Ьу] \Ьх]
в = -У 10
г
2.2. Численный алгоритм
Дискретизация расчетной области по переменным x и y проводилась следующим образом:
Xi = Xi-i + i = 1,..., N1, yj = yj-i + (hy)j, j = 2,..., N2,
здесь (hx)i, (hy)j — массивы шагов сетки по горизонтали и вертикали, N1 — число узлов сетки по горизонтали, N2 — число узлов сетки по вертикали; (hx)i = const, а в y-направлении сетка могла сгущаться в области пограничных слоев.
По переменной t шаг т принимался постоянным, n — номер шага по времени:
fn = f (tn),fn+1 = f (tn + т).
Алгоритм для решения системы уравнений (1)-(3) состоял в следующем. Пусть в момент времени t = tn значения искомых функций (—, T, u, v) известны. Для того, чтобы найти значения этих функций в момент t = tn+1, необходимо выполнить цикл из трех последовательных шагов.
На первом шаге численное интегрирование уравнения теплопереноса (1) осуществляется мажорантной схемой расщепления [20, 21] с аппроксимацией конвективных слагаемых направленными разностями:
Tn+1/2 _ Tn
-т-+ ^ (unTn+1/2) = Л11Т n+1/2 + An,
Tn+1 _ Tn+1 /2
-T- + Л^ VnT "+1) = Л22Tn+1.
Здесь
^(«T)i = 0.5(u - |u|)Tx,i + 0.5(u + |u|)Tx,i, ЛllTi = [ai+1Tx,i - aiTX,i] /hi, a = 0.5(a + ai-1); fn = f (tn,Xi,yj), hi = Xi - Xi-1, hi = 0.5(hi + hi+1);
fx,i = (fi+1 — fi)/hi+1; fx,i = (fi — fi-1)/hi; fx,i = (fi+1 — fi)/hi.
Индекс j при этом фиксирован. Операторы Л2 и Л22 определяются аналогичным образом для переменной y.
Граничные условия для температуры строятся следующим образом. Поскольку внизу области (j = 1)
— = - Q (—1 = -—(3Ti1 - 4Ti2+Ti3),
dy Л м \dyj i , 1 2hy i '1 i '2 i '
получим
Ti , 1 = [2hyQ/Лм + 4Ti ,2 - Ti ,з]/3.
На втором шаге при решении уравнения четвертого порядка с переменной вязкостью для функции тока (2) используется итерационный метод стабилизирующей поправки [22]:
ф8+1/2 _ ф8
ф + Лц (^ПЛПф8+1/2) + Л22 (^гаЛ22„0 + 4Лю (^Люф'
Оп
-Лп(ипЛ22ф*) - Л22(^Лпфв) = ИаА!Тп,
^-О^- + Л22 (- Л2^иПЛ22/ф3) = 0.
Здесь Л11Ф = фхх, Л22Ф = фуу, Л12Ф = 0.5(фуж + фух), Лп(^Лпф) = (иф^ Л22 (V Л22Ф) = ^фуу)уу, ЛП(^Л22ф) = (^фуу)йх, Л22 (V Лцф) = (^фхх )уу;
п — номер шага по времени, 5 — индекс итерационного цикла, оп — итерационный параметр:
п . \ hihi+lhjhj+l\ . . оп = Ш1П ^-„ ^ , г = 1,..., N1; ^ = 1,..., N2.
ij I vn-'J К. "ij
Уравнение (2) решается только в области конвектирующей мантии; в литосфере функция тока полагается равной нулю.
Граничные условия для функции тока реализуются с применением следующих конечно-разностных аппроксимаций для случаев проскальзывания (д2ф/дп2 = 0) или прилипания (дф/дп = 0):
/дф\ 1 /дф\ 1
UJ м = - 2hy t3^1 - 4„i'2 + Лду) N = (3ф^ - 4ф^-1 + ф^-2)'
д 2ф\ 1 (д 2ф\ 1
—2 I = -2 (2ф1'1-5ф1'2+4ф1'3-ф1,^ Л —j 1 = — (2ф^-5„i,N2-1+4„i'N2-2—ф^-з) , дУ / i'1 hy \дУ J i'N2 hy
i = 2, 3,N1 — 1; hy = const.
На третьем шаге компоненты вектора скорости определяются из (3) путем численного дифференцирования функции тока. И, наконец, для сходимости алгоритма (в связи с нелинейной зависимостью между определяемыми величинами) описанные выше три шага необходимо выполнить 2-3 раза.
Численная модель тестировалась путем сопоставления с результатами расчетов, представленных в [13, 23]. Подробно результаты тестирования этой и других численных моделей описаны в [24].
Расчеты выполнялись на сетке 201 х 67 (т.е. hx = 30 км, hy = 10 км) с шагом по времени 0.1 млн лет. Вычисления проводились до значения времени 10 000 млн лет (100 000 шагов).
3. Установление стационарного режима конвекции
В данной работе исследование верхнемантийной конвекции под континентами мы ограничили случаем литосферы постоянной мощности. Конвекция изучалась в области с аспект-ным отношением Л = 9. В рамках нашей модели литосфера была представлена жестким кондуктивным слоем мощностью в 120 км, вещество которого не принимает участия в конвективном мантийном перемешивании. Такая постановка задачи сокращает перепад вязкости в области конвекции до 1-2 порядков. Слой пониженной вязкости (астеносфера) располагается непосредственно под литосферой. Наиболее распространенное объяснение
природы астеносферы состоит в том, что геотерма максимально близко подходит здесь к кривой солидуса базальта. Однако Karato [25], например, связывает понижение вязкости под литосферой с резко повышенной концентрацией водорода (до 100-1000 г/т H/Si), что не требует теплового размягчения или появления расплава для объяснения природы астеносферы. Поэтому в нашей модели мы будем учитывать наличие астеносферы просто зависимостью величины вязкости от глубины безотносительно к ее природе, исходя из современных работ по изучению послеледниковых изостатических движений поверхности [26]. Для континентальных областей параметры астеносферы, удовлетворяющие этим данным, включают величину ее вязкости 1.4 • 1019 Па-с при мощности слоя 70-75 км. Вязкость остальной верхней мантии принята как 1021 Па-с.
Еще одной особенностью континентальных областей по сравнению с океаническими является концентрация долгоживущих радиоактивных изотопов в коре континентов. Тепло, выделяющееся при их распаде, вносит существенный вклад в наблюдаемый тепловой поток, поэтому его необходимо учитывать при моделировании. А. Д. Дучковым с коллегами на основе многолетних наблюдений за тепловым потоком Сибири и лабораторных исследований теплофизических свойств образцов были созданы тепловые модели региона [27, 28]. Средний тепловой поток через поверхность Западно-Сибирской плиты составляет 53 • 10-3 Вт/м2. Генерация тепла радиоактивными изотопами в коре платформ рассчитывалась из двухслойной модели коры. Для литосферы постоянной мощности было использовано постоянное по латерали распределение изотопов, соответствующее параметрам Западно-Сибирской плиты: в первом слое мощностью в h1 = 14 км теплогенерация составляла 1 = 0.7 • 10-6 Вт/м3, в следующем слое, расположенным ниже h2 = 17 км, 2 = 0.3 • 10-6 Вт/м3.
Эволюцию конвективного течения под плитой и выход системы к стационарному состоянию будем исследовать с помощью зависимости числа Нуссельта от времени, приведенной на рис. 2, что является достаточно распространенным приемом [8].
Как уже упоминалось, начальное малое возмущение теплового поля в исследуемой области было выбрано нами в виде суммы десяти синусоид различного периода. Поэтому на начальном этапе эволюции сформировалось шесть конвективных ячеек, аспектное отношение которых варьировалось от 3 до 0.8 (рис. 3).
На графике Nu(t) данному режиму конвекции соответствует отрезок A-B, отличающийся наиболее высокими значениями числа Нуссельта, Nu ~ 11 (рис. 2).
Объяснение этому факту дается в [29], где показано, что число Нуссельта уменьшается с ростом горизонтального размера конвективных ячеек. Максимальному значению Nu соответствует аспектное отношение ячейки 0.8. В нашем случае отрезку A—B отвечает максимальное количество ячеек с минимальным аспектным отношением.
Другой особенностью данного режима является высокочастотные временные вариации числа Нуссельта, обусловленные периодическим возникновением тепловой неустойчивости в граничных тепловых слоях ячеек. Область развития такой неустойчивости в граничном слое характеризуется локальным увеличением мощности слоя, однако при данных параметрах конвекции потенциальной энергии, запасенной в районе неустойчивости, недостаточно для ее отрыва. Поэтому описываемые неустойчивости сносятся ламинарным течением граничного слоя в вертикальные потоки ячейки, где и диссипируют. Причина периодического возбуждения граничных слоев состоит в том, что конвектирующая область и, что важнее, перекрывающий ее кондуктивный слой литосферы еще не прогрелись. В литосфере существуют заметные латеральные вариации температуры, которые, взаимодействуя с вариациями теплового граничного слоя конвективной ячейки, увеличивают
1000 2000 3000 4000 5000 6000 7000 8000 Рис. 2. Зависимость (Т), Ми и от времени.
или уменьшают тепловой поток в литосферу, что и ведет к формированию локальной неустойчивости. Однако с течением времени, как видно из графика частота вариа-
ций данного параметра со временем уменьшается, это особенно заметно в последней трети отрезка А—В. Одновременно с этим уменьшается и число конвективных ячеек в слое. Так, к I ~ 4 000 млн лет ячеек осталось только три — две с аспектным отношением 3 и 2 и одна в центре с горизонтальным размером, в 5 раз превышающим вертикальный (рис. 4).
К этому времени литосфера прогрелась настолько, что уже не порождает заметных латеральных вариаций теплового потока из конвектирующей мантии. Причина возникновения локальных возбуждений верхнего теплового граничного слоя исчезла, поэтому зависимость числа Нуссельта от времени становится более гладкой, а его значение продолжает уменьшаться, поскольку аспектное отношение ячеек растет.
В интервале В—С продолжается увеличение горизонтального размера центральной ячейки. Ее аспектное отношение увеличивается при этом с 5 до 7. Рост данной ячейки осуществляется за счет уменьшения размеров ячеек, расположенных слева и справа от центральной. Эффективность конвекции в выносе тепла уменьшается, что фиксируется монотонным понижением значения числа Нуссельта.
Исчезновение правой ячейки отмечено на графике числа Ми(£) его резким уменьшением. Этому событию соответствует положение буквы С. Поведение Ми на отрезке С—Б отражает процесс исчезновения левой приграничной ячейки, что отвечает минимуму числа Ми в этом интервале (рис. 5).
После I = 7000 млн лет, чему соответствует буква Б на рис. 2, структура течения
Рис. 3. Поле скорости (см/год), функция тока и температура (о0) на момент времени £ = 2 500 млн лет.
представляет собой одну вытянутую ячейку с аспектным отношением 9 (рис. 6).
Этот режим оставался неизменным до конца счета, т.е. в течение более 3 000 млн лет. Следует отметить, что на этом отрезке времени оставались неизменными средняя по области температура (Т), число Нуссельта Ми и среднеквадратичная скорость движения вещества Угт8. Все это позволяет предположить, что здесь было достигнуто стационарное состояние течения.
Время установления стационарного режима является одним из базовых параметров, который характеризует эволюцию тепловой конвекции. В одной из первых современных работ по численному моделированию тепловой конвекции [4] была предложена простая формула, определяющая время выхода конвектирующей системы на стационарное состояние. Авторы [4] предположили, что главное условие стационарного режима — кондук-тивный прогрев всего вещества конвективной ячейки. Прежде всего это касается изотермического ядра, т.е. центральной части ячейки, вращающейся практически как единое целое при режиме развитых пограничных слоев. При таком подходе предполагалось, что начальное распределение температуры по глубине задавалось в виде линейной функции. Время достижения стационарного состояния было оценено ими как Т1 = 42/(кп2). Для верхней мантии Земли 4 = 670 км и время выхода на стационарный режим оценивается тогда величиной т1 = 1700 млн лет. Однако эта формула, как резонно было отмечено В. П. Трубицыным с коллегами [30], не учитывает того факта, что интенсивность прогревания зависит от интенсивности конвекции, которая хорошо описывается числом Нуссельта.
Рис. 4. Поле скорости (см/год), функция тока и температура (о0) на момент времени £ = 4 000 млн лет.
Поэтому они предложили свой вариант оценки т2 _ С2/(&№и), где Ми — число Нуссельта. В этом случае т2 _ 1500 млн лет, поскольку, как следует из выполненных в настоящей работе расчетов, число Нуссельта для верхнемантийной конвекции под жесткой литосферной плитой континента не превышает 10.
В дальнейшем авторы [4] уточнили выражение для т, более строго оценив время кон-дуктивного нагрева изотермического ядра, и в окончательном виде предложили формулу вида т3 _ С2/(кпМи), что в нашем случае соответствует значению т3 _ 500 млн лет. Сравнивая представленные выше методы оценок и сами значения с результатами выполненного в настоящей работе моделирования конвекции, обнаруживается явное несоответствие между результатом оценок и экспериментом.
Для нашей модели время выхода на стационарное состояние не ниже т ~ 7 000 млн лет, что в 5-10 раз выше оценочных значений. Если принять, что главной причиной установления конвекции является прогревание и выравнивание температуры в конвектирующей области, то становятся понятными завышенные значения периода установления для нашей модели. В отличие от моделей, для которых были сделаны приведенные выше оценки т, в нашем случае процесс перераспределения тепла осуществлялся медленнее из-за верхнего кондуктивного слоя — литосферы. Кроме того, наличие данного слоя привело к снижению эффективности передачи тепла верхнемантийной конвекцией более чем в два раза, понизив число Нуссельта до 9.7 по сравнению с общепринятым его значением для верхней мантии Ми=20 [5]. Учитывая этот факт, для моделей верхнемантийной конвекции под жесткой кондуктивной литосферой можно предложить формулу для оценки с более
Рис. 5. Поле скорости (см/год), функция тока и температура (о0) на момент времени £ = 6 000 млн лет.
слабой зависимостью от числа Нуссельта, например в виде т = й2 / (ykNu1/^. Это даст значение т — 7 500 млн лет, сопоставимое с результатом эксперимента.
В процессе численных экспериментов изучались также зависимости от времени и других обобщенных характеристик течения. Изменение во времени средней квадратичной скорости течения вещества Кт8(£) повторяет в общих чертах поведение числа с
единственным исключением. В отличие от Ми(£), среднее значение Кт8(£), около которого варьировались мгновенные значения скорости, оставалось неизменным — порядка 1.1 см/год (см. рис. 2). Локальные временные вариации скорости обусловлены формированием тепловых неустойчивостей в граничных слоях, а постоянство усредненного значения скорости Угт8 можно связать с постоянством граничных условий и отсутствием диссипа-тивных членов в уравнениях, описывающих движение. Временные вариации средней по объему температуры (Т)(£) имеют локальные экстремумы, связанные, вероятно, с режимами перестройки структуры течения. На заключительном этапе эксперимента средняя температура, как и другие характеристики, приблизилась к своему постоянному значению (Т) - 1150 оС.
Аспектное отношение стационарной структуры течения в эксперименте оказалось равным 9, что соответствует горизонтальному размеру ячейки около 6000 км. Подобный размер был получен в экспериментах [31]. Это говорит о том, что введение в рассмотрение верхнего кондуктивного слоя подобно заданию на верхней границе конвектирую-щей области граничного условия в виде постоянного теплового потока. Горизонтальному
Рис. 6. Поле скорости (см/год), функция тока и температура (о0) на момент времени £ = 8 000 млн лет.
вытягиванию ячейки способствовал также тот факт, что на нижней границе также задавался постоянный тепловой поток, а не температура. Предпочтение такому граничному условию (как и в [32]) было отдано потому, что связь между верхней и нижней мантией на глубине 670 км имеет выраженную тепловую, а не механическую природу. Другими словами, из-за различия вязкости в верхней и нижней мантии более чем на порядок конвективные ячейки в мантиях имеют слабую механическую связь, т.е. вращаются без существенного взаимодействия. При таком режиме значение теплового потока в верхнюю мантию обнаруживают слабые латеральные вариации по сравнению с самой температурой, т.е. поток почти не зависит от x-координаты и в первом приближении его можно считать постоянным (Ellsworth, Schubert [33]).
В заключение подчеркнем, что главный результат данного исследования, как нам кажется, состоит в том, что конвекция при Ra = 7 • 105 в расчетной области с аспектным отношением 9, через t ~ 7000 млн лет вышла на стационарный режим перемешивания. Этот вывод, помимо визуального наблюдения за структурой течения, подтверждается поведением эволюционных характеристик (T)(t), Nu(t), Vrms(t).
С одной стороны, этот результат подтверждает вывод [9] о том, что наличие верхнего прочного кондуктивного слоя — литосферы — подавляет возникновение тепловых возмущений в тепловом слое; последнее отодвигает границу возможных стационарных течений в область более высоких значений чисел Рэлея. В данном случае эта область увеличилась почти на порядок. С другой — время счета было ограничено и составляло около
10 000 млн лет. Это не очень много для исчерпывающего анализа эволюции, поскольку, как было показано в [8, 14], полученное нами стационарное состояние может сохраняться 3 000 - 5 000 млн лет, а затем система вновь эволюционирует в новое устойчивое состояние, сохраняющееся опять длительное время.
Объяснение подобному поведению, как уже упоминалось, следует искать в нелинейности самой задачи тепловой конвекции. Времена порядка t ~ n • 1000 млн лет, где n > 5, представляют в основном интерес лишь для гидродинамики явления, а не для геологии, поскольку возраст Земли не превышает 5 000 млн лет. Поэтому, по-видимому, вряд ли можно считать, что структура тепловой конвекции под обширными континентальными областями стационарна (если не принимать во внимание латеральных вариаций мощности литосферы).
Авторы благодарят Г.Г. Черных за полезное обсуждение численной модели.
Список литературы
[1] SoLHEIM L.P., PELTIER W. R. Mantle phase transitions and layered convection. Can. J. Earth Sci., 30, 1993, 881-892.
[2] Ландау Л. Д., Лифшиц е. М. Гидродинамика. Наука, М., 1986.
[3] Krishnamurti R. On the transition to turbulent convection. J. Fluid Mech., 42, 1970, 295-320.
[4] McKenzie D. P., Roberts J. M., Weiss N. O. Convection in Earth's mantle: towards a numerical simulation. J. Flrnd Mech., 62, Pt 3, 1974, 465-538.
[5] Трубицын В. П., НиколАйчик В. В. Режимы тепловой конвекции. Физика Земли, 6, 1991, 3-12.
[6] Добрецов Н. Л., Кирдяшкин А. Г. Глубинная геодинамика. ОИГГМ СО РАН, Новосибирск, 1994.
[7] Travis B., Olson P. Convection with internal heat sources and thermal turbulence in the Earth's mantle. Geophys. J. Inter., 118, 1994, 1-19.
[8] Hansen U., Ebel A. Time-dependent thermal convection — a possible explanation for a multiscale flow in the Earth's mantle. Geophys. J., 94, No. 2, 1988, 181-191.
[9] Davies G. F. Role of the lithosphere in mantle convection. J. Geophys. Res., 93, 1988, 10451-10466.
[10] Travis B., Olson P., Shubert G. The transition from two-dimentional to three-dimentional planform in infinite Prandtl number thermal convection. J. Flrnd Mech., 216, 1990, 71-91.
[11] JARVIS G. T. Time-dependent convection in the Earth's mantle. Phys. Earth Planet. Interiors, 36, 1984, 305-327.
[12] Schubert G., Anderson C.A. Finite element calculations of very high Rayleigh number thermal convection. Geophys. J. Roy. Astr. Soc., 80, 1985, 289-318.
[13] BLANKENBACH B. ЕТ AL. A benchmark comparison for mantle convection codes. Geophys. J. Int., 98, 1989, 23-38.
[14] CHRISTENSEN U. Time-dependent convection in elongated Rayleigh—Benard cell. Geophys. Res. Lett., 14, 1987, 220-223.
[15] Трубицын В. П., Бобров А. М., Кувышкин В. В. Влияние континентальной литосферы на структуру мантийной тепловой конвекции. Физика Земли, №5, 1993, 3-11.
[16] GURNIS M., ZONG S. Generation of long wavelength heterogeneity in the mantle by the dynamic interaction between plates and convection. Geophys. Res. Lett., 18, 1991, 581-584.
[17] ШАРАПОВ В. Н., Миловл Л. В. Динамика гранитизации магмы стационарным потоком флюида при развитии конвективного плавления пород земной коры. В "Динамические модели физической геохимии". Наука, Новосибирск, 1982, 16-19.
[18] CHRISTENSEN U. Convection with pressure- and temperature-dependent non-Newtonian rheology. Geophys. J. Roy. Astr. Soc., 77, 1984, 343-384.
[19] FLEITOUT L., YUEN D. A. Secondary convection and the growth of the oceanic lithosphere. Phys. of the Earth and Planetary Interiors, 36, 1984, 181-212.
[20] ЯНЕНКО Н. Н. Метод дробных шагов решения многомерных задач математической физики. Наука, Новосибирск, 1967.
[21] МАРЧУК Г. И. Методы расщепления. Наука, М., 1988.
[22] HOUSTON M.H., Jr., De BREMAECKER J.Cl. ADI solution of free convection in a variable viscosity fluid. J. of Comput. Physics, 16, 1974, 221-239.
[23] MOORE D.R., WEISS N.O. Two-dimensional Rayleigh-Benard convection. J. Fluid Mech, 58, Part 2, 1973, 289-312.
[24] Chernykh G.G., MOSHKIN N.P., RYCHKOVA E.V., TYCHKOV S.A. Comparison of some numerical algorithms for two-dimensional convection of fluid with nonlinear viscosity. In "Int. Conf. on the Methods of Aerophys. Research: Proc., Pt I", Novosibirsk, 1996, 79-84.
[25] KARATO S. The role of hydrogen in the electrical conductivity of the upper mantle. Nature, 347, 1990, 272-273.
[26] FjELDSKAAR W. Viscosity and thickness of the astenosphere detected from Fennoscandian uplift. Earth Planet. Sci. Lett., 126, 1994, 339-410.
[27] ДУЧКОВ А. Д., СОКОЛОВА Л. С. Геотермические исследования в Сибири. Наука, Новосибирск, 1974.
[28] ДУЧКОВ А. Д., БАЛОБАЕВ В. Т., ВОЛОДЬКО Б. В. и др. Температура, криолитзона и радиогенная генерация в земной коре Северной Азии. ОИГГМ СО РАН, Новосибирск, 1994.
[29] Olson P., CoRcos G. M. A boundary layer model for mantle convection with surface plates. Geophys. J. Roy. Astr. Soc., No. 62, 1980, 195-219.
[30] Трубицын В. П., Васильев П. П., Карасев А. А. Конвекция при неравномерно распределенных источниках тепла. Физика Земли, №7, 1984, 13-21.
[31] Hewitt J.M., McKenzie D.P., Weiss N.O. Large aspect ratio cells in two-dimensional thermal convection. Earth Planet Sci. Lett., 51, 1980, 370-380.
[32] Cherepes L., Rabinowicz M. Gravity and convection in a two-layer mantle. Ibid., 76, 1985, 193-207.
[33] Ellsworth K., Schubert G. Numerical models thermally and mechanically coupled two-layer convection of highly viscous fluids. Geophys. J., 93, 1988, 347-363.
Поступила в редакцию 1 октября 1997 г.
