1. 研究背景与整体思路拆解

1.1 非饱和裂隙土降雨入渗难在哪

做非饱和裂隙土降雨入渗模拟这几年，一个经常被同行问起的问题就是：COMSOL里到底用van Genuchten模型（简称VG模型）还是Brooks-Corey模型（简称BC模型）？这个选择看似只是给软件填两个参数，实际上直接影响湿润锋的推进速度、裂隙周围的优先流形态，以及模型能不能顺利收敛。

裂隙土和非裂隙土最大的区别，在于水的迁移路径极不均匀。基质部分渗透系数很小，水走得慢；裂隙部分开度大、渗透能力强，降雨一上来就沿裂缝往下窜，形成典型的优势流路径。你可以把裂隙想象成高速公路，基质是旁边的乡间小路，降雨就是车流，大部分水会优先走高速路，只有边上的水慢慢渗进小路。这种双重孔隙结构使得饱和度和压力水头在空间上变化剧烈，如果用太简单的模型去描述，极易出现收敛失败或者入渗量偏大的问题。

1.2 为什么用VG和BC做对照

第二个常被问到的问题是：为什么不只用一种模型，非要做VG和BC的对照？原因很简单，国内大部分地勘报告和室内土工试验给出的往往是“饱和渗透系数+孔隙比”，真正能把非饱和土水特征曲线完整测出来的项目并不多。这时候你需要根据手头数据的精度选模型——BC模型参数少、形式简单，适合初算和量级判断；VG模型公式连续、曲线S形过渡平滑，更适合裂隙发育、孔径分布范围宽的土体。把它们同时放进一个COMSOL模型里做对比，是验证参数合理性最直接的手段。

这也决定了下面整篇文章的思路：先从数学模型讲清楚两个模型究竟差在哪，再把COMSOL建模的全流程串一遍，接着给出参数标定和结果差异，最后集中讲我实际运行中踩过的坑。无论是刚入门的研究生，还是已经开始做二维边坡入渗模拟的工程师，应该都能从中找到可以直接复用的内容。

2. 模型内核：VG与Brooks-Corey的数学与参数

2.1 两个模型的表达式与物理解释

先明确一点，非饱和渗流在COMSOL里最常用的是Richards方程接口，核心问题是压力水头h与体积含水率θ之间的关系。VG和BC其实都是描述这条“土水特征曲线”（SWCC）的经验公式，区别在于曲线形态的表达方式。

VG模型用有效饱和度Se作为中间变量：

Se = (θ - θr) / (θs - θr) = [1 + (α |h|)^n]^(-m) m = 1 - 1/n kr = Se^(1/2) * [1 - (1 - Se^(1/m))^m]^2

其中θs是饱和含水率，θr是残余含水率，α是进气压力倒数（单位1/m），n是孔径分布指数。α决定了曲线在压力水头坐标下的水平位置，n决定了曲线过渡段的陡峭程度。VG模型的特点是整条曲线处处连续、导数连续，数值上非常友好。

BC模型则是分段函数：

|h| >= hb 时: Se = (hb / |h|)^λ |h| < hb 时: Se = 1 kr = Se^(3 + 2/λ)

hb是空气进入值，理解为土体开始排水对应的吸力值；λ是孔径分布指数。BC模型在hb这个分界点处导数不连续，这是它数值上比VG难处理的原因之一。

从物理直觉上讲，孔径分布越均匀（比如干净砂土），BC模型越合适；孔径分布越宽、大小孔隙混杂（比如裂隙粘性土），VG的S形曲线更贴近真实。裂隙本身相当于给土体增加了一批“超大孔隙”，孔径分布被拉得很宽，这是裂隙土优先选VG的现实依据。

2.2 关键参数对照与取值经验

我整理了一张参数对照表，方便在COMSOL里设置参数时对照查看：

经验层面，裂隙粘性土的VG参数常见范围：α大约0.5~5（1/m），n大约1.3~2.5，θr在0.02~0.08，θs在0.35~0.45。BC模型的hb对粘性土往往在0.2~2 m水头，λ在1.5~4之间，砂土则更小。

