处理地下水模型单元疏干-湿润的两种算法对比研究

陆文， 陆垂裕， 何鑫， 孙青言， 张博， 贾仰文

摘要：在使用网格单元中心差分格式的地下水模型中，对地下水网格单元“疏干（干）-湿润（湿）转化”的模拟极易引发模型迭代不收敛等异常情况，很大程度上影响模型的应用。本研究使用理想案例和丹麦应用实例，综合比较了MODFLOW 模型的试算法与COMUS 模型的全有效单元法对网格单元“干-湿转化”问题的模拟能力及算法特性。结果表明：（1）试算法的参数组合选取对模拟的收敛性和模拟结果都有明显影响，使用试算法时需要不断优化参数组合以避免模型迭代不收敛或模拟失真等异常情况，很大程度上增加了用户使用模型的难度和时间成本；（2）全有效单元法的模拟结果比试算法的模拟结果更具可靠性，全有效单元法的作用等同于理论上最优的试算法参数组合，使用全有效单元法时用户无需进行复杂的参数组合工作，因此该方法能有效降低模型的使用难度与模拟结果的不确定性；（3）全有效单元法中单元间水平向水力传导度算法实现了可以与经典调和平均法相比较的数值计算精度，说明全有效单元法在不涉及网格单元“干-湿转化”问题的地下水模拟中同样具有应用潜力。综上所述，全有效单元法更适用于处理地下水模型单元的疏干-湿润转化问题，并且有望在地下水数值模拟领域中得到更为广泛的应用。

关键词：地下水模型；网格单元疏干-湿润；全有效单元法；试算法；饱和-非饱和

在类似MODFLOW 等采用网格单元中心差分格式模拟地下水饱和流动的数值模拟模型中（含水层上部可能耦合土壤水的饱和-非饱和流模拟，但地下水部分只模拟饱和流动），疏干单元指网格单元的模拟水头低于网格单元的底板高程，此时模拟水头只有计算意义而无水量意义，意味着网格单元中不存在重力水，由于不模拟非饱和流，这也意味着网格单元中没有任何水分存在；湿润单元指网格单元的模拟水头高于网格单元的底板高程，此时模拟水头具有水量意义，也即网格单元中存在重力水。网格单元一旦疏干，基于经典调和平均法计算的单元间水平向水力传导度为零，意味着地下水的侧向流动在该疏干单元处中断，这极有可能导致地下水运动差分方程组无法求解，即模型无法迭代收敛，使得模型模拟失败[1—2]。为保证模型的可用性，疏干单元将被转化为无效网格单元，等同于隔水边界，不再模拟其水头，但也意味着作用于该网格单元之上的源汇项全部丢失[3]。另一方面，当周边单元的水头高于疏干单元的底板高程时，又希望疏干单元能够重新湿润成为有效单元，模拟其水头和作用于其上的源汇项，以保证模拟的准确性，然而这一看似简单的问题在过去很长一段时间内都未能得到解决[4]。总之，模拟“干-湿转化”容易引发模型迭代不收敛，而不模拟则会引发模拟结果失真等异常情况，使地下水模型的应用和准确性受到很大限制[5—6]。

当前主要有2 种应对“干-湿转化”问题的思路。（1）类似HYDRUS-2D、VS2DT 模型完整模拟地下水的饱和-非饱和流动，此时不存在“干-湿转化”问题。这一思路缺陷在于饱和-非饱和渗流方程的求解效率较低，难以适用于区域尺度的现实复杂三维模拟问题。（2）仍然在饱和地下水流数值模拟模型中对“干-湿转化”问题进行考虑，目前这一思路发展出的方法总结起来主要有3 类。第1 类方法不具备物理意义，包括调大模型的收敛容许误差使得模型容易收敛；调低关键区域的底板高程[7]、调整模型的垂向分层结构[8]、为疏干单元赋以极小的饱和厚度[9—10]、将接近疏干的非承压域转换为承压域模拟等[11]，从而避免模拟域疏干。这类方法在一些情况下能够获得可以接受的模拟效果，但不可避免会引入模拟失真等问题[12]。第2 类方法具有一定的物理意义，以MODFLOW-2005采用的试算法（empirical trial method, ET）为代表。该方法根据无效疏干单元周边相邻湿润单元的水头是否达到预设湿润条件，反复试算确定疏干单元能否被重新湿润 [13]，然而试算法不能保证模拟结果的稳定性[4,12]。尽管如此，它仍然是当前最常用的“干-湿转化”模拟算法。第3 类方法具有坚实的物理机制，主要包含MODFLOW-NWT 模型提出的上游加权法[14]和笔者近期在基于网格单元中心差分格式自主研发的COMUS 模型[15]（C++ objected-oriented model for undergroundwater simulation）中提出的全有效单元法[16]（always active cell method，AAC）。这类方法能够从地下水动力学方程本身出发解决网格单元的“干-湿转化”模拟问题，主要特点在于模拟区域内的单元无论是疏干还是湿润，始终都是有效的单元，都考虑其水量平衡，建立其地下水运动差分方程，并纳入整体矩阵方程计算，网格单元“干-湿转化”能够无缝切换。然而由于上游加权法本身数值精度较低，因此目前MODFLOWNWT 用户较少。

