一种高阶保界的有限差分多分辨率WENO方法，用于高维理查兹方程

《Advances in Water Resources》：A high-order bound-preserving finite difference multi-resolution WENO method for high-dimensional Richards equation

　　吴新宇|程娟|舒志旺

摘要
在本文中，我们提出了一种新的高阶保守型有限差分多分辨率加权本质无振荡（WENO）方法，用于求解非线性混合形式的Richards方程，以模拟非饱和渗透流。时间离散化采用三阶强稳定性保持（SSP）龙格-库塔（RK）方法，具有较低的计算复杂度。作为一个关键创新点，我们引入了一种参数化的通量限制技术，以在保持计算稳定性和鲁棒性的同时，将数值解保持在物理范围内。这种新方法在有效抑制 steep 润湿前沿附近的非物理振荡的同时，实现了真正的高阶精度。该算法框架易于构建，在三维（3D）空间中具有高效性。通过针对多个维度中的各种水力模型进行了广泛的数值实验，证明了所提方法的良好性能和实用性。

引言
在非饱和土壤中预测和模拟渗透流是环境工程、农业灌溉、石油回收和水资源等多个科学和工程领域的重要课题。用于描述非饱和水运动的经典模型是Richards方程，该方程由广义达西定律和质量守恒方程结合而成。由于复杂的地质介质导致的高度非线性，解析解几乎无法获得。因此，迫切需要开发具有先进数值方法的高效且精确的模拟器。
Richards方程通常有三种表达形式，其中一种是混合形式：
?θ?t???(K(h)?h)+?K(h)?z=0
其中 Hz 是含水量，h 是水头压力，K 是水力传导率，z 是向下取正的垂直坐标，t 是时间。通过特定的流体含量函数 C(h)=?θ/?h，基于 h 的形式表示为
C(h)?h?t???(K(h)?h)+?K(h)?z=0。
定义流体扩散率函数 D(θ)=K(θ)?h/?θ=K(θ)/C(θ) 可得到基于 θ 的形式：
?θ?t???(D(θ)?θ)+?K(θ)?z=0。
与基于 h 的形式相比，基于 θ 的形式在建模饱和-非饱和区域时更加灵活，但表现出强烈的非线性。由于它不满足质量守恒的形式，因此设计保守的数值方法极为困难。因此，许多当前的努力都集中在减少质量守恒误差（Celia 等人，1990年；Milly，1985年；Sadegh Zadeh，2011年），以确保数值结果的可靠性。相比之下，基于 θ 的形式具有适度的非线性和完美的守恒公式。然而，在饱和区域由于扩散率变为无穷大而面临计算挑战，并且对于复合的、分层的或真实的非均质多孔介质则不适用。混合形式通过满足质量守恒、避免数值奇异性并适用于真实的非饱和区域，解决了上述问题。因此，本文的重点是开发一种用于混合形式Richards方程（1.1）的高分辨率保守数值方法。值得注意的是，所提出的方法也适用于保守的基于 θ 的形式（1.3），但需要对饱和引起的奇异性给予额外关注。

