Thermal convection modeling of the evolution of the earth core

Введение

Начало кристаллизации и рост твёрдого ядра $–$ одно из важнейших событий эндогенной геологической истории нашей планеты. Оно определило смену интенсивности магнитного поля Земли [1, 2], величину теплового потока на границе жидкого ядра и мантии (СМВ) [6], а также, вследствие фракционирования лёгких элементов между твёрдым и жидким ядром на их границе (ICB) [5], характер взаимодействия жидкого ядра с веществом мантии [4, 6].

Современные оценки возраста (и, соответственно, скорости роста) твёрдого ядра сильно разнятся от >2.5 млрд лет [2, 11] до 0.5 млрд лет [3], в зависимости от метода оценки (палеомагнитные данные об интенсивности магнитного поля Земли или термодинамические расчёты), от граничных условий (адиабатическая Т, постоянный тепловой поток на СМВ) и значений физических параметров (в первую очередь, теплопроводности ядра), используемых в расчётных моделях.

Большая часть работ, посвящённых процессам в ядре, связана с моделированием генерации магнитного поля [1]. Для этого используется система уравнений термической конвекции в приближении Буссинеска с учётом силы Кориолиса, обусловленной вращением Земли, и магнитной индукции (например, [14] и ссылки там), взаимодействующих через силу Лоренца. Современные модели в полную систему уравнений включают дополнительно химическую составляющую конвекции, обусловленную фракционированием элементов на границах СМВ и ICB [6]. Очевидно, что происходящие в ядре и на его границах процессы взаимосвязаны и взаимодействуют сложным нелинейным образом. При этом гидродинамика играет в них главную роль. Поэтому для адекватного понимания сути различных процессов целесообразно в первом приближении исследовать конвекцию в “чистом виде”, ограничившись только главными факторами, существенно влияющим на течение, а именно $–$ термической конвекцией и вращением (т.е. без магнитного поля и химической конвекции). Решению этой задачи посвящена предлагаемая работа.

Модель термической конвекции с кристаллизацией в приближении Буссинеска

Выберем характерные масштабы: $L,v_0,T_0,{\rho }_0,p_0,g_0,\nu ,\kappa ,\alpha ,{\Omega }^E,C_V,\text{H}$

$L,v_0,T_0,{\rho }_0,p_0,g_0,\nu ,\kappa ,\alpha ,{\Omega }^E,C_V,\text{H}$ (таблица).

Таблица. Физические величины, использованные в модельных расчётах

Запишем модель термической конвекции в безразмерном виде, используя принятые в гидродинамике критерии подобия. В приближении Буссинеска [20] модель термической конвекции включает:

уравнения Навье$–$Стокса

 $dV/dt=-eF{r}^{}{}^{-1}T-Eu\nabla p/r_0+R{e}^{-1}DV+R{o}^{-1}\left[e_{z}xV\right] +Cfr$ 

 $dV/dt=-eF{r}^{}{}^{-1}T-Eu\nabla p/r_0+R{e}^{-1}DV+R{o}^{-1}\left[e_{z}xV\right] +Cfr$ (1)

уравнение теплопроводности

 $dT/dt=P{e}^{-1}\Delta T- P_{0}divV+CrdQcr/dt$(2)

и уравнение неразрывности

 $V_{r}d{r}_0/dr+r_{0}divV=0$. (3)

В этих уравнениях $\rho { }_0\left(r\right)$ и $P_0(r)$ $–$ известные распределения плотности и давления, модель PREM [8], $V_r$ $–$ радиальная компонента скорости. Учтено также, что в приближении Буссинеска малые возмущения плотности ( $\epsilon <<1$ ) линейно выражаются через температуру

 $\rho \approx {\rho }_0\left(r\right)(1-\epsilon T)$. (4)

Полное давление $P=P_0+p$, где $P_0>>p$, состоит из гидростатической _{$P_0$} и динамической $p,$ частей. Уравнения движения (1) записаны во вращающейся системе координат ( _{$e_z$} $–$ единичный вектор оси z, направленной по оси вращения Земли), они учитывают силу тяжести, характеризуемую числом Фруда $Fr=v_0{}^2/g_{0}L$, $\epsilon =\alpha \delta T$; силы давления, характеризуемые числом Эйлера $Eu=p_0/{\rho }_0{v}_0{}^2$; вязкие силы, характеризуемые числом Рейнольдса $Re=v_{0}L/\nu$; силы Кориолиса, характеризуемые числом Росби $Ro=v_0/2{\Omega }^{E}L$, и центробежные силы, характеризуемые коэффициентом $Cf=[{\Omega }^E]2L^2/v_0{}^2$.