实际项目中很少直接套用文献值。我之前帮某课题组的边坡模型做参数校核时，他们最初从一篇论文里抄了一组VG参数，α取0.1，结果湿润锋几个月都不往下走，明显不合理。后来换了室内试验数据重新拟合，α变成了2.7，入渗结果才正常。参数这东西，看着是小事，差一个量级模型表现就完全不同。

3. COMSOL实操建模：几何、边界与求解器

3.1 几何处理：裂隙的薄层等效

COMSOL建模的第一步是几何，裂隙怎么处理决定后续所有环节的走向。目前主流做法有两种：一是离散裂隙网络（DFN），把每条裂隙画成实体的薄层区域；二是等效连续介质，把裂隙的影响折算到渗透系数场里。

我个人在COMSOL里更推荐第一种做法，因为非饱和裂隙土模拟的核心就是要看裂隙周围的局部饱和度和压力变化，如果一开始就把裂隙平均掉，优势流就没法在结果里体现。具体操作时，把裂隙画成宽度1~10 mm的薄矩形域，材料单独赋予裂隙的渗透系数和孔隙率。裂隙渗透系数可以用立方定律近似计算：

k_f = b^2 / 12

其中b是裂隙开度（单位m），k_f的单位是m²，再乘以水的密度和重力加速度除以动力粘滞系数，换算成m/s。如果你只关心宏观入渗量，工程上也可以简单取k_f为基质渗透系数的100~1000倍，然后做敏感性分析。

这里有个细节：薄层区域的长宽比很容易拉爆网格质量。裂隙宽度1 mm，长度5 m，长宽比5000：1，COMSOL自动网格基本没法处理。我的做法是先用一个全局“最大单元尺寸”约束把网格加密，再在裂隙处加“边界层网格”。实测下来这种组合比直接给裂隙区域设超细网格更稳，也不会把计算量撑爆。

3.2 边界条件：降雨强度与积水压力的动态切换

降雨入渗的边界条件，新手第一次做就会踩坑。降雨强度是流量边界条件，在COMSOL里用通量边界实现，向下的入渗通量等于降雨强度I（单位m/s），即达西通量q = -I。

但降雨入渗的本质是：如果土体入渗能力大于降雨强度，边界保持通量控制；一旦地表达到饱和，多余的水不能立刻渗进去，边界会积水，此时边界控制条件从通量切换为压力边界，地表孔压为0。这个切换如果处理不好，模型会自动让地表压力水头变成正的大数值，物理上就失真了。

我在实际建模中的处理方式是分两步：

1. 先给定通量边界q = -I，计算到地表压力接近0的时间点；
2. 从这个时间点开始，把边界改成压力水头h = 0的定压力边界。

如果你不希望手动切，也可以在COMSOL里利用事件接口（Events）设置一个阈值：当地表压力水头>=0时，自动切换到定压力边界。实测下来事件接口需要注意时间步长的控制，如果不限制最大步长，事件检测时刻和实际物理时刻会偏差。

3.3 网格与求解器配置

COMSOL的Richards方程求解是非线性瞬态问题，最容易让人放弃的就是求解器设置。几个关键点：

初始条件不要给得太“干”。很多模型把初始压力水头统一设为-100 m，这在COMSOL里收敛非常困难。强烈建议先用一个较小吸力（比如-10 m）计算到准稳定状态，再把结果作为后续正式计算的初始值。这叫“分阶段初始化”，比直接硬算省太多时间。

时间步长建议设置初始步长1e-5 s级别，最大步长根据降雨时长控制。比如模拟24小时降雨，最大步长可以设600 s；模拟10天甚至更长，最大步长可以放宽到3600 s。太大会漏掉湿润锋推进的细节，太小则白算。

绝对容差我通常设为0.001或者0.01，单位是压力水头（m）。这个值不是越小越好，因为Richards方程在大吸力区间对数值误差非常敏感，容差过小会导致求解器在每个时间步反复迭代，计算时间膨胀十倍。

如果你用的是COMSOL内置的“Richards方程”接口，它默认的因变量是压力p（单位Pa）。定义VG表达式时，最好先把压力换算成压力水头：