Richards方程的强非线性源于复杂的土壤本构关系（即 K(h) 和 h(θ) 之间的依赖性），使其成为所有水科学领域中最具有挑战性的方程之一（Farthing 和 Ogden，2017年）。只有通过做出一些简单而不切实际的假设（Basha，1999年；Gardner，1958年；Hayek，2015年；Hayek，2016年；Tracy，2006年；Tracy，2007年），才能获得解析解，因此数值方法成为解决该方程的唯一可行方法。在过去的四十年中，科学家们提出了许多模拟非饱和流动问题的方法，特别是有限元方法（FEM）（Abriola 和 Lang，1990年；Celia 等人，1987年；Neuman，1973年；Simunek 等人，2006年）、有限体积方法（FVM）（Horgue 等人，2015a年；Svyatskiy 和 Lipnikov，2017年；Manzini 和 Ferraris，2004年；Su 等人，2022年）以及有限差分方法（FDM）（Dogan 和 Motz，2005年；Celia 等人，1990年；Haverkamp 等人，1977年；Romano 等人，1998年）。还有一些其他方法，如欧拉-拉格朗日方法（Huang 等人，1994年）和混合有限元方法（Farthing 等人，2003年；Bause 和 Knaber，2004年；Radu 等人，2004年），也被应用于求解Richards方程，实现了几乎无振荡的行为和质量守恒等优点。此外，还投入了大量精力设计高效的自适应算法（Miller 等人，2006年；Kavetski 等人，2001年）和研究主要变量切换技术（Diersch 和 Perrochet，1999年；Wu 和 Forsyth，2001年；Sadegh Zadeh，2011年），以适应各种区域特征。
Richards方程的解通常表现出陡峭的梯度和移动的润湿前沿。此外，当对流项占主导时，还会出现双曲行为，表现为几乎不连续的尖锐前沿和上风特性。传统方法经常受到非物理振荡、数值扩散和质量守恒误差的困扰。此外，大多数当前的工作仅具有二阶精度，对于Richards方程的高精度空间离散化方法的研究很少。Zambra 等人（2012年）通过使用有限体积WENO方法进行空间离散化和空间及时间的任意导数（ADER）框架进行时间离散化，实现了三阶精度。Boada 等人（2020年）采用高阶仿射有限差分算子方法解决了混合形式Richards方程的数值问题，达到了四阶精度。但这些方法仅在一维（1D）空间中实现，不适用于复杂的高维模拟。此外，Richards方程的解满足严格的界限保持（BP）性质，即含水量 θ∈[θr,θs] 和非饱和水头 h≤0。由于水力性质的强非线性，这一约束对于数值稳定性至关重要。一旦解超出物理限制，模拟将失败。不幸的是，我们不知道有任何研究关注Richards方程的界限保持技术。虽然直接截断可以强制执行解的边界，但会破坏质量守恒和高精度，使得数值结果缺乏可靠性和可信度。
为了准确捕捉润湿前沿的动态并实现真实的3D模拟，需要一种高阶、保持界限、高效且保守的数值方法。WENO方法是一类具有高精度的高阶方案，适用于求解双曲守恒律和对流占主导的对流-扩散偏微分方程（PDEs），具有诸如守恒性、在平滑区域的 high 阶精度以及在前沿附近无振荡等理想属性。WENO方法基于本质无振荡（ENO）方法（Harten 等人，1987年；Shu 和 Osher，1988年；Shu 和 Osher，1989年），不同之处在于使用加权模板而不是仅选择最平滑的模板来重构解。第一个WENO方案是由Liu等人（1994年）设计的三阶有限体积方案。后来，Jiang和Shu（Jiang 和 Shu，1996年）通过优化平滑度指标和非线性权重的计算，提出了一个更高效和通用的WENO方案，从而提高了鲁棒性和更高阶精度。最近在Zhu和Shu（2018年）的工作中，他们基于多分辨率思想设计了一类新的WENO方案，在这种方案中，模板大小不等，线性权重可以是任意正数，仅要求它们的总和等于1。多分辨率WENO方法可以轻松实现任意高的精度和更高维度，已广泛应用于求解双曲问题（Tsybulnik等人，2023年；Wang和Zhu，2020年）和非线性退化抛物方程（Jiang，2021年；Xu和Zhang，2025年）。为了在保持高精度的同时保留数值解的物理界限，Jiang和Xu（Jiang 和 Xu，2013年）为高阶有限差分龙格-库塔（RK）WENO方法设计了一种参数化的最大原理保持（MPP）通量限制器。这种方法直接修改高阶数值通量，并适用于保守形式的多分辨率WENO方案。
在本文中，我们开发了一种新的高阶有限差分多分辨率WENO方案来求解混合形式的Richards方程。有限差分框架允许逐维度进行计算（而不是维度分割，因为维度分割会降低精度），从而便于编码并提高效率，同时保持高阶精度，特别是在3D空间中。采用三阶SSP RK方法进行时间离散化，以保持高阶时间精度和强稳定性，同时避免非线性迭代。为了增强算法的稳定性和鲁棒性，结合了BP通量限制技术，以确保数值解保持在物理上现实的范围内。所提出的算法流程可以总结如下：
(1) 为混合（或基于 θ 的）形式的Richards方程构建一个半离散保守有限差分方案。通过对流项使用多分辨率WENO重构（Zhu和Shu，2018年），获得高阶数值通量，不同之处在于使用Roe速度而不是通量分割来确保正确的上风偏置和稳定性。对于扩散项，依次执行两次具有中心模板的多分辨率WENO重构。
(2) 使用三阶SSP RK方法对时间导数进行离散化。
(3) 在RK离散化的最后阶段应用通量限制器来修改高阶数值通量，从而在不牺牲精度的情况下强制执行解的界限。
我们对三种不同的水力模型——Van Genuchten–Mualem、Brooks-Corey和Gardner模型进行了数值模拟，证明了所提方法的广泛适用性和普适性。展示了四阶和六阶精度的数值结果，并且使用更宽的模板可以获得任意更高的精度。
本文的其余部分组织如下。第2节描述了三种常用的Richards方程水力模型，并详细介绍了所提出的有限差分多分辨率WENO方案。第3节详细介绍了参数化的界限保持通量限制技术的细节。第4节进行了一组全面的数值实验，以证明该方法的有效性、精度、稳定性和一致性。最后，在第5节得出结论。

