Почвоведение, 2021, № 6, стр. 715-724

Функциональное представление водоудерживающей способности и относительной гидравлической проводимости почвы с учетом гистерезиса

В. В. Терлеев a*, Р. С. Гиневский a, В. А. Лазарев a, А. Г. Топаж b, Е. А. Дунаева c

a Санкт-Петербургский политехнический университет Петра Великого
195251 Санкт-Петербург, Политехническая ул., 29, Россия

b ООО “Бюро Гиперборея”
193312 Санкт-Петербург, ул. Подвойского, 40-2, Россия

c Научно-исследовательский институт сельского хозяйства Крыма
295543 Симферополь, ул. Киевская, 150, Россия

* E-mail: Vitaly_Terleev@mail.ru

Поступила в редакцию 10.07.2020
После доработки 03.11.2020
Принята к публикации 20.11.2020

Полный текст (PDF)

Аннотация

Целью исследования является моделирование водоудерживающей способности и относительной гидравлической проводимости почвы как капиллярно-пористой среды, а также верификация предлагаемых моделей в сравнении с наиболее известными мировыми аналогами. Указанная цель достигается решением следующих задач: 1) описанием гидрофизических свойств почвы в виде трех систем функций с соответствующими наборами общих параметров; 2) верификацией этих систем путем оценивания относительной гидравлической проводимости, а также сканирующих ветвей водоудерживающей способности с использованием параметров, идентифицированных по данным из литературного источника о главных ветвях гистерезиса водоудерживающей способности опесчаненного суглинка; 3) применением равенства значений экспоненциального параметра в вычислениях ветвей иссушения и увлажнения гистерезиса водоудерживающей способности почвы для устранения нежелательного искусственного “эффекта помпы”; 4) исследованием влияния аддитивного параметра на погрешности точечной аппроксимации данных о главных ветвях водоудерживающей способности, а также на погрешности оценок относительной гидравлической проводимости и сканирующих ветвей гистерезиса водоудерживающей способности почвы, а также 5) выявлением достоверных различий между погрешностями этих оценок по критерию Вильямса–Клута для выбора лучшей системы функций. В поисках решения проблем точного ирригационного земледелия, таких как прогнозирование влагообеспеченности сельскохозяйственных культур и расчет прецизионных норм орошения, применение предлагаемых авторами моделей представляется наиболее предпочтительным.

Ключевые слова: математическая модель, точечная аппроксимация, “эффект помпы”, гидрофизические свойства почвы

ВВЕДЕНИЕ

Насыщенность почвы водой, имеющей свойства жидкости, оценивается эффективным влагонасыщением Se = (θ – θr)/(θs – θr), где θ (см3 · см–3) – объемная влажность; θs (см3 · см–3) – объемная влажность насыщения; θr (см3 · см–3) – остаточная влажность, соответствующая минимальному удельному объему воды как жидкости в почве. К гидрофизическим свойствам почвы относятся ее водоудерживающая способность и гидравлическая проводимость. Водоудерживающая способность описывается в виде зависимости Se или θ от капиллярного давления (капиллярно-сорбционного потенциала) почвенной влаги ψ (см вод. ст.) (ψ < 0 в почвах, не насыщенных водой). Для описания гидравлической проводимости используется зависимость коэффициента влагопроводности почвы k (см · сут–1) от величин Se, θ или ψ. Максимальное значение k равно коэффициенту фильтрации влаги в почве ks (см · сут–1). Отношение k/ks называется относительной гидравлической проводимостью почвы. Неоднозначный характер зависимостей, описывающих гидрофизические свойства почвы в виде функций, аргументами которых является ψ, обусловлен феноменом гистерезиса [22, 25, 27, 28].

В литературе приведены примеры описания зависимости θ(ψ) для различных интервалов значений ψ на основе физических представлений, например – экспоненциальная модель [5, 6, 8]. Кривая, которая графически представляет экспоненциальную модель, не имеет S-образной формы в общем диапазоне пленочной и капиллярной влаги; для этого диапазона почвенной влаги к настоящему времени не предложено исчерпывающего обоснованного математического описания гидрофизических свойств почвы в виде обобщенной системы функций [10]. Поэтому при формулировании зависимостей Se(ψ) и k(Se) обычно используют регрессионные и другие эмпирические функции [2, 3, 7]. В большинстве случаев данные функции применяются для точечной аппроксимации данных прямых измерений, например – с применением алгоритма Левенберга–Марквардта для минимизации квадратичных отклонений расчетных значений от опытных данных [21, 23]. Как известно, эти измерения являются весьма трудоемкими. Их трудоемкость особенно возрастает при измерении сканирующих ветвей гистерезиса. Вместе с тем для решения ряда задач гидрофизики почвы требуются не только данные прямых измерений, но и функциональное представление зависимостей Se(ψ) и k(Se). К числу таких задач относятся, например, прогнозирование влагообеспеченности сельскохозяйственных культур и расчет прецизионных норм орошения.

Для решения первой задачи применяется дифференциальное уравнение в частных производных параболического типа (уравнение Ричардса) [32], одним из коэффициентов которого является функция приведенной дифференциальной влагоемкости почвы dSe/dψ. Использование регрессионных зависимостей (в том числе степенных полиномов) и других эмпирических функций [2, 8], которые аппроксимируют опытные данные, нередко приводит к физически абсурдным результатам, например, при вычислении производной dSe/dψ, поскольку дифференцирование не является устойчивой операцией по отношению к аппроксимациям. Эта проблема особенно обостряется, если для решения уравнения Ричардса применяются численные методы с итерационными процедурами [30].

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

