当前位置：首页 > news >正文

MATLAB一维土壤水分运动模拟脚本：Richards方程差分求解器

news 2026/6/12 18:55:38

本文还有配套的精品资源，点击获取

简介：这个MATLAB资源包提供一个可直接运行的一维Richards方程数值求解工具，专注非饱和带土壤水分运移过程模拟。核心脚本solve_one_dimension_Richards.m采用有限差分法离散控制方程，支持显式与隐式两种时间推进格式，用户通过修改脚本内变量即可设定初始含水率分布、上下边界类型（定水头、给定通量或自由排水）以及关键土壤参数——包括水分特征曲线和非饱和导水率函数。不需要额外配置文件或接口，所有输入均在脚本头部集中定义。程序输出为各时间节点下空间网格点的体积含水率序列，结果以数组形式返回，便于后续绘图（如含水率随时间/深度变化曲线）或接入其他模型做耦合计算。配套包含示例土壤数据文件Hor_SL_L_0.5m.txt，可用于快速验证和教学演示。代码结构清晰，关键离散步骤附有物理含义注释，适合算法理解、课堂演示及中小尺度工程估算场景。

1. 项目概述：为什么一个“能跑通”的Richards方程求解器比想象中更难写

你有没有试过在MATLAB里敲下第一行theta = theta0;，然后盯着屏幕等它算出土壤里那滴水到底往下渗了多深？我干过——而且不止一次。十年前第一次用有限差分法解Richards方程时，我花整整三天才让程序不报错；又花了两周，才确认它输出的不是数学幻觉，而是物理上说得通的含水率变化曲线。这不是因为公式太复杂（其实核心就那一行偏微分方程），而是因为非饱和带水分运动本身就是一个处处藏着“非线性陷阱”的系统：导水率K(θ)随含水率θ指数级变化，水分特征曲线h(θ)在干燥端陡得像悬崖，而时间步长一选大，显式格式立刻发散；一选小，算一天才推进两小时。更别提边界条件切换时的数值震荡——比如上边界从定水头突然切到自由排水，模型常会“打个嗝”，在表层生成一段虚假的干裂带。

这个脚本solve_one_dimension_Richards.m，就是我从那些坑里爬出来后，用最朴素的方式重写的“防坑版”求解器。它不追求工业级精度（比如自适应网格或全隐式牛顿迭代），但保证三点：**第一，所有参数都在脚本开头50行内集中定义，改个初始含水率、换种土壤类型，不用翻三四个配置文件；第二，显式与隐式两种格式共存，你可以用显式快速试算、调参，再切到隐式做正式模拟；第三，每个差分离散步骤旁边都标着物理含义——比如K_half(i) = 0.5*(K(theta(i)) + K(theta(i+1)))这行，注释直接写“节点i与i+1间界面处的调和平均导水率，体现达西定律在离散界面上的连续性要求”。它不是黑箱，而是一张摊开的演算草稿纸。

关键词里提到的“Richards方程”“土壤水分模拟”“MATLAB差分”，恰恰对应三个层次的需求：科研人员要验证新提出的K(θ)函数是否合理，农业工程师想估算某次灌溉后水分入渗深度，而高校教师需要一个学生能在两节课内看懂、修改、并画出动态剖面图的教学工具。这个脚本就是为这三类人写的——它不替代HYDRUS，但当你只需要知道“第3天时40cm深处含水率是不是降到萎蔫点以下”，它比打开一个带GUI的商业软件快五倍。配套的Hor_SL_L_0.5m.txt是真实砂壤土的实验室测定数据，不是合成的光滑曲线，你拿它跑一遍，看到含水率前沿在前6小时推进缓慢、之后加速，再在24小时后趋于平缓——那种“啊，这就是毛管势在起作用”的实感，是任何理论推导给不了的。

2. 核心原理与方案设计：为什么必须同时提供显式与隐式格式？

2.1 Richards方程的物理本质与数值困境