本文通过理想案例和丹麦应用实例，比较传统的试算法和全有效单元法对“干-湿转化”问题的模拟能力和算法特性，检验全有效单元法的可靠性和实用性，并对全有效单元法的特点与适用性进行讨论。需要说明的是，COMUS 模型同样支持采用调和平均法计算单元间水平向水力传导度，也支持使用试算法模拟网格单元的“干-湿转化”过程；本研究中试算法和全有效单元法的模拟结果对比是基于COMUS 模型平台。

1 试算法原理

对于疏干单元，模型将在每次迭代过程中检查该单元下侧、前侧、后侧、左侧和右侧5 个相邻单元的水头情况，见图1，若发现某个相邻湿润单元的水头高于或等于湿润判断水头阈值，则在下次迭代时对该疏干单元进行重新湿润尝试[17]。湿润判断水头阈值（hTurnon）的计算公式如下：

图1 试算法概念图
Fig.1 Schematic diagram of the ET method

式中：Bot——疏干单元的底板高程/m；

WETDRY——用户指定的湿润特征参数。

WETDRY 的绝对值（|WETDRY|）表示湿润饱和厚度阈值/m。WETDRY＜0 表示若下侧湿润网格单元的水头值高于hTurnon，则认为疏干单元能在下次迭代时被重新湿润；WETDRY＞0 表示若5 个相邻单元中有任一湿润网格单元的水头值高于hTurnon，则认为疏干单元能在下次迭代时被重新湿润，通常WETDRY＜0 有助于模型收敛；WETDRY=0 表示疏干单元不可被重新湿润。为了避免刚被重新湿润的网格单元下次迭代又被疏干，还允许用户设置NWETIT 参数指定间隔多少次迭代才对疏干单元进行重新湿润尝试（例如NWETIT=2表示间隔两次迭代才对疏干单元进行重新湿润尝试），以提高成功率。

在判断疏干单元能够被重新湿润后，下次迭代时该疏干单元将被赋予高于其底板高程的试算水头值（hTry），从而再次成为有效网格单元参与模拟计算：

式中：WETFCT——用户指定的调整系数，一般取0～1之间的值[17]；

hn——驱动疏干单元重新湿润的相邻湿润网格单元的水头值/m；

IHDWET——指定试算水头值计算方式的标志参数。

IHDWET=1 表示试算水头值依据hn 计算，IHDWET=2 表示试算水头值依据湿润饱和厚度阈值计算。

疏干单元重新湿润后基于调和平均法计算的单元间水平向水力传导度将不为0，地下水的流动在该单元处恢复，该单元重新转变为有效网格单元，该单元的水头和作用于该单元之上的源汇项将再次被模拟。由试算法原理可知，试算法的参数组合对模拟结果有影响，并且理论上|WETDRY|越小，单元间地下水流动受到的人为限制越小，模拟结果越合理。

2 全有效单元法原理

全有效单元法中所有网格单元在整个模拟期内无论其是否疏干都是有效网格单元，都需要建立单元的地下水运动差分方程，模拟水头和源汇项。对任一网格单元（i, j, k）地下水运动差分方程为：

