Одномерная модель вертикального переноса химических составляющих в атмосфере Марса вплоть до высот термосферы

Full Text

ВВЕДЕНИЕ

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

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

Одной из основных трудностей, возникающих в моделировании аэрономических процессов в планетной атмосфере, является выбор достаточно обоснованных значений такой ее важной динамической характеристики, как коэффициент $K$ вертикальной турбулентной диффузии для компонент. При изучении работ по этой тематике были определены ключевые модели, описывающие вклад турбулентной диффузии в стратосфере и термосфере Марса.

Основной мотивацией для разработки модели атмосферы Марса, работающей вплоть до высот термосферы, является необходимость понимания трансформации химической составляющей атмосферы Марса в верхней атмосфере, включающей заряженную компоненту (ионосферно-тропосферная система, известная как I/T-система). Для решения такой задачи требуется модель, хорошо согласующаяся с I/T-системой, в которой определяющими являются внутренние гравитационные волны. Поэтому для параметризации турбулентной вязкости естественно использовать модель, в которой турбулентность и турбулентная вязкость связаны с нестабильным разрушением приливов и с нелинейным взаимодействием внутренних гравитационных волн. Представляется также необходимым провести сравнение с традиционными моделями с турбулентной диффузией.

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

ИСХОДНЫЕ УРАВНЕНИЯ. ОПИСАНИЕ ИЗВЕСТНЫХ МОДЕЛЕЙ

С момента получения первых данных об атмосфере Марса были предложены различные способы построения теоретических моделей вертикального переноса химических составляющих. Рассмотрим основные из них.

В качестве исходной используется система уравнений неразрывности (Куликов, 2018), описывающая процесс вертикального массопереноса в приближении горизонтально однородной атмосферы. Она имеет следующий вид:

$\begin{array}{c}\left[\frac{\partial }{\partial t}+\frac{\partial }{\partial z}{w}_i\left(z,t\right)\right]n_i\left(z,t\right)=Q_i\left(z,t\right)\end{array}.$ (1)

Здесь t – время, z – высота над средним уровнем поверхности, $w_i$ и $n_i$ – вертикальная составляющая скорости и концентрации i-го компонента химического состава, $Q_i$ – скорость образования и уничтожения i-го компонента в результате всех химических реакций, включая процессы фотодиссоциации. Приведем явный вид “источникового” члена $Q_i$, связанного со скоростью появления и исчезновения частиц сорта $i$ в химических реакциях, включая процессы диссоциации (Хайрер, Ваннер, 1999). Существуют два типа реакций: одна, в которой возникает i-компонент, и вторая, в которой он расходуется ($n_{+}, n_{-}=n_i$, просто стоят в разных частях):

$k_a\times a+k_b\times b=k_i\times n_{+}+k_d\times d,$

$k_i\times n_{-}+k_d\times d=k_a\times a+k_b\times b.$

Тогда $Q_i$ может быть выражена как

$\begin{array}{c}{Q}_i\left(z,t\right)=\frac{\text{d}{n}_{i_{+}}}{\text{d}t}-\frac{\text{d}{n}_{i_{-}}}{\text{d}t}=k_i\times v_{+}-k_i\times v_{-}=\\ =\sum J_{ab}\left(T\right)n_a^{k_a}{n}_b^{k_b}-\sum J_{id}\left(T\right)n_i^{k_i}{n}_d^{k_d},\end{array}$

где $n_a, n_b, n_i, n_d$ – концентрации компонент в реакции, $a$ – доля компоненты $n_i$ в продуктах реакции, $J_{ab}\left(T\right)$ – константа реакции, приведенная в таблице.

 

Система химических реакций, рассматриваемая в данной работе

