1 简介
计算地球物理学的一个中心问题是如何求解以下形式的守恒方程
∂φ∂t= -1VS∂F∂z
(1)其中φ是感兴趣的状态变量，F是 z 方向（此处视为垂直）的通量，t是时间，C是常数。例如，如果φ是温度，C是体积热容，则 (1) 是众所周知的热方程（能量守恒定律），热通量由下式给出
F( z）=−K _∂φ∂z
(2)其中K是热导率。或者，如果液体）土壤湿度且C = 1，则 (1) 是水守恒定律。对于不饱和流，这是理查兹（引用1931）方程，可以写出水分通量（例如Verseghy，引用1991）作为
F( z) = λ ( φ ) [−bψ ( φ ) _ _φ∂φ∂z+ 1 ]
(3)式中，λ为导水率，ψ为土壤吸水力，b为土壤质地参数。在这两种情况下，（1）都可以使用简单的显式（或隐式）有限差分方案进行数值求解。一种常见的方法是采用交错垂直网格，其中在层界面处计算通量（标记为i = 0,1,2,3 …），并在层中点处计算状态变量（标记为i = 1/2, 3/2, 5/2 ，……）。因此 (1) 可以明确地离散为
φd + 1我+ 1 / 2=φj我+ 1 / 2-1VS我+ 1 / 2[Fj我+1 _-Fj我]Δt_ _Δz_ _
(4)其中j是时间索引。给定适当的层厚度Δz和时间步长Δt以及 (2) 或 (3) 的有限差分形式，这可以在均匀网格上轻松解决。当然，如果包括所有适当的源项和汇项，例如横向通量、水相变化引起的加热、地表水渗透等，则（1）和（4）只是一个守恒方程，尽管这些可以很容易地表示为作为附加条款添加。本文解决的问题是如何最好地离散方程中的F。（3）。值得注意的是（3）中的导水率和土壤吸水力都是高度非线性的，并且要计算通量F土壤湿度 φ 及其梯度必须在层界面处进行估计。通过使用足够密集、均匀的垂直网格（如果需要，还可以使用隐式数值方法）可以在一定程度上改善这些问题。

然而，许多陆地表面模型（尤其是用作数值天气预报或气候建模的组成部分时）采用相对稀疏且通常高度不规则的垂直网格。问题是计算效率和准确性之间的平衡之一。一般来说，大部分力发生的地表附近的土层很薄，但在更深的土层中可以变得更厚，以便在没有过多层数的情况下模拟深层土柱。有时，层厚度随深度的增加是系统性的，例如几何增加（例如 de Rosnay 等人，引用2000年,引用2002年；海斯等人，引用2003）。然而，层间距通常非常不规则，并且在文献中很少合理化（表格1）。

加拿大陆地表面计划（类别：Verseghy，引用1991年,引用2017）和土壤、植被和雪模型（SVS：Alavi 等人，引用2016年；侯赛因等人，引用2016）还采用了不均匀的垂直网格。为了在这些模型中离散化 (3)，必须在层界面处估计土壤湿度及其梯度，并且这两个模型都遵循 Verseghy 中使用的方法（引用1991），概述如下。我们在这里认为，这种方法在不规则网格上是有问题的——特别是当层厚度突然增加时——并且提出了一种替代方法，可以系统地离散任意、交错的一维网格上的一维和更高阶的通量。讨论了比较这两种方法的数次多年模拟的结果。

2 离散化Richards方程
CLASS 和 SVS 中的传统方法
如图。1(a) 说明了 CLASS 和 SVS 用于土壤湿度的非均匀交错网格类型的一部分。在此描述中，z i表示第i层底部的深度，因此z i – z i-1是第 i层的厚度。对于 CLASS 的传统应用，包括加拿大地球系统模型（CanESM5：Swart 等人，引用2019），这些水平是： z 1  = 0.1 m、 z 2  = 0.35 m 和z 3  = 4.1 m，但其他配置也是可能的，包括许多规则或不规则间隔的水平（例如 Elshamy 等人，引用2020 年；梅尔顿等人，引用2019）。SVS 的典型应用使用 7 个级别。例如 Maheu 等人。（引用2018）和达布尔等人。（引用2019）使用 0.05、0.10、0.20、0.40、1.0、2.0 和 3.0 m。土壤湿度φ代表平均层值，被认为发生在层中点。因此，对于第 i层，土壤湿度φ i-1/2为z i-1/2  = ( z i + z i-1 )/2。任意土壤湿度剖面也如图所示如图。1(a)（实心蓝色圆圈）。第 i层和第 i+1层之间的水分通量F i由z i处的土壤湿度梯度决定，但也由z i处土壤湿度本身的值（以及其他变量，如电导率和吸力）决定，必须对其进行插值。CLASS 和 SVS 都基于简单平均值 将土壤湿度插值到z i ：
φ我= 1 / 2 (φ我+ 1 / 2+φ我- 1 / 2）
(5)

（如图。1(a) 红色空心圆圈），这样z i处的湿度梯度就成为斜率的平均值：
∂φ∂z|z我=φ我+ 1 / 2-φ我z我+1 _-z我+φ我-φ我- 1 / 2z我-z我-1 _
(6)图示于如图。1(a)（红实线）。在均匀网格上，这表示中心一阶有限差分。然而，在不均匀的网格上，较薄层中的梯度权重大于较厚层中的梯度。韦尔塞吉 (引用2017）认识到这个问题，但表明这模拟了指数衰减的土壤湿度剖面。在本研究中，我们将此插值称为传统方法或 V91 方法。下面我们将证明这种方法可能会导致不必要的数值错误。

b 拉格朗日插值多项式
作为替代方案，我们可以对φ和∂φ / ∂z
到z i基于Φ附近值的更广泛线性组合（例如 Φ i-1/2、Φ i+1/2、Φ i-3/2等）。实现这一点的系统方法是利用拉格朗日插值多项式。回想一下，对于N点插值，拉格朗日多项式L(z)是通过每个N （任意间隔）点的唯一N-1阶多项式。这可以根据拉格朗日基多项式l(z)的线性组合轻松构建，如下所述。

N =  2：线性插值

对于N  = 2，基于 2 个最近点（φ i-1/2，φ i+1/2 ）的z i处土壤湿度的拉格朗日插值可以写成
Li(z)=∑j=i−1/2i+1/2φjlj(z)=φi−1/2li−1/2(z)+φi+1/2li+1/2(z)
(7)其中
li−1/2(z)=(z−zi+1/2)(zi−1/2−zi+1/2)(8a)
(8a)
li+1/2(z)=(z−zi−1/2)(zi+1/2−zi−1/2)(8b)
(8b)

很容易看出，L i ( z ) 是一个连续的线性函数，它同时经过 ( z i-1/2 , φ i-1/2 ) 和 ( z i+1/2 , φ i+1/ 2），即
Li(zi+1/2)=φi+1/2
和
Li(zi−1/2)=φi−1/2

结合(7)和(8)我们发现
Li(z)=1zi+1/2−zi−1/2[(z−zi−1/2)φi+1/2−(z−zi+1/2)φi−1/2]
(9)

为了估计z i处的土壤湿度，我们只需评估
φi≈Li(zi)=1zi+1/2−zi−1/2[(zi−zi−1/2)φi+1/2−(zi−zi+1/2)φi−1/2]
(10)

这在如图。1(b)（橙色空心圆圈）。请注意，在层厚 Δz 的均匀网格上，这变为
φi≈Li(zi)=1Δz[(Δz2)φi+1/2−(−Δz2)φi−1/2]=1/2(φi+1/2+φi−1/2)
这是 CLASS 和 SVS 中使用的常规表达式，尽管在它们的默认配置中都不使用统一网格。

评估∂φ/∂z
在z i处，我们对 z 求导 (9) 并评估z i处的结果，得到
∂φ∂z|zi≈dLidz|zi=1zi+1/2−zi−1/2[φi+1/2−φi−1/2]
(11)这是橙色曲线的斜率如图。1（二）。这是层界面土壤湿度梯度最简单的表达方式，并用于求解理查兹方程（基于湿度的形式）的其他地表模型（例如de Rosnay 等人，引用2002年；牛等人，引用2011）。有趣的是，对于上层（在水文上）比下层厚的特殊情况，CLA​​SS 确实使用了等价于（10）和（11）的表达式，当下层存在基岩而使其渗透深度较薄时，就会发生这种情况。韦尔塞吉 (引用2017）表明在这种情况下（5）不再合适。当一阶精度足够时，我们认为没有理由不将 (10) 用于所有情况。

例子

一个具体的例子将有助于阐明方法并巩固概念。考虑CLASS 的标准3 级配置（例如，在CanESM5 中实现），其中z 1  = 0.1 m、z 2  = 0.35 m 和z 3  = 4.1 m。我们感兴趣的是新的离散化对F 2（第二层和第三层之间的通量）的影响。我们可以预见，鉴于第三层比第二层厚一个数量级：3.75 m 与 0.25 m，上述方法可能会产生显着不同的结果。我们将特别关注湿度梯度

因为这是（3）的重要组成部分，将受到离散化的强烈影响。

3 模拟
为了展示新离散化的影响，开源继承了 CLASS（CLASSIC：Melton 等人，引用2020）已配置为对五个对比 FLUXNET 进行多年模拟（Pastorello 等人，引用2020）网站。地理学详细信息见表2。对于每个站点，通过每半小时的气象测量（传入的短波和长波辐射、气温、降水率、气压、比湿度和风速）进行四次不同的模拟。模拟持续时间受到数据可用性的限制，根据地点的不同，模拟持续时间从 3 年到 11 年不等（表2）。CLASSIC 通过气象强迫循环而旋转，直到碳库达到平衡状态（例如Melton 等人，引用2020）。两次模拟均使用了传统的 CLASS 3 土层配置；一个使用传统的V91插值（方程（5）和（6）），另一个使用建议的一阶拉格朗日（即线性）插值（方程（10）和（11））。下面分别将它们标记为 L3_V91 和 L3_线性。根据 CLASSIC 中使用的标准 20 个土层配置了另外两个模拟（请参阅表格1），分别标记为 L20_V91 和 L20_线性。

GF-Guy
法属圭亚那热带雨林Guyaflux遗址（GF-Guy；北纬5.28°，西经52.92°）拥有深厚的沙质土壤，透水深度为37 m。年总降水量很大，但该地区的 8 月至 11 月有明显的气候旱季（Bonal 等人，引用2008）。模拟中植被表示为常绿阔叶树，覆盖率为 100%。

GF-Guy 2004-2014 年平均水平衡如下所示：表3。值得注意的是，所有 4 个模拟中水平衡的残差都很小：理查兹方程的新（线性）离散化和传统（V91）离散化都不违反水守恒定律（对于本研究中检查的所有五个地点都是如此）学习）。许多附加特征是显而易见的。在年度时间尺度上：

两个 20 层模拟在年度水平衡方面几乎相同

L3_线性水平衡与20层结果相似：ET比20层结果大约11毫米，基流小约11毫米，约占ET的1%。

另一方面，L3_V91的水平衡则显着不同。在这种情况下，20 层结果中的 ET 过量约为 69 毫米，由约 70 毫米的基流赤字（进入径流的差值）平衡。

表3表明 L3_V91 是这些实验中的异常值。通过检查年平均水平衡周期（如图。2）。如图。2(a)（重复于如图。2(e)) 显示了 2004-2014 年模拟期间的平均月降雨量。雨季从 12 月至 2 月以及 4 月至 7 月持续，主要由热带辐合带的运动驱动（Bonal 等人，引用2008）。此期间第一层土壤水分通量以入渗为主，模拟结果差异较小（如图。2(b))。然而，在旱季（8 月至 11 月），L3_V91 和 L3_线性之间可以看到显着差异，线性方案比传统方法更忠实地代表了观察到的第一层土壤湿度的下降。这对 ET 有直接影响，在此期间，L3_线性中的 ET 也减少了（如图。2（与））。20层模拟的相应结果如图所示如图。2（f，g）。如上所述，年度水平衡总量（表3）两个 20 层模拟对于第一层（还有更深的层，未显示）土壤湿度产生几乎相同的结果（如图。2(f)) 和 (如图。2（G））。这些模拟的土壤湿度和蒸散结果与 L3_线性的相应结果类似（如图。2(b,c) 分别) – 所有这些都与旱季 L3_V91 结果显着不同。这表明 L3_V91 确实是一个异常值，并表明如果限制在该具有深层土壤（渗透深度 37 m）的场地使用传统的 3 层 CLASS 土层，则使用线性插值在很大程度上会导致类似的结果到完整的20层模型。当然，例外的是基流，它受到 3 层配置中 4.1 m 土柱底部的强烈影响（如图。2(d,h))。虽然年度总流量相似，但较深的土柱往往会抑制基流年度循环（如图。2(h)) 而截断的土柱往往会促进闪光（如图。2(d))。这不是本文讨论的离散化方案的函数（V91 和线性曲线）如图。2(d) 类似），但传统 3 层配置中最大深度 4.1 m 的限制。

