CC BY
152
УДК 551.584.32
Оценка в реальном времени
потока тепла в почву
по данным почвенных термометров
И. Н. Русин, Е. А. Куканова
Описан метод расчета мгновенных значений потока тепла от поверхности вглубь деятельного слоя почвы, а также коэффициента температуропроводности почвенного слоя непосредственно в момент наблюдений. В предлагаемом варианте он позволяет использовать произвольное размещение почвенных термометров по глубинам и любые интервалы регистрации данных по времени. Пригоден для проведения оценок как в условиях стационарных наблюдений, так и в экспедиционных условиях.
Ключевые слова: поток тепла в почву; коэффициент температуропроводности.
Под тепловым балансом деятельной поверхности подразумевается динамическое равновесие потоков энергии, проходящих через эту поверхность. В результате этого равновесия формируются микроклиматические условия жизнедеятельности. Контроль над состоянием энергообмена атмосферы с подстилающей поверхностью является одной из важнейших задач метеорологии. Состояние динамического равновесия потоков энергии на подстилающей поверхности описывается уравнением теплового баланса этой поверхности. Регистрация входящих в уравнение теплового баланса потоков энергии современными измерительными комплексами по всей территории России требует очень больших затрат. Поэтому, несмотря на значительный прогресс в об-
ласти метеорологических наблюдений, наиболее подробные сведения о пространственном распределении составляющих теплового баланса подстилающей поверхности до настоящего времени получают не прямыми измерениями, а косвенным путем, производя расчеты на основе наблюдений сравнительно простыми приборами на сети специализированных метеорологических станций.
Одной из задач регламента теплобалансовых наблюдений является оценка потока тепла в почву. Для этого используются наблюдения за температурой на нескольких уровнях в слое почвы, прилегающем к поверхности. В России эти наблюдения ведутся в синоптические сроки с помощью коленчатых термометров Савинова, устанавливаемых на глубинах 5, 10, 15 и 20 см. Обработка результатов этих наблюдений в соответствии с руководством [Руководство ..., 1977] позволяет получить оценки характеристик почвы и потока тепла в нее.
В связи с давностью издания этого документа [Руководство ..., 1977] и последующим внедрением методов компьютерной обработки данных полезно сравнить рекомендуемые в нем методы с теми, которые разработаны к настоящему времени, и усовершенствовать процедуру расчета потока тепла в почву, учитывая, например, то, что в экспедиционных условиях (также и в других странах мира) сроки наблюдений и глубины установки термометров могут быть различными, а термические характеристики почвы необходимо получать в реальном времени. Рассмотрению этих вопросов и посвящена данная раб^а.
Анализируя современные методы расчета потока тепла в почву, описанные в работах [Gao, 2005; Gao et al., 2009; Surface energy ..., 2004; Understanding ..., 2010], следует отметить, что общепризнанной концепцией описания распространения тепла в почве была и остается модель теплообмена Фурье с постоянным коэффициентом температуропроводности, описываемая уравнением:
к д2Т
(1)
dt dz
Здесь Т — температура почвы, 2 — вертикальная координата, направленная от поверхности почвы вниз, t — время , k — коэффициент температуропроводности.
Конечно, при конструировании математических моделей тепла и влагообмена в почве учитывается возможность существенного изменения этого коэффициента во времени и по глубине. Однако в процессе наблюдений пока нет возможности измерить с нужной для моделирования точностью все необходимые физические характеристики почвы, так что концепция теплообмена Фурье пока остается наиболее приемлемой для оценки потоков по данным стандартных наблюдений.
Наиболее распространенным дополнительным предположением для расчета к является использование частного решения уравнения (1) в форме суточной тепловой волны, распространяющейся от уровня z1 до уровня z2:
В работах [Gao, 2005; Gao et al., 2009] приведено несколько способов получения коэффициента температуропроводности с помощью этих формул. Например, приняв P равным суткам и рассчитав из наблюдений на уровне zt среднюю суточную температуру T а также амплитуду суточного хода температур на обоих уровнях как полуразность максимальной и минимальной температур, можно получить оценку температуропроводности почвы по уменьшению амплитуды и запаздыванию фазы температурной волны с глубиной. Для этого можно использовать либо амплитудное соотношение
Т(z, >t) = Tm + A, sin(® •t + 9,), i = 1,2
к = ^1 - ¿2 )2
2Ь (Д/А2)
(3)
либо фазовое соотношение
, = д(г1 - ¿2)2 " 2(^1 -%)2
(4)
Несмотря на широкое распространение, метод тепловой волны имеет ряд недостатков, которые были указаны еще в первой из посвященных данной теме работ, выполненной в нашей стране Д. Л. Лайх-тманом [Лайхтман, 1947]. Он отмечал, что аппроксимация суточного хода температуры функцией (2) является грубым приближением. Это связано и с нерегулярностью изменений облачности в течение суток, и с синоптическими процессами, которые делают суточные изменения температуры непериодическими. Кроме того, частное решение (2) уравнения теплопроводности почвы неявно предполагает, что коэффициент температуропроводности постоянен не только по вертикали, но и во времени, по крайней мере, в течение суток, что исключает возможность оценки мгновенных значений потока тепла в почву. Для ликвидации указанных недостатков Д. Л. Лайтхман предложил метод, который в дальнейшем был усовершенствован в работе [Цейтин, 1953], проверен и уточнен полевыми испытаниями в работе [Русин Н. П., 1955] и внедрен в оперативную практику [Руководство ... , 1977]. Однако этот метод был жестко привязан к фиксированным глубинам наблюдений и позволял оценивать только средние за некоторый период значения искомых величин. Ниже приведены в усовершенствованной форме основные этапы преобразований уравнения (1) и обновленные рабочие формулы для расчета текущих значений потока тепла в почву и коэффициента температуропроводности в режиме реального времени.
Умножая обе части уравнения (1) на (^-£), где £ — некоторая фиксированная глубина, и используя равенство:
/(гА = <^ +[7(/,0) - Т(1,Ц)],
&
(5)
получим интегральное соотношение, в которое явно входит поток тепла на поверхности почвы Q:
£
ы
(6)
ОН) = ~к
ат
дг
г=0
Выберем два значения глубин £ = h и £ = Н, удовлетворяющих неравенству h<H, и запишем для них полученные интегральные соотношения:
^ ЯТ
ИО, + к[ту,0) - Т0, А)]= Г(г - А)—^
о ^
■ ДО + к [7-((,0) - Т (I, Н )]=/(; - Н )|^<£
(7)
Для получения рабочего алгоритма удобно получить из этой системы линейных относительно Q и к уравнений сначала выражение для потока тепла:
О = ] т( г) ^с1: + [Т ((, А) - Т (1, Н)],
т{ г) =
1,2 < И
{Н - г) 9 г > И (Н - И)
(8)
а затем после небольших преобразований, которые подробно описаны в [5], получить и выражение для k:
к (Н - Н)[т(/,0) - Т(/, К)]- ь[т(/, К) -Т(/,Н)]
Г . ч о! ,
(9)
к =
Формулы (8) и (9) несколько отличаются от выведенных в работе [Цейтин, 1953]. Они предназначены для вычисления мгновенных значений величин, а не осредненных по времени.
Применяя формулы (8) и (9), следует различать случаи, когда расчет потоков производится по полностью сформированному ряду наблюдений и когда обработка ведется в реальном времени. В первом случае предпочтительно использовать прошлые и будущие значения температур и пользоваться симметричными центрально-разностными аппроксимациями. А во втором приходится использовать только прошлые значения и менее точные направленные разностные формулы. В этой работе использованы формулы направленных разностей первого, второго и четвертого порядка аппроксимации для равноотстоящих моментов:
(10)
д/
3/(О - 4/6 ) + / (?-)
(11)
2Д
д/
25/+1) - 48/(;г) + 36/) -16/(/,) + з/(;г )
12А
Для вычисления интегралов по глубине удобно применить метод, предложенный в работе [Цейтин,1953], т. е. воспользоваться интерполяционным полиномом, построенным по наблюдениям за температурой на глубинах. Поскольку в упомянутой работе этот метод был привязан к фиксированным глубинам 5, 10 ,15 и 20 см, то ниже он обобщен так, чтобы использовать любые глубины или моменты времени.
Предположим, что наблюдения за температурой проведены на поверхности почвы ^=0) и на п глубинах ..., zn). В качестве примера использовано п = 4. Тогда коэффициенты а1 интерполяционного полинома
дТ ( г)
I ( г ) =
дг
= 1
а. ■ г
(13)
2=0
можно представить в виде решения линейной системы уравнений:
М • а = /
М =
"1 *1 21 21 (а > и 0 ( / ( 2г) ^
1 2 2 2 2 2 2 2 2 а1 / ( г 2)
1 2 3 2 ¿3 3 гз 4 23 а = а2 7 = / ( 2з)
1 2 4 2 * 4 3 2 4 4 2 4 а3 / ( г 4)
1 2 5 2 ^5 3 ^5 4 25 _ V а4 1 / ( ^ з) J
(14)
к=5
а = М-х • / ^аг =2 тик/(г, )
к=1
м _1 = К* }
Здесь, помимо очевидных рабочих обозначений, введено обозначение т,к для элементов матрицы, обратной матрице Вандермонда, построенной по заданным глубинам наблюдений. Использование этих коэффициентов позволяет получить интегралы по глубине, входящие в формулы (8) и (9) в виде сумм:
к=1
Ztt
dt
1 -
н
Н - h i=0 (i +1)(/ + 2) íи (z )f^z = Ё D*(z*)
к=1
d = —— y 'w
k H - h ¿(/ +1)(/ + 2)
1 -
h<
v i+1
(15)
Формулы (8), (9) и (15) позволяют произвести все необходимые расчеты для произвольного количества глубин наблюдений, расположенных на разных расстояниях друг от друга и с разными шагами по времени. Если число точек по вертикали отличается от пяти, то соответственно следует изменить порядок полинома и матрицу Вандермонда.
Для иллюстрации применения метода на рис. 1 приведены результаты расчетов потока тепла в почву по программе наблюдений, отличающейся от принятой в России. Наблюдения проведены в штате Вермонт, США [Peck, Fiori, 1992]. Они относятся к июлю 1992 года. Интервалы между сроками наблюдений — 30 минут. Использованы данные о температуре поверхности почвы на принятых в США глубинах (7.5, 15, 22.5, 30, 37.5, 45, 52.5 и 60 см).
Во всех вариантах, кроме Г, расчет потока проведен по формулам (8), (9) и (15). Варианты А, Б и В различаются способом оценки производных: A — по формуле (12), Б — по формуле (11), В — по формуле (10). В этих вариантах в качестве нижней границы H выбрано 30 см, а в качестве h — середина слоя, т. е. 15 см. Вариант Г вычислен способом, использованным в работах [Gao, 2005; Gao et al., 2009], то есть на основе определения Q по Фурье при коэффициенте темпе-
ратуропроводности, полученном по затуханию фазы суточной волны, по формулам (12—14) и оценке вертикального градиента с помощью разностей температуры на поверхности и на глубине 7.5 см. Вариант АИ аналогичен А, но при h = 7.5 см. Вариант АН аналогичен А, но при Н = 7.5 см.
150
Рис. 1. Потоки тепла в почву, вычисленные по данным [9].
А — базовый вариант с применением формулы (14),
Б — с применением формулы (13), В — с применением формулы (12),
Г — рассчитанный по методике [6; 7], АИ — базовый вариант при И=7.5 см,
АН — АИ — базовый вариант при Н=22.5 см
Рисунок 1 позволяет убедиться, что результаты расчетов потока тепла по этому методу несущественно реагируют на замену формул аппроксимации производных по времени. Среднее значение разностей между А и Б составляет 0.17 Вт/м2 при среднем квадрате разностей 16.4 Вт/м2, а для разностей А и В соответственно — 16.7 и 23.2 Вт/м2. Это показывает, что использование формулы (10) заметно снижает точность расчета, а различие между формулами (11) и (12) для расчета потока тепла в почву несущественно.
Влияние выбора h и H на значения потоков несколько сильнее: Среднее значение разностей между А и Ah составляет 8.2 Вт/м2 при среднем квадрате разностей 81.1 Вт/м2. Это означает, что глубину h следует выбирать вблизи средней глубины установки термометров. Среднее значение разностей между А и AH составляет 16.1 Вт/м2 при среднем квадрате разностей 197.2 Вт/м2 Причина в том, что выбор такого значения H уменьшил общее количество глубин наблюдений с четырех до трех. В стандарте наблюдений это не допускается, и число измерений температуры на глубинах не может быть меньше четырех. Однако даже в таком искусственном случае фаза и амплитуда суточного хода потока тепла в почву передается верно.
Такого утверждения нельзя сделать, сравнивая варианты А и Г. Можно видеть, что при замене производной по вертикали конечной разностью, как это предполагается в методе [Gao, 2005; Gao et al., 2009], фаза потока тепла смещается и приближается к фазе суточного изменения температуры поверхности почвы. Это понятно, потому что изменчивость температуры поверхности значительно сильнее изменчивости температуры даже на самой малой из глубин измерения. Но физический анализ и известные данные наблюдений [Лайхтман, 1947; Русин Н. П., 1955] показывают, что поток тепла в почву опережает по фазе изменение поверхностной температуры. Это и приводит к большим средним разностям значений между А и Г; в нашем случае разность равна 37.8 Вт/м2 при средней квадратичной разности более 3000 Вт/м2. Отсюда следует, что метод [Gao, 2005; Gao et al., 2009] пригоден для расчета потока тепла по наблюдениям за температурами, только если устанавливать один из почвенных термометров на глубинах 1—2 см, что весьма затруднительно на практике.
К достоинствам рассматриваемого здесь метода относится также и возможность определять текущие значения коэффициента температуропроводности почвы. Это создает предпосылки для изучения влия-
ния смены погоды и условий увлажнения на его значения. Однако, как показано на рис. 2, значения к весьма чувствительны к выбору способа оценки производной по времени. На этом рисунке показаны результаты расчетов по формулам (9) и (15) при двух способах конечно-разностного дифференцирования по времени. Вариант А соответствует пятиточечной формуле (12), вариант Б — менее точной формуле (11).
Основной вывод, который позволяет сделать рис. 2, состоит в том, что при достаточной точности дифференцирования по времени (А) коэффициент к не меняется в течение суток, как это и предполагается при использовании метода суточной тепловой волны (2) — (4). Однако это лишь в рассмотренные сутки. Поскольку рассматриваемый метод не предполагает постоянство коэффициента в течение суток, он позволяет изучать его изменчивость в зависимости от изменения влажности почвы при выпадении осадков или продолжительном испарении. Это сильная сторона метода.
Сравнивая случаи А и Б (рис. 2), можно видеть, что для получения качественной информации о коэффициенте температуропроводности следует применять формулы численного дифференцирования по времени повышенной точности. Физические условия для этого обсуждены выше, при выводе критерия (11). Тот факт, что применение формул повышенной точности повышает качество расчета k, означает, что существуют моменты суток, когда близко к нулю изменение теплосодержания всего деятельного слоя почвы. Это значит, что изменения температуры в верхней и нижней части деятельного слоя противоположны по знаку друг другу. Эти моменты близки ко времени рассвета или заката, когда температура почвы мало меняется по вертикали, что приводит к росту погрешности знаменателя формулы (9) и понижает точность оценок. После прохождения этих моментов в утренние часы нижняя часть деятельного слоя охлаждается, а верхняя нагревается, что обусловливает резкий рост вертикального градиента температуры
почвы и, соответственно, потока тепла вглубь почвы, а значит, быстрый прогрев. Вечером происходит противоположный процесс, почва быстро охлаждается. Этим объясняется и опережение — тот факт, что поток тепла в почву опережает по фазе изменения температуры ее поверхности.
Рис. 2 Коэффициенты температуропроводности, полученные в процессе вычисления потоков тепла в почву. Рассчитаны для тех же суток, что и потоки тепла, значения которых приведены на рис. 1.
В качестве примера применения метода в разных условиях наблюдений были произведены расчеты потока тепла для трех выборок. Первый массив наблюдений получен в штате Вермонт, США [Peck et al., 1992]. Он описан выше. Второй массив — данные наблюдений, выполненных в экспедиции Главной геофизической обсерватории имени А. И. Воейкова (ГГО) летом 1961 года в Днепропетровской области [Лайхтман, Орленко, 1963]. Наблюдения проведены с интервалом в два часа. Глубины установки термометров: 0, 5, 10, 15, 20, 40, 50 см. Третий массив получен во время учебно-научной практики студентов кафедры метеорологии и климатологии географического факультета
МГУ летом 2011 года недалеко от города Кировска в Хибинах. Наблю -дения проводились каждые 4 часа с помощью термометров, установленных на стандартных глубинах (5, 10, 15 и 20 см), и поверхностного термометра. При обработке наблюдений глубина Н выбрана так, чтобы изменения температуры на ней превышали погрешность измерения (в данном случае более 0.1ОС). Для первого массива такой глубиной является 30 сантиметров. Для второго и третьего — 20 см. Результаты расчетов приведены далее (рис. 3).
Рис. 3. Потоки тепла в почву, рассчитанные по трем выборкам данных, полученных в разных физико-географических условиях и при разных режимах наблюдений: А — Вермонт, США, Б — Днепропетровская область, В — Хибины.
По оси x — часы, отсчитываемые от начала наблюдений
Рисунок демонстрирует большие значения потока тепла в почву, его выраженный суточный ход и междусуточную изменчивость. Поэтому можно заключить, что поток тепла в почву может играть важную роль в энергообмене атмосферы с подстилающей поверхностью и требует систематического мониторинга во всех микроклиматических экспериментах.
Рассмотренный метод позволяет в режиме реального времени получать данные о потоке тепла в почву и коэффициенте температуропроводности не только на специализированных метеорологических станциях, но и в нестандартных условиях экспедиционных наблюдений.
Литература
1. Лайхтман Д. Л. Новая формула для вычисления теплового потока в почве по экспериментальным данным / Д. Л. Лайхтман // Труды НИУ ГУГМС. Сер. 1. - 1947. - Вып. 39. - С. 23-26.
2. Лайхтман Д. Л. Программа и методика экспедиционных наблюдений. Приложение 4 / Д. Л. Лайхтман, Л. Р. Орленко // Труды ГГО им. А. И. Воейкова. - 1963. - Вып. 144. - С. 3-8, с. 234-235.
3. Руководство по теплобалансовым наблюдениям. - Ленинград : Гидрометеоиздат, 1977. - 146 с.
4. Русин Н. П. Об определении теплообмена в почве на гидрометеорологических станциях / Н. П. Русин // Труды ГГО им. А. И. Воейкова. -1955. - Вып. 52(114). - С. 33-40.
5. Цейтин Г. Х. К вопросу об определении некоторых тепловых свойств почвы / Г. Х. Цейтин // Труды ГГО им. А. И. Воейкова. - 1953. -Вып. 39(101). - С. 201-213.
6. Gao Z. Q. Determination of soil heat flux in a Tibetan short-grass prairie // Boundary-Layer Meteorology / Z. Q. Gao. - 2005. - Volume 114. -Issue 1. - P. 165-178.
7. Gao Z. Comparison of six algorithms to determine the soil thermal diffusivity at a site in the Loess Plateau of China / Z. Gao, L. Wang; R. Horton // Hydrology and Earth System Sciences Discussions. - 2009. -Volume 6. - Issue 2. - P. 2247-2274.
8. Peck L. 1992. Soil temperatures at South Royalton, Vermont, USA. Boulder, CO : National Snow and Ice Data Center / World Data Center for Glaciology. Digital Media / L. Peck, J. Fiori. - Режим доступа : http://nsidc. org/data/ggd489.html.
9. Surface energy balance closure in an arid region: role of soil heat flux / B. G. Heusinkveld, A. F. G. Jacobs, A. A. M. Holtslag, S. M. Berkowicz // Agricultural and Forest Meteorology. - 2004. - Volume 122. - P. 2137.
10. Understanding heat transfer in the shallow subsurface using temperature observations / M. M. Rutten, S. C. Steele-Dunne, J. Judge, N. van de Giesen // Vadose Zone Journal. - 2010. - Volume 9. - No. 4. -P. 1034-1045.
© Русин И. Н., Куканова Е. А., 2012
Real-Time Estimation of Heat Flow in Soil According to Soil Thermometer Readings
I. Rusin, E. Kukanova
The article gives a description of the calculation method for instantaneous values of the heat flow from the surface down to the active soil layer as well as for the coefficient of soil layer temperature conductivity right at the moment of observation. In the provided variant it allows random depth placement of soil thermometers and using any time intervals of data recording. This method can be used to carry out estimations in conditions of both stationary observation and field research.
Key words: heat flow in soil; temperature conductivity coefficient.
Русин Игорь Николаевич, доктор географических наук, профессор кафедры климатологии и мониторинга окружающей среды, факультет географии и геоэкологии, Санкт-Петербургский государственный университет (Санкт-Петербург), inrusin@mail.ru.
Rusin, I., PhD in Geography, professor, Department of Climatology and Environmental Monitoring, Faculty of Geography and Geo-Ecology, Saint Petersburg State University (Saint Petersburg), inrusin@mail.ru.
Куканова Евгения Андреевна, студентка кафедры метеорологии и климатологии, географический факультет, Московский государственный университет им. М. В. Ломоносова (Москва), vivera5@ya.ru.
Kukanova, E., student, Department of Meteorology and Climatology, Faculty of Geography, Lomonosov Moscow State University (Moscow), vivera5@ya.ru.