水力模型
求解Richards方程需要了解水力传导率K和含水量θ与水头h之间的功能关系。这些功能关系被称为多孔介质的水力性质。已经实现了几种常用的水力模型，其中之一是Van Genuchten–Mualem模型（Vangenuchten，1980年）：
θe=θ?θrθs?θr=1+(α|h|)n?m，
K=Ksθel1?1?θe1/mm2，
其中θr和θs分别是残余含水量和饱和含水量，θe是...

参数化的界限保持通量限制器
1D Richards方程（2.4）的一个关键性质是解θ(z,t)满足最大原理，因为K(h)≥0，这意味着解被初始条件θ(z,0)=θ0(z)和空间边界所限制。具体来说，对于周期性或均匀的诺伊曼（无通量）边界条件，θ(z,t)∈[θm,θM]，t>0，
其中θM=maxzθ0(z)，θm=minzθ0(z)。
某些施加非零通量（非均匀诺伊曼条件?h?z≠0）或固定值的边界条件...

数值实验
在本节中，我们进行了数值测试，以证明所提出的BP有限差分多分辨率WENO方法对Richards方程的性能。对于下面介绍的WENO-2k方案，扩散项以2k阶精度离散化，以降低其较高的计算成本，而对流项采用更高一阶（2k+1）近似，以降低其边际成本。时间步长取为Δt=0.4c/Δz+b/Δz2（对于1D问题），Δt=0.4c/Δz+b/Δx2+b/Δz2（对于...

结论
在本文中，我们开发了一种新颖且高效的高阶界限保持有限差分多分辨率WENO方法，用于求解保守的Richards方程，包括混合形式和基于θ的形式。凭借其自然的守恒性、简单性和高阶精度优势，该方法有效地处理了Richards方程的强非线性，并准确捕捉了陡峭的润湿前沿，没有非物理振荡。我们工作的一个关键创新点...

作者贡献声明
吴新宇：撰写 - 审稿与编辑、撰写 - 原始草案、软件、方法论、形式分析、概念化。
程娟：撰写 - 审稿与编辑、撰写 - 原始草案、监督、方法论、资金获取、形式分析、概念化。
舒志旺：撰写 - 审稿与编辑、撰写 - 原始草案、监督、方法论、资金获取、形式分析、概念化。

未引用的参考文献
表D.2

利益冲突声明
作者声明他们没有已知的财务利益或个人关系可能影响本文报告的工作。