两年前我接了一个雨养边坡的渗流场评估项目最初模型只在饱和带范围里跑地下水位以下的流网画得干干净净专业上也没人挑出毛病。可甲方的核心关注点是雨季滑坡风险这就逼着我把注意力从地下水位挪到地表以下那几米才意识到自己绕开了整个项目最关键的环节——非饱和区处理逻辑。这篇文章打算把这一块涉及的内容完整梳理一遍非饱和区的水分运动机制、Richards方程和土水特征曲线在其中的地位、数值模拟里边界条件和收敛性的处理逻辑再附一个完整的降雨入渗模拟案例复盘。适合水文地质、环境岩土、农业水文以及做渗流数值模拟的同行参考也适合刚接触这类模型的初学者建立整体概念。1. 非饱和区是什么包气带里的水到底在做什么运动1.1 先看清楚包气带的结构位置非饱和区也就是包气带是地面以下、潜水面以上的那一段地层。很多人做地下水数值模拟时习惯把模型顶面直接设在潜水面把潜水面以上当成一个“黑箱”只在顶部给个入渗补给量。这个简化在很多区域尺度问题里能用但一旦问题涉及降雨入渗、蒸发、植物根系吸水、污染物下渗这个黑箱就不能再黑了。包气带从结构上可以分成三个亚带最上面是土壤带厚度通常零点几米到几米是根系活动和大部分生物过程发生的区域中间是中间带厚度取决于地下水位埋深可以从几米到几十米最下面是毛细带紧贴潜水面水在毛细张力作用下被抬升到潜水面以上一定高度含水量相对较高。这三个亚带的水力特性差异很大建模时如果从土壤库选一个“平均参数”从上打到下湿润锋的推进速度、补给到达地下水的时间都会算偏。至少应按土壤分层和埋深条件做纵向分层。1.2 非饱和不等于干三相体系的物理现实非饱和区最容易被误解的地方是“非饱和”听起来像“没有水”。实际上包气带里的水不仅存在而且在部分地层还相当可观。孔隙中水和气同时存在就是典型的三相体系土壤颗粒为固相孔隙溶液为液相空气为气相。定量描述非饱和区含水量最常用的参数是体积含水量θ即单位体积土体中水所占的体积比。θ的上限是饱和含水量θs工程上近似等于孔隙度θ的下限通常不取0而是取残余含水率θr。残余含水率指的是土体再干也很难释放的那部分水以薄膜水和封闭孔隙水等形式存在。换句话说一块看着干透了的土里面其实还有百分之几的水。描述非饱和状态除了含水量还需要一个更关键的量——负压水头ψ也就是基质吸力。饱和区压力水头一般不低于0而非饱和区因为毛细作用和颗粒表面吸附的作用压力水头是负值。负压绝对值越大土越干。打个比方非饱和区的土就像一块被拧过水的海绵湿海绵轻轻一压就出水干海绵你使劲压也挤不出多少水这个“越干越难挤”的感觉在力学上就是基质吸力。1.3 为什么非饱和处理比饱和区麻烦一个量级做饱和带模拟导水率K和贮水系数通常是常量方程性质相对温和。到了非饱和区两个核心参数全部随含水量变化而且变化幅度极大。土壤在接近饱和时导水率可能是每小时几十厘米到了很干的阶段可能降到每小时10的负六次方厘米量级跨越七八个数量级。这种强非线性是数值模型发散的最常见根源。第二个麻烦是干湿交替带来滞后效应。土水特征曲线在吸湿和脱湿两个方向上并不是同一条线这会导致同样一个含水量对应的吸力不同。后面我会专门讲这个问题。第三个麻烦是优先流。包气带里存在动物通道、植物根孔、干缩裂缝水会沿着这些通道快速绕开基质向下运动。传统模型假定土体均质各向同性在野外经常算出比实测偏慢的下渗速度。这三个问题叠加在一起决定了非饱和区不是“给套参数就能跑”的模块而是需要专门处理逻辑的复杂系统。2. Richards方程与土水特征曲线非饱和区的两条支柱2.1 达西定律的非饱和扩展与Richards方程的由来非饱和区水流运动仍然服从达西定律只不过导水率不再是常数。饱和状态的达西公式是q -Ks·∂H/∂l到了非饱和状态写法变成q -K(ψ)·∂H/∂l。K从Ks变成了随负压水头ψ变化的函数。水头H仍然等于位置水头加压力水头只是这里的压力水头在非饱和区是负值。把达西扩展式代入连续性方程就得到描述非饱和水分运动的基本方程——Richards方程。以压力水头为自变量的一维垂直形式可以写成C(ψ)·∂ψ/∂t ∂/∂z[K(ψ)·(∂ψ/∂z 1)]这里z取向上为正C(ψ)叫容水度是含水量对吸力的变化率∂θ/∂ψ。容水度本质上描述的是土体在单位吸力变化下能释放或吸纳多少水。饱和区容水度在理论上趋于无穷所以直接用ψ形式处理饱和—非饱和耦合时数值上容易振荡。实际工程软件里很多模型选择保留∂θ/∂t项的混合形式。这么做的好处是在饱和区和非饱和区交界处质量守恒更好不会因为C(ψ)的剧烈变化产生额外的数值误差。这也解释了为什么在Hydrus这类软件中模型内部迭代总是会明确告诉你“当前用的是基于含水量的方程组还是基于压力水头的方程组”这个选项不是性能参数而是直接关系收敛性的物理考虑。2.2 土水特征曲线是打开非饱和问题的钥匙Richards方程看似简单真正难的地方在材料关系也就是土水特征曲线SWCC。SWCC建立的是吸力ψ与含水量θ之间的一一对应关系它是把“应力状态”和“储水状态”连接起来的桥梁。配合非饱和导水率函数整套模型才能闭合。工程上最常用的解析表达式是van Genuchten模型。含水量表示为θ(ψ) θr (θs - θr) / [1 (α|ψ|)^n]^m其中m 1 - 1/n非饱和导水率表示为K(ψ) Ks·Se^L·[1 - (1 - Se^(1/m))^m]^2Se是标准化含水量也叫有效饱和度。α、n是两个拟合参数。α的倒数大致对应进气值进气值可以理解为土开始排水时的临界吸力砂土的进气值通常很小因为大孔隙不擅长持水黏土的进气值较大。n控制SWCC的形状陡峭度n越大表示孔径分布越均匀含水量随吸力的变化越集中在一个窄区间砂土n可以达到2到3以上黏土n往往只有1.1到1.3。下面这组参数来自常用土壤数据库适合作为没有实测数据时的初始参考注意它只是经验平均值绝对不能替代场地实测。土壤类型θrθsα (1/cm)nKs (cm/h)砂土0.0450.430.1452.6829.70壤土0.0780.430.0361.561.04黏土0.0680.380.0081.090.24我自己习惯把SWCC做成一个Excel模板A列输入一系列负压ψB列按公式算SeC列算θD列算K。这样拿到一组实验室数据后先把数据散点图拉出来根据曲线拐点初估α和n再用规划求解或专门的拟合工具标定。比起一上来就丢给软件自动拟合“肉眼先看曲线再让机器优化”的方式更容易发现异常数据点。2.3 没有实验数据时怎么干活经验参数与敏感性意识很多项目时间紧张根本没条件做完整的压力板试验。这时参考值库确实是唯一选择但要清楚低估了什么。非饱和参数最大的特点是“异参同效”也就是好几组完全不同的α、n组合在拟合SWCC观测数据时都能给出几乎相同的曲线但对预测目标的影响却不同。比如预测降雨后是否产生地表积水Ks影响最大预测污染物多久到达潜水面n和θs的影响更大。所以做敏感性分析比盯着某一组参数本身更有意义。我的经验是先按质地类型选用一套基础参数然后围绕Ks上下浮动半个数量级、围绕n上下浮动20%各跑一遍模型观察关键输出变量变化幅度。如果输出变化大说明这个参数值得投入资源去实测如果输出变化小那就可以放过不要把预算浪费在对结果影响不大的参数上。3. 数值模拟里的关键处理逻辑边界条件、网格与收敛3.1 边界条件选择降雨不能简单套流量边界非饱和区模拟的边界条件处理是最容易暴露“处理逻辑”问题的地方尤其是降雨入渗的上边界。很多人想当然认为降雨就是一个流量边界把降雨强度直接当作表面入渗通量。这在降雨强度小于土体入渗能力时可以成立但一旦降雨强度超过入渗能力地表就开始积水这时表面不再是纯流量边界而是变成一个压力边界压力水头约为0允许积水时是积水深度。很多代码在处理这个问题时采用了一种动态切换机制即“大气边界”逻辑。大气边界是这样工作的在未积水状态下施加给定的降雨或蒸发通量当表层负压水头达到某个阈值比如-1cm到0cm表面自动切换成积水条件此时实际入渗量由下部土体能传导多少决定多余的水表现为地表积水或径流。反过来积水消退后边界又会切回流量控制。我在这个环节踩过很实在的坑。有一次图省事直接给上边界恒定通量降雨强度取2cm/h土层Ks只有1cm/h左右。模型为了守恒把表层负压不断往负向拉大计算结果里地表含水量始终偏低根本看不到积水而实际上这种雨强下早就该积水了。把边界改成大气边界之后积水时间、表层含水量剖面才与实测吻合。所以遇到降雨入渗问题先问一句我的软件支持这种通量—压力动态切换吗不支持的话就得自己判断什么时候切换边界。3.2 网格与时间步长为什么干湿交界处容易发散Richards方程是强非线性方程网格过渡区和初始条件对收敛性的影响极大。以入渗问题为例湿润锋是一个含水率从0.2跳到0.4甚至更高的陡峭前缘如果这条前缘的横向跨度只有一个网格计算就会非常勉强。处理逻辑通常是在湿润锋预计经过的深度加密网格特别是表层。一般来说地表第一个节点控制在5到10cm以内往下逐渐放松。一个5m深剖面可以把0到0.5m加密到2到5cm间距0.5m以下用10到20cm间距这样兼顾精度和算力。初始条件也不能太干。很多初学者把初始负压设到-1000m相当于若干旱土壤结果程序跑第一步就把时间步长压到10的负10次方小时量级基本算不出去。物理原因是太干的土壤导水率极低湿润锋处水力梯度又大数值处理难以同时满足局部流量平衡。更稳妥的做法是先做一次稳态求解让含水率剖面与底部水位平衡再以此为初始场叠加降雨。收敛还跟迭代算法有关。混合形式方程配Picard迭代是经典组合但在干湿交替剧烈时松弛因子放太大会振荡。我在项目里习惯把松弛因子控制在0.7到0.8而不是默认1。代价是收敛慢一点好处是稳健得多。最后必须盯质量平衡误差可靠模型一般应该小于1%超过2%就该回去检查网格、步长或者边界条件。3.3 滞后效应一个被低估的误差来源土水特征曲线并不是独一无二的。同样一个含水量土在湿润过程中对应的吸力与在干燥过程中对应的吸力并不相同。这个现象叫滞后或回滞。主要成因包括“墨水瓶效应”即孔隙截面忽大忽小在吸湿时先经过窄颈才能填满宽腔在脱湿时宽腔要先排空才能释放窄颈中的水以及接触角、截留气泡等因素。滞后在实际项目里影响多大单次降雨入渗过程影响通常可控因为整个过程主要在吸湿分支上进行。但如果是长期水量平衡模拟比如连续几年计算蒸发、降雨交替条件下覆盖层的水分动态忽略滞后可能让表层含水量的预测误差达到0.05到0.1体积比这对边坡稳定性或植被蒸腾计算就是不可接受的量级。如果软件支持滞后选项那就需要提供主吸湿曲线和主脱湿曲线两组参数软件内部再插值生成扫描线。不支持滞后或没有足够数据时至少要知道自己的结果是偏向吸湿分支还是脱湿分支的用单曲线估算时要说明取舍。4. 一个完整案例复盘降雨入渗模拟从头跑到底4.1 场景设定与模型搭建用一个典型场地来说明整个非饱和区处理逻辑。地层我取一层均质壤土厚度3m底部设为潜水面采用自由排水边界。降雨强度2.0cm/h持续8小时随后停止降雨并允许地表蒸发总模拟时长24小时。初始条件按静水平衡设置即底部压力水头为0向上逐步递减3m顶部位置压力水头约为-300cm。网格在上部加密0到0.5m间距2cm其余部分间距5到10cm竖向共约100个网格。这种设定和现实中很多垃圾填埋场覆盖层、浅层边坡的渗流分析场景接近。这里有一个容易被忽略的逻辑底部边界该用自由排水还是定水头。如果地下水位恒定且深自由排水边界更合理因为下边界允许土体在重力作用下把水排走不会人为抬高底部负压。如果项目关心的是潜水补给量那么底部可以考虑用单位梯度边界也就是让流出通量只受重力影响。边界的选择直接决定底部的渗漏速率这一步选错后面的结果都白做。4.2 参数选取与初始条件标定壤土参数采用上一节表中的参考值θr0.078θs0.43α0.036 per cmn1.56Ks1.04cm/h。初始剖面由静水平衡给出表层含水量对应负压-300cm时约等于0.16左右。这个数值在壤土上很典型。有人可能想直接给均匀含水量0.16而不是压力平衡剖面这会导致模拟前几个小时出现人为的排水或吸气过程掩盖真实降雨响应。所以初始条件不要图省事尽量用稳态解。把参数输入模型后先不急着开降雨先让程序在没有外荷载条件下跑一段时间确认含水率剖面没有明显漂移。这一步在各类数值模拟中叫“初始平衡检验”做透了后续结果才经得起推敲。4.3 结果判读与关键输出模拟结果大概表现出这么几个阶段降雨初期表层入渗能力很强因为表层负压大水被基质吸力牢牢吸住头一两个小时降雨基本全部入渗。随着上层含水量升高导水率逐渐接近饱和值入渗能力下降当降雨强度2.0cm/h超过Ks的1.04cm/h时表层很快达到饱和并积水。在Hydrus这类软件里可以在输出结果中看到表层通量在某一时刻从降雨强度自动跳变为略低于降雨强度的入渗值多余的降雨被程序计入积水或径流损失。湿润锋大约在8小时降雨结束前推进到0.6到1.0m深度下部土体仍然保持初始低含水量状态。停雨后表层开始排水和蒸发含水量回落到0.25到0.30附近但底部渗漏通量还会继续一段时间。这个“蓄水池加节流阀”效应正是非饱和区的典型行为——即使降雨停止了深层土壤仍能缓慢向下输送水分而不是像管道一样立刻截断。我通常会输出三个检查项质量平衡误差、表层含水率过程线、不同深度含水率剖面。质量平衡误差在1%以内说明数值方案可信表层含水率在雨后应当回落但不能出现负值或大于θs的异常值含水率剖面应保持单点连续不出现锯齿状跳变。4.4 我在这个案例里踩过的几个坑第一个坑就是前面说的恒定通量边界。早期我为了省事直接把降雨强度作为上边界通量结果地表始终不能积水而实际现场降雨后几小时就出现了地表径流。后来才明白降雨边界是具备切换逻辑的不是单纯的给定流量。第二个坑是初始条件设得太干。有一次我是想模拟极端干旱后骤降暴雨把初始负压设成-1000cm结果程序时间步长被压到几乎停滞而且湿润锋处的迭代收敛极慢。处理这个问题的方法是把初始负压设在-200到-400cm量级并允许程序使用自适应时间步长让湿润锋逐步推进而不是一开始就形成剧烈的干湿交界面。第三个坑是忽视了蒸发期的质量平衡。停雨后表层土壤迅速变干上部网格如果太粗蒸发吸力会造成局部负压极大值迭代容易发散。这时需要加密表层网格并适当限制最大负压范围。我记得有一次跑出来的表层含水量出现负值第一反应以为是bug后来发现是蒸发期网格太粗基质吸力在单个网格内被严重高估。加密后问题自然消失。5. 从包气带走向工程决策三个典型应用场景的延伸思考5.1 降雨型滑坡与基质吸力非饱和区处理逻辑在工程中一个非常重要的应用是降雨型滑坡分析。浅层滑坡通常发生在土体饱水后抗剪强度骤降的位置而湿润锋推进就是触发机制。非饱和土力学里有一个实用概念基质吸力为正时能增加有效应力相当于额外提供了抗剪强度。降雨使表层负压减小甚至转为正值吸力贡献消失边坡安全系数随之下降。正确的分析链路是先做非饱和渗流算出不同时刻的负压场和饱和度场再把这些结果传递到稳定性分析中用考虑吸力的抗剪强度公式去评估安全系数。难点在于稳定分析关心的是局部往往滑动面就在地下0.5到3m正好是包气带。如果模型还在饱和区边界上打转整个稳定性评估就失去了基础。5.2 渗滤液或污染物在非饱和区的迁移逻辑非饱和区是污染物到达地下水的最后屏障。传统评估经常简单地把地表污染源通量直接折算成地下水补给量这就跳过了包气带的缓冲作用。实际上包气带既能通过基质吸附和降解作用延缓污染物迁移也可能通过优先流通道快速突破。污染物模拟需要在Richards方程基础上耦合对流弥散方程再考虑吸附项和降解项。以垃圾填埋场渗滤液为例如果只做饱和带模拟计算结果往往会高估地下水浓度峰值出现的时间。把非饱和区加进去后渗滤液在覆盖层和包气带中的迁移速度受Ks控制深层污染到达潜水面的时间会明显推迟。但另一个方向也要小心如果场地存在裂缝或大孔隙优先流会让污染物走捷径这时候标准Richards方程的均质假设反而不安全。5.3 根系吸水与蒸散发生态水文模型植被与水的相互作用同样发生在非饱和区。根系吸水在Richards方程里表现为一个汇项经典的Feddes模型把吸水强度与根系深度、土壤负压关联起来吸力过大时根系吸水困难接近饱和时缺氧也会降低吸水。这个逻辑决定了农田灌溉需要什么时候灌、灌多少也决定了生态修复工程里植物能不能在旱季活下来。从这个角度看非饱和区处理逻辑不只是工程师和科研人员的专属话题农业水文学家、生态修复从业者每天都在跟它打交道。如果模型只关心饱和带水位不关注根区含水量变化就完全可能给出错误的灌溉决策或生态补水建议。最后说一个我在实际项目中始终坚持的小习惯无论软件跑出多么精细的云图我都会先用简单的入渗公式估算一下湿润锋深度、积水时间和总下渗量再和数值结果对量级。这个习惯帮我抓住过好几次参数标定偏离的问题。非饱和区处理逻辑的难点从来不在公式多复杂而在每一步取舍得是否清醒——边界条件怎么选、参数可靠性到什么程度、滞后和优先流要不要考虑这些决策决定了模拟结果到底是在解决问题还是在生产误差。