Hp = p / (rho * g) Se = (1 + (alpha * abs(Hp))^n)^(-(1 - 1/n))

否则α的单位不是1/m就会出现一个很隐蔽的“参数对不上”的问题——COMSOL本身不报错，但结果完全不对。

4. 参数标定：把实测土水特征曲线转化为模型输入

4.1 从试验到参数拟合的流程

VG和BC参数不是拍脑袋定的，标准的来源是压力板试验或离心机试验，得到若干个“吸力-含水率”数据点，然后拟合出整条SWCC曲线。试验数据到我手里的通常只有十几个点，分散在吸力0.1 m到100 m的范围内。

拟合我一般用Python的scipy.optimize.curve_fit，不用在COMSOL里做拟合，因为COMSOL的优化模块做这类参数辨识也可以，但设置起来要额外花时间。流程是这样的：

1. 把压力板试验数据整理成两列：吸力h（m）、体积含水率θ；
2. 定义VG公式作为拟合函数，初值从文献范围猜；
3. 最小二乘拟合，得到α、n、θr；
4. 把拟合曲线和实测点画在一起，看R²和残差分布。

曾经有一组数据怎么拟合都不收敛，后来发现是压力板试验在低吸力段的数据点太少，整条曲线的形状被中间几个点“牵着鼻子走”。解决办法是给θs一个上下限约束，直接用实测饱和含水率固定θs，只拟合α、n、θr三个参数，瞬间稳定了。

4.2 拟合中的单位陷阱与裂隙土特有标定问题

单位问题我专门写一小节，因为太多人在这里翻车。压力板试验的吸力通常用kPa表示，而COMSOL的Richards方程里压力水头习惯用m。换算公式是：1 m水头等于约9.8 kPa。也就是说，一个进气值hb = 3.5 kPa的土样，按水头算就是约0.36 m。如果你把3.5当米输入，整个SWCC曲线相当于向右平移了一个量级，模型结果能对吗？

裂隙土还有一个特殊现象：裂隙开度不是固定的，它随含水率变化。干燥时裂隙张开，湿润后土体膨胀、裂隙闭合。这意味着同一块土体，降雨开始前和降雨过程中的渗透系数是完全不同的。如果你的研究对象正好是膨胀性裂隙土，单纯用恒定VG参数是不够的。我在某次模拟里把裂隙渗透系数设为常数，结果算到一半发现裂隙周围压力持续偏高，实际上如果裂隙闭合，雨水根本下不到那么深。

这种问题我建议做一步参数折减：让裂隙渗透系数k_f随有效饱和度Se单调递减，用表达式实现：

k_f = k_f0 * (0.1 + 0.9 * (1 - Se))

这样模拟开始后裂隙逐渐“自我封堵”，和实际土体表现更接近。当然这只是一种简化处理，真正的裂隙-基质耦合模型需要引入双孔隙/双渗透模型，那是另一个量级的课题。

5. VG与BC的模拟结果差异与选型建议

5.1 结果差异的机理分析

很多人在同一个模型里分别跑VG和BC，出来的云图看着差不多，就把这个步骤跳过了，其实数据层面差异非常大。我做过一个典型算例：二维土柱，中间有一条垂直裂隙贯穿，降雨强度恒定，模拟48小时。两组参数都做了土水特征曲线拟合，所以“宏观上”曲线重叠度很高。

实际结果差异体现在两处。第一，湿润锋的形态：VG模型算出的湿润锋过渡带更宽，裂隙两侧的饱和度分布平滑；BC模型因为进气值对应的导数突变，湿润锋在吸入段会略微“拖尾”，看起来就像有一层干土突然被水侵入，锋面更陡。第二，积水开始时间：BC模型由于相对渗透系数公式更简单，在低饱和度区域的渗透阻力衰减更快，地表更容易提前积水。