Richards方程描述的是非饱和土壤中水分在重力与基质势梯度共同驱动下的运动，其一维垂直形式（以体积含水率θ为因变量）为：

$$
\frac{\partial \theta}{\partial t} = \frac{\partial}{\partial z} \left[ K(\theta) \left( \frac{\partial h}{\partial z} + 1 \right) \right]
$$

这里h是压力水头（负值表示吸力），K(θ)是非饱和导水率，z轴正向向上。方程右边括号内是达西通量q = -K(∂h/∂z + 1)（负号被吸收进坐标定义，此处按常规取向下为正）。关键难点在于：K和h都是θ的强非线性函数，且通常只能通过经验模型（如van Genuchten模型）给出解析表达式。这意味着方程本质上是一个关于θ的非线性抛物型偏微分方程。

数值求解时，非线性会放大两个经典问题：稳定性与收敛性。显式格式（如前向欧拉）计算快、每步只需矩阵乘法，但时间步长Δt受CFL条件严格限制：Δt ≤ Δz² / (2·max(K))。对砂土（K可达10⁻⁴ m/s），若空间步长Δz=0.01m，Δt不能超过5秒——模拟7天就得算12万步，内存和耗时爆炸。而隐式格式（如后向欧拉）无条件稳定，Δt可设为1小时甚至更长，但每步需解非线性方程组，传统做法是线性化（如Picard迭代），而线性化初值选不好，迭代就发散。

2.2 本脚本的双轨制设计逻辑

solve_one_dimension_Richards.m没有在“快”与“稳”之间妥协，而是把两者做成可切换的平行轨道：

显式轨道（scheme = 'explicit'）：用于参数敏感性分析与快速原型验证。比如你想测试不同初始含水率θ₀对入渗速率的影响，显式格式让你在30秒内跑完10组对比；或者当土壤参数不确定（如K(θ)函数有多个候选模型），先用显式扫一遍参数空间，快速定位合理区间。
隐式轨道（scheme = 'implicit'）：用于正式模拟与结果交付。它采用半隐式线性化策略：将方程右边的K(θ)在当前时间层n取值，而∂h/∂z中的h(θ)则用θⁿ⁺¹的线性近似。具体离散后得到形如A·θⁿ⁺¹ = b(θⁿ)的线性方程组，其中系数矩阵A只依赖于θⁿ，无需迭代即可直接求解。这比全隐式牛顿法简单，又比纯显式稳定得多——实测中，对典型壤土，Δt=3600s（1小时）仍能保持误差<2%。

提示：脚本中if scheme == 'implicit'分支下的A矩阵构建是核心。它不是简单的三对角阵，因为K(θ)的非线性导致A的非对角元包含∂K/∂θ项，但本脚本用K_half的调和平均替代了该导数计算，既避免求导误差，又保持矩阵对称正定——这是多年调试后找到的平衡点。

2.3 边界条件的物理实现与数值处理

边界条件不是数学符号，而是物理过程的翻译。脚本支持三类边界，每种都对应真实场景：

定水头边界（bc_top = 'head',h_top = -0.1）：模拟地表积水10cm，此时h(z=0,t) = -0.1m。数值上，它直接固定顶层节点的h值，再通过h(θ)反解出θ，作为Dirichlet条件代入方程。
给定通量边界（bc_top = 'flux',q_top = -5e-7）：模拟降雨入渗（负号表示向下），此时q(z=0,t) = -5×10⁻⁷ m/s。数值上，它转化为顶层界面的K·(∂h/∂z + 1)值，构成Neumann条件，需在离散时外推至边界。
自由排水边界（bc_top = 'free_drainage'）：模拟无积水时的自然蒸发/入渗，此时∂h/∂z = 0（即h在表层梯度为零）。这等价于q = K(θ)，意味着表层通量完全由当地导水率决定——代码中通过设置K_half(1) = K(theta(1))实现，物理意义是：水分只靠自身导水能力排出，不受外部压力驱动。

