饱和与非饱和多孔介质水过滤建模的确定性与随机性方法：一篇综述

《Hydrology》：Deterministic and Stochastic Approaches to Modeling Water Filtration Through Saturated and Unsaturated Porous Media: A Review

　　

水流通过饱和与非饱和多孔介质的过程是地下水评估、污染物运移、农业排水及坝体渗流分析的基础，然而现有综述通常仅涉及该链条中的单一环节：单一控制方程、单一流态或单一建模范式。本综述将这些环节串联起来。本文重新审视了过滤的水动力学基础，从Darcy的原始实验到稳态与瞬态饱和流方程，并通过Richards方程及土壤水分保持模型将处理范围扩展至变饱和条件。随后，本文直接比较了确定性与随机性建模范式，并通过一项原创数值案例研究使该比较具体化：针对砂滤柱中的积水入渗求解Richards方程，先进行确定性求解，再对空间相关的水力传导度场进行经过收敛性检验的Monte Carlo系综（最多500次实现）求解，结果表明，仅现实中的非均质性即可使预测的湿润锋到达时间产生超过三倍的离散度。本文还综述了数值求解器、高性能计算及机器学习替代模型的最新进展，以展示两种范式如何融合为混合的、物理信息驱动的框架，并明确将本综述与2017至2025年间发表的七项相关研究进行了对比定位。本文指出了非均质性表征及混合模型验证方面持续存在的挑战，并概述了未来研究方向。

主体内容总结

1. 引言

本文指出，液体在饱和多孔介质（土壤、含水层及工程滤床）中的运移控制着地下水补给、污染物运移、土壤—含水层途径的饮用水处理、农业排水以及坝体和堤防的渗流稳定性等一系列过程。自Henry Darcy于1856年对砂柱水流开展实验以来，该水动力学规律已被研究近两个世纪，且因地下水流精确预测对水资源管理、环境保护和岩土安全至关重要而持续吸引着活跃研究。饱和区中所有孔隙空间均被水占据，经典水动力过滤定律直接适用；其上的非饱和（包气带）带则在重力和毛细力共同作用下输水，需要扩展的理论处理。确定性方法从给定的介质参数、源项和边界条件求解控制偏微分方程以获得单一可复现的解场；随机方法则将介质的水力性质——最重要的是水力传导度——视为空间相关的随机变量，并将该不确定性传播为可能的流态系综。由于天然土壤和含水层很少呈现纯确定性模型所假设的均质、充分表征的条件，实际中最稳健的流模型通常结合两种传统。

2. 数学基础

本节以压缩形式总结了饱和与非饱和多孔介质流的经典控制方程，为全文的确定性—随机性比较和数值案例研究提供物理基础。值得注意的是，并非每种多孔、裂隙或混合介质都需要完全独立的数学处理：近期已发展出统一的积分微分及非局部（近场动力学和分数阶导数）公式，可在单一数学框架内描述跨越裂隙、奇异性和强非均质夹杂的流动。

2.1 饱和流：Darcy定律与控制方程

过滤的数学描述由法国工程师Henry Darcy基于垂直砂柱过滤实验建立。Darcy定律以过滤（Darcy）速度表达：v = Q/A = ?K dh/dL，其中Q为液体流量[m3/s]，A为流动截面积[m2]，K为饱和介质的过滤系数即水力传导度[m/s]，dh/dL为水力梯度[–]，h为水头[m]，L为液柱长度[m]。水力传导度同时取决于过滤介质结构和渗流液体性质（尤其黏度）。Darcy定律严格而言是唯象近似，仅在理想化平行毛细管束情形下精确，且在有限条件范围内成立：在高Reynolds数的岩溶或强裂隙岩石中失效，在极低渗透性黏土和页岩中低于阈值水力梯度时则进入所谓前Darcy流区域。对各向异性介质，Darcy定律以微分形式表达为v_x = ?K_x?h/?x等。将Darcy定律与单元土体的质量守恒结合，依次得到饱和流控制方程：稳态条件下为?(ρv_x)/?x + ?(ρv_y)/?y + ?(ρv_z)/?z = 0；对不可压缩常密度流体简化为?v_x/?x + ?v_y/?y + ?v_z/?z = 0。结合Darcy定律得到饱和多孔介质流动控制方程?/?x(K_x?h/?x) + ?/?y(K_y?h/?y) + ?/?z(K_z?h/?z) = 0；对各向同性均质介质简化为Laplace方程。非稳态条件下，质量守恒要求进入单元体的流体质量等于其内部储存质量的时间变化率，结合储水关系得到瞬态流动方程?/?x(K_x?h/?x) + ?/?y(K_y?h/?y) + ?/?z(K_z?h/?z) = S_s?h/?t，其中S_s = ρg(α + φβ)为比储水系数。求解该方程需预先知道K、S_s等参数，而这些参数在整个流动域内很少能确定已知——这是引入随机处理的核心动机。

