МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ПРИРОДООБУСТРОЙСТВА

 

 
«РОЛЬ ПРИРОДООБУСТРОЙСТВА СЕЛЬСКИХ ТЕРРИТОРИЙ В ОБЕСПЕЧЕНИИ УСТОЙЧИВОГО РАЗВИТИЯ АПК»
 
(МАТЕРИАЛЫ МЕЖДУНАРОДНОЙ НАУЧНО-ПРАКТИЧЕСКОЙ КОНФЕРЕНЦИИ)
 
Москва 2007


УДК 631.6
АНАЛИЗ ГИДРОФИЗИЧЕСКИХ ФУНКЦИЙ В ПРИЛОЖЕНИИ К ПРОГНОЗАМ ВЛАГОПЕРЕНОСА  В ПОЧВАХ

Ф.В. Серебренников
ФГОУ ВПО «Московский государственный университет природообустройства»,
г. Москва, Россия



Both well-known and new  HCR relations are analyzed and their dependence on the capillary-sorbing potential is estimated using Mualem model. The approach on moisture transfer parameters identification is suggested on the base of experimental results.

Для прогнозирования движения влаги в зоне аэрации необходимо обладать достоверными сведениями о гидрофизических свойствах почвы, тем самым требуется знание ряда функций, которые описывают  способность почв удерживать и проводить почвенную влагу под действием термодинамических сил и их градиентов. К числу таких функций относятся: функция водоудерживающей способности почв или основная гидрофизическая характеристика (ОГХ) и функция влагопроводности  ненасыщенных почв. В последнем случае различают капиллярную характеристику влагопроводности (КХВ), связывающую коэффициент влагопроводности почвы и капиллярно-сорбционный потенциал, и влажностную характеристику влагопроводности (ВХВ), связывающую коэффициент влагопроводности почвы и ее влажность. Как ОГХ, так и ВХВ (или КХВ) являются отражением строения порового пространства почвы и поэтому генетически связаны друг с другом. На основе целого ряда представлений о строении порового пространства и феноменологии передвижения влаги в порах почвы были предприняты многочисленные попытки установить аналитически связь между ОГХ и ВХВ (КХВ). Наиболее поздней по времени и в известной степени обобщающей результаты предыдущих исследований является так называемая модель Муалема [1].
В связи со сказанным следует отметить, что определение ОГХ в лабораторных исследованиях выполняется относительно просто, так как соответствующие методики детально разработаны и создано необходимое лабораторное оборудование. Определение ВХВ требует постановки значительно более трудоемких опытов. Поэтому представляется вполне логичной мысль определять ВХВ (аналитически) на основе экспериментальных данных об ОГХ. Для этого требуется подобрать для описания ОГХ некоторую зависимость, отвечающую существу процесса, а затем по модели Муалема найти выражение для ВХВ. Тем не менее, этим путем не так легко воспользоваться, так как модель Муалема зачастую приводит к сложным зависимостям в виде комбинаций специальных функций. Формулы Ван-Генухтен, (в свое время широко рекламируемые), позволяющие осуществлять именно такой переход от ОГХ к ВХВ, как это будет показано ниже, имеют ограниченное применение, несмотря на их простую аналитическую запись.
По этим и другим причинам (например, из-за сложившихся традиций) при прогнозировании водного режима почв очень часто для описания ОГХ и ВХВ используют формулы, не связанные между собой какой-либо моделью строения порового пространства.
Вот почему представляет интерес выполнить анализ получивших к настоящему времени распространение моделей ОГХ и выявить среди них наиболее перспективные как с точки зрения отражения существа процесса влагопереноса в почве, так и с точки зрения удобства организации вычислительного алгоритма.
Ограничимся случаем неизотермического передвижения почвенной влаги, которое с достаточной полнотой описывается уравнением Ричардса-Клютта. Кроме этого, ограничимся рассмотрением ряда самых известных функций водоудерживающей способности почв, привлекших наибольшее внимание специалистов.

В работе [2] для функций водоудерживающей способности почв предложена следующая общая структура их математической записи


,                                                    (1)


где некоторые эмпирические параметры, причем параметр  может быть численно близок значению предельной влагоемкости (W1 по С.Ф. Аверьянову), а параметр  может быть численно близок к значениям максимальной молекулярной влагоемкости (Wr); соответственно, влажность почвы и капиллярно-сорбционный потенциал.
Наиболее полный обзор функций водоудерживающей способности и влагопроводности почв сделан в [1]. Кроме этого, наиболее используемые зависимости ОГХ представлены в работах [2, 3]. Именно эти функции большей частью и применяются в расчетах влагопереноса в зоне аэрации при научном обосновании мелиораций и проектировании мелиоративных систем.
Ниже дается подборка функций, получивших наибольшее распространение
функция Вейбулла [2]


  ;                               (2)