注意：底部边界统一设为free_drainage，但实际应用中若模拟深层地下水顶托，可轻松扩展为head类型。脚本预留了bc_bottom变量，只是默认未启用——这种“够用就好，留有余地”的设计，正是工程脚本与学术代码的本质区别。

3. 关键参数配置与土壤水力函数实现

3.1 参数集中配置区：50行内掌控全局

打开solve_one_dimension_Richards.m，前50行就是你的“控制台”。这里没有XML、没有JSON、没有.mat加载，所有参数以MATLAB变量直白呈现：

% === 网格与时间设置 === L = 1.0; % 土壤柱总长度 (m)，正方向向下 Nz = 101; % 空间节点数（含上下边界） dt = 3600; % 时间步长 (s)，隐式下可大胆设为3600 nt = 168; % 总时间步数（7天） % === 初始与边界条件 === theta0 = 0.25 * ones(Nz,1); % 初始含水率均匀分布 bc_top = 'head'; h_top = -0.05; % 上边界：地表积水5cm bc_bottom = 'free_drainage'; % 下边界：自由排水 % === 土壤参数（van Genuchten模型） === alpha_vg = 1.5; % (1/m) 形状参数 n_vg = 2.0; % 无量纲 形状参数 theta_s = 0.42; % 饱和含水率 theta_r = 0.06; % 残余含水率 Ks = 1.2e-5; % (m/s) 饱和导水率

这种设计消灭了“配置地狱”。你想换土壤？改alpha_vg、n_vg、Ks四行；想模拟干旱？把theta0改成0.12*ones(Nz,1)；想看暴雨响应？把bc_top换成'flux'并设q_top = -1e-6。所有改动实时生效，无需重启、无需编译、无需理解路径规则。

3.2 van Genuchten模型的MATLAB向量化实现

土壤水力函数是Richards方程的“心脏”。脚本采用最经典的van Genuchten (1980) 模型，但实现上做了关键优化：

水分特征曲线h(θ)：
matlab h_vg = @(theta) -1./alpha_vg .* ((theta-theta_r)./(theta_s-theta_r)).^(-1./(1-1./n_vg)) ... ./ (1 - ((theta-theta_r)./(theta_s-theta_r)).^(n_vg./(n_vg-1))).^(1./n_vg);
这里用向量化匿名函数，输入theta向量，直接输出h向量。重点在分母的(1 - ...)项——它确保当θ→θₛ时h→0，当θ→θᵣ时h→-∞，完美复现物理极限。
非饱和导水率K(θ)（Mualem模型耦合）：
matlab Se = @(theta) ((theta-theta_r)./(theta_s-theta_r)).^(1./n_vg); % 有效饱和度 K_vg = @(theta) Ks .* Se(theta).^0.5 .* (1 - (1 - Se(theta).^(n_vg./(n_vg-1))).^((n_vg-1)./n_vg)).^2;
Se是有效饱和度，K_vg中Se^0.5来自Mualem假设，后半段是van Genuchten的湿润锋修正。整个函数在θ=θᵣ处K=0，在θ=θₛ处K=Kₛ，且全程光滑可导——这对隐式格式的数值稳定性至关重要。

实操心得：我在Hor_SL_L_0.5m.txt里放的实测数据，theta从0.07到0.41，h从-10⁵Pa到-10Pa。用van Genuchten拟合时，alpha_vg和n_vg必须协同调整：alpha_vg大则干燥端更陡，n_vg大则过渡带更宽。脚本默认值alpha_vg=1.5, n_vg=2.0对应典型粉壤土；若你用砂土，建议alpha_vg=5.0, n_vg=1.8——这个经验值，是我对比27组实测数据后总结的。

3.3 空间离散与网格设计的物理考量

空间网格不是越密越好。脚本采用等距网格（dz = L/(Nz-1)），原因很实在：非饱和带水分运动的前沿（wetting front）通常集中在表层0~30cm，此处梯度最大。但若为此在表层加密网格，会导致dz不一致，K_half计算复杂化，且对教学演示不友好。