2.2 非饱和流：Richards方程

在地下水位以上非饱和（包气带）带中，含水量θ和水力传导度K均随局部基质吸力水头ψ变化。将Darcy定律推广至K(ψ)或K(θ)并与质量守恒结合，得到Richards方程——1931年首次提出，至今仍是变饱和流的标准物理描述；其一维垂直混合形式为?θ/?t = ?/?z[K(θ)(?ψ/?z + 1)]，其中θ为体积含水量[–]，z为高程[m]，t为时间[s]，ψ为基质吸力水头[m]，K(θ)为非饱和水力传导度[m/s]。该方程强非线性，一般无闭式解，由此产生的数值挑战（质量平衡误差、接近饱和时的收敛问题、伪振荡）已有综述讨论。闭合该方程所需的本构关系θ(ψ)和K(θ)最常以van Genuchten–Mualem模型表示：S_e = (θ ? θ_r)/(θ_s ? θ_r) = [1 + |α_vGψ|ⁿ]^?m，m = 1 ? 1/n；K(S_e) = K_sS_e^L[1 ? (1 ? S_e^1/m)^m]²，其中S_e为有效饱和度，θ_r和θ_s为残余和饱和含水量，α_vG和n为拟合土壤水分保持数据的曲线形状参数，K_s为饱和水力传导度，L为孔隙连通参数。这些参数对土壤质地和结构高度敏感，是非饱和流建模预测不确定性的主要来源之一。

3. 确定性与随机性建模方法

3.1 确定性方法

从确定性角度建立水动力过滤问题遵循既定步骤序列：确定研究域几何范围；为过滤介质参数赋值；指定流动域的初始和边界条件；求解控制方程以获得水头或含水量场。确定性水动力问题进一步分为正问题和反问题：正问题从已知的地过滤参数和边界条件定量求解地下流的动态函数；反问题则从观测的动态流函数值确定含水层或过滤层的水文地质参数，通常通过水头或流量实测值校准实现。

3.2 随机方法

从统计角度处理过滤问题需要识别水文地质问题并构建初始概念模型、对研究对象进行概化、选择建模参数并构建因变量初始矩阵。地统计方法将水力传导度等介质性质表示为空间相关随机场，以均值、方差和由稀疏野外测量推断的空间相关结构为特征；流和运移方程在给定这些随机场的条件下求解，更常见的是通过Monte Carlo模拟在等可能的介质实现系综上进行。该系综方法将介质表征的不确定性直接传播到预测的水头或流量场中。近期水文研究证实，区域尺度反演建模比较显示随机（地统计）反演能比确定性分区更忠实地再现水力传导度场的大尺度趋势和局部变异性，但代价是计算量大幅增加且需要更大的观测数据集来约束随机场参数。

3.3 比较与混合方法

实践中两种建模哲学均不能自足。最有效的土壤和含水层液体行为模型通常是随机与基于物理方法的结合：用地统计方法表征控制参数K、S_s、θ及保留参数的统计变异性，所得参数场或由之构建的初始非线性回归模型再用于约束基于物理的偏微分方程求解器，使基于物理和随机成分联合表征所研究的过滤现象。近期地下水建模综述也将确定性与随机方法视为单一建模工作流中的互补层次而非竞争替代方案，二者间的适当平衡取决于数据可用性、问题尺度及模型所支撑的决策。

4. 说明性案例研究：砂滤柱中的积水入渗

4.1 数值模型与模拟设置

为将前述比较具体化，本节将确定性与随机建模哲学应用于单一物理动机明确的问题：垂直砂滤柱中的积水入渗，直接类比水处理中的快砂滤池和砂质包气带的浅层补给。采用一维混合形式Richards方程，对1.2 m深柱体离散为121个节点，使用Celia、Bouloutas和Zarba的质量守恒修正Picard有限差分格式，后向Euler时间步进（Δt = 20 s），节点间传导度取几何平均。van Genuchten–Mualem本构关系采用Carsel和Parrish报告的砂类参数。柱体初始为均匀较干压力水头（ψ = ?1.0 m）；顶部为Dirichlet条件，ψ_top在前60 s内从初始值线性升至0 m以避免t = 0处的人工数值不连续；底部为自由排水（单位水力梯度）条件。确定性运行中K_s取单一均值均匀赋值；随机运行中K_s(z)表示为空间相关对数正态随机场，具有相同几何均值、匹配Carsel和Parrish报告的变异系数（σ_lnK ≈ 0.49）及相关长度为0.10 m的指数空间相关结构，通过协方差矩阵Cholesky分解生成50个独立实现构成Monte Carlo系综。

