CC BY
58
Вычислительные технологии
Том 9, № 1, 2004
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕПЛОВЛАГООБМЕНА В СИСТЕМЕ АТМОСФЕРА - ПОЧВА В ЗАСУШЛИВЫЙ
ПЕРИОД*
В. А. Шлычков Институт водных и экологических проблем СО РАН,
Новосибирск, Россия e-mail: slav@ad-sbras.nsc.ru
A mathematical model for heat and moisture transfer between the atmosphere and the soil active layer in summer period is presented. Data of observations help to identify the governing parameters of the problem. Since the problem was integrated over the month period, the model can be interpreted as a climatic model. Sensitivity of the climatic parameters to variations of some physical parameters was investigated.
Введение
Основными физическими механизмами, определяющими характер термодинамических процессов вблизи поверхности земли, являются турбулентный тепловлагоперенос в приземном слое атмосферы и массоэнергообмен в деятельном слое почвы. Теоретические исследования показали высокую восприимчивость атмосферы к испарению с земной поверхности. При этом наиболее интересными и физически содержательными оказываются уравнения баланса тепла и влаги на подстилающей поверхности, описывающие тонкие эффекты трансформации радиационных и турбулентных потоков.
Современные математические модели испарения с поверхности почвы [1,2] ориентированы на воспроизведение гидрологического цикла крупных речных бассейнов. В практических расчетах водосборных территорий исследователи, как правило, прибегают к ряду упрощений, схематизируя сложную физику испарения с помощью тех или иных параметризаций. Чувствительность расчетных характеристик к параметрическому замыканию зависит не только от качества математической модели, но и от типа физических процессов, которые она описывает. Весьма трудными для описания и слабо изученными являются процессы естественного испарения с почвы в засушливые годы, когда сумма осадков значительно ниже нормы. Расчет засушливых сезонов особенно актуален для условий Сибири, где вероятность развития экстремальных погодных явлений весьма высока.
* Работа выполнена в рамках проекта 13.10 фундаментальных исследований Президиума РАН и при поддержке Российского фонда фундаментальных исследований, грант № 03-05-96825.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.
Цель данной работы заключается в построении математической модели взаимодействия приземного слоя атмосферы и почвы и ее апробации в задаче исследования чувствительности локальных климатических параметров в летний сезон. Подобная задача без учета суточных вариаций метеорологических величин была рассмотрена в [3] в аспекте изучения режимов просыхания почвы.
1. Постановка задачи
Математический аппарат, используемый при построении модели, основан на детерминированном описании процессов в атмосфере и почве [4]. Так, восстановление вертикальных профилей в атмосферном слое постоянных потоков (СПП), примыкающем к подстилающей поверхности, проводится на основе теории подобия Монина — Обухова с универсальными функциями Бусинжера [5]. Основные соотношения модели без учета растительности имеют вид
u(z) = ^fu(C), û(z) = во + -fe(Z), q(z) = qo + ^fe «). (1)
к к к
Здесь z — вертикальная координата, направленная вверх (начало координат совмещено с подстилающей поверхностью); и — скорость ветра; в — температура воздуха; q — удельная влажность воздуха; во, q0 — температура и влажность на уровне шероховатости ze; и*, в*, q* — масштабы скорости, температуры и влажности; к — постоянная Кармана; fi — универсальные функции; Z = z/L, L — масштаб турбулентности в СПП, для определения которого используется уравнение
Ri = Zfe/f2 при Z = h/L, (2)
где Ri = \(въ, — в0)/u2h — интегральный критерий устойчивости; Л — параметр плавучести; ин,вн — значения скорости ветра и температуры на верхней границе СПП z = h.
Уравнения взаимосвязанного переноса тепла и влаги в почве для неизотермических условий имеют вид
ÔTS = Я dt dz dz
дп д ^д'Ф дК д дТ3 е3 дЬ дг дг дг дг ш дг рш'
ду3 _ ^^ -
дЬ дг ь дг
Здесь Т3, п, Ф, ^ — искомые температура почвы, объемное влагосодержание, потенциал почвенной влаги и удельное содержание водяного пара в почвенном воздухе соответственно; е3, Х3 — объемная теплоемкость и коэффициент кондуктивной теплопроводности почвы; — скрытая теплота конденсации; — скорость фазовых переходов; К — гидравлическая проводимость; Хш — коэффициент термовлагопроводности; ^ — коэффициент диффузии пара; — плотность воды. В отличие от традиционной модели Ричардса [6] уравнение переноса жидкой влаги в (3) содержит слагаемое, описывающее неизотермический влагообмен. Механизм термопереноса может оказаться ведущим в условиях иссушения и высоких температурных градиентов в почве, связанных с радиационным нагреванием поверхности.
Перейдем к постановке краевых условий. На верхней границе СПП будем считать известными поля скорости, температуры и влажности
и = пн, в = вн, д = дн при г = к. (4)
Граничные условия на подстилающей поверхности вытекают из анализа балансовых соотношений задачи в полной постановке [4]. Уравнение баланса тепла имеет вид
Н + Е + Я = В при г = 0, (5)
где Н = раврсвин(вн — во) — турбулентный поток тепла; ЬшЕ = раЬшсвпн(дн — до) — поток скрытого тепла; Я — радиационный баланс; В = Л3дТ3/дх + Бду3/дх — поток тепла в почву; ра — плотность воздуха; ср — удельная теплоемкость воздуха; вв — интегральный коэффициент переноса тепла. Уравнение баланса влаги запишем в виде
^Т = -— Р^ при г = 0, (6)
аЬ рт
где Р^а) = ввин(дн — д0) + идт — атмосферный поток влаги; идт — интенсивность осадков; Р«з) = К • дф/дх + К + • дТ3/дх + Б • ду3/дх/рш — поток влаги в почву; — глубина "лужи" на поверхности (горизонтальный отток воды по склону в данной постановке не учитывается). При > 0 поверхностный слой полагаем полностью насыщенным
ф = при х = 0. (7)
Если интенсивность осадков не превышает впитывающую способность почвы, то уравнение (6) при кш = 0 служит для расчета поверхностного влагосодержания.
При определении внутрипочвенного испарения будем использовать гипотезу о термодинамическом равновесии между жидкой и парообразной фазами воды и положим
у3 = ра (П — п)д(Т3)ехр(Лф/Т3), (8)
1
где П — пористость; 7 — объемный вес почвы; д — насыщающая влажность; Л — константа;, ф — термодинамический потенциал почвенной влаги. Уравнение для у3 в (3) превращается в этом случае в соотношение для расчета скорости фазовых переходов.
На нижней границе деятельного слоя почвы Н3 будем считать известными температуру и влажность
Тз = Тн, п = Пн при х = Н3. (9)
Для описания перепадов температуры и влажности в вязко-буферном подслое, примыкающем к поверхности почвы и обусловленном молекулярными эффектами, будем использовать полуэмпирические параметризационные формулы вида [7]
во — Тз = ^вв*(и*хв/^а)0'45, до — дз = ^д(дн — до), (10)
где ^в, ^д — безразмерные параметры, подлежащие определению; иа — молекулярная вязкость воздуха. Способ задания ^в указан ниже, а величина ^д служит для удовлетворения уравнения баланса влаги в условиях иссушения поверхности. Как известно, объемная влажность почвы в естественных условиях не может опуститься ниже значения максимальной гигроскопичности пмг. Выполнение условия п > Пмг контролируется в ходе
интегрирования на каждом временном шаге и при необходимости проводится коррекция значений п путем принудительного задания п = Пмг- Это приводит к нарушению баланса потоков в (6). Однако, считая параметр свободным, удовлетворим уравнению (6) за счет специального выбора в виде
В сущности, введением вязко-буферного подслоя мы ослабляем градиент приземной влажности так, чтобы скорость испарения находилась в согласии с ограниченной скоростью поступления почвенной влаги к поверхности. В условиях свободного испарения при достаточной влажности поверхности будем полагать = 0.
В процессе испарения, когда отсутствует дефицит почвенной влаги и имеет место значительный скачок температуры в вязко-буферном подслое, может оказаться, что искомая влажность с0, равная по условию влажности почвы при г = 0, превышает свое предельное (насыщающее) значение с/, соответствующее температуре на уровне шероховатости. В этом случае будем полагать с0 = с/, а разность с3 — с/ отнесем к вязкому подслою в поле влажности. При этом излишнюю влагу, фиктивно испарившуюся в результате завышения предельного значения с, "сбросим" обратно в почву путем коррекции потока атмосферной влаги на величину упомянутой разности, сохранив тем самым водный баланс в системе.
Начальные условия для температуры почвы выбирались из данных наблюдений с последующей сплайн-интерполяцией по глубине, а начальный профиль влажности принимался постоянным по вертикали и равным полевой влагоемкости для рассматриваемого типа почвы.
2. Задание связей и идентификация параметров
Радиационный баланс в (5) рассчитывался по соотношению
где I — приходящая коротковолновая радиация; Г — эффективное длинноволновое излучение. При расчете I, Г использовались соответственно формулы Альбрехта и Брента, обычно применяемые в метеорологии [3].
Для замыкания уравнений (3) необходимо задать влажностные характеристики почвы в виде функциональных зависимостей п, К от ф. Соответствующие формулы взяты из монографии [6], почвенно-гидрологические константы подбирались для условий станции Огурцово (г. Новосибирск). Считалось, что уровень грунтовых вод расположен достаточно глубоко (более 8 м).
Задание коэффициента Лад, характеризующего влияние градиента температуры на миграцию влаги в почве, затруднено по причине отсутствия экспериментальных данных по изучению термопереноса влаги. Определение Лад проводилось на основе принципа взаимности Онсагера, который требует симметрии матрицы феноменологических коэффициентов. При некоторых упрощающих предположениях условие симметричности можно записать в виде уравнения
^я = 1 + (Сз — Со )св Пн/Р^-
Я = I — Г,
(11)
+ Р'Ш фК = 0,
связывающего искомую величину с известными параметрами.
Шаг интегрирования по времени составлял 30 мин, период физического моделирования — 30 сут., дискретизация системы (3) по вертикали осуществлялась с шагом 1 см. Принято П = 0.5, 7 = 1.2 г/см, ги = = 2 см, И = -80 см. Численные значения остальных параметров взяты равными или близкими к используемым в [3].
В рамках поставленной задачи разработана схема идентификации параметров и верификации модели. Для этого использовались данные наблюдений на теплобалансовой станции Огурцово, включающие шестисрочные градиентные измерения ветра, температуры и влажности на уровне 0.5 и 2 м над оголенной поверхностью почвы. Кроме того, в расчетах использованы информация о нижней и общей облачности, данные по температуре почвы на поверхности и четырех нижележащих горизонтах, а также данные об осадках. В целях калибровки механизма испарения моделирование проводилось для июля 1989 г., который характеризовался засушливой погодой (объем осадков составил 22 мм).
Рассмотрим задачу определения параметра ^, характеризующего вязкий температурный подслой. Заметим, что все величины, входящие в соотношение (10), известны из данных наблюдений и в результате решения задачи СПП. Это позволяет рассчитать оптимальное значение ^ с помощью статистического моделирования на основе требования минимизации отклонений расчетной температуры от фактической. Алгоритмическая база оптимизации параметра заимствована из [8]. В результате проведенных расчетов выяснилось, что значения ^ существенно зависят от стратификации приземного слоя и сильно различаются в дневной и ночной периоды суток. В связи с этим оценка этой величины проводилась раздельно для устойчивых (в* > 0) и неустойчивых (в* < 0) режимов. После применения статистической модели получены следующие значения: ^ = 0.18 для в* > 0 и = -0.02 для в* < 0.
Проведем оптимизацию параметров длинноволнового излучения Г, пользуясь фактической информацией по энергетическому бюджету подстилающей поверхности и степени облачного покрытия. Интенсивность излучения удобно оценивать в ночной период, когда I = 0 в (11) и Я = Г; при этом левую часть последнего равенства зададим из актиномет-рических наблюдений, а правую будем вычислять по формуле Брента
Г = а(а - Ъу/ё2)Т4с, (13)
где а — постоянная Стефана — Больцмана; е =1 — 0.03^н — 0.05^о, Пн, По — балл нижней и общей облачности; е2 — упругость водяного пара при г = 2 м. Имеющаяся информация позволяет рассчитать все компоненты правой части в (13), за исключением агрегата а — Ъ^/е2, в котором коэффициенты а, Ъ будем полагать неизвестными и оптимизировать их значения калибровкой по данным наблюдений. Методическая основа оптимизации заключается в вариации коэффициентов вблизи традиционно применяемых значений [3] и минимизации средней абсолютной ошибки. Были получены значения а = 0.2, Ъ = 0.028.
3. Результаты расчетов
Базовым вариантом модели будем считать задачу (1)—(10) с фиксированным набором параметров, обеспечивающих наиболее близкое к реальному описание метеорологического режима по данным за июль 1989 г. (численный эксперимент 1). Согласно постановке задачи температура и влажность задаются на уровне г = к = 2 м, а на поверхности и в почве определяются в результате расчетов. Учитывая, что температура поверхности была использована как средство оптимизации соотношения (10), оценку точности модельных
6,°С
Фактические (сплошные линии) и расчетные (пунктир) значения температуры (а) и удельной влажности (б) на первые 10 суток интегрирования.
экспериментов проведем путем сопоставления расчетных и фактических величин в, д на уровне г = 0.5 м. На рисунке представлен ход реальной и расчетной температуры и удельной влажности над оголенной (без растительности) поверхностью почвы с дискретностью 3 ч, соответствующей цикличности сбора данных. Видим, что кривые расположены достаточно близко, что свидетельствует о возможности модели правильно воспроизводить структуру приземного слоя. Интегрирование проводилось на срок 30 сут. Средняя абсолютная ошибка в температуре составила 0.8 °С, относительная ошибка в поле влажности равняется 5.4%.
Далее будем анализировать среднемесячные величины. При этом удобной характеристикой является величина суммарного испарения Е, отражающая как термический, так и влажностный режимы поверхности.
На первой стадии (1-7 сут.) скорость испарения регулируется целиком процессами в атмосфере, поскольку начальный запас почвенной влаги в приповерхностном слое не ограничивает потребный поток влажности, обусловленный метеорологическим режимом СПП. По мере просыхания почвы снабжение поверхности водой падает ниже уровня, необходимого для свободного испарения в данных атмосферных условиях. Тогда наступает вторая стадия (£ > 7 сут.), в которой интенсивность испарения ограничена условиями в почве, а именно скоростью передвижения воды с нижних горизонтов к поверхности за счет механизма капиллярной проводимости. На этой стадии влажность поверхности почвы минимальна и близка к гигроскопическому значению, а перенос влаги становится чувствительным к температурным градиентам.
На основе откалиброванной модели была проведена следующая серия численных экспериментов.
Эксперимент 2. Температура воздуха на уровне 2 м увеличивалась на 2 °С по отношению к натурным данным. Цель эксперимента — изучить влияние потепления на температурно-влажностные параметры климата в засушливый период.
Эксперимент 3. Скорость ветра увеличивалась на 50 % от данных. Согласно (5) интенсивность поверхностного испарения пропорциональна скорости ветра в условиях достаточного количества влаги. Поэтому представляет интерес выяснить влияние ветра на испарение в засушливых условиях.
Эксперимент 4- Балл общей и нижней облачности уменьшался на 50 % от фактического. Расчет выполнен для выяснения зависимости решения от интенсивности прямой солнечной радиации и неинструментально наблюдаемой характеристики облачного покрытия. Ненадежность визуального определения облачности может сказаться на итоговых результатах; вследствие этого необходимо изучить чувствительность системы к данному фактору.
Эксперимент 5. Коэффициент гидравлической проводимости почвы увеличивался в 10 раз по сравнению с экспериментом 1. Гидрологические характеристики почвы обладают значительной изменчивостью, и задать их точное значение не представляется возможным. Что касается гидравлической проводимости, то в зависимости от вида почв ее величина может меняться на 4-5 порядков. В связи с этим изменение К на один порядок представляется вполне допустимым.
Эксперимент 6. Коэффициент гидравлической проводимости почвы уменьшался в 10 раз.
Представленные варианты расчетов служат для исследования устойчивости модельного гидрометеорологического режима к вариациям внешних параметров. Однако их можно интерпретировать также как некоторые сценарии реакции климатической системы на изменение природных факторов. Так, например, эксперимент 2 может быть истолкован как попытка климатического прогнозирования влияния повышения глобальной температуры на характеристики локального климата.
Результаты расчетов, осредненные за весь период интегрирования (суммарное испарение не осреднялось), для описанных вариантов представлены в таблице, где По,П5 — влажность почвы на поверхности и на глубине 5 см.
Анализ таблицы показывает, что вариации параметров приземного слоя слабо сказываются на средних характеристиках (эксперименты 1-4). Это вполне объяснимо с физической точки зрения: процесс испарения, протекающий в характерных для выбранного периода условиях дефицита почвенной влаги, не должен зависеть от атмосферного режима, а должен полностью определяться скоростью транспорта почвенной воды к поверхности. Подтверждение этому находим в результатах экспериментов 5, 6. При увеличении влагопроводящей способности почвы (эксперимент 5) суммарное испарение увеличилось почти в пять раз, что связано с повышенной скоростью подвода влаги к поверхности и нелимитированным испарением. Интенсификация испарения вызвала понижение средней температуры на 4 °С вследствие увеличения затрат тепла. В эксперименте 6 картина меняется на противоположную: затрудненный приток влаги к подстилающей поверхности существенно снижает норму испарившейся воды и в некоторые периоды времени даже меняет направление потока влажности — из атмосферы в почву, соответственно увеличивается приповерхностная температура.
Варианты расчетов
Эксперимент 0о, °С Е, см По П5
1 22.7 4.34 0.017 0.109
2 23.7 4.30 0.013 0.098
3 22.3 4.33 0.011 0.104
4 24.8 4.05 0.014 0.105
5 18.8 21.1 0.077 0.115
6 23.6 0.78 0.011 0.167
Представленные результаты дают основания говорить о приемлемости численной модели для описания процессов испарения с поверхности почвы в условиях иссушения. В качестве следующего этапа предполагается провести испытание модели для влажного сезона и сезона со средней нормой осадков. Развитие исследований будет проводиться также путем обобщения параметризации СПП для описания обменных процессов в растительном покрове с учетом транспирации. Данный алгоритм предполагается использовать как компонент моделирующей системы для воспроизведения гидрологического цикла в бассейнах рек Западной Сибири.
Список литературы
[1] Гусев Е.М., Насонова О.Н. Опыт моделирования процессов тепловлагообмена на поверхности суши в региональном масштабе // Водные ресурсы. 2000. Т. 27, № 1. С. 32-47.
[2] ВолодинА Е.Е. Численное исследование чувствительности гидрологических характеристик к вариациям физических параметров системы почва — растительность — снег // Изв. РАН. ФАО. 2001. Т. 37, № 5. С. 700-711.
[3] ПАлАгин Э.Г. Математическое моделирование агрометеорологических условий перезимовки озимых культур. Л.: Гидрометеоиздат, 1981. 191 с.
[4] Пушистов П.Ю., Шлычков В.А. Вывод балансовых соотношений для мезомас-штабной климатической модели тепловлагообмена в системе атмосфера — гидросфера — растительность — почва // Тр. СибНИГМИ. 1992. Вып. 100. С. 129-146.
[5] Лыкосов В.Н. Параметризация пограничного слоя атмосферы в моделях крупномасштабной циркуляции. // Вычисл. процессы и системы. Вып. 10 / Под ред. Г.И. Марчука. М.: Наука, 1993. С. 64-95.
[6] Кучмент Л.С., Демидов В.Н., Мотовилов Ю.Г. Формирование речного стока. М.: Наука, 1983. 216 с.
[7] Братсерт У.Х. Испарение в атмосферу. Теория, история, применения. Л.: Гидро-метеоиздат, 1985. 352 с.
[8] Дрейпер Н., Смит Г. Прикладной регрессионный анализ. М.: Статистика, 1973. 391 с.
Поступила в редакцию 30 июня 2003 г.