其中，下标（i, j-1, k）、（i, j+1, k）、（i-1, j, k）、（i+1, j, k）、（i, j, k-1）、（i, j, k+1）分别代表网格单元（i, j, k）沿行（水平向）、列（水平向）、层（垂向）方向上的6 个相邻网格单元，见图1；CRi,j-1/2,k、CRi,j+1/2,k、CCi-1/2,j,k、CCi+1/2,j,k、CVi,j,k-1/2、CVi,j,k+1/2 分别为网格单元（i, j, k）和周围6 个网格单元之间的水力传导度/（m2·d—1）；分别代表各网格单元第m 个模拟时段末第n 次迭代待求解的水头值/m；代表网格单元（i, j, k）在第m 个模拟时段初的水头值/m；Pi,j,k为作用于网格单元（i, j, k）之上的与水头有关的源汇项相关系数/ （m2·d—1）；Qi,j,k 代表作用于网格单元（i, j,k）之上的流量源汇项/（m3·d—1）；SSi,j,k 代表网格单元（i, j,k）的给水度或贮水系数；ΔRj 和ΔCi 分别代表网格单元（i, j, k）的行宽和列宽/m；tm、tm-1 分别代表第m 个模拟时段末和初的时刻/d。

2.1 单元间水平向水力传导度计算

全有效单元法采用不同于调和平均的新方法计算单元间水平向水力传导度，即CR、CC 项。以沿行方向上任意2 个相邻网格单元（i, j, k）、（i, j+1, k）为例，沿列方向类同：

式中：MaxB——沿行方向上2 个相邻网格单元相对较高的底板高程/m；

Boti,j,k、Boti,j+1,k——相邻网格单元的底板高程/m；

KRave——相邻网格单元之间的平均渗透系数/（m·d—1）；

HTave——相邻网格单元之间的平均有效饱和厚度/m；

ΔRj+1/2——相邻网格单元中心点之间的距离/m；

hup——相邻网格单元中水头相对较高的网格单元（即上风单元）的水头/m。

式（4）中单元间的平均渗透系数KRave 采用加权调和平均法计算：

式中：KRi,j,k、KRi,j+1,k——网格单元（i, j, k）、（i, j+1, k）沿行方向上的渗透系数/（m·d—1）；ΔRj+1——网格单元（i, j+1, k）的行宽/m。

下风单元为相邻网格单元中水头相对较低的网格单元。网格单元之间的平均有效饱和厚度HTave通过上风单元和下风单元的有效饱和厚度线性加权（图2）计算如下：

图2 全有效单元法下同层相邻网格单元之间的平均有效饱和厚度计算示意图
Fig.2 Schematic diagram of calculating the average effective saturated thickness between adjacent grid cells in the same layer using the AAC method

式中：HTup、HTdown——上风单元、下风单元的有效饱和厚度/m；

hdown——下风单元的水头/m；

ΔRup、ΔRdown——上风单元、下风单元的行宽/m；

Topup、Topdown——上风单元、下风单元的顶板高程/m；

Botup、Botdown——上风单元、下风单元的底板高程/m。

在全有效单元法中，疏干单元的水头只有计算意义而无水量意义。图3 展示了式（4）规定的水平向同层相邻网格单元在不同底板高程和水头状态下的8种地下水流动关系，根据式（4）和图3 可知，尽管疏干单元有计算意义的水头，但不存在任何地下水流出疏干单元，见图3（a）（b）（e），因此不会破坏水量平衡原则；同时地下水可以从同层相邻的湿润单元向疏干单元流入，见图3（c）（g），即地下水的侧向流动不会在疏干单元处中断，网格单元能够在疏干-湿润状态之间无缝变化。据此可将疏干单元保持为有效单元进行模拟，也就避免了丢失源汇项。

图3 全有效单元法下同层相邻网格单元之间的水分流动关系
Fig.3 Groundwater flow between adjacent grid cells in the same layer using the AAC method

2.2 垂向水力传导度计算

多层含水层模拟时涉及到网格单元间垂向水力传导度的计算。全有效单元法根据网格单元的垂向渗透系数和饱和厚度采用调和平均法动态计算越流系数和垂向水力传导度。用户还可以额外输入上、下2 层网格单元之间低渗透性介质的厚度和垂向渗透系数，从而模拟“拟三维”的情况：

式中：HSati,j,k——上层网格单元的含水层饱和厚度/m，为水头减去底板高程的值，当上层网格单元承压时，该值等于该网格单元的含水层厚度；

TKi,j,k+1——下层网格单元的含水层厚度/m，全有效单元法规定下层网格单元无论是否疏干，其含水层厚度都需要参与垂向水力传导度的计算，从而可避免由垂向水力传导度计算不连续导致模型收敛失败的情况；