从数学上理解这件事，关键在相对渗透系数kr对饱和度的敏感度。VG的kr公式含有一项(1 - Se^(1/m))^m，在Se接近0或接近1时两端都趋于平缓；BC的kr是纯指数函数，在Se中等区间变化极快。裂隙土恰恰是基质和裂隙饱和度差异巨大的场景，基质区域的Se可能只有0.3，裂隙附近已经接近0.9，这种情况下两个模型的kr几乎差一个数量级。

5.2 对照思路与结果判断标准

那我建议怎么用这个对照？不要只盯着云图，我习惯把两种模型的结果放到同一组探针点数据里比较。COMSOL里在裂隙中点和基质深处各放一个点探针，输出压力水头或饱和度随时间的变化曲线。如果两条曲线趋势一致但峰谷幅度差异明显，说明你的模型对土水特征曲线的形态很敏感，参数标定工作需要重点复核。

另一种判断标准是质量平衡。COMSOL的“Derived Values”里可以积分边界入渗总量，再减去域内储水变化量，两者之差如果超过2%，结果就不可信。BC模型因为分段函数的导数不连续，在裂隙附近特别容易出现质量不平衡，一旦发现优先检查时间步长是否过大。

选型上我的经验是：正式研究尽量用VG，工程初算或者缺少残余含水率数据时用BC。另外BC模型在粗颗粒土、单一孔径介质里表现并不差，如果你模拟的是砂土基质的裂隙网络，BC反而简洁可靠，参数少意味着需要标定的试验点更少，这对实际项目来说是很现实的优势。

6. 常见问题与排查实录

6.1 收敛失败的典型原因与解法

我整理了实际运行中遇到频率最高的几个问题，做成速查表：

最常见的是第一类。我之前做深埋裂隙土模型，初始吸力设为-80 m，直接让COMSOL瞬态计算，结果连初始步长都过不去。后来我把问题拆成三步：第一步初始吸力-5 m计算到稳态，第二步把稳态结果作为初值，用-20 m吸力快速计算，第三步才接正式降雨边界。三维模型从完全跑不动变成半小时完成，这里的分阶段启动其实就是给非线性求解器一个相对接近真实解的初值，道理和用牛顿法迭代必须给好初值一模一样。

6.2 怎么判断模拟结果是不是可信

最后这部分是经验之谈，如果你做完一个模拟，只知道看云图颜色好不好看，那基本等于白算。我每次拿到结果会按顺序检查三件事。

第一是质量平衡。在COMSOL里定义全局积分算子intop1，算降雨边界总入渗量，再算整个域内含水率总量的变化。两者之差除以入渗总量就是相对质量误差。这一步在瞬态模拟里等于照妖镜。

第二是探针曲线的平滑性。正常情况下压力水头随时间变化的曲线应当是单调光滑的，如果某段出现高频抖动，说明时间步长在该阶段偏大，或者边界条件切换的瞬间数值算法在振荡，这时候我一般会回退到切换时刻附近，缩小最大步长重算。

第三是局部物理合理性。裂隙土模拟最常见的一个低级错误是裂隙区域饱和度始终处于0.9以上，而基质区域完全没有响应。如果模拟的是降雨强度很小的入渗过程，裂隙确实会优先饱和，但基质区域不可能完全无响应，至少应该有缓慢上升的趋势。如果基质区域始终一条直线，往往是裂隙渗透系数设得过大，或者基质区域的体积含水率参数设定有问题。

最后再说点个人的实际操作体会

分享一个我后来固定下来的工作模板。做任何一个非饱和裂隙土模型，我都会先在COMSOL里用同一个几何跑一遍BC模型，目的不是为了得到最终结果，而是用它来快速测试边界条件和网格设置是否合理。BC模型参数少，调试方便，哪怕中途报错也容易定位。等整个模型框架稳定了，再切到VG模型做正式计算。这个“先BC调试、后VG精算”的思路，帮我省了不知道多少迭代排查的时间。

另外提一句，如果预算允许，把COMSOL的土壤水分特征曲线直接和室内试验数据做交叉验证，远比从文献里拿近似参数可靠。若干年后回看这些计算结果，你会发现最花时间的往往不是模型本身，而是参数背后那些看似不起眼的土工试验细节，但它们恰恰决定了整个模型的成败。