№	Реакция			Константа реакции J
1	CO2 + hν(40 − 2062)	→	CO + O	4 × 10−7
2	O2 + hν(40 − 2062)	→	O + O	1.8 × 10−6
3	H2O + hν(1515 − 2000)	→	OH + H	1.46 × 10−6
4	HO2 + hν(2000 − 3600)	→	O + OH	5.14 × 10−4
5	H2O2 + hν(2000 − 3600)	→	OH + OH	7.12 × 10−5
6	CO + O + CO2	→	CO2 + CO2	1.0 × 10−32exp(−2190/T)
7	CO + OH	→	CO2 + H	2.1 × 10−32exp(−116/T)
8	O + O + CO2	→	O2 + CO2	3.0 × 10−33exp(−300/T)2.9
9	O + O2 + CO2	→	O3 + CO2	1.4 × 10−33exp(−300/T)2.5
10	O + O3	→	O2 + O2	2.0 × 10−11exp(−2280/T)
11	O + H2	→	OH + H	7.0 × 10−11exp(−5000/T)
12	O + OH	→	O2 + H	4.0 × 10−11
13	O + HO2	→	OH + O2	8.0 × 10−11exp(−500/T)
14	O + H2O2	→	H2O + O2	2.3 × 10−11exp(−3200/T)
15	O + H2O2	→	HO2 + O2	2.3 × 10−11exp(−3200/T)
16	O(1D) + CO2	→	O(3p) + CO2	1.8 × 10−10
17	O(1D) + H2O	→	OH + OH	3.0 × 10−10
18	O(1D) + H2	→	OH + H	3.0 × 10−10
19	H + HO2	→	OH + OH	4.2 × 10−10exp(−950/T)
20	H + HO2	→	H2O+ O	8.3 × 10−11exp(−500/T)
21	H + HO2	→	H2 + O2	4.2 × 10−11exp(−350/T)
22	H + O2 + CO2	→	HO2 + CO2	2.0 × 10−31exp(−273/T)1.3
23	H + H + CO2	→	H2 + CO2	1.0 × 10−32exp(−273/T)0.7
24	H + O3	→	OH + O2	2.6 × 10−11
25	OH + OH	→	H2O + O	1.0 × 10−11exp(−550/T)
26	OH + H2	→	H2O + H	3.6 × 10−11exp(−2590/T)
27	OH + O3	→	HO2 + O2	1.3 × 10−12exp(−956/T)
28	OH + HO2	→	H2O + O2	8.3 × 10−11exp(−500/T)
29	OH + H2O2	→	H2O + HO2	4.1 × 10−11T0.5exp(−600/T)
30	HO2 + HO2	→	H2O2 + O2	1.7 × 10−11exp(−500/T)
31	CO2++ H2	→	CO2H++ H	1.4 ×10−9
32	CO2H++e	→	CO2 + H	3.5 × 10−7

 

Для паров (Courtney, 1962):

$Q{V}_{\text{v}}=k_{\text{V}}{n}_{\text{C}}\left(n_{\text{VS}}-n_{\text{V}}\right),$

где $k_{\text{v}}=3\alpha m_{\text{w}}V/4r{\rho }_{\text{ice} }$ – “константа” скорости конденсации пара; $m_w$ – масса молекулы воды; $\alpha$ – коэффициент конденсации паров, называемый также коэффициентом аккомодации для процесса конденсации; $r$ – средний радиус частиц конденсата; ${\rho }_{ice}$ – плотность водяного льда; $n_C$ – концентрация конденсата, $n_V$ – концентрация паров воды. Концентрация насыщенного водяного пара $n_{VS}$ – над плоской поверхностью конденсата рассчитывалась по формуле Гоффа–Грэтча (Куликов, 2018).

Усредненный по времени поток i-го компонента в вертикальном направлении

$\begin{array}{c}{\psi }_i=n_i{w}_i+\overline{{n_1}^{\text{'}}{w_1}^{\text{'}}}\end{array}$ (2)

состоит из двух частей: молекулярная диффузия ${\psi }_i^m=n_i{w}_i$ и турбулентная диффузия ${\psi }_i^T=\overline{{n_i}^{\text{'}}{w_i}^{\text{'}}}$. Здесь $w_0$ среднемассовая скорость движения атмосферы, $w_i$ диффузионная скорость, а $n_i\text{'}$ и $w_i\text{'}$ – значения величины пульсаций концентрации и скорости i-го компонента, возникающих за счет турбулентности.

Чтобы описать каждую из величин ${\psi }_i^T$, можно воспользоваться соотношениями Стефана–Максвелла, которые хорошо подходят для описания газа, состоящего из многих компонент. Однако, как правило, ввиду того факта, что углекислый газ занимает почти 95% от всего химического состава Марса, используется следующее, так называемое приближение диффузии малой составляющей:

${\psi }_i^m=-D_i\left[\frac{\partial n_i}{\partial z}+\left(\frac{1}{H}+\frac{1+a_i}{T}\frac{\partial T}{\partial z}\right)n_i\right],$

${\psi }_i^T=-K\left[\frac{\partial n_i}{\partial z}+\left(\frac{1}{H}+\frac{1}{T}\frac{\partial T}{\partial z}-\frac{1}{\mu }\frac{\partial \mu }{\partial z}\right)n_i\right].$ (3)

Здесь $D_i$ коэффициент диффузии i-го компонента в $C{O}_2$, $H_i=RT/{\mu }_{i}g$ – высота однородной атмосферы, $R$ – универсальная газовая постоянная, $T$ – температура газа, $g$ – ускорение свободного падения, µ_i и ${\alpha }_i$ – молекулярный вес и коэффициент термодиффузии i-го компонента, H шкала высот однородно перемешанной атмосферы, а µ – ее средний молекулярный вес.

Основное отличие моделей турбулентной диффузии состоит в определении коэффициента K. Например, согласно выводам работы (Изаков, 1976), для удовлетворительного соответствия рассчитанного состава атмосферы Марса и данных измерений необходимо, чтобы $K=4\times {10}^{8 }с{м}^2{с}^{-1}$. Столь высокое значение $K$ требуется, с одной стороны, для того, чтобы обеспечить быстрый перенос в нижнюю атмосферу атомарного кислорода и ограничить тем самым скорость образования $O_2$ на уровне $~30$ км в реакции $\text{O}+\text{OH}\to {\text{O}}_2+\text{H.}$. С другой стороны, большие величины $K$ необходимы для достаточно быстрого переноса вниз атомарного водорода, образующегося в ионосфере, с тем чтобы обеспечить поддержание концентрации H на экзобазе ( 225 км) на уровне, не превышающем $3\times {10}^4{\text{ см}}^{-3}$. В работе (Hunten, 1974) предполагалось, что коэффициент $K$ зависит от высоты и может изменяться в пределах $\left(3\times {10}^{6 }{\text{см}}^2{\text{с}}^{-1}\right)<K<\left(2.8\times {10}^{10 }{\text{см}}^2{\text{с}}^{-1}\right)$. На основании результатов своих расчетов авторы (Hunten, 1974) пришли к выводу, что величина K на высотах 30–40 км является критическим параметром модели и должна составлять здесь не менее $7\times {10}^7 {\text{см}}^2{\text{с}}^{-1}$. Влияние высотной зависимости K на состав атмосферы Марса исследовалось также в (Изаков, 1976). Согласно результатам этой работы, уменьшение K приводит к падению концентраций $\mathrm{H}$ и ${\mathrm{H}}_2$ на уровне экзобазы.

Важный прогресс был достигнут в работе (Куликов, 2018), где использовалась следующая модель:

$\begin{array}{c}K\left(z\right)=K_1\text{exp}\left(S_1\left(z-z_m\right)-S_2{\left(z-z_m\right)}^2\right),\end{array}$ (4)

где $S_2=0$ при $z\le z_m$. Значения параметров были взяты следующие: $K_m={10}^{7.5} с{м}^2{\text{с}}^{-1}$, $z_m=107$ км, $S_1=5.83\times {10}^{-2} к{м}^{-1}$ и $S_2=3.89\times {10}^{-3} {\text{км}}^{-2}$.

В работе (Куликов, 2018) исследуются следующие нейтральные компоненты: CO₂, CO, O(³ P), O(¹ D), O₃, H₂O, OH, H, HO₂, H₂O₂, H₂, ионы ${\text{CO}}_2^{+} \text{\hspace{0.33em}è} {\text{\hspace{0.33em}CO}}_2{\text{H}}^{+}$, а также исследуется аэрозольные составляющая, образованная частицами конденсата водяного пара. В нашей работе мы будем учитывать компоненты в рамках той же самой аэрономии атмосферы Марса, с тем чтобы понять важную роль модели турбулентной диффузии в средней и верхней атмосфере. Схема химических реакций приведена в таблице.