В уравнение теплопроводности входит число Пекле $Pe=v_{0}L/\kappa$ и коэффициент, характеризующий тепло фазового перехода $Cr=H/\left(C_V{T}_0\right)$. Полагая $p_0={\rho }_0{g}_{0}L$, получим $Fr·Eu=1$ и исключим число Эйлера. Заметим, что у нас ${\rho }_0,g_0,L$ такие, что $p_0\approx p_{max}=360$ ГПа равняется давлению в центре Земли. Число Рэлея выражается через используемые критерии подобия $Ra=\epsilon F{r}^{-1}RePe$.

Отметим, что при большом линейном размере задачи, несмотря на неопределённость коэффициентов вязкости и теплопрово дности, числа Рейнольдса $Re$ и Пекле $Pe$ очень велики, больше 10¹⁰, так что в жидкой части ядра происходит интенсивная турбулентная конвекция фактически невязкой и нетеплопроводной жидкости. Вязкость и теплопроводность играют роль только в тонких пограничных слоях около твёрдых поверхностей. Заметим, что при численном моделировании появляется схемная диффузия, зависящая от шага сетки, поэтому мы будем вынуждены ограничиться числами Рейнольдса и Пекле $Re=Pe={10}^6$. Остальные коэффициенты, входящие в наши уравнения, равны: $Eu=F{r}^{-1}=g_{0}L/v_0{}^2=1.4\cdot {10}^6$, $\epsilon =\mathrm{0.01,} R{o}^{-1}=2{\Omega }^{E}L/v_0=102$.

Начальные значения и граничные условия

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

Условия на границе с мантией $—$ это условие прилипания для вектора скорости и условие для температуры. Будем считать температуру одинаковой во всех точках внешней границы и медленно убывающей по экспоненциальному закону со временем:

 $T_{\text{CMB}}\left(t\right)=T_{\text{CMB}}\left(0\right)\cdot \text{exp}(-{\lambda }_{M}t)$. (10)

Тогда темп охлаждения будет регулироваться показателем ${\lambda }_M$. Ясно, что мантия, оправдывая своё название, экранирует выход тепла из ядра и замедляет остывание планеты, поэтому показатель ${\lambda }_M$ должен быть малой величиной.

Чтобы обойти проблему начальных условий, сначала проведём вспомогательный расчёт. Зададим начальное состояние покоя с малыми случайными возмущениями температуры, тогда в силу неустойчивости такого состояния в системе возникает движение, которое начнет раскручиваться по законам термической конвекции. Проведём расчёт до стабилизации среднестатистического течения, т.е. до того состояния, когда начальные условия забудутся. Полученные таким образом распределения скоростей $V\left(r\mathrm{,0}\right)$ и температуры $T\left(r\mathrm{,0}\right)$ примем за начальное состояние для нашего моделирования.

Для моделирования кристаллизации ядра требуется задать зависимость температуры плавления вещества ядра от гидростатического давления. Дефицит плотности и скорости продольных волн ядра, установленный сейсмологией [8] указывает на присутствие в его составе лёгких элементов. Одним из наиболее часто предлагаемых лёгких элементов в ядре Земли является водород [10, 12, 19]. Температура плавления водородсодержащего железа сильно зависит от состава FeH_x, где х $–$ количество атомов Н в формульной единице ( $0<x<1.2$ ). При х = 1 (состав FeH) кривая плавления (Т $–$ температура К, Р $–$ давление, Гпа) описывается уравнением [10]:

$T{1}_{melt}=T_0((P-P_0)/a+1)/c$, (11)

где

$a=\mathrm{24.6,}c=\mathrm{3.8,}{T}_0=1473\text{K}$ и $P_0=9.5$ Гпа.

Уравнение кривых плавления при переменном содержании Н получено нами путём интерполяции между кривой при х = 1 (уравнение 11) и кривой плавления чистого Fe по [13]:

$T_{\text{melt}}=T{1}_{\text{melt}}+\left(d+e\cdot P\right)\cdot (1-х),\text{где}d=\mathrm{521.7391,}e=4.7826$ (12)