TKCB——上、下层网格单元之间低渗透性介质的厚度/m；

VKi,j,k、VKi,j,k+1、VKB——上层网格单元、下层网格单元、低渗透性介质的垂向渗透系数/（m·d—1）。

2.3 贮水项计算

全有效单元法中有效网格单元可以分为湿润-承压、湿润-非承压、疏干3 种状态。在非稳定流模拟时，网格单元在3 种状态下将分别采用贮水系数、给水度、零值贮水系数计算贮水项，见图4，网格单元的水头在不同状态之间切换时还需要考虑贮水项的连续性处理。假定网格单元模拟时段初处于疏干状态，模拟时段末处于湿润承压状态，时段内网格单元贮水项可表达为：

图4 全有效单元法不同贮水系数的适用范围
Fig.4 Model cell for which storage factors are used during different states under the AAC method

式中：SYi,j,k、SCi,j,k——网格单元的给水度、贮水系数。

2.4 源汇项处理

全有效单元法中疏干单元是有效网格单元，可以接受源汇项。补给性质的源项作用于疏干单元不存在计算问题，然而排泄性质的汇项作用于疏干单元会违反水量平衡原理，因此对汇项有针对性处理。

（1）与水头有关的汇项，此时仅需要进行数据合理性的检查。以常年性河流为例，在河水位低于地下水头时地下水向河流排泄；若网格单元有底板高程限制条件，此时只要用户输入的河水位参数合理，即河水位始终不低于网格单元的底板高程，就不会产生计算问题；若用户输入的河水位参数不合理，全有效单元法将会报错。类似的汇项还有通用水头、潜水蒸发、排水沟、季节性河流、湖泊等。

（2）与地下水头无关的汇项，此时需要对汇项进行修正处理。这类汇项主要是开采井，全有效单元法根据网格单元饱和厚度阈值对开采量进行线性衰减处理：

式中：Φ——开采井网格单元的饱和厚度阈值/m，为用户输入值。一般 Φ值取网格单元厚度的0.10～0.25 倍比较合适，如果过小迭代计算不容易收敛；

QWell——用户在网格单元指定的开采量/（m3·d—1）；

h——网格单元当前的模拟水头/m；

QPump——模型在网格单元实际应用的开采量/（m3·d—1）。

除模拟时段初和模拟时段末网格单元分别处于疏干-湿润承压状态，还可能处于疏干-疏干、疏干-湿润非承压、湿润非承压-疏干、湿润非承压-湿润非承压、湿润非承压-湿润承压、湿润承压-疏干、湿润承压-湿润非承压、湿润承压-湿润承压状态，都需要对贮水项进行类似的连续性处理，详细过程可参考COMUS模型说明书[16]。

3 理想案例

理想案例用于验证全有效单元法模拟结果的可靠性。本案例由20 层、21 列网格单元组成，网格的长度、厚度均为5，2 m。模拟周期共20 d，等分为20 个模拟时段。在第20 层、第11 列网格单元处有一开采井，模拟周期前10 天开采量为30 m3/d，后10 天开采量为10 m3/d，用于模拟水位下降到恢复的过程。左、右侧2 列为40 m 定水头网格单元，定水头和井流边界共同作用使得网格单元反复疏干-湿润。网格单元的初始水位为40 m，给水度取0.1，第11 列网格单元的渗透系数取100 m/d，其他网格单元的渗透系数取1 m/d，分别使用试算法方案、全有效单元法方案和基于有限差分的饱和-非饱和二维水流模型VS2DT（variably saturated two-dimensional flow and transport model，由USGS 开发）对本案例进行模拟。VS2DT 模型无需面对网格单元的“干-湿转化”问题，可以作为参比标准比较试算法和全有效单元法模拟结果的可靠性。试算法方案参数为WETDRY = —0.001 m, NWETIT = 1,WETFCT = 1, IHDWET = 2，由于WETDRY 取正值难以收敛，因此取负值，同时|WETDRY|很小，因此是试算法中较为理想的参数方案。COMUS 模型的水头控制精度和流量控制精度分别取10—5 m 和10—4 m3/d，采用预处理共轭梯度求解器（preconditioned conjugate gradient solver, PCG）求解矩阵方程。