4.2 结果：确定性对随机性预测

湿润锋深度定义为有效饱和度S_e首次达到0.5处的深度。确定性模拟预测相对较快的前锋：湿润锋在0.62 h（37 min）到达0.5 m深度，1.5 h到达柱底，此时底部出流达到完全饱和通量（K_s = 297 mm/h）。随机系综均值在大部分模拟时间内滞后于确定性曲线，5th–95th百分位带很宽：在500成员系综中，湿润锋到达0.5 m所需时间从0.42 h（25 min；5th百分位）到1.43 h（95th百分位），系综均值为0.79 h（47 min）——仅由K_s的假定空间变异性产生的超过三倍离散度。底部排水通量方面，确定性模拟显示通量在1.5 h前锋到达柱底前保持为零，之后急剧上升至饱和水力传导度；随机系综则显示逐渐弥散上升，系综均值出流通量远低于确定性值，5th–95th百分位范围从零（低传导度层尚未克服的实现）到319 mm/h（已达稳态排水的实现）。在500个随机实现中，仅312个（62.4%）在模拟的2 h窗口内到达柱底，其余188个（37.6%）在t = 2 h时尚未开始排水——超过确定性到达时间的两倍。该右删失结果本身就是非均质性效应的直接定量例证。

4.3 数值验证：质量平衡、网格收敛及与Green–Ampt解的比较

全局质量平衡误差定义为最终模拟时刻（t = 2 h）顶部累计入渗、底部累计排水与柱内储水量净变化之间的相对偏差，以累计入渗的百分比表示。在固定时间步长下，质量平衡误差不随网格细化单调递减：将Δz从2.0 cm细化至1.0 cm使误差从1.07%降至0.34%，但进一步细化至0.5 cm又增至0.56%——这是空间与时间截断误差量级相当时的质量守恒格式的预期特征。网格细化序列显示到达0.5 m湿润锋深度的预测时间趋于稳定值，从61节点的1.39 h到121节点的0.62 h再到241节点的0.58 h。第三，数值累计入渗F(t)与Green–Ampt解析入渗方程比较，有效湿润锋吸力G_c = 3.8 cm由van Genuchten–Mualem参数经Morel-Seytoux和Khanji的理论关系独立计算而非拟合得到。两条曲线在早期差异最大，至模拟期末收敛至7%（相对）以内，RMSE = 39.3 mm，归一化RMSE为模拟期F值范围的6.5%。该一致性支持Richards方程实现用于随机系综的物理合理性。

4.4 扩展验证、敏感性分析与可复现性

湿润锋定义为S_e首次达到0.5处的深度，等价于θ = 0.237。对θ_s、θ_r和α分别扰动±1个标准差：扰动θ_s使到达0.5 m湿润锋的时间从0.62 h基线移至0.51 h和0.73 h（约35%摆动）；α产生约12%摆动；θ_r影响最小约6%。时间步长敏感性方面，Δt = 10、20和40 s下前锋到达时间和质量平衡误差均以与网格细化残余变化相当的幅度变化，且二者均不随Δt单调变化，作者明确披露该观测敏感性而非将其呈现为形式化收敛阶检验。Monte Carlo系综规模收敛性方面，500成员系综中4个实现（0.8%）在2 h窗口内未到达0.5 m湿润锋深度而被排除于该统计之外；50、100、250和496个有效实现的子集统计中，均值和95th百分位在约3%内一致，5th百分位作为极尾统计量噪声较大，变化达约6%。到达0.5 m时间的经验分布适度右偏（偏度1.13），排水通量和累计入渗接近对称。可复现性声明披露了Python 3.11实现、非线性迭代容差（max|Δψ| < 10^?7 m，每步最多60次迭代，不收敛时递归减半时间步长）、随机种子（5000起顺序编号）、σ_lnK = 0.49由Carsel和Parrish报告的砂类算术均值和标准差经对数正态矩关系推导而非假定，以及0.10 m相关长度代表砂质沉积物中常见的垂直相关尺度。补充的非平稳重尾敏感性运行（n = 100）结合了深度方向ln(K_s)均值线性趋势（跨度±1个标准差）和4自由度的多元Student t空间场：到达0.5 m深度的时间均值从匹配平稳基线的0.77 h移至0.51 h，偏度从1.17增至3.67，近地表0.2 m平均ln(K_s)与到达时间负相关（r = ?0.67），表明下细上粗趋势使近地表传导度升高从而加速早期入渗。