В работе (Chaffin и др., 2017) предсказан повышенный выброс атомарного водорода из атмосферы Марса, вызванный присутствием воды на большой высоте. При моделировании процессов переноса газа в нижней атмосфере также учитывались химические реакции, показанные в таблице. До высоты  км коэффициент турбулентной диффузии был константой и численно равен $1.0\times {10}^{6 }с{м}^2{с}^{-1}$, после чего менялся с высотой следующим образом (Краснопольский, 1982):

$\begin{array}{c}{K}_2\left(z\right)=K_0{\left(\frac{P\left(z\right)}{k_{b}T\left(z\right)}\right)}^{-\frac{1}{2}},\end{array}$ (5)

где $K_0={10}^{13}{\text{\hspace{0.33em}см}}^2{\text{c}}^{-1}.$

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

ГРАВИТАЦИОННАЯ ТУРБУЛЕНТНОСТЬ. ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ МОДЕЛИ ВЕРТИКАЛЬНОГО ПЕРЕНОСА

Мы предлагаем использовать модель турбулентности Линдзена (Lindzen, 1981), которая хорошо описывает турбулентную диффузию для земной атмосферы, для улучшения результатов моделирования марсианской атмосферы на средних и больших высотах. Предполагается, что механизмом, вызывающим турбулентную диффузию, является нелинейное взаимодействие внутренних гравитационных волн различных типов. Известны два механизма разрушения внутренних гравитационных волн на больших высотах атмосферы Марса. Один из механизмов представляет собой увеличение амплитуды внутренней волны вследствие уменьшения давления атмосферы с высотой, а другой механизм – это резонансное взаимодействие ветра и внутренней волны при равенстве скорости ветра и групповой скорости внутренней волны. Поэтому такая модель может рассматриваться как адекватная альтернатива для предсказания трансформации химических компонент атмосферы Марса на высотах, включающих ионосферу. Разрыв волны, по-видимому, происходит на высоте около 85 км, что приводит к слою усиленной вихревой диффузии и вызванного волной ускорения, простирающемуся в области между 85 км и примерно 108 км, выше которой преобладает молекулярный перенос.

Ситуация в средних и высоких широтах более сложная. Гравитационные волны доминируют на этих высотах, и если такие волны возникают в тропосфере (что вероятно), то их фазовые скорости обычно будут варьироваться от нуля (что соответствует горным волнам и т. д.) до типичных скоростей тропосферного потока.

Для таких низких фазовых скоростей влияние среднего зонального распределения потока будет иметь решающее значение. Распределение потока будет эффективно определять, какие гравитационные волны (в зависимости от фазовой скорости) могут достигать мезосферы, и соответственно, амплитуды мезосферных гравитационных волн и их предельную высоту. Этот эффект, например, подразумевает, что высота, до которой учитывается гравитационная волна летом (70 км), будет большей, чем зимой (50 км). Кроме того, турбулентная диффузия и вызванные волнами ускорения больше не распространяются до некоторого уровня, где доминируют молекулярные процессы. Вместо этого они ограничены слоями, простирающимися до некоторого приблизительного критического уровня (где среднее значение скорости зонального ветра равно или почти равно фазовой скорости волны). Между критическими уровнями и термосферой ожидается резкий минимум вихревой диффузии. Вызванные волнами ускорения обеспечивают явный источник “трения”, необходимого для обращения вспять мезосферных сдвигов и, соответственно, для изменения температурного градиента от полюса к полюсу в мезопаузе.

Данная концепция приводит к следующему виду турбулентной диффузии на высоте 70–80 км:

$\begin{array}{c}{K}_3\left(z\right)=\frac{k{c}^4}{N^3{\left(1+l^2/k^2\right)}^{\frac{3}{2}}}\frac{1}{2H},\end{array}$ (6)

где k = 3 × 10^–6 м^–1, с = 466 м/с,

$\begin{array}{c}{\left(1+l^2/k^2\right)}^{\frac{1}{2}}~\mathrm{1.54,} H=RT/g, \\ N^2=g/T\text{Ã}, \text{Ã}=\left(\partial T/\partial z\right)+\left(g/c_p\right).\end{array}$

Целью данной работы является построение глобально усредненной аэрономической модели Марса и определение на ее основе приближенного вертикального распределения величины K.