bCA-TPD
土耳其角成熟落叶林（CA-TPD；北纬 42.64°，西经 80.56°）位于加拿大安大略省，也有深厚的沙质土壤（表2），透水深度为39 m。植被由 95% 阔叶落叶冷乔木和 5% 针叶常绿乔木 (Arain & Restrepo-Coupe,引用2005）。

CA-TPD 2012-2014 年平均水平衡如下所示：表4。在年度时间尺度上：

同样，两个 20 层模拟在年度水平衡方面几乎相同

L3_线性水平衡与 20 层结果相似，但不如 GF-Guy 那样接近：ET 比 20 层结果大约 21 毫米，基流小约 21 毫米，约占 ET 的 6%。

L3_V91水平衡再次存在较大差异。在这种情况下，20 层结果中的 ET 过量约为 65 毫米，由约 64 毫米的基流赤字（进入径流的差值）平衡。

年平均水平衡周期见如图。3。该站点年降水量适中（460毫米），与之前站点相比，全年分布更加均匀（如图。3(a,e))。正如 GF-Guy 所指出的，两个 20 层模拟在水平衡方面产生了几乎相同的结果（如图。3(f,g))，L3_线性结果与这些类似(如图。3(b,c)) 虽然不如 GF_Guy 那么接近。在这种情况下，我们建议使用线性插值是对 3 层配置中 V91 的改进，尽管这种改进可能还不够，并且可能需要完整的 20 层配置。请注意，L3_线性和 L20 模拟 (21 毫米) 之间的 ET 差异约为年降水量的 5%，但 L3_V91 的差异为 14%，表明水平衡存在显着的不确定性。