5. 近期进展：数值、计算与数据驱动方法

三项发展重塑了过滤和地下水流建模实践。第一，数值方法和计算硬件方面的进展使高分辨率三维瞬态模拟成为常规而非例外，下一代计算架构（GPU加速、百亿亿次和云计算）被认为需要相应新一代地下水模型加以利用，特别是用于本文第4节所用系综运行及随机和Monte Carlo方法普遍要求的那些运行。第二，大量案例研究的系统综述确认模型可靠性仍更多取决于用于参数化和校准的野外数据质量与密度，而非数值方案的 sophistication，这一挑战为确定性和随机方法所共有。第三，机器学习和深度学习方法日益被用作独立预测工具和嵌入基于物理框架的替代模型。机器学习模型在训练所覆盖的历史记录上可匹配或超过基于物理模型的预测精度，但其在观测条件之外的 extrapolation 能力以及与质量守恒原理的一致性仍是开放问题。这推动了物理信息和混合模型的兴趣增长，其中机器学习组件受控制过滤方程约束或与之联合训练而非作为无约束黑箱使用——例如物理信息神经网络已被用于直接求解Richards方程以估计非饱和水力传导度。这三条线索——更快更可扩展的数值求解器、更丰富的地统计非均质性表征、数据驱动替代或混合模型——应理解为互补而非竞争。

6. 讨论

6.1 与现有综述的定位比较

本文通过与2017至2025年间发表的七项相关研究对比来定位自身贡献。数值导向的综述如Farthing和Ogden以及Stankevi?ius等人的数值对深度学习求解器基准研究聚焦于固定控制方程的求解方法，未涉及确定性对随机性区分或从饱和到非饱和流的扩展；Govindarajan等人和Teng等人对Darcy定律本身的综述关注石油和低渗透情境中线性通量—梯度关系的有效性极限；Lozano Hernández等人的更广泛地下水建模综述定性地讨论确定性和随机实践，但作为文献调查而非算例，未量化特定问题的预测不确定性；Boo等人的机器学习综述和Osame等人的传感器综述各覆盖数据驱动预测和实验表征的补充角落。据作者所知，尚无现有综述同时结合饱和流理论、非饱和流理论、明确的确定性对随机性比较及原创案例研究这四要素，这正是本文填补的具体空白。

6.2 持续挑战与未来研究方向

表征地下非均质性仍是饱和与非饱和流模型中预测不确定性的最大单一来源，因为水力传导度和保留参数的空间相关结构很少以高保真数值模型所需的分辨率已知。案例研究表明，对单柱入渗问题，仅合理的适度非均质性即可使预测湿润锋到达时间偏移超过三倍。模型验证仍受地下监测成本高和稀疏性的制约，尤其在非饱和带。机器学习替代或混合模型与物理方程的结合仍处早期阶段，质量守恒的强制实施和纯数据驱动预测的物理可解释性尚存开放问题。随机处理还有三个更基本的局限：K_s(z)场按构造为对数正态，但一些沉积含水层的野外调查显示非高斯重尾分布和多尺度非平稳空间结构，单一对数正态指数相关场不能捕捉，因此预测结果包络宽度应视为真实不确定性的下界；Richards方程本身是基质流公式，不表征大孔隙、根通道或不稳定湿润锋指进中的优先流，其省略是模型本身的结构性局限；Monte Carlo结果仅量化通过假定结构正确的控制方程传播的参数不确定性，模型结构误差和先验约束不良的敏感性不由系综离散度本身解决。

7. 结论

本综述追溯了液体通过饱和多孔介质的水动力学理论，从Darcy原始实验及由此产生的稳态和瞬态流方程，经Richards方程扩展到变饱和条件。确定性和随机建模传统被证明是互补而非竞争的：确定性公式提供物理控制方程，随机和地统计方法提供表征和传播天然土壤和含水层固有参数不确定性的原则性手段。案例研究为这一比较提供了具体可复现的形式：对砂滤柱积水入渗，纯确定性模拟预测单一相对较快湿润锋到达，而由同一控制方程构建的经收敛性检验的随机系综（最多500次实现）表明，仅水力传导度的现实空间变异性即可产生超过三倍的到达时间离散度和可比宽度的排水通量范围。数值方法、高性能计算和机器学习的最新进展正与这一经典框架融合为混合的物理信息模型。统计与确定性方法的结合仍是模拟土壤剖面或含水层饱和与非饱和带液体流动的核心，未来工作的主要机遇在于更好地表征地下非均质性、将本文所示混合模型对野外数据验证，以及将数据驱动组件嵌入而非脱离本综述推导的守恒原理。