где

$T_{\text{melt}}=T{1}_{\text{melt}}+\left(d+e\cdot P\right)\cdot (1-х),\text{где}d=\mathrm{521.7391,}e=4.7826$.

Согласно оценке [7], температура на современной границе СМВ с вероятностью 95% лежит в диапазоне 3470$–$3880 °К. Поэтому для моделирования процесса кристаллизации взята кривая плавления с содержанием водорода х = 0.5 (0.9 мас. % Н), попадающая в указанный диапазон $T_{\text{CMB}}$, а также хорошо согласующаяся с оценкой плотности и скорости продольных волн в ядре [12]. Отметим, что на этом этапе моделирования мы не рассматривали фракционирование лёгкого элемента между твёрдым и жидким ядром.

О численном методе и 2D-модели

В качестве инструмента нашего исследования используем 2D-вариант нашей модели термической конвекции в “плоском ядре”, плоскость которого ортогональна оси вращения Земли. Конечно, для реального моделирования, в частности генерации магнитного поля, необходимы полномасштабные 3D-расчёты. Но в данной работе для установления основных свойств термической конвекции используется более простая и понятная, благодаря полной визуализации, 2D-модель, к тому же более быстрая и точная в вычислительном отношении.

Численное моделирование осуществлялось конечно-разностным методом с помощью аппроксимации второго порядка точности уравнений в частных производных на равномерных декартовых сетках. Расчёты выполнены на сетках размером 512×512 и 1024×1024 узлов, что достаточно для прямого численного моделирования турбулентных режимов и учета больших чисел Рейнольдса и Пекле.

Конвекция в полностью жидком ядре. В геофизической литературе принято исследовать и определять радиальные распределения параметров в Земле в её современном состоянии, соответствующем модели PREM. Для температуры ядра вычисляется и приводится её адиабатическое (или изоэнтропийное) распределение. Например, в классической монографии В.Н. Жаркова [9] используется соотношение $T_{\text{ad}}={\rho }^{\gamma }$ с параметром Грюнайзена $\gamma \approx 1.45$. Так как при таком вычислении температуры конвекция в ядре не учитывается, то остаётся вопрос о её реальном распределении. Получаемая в нашей конвективной 2D-модели температура ядра зависит от времени и координат; в полярных координатах это функция $T(t,r,\varphi )$. Для получения распределения температуры в принятом радиальном виде достаточно осреднить её по угловой координате: $T_{ñð}\left(t,r\right)=⟨T(t,r,\varphi )⟩$.

На рис. 1 в естественной цветовой шкале слева показаны отклонения температуры от её среднего распределения, а справа Z-компонента завихренности. Оттенки красного цвета здесь и далее представляют положительные значения отклонения, синего цвета $–$ отрицательные, белый цвет соответствует около нулевым значениям (при превышении максимальных значений цвета зацикливаются).

Рис. 1. Температура (Т) и завихренность (Ω) в ядре до начала кристаллизации

Видно, что на границе с мантией образуются тонкие пограничные слои $–$ холодный тепловой и вязкий гидродинамический. Тяжёлый холодный слой часто отрывается и начинает тонуть в виде множества мелких струй. Тонкие струйки сливаются и укрупняются, и это повторяется каскадно, в итоге остаётся небольшое число быстро погружающихся струй, которые по инерции проскакивают область невесомости в центре и снова движутся к мантии, при этом максимальные скорости достигаются в центральной части (видео 1 в приложении; дополнительные материалы размещены в электронном виде по DOI статьи и на сайте редакции). На обеих панелях видна характерная особенность конвекции в полностью жидком ядре: вещество не скапливается в центральной части, а по инерции проходит её и снова выходит в периферийную область.

Мелкие вихри, сливаясь, укрупняются до размера, сопоставимого с радиусом ядра. В связи с этим следует отметить существующие в литературе расхождения: некоторые исследователи (например, [4]) отмечают преимущественно крупномасштабный организованный ламинарный характер конвекции, а другие [6, 15] $–$ мелкомасштабный турбулентный. Результаты нашего моделирования показывают, что оба режима реализуются одновременно. Формирующиеся при этом крупные округлые вихри являются аналогами вихревых столбов Тейлора, с образованием которых связывается генерация дипольного магнитного поля [15].