我的折中方案是：用足够细的等距网格捕捉前沿，再通过theta0的初始分布引导数值聚焦。例如，设Nz=101（dz=0.01m），虽全柱1m都用1cm步长，但初始theta0设为[0.12*ones(30,1); 0.25*ones(71,1)]，即表层30cm干燥、下层湿润。这样，数值解自然在干燥-湿润交界处产生陡峭梯度，无需手动网格加密。实测表明，对入渗模拟，dz=1cm已足够分辨毫米级的含水率变化；若研究蒸发，可将Nz提到201（dz=0.5cm），耗时仅增40%，精度提升显著。

4. 核心算法实现与差分步骤详解

4.1 显式格式：一步到位的物理直觉

显式格式的核心思想是：用已知的θⁿ时刻所有信息，直接算出θⁿ⁺¹。脚本中关键离散步骤如下（以内部节点i=2:Nz-1为例）：

% 计算节点i与i+1间界面的调和平均导水率 K_half(i) = 2 ./ (1./K_vg(theta(i)) + 1./K_vg(theta(i+1))); % 计算节点i处的压力水头（用h_vg反解） h_i = h_vg(theta(i)); % 构建通量差分：q_{i-1/2} - q_{i+1/2} % 其中 q_{i+1/2} = K_half(i) * ( (h_{i+1}-h_i)/dz + 1 ) flux_diff(i) = (K_half(i-1) * ((h_i - h_prev(i-1))/dz + 1) ... - K_half(i) * ((h_next(i+1) - h_i)/dz + 1));

这里h_prev和h_next是上一时间步和下一时间步的h向量。注意K_half用调和平均而非算术平均，这是达西定律在异质介质界面的正确处理方式——它保证通量连续，避免因K突变产生的虚假源汇。

显式更新公式为：

theta_new(i) = theta(i) + (dt/dz) * flux_diff(i);

即∂θ/∂t ≈ (q_{i-1/2} - q_{i+1/2})/Δz。整个过程无需矩阵运算，纯向量化循环，Nz=101时单步耗时<0.5ms。

踩过的坑：早期版本用算术平均K_half = 0.5*(K1+K2)，结果在干燥-湿润交界处出现θ振荡（数值色散）。换成调和平均后，振荡消失，且K_half在θ→θᵣ时自动趋近于0，物理意义更坚实。

4.2 隐式格式：稳定性的代价与收益

隐式格式将θⁿ⁺¹作为未知量，方程变为：
$$
\frac{\theta_i^{n+1} - \theta_i^n}{\Delta t} = \frac{1}{\Delta z} \left[ K_{i-1/2}^{n+1} \left( \frac{h_{i}^{n+1} - h_{i-1}^{n+1}}{\Delta z} + 1 \right) - K_{i+1/2}^{n+1} \left( \frac{h_{i+1}^{n+1} - h_i^{n+1}}{\Delta z} + 1 \right) \right]
$$

直接求解此非线性方程组极难。脚本采用半隐式线性化：K_{i±1/2}仍用θⁿ计算（即K_half基于θⁿ），而h^{n+1}用θ^{n+1}的线性近似h_i^{n+1} ≈ h_i^n + (∂h/∂θ)_i^n · (θ_i^{n+1} - θ_i^n)。经整理，得到线性系统A·θ^{n+1} = b。

A矩阵构建代码精炼如下：

% 对每个内部节点i (2:Nz-1) dhdtheta_i = dh_dtheta_vg(theta(i)); % h对θ的导数，预计算 A(i,i-1) = -K_half(i-1) / dz^2 * dhdtheta_i; A(i,i) = 1/dt + (K_half(i-1)+K_half(i))/dz^2 * dhdtheta_i; A(i,i+1) = -K_half(i) / dz^2 * dhdtheta_i; b(i) = theta(i)/dt + (K_half(i-1)-K_half(i))/dz * dhdtheta_i ... + K_half(i-1)/dz^2 * (h(i)-h(i-1)) + K_half(i)/dz^2 * (h(i)-h(i+1));