3 种方案第10 天末和第20 天末的模拟水位分布如图5 所示，可以看出不论在水位下降阶段还是在水位上升阶段，试算法方案、全有效单元法方案的模拟水位都低于VS2DT 模型的模拟水位，但全有效单元法方案的模拟结果与VS2DT 模型更为接近，同时在水位上升阶段3 种方案的模拟结果差别很小。水量平衡对比结果（表1）与水位对比结果一致，全有效单元法方案的模拟结果与VS2DT 模型非常接近，二者模拟结果都显示模拟域处于负平衡状态（蓄变量为负），而试算法方案模拟结果显示模拟域处于正平衡状态。使用试算法方案时，在模拟周期前10 天的水位下降阶段，一旦网格单元在某模拟时段内由时段初的湿润状态转化为时段末的疏干状态，该疏干单元立即被转变为无效网格单元，不再参与计算，该单元在模拟时段初处于湿润状态时的贮水量也不会被模型统计，因此产生了水位总体下降但是蓄变量为正的悖论。综上所述，全有效单元法模拟结果的可靠性优于试算法。

表1 理想案例水量平衡模拟结果对比
Table 1 Comparison of water balance in the ideal case /m3

图5 理想案例水位模拟结果对比
Fig.5 Comparisons of simulated groundwater heads in the ideal case

4 丹麦应用实例

本案例的目的在于检验全有效单元法的实用性，同时全面比较试算法和全有效单元法对于网格单元“干-湿转化”问题的模拟能力和算法特性。

4.1 流域概况和模型配置与模拟方案设置

汉普（Hampen）湖流域位于丹麦日德兰半岛中部斯凯恩河上游地区，临近日德兰半岛南北走向的主要地下水分水岭，见图6，总面积23.7 km2。流域地势东高西低，地下水大致自东向西流出本流域，见图7（a），流域内土地利用类型主要为森林和草地，见图7（b），年均降水量888 mm（1970—2007 年）。流域分为4 个水文地质单元，见图7（c），其中黏土层之上为非承压含水层，之下为承压含水层[18—20]。Hampen 湖总面积约为0.76 km2，湖泊水位与周围地下水位密切相关[18]，湖水可通过河道向下游排泄。

图6 汉普湖流域概况[18]
Fig.6 Information of the Lake Hampen basin in the practical simulation in Denmark [18]

图7 Hampen 湖流域模型配置
Fig.7 Model configuration of the Lake Hampen basin

基于试算法的模型配置已被证明在稳定流状态[19,21]和非稳定流状态[18]下都能有效模拟Hampen 湖流域的“河-湖-地下水”交互过程，本文简要介绍。Hampen 湖流域分为72 行、84 列、8 层共36 712 个有效网格，流域外围网格大小为100 m×100 m，在湖泊附近细化为50 m×50 m，见图7（a）（c）, 细致的含水层划分和水位升降使得大量网格单元反复疏干-湿润。流域西部概化为通用水头边界，北、东、南部为隔水边界，见图7（a）。根据植被覆盖类型和非饱和带厚度将流域分为4 个地下水补给区，见图7（b），补给率由HYDRUS-1D 模拟得到。每个网格平均意义上的水文地质参数由水文地质单元的水文地质参数求得，见图7（c）。湖泊采用Lu 等[22]提出的“倾斜湖底”算法模拟，河流采用Streamflow-Routing 模块模拟。模拟周期为2003 年12 月10 日至2012 年1 月7 日，共2 950 d。

模拟方案设置如表2 所示, 共包含7 个试算法模拟方案和一个全有效单元法模拟方案。Hampen 湖流域模型原始试算法配置如试算法5 方案所示[18—19,21]，本研究额外设置了6 个试算法模拟方案，共同说明试算法参数组合对模拟收敛性和模拟结果的影响。其中试算法1、2、3、4 方案的WETDRY 都相同，其它参数有所不同；试算法5、6、7 方案只有WETDRY 不同，其它参数都相同。COMUS 模型水头控制精度和流量控制精度分别取10—5 m 和10—4 m3/d，采用PCG 求解器求解矩阵方程。

表2 Hampen 湖流域模型模拟方案设置与模拟收敛性
Table 2 Simulation scheme settings and convergence of the Lake Hampen basin model

4.2 模拟收敛性的比较与分析