cIT-Noe
意大利灌木区 Arca di Noe（IT-Noe；北纬 40.61°，东经 8.15°）位于撒丁岛西北部的 Le Progionette 自然保护区，覆盖率为 70% 的常绿阔叶灌木（表2）。该地区属半干旱地中海气候，夏季温暖，冬季温和。降水主要发生在 9 月至 5 月期间，6 月至 8 月降水很少（Marras 等，引用2011年；帕帕莱等人，引用2015）。该场地土壤非常浅，透水深度仅为0.4 m。这意味着在 3 层配置中，L3 在水文学上仅 5 厘米厚，比 L2 (25 厘米) 薄，因此 V91 插值方案仅对第一层和第二层之间的通量有效。此配置中的渗透层厚度为 10、25 和 5 厘米，这与 20 层配置中四个等距的 10 厘米（渗透）层形成对比。因此，我们预计四个模拟之间的差异将主要由分层差异主导，而插值造成的影响较小。这明显地体现在表5其中 L20 模拟彼此相似，L3 模拟彼此相似。层配置最重要的影响是对径流和基流之间水的分配：L3 流道产生约 25 毫米的径流和约 31 毫米的基流。尽管如此，L3_V91 和 L3_线性之间的差异是明显的，与前两个站点一致：L3_V91 中过多的 ET (8 毫米) 在很大程度上被基流中的不足 (5 毫米) 所平衡。