dh_dtheta_vg是h_vg的解析导数，避免数值微分误差。A是三对角矩阵，用MATLAB的\运算符高效求解。

实操心得：dhdtheta_i在θ接近θᵣ时极大（因h→-∞），这会使A(i,i)主导，抑制θ剧烈变化——这正是物理上“干燥土壤导水率极低，含水率难以下降”的体现。若忽略此项，隐式解会过度平滑前沿，失真严重。

4.3 边界节点的特殊处理

边界是数值误差高发区。脚本对三类上边界分别编码：

定水头边界（bc_top = 'head'）：
matlab % 固定h(1) = h_top，反解theta(1) theta(1) = fzero(@(th) h_vg(th) - h_top, [theta_r, theta_s]); % 在A矩阵中，将第一行设为[1,0,0,...]，b(1)=theta(1)，强制θ₁不变
给定通量边界（bc_top = 'flux'）：
matlab % 通量q_top = K_half(1) * ( (h(2)-h(1))/dz + 1 )，解出h(1) h1_est = h(2) + dz * (q_top/K_half(1) - 1); % 再用h1_est反解theta(1) theta(1) = fzero(@(th) h_vg(th) - h1_est, [theta_r, theta_s]);
自由排水边界（bc_top = 'free_drainage'）：
matlab % 设∂h/∂z = 0 at z=0，即h(1) = h(2)，故K_half(1) = K(theta(1)) % 通量q_top = K(theta(1))，直接赋值

底部边界同理，但仅实现free_drainage（h(end) = h(end-1)）。这种边界处理不追求数学完美，但确保通量守恒——我用质量平衡检查过：sum(theta_new - theta)*dz（总水量变化）与边界净通量积分误差<0.1%，满足工程精度。

5. 输出结果解析与可视化实践

5.1 结果数据结构：为后续分析而生的设计

脚本输出theta_history是一个Nz × nt的矩阵，theta_history(i,j)表示第j个时间步、第i个空间节点（深度z(i)）的体积含水率。这种列优先存储（MATLAB默认）便于绘图：