采用试算法1、2、3 方案模拟时，模型迭代不收敛，模拟失败；采用试算法4、5、6、7 方案时模型迭代收敛，模拟成功，见表2。试算法参数组合选取对“干-湿转化”模拟的收敛性有明显影响，这与前人的研究结论相符[4,12]。采用全有效单元法的模拟方案能够收敛，后续模拟结果的对比只在试算法4、5、6、7 方案和全有效单元法方案之间开展。

4.3 湿润地下水网格单元数量的比较与分析

各方案识别的湿润地下水网格单元数量对水位和水量平衡模拟结果有重要影响，各方案下逐日湿润单元数量变化见图8。全有效单元法方案识别的湿润单元数量最多，试算法下各方案识别的湿润单元数量随着|WETDRY|的升高而减少。除模拟初期外，|WETDRY|最小的试算法7 方案其余时段识别的湿润单元数量与全有效单元法方案基本一致。试算法下|WETDRY|越大，疏干单元越难重新湿润，而全有效单元法网格单元可以在疏干-湿润状态间无缝转换，可以认为|WETDRY|=0，因此全有效单元法方案和试算法7 方案识别的湿润单元数量较多是可以预料的。试算法4、5、6 方案识别的湿润单元数量在年内比较稳定，年际间呈现阶梯状分布，全有效单元法方案和试算法7 方案识别的湿润单元数量变化与作为当地地下水主要补给来源的面上补给量的动态之间存在良好的一致性关系，更好地体现了补给量变化引起的含水层连续“疏干-湿润”过程，表明全有效单元法方案和试算法7 方案的模拟结果更加合理。

图8 试算法方案与全有效单元法方案识别的湿润地下水网格单元数量对比
Fig.8 Comparisons of the numbers of the wet cells identified by the ET schemes and the AAC scheme

4.4 模拟水位与水量平衡的比较与分析

图9 为Hampen 湖和位于非承压含水层的地下水位观测井G1—G6（图7a）在2008 年全年和2011 年部分时段的各方案水位模拟误差对比，绝大多数水位模拟误差不超过0.5 m。|WETDRY|越小，试算法方案与全有效单元法方案的模拟结果越相似，主要是因为全有效单元法可以看作是|WETDRY|=0 的情况。同时也可以注意到，在任意一个观测站点，全有效单元法方案的模拟结果既存在误差大于试算法方案的情况，也存在误差小于试算法方案的情况（图9），说明采用全有效单元法模拟时，并未导致水位模拟精度出现系统性降低，全有效单元法的单元间水平向水力传导度算法实现了可以与调和平均法相比较的数值计算精度。

图9 试算法方案与全有效单元法方案的水位模拟误差对比
Fig.9 Comparisons of head errors simulated by the ET schemes and the AAC scheme

注：模拟水位误差指实测值与模拟值之差。

各方案模拟的流域逐日平均潜水位与湖泊水位如图10 所示。除模拟初期全有效单元法方案的模拟潜水位明显较高之外，其余时段2 种方案模拟的潜水位升降趋势基本一致，见图10（a）—（d）。随着|WETDRY|降低，两种方案的模拟结果愈发一致，|WETDRY|最小的试算法7 方案和全有效单元法方案模拟结果间的平均绝对误差（MAE）仅为0.013 m，见图10（d）。对比图8 和图10（d）可知，试算法7 方案与全有效单元法方案识别的湿润单元数量差异与模拟潜水位差异之间具有高度相似性，因此模拟期初两种方案较大的潜水位模拟偏差可能主要与识别的湿润单元数量差异有关。两种方案模拟的湖泊水位升降趋势基本相似，见图10（e）—（h）。特别地，试算法7 方案、全有效单元法方案模拟的湖泊水位几乎完全一致，R2=1，MAE仅为0.003 m，见图10（h）。总之，全有效单元法方案可以产生与理论上最为合理的试算法7 方案类似的水位模拟结果。

图10 试算法方案与全有效单元法方案逐日平均潜水位和逐日湖泊水位模拟结果对比
Fig.10 Comparison of simulated daily average water table and daily lake water level in the ET schemes and the AAC scheme

水量平衡方面，|WETDRY|相对较小的试算法4、5、7 方案和全有效单元法方案模拟的地下水系统水平衡状态总体一致，见表3，即研究区地下水的主要补给来源为面上补给，湖泊既向地下水系统渗漏，也接受来自地下水系统的排泄，存在双向联系，地下水系统处于正平衡状态，贮水量增加。同时随着|WETDRY|逐渐降低，2 种方案模拟的地下水系统蓄变量愈发接近，也即识别的水平衡状态愈发一致。|WETDRY|最大的试算法6 方案模拟结果显示地下水系统处于负平衡状态，主要原因是该方案模拟的湖泊渗漏补给量较少，可能是因为等同于隔水边界的大量无效疏干单元阻碍了湖泊向地下水系统的渗漏。