d MY-PSO
位于马来西亚中南部的热带雨林帕索森林保护区（MY-PSO；北纬 2.97°，东经 102.31°）由 100% 常绿阔叶树组成（表2; 小杉等人，引用2008）。该地区全年都有降水，三月和四月以及十月和十一月是降水高峰。该场地还有非常浅的土壤，透水深度为 1.0 m。因此，在 3 层配置中，水文活动层的厚度为 10、25 和 65 厘米，而在 20 层配置中，对 10 层相同厚度 (10 厘米) 进行建模。因此，所有层的 V91 插值都被激活（使用时），但渗透层厚度的跳跃永远不会超过 3 倍，与检查的前两个位置相反，在 3 层配置中，第三层是 15 倍比第二个厚。年平均水平衡结果列于表6。

在这种情况下，所有四个模拟都是可比较的，并且使用线性插值方案的影响很小。尽管如此，L3_V91 和 L3_线性之间的差异与之前的所有站点一致：L3_V91 中过多的 ET (3 mm) 由基流中的不足 (3 mm) 来平衡。

e US-Wkg
Walnut Gulch Kendall 草原场地（US-Wkg；北纬 31.74°，西经 109.94°）植被稀疏，C 4草覆盖率为 35%（表2; 斯科特等人，引用2010）。该地点非常干燥，大部分降水发生在 7 月至 9 月，没有明显的径流或基流，年平均降水量完全由蒸散平衡。年平均水平衡（表7）以及平均年水平衡周期（未显示）表明，四个模拟实际上是相同的：插值方案和分层配置对模拟水平衡结果没有任何影响。