В настоящей работе мы предлагаем новую формулу турбулентной вязкости для численного моделирования переноса газовых составляющих атмосферы Марса, в которой вместо формулы (5) используется выражение Линдзена (Lindzen, 1981) для коэффициента турбулентной диффузии (6). В работе проведено численное моделирование на основе уравнения (1) с использованием двух моделей турбулентной диффузии (4) и (6) для атмосферы Марса на высотах 0–200 км. Проведено подробное сравнение результатов моделирования на основе обеих моделей турбулентной диффузии. Построена разностная схема и заданы граничные условия для каждой из рассматриваемых компонент, а именно: CO₂, CO, O₃, H₂O, OH, H, HO₂, H₂O₂, H₂, Ar, N₂.

Уравнения химических реакций берутся из таблицы. Запишем уравнение (1) для двух моделей турбулентной вязкости:

$\begin{array}{c}\left\{\begin{array}{c}\frac{\partial }{\partial t}{n}_i +\frac{\partial }{\partial z} Ô=Q_i\\ Ô=-\left(D_i+K\left(z\right)\right)\left[\frac{\partial n_i}{\partial z}+\left(\frac{1}{H_i}+\frac{1}{T}\frac{\partial T}{\partial z}\right)n_i\right]-D_i\frac{{\alpha }_i}{T}\frac{\partial T}{\partial z}-K\left(z\right)\frac{1}{\mu }\frac{\partial \mu }{\partial z}.\end{array}\right.\end{array}$ (7)

Нам необходимо найти установившейся средний профиль концентраций химических компонент от высоты, соответственно: $\partial n_i/\partial t=0$. В итоге запишем уравнение (7) в следующем виде:

$\begin{array}{c}{Q}_i=-\left(\frac{\partial D_i}{\partial z}+\frac{\partial K}{\partial z}\right)\left[\frac{\partial n_i}{\partial z}-\left(\frac{1}{H_i}+\frac{1}{T}\frac{\partial T}{\partial z}\right)n_i\right]-\\ -\left(K+D_i\right)\left[\frac{{\partial }^{2}T}{\partial z^2}+\left(-\frac{1}{H_i^2}\frac{\partial H_i}{\partial z}+\frac{1}{T^2}{\left(\frac{\partial T}{\partial z}\right)}^2+\frac{1}{T}\frac{{\partial }^{2}T}{\partial z^2}\right)n_i+\left(\frac{1}{H_i}+\frac{1}{T}\frac{\partial T}{\partial z}\right)\frac{\partial n_i}{\partial z}\right]-\\ -\frac{\partial D_i}{\partial z}\frac{{\alpha }_i}{T}\frac{\partial T}{\partial z}+D_i\frac{{\alpha }_i}{T^2}{\left(\frac{\partial T}{\partial z}\right)}^2-D_i\frac{{\alpha }_i}{T}\frac{{\partial }^{2}T}{\partial z} -\frac{\partial K}{\partial z}\frac{1}{\mu }\frac{\partial \mu }{\partial z}\begin{array}{c}+D_i\frac{1}{{\mu }^2}{\left(\frac{\partial \mu }{\partial z}\right)}^2-K\frac{1}{\mu }\frac{{\partial }^2\mu }{\partial z^2}.\end{array}\end{array}$ (8)

Зависимость коэффициента молекулярной диффузии от температуры и, как следствие, от высоты для каждой компоненты вычисляется по формуле Чепмена–Энскога:

$\begin{array}{c}{D}_i=\frac{A{T}^{\frac{3}{2}}\sqrt{\frac{1}{{\mu }_0}+\frac{1}{{\mu }_i}}}{P\left(z\right){\sigma }_{\mathrm{1,2}}^2\Omega }.\end{array}$ (9)

Здесь 1, 2 – индексы двух видов молекул, присутствующих в газовой смеси, µ_i – молярные массы компонент (в нашей модели рассматривается диффузия компонент в газе СО₂ ввиду его преобладания в сравнении с остальными газами), коэффициент А = 1.859 × 10^(–3) атм · Ų · см²  · K^(–3/2) / с, P(z) – профиль давления (атм), ${\sigma }_{\mathrm{1,2}}=\frac{1}{2}\left({\sigma }_1+{\sigma }_2\right)$ – эффективный диаметр столкновения в Å (значения заданы в форме таблицы в [13]), $\Omega ~1$ – безразмерное значение интеграла столкновений.

На рис. 2 изображены используемые профили коэффициента молекулярной диффузии для некоторых компонент.

 

Рис. 2. Коэффициенты диффузии нескольких химических компонент атмосферы Марса, вычисленные по формуле Чепмена–Энскога. Использованные профили температур изображены на рис. 1. Параметры самих газов были взяты из (Hirschfelder и др., 1964).

 

Выражения $H_i\left(z\right)$ вычислялись для каждой компоненты отдельно по профилю температуры $T\left(z\right)$. Необходимые распределения температуры и давления взяты из (Chaffin и др., 2017; Краснопольский, 1982). Соответствующие профили изображены на рис. 1. Изменением средней молярной массы в (Lindzen, 1981) в первом приближении пренебрегаем, как уже упоминалось, ввиду того факта, что во всей нижней атмосфере Марса концентрация углекислого газа на порядок и больше превышает остальные показатели.

 

Рис. 1. Усредненные среднесуточные профили температуры (а) и логарифма давления (б) для среднеширотной атмосферы Марса. Данные взяты из статей (Bougher и др., 2015; 2017; Chaffin и др., 2017).

 

Система дифференциальных уравнений (1), в которой присутствуют быстрые и медленные текущие химические реакции, относится к типу жестких дифференциальных систем уравнений в форме Коши (Маров, Колесниченко, 1987). Такие системы дифференциальных уравнений сталкиваются с принципиальными трудностями при их численном решении. Прежде всего требуется резкое увеличение числа вычислений. Отметим также резкое возрастание погрешности вычислений при недостаточно малом шаге. Поэтому в качестве метода численного решения мы используем неявные методы, которые обеспечивают результат, несравненно более надежный, чем явные методы (Маров, Колесниченко, 1987). Фактически речь идет об обеспечении устойчивости численного алгоритма при использовании неявных схем. Использование неявных схем в задачах планетной аэрономии и в химии планетных атмосфер имеет определяющее значение, поскольку одни и те же системы уравнений при разных параметрах могут быть как жесткими, так и нежесткими.

Все производные в (8) аппроксимированы c помощью центральной разностной схемы:

$\begin{array}{c}\left(\frac{\partial f}{\partial z}\right)\left(z_k\right)=\frac{f_{k+1}-f_{k-1}}{2h},\end{array}$ (10)

$\begin{array}{c}\left(\frac{{\partial }^{2}f}{\partial z^2}\right)\left(z_k\right)=\frac{f_{k+1}-2f_k+f_{k-1}}{h^2}.\end{array}$ (11)

Здесь $h=z_{k+1}-z_k$. Таким образом, по полученным профилям давления и температуры мы находим соответственно зависимость от высоты следующих величин:

$\frac{{\partial }^{2}T}{\partial z^2},\frac{\partial T}{\partial z},\frac{\partial D_i}{\partial z},\frac{\partial K}{\partial z},\frac{\partial H_i}{\partial z}$.

В дальнейшем для удобства производные функции будем просто обозначать $f\text{'},f\text{'}\text{'}$. В результате получим уравнение непрерывности для i-й компоненты в разностной форме:

$\frac{n_{k+1}-2n_k+n_{k-1}}{h^2}\underset{A_k^2}{\underbrace{D_k^0}}+n_{k+1}-\underset{A_k^1}{\underbrace{\frac{n_{k-1}}{2h}\left[\left(\frac{1}{H_k}+\frac{{T_k}^{\text{'}}}{T_k}\right)D_k^0 + {D_k^0}^{\text{'}}\right]}}+\underset{A_k^0}{\underbrace{n_k\left[\left(\frac{1}{H_k}+{T_k}^{\text{'}}{T}_k\right){D_k^0}^{\text{'}} +D_k^0\left(-\frac{{H_k}^{\text{'}}}{H_k^2}-\frac{{T_k}^{\text{'}}}{T_k^2}+\frac{{T_k}^{\text{'}\text{}\text{'}}}{T_k}\right)\right]}}=\underset{B_k}{\underbrace{Q_k}},$

где $D_k^0=K_k+D_k^i$.

Обозначая соответствующие коэффициенты перед производными концентрациями $A_k^0,A_k^1,A_k^2,$ получаем следующую систему уравнений:

$\begin{array}{c}\begin{array}{c}{n}_{k+1}\underset{{\alpha }_1}{\underbrace{\left(A_k^2+\frac{h{A}_k^1}{2}\right)}}+n_k\underset{{\alpha }_2}{\underbrace{\left(-2A_k^2+h^2{A}_k^0\right)}}+\\ +n_{k-1}\underset{{\alpha }_3}{\underbrace{\left(A_k^2-\frac{h{A}_k^1}{2}\right)}}+B_k{h}^2=0.\end{array}\end{array}$ (12)

Откуда следует:

$\begin{array}{c}{n}_{k+1}=-\frac{1}{{\alpha }_1}\left(n_k{\alpha }_2+n_{k-1}\right).\end{array}$ (13)

Граничные условия на поток и концентрацию i-й химической компоненты при z = 0, брались из (Chaffin и др., 2017). В местах разрывов производной функций температуры, коэффициентов диффузии и других параметров граничные условия выбирались такими, чтобы соблюдалась непрерывность концентраций и потоков.

Выбранные граничные условия для некоторых компонент: распределение концентрации CO₂ определялось из условия гидростатического равновесия; у CO и O₂ потоки равны 0 (диффузионное равновесие). У H₂ и H₂Oи других паров начальные потоки вычислялись по методу из (Куликов, Рыхлецкий, 1983). На верхней границе (граница турбопаузы 225 км) принималось, что для всех компонент, кроме водородосодержащих, диффузионные потоки равны 0.

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

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

В работе проведено сравнение двух моделей, которые наиболее точно среди всех описывают итоговые распределения. Соответственно, в первой модели коэффициент турбулентной диффузии вычислялся по формуле (4), а во второй – по формуле (6). Профили коэффициентов диффузии изображены на рис. 3. Соответственно, можно видеть, что на высоте 120 км коэффициент молекулярной диффузии начинает преобладать.

 

Рис. 3. На графиках сплошной линией изображены коэффициенты турбулентной диффузии, вычисленные по формуле (4) для (а) и по формуле (6) для (б). Пунктирной линией изображены вычисленные ранее коэффициенты молекулярной диффузии для компоненты O. Штрихпунктирной линией обозначен суммарный коэффициент диффузии D₀ = D + K.

 

Результаты расчетов изображены на рис. 4 и 5. Главное отличие состоит в высотных профилях распределений водородсодержащих компонент.

 

Рис. 4. Полученный по результатам расчетов по первому методу высотный среднесуточный профиль концентраций химических компонент атмосферы Марса для средних широт и условия равноденствия.

 

Рис. 5. Полученный по результатам расчетов по второму методу высотный среднесуточный профиль концентраций химических компонент атмосферы Марса для средних широт и условия равноденствия.

 

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

Главными недостатками данных подходов является пренебрежение горизонтальным переносом компонент, особенно в турбопаузе, однако данные поправки будут учтены в будущих работах.

ЗАКЛЮЧЕНИЕ

Проанализированы различные модели описания турбулентной диффузии газовых компонент атмосферы Марса. Среди всех вариантов были выбраны две теории, которые в дальнейшем использовались при расчетах.

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

Осуществлено решение полученной системы уравнений с использованием модернизированного метода Ньютона. Модели показали существенные различия в распределении малых составляющих атмосферы, в частности, водородсодержащих молекул, что указывает на важность выбора модели для описания турбулентной диффузии при построении одномерной фотохимической модели атмосферы.

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

Авторы признательны Д.А. Климачкову за полезные обсуждения и помощь в работе.

Работа поддержана темой “Планета” госзадания Минобрнауки РФ.

About the authors

Е. О. Киливник

Институт космических исследований РАН; Московский физико-технический институт (национальный исследовательский университет)

Author for correspondence.
Email: Kilivnik.e@phystech.edu
Russian Federation, Москва; Москва

А. С. Петросян

Институт космических исследований РАН; Московский физико-технический институт (национальный исследовательский университет)

Email: Kilivnik.e@phystech.edu
Russian Federation, Москва; Москва

А. А. Федорова

Институт космических исследований РАН

Email: Kilivnik.e@phystech.edu
Russian Federation, Москва

О. И. Кораблев

Институт космических исследований РАН; Московский физико-технический институт (национальный исследовательский университет)

Email: Kilivnik.e@phystech.edu
Russian Federation, Москва; Москва

References

Supplementary files