表3 试算法方案与全有效单元法方案模拟的水量平衡对比
Table 3 Comparison of simulated water balance by the ET schemes and the AAC scheme/104 m3

5 讨论

地下水模型中对于网格单元“干-湿转化”过程的模拟长久以来是困扰模型应用的一个难点[4,7,12]。目前被广泛应用的“干-湿转化”模拟算法是经验性的试算法，然而试算法参数组合对模拟的收敛性和模拟结果都有很大影响。例如在丹麦应用实例中试算法3、4 方案只有NWETIT 参数不同，但模拟的收敛性完全相反。再比如拥有最大|WETDRY|的试算法6 方案虽然能收敛，但模拟结果显示Hampen 湖流域地下水系统处于负平衡状态，其他拥有相对较小|WETDRY|的收敛方案都显示Hampen 湖流域地下水系统处于正平衡状态。使用试算法时需要反复调整参数，在保证模型收敛的前提下，寻找|WETDRY|最小的参数组合以保证模拟的相对准确性和合理性，很大程度上增加了用户使用模型的难度和时间成本。在过去数十年，包括人口增长和气候变化在内的多种因素导致全球干旱和半干旱地区地下水过度开采，进而使得含水层水位下降甚至枯竭[23]，具有鲁棒性的“干-湿转化”模拟算法对于科学合理地管理地下水系统的补给与排泄、恢复地下水位具有重要支撑作用[7,24—25]，目前基于试算法的模拟结果及衍生出的相关管理政策具有明显的不确定性。

全有效单元法的优点在于将疏干单元作为有效单元，从而能够从地下水动力学方程本身出发解决网格单元的“干-湿转化”模拟问题。理想案例模拟结果表明，全有效单元法模拟结果的可靠性优于试算法。通过丹麦应用实例研究发现，随着|WETDRY|的降低，试算法方案和全有效单元法方案的湿润单元数量、水位和水量平衡模拟结果愈发趋近，可以认为全有效单元法的作用等同于|WETDRY|趋近于零的最优试算法参数组合。使用全有效单元法时无需进行复杂的参数组合工作，很大程度上降低了模型的使用难度和模拟结果的不确定性。同时，全有效单元法方案和|WETDRY|最小的试算法7 方案模拟结果之间的高度相似性也说明全有效单元法中单元间水平向水力传导度算法实现了可以与经典调和平均法相比较的数值计算精度，这意味着全有效单元法在不涉及网格单元“干-湿转化”问题的地下水模拟中同样具有应用潜力。目前在MODFLOW-2005 等使用调和平均法计算单元间水平向水力传导度的模型中，是否启用试算法模拟网格单元的“干-湿转化”过程仍然被作为模型的可选项，全有效单元法的全场景应用潜力有助于弥合当前模型的不足，但这一潜力仍需要更多的理论案例和应用实例来验证，这也是下一步的研究方向。

最后需要说明的是，MODFLOW-USG[26]等基于非结构化网格的地下水模型同样面临着由“干-湿转化”模拟引发的各种异常情况，尽管全有效单元法目前集成于基于结构化网格的COMUS 模型之中，但该方法也可为基于非结构化网格的地下水模型应对“干-湿转化”问题提供有益借鉴。

6 结论

（1）试算法的参数组合选取对模拟的收敛性和模拟结果都有明显影响，使用试算法时需要不断优化参数组合，在实现模型收敛的前提下寻找|WETDRY|最小的参数组合以保证模拟的相对准确性和合理性，很大程度上增加了用户使用模型的难度和时间成本。

（2）相比于试算法，全有效单元法的模拟结果更为准确；全有效单元法的作用等同于理论上最优的试算法参数组合，使用全有效单元法时用户无需进行复杂的参数调优工作，因此该方法能有效降低模型的使用难度与模拟结果的不确定性。

（3）全有效单元法中单元间水平向水力传导度算法实现了可以与经典调和平均法相比较的数值计算精度，因此全有效单元法在不涉及网格单元“干-湿转化”问题的地下水模拟中同样具有应用潜力。