4 讨论与结论
在交错垂直网格上计算基于湿度的理查兹方程的地表模型必须估计土壤湿度及其在层界面的梯度，这需要插值。在足够密集、均匀的网格上，这通常不会有问题，但许多模型使用稀疏且非常不规则的网格，需要更加小心。韦尔塞吉 (引用1991）将层界面处的（体积）土壤湿度插值为两层的算术平均值（方程 5），即使各层的厚度差异很大（如图。1(a)，空心红色圆圈）。然后将界面处的土壤湿度梯度估计为界面上方和下方梯度的平均值（方程 6；如图。1(a) 红实线）——优先加权较薄层中的梯度。我们在本研究中将这种方法标记为 V91，至今仍在许多加拿大建模系统中使用，包括 CLASSIC（Melton 等人，引用2020），SVS（Alavi 等人，引用2016），MESH（Pietroniro 等人，引用2007），CanESM（Swart 等人，引用2019），CanRCM（Scinocca 等人，引用2015）和CRCM（Martynov 等人，引用2013）。我们已经表明，与基于线性插值的方法相比，通过构造，这总是会夸大层界面处的土壤湿度梯度（当上层比下层薄时），即当界面土壤湿度是两层的加权平均值时（等式10，如图。1(b) 空心橙色圆圈），土壤湿度梯度是简单线性一阶有限差分（方程 11，如图。1(b) 橙色实线）。对于传统的 CLASS 3 层配置，我们发现第二层和第三层之间的湿度梯度大 4.27 倍（即当第三层完全渗透时）。当第二层比第三层干燥时，土壤水分扩散率（16）也过大，因此在这种情况下水分通量误差可能很大。在这些条件下，过多的吸力往往会增加第二层土壤湿度，其中更多的水分可用于蒸散，而较少的水分可用于基流。事实上，我们的模拟证实了这一点。

我们使用 CLASSIC 陆地表面模型检查了五个对比地点的模拟，该模型配置有 CLASS 通常使用的标准 3 层和 20 层，以及传统的 V91 方案 (5)、(6) 或线性方案 ( 10）、（11）。使用线性方案的重要性取决于土层配置、场地特征和气候。对于 20 层配置，土层厚度的跳跃很小（不超过 3 倍），使用线性插值对水平衡影响很小（尽管非零）。在对北美地区进行的一系列 SVS 模拟中也发现了这一普遍结论（Étienne Gaborit，个人通讯）。SVS 通常也使用厚度随深度增加相对缓慢的层（表格1）。这并不是说这里提出的问题对于 CLASSIC 和 SVS 建模研究不重要。在这两种模型中，用户可以自由（并且经常）指定他们想要的任何土壤层；很少考虑理查兹方程的离散化。

另一方面，对于 CLASS 的传统 3 层配置，从第二层到第三层的层厚度跳跃很大 - 15 倍，§2 中的分析表明这可能对土壤水分通量产生很大影响这些层之间，从而实现整体水平衡。我们的模拟表明，与 L20 模拟相比，V91 插值的影响一般是增加土壤湿度，从而增加蒸散量，但会牺牲基流（水平衡的其他组成部分相对不变）。在某些情况下，这些差异的幅度很小——特别是对于干燥地点或所检查的土壤浅层的地点。对于土壤较深的地点，差异要大得多。在加拿大站点（CA-TPD）蒸散量过多（即相对于 20 层模拟）是年降水量的 14%。在许多情况下，L3_线性模拟与相应的 L20 模拟非常接近，这表明即使受限于使用 CLASS 的传统 3 层配置，与 V91 相比，线性插值的使用也能改善水平衡模拟。

现代陆地表面模型代表了无数复杂的、通常约束很差的过程，并且很难将模拟错误归因于代码的任何特定方面。在本研究中，我们将注意力集中在一个问题上：理查兹方程的离散化。目前，土壤湿度及其梯度在 CLASS 和 SVS 中基于（5）和（6）进行离散化，我们在此建议由（10）和（11）给出更好的一阶精度离散化。如果需要，可以轻松且系统地构建更高阶的方案。

这对实验设计来说是一个潜在的强大帮助。然而，对这种可能性的进一步讨论超出了本研究的范围。一阶拉格朗日插值可能足以满足典型的水文学应用。