摘要

变饱和多孔介质中水流的Richards方程解对水资源评估和环境管理越来越有用。

除了解的准确性，对所需的计算工作量也颇为关注，特别是在涉及高度非线性土壤水力特性和干燥初始条件时。

本文评估了在使用修正的Picard迭代方法解混合形式Richards方程时不同收敛准则的性能。

结果以计算机处理（CPU）时间和迭代次数进行比较。

通过使用含水量的Taylor级数展开导出的新的非线性收敛准则被实施到混合形式的数值算法中。

新准则的计算效率针对不同的土壤类型、边界条件、初始条件和分层土壤进行了评估。

尽管三个准则在计算的含水量、压力水头和水通量分布方面产生了几乎相同的结果，且质量平衡误差微乎其微，但所需的CPU时间明显不同。

总体而言，新的非线性收敛准则在计算效率上远远优于其他两个准则。

在高度非线性流问题中，其他两个收敛准则失败时，新的准则也更为健壮（即解保持收敛）。

本研究结果表明，当新的收敛准则应用于混合形式Richards方程的修正Picard解时，可以产生一种非常高效和准确的方法，用于模拟土壤中的变饱和水流。

1. 引言

在环境研究中，数值模型是评估非饱和带释放化学物质可能导致地下水污染风险的重要工具

地质材料的极端变异性和复杂性、干燥的初始条件以及不断变化的边界条件可能使得在可接受的准确性和计算工作量范围内解决流动和传输问题变得困难。

当前大多数可用的数值方法在准确性或计算效率之间做出了妥协。

例如，对于相对粗纹理的土壤的渗透通常难以建模，因为其具有高度非线性的水力特性。

干燥的初始条件和高度非线性的水力特性通常需要使用非常细致的空间和时间离散化来避免数值不稳定性。

这些条件导致数值算法变得CPU密集型，特别是当需要模拟长期和/或多维问题时。

因此，追求最佳数值算法的努力应该涉及优化方案的准确性和健壮性，以及最小化所需的计算时间。

 

本文的目标是提高Celia等人（1990）的混合形式算法用于解决变饱和流问题的计算效率。

通过引入新的非线性收敛准则，改进了算法的收敛速度。

所提出的准则将通过与其他两个广泛使用的收敛准则进行比较来评估其性能。

2. 背景

在等温条件下，可变饱和刚性土壤中的一维垂直水流通常用Richards方程（Richards，1931）描述，该方程可以用压力头或水含量的形式表示，即

或

其中，为压力头L，

为体积含水量（），

为时间（），

表示距离土壤表面向下的垂直距离[L]，

为比水含量容积（），

为水力导度（），

为土壤水非饱和扩散率（），

为源/汇项（）。

直到最近，大多数数值研究使用Richards方程的压力水头或含水量形式来描述可变饱和土壤中的流动（例如，Davis和Neuman，1983；Huyakorn等，1983，1989；Hills等，1989；Kool和Van Genuchten，1991；Kirkland等，1992）。

使用方程（1）的基于含水量的方案。

(2) 可以被写成保守质量的形式，因此在大多数情况下，无论时间步长和网格间距如何，都应在计算域内保持质量（Hills等人，1989）。

Huyakorn和Pinder（1983）表明，对于最初是干燥均质土壤，使用这样的方案是有优势的。

基于含水量的公式的一个局限性是，这种形式不能用于描述饱和带的流动，而且分层土壤中的流动也不容易模拟。

此外，即使这种公式在流动系统内部准确保持质量，基于含水量的算法在边界处可能出现质量平衡错误。

基于压力头的公式（1）被认为对涉及分层或空间异质土壤的实际问题以及变饱和流问题更有用。

不幸的是，使用基于压力头的公式模拟干燥和/或高度非线性土壤的入渗通常在保持质量方面面临困难。

有关基于压力头和基于水含量的Richards方程形式的相对优势和劣势的更详细讨论已由Hills等人（1989）提供。

一些研究人员探索了解决基于压力头或基于水含量形式的Richards方程的替代数值技术。

这些研究的目标包括获得更稳定的数值算法，加速计算，最小化质量平衡误差，并针对不同的土壤类型或初始和边界条件实现更精确的解。

例如，Milly（1985）提出了一种质量保守的解决方案，其中使用了有效元素土壤水容量项。

这种方法，结合质量集中法（Neuman，1973），有效地确保了基于压力头的方程的整体质量平衡。

Rathfelder和Abriola（1994）通过扩展和离散化土壤水容量开发了类似基于压力头方程的质量保守解决方案。

其他人使用了Kirchhoff类型的转换（例如，Ross和Bristow，1990）或替代函数（Ross，1990；Pan和Wierenga，1995）以简化非线性流的描述。

在不同的方法中，Gottardi和Venutelli（1992）使用了一种移动有限元方法，其中网格点沿润湿前沿移动，从而允许较少的节点而不损害数值精度。

然而，移动网格法在应用于分层系统或用于时变边界条件时存在一些局限性。该方法也被发现不如传统的固定网格公式具有较差的质量保守性。

El-Kadi和Ling（1993）提出了Peclet和Courant数标准，用于空间和时间离散化，以描述解决Richards方程的数值方案的准确性和效率。

通过引入源项，Hills等人（1989）成功地用基于含水量的算法解决了涉及一维水流进入分层土壤的问题。

Kirkland等人（1992）随后开发了一个基于含水量的算法，通过变量转换模拟了二维变饱和流。

最近，Huang等人（1994）提出了一种基于特征的粒子跟踪技术，用于解决基于压力头的Richards方程，用于高度非线性的入渗问题。

察觉到现有基于压力头和基于水含量的Richards方程解决方案的缺点，许多人尝试将这两种方法的优势结合起来。

Richards方程的混合形式被认为保持了基于水含量方程固有的质量守恒特性，同时提供了以压力头h表示的解。

Richards方程的混合形式表达为：

Brutsaert（1971）是最早使用混合形式Richards方程来解决饱和-非饱和流问题的人之一。

他将混合形式方程的有限差分逼近与牛顿迭代方案结合起来，以有效处理陡峻的润湿前缘。

Allen和Murphy（1985，1986）在其共轭有限元算法中使用了Richards方程的混合形式，并采用了“拟牛顿”迭代方法。

最近，Celia等人（1987）和Celia等人（1990）提出了一种质量保守的数值方案，使用“修正的Picard迭代”解混合形式Richards方程（3）。

由于完美的质量平衡，修正的Picard迭代技术被证明是对早期Picard方法的重大改进。

修正的Picard迭代方法在模拟具有陡峭润湿前缘的非饱和流方面也表现出很大的潜力（Celia等人，1990；Celia和Bining，1992）。

Ray和Mohanty（1992）随后重新审视了混合形式算法，并通过几个数值实验展示了它相对于基于压力头的方案的优势。

与基于压力水头和基于含水量的算法类似，Celia等人（1990）的修正Picard迭代方案基于Richards方程（3）的全隐式（向后Euler）时间逼近，具体如下：