📚 课程目录第3章土壤-作物-大气连续体

3. 土壤-作物-大气连续体：水分平衡方程、土壤水分胁迫系数与根系吸水模型

各位同学，欢迎来到第三章。这一章我们聊点实在的——水。

搞作物模型，说白了就是在模拟“土壤-作物-大气”这个连续体里，水是怎么跑的。我刚开始做这个方向时，总觉得不就是水从土里到根里再到天上吗？后来发现，这里面的坑，一个比一个深。

今天我们就拆开来看三个核心模块：水分平衡方程、土壤水分胁迫系数、以及根系吸水模型（Feddes模型）。这三块是连在一起的，你搞懂了它们，作物模型的水分模块基本就拿下了。

3.1 水分平衡方程：水账本怎么记？

水分平衡方程，说白了就是给土壤水分记一本流水账。收入多少，支出多少，最后剩多少。

我习惯用这个形式：

ΔS = P + I - ET - D - R

其中：

ΔS：土壤储水量的变化量（mm）
P：降水量（mm）
I：灌溉量（mm）
ET：蒸散量（mm）——这是大头
D：深层渗漏量（mm）
R：地表径流量（mm）

你想想看，这个方程其实特别像你记账。收入是P和I，支出是ET、D和R。如果收入大于支出，ΔS就是正的，土壤变湿；反之变干。

我的经验：在实际项目中，D（深层渗漏）是最容易被忽略的。我曾经在一个水稻田项目里，只算了P和ET，结果模拟出来的土壤含水量一直偏高。后来一查，原来是渗漏没算进去。水稻田的渗漏量可不小，一天能跑掉2-5mm。

水分平衡方程的时间步长通常取日尺度。为什么？因为气象数据大多是日值，而且作物生长过程在日尺度上已经能描述得比较好了。

3.2 土壤水分胁迫系数：作物到底渴不渴？

有了水分平衡方程，我们知道了土壤里有多少水。但问题来了：土壤有水，作物就一定能吸到吗？

答案是否定的。土壤水分太低，根系吸不动；水分太高，根系缺氧也吸不动。这就引出了土壤水分胁迫系数的概念。

我一般用Ks来表示这个系数，取值范围0到1：

Ks = 1：水分充足，无胁迫
Ks = 0：完全胁迫，作物停止吸水
0 < Ks < 1：部分胁迫，吸水受限

怎么算Ks？最经典的是分段线性函数：

如果 θ ≥ θfc：Ks = 1
如果 θpwp < θ < θfc：Ks = (θ - θpwp) / (θfc - θpwp)
如果 θ ≤ θpwp：Ks = 0

其中：

θfc：田间持水量（Field Capacity）
θpwp：永久萎蔫点（Permanent Wilting Point）
θ：当前土壤含水量

关键点：这个线性模型虽然简单，但很多商业软件（比如DSSAT、AquaCrop）都在用。为什么？因为它抓住了主要矛盾——土壤水分在有效范围内时，胁迫程度与含水量基本成正比。

不过，我遇到过一个问题：有些作物对水分胁迫特别敏感，比如玉米在抽穗期。这时候线性模型就不太准了。我建议在关键生育期，可以改用非线性函数，比如指数型：

Ks = (θ - θpwp) / (θfc - θpwp) ^ 0.5

这样胁迫系数下降得更快，更符合实际情况。

3.3 根系吸水模型：Feddes模型

好了，现在我们知道土壤里有多少水，也知道作物受不受胁迫。但根系到底怎么吸水？这就需要一个根系吸水模型。

Feddes模型是我用得最多的，也是目前作物模型里最主流的。它把根系吸水看作一个“源汇项”，直接加到水分平衡方程里。

Feddes模型的核心公式：

S(h) = α(h) × Smax

其中：

S(h)：实际根系吸水速率（cm³/cm³/day）
Smax：最大潜在根系吸水速率
α(h)：水分胁迫响应函数（0到1之间）

这个α(h)就是Feddes模型的精髓。它把土壤水势（h）和根系吸水能力关联起来。我画了一张图，你看一眼就明白了：