Моделирование кристаллизации. Моделирование процесса кристаллизации на мелких сетках проводилось следующим образом. Для каждого узла известна температура плавления, зависящая через гидростатическое давление от радиуса. Узлы сетки разделялись на “жидкие” и “твердые” и для описания их текущего состояния вводилась булевская функция. В начале моделирования все узлы $–$ “жидкие”. Затем в каждый момент времени в каждом узле осуществляется проверка и, если в жидком узле температура опустилась ниже температуры плавления (с небольшим переохлаждением), то этот узел переходит в категорию твёрдых узлов. Аналогично, если в твёрдом узле температура повышается выше температуры плавления (с небольшим перегревом), то этот узел возвращается в категорию жидких узлов. Каждый переход сопровождается выделением или поглощением энергии фазового перехода и скачком плотности. Увеличение плотности твёрдой фазы способствует её перемещению и сосредоточению в центральной части ядра. Кроме того, движение твёрдой фазы должно быть вращательно-поступательным и связанным с течением жидкой фазы условиями прилипания.

На рис. 2 зелёным цветом последовательно показаны стадии кристаллизации ядра. Видно, что кристаллизация начинается в центральной области с появления отдельных центров кристаллизации, вокруг которых растут кристаллизующие участки, затем сливающиеся друг с другом, образуя сплошную затвердевшую область.

Рис. 2. Последовательные стадии кристаллизации ядра от ранних (панель а) к поздним (панель м). Современному размеру твёрдого ядра соответствует панель и

Рис. 3. Температура (Т) и завихренность (Ω) в жидком ядре современной конфигурации

Более детально процесс кристаллизации, соответствующий нашей численной модели, показан в видео 2 Приложения. Кристаллизующееся вещество становится не абсолютно твёрдым, а только очень вязким. Его плотность, а вместе с ней сила тяжести, также немного увеличиваются. Поскольку в центральной области имеет место невесомость, кристаллизующееся вещество по инерции может пересечь область кристаллизации, выйти из неё в область низкого давления и снова расплавиться. Однако по мере дальнейшего остывания ядра радиус зоны кристаллизации увеличивается, и затвердевающие фрагменты начинают быстро накапливаться в центральной области. При этом между ними остаются жидкие прослойки, так что на этом начальном этапе кристаллизации в центральной области наблюдается рыхлая пористая структура [21]. Позже жидкие прослойки также кристаллизуются, но происходит это в условиях “невесомости”. По мере роста твёрдого ядра происходит перестройка структуры конвекции, большие струи не могут проходить через центр, они разворачиваются, в том числе и из-за тепла, выделяющегося при кристаллизации и увеличивающего силу плавучести. Вихревые структуры, закручиваемые струями, уменьшаются в размерах, вследствие чего число струй и вихревых структур постоянно возрастает. Дальнейший рост внутреннего ядра происходит неравномерно в тех местах, где подходят холодные струи. Так как конвекция хаотичная, то струи подходят в разных местах, и форма кристаллической части ядра становится всё более круглой. Ко времени, когда внутреннее ядро достигает современного размера R = 0.35 (панель и на рис. 2), его форма становится почти круглой (панели к, л, м на рис. 2; видео 2 приложения).

Конвекция в ядре современной конфигурации. Конвекция в ядре современной конфигурации (современный радиус твёрдого ядра 1221.5 км составляет 0.35 от радиуса всего ядра) исследуется наиболее часто. В частности, имеются аналогичные исследования в чисто термической постановке, без учёта магнитного поля [17]. Результаты моделирования термической конвекции в современном ядре показаны на рис 3; видео 2 приложения. В отличие от полностью жидкого ядра (см. рис. 1), к холодным нисходящим потокам добавляются горячие восходящие потоки. Видно, что максимальный размер вихрей уменьшается.

Интегральные результаты моделирования

2D-моделирование термической эволюции ядра, с условием экспоненциально убывающей температуры на границе ядро/мантия, позволяет находить полные распределения температуры, скоростей и конфигураций фаз во все моменты времени. Видеозаписи численных экспериментов даны в приложении (см. электронное приложение). Для анализа результатов, получаемых в каждый момент времени, дополнительно вычислялись осреднённая по пространству скорость конвекции $V_{av}$, осреднённое вдоль CMB число Нуссельта $Nu$, представляющее безразмерный тепловой поток, и приведённая к радиусу круга $R_{ic}$ площадь твёрдой фазы (рис. 4 а). Распределения приведены в относительных переменных, реперной точкой является современное состояние t = 0, $R_{ic}$ = 0.35.