функция А.Я. Почепского  (или Брутшаэра) [2], приведенная по своей структуре к записи (1)


;                                                     (3)


функция Ван-Генухтена [2]


;                                        (4)


функция Фарелла и Янсона [1, 3]


.                                            (5)


Введем обозначение Q (нормирование влажности в интервале от 0 до 1).
Тогда зависимости (2)-(5) можно записать в виде функции Y = f (Q).
Тем самым формулы (2)-(5) преобразуются к виду:


,                                                         (6)
,                                                          (7)
,                                                       (8)
.                                                    (9)


Все эти функции были подвергнуты детальному аналитическому исследованию, а именно, для каждой функции в декартовых координатах строился график для ряда характерных значений параметров  (при В = 1, так как по сути дела В представляет собой масштабный коэффициент). Затем бралась первая производная на предмет изучения поведения функция в окрестности точки (Q = 1, Y = 0). Далее бралась вторая производная для выяснения наличия на кривой точки перегиба. Если таковая имелась, то строился график функции , полученный путем приравнивания второй производной нулю и характеризующей область существования точки перегиба на графике функции   .
При исследовании функции (6) было установлено:
при график функции  имеет точку перегиба, координаты ищутся по зависимости


                                                          (10)


и далее по (6), причем:
при  график функции  вогнутый (в сторону начала координат) и не имеет точки перегиба;
все кривые графики функции   независимо от значения параметра  пересекаются в одной точке с координатой Q = 0,367879... и ; это свойство может быть использовано для первого приближения значений параметров  при наличии соответствующих экспериментальных данных;
точка перегиба (при  смещается по графику в интервале от Q = 0,367879, ;
при  тангенс угла касательной в окрестности точки Q = 1, Y = 0 всегда равен бесконечности;
при   тангенс угла касательной в той же точке всегда равен нулю;
при  график функции вырождается в прямую  в разомкнутом интервале от Q = 0 до Q = 1 (на концах интервала функция не существует).
При исследовании функции (7) было установлено:
эта функция получается из (6) после элементарных преобразований, если воспользоваться первым членом ряда для логарифма в (6)  [4, с. 120, п. 601.5];
при  график функции  имеет точку перегиба; координаты точки перегиба находятся по зависимости  и далее по (7);
при  график функции  вогнутый (в сторону начала координат) и не имеет точки перегиба;
все кривые графика функции пересекаются в одной точке с координатами Q = 0,5 и Y = В, как и ранее, это свойство может быть использовано для первого приближения значений параметров  при наличии соответствующих экспериментальных данных;
точка перегиба (при ) смещается по графику функции   в интервале от Q = 0.5, Y = В и Q = 1, Y = 0;
при  и   тангенс угла касательной в окрестности Q = 1, Y = 0 ведет себя так же, как и для функции (6);
при  график функции  вырождается в прямую Y = В с тем же поведением на концах интервала Q = 0 до Q = 1.
При исследовании функции (8) было установлено:
график функции имеет физический смысл при с > 1;
при с > 1 график функции всегда имеет точки перегиба; координаты точки перегиба находятся по зависимости


                                                            (11)


и затем по (8);
единая точка пересечения кривых графика функции  при В = const отсутствует;
точка перегиба смещается по графику функции Y = f(Q) в интервале  Q = 0.5 , Y = В и Q = 1, Y = 0 (интервал незамкнутый);
при Q ® 1 график функции Y = f (Q) не достигает координаты Q = 1, Y = 0;
тангенс угла наклона касательной в окрестности точки Q = 1, Y = 0 всегда стремиться по модулю к большим значениям (теоретически к бесконечным), но сама производная в этой точке не существует;
при с ® ?  график функции Y = f (Q) вырождается в прямую Y = В  в разомкнутом интервале от Q = 0 до Q = 1 (на концах интервала функция не существует).
При исследовании функции (9) было установлено:
функция существует при 0 < С < 1 (имеет физический смысл);
график функции всегда вогнутый по отношению начала координат;
при Q ® 1 кривые графика функции Y = f (Q) при различных значениях параметра   (С > 0) приближаются к точке с координатами Q = 1, Y = В;
при Q ® 1 первая производная приближается к значению ;
при с ® ? график функции Y = f (Q) вырождается в прямую Y = В  в интервале от Q = 0 до Q = 1 (на концах интервала функция существует).
Представляет интерес оценить свойства функций (6) – (9), сравнивая их поведение с имеющимся экспериментальным материалом [8], а также высказать соображения о физической картине распределения влаги в почве над грунтовыми водами (зона неполного насыщения). С этой целью воспользуемся моделью Муалема [1]


,                                         (12)


где    Q – в данном  случае степень эффективного насыщения;; W – объемная влажность почвы Wr – максимальная гигроскопическая влажность   (доли ед.); WS – влажность полного насыщения (доли ед.);    – коэффициент фильтрации почвы, м/сут.;  kW – коэффициент влагопроводности, м/сут.
Если воспользоваться функциями (6) и (7) и подставить их в (12) для получения аналитической связи kW = f (Q), то возникают определенные затруднения чисто вычислительного свойства. В первом случае функция kW = f (Q) представлена комбинацией неполных гамма-функций, во втором – комбинацией неполных бета-функций [2]. Следовательно, вести вычисления можно только с использованием рядов, которые медленно сходятся и, главное, для таких рядов пока не разработано общедоступных стандартных программ для компьютера. 
Для функции (8) решение существует и довольно простое, так как Ван Генухтен специально подбирал вид функции, который допускал интегрирование по модели (12). Тем не менее, соответствующая формула для расчета коэффициента влагопроводности применима лишь для относительно редкого случая, когда кривая ОГХ обладает точкой перегиба. Что же касается функции (4), то она также позволяет получить простую зависимость для  kW, но вызывает сомнение соответствие функции (9) физике явления: при           Q = 1 (полное насыщение)  y = В, а должно быть y = 0.
Необходимо остановиться еще на одном существенном обстоятельстве.
Когда речь заходит о функциях водоудерживающей способности почв (основной гидрофизической характеристики почв – ОГХ), приходится сталкиваться с широко распространенным мнением, согласно которому вид соответствующих зависимостей отвечает кривой  на рисунке (кривые построены в обычных декартовых координатах). Этот рисунок заимствован из работы  [5, с.38]. Причем кривая  дает представление о принципиальном виде функций водоудерживающей способности почв, однородных по своему механическому составу («хорошо отсортированных», как отмечается в работе [5, с. 39]); кривая  – то же для почв неоднородных по механическому составу («плохо отсортированных», там же). Для последующего изложения очень важным является четкое представление о поведении производной  в точке пересечения кривых ОГХ с координатной осью Q (ось абсцисс): при y ® 0 для кривой  (градиент) существенно возрастает, а для кривой  производная стремится к некоторым минимальным значениям, что подчеркивается в [5, с. 224], а также в [1, с. 63, 64], причем ни бесконечность, ни ноль в реальных условиях не достигаются. Тем самым, формулы (6) и (7)  в точке , равно как и (8), находятся в определенном противоречии с физической картиной ОГХ.
Между тем, в работе [1, с.186, рис. 4.1.1], равно как и во многих других публикациях [10, с.148] в качестве наиболее общей формы ОГХ представлен график y = f (W), вид которого соответствует кривой  на рисунке. Затронутый вопрос важен и с точки зрения оценки характера распределения влаги в зоне аэрации над зеркалом грунтовых вод, поскольку функция ОГХ описывает равновесную эпюру влажности над УГВ.
Чтобы разобраться, насколько для реальных почв типичен тот или иной вид функции ОГХ, нами были привлечены соответствующие обширные экспериментальные материалы, опубликованные в работах [6, 7, 8], а также находящиеся в фондах НПО «ВолжНИИГиМ».
Так, например, А.Н. Морозов [6] для почв среднеазиатского региона приводит многочисленные данные о связи y = f (W). Аналогичные данные для Европейской части России приводит  С.А. Лавров [7]. Такие же сведения содержатся в упомянутой ранее работе [8, с. 58-60].
Для выяснения реального вида графиков функций ОГХ скептикам рекомендуется построить соответствующие кривые по экспериментальным данным не в виде традиционных (что полезно, но может ввести в заблуждение о фактическом распределении напоров y в капиллярной кайме над УГВ), а в обычных декартовых координатах.  Нами такая работа была проделана, что позволило констатировать очевидное: во всех рассматриваемых случаях графическое изображение функций ОГХ соответствует кривой  рисунке (за исключением функции ОГХ для песка, данные С.А. Лаврова).
Следовательно, есть все основания утверждать, что кривая y = f (W) или y = f (Q) преимущественно имеют вогнутый характер (в сторону начала координат) и было бы неправильно недооценивать этот факт.
Высказанные соображения свидетельствуют  о том, что есть основания рассмотреть другие функции ОГХ, которые отвечали бы случаю движения влаги в почвах неоднородного состава.

Кривые водоудерживающей способности

а – почвы однородного механического состава;
б – почвы разнородного механического состава

В качестве таких функций ОГХ рассмотрим следующие, предложенные ранее автором данной работы [ВолжНИИГиМ, 1994, МГУП, отчет о НИР, 1999, 2000, кафедра КИВР]:


,                 (13)
,              (14)


Исследование этих функций проведем по схеме, изложенной ранее.
При исследовании функции (13) было установлено:
функция имеет физический смысл при С > 0;
график функции всегда вогнутый по отношению начала координат;

при Q - 0 кривые графика функции

y = f (Q) при различных значениях параметра С (С > 0) приближаются к точке с координатами Q = 1, y = В и совпадает с ней при Q = 1.
при  Q ® 1 первая производная стремится к значению – В/С (точка Q = 1, y = 0).
При исследовании функции (14) было установлено, что она обладает аналогичными cвойствами, что и (13), несмотря на различия в их записи.
На основании выполненного анализа для определения аналитического выражения kW = f (Q) по модели (12) принимаются зависимости (13) и (14) как отвечающие существу рассматриваемого явления – графического выражения ОГХ.
В модели Муалема (12) обозначим через  J1 неопределенный интеграл, а через J2 – определенный. Тогда на основании функции (12) получаем


,                                     (15)
.                                          (16)


В качестве J2 используются зависимости (15) и (16) стой лишь разницей по сравнению с формулой (1), что под   понимается  ,  где всегда WA < A. При выполнении этого условия отношение интегралов в (12) всегда численно определено.
Таким образом


.                                                  (17)


Если воспользоваться функцией (14), то по модели Муалема получаем


,                                                      (18)


то есть это часть функции (15) без знакопеременного ряда. Для представления выражения kW    в аналитическом виде справедливы все предыдущие рассуждения.
По зависимостям (15)-(18) были проведены вычислительные эксперименты по опытным данным, помещенным в [2, 6, 8]. Сравнение расчетных и фактических значений kW  показало достаточно удовлетворительную точность предложенного подхода для вычисления функции влагопроводности ненасыщенных почв (ВХВ) при наличии экспериментальных данных об основной гидрофизической характеристики (ОГХ). В [9] также было получено подтверждение правомерности использования зависимостей (13)-(18).
Модель Муалема замечательна тем, что она обратима: с ее помощью можно получить и функцию , если известна функция
Так, например, зависимость (12) после очевидных преобразований можно представить в виде дифференциального уравнения    

 
.                                    (19)


Чтобы в качестве демонстрационного примера решить уравнение (19), обратимся к самому простому случаю, а именно воспользуемся уравнением С.Ф.Аверьянова  [10]


kW = ? q 3,52.                                                         (20)


В данном случае формула (20) связывает коэффициент влагопроводности с влажностью, следовательно, (20) представляет собой влажностную характеристику влагопроводности (ВХВ). Уместно отметить, что при расчетах влагопереноса по модели Ричардса-Клютта необходимо знать также функцию водоудерживающей способности почвы, то есть основную гидрофизическую характеристику (ОГХ). Зачастую при построении вычислительного алгоритма используют наряду с формулой С.Ф. Аверьянова в качестве ОГХ уравнения, предложенные А.Я. Поченским, А.М. Якиревичем и др. Разумеется, что параметры в моделях ОГХ, равно как и в ВХВ, должны находиться путем аппроксимации экспериментальных данных, что гарантирует достоверность результатов расчетов, даже если эти зависимости не согласуются с моделью (12). Однако следует отдать предпочтение случаю, когда связь между ОГХ и ВХВ определена в рамках модели Муалема, которая является в значительной мере обобщением результатов предыдущих исследований по данной тематике.
Для общности предадим формуле (20) вид: kW=?q N /Ведерников В.В., 1976/, что позволяет записать , после чего на основании решения уравнения (19) получаем следующую функцию водоудерживающей способности почвы (ОГХ)


.                                                      (21)


Необходимо заметить, что зависимость (21) имеет физический смысл при N ? 3,5.


На основании (21)  можно зависимость (20) записать в виде капиллярной характеристики влагопроводности (КХВ)
.                                             (22)


Учитывая, что N ? 3,5, для большего удобства (21) и (22) записываются как


,                                                   (23)
и
.                                           (24)


Если рассмотреть поведение функции (23) на интервале 0 ? q ?1, то выяснится, что кривая Y=¦(q) представляет собой вогнутую линию по отношению начала координат, что отвечает случаю движения влаги в почвах, неоднородных по механическому составу. Вместе с тем при q  = 1 (полное насыщение) Y = const, что противоречит существу процесса.
Вот почему представляет интерес придать формуле (23) вид, обеспечивающий равенство нулю сорбционно-капиллярного потенциала при полном насыщении почвы влагой, то есть Y = 0 при q = 1. Наиболее просто этого добиться, если кривую Y=¦(q) перенести параллельно оси ординат на шаг, равный , что позволяет записать (6) в виде


.                                                    (25)


В случае (25) при q ® 0, Y ® ?, а при q = 1, Y = 0.
Далее для ОГХ в записи (25) следует найти ВХВ с помощью модели Муалема. Тем самым, речь идет о решении следующего интеграла, пределы интегрирования которого должны быть приняты такими же как и в модели (12)


 или  .                                         (26)


Анализ интеграла (26) показывает, что для ряда частных значений N существуют соответствующие аналитические решения этого интеграла.
Например:


N = 3,5
;                            (27)

N = 4,5
;                                    (28)
N = 6,5
.                                      (29)


В общем случае интеграл (26) выражается с помощью ряда [4]. Для удобства записи в модели Муалема (12), как и ранее  обозначают через J1 неопределенный интеграл (числитель), а через J2 – определенный (знаменатель). Тогда, согласно [4, с.67]


;                         (30)
.                            (31)


В (30) и (32) для удобства записи принято  m = (2N – 5)/4.
Для нескольких частных случаев с помощью (27)-(29) была проверена сходимость рядов (30) и (31), которые присутствуют в (17). Оказалось, что достаточно учесть 50 членов каждого ряда, чтобы вычисления по сравнению с аналитическим решением не отклонялись более чем на 1,0-1,5%.
В заключение уместно обратить внимание на следующее обстоятельство. На основании уравнения (19) для достаточно простого случая  было получено соответствующее аналитическое решение, что в общем случае является исключением, а не правилом, так как запись ВХВ обычно при подстановке ее в (19)  не позволяет решить  дифференциальное уравнение традиционными методами. Однако это не может служить препятствием, если иметь в виду возможности современной вычислительной техники. Действительно, если задаться конкретной функцией , то можно, исходя из уравнения (19), восстановить зависимость , например, с помощью метода Ньютона и далее вести расчеты влагопереноса для каждого шага по времени в рамках обычной вычислительной процедуры.
И еще одно замечание: если вместо записи для ВХВ воспользоваться КХВ, то в уравнение (19) следует заменить переменные, исходя из обычных преобразований, принятых в дифференциальном исчислении, а далее воспользоваться тем же методом Ньютона или аналогичным, исходя  из существа рассматриваемой задачи.

Библиографический список

1.   Глобус А.М. Почвенно-гидрофизическое обеспечение агроэкологических математических моделей. – Л.: Гидрометеоиздат, 1987. 428 с.
2.   Зейлигер А.М. Сопоставление моделей водно-физических характеристик почв с экспериментальными данными. /В кн.: Оптимизация процессов комплексного мелиоративного регулирования. – М., 1985. С. 61-72.
3.   Заславский Б.Г., Опарина И.В., Терлеев В.В. Диалоговая система формирования банка гидрофизических характеристик почв. //Докл. ВАСХНИЛ. 1988. № 11. С. 40-42.
4.   Двайт Г.Б. Таблицы интегралов и другие математические формулы. – М.: Наука, 1964.  228 с.
5.  Бэр Я., Заславски Д., Ирмей С. Физико-математические основы фильтрации воды.                  – М.: Мир, 1971. 452 с.
6.  Морозов А.Н. Исследование возможности полива хлопчатника минерализованными водами (отчет НИР, 1984 г., фондовые материалы института «Средазгипроводхлопок», Ташкент).

7. Лавров С.А. Определение основной гидрофизической характеристики по данным почвенно-гидрологических констант. /В кн.: Вопросы гидрофизики. – Л.: Гидрометеоиздат, 1986. Вып. 308. С. 39-45.

8. Моделирование процессов засоления и осолонцевания почв. – М.: Наука, 1980. 264 с.   

9. Холуденева О.Ю. Автоматизированные технологии ведения комплексного мониторинга орошаемых агроландшафтов Поволжья. Автореф. дис…. канд. техн. наук.                               – Саратов, 2002.

10. Шеин Е.В. Курс физики почв. Учебник. – М.: Изд-во МГУ. 2005. 432 с.