По причине высокой трудоемкости измерения сканирующих ветвей в ирригационном земледелии нередко применяется методика расчета норм орошения, основанная на разности между наименьшей влагоемкостью и предполивной влажностью почвы: при этом наименьшая влагоемкость обычно определяется главной ветвью иссушения Se(ψ). На этой ветви при критическом значении ψ, которому соответствует наибольший запас влаги, удерживаемый почвой, величина Se достигает максимального значения в сравнении с любой другой ветвью гистерезиса. Однако смене состояний почвенной влаги при увлажнении почвы соответствует не главная ветвь иссушения, а определенная сканирующая ветвь увлажнения [35]. Поэтому при расчете нормы орошения по данной методике результат может оказаться завышенным. Для получения более точного результата следует использовать определенное (более низкое) значение Se на соответствующей сканирующей ветви увлажнения. Попутно отметим, что для определения этого (более низкого) значения Se может быть использована эмпирическая зависимость Воронина [1].

Наряду с трудоемкостью измерения сканирующих ветвей Se(ψ), существует проблема непредсказуемости выпадения атмосферных осадков в предстоящем вегетационном периоде (в условиях сельскохозяйственного поля): невозможно предугадать, какие именно сканирующие ветви увлажнения потребуются для вычисления прецизионных норм орошения. В данном случае применение математической модели гистерезиса Se(ψ) не имеет альтернативы. Поэтому в значительной мере актуальными являются исследования, направленные на разработку методов, которые позволяют уменьшить объем прямых измерений гидрофизических свойств почвы [40], а также на построение физически обоснованных математических моделей, которые необходимы для прогнозирования влагообеспеченности сельскохозяйственных культур и расчета прецизионных норм орошения.

Цель исследования – функциональное описание гидрофизических свойств почвы как капиллярно-пористой среды в форме математических моделей, его верификация и сравнение с аналогами на примере данных из авторитетного литературного источника [24]. Указанная цель достигается решением следующих задач:

• описанием зависимостей Se(ψ) и k(Se)/ks в виде трех систем функций с соответствующими наборами общих параметров (применение общих параметров позволяет оценивать значения функции k(Se)/ks по опытным данным Se(ψ) и уменьшить объем прямых измерений);

• учетом гистерезиса и представлением зависимости Se(ψ) в виде функции, которая имеет два набора параметров: один набор – для ветвей иссушения, другой – для ветвей увлажнения (применение таких наборов параметров позволяет оценивать сканирующие ветви по данным о главных ветвях гистерезиса Se(ψ) и уменьшить объем прямых измерений);

• идентификацией параметров трех систем функций путем точечной аппроксимации данных из литературы о главных ветвях гистерезиса Se(ψ);

• верификацией трех систем функций путем оценивания значений k(Se)/ks и сканирующих ветвей Se(ψ) с использованием параметров, идентифицированных по данным о главных ветвях гистерезиса Se(ψ);

• выявлением достоверных различий между погрешностями оценок по критерию Вильямса–Клута [4] для выбора лучшей модели.

ОБЪЕКТЫ И МЕТОДЫ

В настоящее время для описания гидрофизических свойств почвы широко применяются функции водоудерживающей способности и относительной гидравлической проводимости почвы, предложенные Ван Генухтеном [39]:

(1)
${{S}_{e}}(\psi ) = \left[ \begin{gathered} {{(1 + {{( - \alpha \psi )}^{n}})}^{{ - (1 - {1 \mathord{\left/ {\vphantom {1 {n)}}} \right. \kern-0em} {n)}}}}},\,\,\,\psi < 0; \hfill \\ 1,\,\,\psi \geqslant 0, \hfill \\ \end{gathered} \right.$
(2)
${{k({{S}_{e}})} \mathord{\left/ {\vphantom {{k({{S}_{e}})} {{{k}_{s}}}}} \right. \kern-0em} {{{k}_{s}}}} = \left[ \begin{gathered} \sqrt {{{S}_{e}}} {{(1 - {{(1 - S_{e}^{{{1 \mathord{\left/ {\vphantom {1 {(1 - {1 \mathord{\left/ {\vphantom {1 n}} \right. \kern-0em} n})}}} \right. \kern-0em} {(1 - {1 \mathord{\left/ {\vphantom {1 n}} \right. \kern-0em} n})}}}})}^{{({{1 - 1} \mathord{\left/ {\vphantom {{1 - 1} {n)}}} \right. \kern-0em} {n)}}}}})}^{2}},\,\,\,\theta < {{\theta }_{s}}; \hfill \\ 1,\,\,\,\theta = {{\theta }_{s}}, \hfill \\ \end{gathered} \right.$
где n > 1 и α (см вод. ст.–1) – формальные параметры.

Метод оценивания значений функции относительной гидравлической проводимости по данным о водоудерживающей способности почвы с использованием соотношений (1) и (2) называется методом Муалема–Ван Генухтена [26, 39]. Для функций, которые описываются соотношениями (1) и (2), воспользуемся обозначениями соответственно: WRC-VG (WRC – water retention capacity) и RHC-MVG (RHC – relative hydraulic conductivity). Важным достоинством функций WRC-VG и RHC-MVG является то, что они имеют общие параметры. Но формальный характер параметров n и α является причиной ряда недостатков этих функций. Первым недостатком является возможность возникновения упомянутой выше проблемы дифференцирования аппроксимации. Второй недостаток состоит в проблематичности оценки RHC-VG за пределами интервала ψ, на котором измерена зависимость θ(ψ). Третий недостаток заключается в ограничительном условии n > 1, которое приводит к тому, что при малых значениях параметра n имеет место значительное возрастание погрешности оценок значений функции относительной гидравлической проводимости по данным о водоудерживающей способности почвы с использованием формул (1) и (2), отмеченное Ван Генухтеном в отношении аллювиальной глинистой почвы 1006 Beit Netofa clay из каталога Муалема (ks = 9.5 × 10–7 cм · с–1) [24, 31].