Рис. 4. (слева) – рост внутреннего ядра ( $R_{ic}$, зелёный цвет), тепловой поток из ядра в мантию ( $Nu$, синий) и средняя скорость конвекции ( $V_{av}$, коричневый); (справа) – усреднённый профиль температуры в ядре ($T_{conv},$ сплошная красная кривая) в сравнении с адиабатическим по [9] ($T_{ad},$ пунктир)

В наших расчётах принято, что температура CMB экспоненциально убывает со временем (чёрная кривая на рис. 4 слева). Выход тепла из ядра, то есть времена его остывания и кристаллизации, контролирует мантия. Для того чтобы перевести результаты в размерное время, нужна дополнительная, не связанная с ядром, информация о теплоизолирующем влиянии мантии. Например, в работе [18] утверждается, что постепенное охлаждение Земли составляет около 100 °С за миллиард лет. В работе по моделированию мантийной конвекции [16], где первые 0.5 млрд лет отводятся на кристаллизацию самой мантии, показано, что за последующие 4 млрд лет геологической эволюции температура ядра уменьшилась на 12.5%, что созвучно оценке [18]. Значит 4 млрд лет тому назад она равнялась 4434 °К. Исходя из такого темпа остывания $T_{\text{CMB}}$, согласно рис. 4 получаем, что для достижения современного радиуса твёрдого ядра его кристаллизация должна была начаться примерно 0.5 млрд лет тому назад.

Представленные на левом рис. 4 результаты численного эксперимента показывают, что с момента появления твёрдого ядра тепловой поток из ядра в мантию увеличивается вследствие выделения тепла кристаллизации, а скорость конвекции начинает убывать, но происходит это постепенно по мере роста твёрдого ядра, служащего препятствием для конвекции. Также виден хаотически осциллирующий характер процессов тепломассопереноса.

Из моделирования термической конвекции в жидком ядре следует, что развитая турбулентная конвекция происходит с очень высокими скоростями, ориентировочно v₀ ~5$–$7 м/с и сопровождается коротковолновыми осцилляциями как момента инерции ядра, так и его вращательного момента. Эти осцилляции должны выливаться в компенсирующие осцилляции угловой скорости вращения мантии, которые регистрируются на поверхности планеты [22].

Распределение осреднённой по угловой координате температуры в современном жидком ядре показано на рис. 4 справа. Видно, что в результате интенсивной конвекции распределение температуры более пологое по сравнению с адиабатическим, которое чаще всего используется в литературе [9]! И эта температура убывает со временем.

Выводы

Проведённое моделирование чисто термической конвекции позволило выявить некоторые важные особенности эволюции процессов, происходивших в ядре Земли на фоне остывания планеты.

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

2. Быстрый хаотичный рост твёрдого ядра в начальной стадии кристаллизации.

3. Бесформенная конфигурация ядра на начальной стадии кристаллизации и его рыхлость естественным образом объясняются отсутствием в центре силы тяжести.

4. При появлении твёрдого ядра, перекрывающего конвективные потоки через центр, начинается перестройка структуры конвекции, средняя скорость конвекции уменьшается. Но тепловой поток из ядра в мантию при этом увеличивается из-за выделения тепла кристаллизации.

5. Усреднённый профиль температуры в жидком ядре отличается от адиабатического.

Предложенная в нашей работе модель не включает ряд важных процессов, которые могут существенно влиять на тепломассоперенос и, соответственно, характер конвекции в ядре. Это, в первую очередь, фракционирование лёгкого элемента (водорода) между твёрдым и жидким ядром, которое порождает важную химическую составляющую конвективных потоков. Во-вторых, наше моделирование не учитывает электромагнитые силы Лоренца. Кроме того, оно, как и все предыдущие известные нам работы, проведено в предположении, что полный размер ядра соответствует современному, т.е. не менялся во времени, т.е. не учитывает возможность химического обмена металлом и лёгким элементом на границе внешнее ядро$–$мантия. Решение этих проблем $–$ задача будущих исследований.

Источники финансирования

Исследования выполнены в рамках гос. задания ИГЕМ РАН.