ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ, 2021, том 47, № 7, с. 525-534
ПАРАМЕТРИЧЕСКАЯ ТЕПЛОВАЯ МОДЕЛЬ ЭВОЛЮЦИИ ЗЕМЛИ
© 2021 г. М. Ю. Решетняк1,2*
1Институт физики Земли РАН, Москва, Россия
2Институт земного магнетизма и распространения радиоволн РАН, Москва, Россия
Поступила в редакцию 02.06.2021 г.
После доработки 18.06.2021 г.; принята к публикации 29.06.2021 г.
Рассмотрена модель совместного остывания ядра и мантии Земли на интервале времени с момента
аккреции, 4.5 млрд лет назад, и до 1.5 млрд лет в будущее. Модель позволяет получить современный
радиус твердого ядра, близкий к наблюдаемому тепловой поток на поверхности планеты, удовлетво-
рить ограничениямпо температуре и вязкости мантии, параметрам слоя D′′ и литосферы. Исследованы
сценарии эволюции Земли, при которых тепловой поток на границе ядро-мантия не претерпевал
значительных изменений. Приведены следствия для геодинамо.
Ключевые слова: эволюция Земли, ядро, мантия, обратная задача, геодинамо.
DOI: 10.31857/S032001082107007X
ВВЕДЕНИЕ
Qb обосновано. Важность данного вопроса обу-
словлена ограничениями, вносимыми геодинамо.
После завершения процесса аккреции и появ-
Согласно палеомагнитным наблюдениям, нет ука-
ления жидкого ядра (4.5 млрд лет), длившиеся
заний на существенное увеличение напряженности
∼50-100 млн лет (Дрейк, 2000; Соломатов, 2007),
геомагнитного поля (см. обзор Решетняка, Павло-
ядро Земли, имевшее температуру порядка T◦ =
ва, 2016), напрямую связанное с увеличением Qb.
= 6000◦K (Лаброзе и др., 1997), начало осты-
Для ответа на вопрос о постоянстве Qb необходимо
вать, отдавая свое тепло мантии. Насколько про-
решить совместные уравнения теплообмена в ядре
цесс остывания был монотонным, сказать сложно,
и мантии. Этот подход интересен еще и тем, что
поскольку конечный ответ зависит от начальных
полученное решение может быть проверено на
условий и точной оценки концентрации радиоак-
соответствие некоторым известным величинам
тивных элементов. Тем не менее на характерных
для мантии, например, температуры, вязкости
временах порядка миллиарда лет температура ядра
(Трубицын, Трубицын,
2020; Трубицын,
2016),
уменьшалась, со временем появилось твердое ядро
потока Qs, чисел Рэлея, которые также могут быть
и, возможно, небольшой сферический слой вблизи
использованы для проверки работоспособности
границы ядро-мантия с субадиабатическим гради-
моделей остывания Земли. Точность определения
ентом температуры.
этих величин хоть и не высока, но в ряде случаев
За последние два десятилетия выработался
может быть выше точности величин, связанных
непосредственно с ядром.
некоторый стандарт в моделировании тепловой
эволюции ядра (Лаброзе и др., 1997; Лаброзе,
Первые параметрические модели остывания
2003). Основным параметром, регулирующим
мантии начали появляться еще в 80-е гг. (см. обзор
остывание ядра в моделях, является заданный
в Шуберт и др., 2001), но в то время модели ядра
тепловой поток на границе ядро-мантия Qb,
были развиты значительно хуже, чем сейчас, в
как правило убывающий с момента аккреции к
частности, не учитывались ограничения, связанные
настоящему времени на величину порядка 20%. В
с генерацией геомагнитного поля (Лаброзе, 2003),
то же время мощность радиоактивных элементов
только-только стали появляться скейлинговые
в мантии за 4.5 млрд лет уменьшилась приблизи-
законы, связывающие параметры модели, осно-
тельно в 4.5 раза, а величина теплового потока на
ванные на результатах трехмерного моделирования
поверхности планеты Qs — на порядок (Шуберт
тепловой конвекции и динамо. Неудивительно,
и др., 2001). В этой связи возникает вопрос,
что с тех пор некоторые параметры, в частности
насколько предположение о медленном изменении
входящие в зависимости вязкости от температуры,
претерпели изменение. В целом же модели остыва-
*Электронный адрес: m.reshetnyak@gmail.com
ния ядра не потеряли своей актуальности, и могут
525
526
РЕШЕТНЯК
быть использованы и в настоящее время, хотя и
(рис. 1). Радиус твердого ядра c увеличивался со
требуют некоторого переосмысления. Их ценность
временем. Внутри твердого ядра 0 ≤ r ≤ c (I) рас-
связана еще и с тем, что современные трехмер-
пределение температуры определяется кондуктив-
ные модели тепловой конвекции не могут быть
ным теплообменом и радиоактивными источника-
использованы для описания тепловой конвекции на
ми. В жидком ядре r ≤ rb (II) температура распре-
ранних этапах эволюции Земли, когда числа Рэлея
делена по адиабате. Модель описывает эволюцию
были еще слишком велики. В этом случае исполь-
температуры, тепловых потоков в ядре и увеличе-
зование параметрических моделей, основанных на
ние радиуса ядра c во времени t при заданном Qb.
скейлинговых зависимостях между безразмерными
Соответствующим подбором параметров удается
числами, оказывается очень удобным средством.
получить реалистичное значение современного ра-
Ниже мы рассмотрим одну из простейших пара-
диуса твердого ядра, c = 1228 км, и требуемую
метрических однослойных моделей остывания ман-
мощность для генерации магнитного поля (Решет-
тии и ядра (Стивенсон и др., 1983), которая пред-
няк, 2021). Напомним основные уравнения этой
сказывает значения Qb сразу после аккреции на два
модели.
порядка большие, чем сейчас, и покажем, как опре-
Распределения плотности ρ(r), давления P (r) и
деленным выбором параметров уменьшить Qb в
ускорение свободного падения g(r) удовлетворяют
древности, удовлетворив требованиям палеомагни-
гидростатическому балансу, задаваемому соотно-
тологов, наблюдающих геомагнитное поле, сравни-
шениями
мое с современным по напряженности, 3.5 млрд лет
∇P = -ρg,
(1)
назад (Мак-Элхини, 1980; Тардуно и др., 2010).
Модель, которая покрывает интервал 6 млрд лет
(с окончания процесса дифференциации ядра и
r
∫
мантии 4.5 млрд лет в прошлом и на 1.5 млрд лет
4πG
g(r) =
ρ(u)u2 du.
(2)
в будущее), позволяет проследить появление твер-
r2
дого ядра, эволюцию тепловых потоков на грани-
0
цах жидкого ядра, поверхности планеты, сделать
Здесь G — гравитационная постоянная. Замыкает
некоторые выводы об изменении магнитного поля
систему уравнений для трех переменных (P, ρ, g)
на геологических временах. Поскольку интервал
логарифмическое уравнение состояния (Пурье, Та-
моделирования выходит за пределы геомагнитных
рантола, 1998)
наблюдений (3.5 млрд лет), появляется возмож-
ность предсказания поведения магнитного поля на
ρ
ρ
P =K◦
ln
,
(3)
ранней стадии эволюции планеты. Данные подходы
ρ◦
ρ◦
могут быть использованы для описания тепловой
эволюции планет.
где K◦ — модуль объемной упругости, ρ◦ — плот-
ность при нулевом давлении.
Для учета скачка на границе ядро-мантия вво-
ОСТЫВАНИЕ ПЛАНЕТЫ
дится поправка
{
Рассмотрим остывание Земли радиуса rs после
ρ(r) + δρ, если r ≤ c
окончания формирования жидкого ядра радиуса rb
ρ(r) =
(4)
(рис. 1). К этому моменту в планете произошла
ρ(r), если r > c.
дифференциация вещества: ядро стало в основном
железным, а мантия — силикатной. Граница ядро-
Зная распределения (P, ρ, g), можно найти
мантия является переходом от одного химического
адиабатическое распределение температуры
состава к другому. Предполагается, что радиусы rb
∫
- α(u)g(u)
и rs уже через 50-100 млн лет после образования
Cp
Tad(r) = Tc(c) e c
du.
(5)
Земли оставались неизменными во времени. Осты-
вание планеты происходит за счет теплового потока
Здесь Tc(c) — температура на границе r = c, ко-
Qs на rs. Задача состоит в описании совместного
эффициент объемного расширения α задан соотно-
процесса охлаждения Земли: ее мантии и ядра.
шением
γCpρ◦
α(r) =
(
),
(6)
Охлаждение ядра
K◦ 1 + lnρ
ρ◦
Согласно модели, первоначально ядро радиуса
где Cp — удельная теплоемкость при постоянном
rb было полностью расплавлено (Лаброзе и др.,
давлении, γ — коэффициент Грюнайзена.
1997; Соломатов, 2007). Далее, по мере остывания
за счет теплового потока Qb на границе ядра, в
Если твердое ядро еще не появилось, c = 0, то
центре Земли появился зародыш твердого ядра
Tc(c) = Tc(0) = T◦, где T◦ — температура в центре
ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 47
№7
2021
ПАРАМЕТРИЧЕСКАЯ ТЕПЛОВАЯ МОДЕЛЬ ЭВОЛЮЦИИ ЗЕМЛИ
527
m
2
m
1
c
rb
rs
Рис. 1. Схематичное строение Земли радиуса rs, включающее твердое ядро радиуса c и границу ядро-мантия rb. На
m
границах мантии возникают тепловые пограничные слои δm1, δ
2
Земли. Значение T◦ находится из уравнения тепло-
С левой стороны в (10) стоит разность тепловых
вого баланса
потоков Qb, Qc за единицу времени, входящих в
∫rb
адиабатическую область II, справа — скорость из-
∂Tad
∂T◦S
менения энергий. Во многих работах принимается,
-4π ρCp
r2 dr = -
=Qb
(7)
∂t
∂t
что Qb = Q◦b - γt. Коэффициент γ выбран так, что-
0
бы обеспечить уменьшение потока за 4.5 млрд лет
и
на 20%.
∫r
Источник, связанный с латентной теплотой
1
rb
-
α(a)g(a) da
в (11), равен
∫
Cp
4π
0
S =
e
ρ r2 dr.
(8)
PL(c) = 4πρ(c)δSTs(c),
(12)
Cp
0
где δS — удельная энтропия кристаллизации.
Рост твердого ядра начинается, когда темпера-
Оценка изменения гравитационной энергии при
тура станет меньше температуры кристаллизации
росте твердого ядра (Лопер, 1984) имеет вид
)2(γ-1/3)
[
]
( ρ(r)
2π
c3
(c)2
Ts(r) = T◦
,
(9)
QG(c) =
GM◦δρ
1-
,
(13)
s ρ(c)
5
rb
rb
где T◦s — температура кристаллизации в центре яд-
∫rb
где M◦ =43 π
ρ r2 dr — масса ядра, постоянная
ра Земли. Кристаллизация начинается в центре яд-
0
в модели.
ра, т.е. появление твердого ядра сводится к выпол-
нению условия Tc = T◦ = Ts(r) в центре r = c = 0.
Дифференцируя (13) по времени, получаем
[
]
Далее, при росте твердого ядра, c > 0, температура
2π GM◦δρc2
(c)2
кристаллизации дает значение для адиабаты на
PG(c) =
3-5
(14)
границе твердого ядра: Tad(c) = Tc(c) = Ts(c).
5
rb
rb
Положение границы твердого ядра c может быть
Изменение теплоты, связанное с адиабатическим
найдено из баланса энергии
охлаждением в (10), дает следующий вклад:
Qb - Qc = QL + QG + QC,
(10)
rb
∫
∂Tad
где
PC = -Cp ρ
r2 dr.
(15)
∂t
QL = ċPL, QG = ċPG, QC = ċPC.
(11)
c
ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 47
№7
2021
528
РЕШЕТНЯК
Таблица 1. Параметры модели ядра
Название
Обозначение
Значение
Гравитационная постоянная
G
6.6873 × 10-11 м3 (кг с2)-1
Коэффициент температуропроводности
κ
7 × 10-6 м2 с-1
Параметр Грюнайзена
γ
1.5
Радиус жидкого ядра
rb
3.48 × 106 м
Удельная энтропия кристаллизации
δS
118 Дж (кг K)-1
Удельная теплоемкость
Cp
860 Дж (кг K)-1
Плотность при нулевом давлении
ρ◦
7500 кг м-3
Скачок плотности на границе твердого ядра
δρ
500 кг м-3
Плотность ядра
ρ
1.2 × 104 кг м-3
Модуль объемной упругости
K◦
4.76 × 1011 Па
Современный радиус твердого ядра
ĉ
1.22 × 106 м
Уравнения (10)-(15) могут быть решены отно-
первого. Другим важным отличием мантии от яд-
сительно производной по времени ċ и проинтегри-
ра является наличие ярко выраженных тепловых
рованы по времени1.
пограничных слоев, соответствующих слою D′′ и
литосфере, связанных с большими числами Рэлея
Считая температуру непрерывной на границе c,
и сильной зависимостью вязкости от температуры в
получаем, что значение температуры Ts(c) является
мантии (рис. 1). Обычно тепловыми пограничными
граничным условием для задачи теплопроводности
слоями в жидком ядре пренебрегают (Лаброзе и
внутри твердого ядра 0 < r < c (I) с движущейся
др., 1997).
границей c(t):
Для того чтобы описать скачки температуры в
∂T
= κΔT,
(16)
пограничных слоях, используют следующий эле-
∂t
гантный подход (Стивенсон и др., 1983; Шуберт
где κ — коэффициент температуропроводности.
и др., 2001), записав уравнение остывания для
В центре r = 0 имеем второе граничное условие
средней температуры Tm мантии массой Mm с
∂T
теплоемкостью Cmp:
= 0. Совместное решение уравнений (1)-(16)
∂r
∂Tm
дает распределение физических полей в областях I,
CmMm
= Qb - Qs + H(t),
(17)
p
∂t
II (см. подробнее значения параметров в табл. 1).
где Qs — тепловой поток на поверхности пла-
4
(
)
t
Охлаждение мантии
неты, H =
πCmpρm
r3s - r3b
Cmee-
e
— вклад
3
Остывание мантии происходит за счет теплово-
радиоактивных источников в мантии, убывающих
го потока в околоземное пространство. В моделях
во времени, Cme — удельная мощность радиоизо-
обычно принимается, что температура на поверх-
топов, τme — характерное время распада изото-
ности Земли Ts равна фиксированному значению.
пов. Температура Tm осреднена по конвективной
Насколько это упрощение, конечно же, неприме-
области мантии, находящейся, согласно модели,
нимое на ранней стадии эволюции планеты (Абе,
между двумя пограничными слоями δ1, δ2: Tm =
1993), влияет на изменение теплового режима,
3
∫ rs-δ2 T(r) r2dr.
=
r3s-r3
rb+δ1
обсудим далее, а пока отметим, что граничные
b
условия для задачи ядра и мантии различны. Для
В правой части (17) два неизвестных тепло-
ядра — это условие второго рода, а для мантии —
вых потока, которые требуется выразить через Tm.
Скачки температуры удовлетворяют соотношени-
1Для вычислений удобно ввести переменную f = c3/3,
ям δT1 = Tm - Tb, δT2 = Ts - Tm, δT1 + δT2 = δT,
тогда
˙f
= ċc2.
δT = Ts - Tb.
ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 47
№7
2021
ПАРАМЕТРИЧЕСКАЯ ТЕПЛОВАЯ МОДЕЛЬ ЭВОЛЮЦИИ ЗЕМЛИ
529
Поток Qs и толщина пограничного слоя δ2
о структуре слоя D′′, для которого δT1 ∼ 102 K и
выражаются через скачок температуры в мантии
δ1 ∼ 102 км.
δT
δT : Qs = Sskm
Nu, δ2 = Nu-1 (rs - rb),
Как показывают численные эксперименты, для
rs - rb
нахождения δ1 удобнее использовать эмпириче-
(
)
β
Ra
ский прием (Стивенсон и др., 1983; Шуберт и др.,
где Nu =
— число Нуссельта, k =
Racr
2001), основанный на физических экспериментах
= κρmCmp
— коэффициент теплопроводности,
(Букер, Стенгел, 1978; Рихтер, 1978) по нахожде-
нию критического числа Рэлея для пограничного
β — константа, Ss = 4πr2s — площадь поверх-
слоя с переменной вязкостью:
ности Земли. Число Рэлея задано в виде Ra =
3
gαδT (rs - rb)
gαδT1δ31
=
, где α — коэффициент объ-
Racrb =
≈ 2 × 103.
(21)
κν
κν1
емного расширения, g — ускорение свободного
Тогда
падения,
)1/3
A◦
( Racrbκν1
ν =ν◦e
Tm
(18)
δ1 =
,
(22)
gαδT1
– коэффициент гидродинамической вязкости, ν◦,
δT1
A◦ — константы. Предполагается, что Ra ≫ Racrit,
δT1 = δT - δT2, Qb = Sbkm
,
δ1
где Racrit — критическое значение, соответству-
ющее началу тепловой конвекции. Если Qb = 0,
где Sb = 4πr2b — площадь жидкого ядра. После
то уравнение (17) может быть проинтегрировано
этого уравнения (17)-(22) могут быть решены от-
по времени при заданных начальных условиях.
носительно Tm при заданном δT .
Во многих случаях так и поступали, предполагая,
Рассмотрим совместное решение задач для ядра
что тепловым влиянием ядра на мантию можно
и мантии. При постоянном Ts скачок δT зависит
пренебречь (Шуберт и др., 2001). Для описания
только от температуры Tb, которая изменяется при
же эволюции ядра знание Qb крайне важно, и его
остывании ядра. Условие сопряжения остывания
необходимо найти.
ядра и мантии следует из непрерывности темпера-
Перейдем к вычислению Qb. Возможны различ-
туры на границе rb. Пусть на начальной итерации
ные подходы, например, использование большо-
известен поток на границе ядро-мантия Qb. Этого
го числа эмпирических закономерностей (Сотин,
(в совокупности с начальными условиями) доста-
Лаброзе, 1999) как для верхней границы мантии,
точно, чтобы решить уравнения (1)-(16) для яд-
так и для нижней. Универсальность данного подхо-
ра и найти значение адиабатической температуры
да может скрадываться недостаточной точностью
Tad(rb) на внешней границе ядра. В нашем постро-
самих параметров, полученных для чисел Рэлея,
ении, основанном на пренебрежимо малом скачке
много меньших значений сразу поcле аккреции.
температуры со стороны ядра на rb, это и будет тем-
Другая возможность основана на гипотезе
пература Tb, определяющая перепад температуры
(Джарвис, 1993; Вангелов, Джарвис, 1994), со-
в мантии δT . Зная δT , можно решить (17)-(22) и
гласно которой, для нижнего и верхнего погра-
найти Tm и Qb. Итеративный процесс, включающий
ничных слоев должно выполняться равенство
расчет физических полей для ядра и мантии на
локальных чисел Рэлея:
текущем шаге по времени, заканчивается после
gαδT1δ31
gαδT2δ32
выполнения заданного критерия точности. Детали
Ra1 =
,
Ra2 =
(19)
κν1
κν2
численной реализации задачи по охлаждению ядра
и нахождению оптимальных параметров модели на
Поскольку изменение вязкости ν по глубине много
многопроцессорных компьютерах, можно найти в
больше изменений α, g, κ, последние можно счи-
работах (Решетняк, 2019, 2020).
тать постоянными. После чего получаем выраже-
ние
)1/3
(ν1 δT2
РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ
δ1 =
δ2.
(20)
ν2 δT1
Как следует из результатов работы (Стивенсон
Значения ν1, ν2 берутся при температурах (Tb +
и др., 1983), на начальном этапе охлаждения (пер-
+ Tm)/2, (Ts + Tm)/2 соответственно. Но, как
вые сотни миллионов лет) потоки на поверхности
следует из численных экспериментов, получаемые
планеты и границе ядро-мантия превосходили со-
скачок температуры δT1 ≪ δT2 и толщина погра-
временные значения на несколько порядков. Если
ничного слоя δ1 < 1 км оказываются крайне малы
для эволюции мантии большой поток Qs является
и не соответствуют современным представлениям
следствием уменьшения вязкости при больших Tm
ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 47
№7
2021
530
РЕШЕТНЯК
и может быть быстро компенсирован быстрым по-
современного радиуса твердого ядра
ĉ. Мини-
следующим увеличением вязкости в ходе охлажде-
мум функции Ψb соответствует, в дополнение к
ния мантии, то большие потоки тепла на границе
предыдущему условию, еще двум требованиям:
ядро-мантия не имеют подобного компенсацион-
приближению по величине современного теплового
ного механизма, поскольку свойства металла слабо
потока на поверхности ЗемлиQs и минимальному
зависят от температуры. В случае появления боль-
среднеквадратичному отклонению σb потока Qb,
ших потоков Qb возникает сложность с твердым
нормированному на его среднее значение Mb на
ядром, которое при возрастании Qb может стать
интервале времени, начиная с td = 200 млн лет
слишком большим (жидкое ядро может и полно-
и до конца интервала моделирования. Величины
стью затвердеть). Как уже отмечалось, палеомаг-
R1, R2
вычисляются для настоящего момента
нитные наблюдения также не отмечают сильного
времени. Значения ĉ,Qs известны из наблюдений
возрастания напряженности магнитного поля B в
(см. табл. 1, 2).
древности, связанного с тепловым потоком B ∼
Начнем с простейшего условия, потребовав
∼Q1/3b. Анализ показывает, что при небольшом от-
воспроизвести в модели современный радиус
ступлении от значений параметров, приведенных в
твердого ядра ĉ, используя для этого штрафную
(Стивенсон и др., 1983), легко попасть на режимы,
функцию Ψa и следующий набор варьируемых
когда Qb действительно становится большим, и яд-
параметров: T◦, T◦s, Cme, ν◦, A◦. Первые два
ро начинает быстро охлаждаться. Другими слова-
параметра отвечают за размер и возраст твердого
ми, вероятность получить полностью твердое ядро
ядра, Cme — за энергобюджет мантии, два по-
весьма высока. В этой связи кажется интересным
следних — за конвективный теплообмен в мантии.
исследовать сценарии эволюции, когда потоки Qb
Результаты моделирования представлены в табл. 3,
не претерпевают больших изменений. Далее мы
режим 1a. Можно заключить, что с поставленной
подберем такие значения A◦, ν◦, которые, с одной
задачей модель легко справилась: ĉ - c = 300 м.
стороны, дадут правдоподобные оценки вязкости
Возраст твердого ядра в данном расчете a = 4.5 -
в мантии, а с другой стороны, уменьшат тепловой
- tb = 2.36 млрд лет. Каких-либо прямых наблю-
поток Qb на начальном этапе развития планеты.
дений, связанных с оценкой возраста твердого
Заранее отметим, что используемые значения A◦,
ядра, нет. При определенных условиях в момент
ν◦ лежат в диапазоне значений, используемых в
образования твердого ядра может появляться
современных моделях (Шуберт и др., 2001).
избыточная энергия, доступная для геодинамо.
Прежде чем перейти к обсуждению результа-
В свою очередь палеомагнитологи не наблюдают
тов моделирования, отметим, что модель содержит
явных свидетельств повышения напряженности
большое число параметров, известных с разной
магнитного поля, предсказываемого теорией.
степенью точности (см. табл. 1, 2). Зависимость
Оценим, насколько a и c чувствительны к воз-
решения уравнений (1)-(22) от параметров может
мущению температуры поверхности Земли Ts. Со-
быть изучена как отдельно для каждого, так и путем
гласно работам (Абе, 1993; Дрейк, 2000; Солома-
нахождения оптимального решения, лучшим об-
тов, 2007), Ts могла быть весьма высока, достигая
разом удовлетворяющего заданным ограничениям
1500 градусов Кельвина и выше.
(Решетняк, 2020). В общем случае под ограниче-
Пусть Ts = T◦s + (Tmaxs - T◦s) exp(-t/Tws), так
нием понимается требование близости той или иной
что Ts ∼ Tmaxs при t ≪ Tws и Ts ∼ T◦s при t > Tws.
величины, входящей в модель (или производной от
Рассмотрим случай, в котором дополнительно к
нее), к заданному значению. В силу нелинейности
параметрам, варьируемым в режиме 1а, добавлены
задачи и большого числа параметров, для поиска
Tmaxs и Tws. Ограничение оставлено то же, что и в
оптимального решения удобно использовать метод
1a. Расчеты показали, что начальные условия су-
Монте-Карло, идеально подходящий для парал-
щественно не влияют на размер и возраст твердого
лельных компьютеров и позволяющий подобрать
ядра (см. табл. 3, режим 1b). Можно заключить, что
параметры модели в заданном диапазоне возмож-
характерное время тепловой памяти системы много
ных значений. Метод также позволят сократить
меньше возраста Земли. Это напрямую связано с
число параметров до минимума, выделив наиболее
экспоненциальной зависимостью вязкости от тем-
важные.
пературы в (18): небольшое увеличение темпера-
Ниже будет рассмотрено два варианта штраф-
туры приводит к быстрому уменьшению вязкости,
ных функций: Ψa = 1 - eR1 и Ψb = 1 -
усилению конвективного теплообмена и ускорен-
|Qs -Qs|
ному остыванию, компенсирующему первоначаль-
-e3 (R1+R2+R3), где R1 =|c-ĉ|
, R2 =
ное увеличение температуры.
ĉ
Qs
σb
В обоих расчетах 1a и 1b наблюдается боль-
и R3 =
. Минимум Ψa соответствует набору
шое отклонение теплового потока на поверхности
Mb
параметров, дающих лучшее приближение для
Земли Qs от современного
Qs. Также в начале
ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 47
№7
2021
ПАРАМЕТРИЧЕСКАЯ ТЕПЛОВАЯ МОДЕЛЬ ЭВОЛЮЦИИ ЗЕМЛИ
531
Таблица 2. Параметры модели мантии
Название
Обозначение
Значение
Радиус Земли
rs
6.371 × 106 м
Плотность мантии
ρm
3.4 × 103 кг м-3
Критическое значение числа Рэлея
Racrit
1200
Удельная теплоемкость
Cmp
1230 Дж (кг K)-1
Ускорение свободного падения
gm
10 м с-2
Коэффициент объемного расширения
αm
3 × 10-5 K-1
Коэффициент температуропроводности
κm
10-6 м2с-1
Численная константа
β
0.3
Значение современного теплового потока
Qs
44 ТВт
Таблица 3. Модельные значения cm, Qs, Qb для настоящего времени, время появления твердого ядра tb и величина
σb/Mb
№
Заданные параметры
cm, км
Qs, ТВт
Qb, ТВт
tb, млрд лет
σ/M
1a
T◦ = 6167 K, T◦s = 5366 K
1219.7
82.8
8.9
2.14
0.54
Ts = 273 K, Cme = 3.79 × 10-14 K с-1
ν◦ = 1.46 × 107 м2/с, A◦ = 7.04 × 104 K
1b
T◦ = 6167 K, T◦s = 5366 K
1210.5
82.7
8.9
2.9
0.53
Ts = 273 K, Cme = 3.79 × 10-14 K с-1
ν◦ = 1.46 × 107 м2/с, A◦ = 7.04 × 104 K
Tmaxs = 1500 K, Tws = 0.2 млрд лет
2
T◦ = 5637 K, T◦s = 5372 K
1180.0
54.4
6.5
1.53
0.16
Ts = 273 K, Cme = 1.16 × 10-14 K с-1
ν◦ = 1.33 × 107 м2/с, A◦ = 7.42 × 104 K
эволюции Qb превышало современное значение в
ли и минимальную дисперсию теплового потока
2.5 раза. В то же время максимальные значения
на границе ядро-мантия. Обратим внимание, что
Qb в наших расчетах уже значительно меньше
при сравнительно том же размере твердого ядра
наблюдаемых в (Стивенсон и др., 1983), где от-
его возраст увеличился на четверть. Нам удалось
личие от современного значения составляло два
существенно уменьшить Qs, хотя оно все еще оста-
порядка, но также может быть еще уменьшено. Для
ется больше требуемого значения. Попытки умень-
поиска оптимального решения была использована
шить Qs еще более приводят к уменьшению ради-
штрафная функция Ψb (см. результаты расчетов
уса твердого ядра c, вплоть до его исчезновения.
в табл. 3, режим 2). Напомним, что для данного
В то же время нам удалось уменьшить отношение
расчета набор варьируемых параметров совпадает
σb/Mb более чем в три раза, уменьшив вариации Qb
с набором в 1a, но дополнительно к ограничению
практически на всем рассматриваемом интервале
на размер твердого ядра добавлены условия на
времени. Небольшие значения Qb в начальный
амплитуду теплового потока на поверхности Зем-
момент времени связаны с выбором однородного
ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 47
№7
2021
532
РЕШЕТНЯК
распределения температуры в мантии в начальный
В заключение остановимся на поведении тепло-
момент времени. Далее мы рассмотрим поведение
вых потоков (рис. 2г). За 4.5 млрд лет мощность
во времени ряда важных величин, продлив интер-
радиоактивных источников сократилась в 4.5 раза,
вал моделирования в будущее.
а тепловой поток Qs — в 6.8 раза. Скорость
убывания Qs в течение первого миллиарда лет
Скорость роста радиуса твердого ядра c замед-
была намного выше скорости убывания H. Это
ляется со временем и за последующие 1.5 млрд лет
было вызвано быстрым увеличением вязкости при
c увеличится лишь на 135 км (т.е. на 11%) (рис. 2а).
уменьшении температуры мантии (рис. 2б, в). Да-
Скорее всего подобное утоньшение жидкого ядра
лее уменьшение Qs стало заметно слабее, чем у H.
не отразится существенно на генерации магнитного
Тепловой поток на границе ядро-мантия
поля. Другое изменение геометрии связно с эволю-
Qb достиг своего максимума в течение первых
цией тепловых пограничных слоев. Для настоящего
100 млн лет и к настоящему времени медленно
времени их толщины составляют δ1 = 83 км и δ2 =
убывает со скоростью порядка 3-5% за 1 млрд лет.
= 137 км. Это близко к современным оценкам —
Для момента времени t = 4.5 млрд лет Qb =
порядка 100 км для обоих слоев (Дзевонский,
= 6.5 ТВт, укладываясь в диапазон значений из
Андерсон, 1981; Стейси, Лопер, 1983). Толщины
работ (Накагава, Тэкли, 2012; Накагава, 2018)
слоев растут. Появление твердого ядра несколько
6-15 ТВт. Как уже отмечалось выше, полученные
замедлило рост δ1, поскольку за счет выделения
значения потока на поверхности планеты Qs
дополнительной энергии в ядре процесс охлажде-
несколько выше наблюдаемых. При уменьше-
ния мантии также замедлился, приведя к излому
нии Qs твердое ядро становится меньше. Для
кривой δ1. На поведении δ2 появление твердого
уменьшения Qs мы несколько уменьшили вклад
ядра не отразилось: процессы в ядре очень сильно
радиоактивного разогрева в мантии H (Стивенсон
экранируются мантией, если не считать наблюда-
и др., 1983; Шуберт и др., 2001), характеризуемого
емое на поверхности Земли магнитное поле. Через
безразмерными числами Ur1 = H/Qs = 0.21 или
1.5 млрд лет δ1, δ2 достигнут значений 90 и 150 км
Ur2 = H/(Qs - Qb) = 0.25.
соответственно.
на всем вре-
Постоянство теплового потока Qb
Скачки температуры δT1, δT2 в пограничных
менном интервале может оказаться крайне важ-
слоях в настоящее время в модели равны 850 и
ным, поскольку палеомагнитологи не видят суще-
2600 K (рис. 2в). Величина первого скачка близка к
ственных изменений напряженности геомагнитного
наблюдаемой в слое D′′ 840 K (Дзевонский, Андер-
поля, начиная с 3.5 млрд лет в прошлом. Этот
сон, 1981). Величина второго скачка включает ска-
интервал времени, скорее всего, включает и мо-
чок в литосфере и мантии. Вычислим температуру
мент появления твердого ядра. Отсутствие резкого
границы ядро-мантия Tb = Ts + δT1 + δT2 = 273 +
изменения в поведении Qb в момент появления
+ 850 + 2600 = 3723 K, что близко к современным
твердого ядра в рассмотренной выше модели мож-
но расценивать как согласие с палеомагнитными
оценкам ∼3800 K (Трубицын, Трубицын, 2020).
наблюдениями. Оптимизм внушает и тот факт, что
Обратим внимание на разный знак производной по
генерация магнитного поля до и после появления
времени для обоих слоев: δT1 — растет, δT2 —
твердого ядра при заданном монотонном уменьше-
уменьшается во времени. Это связано с тем, что
нии Qb во времени позволяет найти режимы гене-
мантия остывает быстрее, чем ядро, и Ts не зависит
рации, когда появление твердого ядра не приводит
от времени. Скачок температуры в верхнем погра-
к существенному изменению энергии, доступной
ничном слое в три раза больше скачка в нижнем
для генерации магнитного поля (Решетняк, 2021).
слое (и это при том, что площадь нижнего слоя в 4
раза меньше верхнего!). Различие скачков связано
с меньшим значением вязкости на нижней границе
ЗАКЛЮЧЕНИЕ
мантии по отношению к верхней границе. Совре-
Совместное моделирование остывания ядра и
менные значения вязкости, полученные в модели,
мантии — крайне интересная физическая задача.
ν = 1.6 × 1018 м с-2 выше известных оценок 3 ×
Различие физико-химических свойств ядра и ман-
× 1017 м с-2 (Трубицын, 2016). Это является неиз-
тии приводит к неожиданным явлениям. Если для
бежной платой при введении средней температуры
моделирования остывания ядра достаточно знать
Tm по всей толщине мантии. За время эволюции
начальные условия и зависимость теплового пото-
вязкость изменилась более чем на 2 порядка. На
ка на границе ядро-мантия от времени, чтобы вы-
порядок уменьшилась конвективная теплоотдача,
яснить распределение температуры в произволь-
см. график для числа Нуссельта (рис. 2в). Для
ный момент времени, то для совместного реше-
древних эпох это число было больше 100, т.е. Ra ∼
ния уравнений остывания мантии и ядра ситуация
∼2×109, чтоеще достаточносложно для трехмер-
становится намного сложнее. Как мы видели, ха-
ных вычислений.
рактерное время памяти в модели мантии порядка
ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 47
№7
2021
ПАРАМЕТРИЧЕСКАЯ ТЕПЛОВАЯ МОДЕЛЬ ЭВОЛЮЦИИ ЗЕМЛИ
533
(а)
1500
150
c
m
1
m
2
1000
100
500
50
0
0
0
2
4
6
t, млрд лет
(б)
4000
4000
T
m
T
1
T2
3000
3500
2000
3000
1000
2500
0
0
2
4
6
t, млрд лет
(в)
103
1018
1017
m
Nu
2
10
1016
1015
1014
101
0
2
4
6
t, млрд лет
(г)
450
45
Q
s
Qb
H
300
30
150
15
0
0
0
2
4
6
t, млрд лет
Рис. 2. Зависимость от времени (а) радиуса твердого ядра c и толщин тепловых пограничных слоев δm1, δm2;
(б) температуры мантии Tm и скачков температуры в пограничных слоях δT1, δT2; (в) кинематической вязкости мантии
νm и числа НуссельтаNu; (г) тепловогопотока на поверхностиЗемлиQs, потока на границеядро-мантияQb и выделения
тепла радиоактивными источниками H.
ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 47
№7
2021
534
РЕШЕТНЯК
миллиарда лет. На больших временах информа-
6.
Дзевонский, Андерсон (A.M. Dziewonski and
ция о начальных условиях теряется. Происходит
D.L. Anderson), Phys. Earth Planet. Inter. 25, 297
это в первую очередь за счет экспоненциальной
(1981).
зависимости вязкости вещества мантии от темпе-
7.
Дрейк (M.J. Drake), Geochim. Cosmochim. Acta.
ратуры. Нахождение распределений температуры,
64, 2363 (2000).
тепловых потоков в далеком прошлом, т.е. реше-
8.
Кристенсен и др. (U.R. Christensen, J. Aubert, and
ние обратной задачи, становится некорректным:
G. Hulot), Earth Planet. Sci. Lett. 296, 487 (2010).
небольшие изменения в значениях параметров и
9.
Лаброзе (S. Labrosse), Phys. Earth Planet. Inter.
современных наблюдениях приводят к большим
140, 127 (2003).
флуктуациям в начальных условиях, распределе-
10.
Лаброзе и др. (S. Labrosse, J.-P. Poirier, and J.-L. Le
ниях полей в далеком прошлом. Выходом является
Mou ¨el), Phys. Earth Planet. Inter. 99, 1 (1997).
использование как можно большего числа огра-
11.
Ли (C.-H. Lee), Galaxies. 6, 51 (2018).
ничений в моделях. Так, в моделях геодинамо для
12.
Лопер (D.E. Loper), Advanc. Geophys. 26, 1 (1984).
нахождения близкого к земному режима генерации
13.
Мак-Элхини, Сенаяке (M. McElhinny and
магнитного поля использовалось порядка десяти
W. Senanayake), J. Geophys. Res. 85, 3523 (1980).
критериев отбора (Кристенсен и др., 2010). Нечто
14.
Накагава, Тэкли (T. Nakagawa and P.J. Tackley),
подобное потребуется и в модели остывания Земли.
Earth Planet. Sci. Lett. 329, 1 (2012).
Использование методов Монте-Карло для оптими-
15.
Накагава (Т. Nakagawa), Phys. Earth. Planet. Int.
276, 172 (2018).
зации решения крайне полезно. С одной стороны,
иногда можно очень быстро получить выразитель-
16.
Пурье, Тарантола (J.-P. Poirier and A. Tarantola),
Phys. Earth Planet. Inter. 109, 1 (1998).
ный результат. С другой стороны, последователь-
ное увеличение количества ограничений позволя-
17.
Решетняк (M.Yu. Reshetnyak), Magnetohydro-
dynam. 55, 175 (2019).
ет выявить значимые эффекты, что делает метод
18.
Решетняк (M.Yu. Reshetnyak), Russ. J. Earth Sci.
сравнимым по эффективности с аналитикой. Выше
20, ES5007 (2020).
мы увидели, насколько незначимый, с точки зрения
мантии, поток тепла на границе ядро-мантия мо-
19.
Решетняк М.Ю., Докл. РАН. Науки о Земле. 496,
176 (2021).
жет быть важным для выбора параметров вещества
20.
Решетняк М.Ю., Павлов В.Э., Геомагнетизм и
в мантии, если учесть термодинамику жидкого яд-
Аэрономия. 56, 117 (2016).
ра, появление твердого ядра, генерацию магнитного
21.
Рихтер (F.M. Richter), J. Fluid Mech. 89, 553 (1978).
поля. Используемые подходы могут быть полез-
ны для описания процессов эволюции планет как
22.
Соломатов (V.S. Solomatov), Treatise on
geophysics (Ed. G. Schubert) 9, 91 (2007).
Солнечной системы (Броуер и др., 2010), так и
23.
Сотин, Лаброзе (C. Sotin and S. Labrosse), Phys.
экзопланет (Ли, 2018). Отметим, что использова-
Earth Planet. Inter. 112, 171 (1999).
ние различных астрономических объектов крайне
24.
Стейси, Лопер (F.D. Stacey and D.E. Loper), Phys.
плезно и для развития самих моделей.
Earth Planet. Inter. 33, 45 (1983).
Работа выполнена при поддержке гранта РНФ
25.
Стивенсон и др. (D.J. Stevenson, T. Spohn, and
No. 19-47-04110.
G. Schubert), Icarus 54, 466 (1983).
26.
Тардуно и др. (J.A. Tarduno, R.D. Cottrell,
M.K. Watkeys, A. Hofmann, P.V. Doubrovine,
СПИСОК ЛИТЕРАТУРЫ
E.E. Mamajek, D. Liu, D.G. Sibeck, L.P. Neukirch,
1. Абе (Y. Abe), Lithos. 30, 223 (1993).
and Y. Usui), Science 327, 1238 (2010).
2. Броуер и др. (D. Breuer, S. Labrosse, and T. Spohn),
27.
Трубицын В.П., Физика земли 5, 3 (2016).
Space Sci. Rev. 152, 449 (2010).
28.
Трубицын А.П., Трубицын В.П., Докл. РАН. Науки
3. Букер, Стенгел (J.R. Booker and K.C. Stengel),
о Земле 495, 41 (2020).
J. Fluid Mech. 86, 289 (1978).
29.
Шуберт и др. (G. Schubert, D.L. Turcotte, and
4. Вангелов, Джарвис (V.I. Vangelov and G.T. Jarvis),
P. Olson), Mantle convection in the Earth and
J. Geophys. Res. 99, 9345 (1994).
5. Джарвис (G.T. Jarvis), J. Geophys. Res. 98, 4477
planets (Cambridge: Cambridge Univer. Press,
(1993).
2001).
ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 47
№7
2021
