搞非饱和渗流的同行肯定都有这种感受单一 Richards 方程在室内土柱实验里拟合得漂漂亮亮参数一给曲线走得顺顺当当可只要一放到野外坡面、农田耕作层或者人为扰动的填土里实测和模拟就开始“闹离婚”。尤其是下了一场雨深层含水量探头突然猛跳而你用单一方程算出来的结果还在慢吞吞地推进——这种时候脑子里就得蹦出一个词优先流。裂隙、根系孔洞、虫洞、冻融裂缝这些大孔隙里水流的速度比你基质域中那个“精细扩散”过程快几个量级单一 Richards 方程根本描述不了这种“快慢两套系统”同时存在的水分运动。“双 Richards 方程双渗模型”就是冲这个问题来的。它把土体人为分成两个水力区域一个代表高导水能力的大孔隙/裂隙域一个代表低导水能力的基质域两个域各用一套 Richards 方程描述中间再加一个水量交换项用来模拟两域之间的水分交换。这套模型不是给实验室土柱做装饰用的它真正能打的场景是坡地降雨入渗、填埋场覆盖层防渗评估、矿山排土场含水量预测、农田排水这些“结构性强、优先流发育”的现实问题。这篇文章我会从模型初衷讲到数学骨架再到参数率定和实际操作中的坑全程用我自己的项目经验说话适合正在学非饱和渗流模拟的研究生、做水文地质的技术人员以及对优先流模拟有兴趣的工程设计人员参考。1. 为什么单一 Richards 方程不够用双渗模型的初衷先别急着上公式先把“为什么要拆成两个方程”这个问题掰扯清楚。只有理解了物理背景你在给参数的时候才知道自己在干什么。1.1 单域模型的困境基质与大孔隙的流速差了几个量级经典 Richards 方程假设土体是均匀多孔介质水流完全由基质吸力梯度和重力驱动导水率的变化由土水特征曲线决定。这个假设在装填均匀的砂柱、壤土柱里没大毛病但自然界里几乎没有这种“规规矩矩”的土。我做过一个试验场地表层是黏性土夹碎石往下 30 厘米处就有一层老根系通道雨后开挖断面能明显看到染色剂顺着根孔和干缩裂隙一路往下窜而周围的基质部分还是干的。这种时候用单一方程算出来的湿润锋是一个平整的“活塞式推进”界面实际却是水沿着几条高速通道迅速到达深层把深层含水量“顶”上来同时基质域还保持着较低含水量。从数值上看裂隙域或大孔隙域的非饱和导水率在接近饱和时可以比基质域高出 2~3 个数量级。也就是说两个域在同一吸力下水流的“运输速度”差了上百倍甚至上千倍。你把这两套水力行为强行平均到一套参数里结果必然是基质域被高估裂隙域被低估深层响应永远算不出来。更麻烦的是单一方程模拟出的累积入渗量虽然可能做“整体守恒”但它给不出“水优先走哪儿”“基质被晾在一边”这种空间结构性信息。而很多工程问题比如边坡降雨入渗的深层孔压响应、填埋场覆盖层是否会发生短路渗漏恰恰取决于这种“通道流”行为。1.2 双渗模型的本质两套“水路”同时跑还互相换水双渗模型Dual-Permeability Model的思路很直白把同一空间位置上的土体拆成两个重叠的连续介质——域 f快域代表裂隙、大孔隙通道和域 m慢域代表基质、微孔隙。每个域有自己的含水量、自己的压力水头、自己的水力参数各自满足一套 Richards 方程。运动方程在时间项里多了一个体积比例权重 w这 w 就是快域在整个土体体积中的占比。比如 w 0.1意味着这个土体里大孔隙通道占了 10% 的体积剩余的 90% 是基质。两个域之间并不孤立它们通过一个交换项 Γ_w 来传递水分。水势高的一侧会向水势低的一侧流动这个交换速率取决于两域的水势差也取决于接触界面的几何形态和传导能力。注意双渗模型和双孔隙Dual-Porosity模型最大的区别就在这里双孔隙模型把裂隙域当成“管道”只输水、不储水或者只做瞬时平衡而双渗模型里裂隙域自身也是储水介质它自己也能持有大部分水量有自己的含水量变化和入渗过程。这就让它能模拟那种裂隙域先快速入渗、随后缓慢向基质域“反哺”的过程而不是一进裂隙就瞬间全灌给基质。这个“自身储水相互交换”的双通道结构是双渗模型在现实中更可靠的核心原因。当初我在某人工填土边坡项目里第一次用双渗模型就是因为我们在覆盖层里埋了含水量计雨季实测数据始终显示表层的饱和带出现后短时间内深层含水量就有响应但响应幅度又不像裂缝完全贯通时那么剧烈。这种“快又不太快、剧烈又不太剧烈”的特征双孔隙模型给不出来单一 Richards 方程更给不出来只有双渗模型能把“快域主导输运、慢域缓冲储水”的中间态行为描述清楚。2. 双 Richards 方程双渗模型的数学结构与物理含义有了物理图像再上公式就顺理成章了。这一部分我会把方程骨架和交换项的细节拆开讲重点解释参数为什么长成这样、数值上为什么难搞。2.1 两个 Richards 方程加一个交换项长什么样采用最常见的以压力水头 h 为变量的混合形式双渗模型的完整方程组可以写成快域域 f方程[ \frac{\partial \theta_f(h_f)}{\partial t} \frac{\partial}{\partial z}\left[ K_f(h_f)\left(\frac{\partial h_f}{\partial z} - 1\right)\right] - \frac{\Gamma_w}{w} ]慢域域 m方程[ \frac{\partial \theta_m(h_m)}{\partial t} \frac{\partial}{\partial z}\left[ K_m(h_m)\left(\frac{\partial h_m}{\partial z} - 1\right)\right] \frac{\Gamma_w}{1-w} ]其中下标 f 和 m 分别代表快域和慢域θ(h) 是土水特征曲线K(h) 是非饱和导水率函数z 以向上为正Γ_w 是两个域之间的单位体积水量交换速率。总含水量是两者的加权和[ \theta_{total} w\theta_f (1-w)\theta_m ]如果你把两个方程加起来交换项会刚好抵消方程在总量上保持质量守恒。这是双渗模型非常优雅的地方——无论你在两域之间怎么“倒水”全系统的水量都不会凭空增加或减少。每个域的 θ(h) 和 K(h) 通常都用 van Genuchten 模型来描述。这里容易踩一个坑很多人直接把同一套 VG 参数复制给两个域只把饱和导水率 K_s 调大调小这是不对的。快域和慢域的区别不仅是导水速率更是孔隙结构快域对应大孔隙孔径分布集中在一个较宽的尺度范围进气值通常更低曲线形态更“陡”慢域对应微孔隙进气值较高曲线形态更缓。所以两套 VG 参数α_vg、n、θ_r、θ_s、K_s原则上应该独立给定而不是简单乘个系数。实际操作中如果缺乏直接的分域实测数据可以先从土水特征曲线反演一个“全土平均参数”再结合裂隙几何统计来拆分这个我放在后面参数率定部分详细说。方程的数值难点也很明显快域接近饱和时K_f 很大压力水头 h_f 变化一点点水量剧烈变化而慢域这时候可能还处在很干的贫水段导水率低到接近 0两者“脾气”完全不同。两个方程耦合在一起时间步长要是按照慢域的“慢性子”来取快域会浪费大量计算要是按照快域的“急性子”来取慢域的低导水率又会让隐式迭代产生严重的数值弥散。这个矛盾是双渗模型实现时最容易让人崩溃的地方后面我会给到工程上常用的折中策略。2.2 交换项的物理标定从达西近似到经验系数交换项 Γ_w 是整个双渗模型的灵魂也是最难标定的部分。不同文献里 Γ_w 的形式五花八门但主流可以归成两大类。第一类是基于达西近似的“几何物理型”。假设快域和慢域之间的交换发生在特定几何界面上接触面积为 A_s特征长度水流从快域中心到慢域中心的平均路径为 l界面处的有效导水率为两域导水率的某种加权平均于是交换通量可以写为[ \Gamma_w \frac{A_s}{V} K_{int} \frac{h_f - h_m}{l} ]这里 K_{int} 通常取两域导水率的算术平均或几何平均A_s/V 是单位体积土体内的接触面积。这种形式的好处是每一项都有物理含义坏处是 A_s 和 l 很难直接测量尤其是裂隙网络完全随机分布时你根本不知道“两块介质之间的有效距离”到底取多少。不过有一个工程上的近似经验对于常见的裂隙-基质系统l 可以取裂隙间距的一半A_s/V 可以粗略用裂隙的表面积密度来估算。我做过一个试验场地裂隙平均间距约 0.1~0.2 米按上述方法估算出的 A_s/V 落在每米 20~60 平方米区间经过率定后效果不错。第二类是把 Γ_w 写成水势差的线性或非线性函数直接用一个“传递系数” α_a 来打包所有几何和界面因素[ \Gamma_w \alpha_a (h_f - h_m) ]α_a 的量纲是 1/s可以理解为“单位水势差驱动下两域之间每秒交换多少体积水”。这种写法简洁高效实际模拟中很好用但 α_a 完全是个经验参数物理意义不直观只能靠反演率定。常见的率定范围我遇到过从 1×10⁻⁵ 到 1×10⁻² 1/s 不等具体值取决于介质渗透性和裂隙发育程度跨了三个数量级。这说明什么说明你拿一个区域的标定值硬套到另一个区域大概率会翻车必须对每个场地单独标定。我自己的体会是如果目标是做趋势分析或者机理研究用简单的线性交换项就够了如果目标是定量复现某个场的深层含水量峰值那必须给交换项留出足够的自由度甚至可以考虑让 α_a 随含水量变化。但参数越多反演越不稳定这在率定阶段要有心理准备。还有一个实操技巧在初跑模型时先把 Γ_w 暂时设为零让两个域完全独立跑通确认两域的单独行为都合理后再逐步打开交换项从很小的 α_a 往上调。这样可以把“两域自身参数错误”和“交换参数错误”两类问题解耦排查效率至少翻一倍。3. 参数获取与率定这个模型最“脏最累”的环节双渗模型能不能出可靠结果说实话七分靠参数三分靠求解。这一章我讲三件事两套水力参数从哪来、权重系数 w 怎么定、交换系数怎么率定。全程基于我踩过的坑没有教科书式说教。3.1 两套水力参数怎么分室内实验、裂隙几何与反演理想情况下你应该分别对裂隙域和基质域做独立的室内渗流实验拟合各自的 VG 参数。但现实里你很难单独取出一块“纯裂隙”土样因为裂隙是结构性的取样过程往往把它破坏了。退而求其次的做法是用原状土样做全土的土水特征曲线和导水率实验再用重塑扰动土样打碎过筛后压实代表基质域参数两者相减得到裂隙域的贡献。这个方法不完美但工程上比拍脑袋强得多。具体操作上基质域的核心参数 θ_s、θ_r、α_vg、n、K_s 可以通过压力板法或高速离心机法测出。裂隙域的参数里θ_s_f 可以认为是全土饱和含水量减去基质域在接近饱和时的含水量或者直接用裂隙体积率来估计K_s_f 则很难直接测常见做法是用双环入渗试验得到表层综合饱和导水率再根据裂隙面积占比反推快域导水率。比如现场测出综合 K_s 是 1×10⁻⁴ m/s裂隙占比 w 是 0.1基质 K_s_m 是 5×10⁻⁶ m/s那快域 K_s_f 至少要达到约 1×10⁻³ m/s 量级才能撑起这个综合值。这个换算关系虽然粗糙但能帮你判断快域参数有没有给得离谱。另一个被很多人忽略的参数是 n 值孔径分布形状参数。快域的孔隙结构以宽分布的裂隙和大孔隙为主n 通常偏高比如 1.6~2.5基质域微孔隙相对均质n 反而偏低比如 1.2~1.6。我见过有项目把两个域的 n 都设为 1.3算出来快域怎么样都表现不出“快速响应”一查原来是 n 太小气值区间拉得太宽大孔隙的优势被稀释了。改回 2.0 后深层的响应时间直接提前了 40 分钟效果立竿见影。如果你手头没有充足的分域实验数据还有一种实用的过渡方案先做全土 VG 参数反演得到整体参数后按下面的经验规则拆分——快域的 α_vg_f 取整体值的 2~5 倍进气更快快域的 K_s_f 取整体值的 10~50 倍快域的 θ_s_f 等于 w 对应的饱和含水量上限通常 θ_s_f 直接取 1 或者略低于 1 也行因为裂隙在饱和时几乎不残留空气。这套规则是很多文献里的默认做法实测下来作为模型的初始猜测值是够用的。3.2 权重系数 w 和交换系数 α_a 的标定心得w 是双渗模型里最具物理直观性的参数它代表快域体积占比。获取方式主要有三种染色示踪实验、CT 扫描图像分析、渗透实验反演。染色示踪最便宜实用下完雨后开挖剖面统计染色面积占总剖面面积的比例这个比例基本就是 w 的一阶估计。我在某试验场上做过 5 个剖面的染色统计w 稳定在 0.08~0.15 之间和渗流反演的结果非常接近说明这个方法很靠谱。但要注意w 不是一成不变的。表层耕作层根系发达裂隙占比可能 0.2 甚至更高到了深层裂隙密度下降w 可能衰减到 0.05。所以如果模拟剖面对垂直分层敏感w 应该设置成随深度变化而不是全剖面一个固定值。最简单的做法是分成表层、过渡层、深层三段每段给一个 w中间线性插值。我见过有的项目为了省事整个剖面 w 设为一个常数 0.1结果浅层优先流模拟得不错但深层含水量明显滞后这就是权重系数没跟着结构走导致的。交换系数 α_a 的率定就更依赖反演了。我用过手动试错也用过伴随法自动优化实际体会是手动试错时先固定其他所有参数把 α_a 按 10 倍对数步长从 1×10⁻⁵ 扫到 1×10⁻²观察深层含水量时程曲线。你会发现一个明显的现象α_a 太小时深层响应峰值明显偏低两域各跑各的α_a 太大时两域水势迅速拉平模型行为退化成一个等效的单域模型深层响应曲线又变得太“平滑”。真实值一定在这两者之间响应恰到好处是特征峰值时间对峰形不突兀也不拖沓。找到这个范围后再缩小步长做精细标定比一上来就做全局优化高效得多。在自动反演方面我踩过一次大坑把 α_a 和 K_s_f 同时设为待估参数结果反演出来的参数组合虽然拟合精度很高但快域导水率被推到了物理上不合理的 1×10⁻¹ m/s交换系数则被压得很低纯粹是数值补偿。从那以后我坚持“少参数、强约束”每次自动率定最多加两个目标参数并且一定要给 K_s_f 和 α_a 设置物理上下限。宁可拟合优度稍微差一点也不能让参数跑到荒谬区间否则模型一换场景立刻露馅。4. 实操过程与典型场景实现从一个坡面算例说起理论说再多不如实际跑一遍。这一章我用一个模拟的坡地降雨入渗算例把双渗模型从概念到实现的全流程走一遍。场景设定为20 米进深的二维坡面表层 2 米是含碎石裂隙的覆盖层以下为相对完整的低渗透地层模拟一场持续 12 小时、峰值强度为每小时 25 毫米的设计降雨。4.1 从概念模型到可计算的离散模型拿到一个场地后第一步不是建模而是把“概念模型”画清楚哪些区域需要分域分域范围怎么划。我习惯先在剖面图上标出裂隙发育区域和基质主导区域。对双渗模型来说不是全剖面的每个点都要跑两套方程你可以只在覆盖层范围内激活双渗区域下方完整地层用单一 Richards 方程就行。这样能省不少计算量而且物理上也说得通完整地层本来就没有像样的裂隙网络强行给一个低 w 值意义不大还白白增加非线性源项。模型的空间离散方面网格尺寸要照顾到两域的水势梯度差异。快域导水率高局部水势梯度大在降雨入渗前锋附近的网格要适当加密纵向步长我通常控制在 0.02~0.05 米横向根据剖面宽度取 0.2~0.5 米。时间步长的控制我用的是自适应策略初始步长设 1 秒降雨强度大时自动压到 0.1 秒量级降雨间歇期逐步放大到 60 秒。这里要特别提醒一个困扰过我好几个项目的问题交换项导致两域水势差在降雨初期突变显式处理交换项时步长稍大就会出现振荡后文问题排查部分会讲具体应对方法。边界条件的设置也有讲究坡面处按降雨边界处理降雨强度超过当前入渗能力时多余的雨水作为坡面径流移走程序里通常以“允许积水、超过阈值后转为定水头/径流”的方式实现底面按自由排水边界两侧按零通量边界。初始条件建议先用“稳态无降雨”预跑一段时间让两域含水量达到一个平衡状态再加载降雨。如果没有这一步预热初始水势场和交换项之间会有剧烈的“磨合”前几个时间步经常出现含水量跳变这种数值假象会毁掉对早期响应的判断。4.2 一个坡面降雨入渗算例的操作路径具体建模过程我分五步走。第一步定义材料参数。覆盖层内激活双渗按前述方法给出两套 VG 参数和 w。我给一个实际用过的初始参数组合供参考快域 θ_s_f0.38、θ_r_f0.03、α_vg_f3.5 m⁻¹、n_f2.0、K_s_f8×10⁻⁵ m/s慢域 θ_s_m0.42、θ_r_m0.05、α_vg_m1.2 m⁻¹、n_m1.4、K_s_m5×10⁻⁶ m/sw 取 0.12交换系数 α_a 初始取 1×10⁻⁴ 1/s。这个组合的特点是快域进气值低、导水率高、孔径分布宽明显区别于慢域。第二步定义初始条件和降雨序列。预跑 480 小时稳态后两域初始水势保持在剖面顶部附近约为 -1.5 米水柱底部由于排水条件形成自然的水势分布。降雨序列按小时数据输入峰值在开始后第 4 小时雨强 25 mm/h要注意程序里的降雨强度单位一般是 m/s不要被 mm/h 搞混换算的时候直接除以 3.6×10⁶。第三步设收敛控制参数。这个太关键了。双渗模型的非线性强度远高于单域方程我一般把两域含水量迭代的容许误差设为 1×10⁻⁴水势迭代误差设为 0.005 米水柱。更关键的是迭代次数上限不要设得太死给到 40~50 次同时开启线性搜索如果求解器支持。另外建议启用“分域交替迭代”即先固定慢域水势求解快域再固定快域水势求解慢域反复多轮。这种方式比网格式的完全联立求解在稳定性上好很多代价是计算时间稍微增加但换来的是少掉一半以上的中途报错。第四步运行计算与结果提取。运行结束后的第一时间我不是直接看整体含水量云图而是先看一个关键点坡面以下 1.5 米处两域的含水量时程曲线。这个点同时包含了快域快速响应和慢域缓释信号能直观反映模型的“双通道”行为是否真的跑出来了。如果两条曲线几乎重合说明交换项太大或两域参数区分度不够模型退化成了单域这是一眼就能看出来的“模型失效信号”。第五步输出结果做水量平衡检查。查看累积入渗量、累积径流量、剖面蓄水增量的差质量平衡误差应该在 1% 以内才算能接受。我遇到过好多次总体平衡很好但局部振荡剧烈的情况这是因为交换项在局部产生了虚假水量搬运应当回到参数或时间步设置排查。4.3 常见问题与排查技巧实录表格式的排查经验放在这里帮助大家少走弯路。现象可能原因处理办法早期不收敛错误提示集中在降雨初刻降雨边界与初始水势突变交换项振荡减小初始时间步到 0.01 秒增加预跑时间交换项隐式化两域含水量曲线几乎重合α_a 设置过大或两域 VG 参数区分度过小检查 α_a按 10 倍步长下调核实快域 n、K_s 是否明显高于慢域深层含水量峰值太低α_a 过小交换不足或 w 太小快域水通量不够增大 α_a检查 w 是否与染色试验统计相符峰值时间偏早、曲线陡峭异常快域占比 w 偏大或快域 K_s 偏高下调 w复核裂隙占比统计压缩快域 K_s质量平衡误差大时间步长切换突变或积水边界迭代未收敛固定降雨初段时间步长关闭大跨度步长倍增参数反演得到离谱数值同时率定参数过多目标函数与物理约束失配减少待估参数设置上下限分阶段率定除表格里这些常规问题外再分享一个深一点的技巧当模型在降雨结束后的排水阶段反复报收敛错误时不要只在数值上想办法。这个时段两域水势差方向反转交换项从“快域向慢域补给”变成“慢域向快域回流”源项方向突变导致雅可比矩阵阶跃。省时的办法是识别到这个物理过程后人为把该时段的时间步长上限压到降雨阶段的 1/5让数值解器平稳翻越这个“交替换向”的拐点。用这个办法我救回过好几个即将放弃的项目。5. 双渗模型的适用边界与后续扩展模型越复杂越要清楚它什么时候值得用、什么时候是杀鸡用牛刀。这一章说适用边界也说模型结果怎么用。5.1 什么时候该上双渗什么时候别硬上双渗模型不是万能的。我之前做过一个均匀黄土场地的入渗实验没有明显的裂隙和根系孔道污染羽中心迁移完全符合层状介质的活塞推进规律。这种情况下用单一 Richards 方程就够了强行上双渗只会增加参数不确定性反演出的 α_a 也毫无物理意义。所以判断标准其实很简单你的场地存不存在显著的大孔隙/裂隙通道存不存在“快通道输运与基质缓释并存”的现象如果没有就别为了论文“创新点”给自己找罪受。反过来下列情形就是双渗模型的典型适用场景一是裂隙性黏土边坡降雨极易沿干缩裂隙快速下渗诱发深层孔压响应二是农田耕作层根系孔和蚯蚓洞导致的优先流影响农药和养分运移三是垃圾填埋场覆盖层施工压实产生的软弱界面上次生裂隙可能形成渗漏短路四是矿山排土场大型机械堆排的松散土体中大孔隙分布极为广泛。在这些场景里双渗模型给出的结果和单一方程往往是“质”的差别而不只是数值上的微调。我再强调一次不要被“两个方程听起来很高级”冲昏头。双渗模型的参数敏感性比单域模型高一个量级对数据质量的要求也高得多。如果项目工期紧、现场数据少我建议直接用带传递函数的简化双渗形式或者先做单域模拟做灵敏度分析确定裂隙通道确实是控制因素后再引入双渗。这种“先证明必要性再增加复杂度”的思路是我这么多年做水文地质模型不翻车的基本原则。5.2 从模拟结果到决策双渗能多告诉我们什么双渗模型最大的价值是它能告诉你“水走哪条路、各走了多少”。举个例子某项目的边坡设置了深层排水孔单一 Richards 方程算出的渗流量预测比较平均但实际排水孔在雨后短时间里涌出的水量远大于平均流量且水质浑浊说明大部分的水是从裂隙通道中直接涌出的。改用双渗模型后模拟结果显示降雨事件中通过快域到达排水孔的流量占总涌水量的 70% 以上且到达时间集中在雨后 2~4 小时内。这个数字直接影响到了排水系统的设计——如果只按平均渗流量设计排水能力很可能被低估。另一个常见应用是填埋场覆盖层防渗性能评估。双渗模型可以追踪“短路”水流的路径占比从而量化不同压实和覆土方案下优先流触发的可能性。我曾经用这个模型比较过两种覆盖层结构结果显示表层含碎石层的方案下水分绕开低渗透层通过侧向裂隙进入下部的比例降低了近一半。这个结论对方案决策起到了关键作用而这种“优化结构而不是单纯加厚”的思路如果没有双渗模型的通道流量输出根本无从谈起。当然模型输出也只是决策的输入之一。双渗模型算出来的“快域流量占比”这类指标最终还需要和现场监测、工程经验对照印证。我的习惯是每一轮模拟都留出一部分实测数据不参与率定专门用来做验证。如果验证期表现稳定再把它作为方案比选的依据如果验证期翻车就回去查参数。这套流程虽然慢但能保证模型不是一个只能在论文里跳舞的玩具。最后说一点个人体会。做双渗模型这几年我最深的感受是模型越复杂你越不要指望“把参数输进去按一下运行”就出真知。它的价值不在于公式长得好看而在于逼迫你去理解你的场地到底长什么样——裂隙怎么分布水流走哪条路快慢过程怎么相互作用。哪怕有一天你用了更复杂的三相多尺度模型这种“先看清结构、再选择模型、最后才信结果”的方法论都不会过时。如果你正打算在项目里尝试双渗模型我强烈建议你从小规模试验场开始先把参数标定路线走通再往大尺度推广这个顺序反着走你在跑通之前就耗尽了所有耐心。
