CC BY
24
Применение двумерной модели для описания турбулентного переноса CO2 в пространственно-неоднородном растительном покрове
Ю. В. Мухартова1,а, Н.Т. Левашова1, А. В. Ольчев2,6, Н.Е. Шапкина1
1 Московский государственный университет имени М.В. Ломоносова, физический факультет, кафедра математики. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.
2Институт проблем экологии и эволюции имени А.Н. Северцова РАН.
Россия, 119071, Москва, Ленинский просп., д. 33.
E-mail: a [email protected], 6 [email protected] Статья поступила 29.07.2014, подписана в печать 30.09.2014.
Построена численная двумерная модель турбулентного переноса CO2 между пространственно-неоднородным растительным покровом и приземным слоем атмосферы. Модель основана на полуторном замыкании усредненных уравнений гидродинамики и диффузии СО2 в движущемся потоке. С помощью этой модели была выполнена оценка влияния неоднородности лесной растительности на турбулентный режим и потоки СО2 в приземном слое воздуха. Смоделирован перенос СО2 на опушке леса, а также в пределах сплошной вырубки.
Ключевые слова: двумерная модель, приземный слой атмосферы, турбулентный перенос CO2, полуторное замыкание, неоднородный растительный покров.
УДК: 519.63. PACS: 02.30.Jr, 92.60.Fm, 91.62.La.
Введение
Задача определения потоков водяного пара, углекислого газа и других парниковых газов между неоднородным растительным покровом и атмосферой является в настоящее время объектом многочисленных экспериментальных и модельных исследований. Значительный интерес к данной проблеме обусловлен главным образом существенной ролью растительности в формировании баланса парниковых газов в атмосфере. Растительный покров поглощает СО2 в процессе фотосинтеза и удерживает его в связанном состоянии на протяжении значительного промежутка времени, компенсируя таким образом его увеличение в атмосфере, вызванное антропогенной активностью. Согласно существующим оценкам [1], наблюдаемое современное глобальное потепление климата обусловлено именно устойчивым ростом содержания парниковых газов в атмосфере. Существующая современная система мониторинга за потоками парниковых газов между земной поверхностью и атмосферой главным образом сосредоточена на изучении переноса парниковых газов в пространственно-однородном растительном покрове. Исследования процессов переноса газовых компонент между пространственно-неоднородной растительностью и атмосферой проводятся довольно эпизодически, что прежде всего связано с недостатком методических подходов, позволяющих провести адекватное описание процессов переноса на границе раздела растительных сообществ. Эффективным инструментом для решения подобных задач могут быть процессориентированные двух- и трехмерные модели турбулентного переноса тепла, водяного пара и углекислого газа.
В настоящей работе для описания процесса переноса углекислого газа между неоднородным растительным покровом и приземным слоем воздуха была разработана двумерная модель, основанная на предложенной в [4-6] модификации полуторного замыкания усредненных уравнений гидродинамики.
С помощью разработанной модели были проведены расчеты деформации поля скорости ветра и коэффициентов турбулентности на опушке леса в случае наветренного края, а также внутри и над лесной сплошной вырубкой. С использованием модели переноса солнечной радиации в растительном покрове, а также модели фотосинтеза и дыхания растительности и почвы было найдено пространственное распределение концентрации СО 2 в рассматриваемых областях и рассчитаны вертикальные потоки СО2 на разных высотах над земной поверхностью.
1. Описание используемой модели 1.1. Основные искомые величины и система уравнений для них
Предположим, что на рассматриваемой территории рельеф местности и структура растительности не меняются вдоль некоторого направления. Выберем декар-тову систему координат таким образом, чтобы ось г была направлена вертикально вверх, а ось х — перпендикулярно направлению, вдоль которого растительность и рельеф можно принять однородными. Тогда задачу можно рассматривать как двумерную в области х е [-Ь, Ь], г е [0, Н]. Далее будем рассматривать именно такую модель.
Пусть V = {и, м] — скорость ветра, где и — ее горизонтальная составляющая, а да — вертикальная. С учетом характерного масштаба турбулентности в приземном слое атмосферы уравнение неразрывности в моделях турбулентного движения воздуха обычно [2, 3] записывают в виде
du dw div v = — + — = 0. dx dz
(1)
Рассмотрим модель турбулентного движения без учета тепловой конвекции. При этом компоненты скорости должны удовлетворять уравнению Навье-Стокса:
dv , VSp
— + (v, V)v =--- + vAv + f,
dt p0
где 5р — отклонение давления от гидростатического; Ро — плотность воздуха, которая в рассматриваемом приближении считается постоянной; V — коэффициент вязкости; / — массовая плотность внешних сил, обусловленная сопротивлением растительности.
Концентрация с углекислого газа описывается уравнением диффузии:
1 д
дс
— + (v, V)c = YcAc + fc,
(3)
где 7с — коэффициент диффузии, \с описывает мощность источников и стоков С02.
1.2. Усреднение системы уравнений
Предполагается, что процесс турбулентного переноса в приземном слое осуществляется с помощью вихрей. При этом типичный размер I (длина пути перемешивания) самых крупных (энергонесущих) вихрей в приземном атмосферном слое соизмерим с глубиной этого слоя, а их скорость V имеет порядок скорости потока. Поэтому турбулентное число Рейнольдса Л = Vll/v, характеризующее отношение динамических сил к вязким, много больше единицы [2, 3].
В турбулентном потоке происходит диссипация кинетической энергии и ее переход во внутреннюю энергию со средней скоростью е. Энергия крупномасштабных вихрей передается вихрям меньших размеров и так вплоть до минимальных диссипативных вихрей, в которых турбулентная кинетическая энергия переходит во внутреннюю за счет вязкого трения. Для того чтобы численное решение задачи было корректным, расчетная сетка должна иметь шаг, позволяющий учесть турбулентные вихри всех масштабов. Отношение I к типичному размеру п диссипативных вихрей в соответствии с предположением Колмогорова о зависимости п только от величин е и V имеет порядок Л^4 [2]. То есть даже в случае двумерной задачи для расчетов необходимо порядка (1/п)2 ~ Л3/2 ~ 1012 узлов сетки.
Один из возможных подходов к решению проблемы состоит в том, чтобы ограничиться исследованием нестационарного турбулентного течения только в масштабах, превышающих некоторый заданный размер. Масштабы вихрей, для которых прямое решение невозможно, моделируются как подсеточная турбулентность с использованием вихревой вязкости. Структура крупных вихрей получается в результате решения усредненных гидродинамических уравнений [7].
При усреднении системы уравнений используем разложение Рейнольдса: любая искомая функция ф представляется в виде суммы ф = Ф + ф' ее флуктуирующей части ф' и среднего значения Ф = (ф) по пространству и времени.
Итак, используя разложение Рейнольдса, представим все искомые функции в следующем виде: и = и + и', т = V + т', 5р = 5Р + 5р', с = С + с'. Усредняя уравнения (1) и (2) по пространству и времени, получим систему для усредненных компонент и и V скорости V = {и, т}:
ди = о
дх дг '
ди + и ди + V ди =
д1 дх дг
= -- —SP + vAU - — /и'2) - — (w'u') + Fu р0 дх дх \ / dz 4 '
—W дW W дW = д1 дх дz
1 д д _ д /-\
=--—SP + vAW - д (u'wЛ - д U'2) + Fw,
р0 дz дх х ' az \ /
д /—
w ' •
(4)
где Г = {Ри, РV} — средняя массовая плотность силы сопротивления растительности.
Усредняя уравнение (3), получаем уравнение для средней концентрации С02
дС Т7дС тт„ дС . „ д '—д '——, Л + и К + W^ = 7сДС " ^ <и'с > " ^ тс ') + РС,
дz
дх
дz
_ (5)
где FC = (fc) — средняя мощность источников и стоков вещества.
1.3. Аэродинамическое сопротивление растительности
При моделировании процесса турбулентного движения воздушного потока удобно рассматривать растительный покров как сплошную среду с непрерывно распределенными объемными силами сопротивления. Величина силы сопротивления, действующей в некотором объеме, пропорциональна плотности р0 движущегося воздуха, квадрату его средней скорости и элементу площади обтекаемого препятствия. Для массовой плотности F силы вязкого трения в пределах кроны дерева будем использовать выражение F = -cd ■ LAD - \V\ ■ V, где V = {U, W} [6, 8, 9]. В этом выражении cd — коэффициент аэродинамического сопротивления. В модельных экспериментах часто используют некоторое среднее значение этого параметра для всего слоя растительности. Измерения коэффициента cd, проведенные в лабораторных условиях, показали его довольно слабую зависимость от скорости ветра [8]. В нашем исследовании мы приняли cd = 0.2. Величина LAD — это так называемая плотность листовой поверхности, т. е. суммарная площадь поверхности листьев в единице объема. Эта величина зависит от структуры и видового состава растительности, она изменяется с высотой и может быть определена экспериментально. В настоящей работе для параметризации профиля LAD(z) использовано эмпирическое выражение [10, 11], предполагающее максимум фитомассы в верхней части кроны деревьев:
LAD(z) =
LAI ■
2h
- Ы nh-*)) - 2 ^^))
wn V \ «crown / 2 \ «crown JJ
h r
LAI =
LAD(z) dz,
где Ясг0№п — высота кроны, гс — высота крепления кроны. Площадь поперечного сечения стволов деревьев в нашем исследовании была принята равной 0.03 м2.
1.4. Потоки Ш 2 в растительности
Для описания источников и стоков СО2 в растительном покрове в модели был использован комплексный подход, основанный на эмпирической модели
Болла с соавторами [12] в модификации Леннинга [13, 14], уравнении Беера-Ламберта для описания процесса ослабления солнечной радиации в растительном покрове [20], а также алгоритме для описания зависимости устьичной проводимости листьев растений для Н2О и СО2 от приходящей фотосинтетически активной радиации (ФАР) [15, 16]. Эмпирические параметры, необходимые для проведения расчетов скорости фотосинтеза, дыхания и устьичной проводимости деревьев, были подобраны из литературных источников [15].
Функция FC, характеризующая интенсивность источников/стоков CO2 в пределах растительности, имеет вид FC = — LAD -An, где An — скорость потока CO2 между листом и прилегающим к листу воздухом (моль CO2 м~2сВеличина An связана с устьичной проводимостью листьев gs (моль CO2 м~2си концентрацией CO2 у поверхности листьев cs (моль CO2/моль сухого воздуха) соотношением
An = — (gs — go)(cs — г, )(1 + Ds /Do),
где g0 — значение gs в световом компенсационном пункте, Г, — углекислотный компенсационный пункт, Ds — дефицит упругости водяного пара, a1 и D0 — эмпирические коэффициенты [12-14].
Зависимость устьичной проводимости gs от интенсивности приходящей фотосинтетическиактивной радиации PAR имеет вид gs = gmax ■ f (PAR), где f (PAR) — функция отклика gs на приходящую к поверхности листа PAR. В соответствии с [15, 16] для описания зависимости gs от приходящей солнечной радиации внутри растительного покрова было использовано выражение f (PAR) = 1 — e~&rpAR, где ßsi — эмпирический параметр, определяющий наклон световой кривой при PAR ^ 0.
В соответствии с [20] определим величину PAR на высоте z внутри растительности выражением
PAR(z) = PARh ■ exp < —k ■
LAD(z') dz'
Здесь PARН — интенсивность потока фотосинтетиче-скиактивной радиации над кронами деревьев высоты Н, £ — коэффициент затухания солнечной радиации в растительном покрове (примем его равным 0.5, что соответствует ослаблению при сферическом распределении фитоэлементов деревьев по углам наклона).
Эмиссия СО2 с поверхности почвы, обусловленная автотрофным и гетеротрофным дыханием почвенной биоты, в рамках модельных экспериментов была принята постоянной (4.0 мкмоль• м-2-с-1).
1.5. Замыкание системы уравнений
Система уравнений (4) и уравнение (5) содержат неизвестные моменты: (и'2), (и'т')
w
u'c'
{т'с'). Для этих величин могут быть получены свои уравнения, но они будут содержать неизвестные моменты более высокого порядка, и так далее до бесконечности. Это так называемая проблема замыкания усредненных уравнений — уравнения для моментов п -го порядка всегда содержат слагаемые с моментами
(п+1) -го порядка [2, 3]. Замкнуть систему усредненных уравнений без дополнительных предположений принципиально невозможно. Поэтому неизвестные величины обычно выражают через известные, исходя из каких-либо физических соображений и гипотез.
В настоящей работе мы будем использовать так называемое полуторное замыкание. Указанные выше моменты представляют собой турбулентные потоки, за счет которых происходит перемешивание, схожее с диффузионным, но по величине превосходящее его на несколько порядков. Это перемешивание называют турбулентной диффузией. При полуторном замыкании турбулентные потоки выражаются по аналогии с обычной молекулярной диффузией в приближении Бусси-неска [4, 5]:
С
u w
= /dW dUN V dx dz
(u2)
dx
= E — 2K^~, (w'2) = E — 2K
dW dz
где K — коэффициент турбулентной диффузии, E = 0.5-
+ ^т— турбулентная кинетическая
энергия (ТКЭ). Аналогично могут быть выражены турбулентные потоки для С02:
(
u c
= - кдС
— 1vc <-, >
дх
w
= - K — dz
где Кс — коэффициент турбулентной диффузии для переносимого вещества.
Коэффициенты К и Кс полностью определяются свойствами потока, а не вещества, как в случае молекулярной диффузии, и зависят от скорости ветра. В теории турбулентности, основанной на понятии коэффициента турбулентной диффузии, считается, что отношение коэффициентов К и Кс есть постоянная величина Бе = К/Кс, называемая турбулентным числом Шмитда. «Классическое» микрометеорологическое предположение состоит в том, что Бе = 1 [3]. Однако многие эксперименты, как натурные, так и проведенные в аэродинамической трубе, говорят о том, что Бе < 1. В частности, в [17] приводится значение Бе = 0.75, а в работе [18] в результате измерений было получено значение Бе = 0.6. В нашем исследовании мы будем использовать Бе= 0.75.
Поскольку при отличных от нуля скоростях ветра турбулентные потоки существенно превышают потоки, вызванные молекулярной диффузией, в нашем исследовании вклад молекулярной диффузии в уравнения (4), (5) не рассматривается. С использованием сделанных допущений приходим к системе уравнений (6) для средних значений компонент скорости и среднего отклонения давления от гидростатического:
ddU+и du+W ddU =
dt дх dz
= _ _L dLSP-— —
p0 дх дх дх
д / ди) д , ^ + ^ K— + K z z z
дW u—W w—W =
()t дх дz
+
дw
дх
+ Fu
__dE d (KdW
р0 dz dz dx y dx
+
9 ío,sdw\ 9 („ди\ ^
+ 2K— + — K— + Fw •
dz \ dz J dx \ dz
dü dW _ о dx dz
(6)
и следующему уравнению диффузии для средней концентрации CO2:
/ 6С\ d / dC\ п (K dx) + d-z{K dz) + Fc
de T1dc de d / de\ d . „ _ . ä + ü dTx + W dz _ dCx\ + dz) + r C.
(7)
1.6. Вычисление коэффициента турбулентной диффузии
Коэффициент турбулентной диффузии K может быть выражен через кинетическую энергию турбулентного движения E и скорость диссипации турбулентной энергии е по формуле K = C^E2е-1, где C^ — безразмерный коэффициент пропорциональности [3]. Измерения, проведенные в аэродинамической трубе, дают значение C^ = 0.09 [4].
Одним из способов нахождения значения турбулентной кинетической энергии и скорости ее диссипации е является решение для них системы дифференциальных уравнений типа уравнений диффузии в движущемся потоке. В общем случае второе уравнение этой системы формулируется для некоторой вспомогательной функции ф, характеризующей масштаб турбулентности. В качестве функции ф в различных моделях берутся следующие три величины: е, е/E или l. Исследование и сравнительный анализ данных моделей проведены в работах [4, 5]. В частности, показано, что при отсутствии растительности указанные три модели в целом дают схожие результаты. Однако для ф = е и ф = е/E модели позволяют лучше воспроизводить турбулентное поле в случае набегания потока на препятствие, например, границу леса, чем при ф = l . Кроме того, выбор ф = е/E оказывается предпочтительным, так как он позволяет избежать неопределенности с числами Прандтля и Шмидта, характерной для ф = е в пределах растительности. Поэтому в нашей работе мы рассматриваем ф = е/E.
Система уравнений для функций E и ф имеет вид
дА + и дА + w дА =
dt dx dz
д ( K dE\ д ( K dE\ „ = д"1 + д- фдZ) + PE - е,
dx
a'E dx )
dz
^ + ü ^ + W ^ _ dt dx dz
_ d í K d ( K d^> \
dx dx ) dz y av dz )
+
+ I (C^i • Pe - Cv2 • e) - Дф.
Здесь афЕ и аф — число Прандтля для ТКЭ и турбулентное число Шмидта для ф соответственно. Они позволяют связать коэффициенты турбулентной диффузии для функций Е и ф с функцией К: аф = К/КЕ, аф = К/Кф. Величина РЕ представляет собой скорость образования
турбулентной кинетической энергии, и в рассматриваемом двумерном случае она имеет вид
»-"( (£)' + (fí) + k(
. . dü dw\
^z) dz + Ж)
Множители Сф1 и Сф2 являются модельными константами. Слагаемое Дф характеризует увеличение диссипации турбулентной кинетической энергии за счет взаимодействия с растительностью. В работах [4, 5] в предположении, что на расстояниях порядка средней длины листьев увеличение турбулентной кинетической энергии за счет взаимодействия ветра с растительностью почти полностью сбалансировано ее диссипацией за счет перехода в тепло, для этого слагаемого получено следующее выражение:
Дф = ф • (Сф1 - Сф2) • ßd ■ Cd ■ LAD-\V E,
где множитель ßd характеризует величину потерь энергии за счет взаимодействия с препятствиями. На основе экспериментальных данных в [4] для ßd в первом приближении получена оценка ßd « 12СУ2.
1.7. Начальные и граничные условия
В качестве начального условия для горизонтальной компоненты скорости ветра возьмем хорошо известный логарифмический профиль, описывающий распределение скорости в приземном слое воздуха по высоте над горизонтально однородной растительностью [3]:
U|t_o _ -
Ч—) •
\ zvisc /
где иФ = у —и'т' — динамическая скорость; к — постоянная фон Кармана, которую обычно принимают равной 0.4; гу;5с — параметр шероховатости; й — высота слоя вытеснения, определяемая шероховатостью поверхности и свойствами растительности. При расчетах в начальный момент времени вертикальную компоненту скорости ветра и избыточное давление будем считать равными нулю.
В качестве начального распределения по высоте коэффициента турбулентной диффузии К возьмем эмпирическое соотношение, полученное в [19] в пространственно-однородном случае:
3/2
K|t_o _
ки* z
1 + 5z/lmo
i - z
где Ег — высота атмосферного пограничного слоя; ¿мо — масштаб длины Монина-Обухова, значение которой определяет высоту, до которой влияние сил трения на турбулентное движение воздушного потока превышает влияние сил плавучести [3].
Начальные значения для турбулентной кинетической энергии Е и функции ф получим из соотношений, связывающих эти величины с коэффициентом турбулентной диффузии К и масштабом энергонесущих вихрей I [2, 3]:
2
E|t_o _
1
(K)'
I С,Е
t_0
t_0
Известно [2, 3], что I при малых г ведет себя как кг, поэтому примем /|^=0 = кг.
Начальное значение концентрации С углекислого газа положим равным некоторому постоянному фоновому значению С0 = 380 ррт (мкмоль СО2 на моль сухого воздуха).
Предположим, что все неоднородности растительного покрова расположены строго внутри расчетной области, достаточно далеко от ее границ. Тогда можно считать, что вблизи вертикальных границ х = ±Е все искомые функции не зависят от х, и сформулировать граничные условия в следующем виде:
ди
х
дф х
= 0,
х=± Е
х=±Е
дЖ х
= о, #-^
х
= 0, тт-
х=± Е
= 0, тт-
х=±Е
ди
дг
г=И
= 0, ^
г
г=И
дЕ
=0, дЕ
г
г=И
дЕ
дх х= ±Е
дС
х х =±Е
и при г=
= 0, дф
дг
= 0,
= 0.
= 0.
г=И
= 0,
, суЕ
ф|г=0 =
г=0
К
г=0
дг
= 0.
г=0
1 ди + иди + жди =
2 дг дх дг
д
= - (О + ¿(к-> + -'К
д_ дх
- ой • ЬЛЭ •¡У!• и
дг )
д_ дг
К!)
1 дЖ идЖ =
2 дг дх дг
д (КдЖ\ д
дх у дх ) дг - ой • ЬЛЭ •¡У]• Ж,
«) + Нк
,ои\ _
дг)
(8)
если г £ [г/, ¿/+1/2];
1 и = д_( Р Е
2 дг дх у р0
,
1 ж = е_ {р Е
2 дг дг у р0
На верхней границе области влияние неоднородно-стей у поверхности сказывается мало, поэтому можно считать избыточное давление 5Р равным нулю, а концентрацию С переносимого вещества равной фоновому значению С0.
При г = 0 компоненты средней скорости будем считать равными нулю. Для остальных функций граничные условия при г = 0 сформулируем следующим образом: дЕ дг
Предполагая, что на нижней границе расчетной области существует постоянный поток ц углекислого газа, обусловленный автотрофным и гетеротрофным дыханием почвенной биоты, сформулируем для С граничное условие как
К дС
дг г=0
2. Разностная схема
Полученную начально-краевую задачу будем решать численно при помощи метода конечных разностей. Введем в расчетной области сетку хп, гт, г/. Пусть на / -м временном слое компоненты скорости, избыточное давление, кинетическая энергия турбулентного движения и функция ф известны. Для того чтобы найти давление и компоненты скорости на новом (/+1) -м слое, используем разностную схему с расщеплением [7]. Введем вспомогательный слой / + 1/2 и осуществим переход со слоя / на слой / + 1 в два шага:
(9)
если г £ [/+1/2, /+1]. При этом будем считать, что на целых слоях выполняется уравнение неразрывности.
Переход / ^ / + 1/2 осуществим с помощью локально-одномерной схемы через промежуточный слой / + 1/4. При этом будем учитывать знак множителей и и Ж в слагаемых, отвечающих за перенос. Рассмотрим переход со слоя / на слой / + 1/4 для функции и : если и'пт ^ 0, то
и/+1/4 _ и/ и> + Х/А - и/'+1/4
ип,т ип,т . т т/ п,т п—1,т _
г ип,т
Т
= 21х
и/+1/4
хп хп — 1
+ Егх [Ж/] - 0Л ЬЛЭп
У/
г п,т
и,
-/+1/4
где Уп,т = {ип'.т,, ¿хх[и] и Егх [Ж] -разностные аналоги дифференциальных операторов д/дх (К • ди/дх) и д/дг (К • дЖ/дх) соответственно,
и/+1/4 и/ и/'+1/4 _ +
ип,т - ип,т + и/ п+1,т ип,т _ + и п,1
Т
= 21х
и/+1/4
х +1 - х + Егх [Ж/] - еа ЬЛЭп
У/
г п,т
и
/+1/4
если ип,т < 0. Для функции Ж этот переход осуществляется аналогично.
Переход со слоя / + 1/4 на слой / + 1/2:
и/+1/2 _ Т //+1/4
и/+1/2 _ г т/+1/2
п,т + ш/+1/2ип,т и п,т— 1
гт гт—1
Жп,т / ^ 0;
= Ег
и/+1/2
и/+1/2
и/+1/2 т 7/ +1/4 и/+1/2 _ и/+1/2
ип,т - ип,т + ш/+1/2ип,т+1 ип,т = .
+ Жп,т — Ег
Т гт+1 гт
< 0,
где Егг [и/+1/2] — разностный аналог оператора д/дг (К • ди/ дг).
Заметим, что рассматриваемые уравнения (8) по структуре схожи с уравнениями для функций Е, ф и С. Уравнения для этих функций будем решать по сходной схеме с той лишь разницей, что в них осуществляется переход со слоя / на слой / + 1 через вспомогательный слой / + 1/2.
Рассмотрим теперь переход / + 1/2 ^ / + 1. Из системы (9) для избыточного давления 5Р получаем уравнение
\(6Р Р\ = 1 д (ди дЖ\ А VР0 + V = -2ш\ох + ж) •
Аппроксимируем производную по времени в правой части этого уравнения конечной разностью:
4 / + Е) =
р0
ди/+1 дЖ/+1
т |\ дх дг
-
(ди/+-1/2 + дЖ/+1/2\
х
+
}
Так как на целых слоях должно выполняться уравнение неразрывности, то выражение (6Р>+1/ро + Е>) должно быть решением уравнения Пуассона
Д
SPi+l Ро
■Е' I =7
1 faui+1'2 dwi+!/2
дх
дг
(10)
Аппроксимируем уравнение Пуассона (10) на сетке и дополним соответствующую алгебраическую систему уравнениями, получаемыми в результате аппроксимации граничных условий для 6Р1+1. Полученную систему будем решать методом матричной прогонки.
Когда величины Р}^т = (8Р'+1/ро + Е>) найдены, компоненты скорости на слое / + 1 могут быть вычислены по формулам
ц}+1 = [//+1/2 _
' 'а, т ' 'а, т 1
т т Т
р/+1 _ р/+1
п+\,т п—\,т
Хп+1 Хп-1 pj+1 _pj+1
п,т+1 п,т— 1
^ т I ^т I
3. Результаты моделирования турбулентного режима и потоков СО 2 в неоднородном растительном покрове
На основе предложенного алгоритма был проведен ряд расчетов поля скорости ветра и концентрации углекислого газа для разных видов неоднородности растительного покрова: наветренной опушки леса и сплошной вырубки в лесу. Средняя высота леса была взята равной 20 м, листовой индекс LAI = 4, высота крепления кроны гс = 8 м, плот-
ность потока приходящеи на верхнюю границу растительного покрова фотосинтетическиактивной радиации PAR/, = 2000 (мкмоль-м~2-с-1), динамическая скорость и* = 0.4 м/с.
На рис. 1 представлены результаты расчетов, проведенных для случая набегания воздушного потока на опушку леса. В рассмотренном случае ветер дует слева направо, лес расположен в области его грани-
ца обозначена светло-серой линией. Для каждого из рисунков, демонстрирующих распределение компонент скорости ветра, турбулентной кинетической энергии, коэффициента турбулентности и концентрации СО 9, справа от каждого из рисунков приведена шкала, устанавливающая связь между цветом и значением соответствующей величины. На графике, демонстрирующем плотность вертикального турбулентного потока СО 9, сплошная линия соответствует величине (ш'с') на высоте 82 м, пунктирная — на высоте 25 м, точками обозначена величина (ш'с') на высоте 42 м над поверхностью.
На рис. 2 приведены результаты расчетов в случае наличия в лесу сплошной вырубки шириной 200 м, расположенной в области х е [—100, 100]. Как и в предыдущем случае, ветер дует слева направо, граница леса также обозначена светло-серой линией. Смысловое значение цвета на рисунках пояснено на приведенных справа от них шкалах.
Результаты модельных экспериментов показывают, что неоднородная растительность оказывает существенное влияние не только на распределение поля ветра и коэффициентов турбулентности, но также и на
Горизонтальная скорость U
Вертикальная скорость
Турбулентная кинетическая энергия Е
-Г, м
Коэффициент турбулентности К
X, м
Концентрация С02
-Г, м
I
Плотность вертикального турбулентного потока С02
395
1380
-100 0 -Г, М
100 200 300
о m -о
с; о
О
О -о
о
U
40 2-О--2-4 -6 -8-10 -12
\
\ \ \
\ ч \ ч \
ч s \ Ч
-z-25 м z=42 м — z=82 м
0
-Г, м
Рис. 1. Модельный эксперимент, описывающий поле ветра и турбулентности в приземном слое воздуха при
набегании воздушного потока на опушку леса
Горизонтальная скорость U
Вертикальная скорость
Турбулентная кинетическая энергия Е
Коэффициент турбулентности К
Концентрация СО
Плотность вертикального турбулентного потока СОг
Рис. 2. Модельный эксперимент, описывающий поле ветра и турбулентности в пределах вырубки в лесу.
Ширина вырубки принята равной 200 м
вертикальные потоки. Пренебрежение неоднородностью растительного покрова при изучении процессов вертикального переноса может привести к существенным погрешностям в оценках локальных и интегральных потоков.
Заключение
Разработана двумерная модель турбулентного переноса СО9 в неоднородной растительности. Реализован алгоритм численного решения соответствующей начально-краевой задачи, основанный на локально-одномерной разностной схеме. С помощью предложенной модели осуществлен расчет полей скорости ветра, коэффициентов турбулентности, турбулентной кинетической энергии, распределения концентрации и вертикальных потоков СО 9 для двух типов неоднородности растительности: опушки леса и вырубки в лесу. Результаты расчетов показывают существенную деформацию поля ветра и турбулентного режима вблизи препятствий, образованных различными формами растительности. Разработанная модель может быть эффективным инструментом для оценки влияния пространственно-неоднородных растительных сообществ на интегральные потоки тепла, Н9О и СО9 между земной поверхностью и атмосферой. Данные оценки необходимы как для прогнозных климатических, так и для экологических моделей, направленных на планирование и оптимальное управление лесными ресурсами.
Работа выполнена при финансовой поддержке РНФ (грант 14-14-00956).
Список литературы
1. IPCC. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the IPCC. Cambridge, 2013.
2. Wyngaard J.C. Turbulence in the atmosphere. Cambridge, 2010.
3. Garratt J. R. The atmospheric boundary layer. Cambridge, 1992.
4. Sogachev A., Panferov O. // Boundary-Layer Meteorol. 2006. 121, N 2. P. 229.
5. Sogachev A. // Boundary-Layer Meteorol. 2009. 130, N 3. P. 423.
6. Sogachev A., Menzhulin G.V., Heimann M., Lloyd J. // Tellus. 2002. 54B, N 5. P. 784.
7. Белоцерковский O.M. Численное моделирование в механике сплошных сред. М., 1994.
8. Дубов А.С., Быкова Л.П., Марунич С.В. Турбулентность в растительном покрове. JL, 1978.
9. Бояршинов M.F., Горемыкин В.Д. // Матем. модел. 2004. 16, № 7. С. 31.
10. Олъчев А.В., Радлер К. // Изв. Самарского научн. центра РАН. 2009. 11, № 1(7). С. 1538.
11. Olchev A., Radler К., Sogachev A. et al. // Ecological Modelling. 2009. 220. P. 3046.
12. Ball J. Т., Woodrow I.E., Berry J.A. // Progress in Photosynthesis Research / Ed. by I. Biggins. 1987. IV. P. 221.
13. Leuning R. // Australian J. of Plant Physiology. 1990. 17. P. 159.
14. Leuning R. // Plant, Cell and Environment. 1995. 18. P. 339.
15. Oltchev A., Constantin J., Gravenhorst G., Ibrom A. // J. of Hydrology and Hydromechanics. 1997. 45, N 1-2. P. 5.
16. Oltchev A., Ibrom A., Constantin J. et al. // J. Phys. Chem. Earth. 1998. 23, N 4. P. 453.
17. Stull R.B. An introduction to boundary layer meteorology. 19. Brost R.A., Wyngaard J.C. // J. Atmos. Sci. 1978. 35. Dordrecht, 1988. P. 1427.
18. Flesch T.K., Prueger J.H., Hatfield J.L. //Agricult. and Fo- 20. Monsi M, Saeki T. // Japanese J. of Botany. 1953. 14. rest Meteorol. 2002. 111. P. 299. P. 22.
Application of a 2D model for describing the turbulent transfer of CO2 in a spatially heterogeneous vegetation cover
Yu.V. Mukhartoval a, N. T. Levashova1, A.V. Oltchev2 b, N. E. Shapkina1
1 Department of Mathematics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.
2A.N. Severtsov Institute of Ecology and Evolution, Russian Academy of Sciences, Moscow 119071, Russia.
E-mail: a [email protected], b [email protected].
A numerical 2D model of the turbulent transfer of CO2 in the air between a spatially heterogeneous vegetation cover and the atmospheric surface layer was developed. The model is based on the one-and-a-half closure scheme for averaged equations of hydrodynamics and diffusion of CO2 in the air. The model was applied to assess the possible influence of forest canopy heterogeneity on the turbulent regime and CO2 fluxes in the atmospheric surface layer. CO2 exchange was simulated both at the windward forest edge and around a forest clearing.
Keywords: 2D model, atmospheric surface layer, turbulent transfer of CO2, one-and-a-half closure, non-uniform vegetation.
PACS: 02.30.Jr, 92.60.Fm, 91.62.La.
Received 29 July 2014.
English version: Moscow University Physics Bulletin 1(2015). Сведения об авторах
1. Мухартова Юлия Вячеславовна — канд. физ.-мат. наук, ст. науч. сотрудник; тел.: (495) 939-10-33, e-mail: [email protected].
2. Левашова Наталья Тимуровна — канд. физ.-мат. наук, доцент; тел.: (495) 939-10-33, e-mail: [email protected].
3. Ольчев Александр Валентинович — канд. геогр. наук, ст. науч. сотрудник; тел.: (926) 246-13-42, e-mail: [email protected].
4. Шапкина Наталья Евгеньевна — канд. физ.-мат. наук, доцент; тел.: (495) 939-10-33, e-mail: [email protected].