这张图怎么看？横轴是土壤水势（h），越往右越湿；纵轴是α(h)，越大吸水越猛。

h₄到h₃：太湿了，根系缺氧，吸水能力为0
h₃到h₂：水分逐渐变好，吸水能力线性上升
h₂到h₁：最优区间，吸水能力拉满（α=1）
h₁到0：开始干旱，吸水能力线性下降

避坑指南：我曾经在参数化Feddes模型时，直接把h₁、h₂、h₃、h₄从文献里抄过来用。结果模拟出来的蒸散量总是偏低。后来发现，不同土壤类型的这些阈值差异很大。砂土的h₁可能只有-100 cm，而黏土能到-500 cm。一定要根据实际土壤质地来调。

3.4 三者的集成：从方程到模拟

现在我们把三个模块串起来。我习惯的流程是这样的：

第一步：用水分平衡方程算出每天的ΔS，更新土壤含水量θ
第二步：根据θ计算土壤水分胁迫系数Ks
第三步：用Feddes模型计算实际根系吸水速率S(h)
第四步：把S(h)作为源项，反馈到水分平衡方程中，更新下一时间步的θ

你看，这是一个闭环。每一步的输出都是下一步的输入。

我给大家一个简化的Python实现思路：

# 伪代码：水分平衡 + 胁迫系数 + Feddes模型
def daily_water_balance(P, I, ET_pot, theta_prev, theta_fc, theta_pwp, h1, h2, h3, h4):
    # 1. 计算水分胁迫系数
    if theta_prev >= theta_fc:
        Ks = 1.0
    elif theta_prev <= theta_pwp:
        Ks = 0.0
    else:
        Ks = (theta_prev - theta_pwp) / (theta_fc - theta_pwp)
    
    # 2. 计算实际蒸散
    ET_act = ET_pot * Ks
    
    # 3. Feddes模型计算根系吸水
    # 假设当前土壤水势h由theta估算
    h = estimate_h_from_theta(theta_prev)
    
    if h <= h4 or h >= 0:
        alpha = 0.0
    elif h >= h2 and h <= h1:
        alpha = 1.0
    elif h > h3 and h < h2:
        alpha = (h - h3) / (h2 - h3)
    else:  # h > h1
        alpha = (0 - h) / (0 - h1)
    
    S_act = Smax * alpha  # Smax需要根据根系分布计算
    
    # 4. 更新土壤水分
    D = calculate_drainage(theta_prev)  # 渗漏计算
    R = calculate_runoff(P)             # 径流计算
    
    delta_S = P + I - ET_act - D - R
    theta_new = theta_prev + delta_S / soil_depth
    
    return theta_new, Ks, S_act

我的建议：实际项目中，Smax（最大根系吸水速率）不要设成一个常数。我习惯把它和根长密度（Root Length Density）挂钩。根多的地方，吸水潜力自然大。你可以用指数衰减函数来描述根系随深度的分布。

3.5 参数获取与校准

模型搭好了，参数从哪来？我列个表，方便你对照：

参数	含义	获取方式	我的经验值
θfc	田间持水量	土壤测定或文献	砂土0.15-0.25，黏土0.35-0.45
θpwp	永久萎蔫点	土壤测定或文献	砂土0.05-0.10，黏土0.15-0.25
h₁	干旱胁迫点	Feddes模型查表	-500到-1000 cm
h₂	最优吸水上限	Feddes模型查表	-100到-300 cm
h₃	最优吸水下限	Feddes模型查表	-300到-600 cm
h₄	厌氧点	Feddes模型查表	-50到-100 cm

校准技巧：我一般用实测的土壤含水量时间序列来反推参数。先固定θfc和θpwp（这两个容易测），然后调Feddes的四个阈值，让模拟的土壤含水量和实测的RMSE最小。别一次性调所有参数，会过拟合。

好了，这一章的内容就到这里。水分平衡方程是骨架，胁迫系数是开关，Feddes模型是执行器。三样东西凑在一起，你就能模拟出作物每天到底喝了多少水。