Следует подчеркнуть, Ван Генухтен не связывал сомнительный результат для этой почвы с использованием условия n > 1. Однако, как показано в статье [9], именно оно является одной из причин достаточно высокой погрешности метода Муалема-Ван Генухтена для почвы 1006 Beit Netofa clay.

В работах [29, 36, 37] на основе представлений о почве как капиллярно-пористой среде сформулированы гидрофизические функции почвы, которые здесь перепишем в виде следующих соотношений:

(3)
${{S}_{e}}(\psi ) = \left[ \begin{gathered} \frac{1}{2}{\text{erfc}}\left( {\frac{{n\sqrt \pi }}{4}\ln ( - \alpha (\psi - {{\psi }_{e}}))} \right),\,\,\,\psi < {{\psi }_{e}}; \hfill \\ 1,\,\,\,\psi \geqslant {{\psi }_{e}}, \hfill \\ \end{gathered} \right.$
(4)
$\begin{gathered} {{k({{S}_{e}})} \mathord{\left/ {\vphantom {{k({{S}_{e}})} {{{k}_{s}}}}} \right. \kern-0em} {{{k}_{s}}}} = \\ = \,\,\left[ \begin{gathered} \frac{{\sqrt {{{S}_{e}}} }}{4}{{\left( {{\text{erfc}}\left( {{\text{inverfc}}(2{{S}_{e}}) + \frac{2}{{n\sqrt \pi }}} \right)} \right)}^{2}},\,\,\,\theta < {{\theta }_{s}}; \hfill \\ 1,\,\,\,\theta \geqslant {{\theta }_{s}}, \hfill \\ \end{gathered} \right. \\ \end{gathered} $
где erfc(x) = 1 – $\frac{2}{{\sqrt \pi }}\int_0^x {\exp ( - {{t}^{2}})dt} $ – дополнительная функция ошибок; inverfc(erfc(x)) = x; ψe (см вод. ст.) – аддитивный параметр: с учетом гистерезиса ψe – капиллярное давление входа воздуха (давление барботирования) для ветвей иссушения (ψe ≤ 0) и ψe – капиллярное давление входа воды для ветвей увлажнения (ψe ≥ 0); α (см вод. ст.–1) – мультипликативный параметр: α = –1/(ψ0 – ψe), ψ0 (см вод. ст.) – капиллярное давление влаги, при котором плотность распределения вероятностей по значениям нормально-распределенной случайной величины ln((ψ – ‒ ψe)/(ψ0 – ψe)) с нулевым генеральным средним и стандартным отклонением σ достигает максимального значения, ψ0 < ψe; n – экспоненциальный параметр: n = ${4 \mathord{\left/ {\vphantom {4 {(\sigma \sqrt {2\pi } )}}} \right. \kern-0em} {(\sigma \sqrt {2\pi } )}}.$

Для функций, которые описываются соотношениями (3) и (4), воспользуемся обозначениями соответственно: WRC-KT и RHC-MKT. В частном случае при ψe = 0 функции WRC-KT и RHC-MKT сводятся к моделям, предложенным Косуги [1720], которые обозначим соответственно WRC-KT0 и RHC-MKT0.

В статье [9] представлена непрерывная аппроксимация соотношений (3) и (4) в классе элементарных функций:

(5)
${{S}_{e}}(\psi ) = \left[ \begin{gathered} {{(1 + {{( - \alpha (\psi - {{\psi }_{e}}))}^{n}})}^{{ - 1}}},\,\,\psi < {{\psi }_{e}}; \hfill \\ 1,\,\,\,\psi \geqslant {{\psi }_{e}}, \hfill \\ \end{gathered} \right.$
(6)
${{k({{S}_{e}})} \mathord{\left/ {\vphantom {{k({{S}_{e}})} {{{k}_{s}}}}} \right. \kern-0em} {{{k}_{s}}}} = \left[ \begin{gathered} \sqrt {{{S}_{e}}} {{\left( {1 - (1 - S_{e}^{{ - 1}})\exp \left( {\frac{8}{{n\pi }}} \right)} \right)}^{{ - 2}}},\,\,\,\theta < {{\theta }_{s}}; \hfill \\ 1,\,\,\,\theta \geqslant {{\theta }_{s}}, \hfill \\ \end{gathered} \right.$
где ψe, α и n – те же параметры, что и в соотношениях (3) и (4).

Для функций, которые описываются соотношениями (5) и (6), воспользуемся обозначениями соответственно: WRC-НT и RHC-MT. В частном случае при ψe = 0 функция WRC-НT сводится к модели, предложенной Хаверкампом и соавторами [14], которую обозначим WRC-НT0. Для функции RHC-MT при ψe = 0 воспользуемся обозначением RHC-MT0. Стоит отметить, что пара функций WRC-НT0 и RHC-MT0 представляет собой математически корректное решение задачи Ван Генухтена в ее исходной постановке.

В статье [11] отмечается, что вопрос интерпретации мультипликативного параметра функции (1) требует дальнейших исследований. Действительно, для этого параметра в литературе нередко встречается весьма сомнительная интерпретация, в соответствии с которой он представляет собой величину, обратную давлению барботирования. Как известно, под давлением барботирования понимается такое капиллярное давление влаги, с которого начинается вход воздуха при иссушении изначально влагонасыщенной почвы, когда еще выполняется равенство Se = 1. Нетрудно убедиться, что для функции (1) при ψ = –1/α это равенство не выполняется, поскольку: Se = 2–(1 – 1/n) и применяется условие n > 1. Поэтому степень обоснованности такой интерпретации представляется более чем проблематичной. Проблему интерпретации мультипликативного параметра α функции (1) как величины, обратной давлению барботирования, иллюстрируют и опытные данные, приведенные в статье [11]. Отметим, что при точечной аппроксимации этих данных с использованием функций (3) и (5) величину давления барботирования может адекватно учитывать аддитивный параметр ψe.

Значительная часть исследований по моделированию гистерезиса Se(ψ) является развитием двух известных разработок: (i) модели Скотта и соавторов [33], а также (ii) модели Кула и Паркера [16]. В первой модели используется функция WRC-НT0; в основу второй модели положена функция WRC-VG. Воспользуемся следующими обозначениями: Hys-SHT0 – для модели Скотта и соавторов; Hys-KPVG – для модели Кула и Паркера.

Функция WRC-VG, применяемая в модели Hys-KPVG, определена только на интервале ψ ≤ 0. Но поскольку главные ветви иссушения и увлажнения гистерезиса могут смыкаться при ψ > 0 (в области вытеснения воздуха, защемленного в тупиковых порах) [12, 13, 15], постольку функция WRC-VG принципиально не может в полной мере описать феномен гистерезиса.

Наряду с моделями Hys-SHT0 и Hys-KPVG в данном исследовании рассматриваются еще три математические модели гистерезиса. Для этих моделей введем следующие обозначения: Hys-SKT, которая основана на функции WRC-KT, описываемой соотношением (3); Hys-SKT0, которая базируется на функции WRC-KT0, описываемой соотношением (3) при ψe = 0; Hys-SHT, в основу которой положена функция WRC-HT, описываемая соотношением (5).

Во всех пяти обозначенных моделях расчет сканирующих ветвей осуществляется по алгоритму, разработанному Скоттом и соавт. [33]. Для отмеченных пяти моделей гистерезиса существует возможность возникновения нежелательного искусственного (методического) “эффекта помпы”. Эффект заключается в том, что при осцилляции ψ в фиксированном интервале возможно пересечение сканирующей и главной ветвей гистерезиса, а также достижение величины Se физически абсурдных значений. По мнению авторов данной статьи: во-первых, пересекаться могут только сканирующие ветви; во-вторых, от главных ветвей начинаются сканирующие (первичные) ветви, но на главных ветвях не могут заканчиваться сканирующие ветви; в-третьих, между двумя точками пересечения двух соседних в последовательной очередности сканирующих ветвей может образовываться замкнутая петля; в-четвертых, в каждой точке на любой ветви производная dSe/dψ принимает только два значения, которые соответствуют ветвям иссушения и увлажнения. Можно предположить, что наиболее предпочтительным способом предотвращения возникновения “эффекта помпы” является применение условия равенства значений экспоненциального параметра n для ветвей иссушения и увлажнения: nd = nw (здесь и далее индекс d используется для ветвей иссушения, а индекс w – для ветвей увлажнения).

Приведенные выше функции Se(ψ) и k(Se)/ks, а также модели гистерезиса Se(ψ) сгруппируем в три системы:

• система 1 (WRC-VG, RHC-MVG, Hys-KPVG);

• система 2 (WRC-KT, RHC-MKT, Hys-SKT или WRC-KT0, RHC-MKT0, Hys-SKT0 для случая ψe = 0);

• система 3 (WRC-HT, RHC-MT, Hys-SHT или WRC-HT0, RHC-MT0, Hys-SHT0 для случая ψe = 0).

В статье [9] на примере глинистой почвы 1006 Beit Netofa clay выявлено преимущество систем 2 и 3 перед системой 1 в отношении погрешностей оценок функции k(Se)/ks, когда значения экспоненциального параметра n для систем 2 и 3 оказались меньше единицы (напомним, что для системы 1 используется ограничительное условие n > 1). Но остается открытым вопрос: имеют ли преимущество системы 2 и 3 перед системой 1, когда значения экспоненциального параметра n для систем 2 и 3 оказываются больше единицы? Кроме того, предстоит ответить на следующие вопросы:

• какая из сравниваемых систем имеет наименьшую погрешность точечной аппроксимации опытных данных о главных ветвях гистерезиса водоудерживающей способности почвы;

• устраняется ли нежелательный искусственный (методический) “эффект помпы” при выполнении условия nd = nw, влияет ли это условие на погрешности оценок функции k(Se)/ks;

• какая из сравниваемых систем имеет наименьшую погрешность оценок функции k(Se)/ks;

• какая из сравниваемых систем имеет наименьшую погрешность оценок сканирующих ветвей гистерезиса Se(ψ);

• влияет ли использование аддитивного параметра ψe на погрешности точечной аппроксимации опытных данных о главных ветвях гистерезиса водоудерживающей способности почвы, а также на погрешности оценок функции k(Se)/ks и сканирующих ветвей гистерезиса Se(ψ)?

Для получения ответов на эти вопросы далее приведены результаты вычислительного эксперимента с использованием данных из авторитетного литературного источника об опесчаненном суглинке 3501 Rubicon sandy loam из каталога Муалема [24, 38]. Данная почва характеризуется следующим распределением частиц по размерам: песок (0.05–2.00 мм) – 65.2%; глина (0.002–0.05 мм) – 25.9%; ил (<0.002 мм) – 8.9%. Исследуемая почва имеет плотность сложения, равную 1.35 г ⋅ см–3; θs = 0.381 см3 ⋅ см–3; ks = 3.0 × 10–4 см ⋅ с–1.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Табл. 1 содержит параметры трех систем гидрофизических функций, идентифицированные путем точечной аппроксимации опытных данных о главных ветвях иссушения и увлажнения водоудерживающей способности почвы. Для варианта ndnw приведены параметры, рассчитанные с применением разработанных авторами компьютерных программ SoilHydrophysics-v.1.0 и SoilHysteresis-v.1.0. Для варианта nd = nw приведены: параметры системы 1 из статьи [16], а также параметры систем 2 и 3 из статьи [34].

Таблица 1.  

Параметры трех систем гидрофизических функций почвы

Система функций ndnw
θr, θs, ψe,d, ψ0,d, αd, ψe,w, ψ0,w, αw, nd, nw,
см3 · см–3 см3 · см–3 см вод. ст. см вод. ст. см вод. ст.–1 см вод. ст. см вод. ст. см вод. ст.–1 б/р б/р
1 0.1804 0.3810 –86.65 0.0115 –25.78 0.0388 6.266 3.311
2 0.1750 0.3810 –34.38 –90.50 0.0178 26.40 –34.86 0.0163 3.676 3.105
3 0.1728 0.3810 –36.37 –90.91 0.0183 24.19 –35.05 0.0169 3.668 3.117
2, ψe =0 0.1840 0.3810 –89.53 0.0112 –30.77 0.0325 6.566 2.711
3, ψe =0 0.1820 0.3810 –89.66 0.0112 –30.75 0.0325 6.794 2.731
nd = nw
1 0.1700 0.3810 –73.53 0.0136 –26.46 0.0378 3.301 3.301
2 0.1730 0.3810 –29.17 –93.03 0.0157 29.99 –35.58 0.0153 3.421 3.421
3 0.1750 0.3810 –18.05 –89.49 0.0140 26.99 –32.39 0.0168 3.679 3.679
2, ψe =0 0.1886 0.3810 –89.33 0.0119 –26.95 0.0371 4.733 4.733
3, ψe =0 0.1869 0.3810 –89.45 0.0112 –27.66 0.0361 4.724 4.724

С применением параметров из табл. 1 проведено исследование условия nd = nw для предотвращения “эффекта помпы”. В качестве примера на рис. 1 приведены результаты такого исследования с использованием модели Hys-SHT. Стрелками на рис. 1 изображены сценарии варьирования величины ψ. В зоне насыщения почвы влагой видно, что при осцилляции ψ в фиксированном диапазоне происходит последовательное асимптотическое приближение замкнутых петель, образованных сканирующими ветвями. Однако в зоне иссушения на рис. 1, A видно, что сканирующая ветвь иссушения пересекает главную ветвь увлажнения. Выполнение условия nd = nw при идентификации параметров приводит к устранению отмеченного пересечения, что подтверждает рис. 1, Б.

Рис. 1.

Осцилляция капиллярного давления влаги в фиксированных диапазонах с использованием Hys-SHT при nd ≠ nw (А) и nd = nw (Б).

Условие nd = nw, устраняющее нежелательный “эффект помпы”, можно признать обоснованным при отсутствии существенных различий между оценками k(Se)/ks, которые полученны с применением параметров, идентифицированных тремя способами: во-первых, при использовании данных только об одной главной ветви иссушения; во-вторых, при использовании данных только об одной главной ветви увлажнения; в-третьих, при использования данных об обеих главных ветвях гистерезиса Se(ψ). На рис. 2 приведены результаты вычислительного эксперимента по точечной аппроксимации главных ветвей гистерезиса Se(ψ) и оцениванию значений k(Se)/ks, выполненного тремя способами, для каждой из сравниваемых систем гидрофизических функций (при ψe ≠ 0 для систем 2 и 3).

Рис. 2.

Точечная аппроксимация данных о главных ветвях Se(ψ) и оценка значений k(Se)/ks с использованием трех систем функций и параметров nd (А), nw (Б) и при nd = nw (В): 1 – результаты расчета; 2 – данные о водоудерживающей способности, 3 – данные об относительной гидравлической проводимости почвы.

Табл. 2 содержит погрешности точечной аппроксимации опытных данных о главных ветвях гистерезиса Se(ψ), а также погрешности оценок k(Se)/ks и сканирующих ветвей гистерезиса Se(ψ).

Таблица 2.  

Корень квадратный из среднего арифметического квадратов отклонений результатов моделирования от опытных данных (RMSE − root mean square error)

Вариант вычислительного эксперимента ndnw nd = nw
системы гидрофизических функций почвы
1 2 3 ψe =0 1 2 3 ψe =0 ψe =0
2 3 2 3
Идентификация параметров по 74 точкам главных ветвей 0.0117 0.0080 0.0081 0.118 0.0111 0.0204 0.0090 0.0113 0.0185 0.0182
Оценка относительной гидравлической проводимости для 15 точек nd 0.0396 0.0184 0.0205 0.0200 0.0242
0.0262 0.0167 0.0203 0.0347 0.0420
nw
0.0484 0.0273 0.0245 0.0496 0.0436
Оценка сканирующих ветвей увлажнения для 33 точек 0.0154 0.0094 0.0101 0.0164 0.0166 0.0148 0.0057 0.0100 0.0212 0.0194
Оценка сканирующих ветвей иссушения: для 26 точек 0.0255 0.0182 0.0183 0.0267 0.0250 0.0177 0.0181 0.0135 0.0301 0.0284
Оценка всех сканирующих: для 129 точек 0.0205 0.0140 0.0143 0.0216 0.0208 0.0161 0.0127 0.0117 0.0255 0.0238

В дополнительных материалах на рис. S1 в качестве примера представлены результаты вычислительного эксперимента с моделью Hys-SКT по точечной аппроксимации главных ветвей и оцениванию сканирующих ветвей гистерезиса Se(ψ) (сплошные кривые – результаты расчета, круглые точки – данные о водоудерживающей способности почвы). Кроме того, дополнительные материалы содержат результаты выявления достоверных различий между погрешностями сравниваемых систем в отношении: (i) точечной аппроксимации опытных данных о главных ветвях гистерезиса Se(ψ) (табл. 3 ); (ii) оценок k(Se)/ks (табл. 4 ); (iii) оценок сканирующих ветвей гистерезиса Se(ψ) (табл. 5 ).

ЗАКЛЮЧЕНИЕ

Представлено описание зависимостей Se(ψ) и k(Se)/ks в виде трех систем функций. В каждой системе используется соответствующий набор общих параметров. Идентифицированы параметры этих систем путем точечной аппроксимации данных о главных ветвях иссушения и увлажнения гистерезиса водоудерживающей способности почвы. Использованы литературные данные об опесчаненном суглинке 3501 Rubicon sandy loam [24, 38]. Выполнено построение петель гистерезиса при осцилляции ψ в заданном диапазоне значений для вариантов ndnw и nd = nw. Осуществлено оценивание значений функции k(Se)/ks и сканирующих ветвей гистерезиса Se(ψ) с использованием трех систем функций. Для выявления различий между погрешностями оценок применен критерий Вильямса–Клута.

Из анализа результатов исследования вытекают следующие выводы:

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

• в варианте ndnw система 1 достоверно уступает системам 2 и 3 в случае использования параметра ψe, при этом между системами 2 и 3 нет существенных различий; система 3 достоверно превосходит системы 1 и 2, если параметр ψe не применяется (при ψe = 0), при этом между системами 1 и 2 нет существенных различий; использование параметра ψe достоверно уменьшает погрешности;

• в варианте nd = nw система 1 достоверно уступает системам 2 и 3 в случае использования параметра ψe, при этом система 2 превосходит систему 3; сравниваемые системы являются неразличимыми, если параметр ψe не применяется (при ψe = 0); использование параметра ψe существенно уменьшает погрешность.

2. В отношении погрешностей оценок относительной гидравлической проводимости почвы:

• между погрешностями оценок, полученных с использованием параметров, идентифицированных по данным только о главной ветви иссушения, и погрешностями оценок, полученных с использованием параметров, идентифицированных по данным только о главной ветви увлажнения, не выявлены достоверные различия ни для одной из систем;

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

• для систем 1 и 2 оценки, полученные с применением параметров, идентифицированных по данным об обеих главных ветвях, являются более точными по сравнению с оценками, полученными с применением параметров, идентифицированных по данным только о главной ветви увлажнения, если используется параметр ψe; однако между погрешностями полученных таким образом оценок нет достоверных различий, если параметр ψe не применяется (при ψe = 0);

• в варианте nd = nw система 1 достоверно уступает системе 2 в случае использования параметра ψe, при этом между системами 1 и 3, а также между системами 2 и 3 нет существенных различий; система 2 достоверно превосходит систему 3, если параметр ψe не применяется (при ψe = 0), при этом между системами 1 и 2 и между системами 1 и 3 нет существенных различий; использование параметра ψe существенно уменьшает погрешности системы 3;

• отсутствие достоверных различий между погрешностями оценок, полученных с использованием параметров, которые идентифицированы тремя способами, свидетельствует о том, что зависимость k(Se)/ks не является гистерезисной в отличие от сложной функции k(Se(ψ))/ks, в которой зависимость Se(ψ) является гистерезисной;

• погрешности оценок при значениях экспоненциального параметра n систем 2 и 3, превышающих единицу, наряду с погрешностями ранее полученных оценок при значениях этого параметра n систем 2 и 3, которые оказались меньше единицы [9], позволяют утверждать о достоверном преимуществе систем 2 и 3 над системой 1.

3. В отношении погрешностей оценок сканирующих ветвей гистерезиса Se(ψ):

• в варианте ndnw система 1 достоверно уступает системам 2 и 3 в случае использования параметра ψe, при этом система 2 превосходит систему 3; между системами 1 и 3 нет существенных различий, если параметр ψe не применяется (при ψe = 0), при этом системы 1 и 3 точнее системы 2; использование параметра ψe достоверно уменьшает погрешности;

• в варианте nd = nw система 1 достоверно уступает системам 2 и 3 в случае использования параметра ψe, при этом система 3 превосходит систему 2; система 1 точнее систем 2 и 3, если параметр ψe не применяется (при ψe = 0), при этом система 3 превосходит систему 2; использование параметра ψe достоверно уменьшает погрешности;

• нежелательный искусственный (методический) “эффект помпы” устраняется при выполнении условия nd = nw, при этом увеличения погрешностей оценок значений функции k(Se)/ks не отмечается.

4. Мультипликативный параметр α системы 1 не является величиной, обратной давлению барботирования. Это давление описывается аддитивным параметром ψe систем 2 и 3. Использование параметра ψe: во-первых, достоверно уменьшает погрешности точечной аппроксимации опытных данных о главных ветвях гистерезиса Se(ψ), а также погрешности оценок значений функции k(Se)/ks и сканирующих ветвей гистерезиса Se(ψ); во-вторых, позволяет описать феномен гистерезиса во всем диапазоне значений ψ, включая положительную область, в которой, как правило, происходит вытеснение воздуха, защемленного в тупиковых порах, на завершающем этапе насыщения почвы водой и смыкаются главные ветви гистерезиса Se(ψ).

5. Оценивание значений функции k(Se)/ks системы 3 при ψe = 0 с применением параметров модели, которая предложена в статье [14], идентифицируемых путем точечной аппроксимации данных θ(ψ), представляет собой математически корректное решение задачи Ван Генухтена в ее исходной постановке [39]. Достоинства систем 2 и 3 позволяют рекомендовать эти системы к применению для моделирования гидрофизических свойств почвы и решения задач ирригационного земледелия. В системе 3 зависимости Se(ψ) и k(Se)/ks сформулированы в достаточно простом виде с использованием элементарных математических функций; при этом во многих случаях погрешности оценок, полученных с применением систем 2 и 3, являются неразличимыми. Поэтому авторы данной статьи отдают предпочтение системе 3 (WRC-HT, RHC-MT, Hys-SHT при nd = nw) с физически интерпретированным аддитивным параметром ψe.

Список литературы

  1. Воронин А.Д. Основы физики почв. М.: Изд-во Моск. ун-та, 1986. 244 с.

  2. Глобус А.М. Почвенно-гидрофизическое обеспечение агроэкологических математических моделей. Л.: Гидрометеоиздат, 1987. 427 с.

  3. Глобус А.М. Экспериментальная гидрофизика почв. Л.: Гидрометеоиздат, 1969. 356 с.

  4. Кобзарь А.И. Прикладная математическая статистика. Для инженеров и научных работников. М.: Физматлит, 2006. 816 с.

  5. Нерпин С.В., Чудновский А.Ф. Физика почв. М.: Наука, 1967. 584 с.

  6. Семенова Н.Н., Терлеев В.В., Сухорученко Г.И., Орлова Е.Е., Орлова Н.Е. Об одном методе численного решения системы параболических уравнений // Вестник Санкт-Петербургского гос. ун-та. Сер. 1. Математика. 2016. Т. 3(61). Вып. 2. С. 230–240. https://doi.org/10.21638/11701/spbu01.2016.207

  7. Смагин А.В. К термодинамической теории водоудерживающей способности и дисперсности почв // Почвоведение. 2018. № 7. С. 836‑851. https://doi.org/10.1134/S0032180X18070092

  8. Судницын И.И. Новые методы оценки водно-физических свойств почв и влагообеспеченности леса. М.: Наука, 1966. 95 с.

  9. Терлеев В.В., Миршель В., Баденко В.Л., Гусева И.Ю. Усовершенствованный метод Муалема-Ван Генухтена и его верификация на примере глинистой почвы Бейт Нетофа // Почвоведение. 2017. № 4. С. 457–467. https://doi.org/10.7868/S0032180X1704013X

  10. Шеин Е.В. Математические физически обоснованные модели в почвоведении: история развития, современное состояние, проблемы и перспективы (аналитический обзор) // Почвоведение. 2015. № 7. С. 816–823. https://doi.org/10.7868/S0032180X15070096

  11. Шеин Е.В., Позднякова А.Д., Шваров А.П., Ильин Л.И., Сорокина Н.В. Гидрофизические свойства высокозольных низинных торфяных почв // Почвоведение. 2018. № 10. С. 1259–1264. https://doi.org/10.1134/S0032180X18100118

  12. Faybishenko B.A. Hydraulic behavior of quasi-saturated soils in the presence of entrapped air: Laboratory experiments // Water Resources Research. 1995. V. 31(10). P. 2421–2435. https://doi.org/10.1029/95WR01654

  13. Gonc’alves Roger D., Teramoto Elias H., Engelbrecht Bruno Z., Alfaro Soto Miguel A., Chang Hung K., Van Genuchten Martinus Th. Quasi-Saturated Layer: Implications for Estimating Recharge and Groundwater Modeling // Ground Water. 2020. V. 58(3). P. 432–440. https://doi.org/10.1111/gwat.12916

  14. Haverkamp R., Vauclin M., Touma J., Wierenga P.J., Vachaud G. A comparison of numerical simulation model for one-dimensional infiltration // Soil Sci. Soc. Am. J. 1977. V. 41. P. 285–294.

  15. Konyai S., Sriboonlue V., Trelo-Ges V. The Effect of Air Entry Values on Hysteresis of Water Retention Curve in Saline Soil // Am. J. Environmental Sci. 2009. V. 5(3). P. 341–345. https://doi.org/10.3844/ajessp.2009.341.345

  16. Kool J.B., Parker J.C. Development and evaluation of closed-form expressions for hysteretic soil hydraulic properties // Water Resources Research. 1987. V. 23(1). P. 105–114.

  17. Kosugi K. General model for unsaturated hydraulic conductivity for soil with lognormal pore-size distribution // Soil Sci. Soc. Am. J. 1999. V. 63. P. 270–277.

  18. Kosugi K. Lognormal distribution model for unsaturated soil hydraulic properties // Water Resour. Res. 1996. V. 32. P. 2697–2703.

  19. Kosugi K. Three-parameter lognormal distribution model for soil water retention // Water Resour. Res. 1994. V. 30. P. 891–901.

  20. Kosugi K., Hopmans J.W. Scaling water retention curves for soils with lognormal pore-size distribution // Soil Sci. Soc. Am. J. 1998. V. 62. P. 1496–1505.

  21. Levenberg K. A Method for the Solution of Certain Non-Linear Problems in Least Squares // Quarterly Appl. Math. 1944. V. 2. P. 164–168.

  22. Mady A.Y., Shein E. Modelling and validation hysteresis in soil water retention curve using tomography of pore structure // International J. Water. 2018. V. 12(4). P. 370–381. https://doi.org/10.1504/IJW.2018.095403

  23. Marquardt D.W. An algorithm for least-square estimation on non-linear parameters // J. Soc. Ind. Appl. Math. 1963. V. 11. P. 431–441.

  24. Mualem Y. A catalogue of the hydraulic properties of unsaturated soils. Research Project 442. Technion, Israel Institute of Technology, Haifa, Israel, 1976. 100 p.

  25. Mualem Y. A conceptual model of hysteresis // Water Resour. Res. 1974. V. 10(3). P. 514–520.

  26. Mualem Y. A new model for predicting hydraulic conductivity of unsaturated porous media // Water Resour. Res. 1976. V. 12. P. 513–522.

  27. Mualem Y. Modified approach to capillary hysteresis based on a similarity hypothesis // Water Resour. Res. 1973. V. 9(5). P. 1324–1331.

  28. Mualem Y., Morel-Seytoux H.J. Analysis of a capillary hysteresis model based on a one-variable distribution function // Water Resour. Res. 1978. V. 4(4). P. 605–610.

  29. Nikonorov A., Terleev V., Pavlov S., Togo I., Volkova Yu., Makarova T., Garmanov V., Shishov D., Mirschel W. Applying the model of soil hydrophysical properties for arrangements of temporary enclosing structures // Procedia Engineering. 2016. V. 165. P. 1741–1747. https://doi.org/10.1016/j.proeng.2016.11.917

  30. Poluektov R.A., Terleev V.V. Crop simulation model of the second and the third productivity levels // Modelling water and nutrient dynamics in soil-crop systems, Springer, Dordrecht, The Netherlands. 2007. P. 75–89. https://doi.org/10.1007/978-1-4020-4479-3_7

  31. Rawitz E. The influence of a number of environmental factors on the availability of soil moisture to plants (in Hebrew). Ph. D. thesis, Hebrew Univ., Rehovot, Israel. 1965.

  32. Richards L.A. Capillary conduction of liquids through porous mediums // J. Appl. Physics. 1931. V. 1(5). P. 318–333.

  33. Scott P.S., Farquhar G.J., Kouwen N. Hysteretic effects on net infiltration. Proceeding of National Conference on Advances in Infiltration, Publication 11–83, American Society of Agricultural Engineers, St. Joseph, Michigan, 1983. P. 163–170.

  34. Terleev V., Nikonorov A., Ginevsky R., Lazarev V., Topaj A., Dunaieva I., Terleeva A. Estimation of precise irrigation rates taking into account the hysteresis of soil water-retention capacity // IOP Conf. Series: Earth Environ. Sci. 2019. V. 403. Article Number 012239. https://doi.org/10.1088/1755-1315/403/1/012239

  35. Terleev V., Nikonorov A., Togo I., Volkova Yu., Garmanov V., Shishov D., Pavlova V., Semenova N., Mirschel W. Modelling the hysteretic water retention capacity of soil for reclamation research as a part of underground development // Procedia Engineering. 2016. V. 165. P. 1776–1783. https://doi.org/10.1016/j.proeng.2016.11.922

  36. Terleev V., Petrovskaia E., Nikonorov A., Badenko V., Volkova Yu., Pavlov S., Semenova N., Moiseev K., Topaj A., Mirschel W. Mathematical modeling the hydrological properties of soil for practical use in the land ecological management // MATEC Web of Conf. 2016. V. 73. Article Number 03001. https://doi.org/10.1051/matecconf/20167303001

  37. Terleev V., Petrovskaia E., Sokolova N., Dashkina A., Guseva I., Badenko V., Volkova Yu., Skvortsova O., Nikonova O., Pavlov S., Nikonorov A., Garmanov V., Mirschel W. Mathematical modeling of hydrophysical properties of soils in engineering and reclamation surveys // MATEC Web of Conferences. 2016. V. 53. Article Number 01013. https://doi.org/10.1051/matecconf/20165301013

  38. Topp G.C. Soil-water hysteresis measured in a sandy loam and compared with the hysteretic domain model // Soil Sci. Soc. Am. J. 1969. V. 33(5). P. 645–651. https://doi.org/10.2136/sssaj1969.03615995003300050011x

  39. Van Genuchten M.Th. A closed form equation for predicting the hydraulic conductivity of unsaturated soils // Soil Sci. Soc. Am. J. 1980. V. 44. P. 892–989.

  40. Vereecken H., Weynants M., Javaux M., Pachepsky Y., Schaap M.G., Van Genuchten M.Th. Using pedotransfer functions to estimate the Van Genuchten-Mualem soil hydraulic properties: A review // Vadose Zone J. 2010. V. 9. P. 795–820.

Дополнительные материалы

скачать ESM.docx
Рис. S1a. - Рис. S1d.
 
Таблица S1. - Таблица S3.