% 绘制含水率时空演化热图 figure; imagesc(t_vec, z_vec, theta_history'); colorbar; xlabel('Time (h)'); ylabel('Depth (m)'); title('Soil Water Content Evolution'); % 绘制特定深度的含水率时间序列 depth_idx = find(z_vec >= 0.3 & z_vec <= 0.31, 1); % 30cm深处 plot(t_vec/3600, theta_history(depth_idx,:)); xlabel('Time (h)'); ylabel('\theta (m^3/m^3)');

z_vec由linspace(0,L,Nz)生成（z=0为地表），t_vec = (0:nt-1)*dt。所有单位统一为国际单位制（m, s），避免单位混淆——这是我见过最多新手错误：把Ks输成cm/s却忘了转换。

配套的Hor_SL_L_0.5m.txt是真实数据，格式为两列：第一列h_cm（压力水头，cm），第二列theta（含水率）。脚本内置读取函数：

data = load('Hor_SL_L_0.5m.txt'); h_cm = data(:,1); theta_exp = data(:,2); h_exp = h_cm / 100; % 转为米

你可以用lsqcurvefit将h_vg拟合到h_exp和theta_exp，得到专属土壤参数——这才是真正“属于你地块”的模拟起点。

5.2 典型可视化案例与物理洞察

运行脚本默认参数（theta0=0.25,bc_top='head',h_top=-0.05），你会得到三类关键图表：

时空热图（上图）：清晰显示“湿润锋”如何随时间下移。前2小时，锋面移动缓慢（毛管吸力主导）；2~12小时加速（重力作用增强）；12小时后渐趋平缓（基质势梯度减小）。在z=0.4m处，θ从0.25升至0.38仅用8小时，印证了砂壤土的快速入渗特性。
深度剖面图（t=24h）：plot(z_vec, theta_history(:,end))显示典型的“S型”曲线——表层饱和（θ≈0.42），中部陡降，深层渐近于初始值。拐点深度即湿润锋位置，可用find(theta_history(:,end) > 0.35, 1, 'first')快速提取。
通量时间序列：由q_top计算得出（若为head边界，则q_top = K_half(1)*( (h(2)-h(1))/dz + 1 )）。你会发现：初始通量高达1e-5 m/s（10mm/h），1小时后降至2e-6，24小时后仅1e-7——这正是“入渗速率递减律”的数值再现。

小技巧：想快速比较不同土壤？复制脚本，改名solve_sand.m和solve_clay.m，分别设Ks=1e-4（砂土）和Ks=1e-8（黏土），运行后用hold on叠绘三条θ-t曲线。你会发现，黏土的湿润锋24小时仅下移5cm，而砂土已达80cm——这种直观对比，比读10页论文更深刻。

6. 常见问题排查与实操避坑指南

6.1 数值发散：90%的问题出在这里

现象：运行几秒后，theta出现负值或远超θₛ（如θ>0.5），程序报错Matrix is singular或Inf/NaN。

根因与对策：
-时间步长过大（显式）：dt > dz²/(2·Ks)。对策：降低dt，或切换至implicit格式。
-初始含水率过高（接近θₛ）：h_vg在θ→θₛ时导数趋零，dh_dtheta计算失效。对策：检查theta0是否≤0.95·θₛ；若必须模拟饱和，改用h为因变量的Richards方程。
-土壤参数不合理：Ks设为1e-3（粗砂）却用n_vg=1.2（黏土参数），导致K(θ)在中间含水率区异常高。对策：用Hor_SL_L_0.5m.txt数据拟合，或查Carsel & Parrish (1988) 表格选参。

快速诊断：在脚本中加入fprintf('Step %d: max(theta)=%.4f, min(theta)=%.4f\n', n, max(theta), min(theta));，观察发散前θ的范围。

6.2 结果失真：看似收敛，实则错误

现象：theta始终在0.25附近小幅波动，湿润锋不动；或通量曲线平滑无衰减。

根因与对策：
-K_half计算错误：误用算术平均或漏掉+1（重力项）。对策：检查K_half公式，打印K_half(1:5)，确认其随θ增大而增大。
-边界条件误配：bc_top='head'但h_top设为-100（相当于-10m水头，远超土壤持水能力），导致theta(1)被强制为θₛ，但下层θ未响应。对策：h_top应介于h(θ₀)与0之间，例如θ₀=0.25时，h_vg(0.25)≈-100cm，则h_top设-50cm合理。
-dz过粗：Nz=21（dz=0.05m）无法分辨1cm级的含水率梯度。对策：增至Nz=101，观察theta_history在交界处的梯度是否变陡。

6.3 效率瓶颈：如何让长时模拟不卡死

现象：模拟7天（nt=168）耗时>10分钟。

优化手段：
-向量化替代循环：脚本已全部向量化，但检查是否有for i=1:Nz内嵌fzero——fzero在循环中调用极慢。对策：对theta0和边界反解，用arrayfun(@fzero, ...)批量处理。
-预分配内存：确保theta_history = zeros(Nz, nt)在循环前声明，避免动态增长。
-减少绘图频率：默认每步绘图会拖慢百倍。对策：注释掉plot语句，或改为if mod(n,10)==0每10步绘一次。

终极提速：将核心循环编译为MEX函数。我已用C写好richards_solver.c，MATLAB调用mex richards_solver.c后，速度提升8倍——需要者可邮件索取。

6.4 教学演示锦囊：让学生一眼看懂的技巧

“冻结”时间步：设nt=1，dt=3600，运行后用disp(theta)查看单步变化，学生能亲手算出θ₂如何由θ₁和通量差推出。
制造“故障”演示：故意将K_half设为常数Ks，运行后展示平坦的θ剖面——对比真实K_half的陡峭前沿，凸显非线性的重要性。
参数滑块互动：用MATLAB App Designer封装，添加alpha_vg、n_vg滑块，实时刷新热图。学生拖动时，亲眼看到alpha_vg增大如何让干燥端更“硬”。