第 2 章：以块为中心的流程 (BCF4) 包

2.1 概述

MODFLOW 代码利用以块为中心的有限差分方法来求解地下水流方程。流域被离散成行、列和层，使得每个节点代表多孔材料的矩形块，称为单元。生成的有限差分网格中的节点表示无流量、可变水头或恒定水头单元，并且与单元相关的任何水力属性都是相对于相应的节点指定的。

对于无侧限流模拟，原始 MODFLOW 代码（McDonald 和 Harbaugh，1988）的 BCF1 包将去饱和可变头单元转换为无流单元。这会导致模型模拟中永久排除部分流域。如果电位水平恢复，MODFLOW 无法将这些单元恢复为可变头单元，这可能会导致代码产生误导或错误的模拟结果。

McDonald 等人 (1991) 尝试将无流动池重新转换为可变头池，以便在需要时允许池重新饱和。他们的再润湿方案在 BCF2 包（块中心流包，版本 2）中实施。该方案利用相邻单元处的头来确定是否将无流单元重新转换为可变头单元。不幸的是，BCF2 再润湿选项在再润湿期间或当撤回使各个单元变干时容易出现收敛和稳定性问题（McDonald 等人，1991）。这是因为用于重新激活干细胞的临时程序可能严重违反流量和质量守恒定律。此外，当相邻单元处的头部相关透射率减小到零时，代码中使用的谐波透射率平均方案在物理上变得不一致。 Goode 和 Appel (1992) 记录了与 BCF2 再润湿包相关的一系列问题以及一系列缓解方案，他们将 BCF3 包添加到 MODFLOW 中，以提供平均透射率的替代程序。鉴于上述再润湿程序遇到的严重困难，需要一种严格的方法来满足流动连续性要求并允许地下水位在无约束层中自由运动，而不需要任何形式的强制转换。 HydroGeoLogic, Inc. 开发了这样一种创新方法，通过使用 3-D 可变饱和流公式和自动生成的伪土壤保水函数，将非饱和流问题减少到寻求地下水位（即，电池中的压头为零或大气压）。该公式旨在准确描绘地下水位，并捕获无约束系统对抽水和补给的延迟产量响应。

模拟的数据要求没有改变，因为使用了伪土壤关系而不是真实土壤的本构关系。我们方法的主要优点包括：

• 对 3-D 流进行严格处理（即，没有 Dupuit 假设），不需要打开和关闭单元。

• 稳健的数值方案，具有良好的收敛性和稳定性特征以及处理细胞完全去饱和和再饱和的所需能力。

使用有限差分法，将带有伪土函数的可变饱和公式编程到 MODFLOW 现有的 BCF3 包中，新版本称为 BCF4 包。 BCF4 包的新选项将块间电导计算为块水力传导率、相对渗透率和平均流通面积的加权调和平均值的乘积。以下各节讨论了 BCF4 包中具有拟土壤函数的可变饱和公式的使用及其应用。

BCF4 软件包还提供使用轴对称圆柱坐标系执行轴对称分析的选项。此选项避免使用通过笛卡尔 (x,y,z) 坐标的完全 3-D 公式，从而大大节省计算量。轴对称分析选项适用于可接受轴对称几何形状的某些情况（即单井流动）。

MODFLOW 的模块化特性和输入数据准备的结构在 BCF4 包的开发中得以保留。新功能的数据要求与原始 MODFLOW 相同，不需要额外的参数值。为了保持向后兼容性并允许用户比较各种计算方案，BCF4 软件包还合并了以前的 BCF3 和 BCF2 软件包的所有选项。在本章中，通过对公式化、验证和演示示例问题的简要讨论，向用户系统地介绍了 BCF4 包的新功能。

2.2A 使用变饱和流动公式建立无侧限流动方程

水在可变饱和系统中的 3-D 运动可以表示为（Huyakorn 等人，1986）：

其中：

x、y 和 z 是笛卡尔坐标 (L)； 

Kxx、Kyy 和 Kzz 分别是沿 x、y 和 z 轴的水力传导率；

krw 是相对渗透率，是水饱和度的函数；

h为水头（L）； 

W 是每单位体积的体积通量，代表水源和/或汇； 

N 是可排水孔隙度，等于比产量；

Sw是水的饱和度，是压头的函数； 

Ss 为多孔材料的比储存量； 

t 为时间（T）。

对于完全饱和介质（即 Sw = 1.0），相对渗透率为 1，方程 (1) 简化为：

方程（2）是 MODFLOW 开发中使用的基本地下水流方程（McDonald 和 Harbaugh，1988）。因此，可变饱和流量方程简化为地下水位以下和受限系统中的传统地下水流量方程。拟土关系用于定义函数关系：Sw = Sw (R) 和 krw = krw (Sw)，其中 R 是压头，定义为 R = h+z，其中 z 在垂直向下方向为正值。请注意，方程 (2) 的以块为中心的有限差分近似包含块间水力电导率，这些电导率是作为参与单元的节点值的加权调和平均值获得的。

因此，对于 BCF4 包建模选项，用户指定水平饱和水力传导率（Kxx 和 Kyy）以及每个单元的顶部和底部高程。块间饱和水力传导率乘以相对渗透率即可得到块间有效水力传导率。垂直块间泄漏（垂直水力传导率 Kzz 除以节点之间的垂直距离 )z）由用户直接输入，就像在所有以前的 BCF 包中所做的那样。这些泄漏可以使用原始 MODFLOW 文档（McDonald 和 Harbaugh，1988）中讨论的各种方法来计算。使用改进的 Picard 方法处理变饱和流量方程中的非线性（Celia 等，1990）。与之前的 BCF 包类似，用户输入无量纲主存储系数 SF1。等式(1)使用的特定存储(Ss)是通过将存储系数SF1除以块厚度来内部计算的。

总之，BCF4 包中嵌入的具有拟土函数的可变饱和公式包括比以前的包有两个改进的方案： (1) 一种解决无侧限流动问题的新方法，它将去饱和和再饱和视为控制方程的一个组成部分但不需要土壤保水数据； (2) 附录 A 中所示的饱和水力电导率的调和平均比有效（水头相关）透射率的调和平均更合适。新选项不会使干电池失活，而是计算通过非饱和区域传输充电水所需的水头。请注意，新的 BCF4 软件包可用于全 3D、准 3D 和轴对称模拟。

对于轴对称流模拟，BCF4 包利用以下形式的饱和地下水流方程：

其中 r 和 z 是轴对称圆柱坐标系的径向和垂直坐标，Krr 和 Kzz 分别是沿径向和垂直轴的主要水力传导率。

当调用轴对称选项时，BCF4包将径向坐标（r）视为水平x坐标，而垂直坐标（z）不变，并且通过设置行数使水平y坐标无效网格块 (NROWS) 为 1。请注意，任何网格块 ()xi , )zk) 的轴对称配置是通过块绕 z 轴旋转生成的。

2.2B 变饱和水流方程的建立

水在可变饱和地下系统中的 3-D 运动由方程 (1) 表示。为了解决变饱和流动问题，还需要指定相对渗透率与水相饱和度、压头与水相饱和度的关系。使用两种替代函数表达式来描述相对渗透率与水饱和度的关系。这些函数由（Brooks 和 Corey，1966）给出：

和（van Genuchten，1977）：

当 n 和 ( 是经验参数，Se 是有效含水饱和度，定义为 Se=(Sw-Swr)/(1-Swr)，Swr 被称为残余含水饱和度。BrooksCorey 表达式产生与 van 类似的曲线Genuchten 函数，当选择参数 n 使得 n = 1 + 2/( (van Genuchten, 1980)。

按照 Mohanty 等人 (1997) 的建议，通过包含分段连续水力传导函数，进一步增强了模拟能力，以包含系统中优先流动路径的双峰描述。因此，如果系统中的饱和度小于临界值S*，则使用等式(4b)的van Genuchten函数。然而，当饱和度高于 S* 时，将应用指数电导率曲线，该曲线会在接近饱和时快速增加电导率，如下所示：

其中 S* 是临界或断点饱和值，其中流量从毛细管主导变为非毛细管主导，反之亦然，kr * 是对应于 S* 的相对渗透率，* 是代表有效大孔隙度或其他贡献结构特征的拟合参数到非毛细管主导的流动。还提供了表格输入选项，可以容纳其他功能形式。

压头（R = h-z，其中 z = 垂直向上坐标）与水饱和度的关系由以下函数描述（van Genuchten，1977，1980）：

其中α和β是经验参数，hc是毛细管水头，定义为(hap-Ψ)，其中hap，空气中的压头取大气压(=0)，Swr是残余水相饱和度。参数 β 和 γ 之间的关系为 γ = 1-1/β。可以在实验室中测量给定土壤的保水性和相对渗透性特性的 Brooks-Corey 和 van Genuchten 函数。

BCF4 包通过 BCF4 输入文件中的索引 IREALSL 使用实际土壤函数 [方程 (4) 或 (5) 和方程 (6)] 调用可变饱和流模拟。对于 van Genuchten 土壤函数，需要额外输入参数 "、$ 和 Swr。如果需要 Brooks Corey 相对渗透率函数，则还需要 Brooks Corey 参数 n 作为输入。其余参数与讨论的参数相同2.2 节中的 3-D、横截面或轴对称情况。

2.2C 变饱和空气流量方程的建立

可变饱和多孔介质中的空气流量可表示为（Huyakorn 等人，1994）：

其中 xi , xj (i, j = 1, 2, 3) 是笛卡尔坐标，koij 是固有渗透率张量，kra 是对空气的相对渗透率，μa 是空气的动力粘度，Pa 是空气相中的压力，ρa 是空气密度，Wa 是代表空气源和/或汇的每单位体积多孔介质的空气体积通量，Sa (=1-Sw) 是空气相饱和度。

在MODFLOWSURFACT/MODHMS框架下求解式(7)时，可以方便地将方程转化为MODFLOW类型变量。因此，空气相需要压头和水头的等效定义：

其中hap是空气相的压力水头，Ma是空气相的势能，Dw是水的密度，g是重力加速度，haM是空气相的势水头或总水头。进一步地，根据定义，水力传导率Kij可表示为： (7)

其中，hap是空气相的压力水头，Φa是空气相的势能，ρw是水的密度，g是重力加速度，haΦ是空气相的势水头或总水头。进一步地，根据定义，水力传导率Kij可表示为：

因此，等式（7）左边的第一项可以使用等式（8）和（9）来操作并重写为

方程 (7) 中的最后一项可以以与 Bear (1979, p. 91) 对水相所做的类似的方式对空气相进行操作，得到：

其中 Ssa 是空气的特定体积存储量，定义类似：

其中 αs 是土壤基质的压缩性，βa 是空气的压缩性。对于大多数含水层系统来说，Ssa 不是一个有形参数，需要进一步操作才能根据 Ss（多孔材料的特定存储量）提出方程，定义为（Bear，1997）

其中 βw 是水的压缩性。等式（12）和（13）产量。

使用(7)中的方程(10)、(11)和(14)给出了以MODFLOW类型变量和基本参数表示的空气相流方程：

仅求解空气相流方程时，必须对系统中的水相做出假设，因为水被视为被动相。

MODFLOW-SURFACT/MODHMS 中提供了两个选项，用于在求解气相流时处理水相。作为第一个选项，假定水处于静水平衡状态，具有恒定的水头，等于域中任何 x、y 位置处地下水位的初始高程 ZWT (x,y)。

因此，域内任何高度 z 处的水相压头 hwp（=之前使用的 Ψ）如下：

其中水压 hwp 等于地下水位处的环境气压 (=0)。

因此，方程（6）中所需的毛细管水头[hc=hap-hwp]可通过使用（19）和（8）获得：

作为第二个选择，可以从同一域上的水流方程的稳态解获得方程(6)的环境水压。因此，通过在环境空气条件下提供正确的空气和水饱和度，可以在空气相流计算期间适应再充电/放电系统。

系统压力条件下的空气密度由理想气体定律得出：

式中ρstd为大气压条件下的标准空气密度，Patm；Pa为系统内空气的绝对压力；并且假设系统条件下空气的表观分子量与标准条件下的相同。重新整理方程（8）并表达绝对压力条件下的气压可得出。

其中参考压力 Patm 出现在参考高程 z=0（即基准面）处。空气密度对系统势能的影响通常很小，并且无论基准面高程如何，都可以针对平均海平面条件规定 Patm，而不会显着降低精度。

将 (22) 代入 (21) 并重新排列，以 MODFLOW-SURFACT/MODHMS 变量的形式给出空气相密度：

总结可变饱和空气相流的公式，MODFLOWSURFACT/MODHMS 使用方程 (15) 来控制受方程 (6) 的毛细管饱和关系和方程 (16) 或方程 (18) 的相对渗透率行为影响的气流行为。毛细管水头由方程（20）确定，采用水的瞬时静水平衡的被动水相假设，或者以环境 hwp 作为输入（从理查兹方程的稳态解获得）作为 hi = hap-hwp 确定。

最后，根据方程（23）计算空气相密度。

BCF4 包使用 BCF4 输入文件中的索引 IREALSL 调用空气相流公式。除了土壤湿度关系和相对渗透率函数所需的参数外，空气相流模拟还需要输入 Dw、μw、μa、$w、$a、Catm、g 和 ZWT (x,y)。其余的仿真输入参数与 2.2 节中针对 3-D、横截面或轴对称情况讨论的参数相同。

BCF4 包使用 BCF4 输入文件中的索引 IREALSL 调用空气相流公式。除了土壤湿度关系和相对渗透率函数所需的参数外，空气相流模拟还需要输入 ρw、μw、μa、βw、βa、Patm、g 和 ZWT (x,y)。其余的仿真输入参数与 2.2 节中针对 3-D、横截面或轴对称情况讨论的参数相同。

2.3A 无侧限流验证示例

BCF4 包通过将数值结果与 Neuman (1972) 的解析解进行了验证，该解析解针对无承压含水层中的井流，该含水层对抽水表现出显着延迟的产量响应。

所考虑的问题涉及以 3658 ft3/d 的速率从均质且各向同性的无承压含水层中抽水的完全穿透井。含水层的初始饱和厚度为 40 英尺，井在含水层厚度较低的 30 英尺部分进行筛选。含水层性质如下：

Hydraulic conductivity, K = 21 ft/d 

Specific yield, Sy = 0.132 

Specific storage, Ss = 1.2 x 10–4 ft–1

建模方法和结果

上述流动问题可使用 BCF4 软件包中提供的全 3D 和轴对称仿真选项来解决。将获得的数值结果与解析解进行比较。

对于完全 3-D 模拟运行，模型区域是一个尺寸为 2,000 英尺 x 2,000 英尺的正方形，井位于中心。为了利用对称性，我们 (22) (23) 仅离散流域的一个象限并考虑四分之一的井流量 Q。3-D 网格由 20 层组成，每层 2 英尺厚，有 21 行， 21 列。水平网格间距 ()x 和 )y) 从 0.5 英尺到 100 英尺不等。VCONT（垂直导水率除以两个相邻节点层之间的层间距离）的值计算为 10.5。每层的恒定主存储系数（SF1，定义为特定存储时间块厚度）被指定为 2.4 x 10–4。井流量（象限网格为 Q/4）均匀分布到第 6 层到第 20 层第 1 行、第 1 列的 15 个井节点。假设外侧、底部和顶部边界为无流动边界条件。强隐式过程 (SIP) 包用于求解矩阵方程。

对于轴对称仿真运行，模型域是半径为 2,000 英尺、厚度为 40 英尺的圆柱体。该网格由 20 层、2 英尺厚、21 行和 1 列组成。水平（径向）网格间距（△x≡△r ）从0.5英尺到100英尺不等。全井流量Q，均匀分布到第6层至第20层第1行第1列的15个井节点。其余信息与完整 3D 模拟运行给出的信息相同。

对距离井 31 英尺、位于底部、中间和顶部节点层中心的三个选定观测点进行了水位下降与时间关系的比较。图 2.1 中绘制的无量纲回撤 sD 和无量纲时间 ty 定义为：

其中b是初始饱和厚度，其余符号如前定义。

图 2.1 将全 3-D 和轴对称模拟的数值结果与 Neuman 解析解进行了比较。对于底部和中间的观测点，数值解和解析解非常一致。对于顶部观测点，两个数值解非常吻合，但偏离了早期值的解析解。这种偏差可能是由于恒定饱和厚度的简化假设和解析解所使用的近似地下水位边界条件的影响。请注意，我们的数值解基于严格的可变饱和流公式。

2.3B 水不饱和流量验证示例

提供了两个示例问题来演示 BCF4 包中非饱和区水流选项的验证和应用。第一个问题考虑的是非常干燥的土柱中水的瞬时渗透。第二个问题考虑非饱和矩形土板中的二维流动。将这些问题与适用于各自情况的数值解进行比较。

2.3B.1 问题 1 - 非饱和垂直土柱中的瞬时渗透

Huang 等人 (1996) 详细介绍了这个问题，涉及非饱和垂直土柱中的一维渗透，该土柱长 160 厘米，离散为 80 层，厚 2 厘米。整个域的初始水头为-104 cm，代表极其干燥的条件，并且在 t=0 时在顶部应用 50 cm/d 的恒定渗透通量。

对典型的粗质地土壤进行了长达 0.8 天的分析，其水力特性如下：

van Genuchten 函数用于保湿性和相对渗透性。

建模方法和结果

选择 PCG4 求解器来求解收敛公差为 0.05 cm 的矩阵方程组。使用AT04包进行自适应时间步进，初始△t=0.2×10-7 d，最大△t=0.01 d，时间步乘数为1.3，时间步除因数为2.7。图 2.B1 显示了 0.4 d 和 0.8 d 时的预测水头剖面。这些结果与使用 Huang 等人 (1996) 和 HydroGeoLogic (1992) 的有限差分和有限元模型获得的结果非常吻合。这个问题还证明了代码处理高度非线性情况和极其干燥条件的能力。

2.3B.2 问题 2 - 非饱和矩形土板中的二维流动

该问题涉及截面尺寸为水平 15 厘米和垂直 10 厘米的矩形土板中的二维非饱和流动。使用 )x=1 cm 和 )z=1 cm 的恒定节点间距来离散化域。整个域的初始水头为-90 厘米。沿着域的左边缘的顶部4cm被规定为6cm的边界水头，并且域的右边缘被维持在-90cm的规定水头。其余边界（即顶部边界、底部边界和左边缘下 6 厘米）为无流动条件。土壤的水力特性如下：

假设土壤是均匀且各向同性的。 van Genuchten 关系用于保湿曲线，而 Brooks Corey 关系用于定义相对渗透率函数。

建模方法和结果

选择 PCG4 求解器来求解收敛公差为 0.01 cm 的矩阵方程组。初始时间步长 0.01 天扩大了 1.2 倍，最大时间步长为 0.05 天。图 2.B2 绘制了时间值为 0.508 天时沿流动区域底部和顶部的水头计算值。这些结果与使用 HydroGeoLogic (1992) 有限元模型获得的结果非常吻合。

2.3C 气相流模拟验证示例

通过将数值结果与使用综合多相流模拟器获得的结果进行比较，验证了 BCF4 包的空气相流模拟选项（Huyakorn 等人，1994 年）。

Panday 等人 (1994) 详细介绍了这里考虑的问题。选择这个问题来研究空气喷射过程的流动动力学。本研究使用了 Marley 和 Magee (1992) 提供的相关现场信息。相关地点已被包气区和饱和区几个泄漏的地下储罐的汽油污染。地下水位平均深度为地表以下-2.4 m。污染面积-450平方米。通过仅考虑可变饱和气流（即假设被动液相），可以对注入井附近的喷射影响半径和空气速度进行初步评估。考虑 30 m × 20 m × 3 m 的矩形域，并将其离散为 900 个节点，其均匀节点间距为 )x = 2m、)y = 2m 和 )z = 0.5m。

初始条件对应于包气区的主要大气压力 (haM=0)。通过在位于 x = 9 m、y = 9 m 和 z = 0.75 m 的节点处以 10-3m3 s -1 的稳定速率注入空气来创建喷射过程。地下水位位于喷水井下方 0.4 m 处（即 z=0.35 处）。在整个模拟过程中，表面边界保持大气压条件，底部和侧面边界被视为无流动条件。模拟中使用的土壤和流体属性采用 Marely 和 Magee (1992) 给出的值，如下所示。

van Genuchten 函数用于描述空气-水系统的相对渗透率和毛细管压力。

建模方法和结果

针对这个问题进行了瞬态和稳态仿真。

瞬态仿真调用AT04包进行自适应时间步进，初始时间步长为1秒，最小△t=1×10-3s，最大△t=1×105s，时间步乘数为1.4，时间步长缩减因子为 2。选择 PCG4 软件包来求解矩阵方程组，瞬态情况下最多 20 次非线性迭代，稳态模拟最多 100 次非线性方程。瞬态和稳态模拟的头部闭合标准分别为 0.1 和 0.21 m。

图 2.C1 描述了沿穿过井的水平线的模拟气压分布，其中包括两种模拟以及在 t=7 分钟和稳态条件下使用 Huyakorn 等人 (1994) 的有限元代码进行的模拟。预测和观察到的压力分布与这两个规范都相当一致。有趣的是，瞬态模拟仅在 1.16 天内就达到了稳态条件，显示了空气相中发生的快速平衡过程。这些模拟的绝对气压是通过等式（22）从电位头计算的，Da是根据等式（23）计算的。

2.4 应用示例

提供了两个示例问题来演示新 BCF4 包的可变饱和流选项的应用。在这两个示例中，都允许发生足够的抽水以使无承压含水层的部分去饱和。第一个示例是稳态泵送模拟。第二个例子是瞬态模拟，其中在四年的压力期后关闭油井以恢复含水层。该场景也使用 MODFLOW 中现有的润湿选项（McDonald 等，1991）和 STAFF3D 有限元代码（HydroGeoLogic, Inc., 1992）进行建模。下面对不同方法获得的结果进行比较和讨论。

2.4.1 问题 1——无承压含水层抽水的稳态模拟

该测试问题考虑了图 2.2 所示的 300 英尺厚的无承压含水层。

建模域是一个尺寸为 75,000 英尺 x 75,000 英尺的正方形。含水层的顶部和底部海拔分别为 50 英尺和 –250 英尺。该区域受到 3.0 x 10-9 英尺/秒的均匀、连续的垂直补给。图中左（西）边界为恒水头边界，其余边界为无流边界。含水层底部 100 英尺处有 15 口井，每口井的抽水速度为 0.95 英尺3/秒。井位如图2.2所示。含水层参数包括：

建模方法
该域被均匀离散为 3 层、15 行、15 列的网格块（即 △z = 100 ft，△x = △y = 5,000 ft）。因此，VCONT（垂直导水率除以两个相邻节点层之间的层间距离）的值为10-7s-1。顶层和中层西侧规定恒定水头为零。代码所需的初始头部分布在整个域中统一设置为零。 15 口井位于第 3 层。虽然这些井充当汇，但地表补给和西部的恒定水头边界为含水层系统提供补给，从而维持所需的稳定流量条件。

执行两次单独的稳态仿真运行：(a) 所有三层的输入 LAYCON 值为 3，(b) 输入 LAYCON 值为 43。请注意，LAYCON 是一个控制参数，指示特定节点层所需的含水层处理选项。案例 (a) 使用原始 MODFLOW 公式来计算水平块间透射率和电导，而案例 (b) 则激活 BCF4 软件包中包含的新选项。 （参见BCF4包的输入指令。）对于这两种情况，SIP包用于求解矩阵方程。除了各层的输入 LAYCON 值不同之外，BCF4 包和其他包（Basic、Well、Recharge、SIP 和 Output Packages）的剩余输入数据是相同的。用于模拟的输入数据文件列于附录 B 中。

为了保持对流动系统的连续充电，规定的充电速率应用于每个垂直列中的最高活性单元。如果任何接受再充电的活性电池变干，这对于情况（a）来说是必要的，以允许连续渗透。原因是，如果可变水头电池干燥并且不将补给传输到地下水位，则可变水头电池会转换为无流动电池。如果使用新的可变饱和流选项[情况(b)]，则情况并非如此，并且避免了数据输入中的这种复杂性。

结果与讨论

图 2.3 至 2.5 分别对第 1、2 和 3 层由先前 MODFLOW 选项 [案例 (a)] 和新的可变饱和流选项 [案例 (b)] 产生的稳态水头分布进行了比较。尽管域中的所有单元都作为可变头输入，但当使用 MODFLOW [案例 (a)] 中的先前 BCF1 选项进行模拟时，第 1 层中大约四分之一的单元被停用（见图 2.3），因为模拟水头值低于该层的底部高程（对于第 1 层为 –50 英尺）。请注意，可变饱和流选项[情况 (b)] 预测去饱和区域内的水头分布，从而允许将垂直补给传输到地下水位以及来自相邻单元的侧向流分量。

第 2 层和第 3 层的稳态水头分布（分别如图 2.4 和 2.5）表明，随着距去饱和区距离的增加，两种模拟选项之间的一致性会提高。在域的东南角附近，两次模拟的头部分布存在显着差异。对案例 (a) 和案例 (b) 的结果进行仔细检查表明，BCF4 封装的可变饱和流选项中使用的水平块间电导表现适当。原始方案[情况(a)]中使用的与水头相关的透射率的谐波加权会导致更大的下降，因为当下降较大时，透射率会向较低的值加权。从物理角度来看，谐波加权应应用于饱和导水率，而不是与水头相关的有效透射率。 （演示见附录 A。）

2.4.2 问题 2——无承压含水层抽水和采收的瞬态模拟

本例中的瞬态分析还考虑了 300 英尺厚的无承压含水层，如图 2.2 所示。与问题 1 一样，建模域对应于尺寸为 75,000 英尺 x 75,000 英尺的正方形。含水层的顶部和底部分别位于 50 和 –250 英尺的海拔处。整个域中含水层的初始地下水位为 0 英尺。该域受到 3.0 x 10 –9 英尺/秒的均匀、连续补给。图中的左（西）边界是恒定（零）水头边界，其余边界是无流边界。如图 2.2 所示，底层厚度超过 100 英尺，有 15 口井被筛选。

然而，在此问题中，每口井的抽水速率高达 3.85ft3/s 。四年后，井中的抽水将被关闭。 （见图 2.6。）然后允许含水层恢复 12 年。含水层参数为：

建模方法

与问题 1 一样，域被均匀离散为 3 层、15 行和 15 列。 MODFLOW所需的VCONT（垂直导水率除以相邻节点层之间的层间距离）的值计算为10–7s–1。

每层指定恒定的主存储系数（SF1 定义为特定存储时间块厚度）10 –4。顶部和中间离散层的西侧规定恒定水头为零。最初，假设域中的头为零。位于第 3 层的 15 口井在前 4 年允许以 3.85 的速率抽水。在随后的 12 年期间停止抽水，以便恢复含水层。

执行三个独立的瞬态仿真运行。前两个模拟分别利用 (1) 现有的 MODFLOW 再润湿功能和 (2) MODFLOW-SURFACT/MODHMS 的新可变饱和流（BCF4 封装）选项。这两个模拟分别称为情况（1）和（2）。第三次模拟运行 [案例 (3)] 使用有限元 STAFF3D 代码独立执行（HydroGeoLogic, Inc., 1992）。对于基于 MODFLOW 的仿真运行，SIP 包用于求解矩阵方程。除了 BCF4 包输入之外，两种情况使用的其余输入数据是相同的。请注意，情况 (1) 需要额外的数据才能重新润湿干电池。

为了在 MODFLOW 再润湿包 [案例 (1)] 模拟选项中保持连续充电，表面充电应用于每个垂直列中最高的活性电池（请参阅充电包输入）。这一规定是必要的，以允许接收补给的活性单元的连续渗透（如果有的话）变得干燥，因为可变头单元如果干燥则转变为非流动单元，并且不将补给传送到地下水位。

结果与讨论

对于使用原始再润湿能力 [案例 (1)] 和新的可变饱和流选项 [案例 (2)] 执行的模拟运行，第 1 层和第 3 层（代表两种极端饱和条件）的水头分布如图 2.7 所示到 2.10。图 2.7 比较了抽水周期结束时（t=4 年）第一层的水头分布。对于情况 (1)，第 1 层中大约四分之一的细胞因 4 年的连续泵送而失活，但第 3 层仍然保持完全活跃。

请注意，较长的泵浦周期会导致第 3 层去饱和并引发数值解的不稳定。图 2.8a 和 2.8b 比较了时间值为 4 年时第 3 层的头部分布。第 3 层中没有失活的单元。案例 (1) 和 (2) 的水头值之间存在显着差异（高达 20 至 30 英尺），发生在各个井附近。远离井场的两个水头分布非常吻合。图 2.9 和 2.10 分别显示了恢复期结束时（t=16 年）第 1 层和第 3 层的水头比较。一般来说，MODFLOW 的先前（再润湿）选项会低估水头值。两个水头分布之间的差异在流域的东南部更为明显。

对三个选定的观察点进行了水头与时间关系的比较。下面参考图 2.2 给出了这些点的位置。

三个观察点处水头的时间变化（图 2.11）表明，两种情况 [情况 (1) 和 (2)] 非常一致，直到情况 (1) 的方案将可变水头单元转换为无流动单元。经过这些转换后，案例（1）的结果与案例（2）的结果有偏差，但遵循相似的趋势。如前所述，BCF2 包 [案例 (1)] 的临时再润湿程序可能无法产生可靠的数值解。例如，如果在瞬态模拟的第一阶段，井的抽水速率从 3.85 增加到 3.90ft3/s，MODFLOW 的再润湿选项将无法收敛。然而，新的可变饱和流选项会产生收敛解。 BCF2 包中使用的再润湿程序失败是由于迭代之间电池的循环干燥和润湿（翻转）产生的不稳定性。

使用 MODFLOW-SURFACT/MODHMS 可变饱和流选项 [案例 (2)] 获得的结果也与使用 STAFF3D 代码获得的结果（称为案例 (3)）进行了比较。图 2.12 显示了三个观察点的头部时间变化图。在整个模拟期间，情况（2）和（3）之间有很好的一致性。

2.5 以块为中心的流包的其他扩展

对以块为中心的流 (BCF4) 包进行了一些细微的增强，以允许模拟具有更大的灵活性。

第一个修改允许输入网格块的垂直水力传导率，作为输入相邻层网格块之间的泄漏的选项。 BCF4 封装输入的第一行中包含了一个标志 IVHYC，以表明这一点。

如果 IVHYC 设置为 1，则代码需要输入所有层的垂直水力电导率，并且根据每个网格块的垂直电导率及其间隔距离的体积加权调和平均值在内部计算泄漏。如果 IVHYC 设置为零，则代码需要输入层之间的泄漏，如原始 MODFLOW 中一样。

第二个修改允许逐个单元地输入水平各向异性，作为每个模型层具有统一的各向异性值的选项。 BCF4 包输入的第一行中包含一个标志 IANIXY，以表明这一点。如果 IANIXY 设置为 1，则代码需要逐个单元输入各向异性。如果 IANIXY 设置为零，则代码期望像原始 MODFLOW 中那样逐层输入各向异性。

第三种修改允许输入曲线网格，作为使用矩形网格在区域平面中离散化域的选项。 BCF4 包输入的第一行中包含一个标志 ICURVL，以表明这一点。如果 ICURVL 设置为 1，则使用曲线网格，并且代码期望为层内的所有节点输入 DELR 和 DELC。由于层彼此垂直堆叠，因此所有层的 DELR 和 DELC 值都是根据第一层的值在内部计算的。如果 ICURVL 设置为零，则像原始 MODFLOW 中一样使用矩形网格。对于此 MODHMS 情况，仅输入 DELR 的列数 (NCOL)，仅输入 DELC 的网格中的行数 (NROW)。 MODHMS 为非饱和区水流模拟提供了额外的灵活性，允许通过重力分离垂直平衡 (GSVE) 伪土壤函数提供的受限/非受限地下水流模拟功能来表示更深的层。此功能对于非饱和区域流量模拟非常有用，可以避免地下水位附近的过度垂直离散，这需要准确表示无承压水位条件和相关的含水层流量。因此，当 IREALSL=1 或 2（不饱和/饱和水流模拟）时，可以设置顶部几层来求解 LAYCON 值为 43 的理查德方程，而下面较厚的层可用于表示无承压/承压地下水系统。这是通过在 MODHMS 中为这些较低层设置 LAYCON=44 来实现的，这些层将被视为 IREALSL=0。

2.6 BCF4 封装的输入指令

2.6.1 概述

BCF4 包所需的输入数据与之前的 BCF 包类似。 BCF4 包采用了新的可变饱和流公式，通过使用伪土壤函数，通过去饱和和再饱和来严格处理无侧限流问题。出于比较目的，我们还将所有先前记录的方案包含在新包中。

因此，我们的 BCF4 封装包括永久停用干电池的 BCF1 无限制封装（McDonald 和 Harbaugh，1988）、将无流动电池转换为可变头电池的 BCF2 封装（McDonald 等人，1991）以及 BCF3 封装，它使用块间透射率关系的不同选项来模拟地下水流（Goode 和 Appel，1992）。建议有兴趣比较新旧方案的用户熟悉各自的 BCF 包文档。

BCF4 的输入数据主要由属性数组组成，用于描述水平水力传导率/透射率、垂直泄漏、单元的顶部和底部标高以及瞬态模拟的存储系数和/或比产量。这些属性数组取决于每层的输入 LAYCON 值，如下所述。 BCF2 和 BCF3 无限制存储方案所需的附加参数数组是 WETDRY，当该选项的润湿标志处于活动状态时使用该数组。 BCF4 内存分配包在读取某些标志后为这些数组保留空间。

然后根据用户输入初始化这些数组，如下所述。以下页面提供的输入说明适用于 BCF4 包。请注意，BCF4 的设计旨在保持与其前身的整体兼容性。

要激活新的可变饱和流选项，用户将层的输入 LAYCON 指定为 40（如果层受到严格限制）或 43（否则）。严格限制层定义为在整个模拟期间测压头保持在该层顶部标高之上的层。如果没有先验知识可用于层中的水头条件，建议用户使用输入 LAYCON 43。BCF4 中实施的新无侧限流公式消除了输入准备中的主要尴尬以及先前单元失活中固有的困难和再润湿计划。我们建议用户改用这种新方法以获得更可靠的收敛模拟。

BCF4 包的 LAYCON=40 或 43 方案使用相邻网格块的水力传导率的加权谐波平均来计算单元之间的电导。通过设置 LAYCON=50 或 53，可以选择使用相邻网格块的水力传导率的加权算术平均来计算单元之间的电导率（类似于有限元装配的执行方式），作为额外的灵活性。

2.6.2 BCF4 输入

BCF4 包的输入是从 IUNIT (1) 中指定的单元读取的。下一小节将提供输入变量的简要说明。数据在模拟开始时读取一次。由基本包读取的 IUNIT 数组确定 MODFLOW 模拟中使用的选项或包。 （有关详细信息，请参阅原始 MODFLOW 文档的基本包。）请注意，BCF4 所需的额外输入变量显示在阴影框中。

以下二维数组的子集用于描述模型网格每一层的属性。每层所需的阵列取决于每层的输入 LAYCON 值（参见第 2 项）、模拟是瞬态（ISS = 0）还是稳态（ISS 不为 0），以及润湿功能是否处于活动状态（IWDFLG非 0) 对于 BCF2 和 BCF3 的润湿方案。请注意，对于 BCF4 封装中引入的新的可变饱和流量方案，IWDFLG 应设置为零（即输入 LAYCON=43）。

如果不需要数组，则必须将其省略。不存在需要所有数组的情况。首先读取第 1 层所需的数组（第 8-15 项）；然后是第 2 层的数组

限制为 200 层。如果层数在40层以下，则使用1条记录；否则，使用两个或多个记录。将记录中未使用的元素留空。 LAYCON 是可以实现 BCF4 封装新选项的参数，如下所述。

BFC4 封装支持 20 个不同的输入 LAYCON 值。这些值以 I2 格式输入。根据该值，第一个数字（十位）（如果有）根据该数字的值确定块间透射率函数类型的应用。输入 LAYCON 值的第二个数字（个）决定了原始模型中的层特征。第二个数字的含义如下（McDonald and Harbaugh，1988）：

上述描述可以总结如下（Goode and Appel，1992）：

请注意，BCF4 包的新可变饱和流选项是通过模型网格中所有层的输入 LAYCON = 40、43、50、53 或 54 来实现的。此外，LAYCON=44 或 54 的值可用于在其他不饱和区域水流模拟（即 IREALSL …0）中将层表示为无约束/约束层。因此，当输入 LAYCON 值的第二位为 4 时，其含义与上面讨论的 3 相同。此新选项需要 IAPART 的非零值。标志 IAPART 从 BAS 包的第 5 项中读取，指示数组 BUFF 是否与数组 RHS 分开。

由输入 LAYCON 的第一个数字表示的因子 10 存储在 LAYAVG 数组中。第二个数字（个）存储在 LAYCON 数组中。下表给出了输入 LAYCON 值的 LAYCON 和 LAYAVG 的存储值：

因此，输入 LAYCON 值的十位部分存储在数组 LAYAVG 中，输入 LAYCON 值的个位部分存储在 LAYCON 数组中。 LAYAVG 的存储值决定了计算块间透射率的方法（当 LAYAVG = 10、20、30、40 或 50 时）。

TRPY – 是一个一维数组，包含每层的各向异性因子。它是沿列的透射率或水力传导率（以所使用的为准）与沿行的透射率或水力传导率的比率。对于各向同性条件，设置为 1.0。这是一个单数组，每层有一个值。不要为每一层读取一个数组；整个数组仅包含一个数组控制记录。如果所有层的值都相同，则可以通过将 LOCAT 设置为 0 并将 CNSTNT 设置为适用于所有层的值来指定整个数组。 LOCAT 和 CNSTNT 的解释见下文。

注意：轴对称仿真不需要 TRPY 数组 (IAXIS>0)。

TRPYNU – 是一个二维数组，包含每层内逐个单元的各向异性因子。它是沿列的透射率或水力传导率（以所使用的为准）与沿行的透射率或水力传导率的比率。对于各向同性条件，设置为 1.0。为每一层读取一个数组。

注意：轴对称仿真不需要 TRPYNU 数组 (IAXIS>0)。

LOCAT – 是填充相关数组的值位置的索引。

<0，读取包含数组值的未格式化记录 = 0，将所有数组值设置为等于常量 (CNSTNT) >0，读取包含数组值的格式化记录。

CNSTNT – 是一个常量，如果 LOCAT 等于 0，则所有数组值都将设置为该常量；如果 LOCAT 不等于 0，则所有数组值都会乘以该常量。

DELR – 是沿行的单元格宽度。读取每个 NCOL 列的一个值。这是一个数组，每一列都有一个值。对于轴对称模拟 (IAXIS>0)，DELR 表示圆柱形电池的径向宽度 ()r)。

DELRCV – 是曲线网格沿行的单元宽度。这是一个二维数组，对于一层的每个 NROW 行，每个 NCOL 列包含一个值。

DELC – 是沿列的单元格宽度。为每一 NROW 行读取一个值。这是一个数组，每一行都有一个值。

注意：对于轴对称模拟 (IAXIS>0)，不会读取 DELC 数组。

DELCCV – 是曲线网格沿列的单元宽度。这是一个二维数组，其中包含一个层的每个 NROW 行的每个 NCOL 列的值。

ZWT – 是用作气流计算参考的域中的地下水位高程。

PWC – 水相的水头，用作气流计算的参考。

SF1——初级存储系数，可存储性。对于瞬态模拟只读（稳态标志 ISS 为 0）。对于输入 LAYCON 等于 1、11、21 或 31，SF1 将始终为具体良率，而对于输入 LAYCON 等于 2、12、22、3、13、23、33、40 或 43，SF1 将始终为受限存储系数。对于输入 LAYCON 等于 0、10 或 20，SF1 通常是受限存储系数。主存储系数除以块厚度等于等式(1)和(2)的特定存储(Ss)。

注 (1)：输入 LAYCON 值 0、10 或 20 也可用于模拟地下水位条件，其中水位下降预计在任何地方都保持为饱和厚度的一小部分，并且上面没有层或流量从上面看可以忽略不计。在这种情况下，将为 SF1 输入特定的屈服值。

TRAN – 是沿网格行的透射率。 TRAN 乘以 TRPY 即可获得沿列的透射率。对于由以下输入 LAYCON 值之一表示的层，TRAN 是只读的：0、2、10、12、20、22。

HY – 是沿网格行的水力传导率。 HY 乘以 TRPY 即可获得沿柱的水力传导率。对于由以下输入 LAYCON 值之一表示的层，HY 是只读的：1、3、11、13、21、23、31、33、40、43。

VHY – 是每个单元的垂直水力传导率。对于由以下输入 LAYCON 值之一表示的层，VHY 是只读的：1、3、11、13、21、23、31、33、40、43。

BOT——含水层底部的高程。对于由以下输入 LAYCON 值之一表示的层只读：1、3、11、13、21、23、31、33、40、43。

VCONT——垂直导水率除以从一层到下一层的厚度。单元的值是水力传导率除以该单元中的节点与下面单元中的节点之间的材料的厚度。因为底层下面没有层，所以没有为底层指定VCONT。

SF2——是二次储存系数，即比产量。如果模拟是瞬态的（稳态标志 ISS 为 0）并且该层的输入 LAYCON 值为以下之一，则 SF2 为只读：2、3、12、13、22、23、33、43。

请注意，SF2 是求解非饱和区域中的气流或理查德方程时的孔隙率。

TOP——含水层顶部的标高。对于由以下输入 LAYCON 值之一表示的层，TOP 是只读的：2、3、12、13、22、23、33、40、43。

WETDRY——是润湿阈值和标志的组合，用于指示哪些相邻单元可以导致单元变湿。如果 WETDRY < 0，则只有干电池下方的电池才会导致电池变湿。如果 WETDRY > 0，干电池下方的电池和四个水平相邻的电池可能会导致电池变湿。如果 WETDRY 为 0，则电池不能被润湿。 WETDRY 的绝对值是润湿阈值。当干电池的 BOT 和 WETDRY 的绝对值之和等于或超过相邻电池的水头时，该电池被润湿。如果 IWDFLG 不为 0 并且输入 LAYCON 值为以下之一，则 WETDRY 为只读：1、3、11、13、21、23、31、33。

VANAL - 是非饱和土壤的 van Genuchten 参数。仅当图层的输入 LAYCON 值为 43 且使用 IREALSL…0 请求真实土壤湿度函数时，VANAL 才为只读。

VANBT——是非饱和土的van Genuchten参数$。当图层的输入 LAYCON 值为 43 且使用 IREALSL…0 请求真实土壤湿度函数时，VANBT 为只读。

VANSR——是非饱和土的残余饱和水平。当图层的输入 LAYCON 值为 43 且使用 IREALSL…0 请求真实土壤湿度函数时，VANSR 为只读。

BROOK——是计算非饱和土相对渗透率的 Brooks-Corey 指数。当图层的输入 LAYCON 值为 43 时，BROOK 为只读，且实际土壤湿度函数与 BrooksCorey 渗透率关系要求 IREALSL=2 或 4。

S_STAR – 是具有 Mohanty 等人 (1997) 的双峰电导率函数的非饱和土壤的临界饱和参数 S*（参见第 2.2B 节的方程 5b）。当图层的输入 LAYCON 值为 43 且 IREALSL=7 时请求双峰电导率函数时，S_STAR 为只读。

D_STAR – 是使用 Mohanty 等人 (1997) 的生物模态电导率函数时，非饱和土壤的有效大孔隙度拟合参数（参见第 2.2B 节的方程 5b）。当图层的输入 LAYCON 值为 43 且 IREALSL=7 时请求双峰电导率函数时，D_STAR 为只读。

H_VAL – 是保湿和相对渗透率函数表格输入插值点的压头值。当层的输入 LAYCON 值为 43 且 IREALSL=5 或 6 时请求函数的表格输入时，H_VAL 为只读。

S_VAL – 是与 H_VAL 相关值的插值点对应的含水饱和度值。当层的输入 LAYCON 值为 43 且 IREALSL=5 或 6 时请求函数的表格输入时，S_VAL 为只读。

K_VAL – 是活动流动相的相应相对渗透率值，对于具有相关 H_VAL 值的插值点。

当层的输入 LAYCON 值为 43 且 IREALSL=5 或 6 时请求函数的表格输入时，K_VAL 为只读。

2.6.3 输入指令中使用的字段说明

第 3 章：裂缝井 (FWL4) 套件

3.1 概述  

HydroGeoLogic, Inc. 在美国地质调查局模块化地下水流模型 MODFLOW 中添加了新的井包，可将井筒模拟为高导率裂缝管。这些包是 MODFLOW 中现有井包（WEL1 包）的替代方案，被称为 Fracture-Well（FWL4 和 FWL5）包。 Fracture-Well 封装旨在克服与原始 WEL1 封装相关的多个问题。

首先，FWL4 和 FWL5 封装模拟通过多层孔的流动。这是因为裂缝管表示允许将与井相关的含水层单元（节点）与一维有限直径圆柱形井单元连接。此外，还严格考虑了该井规定的总提取率。来自与压裂井相关的每个单独节点的体积通量由代码自动计算，以求和到从井中提取的总流量。另一方面，原始的 WEL1 包仅在代表井的含水层节点处实现源/汇项。

因此，它要求用户确定每个节点对穿透两个或多个节点的物理井的贡献。通常估计这些节点通量贡献，以便根据贡献井节点的透射率来分配规定的提取率。这种估计本身可能是一个严重的错误，从而歪曲了物理系统。此外，对于无约束系统，根据 WEL1 包规定节点通量可能会导致计算困难。然而，FWL4 和 FWL5 软件包会自动将净井提取量（由用户输入）分配给井节点。对于每个时间步长，计算出的节点通量可能会有很大变化，就像各个层去饱和和再饱和时的情况一样。

新FWL4和FWL5套件的另一个特点是，对于超抽无侧限系统，当井中水位达到井底时，总井抽水量会自动调整。在此阶段，底部井节点的水头固定在井底条件下，并计算总抽水量，以对应于该井在这些条件下可以供应的可行数量。如果发生抽水需求小于井在井底压力条件下可供应的情况（即，在特定时期减少抽水或增加补给），井节点处的水头将反弹为了反映这些新的情况，该井将提供全部的开采需求。当水位低于井底时，以前的 WEL1 套件不会物理调整抽水速率。相反，它会继续计算到非物理领域，井节点处的水头值不断下降到井底标高以下。在其他情况下，WEL1 包会在迭代之间循环干燥和润湿泵单元，从而导致系统在计算上变得不稳定。

FWL4 和 FWL5 软件包还实现了井眼存储，这在 MODFLOW 的早期版本中被忽略。 FWL4 和 FWL5 软件包之间的差异在于，FWL4 软件包假设叠加，将井流方程添加到多孔基质的方程中，而 FWL5 软件包使用对数井函数来描述从基质块到井节点的流动。因此，如果用户有兴趣使用 FWL4 软件包准确预测局部水位下降，则仍然需要在每口井附近进行精细的面积离散化。

最后，FWL4 和 FWL5 封装仅适用于 BCF4 封装的严格可变饱和模拟选项（即输入 LAYCON = 40 或 43）。

本章的以下部分进一步阐述了 FWL4 和 FWL5 包的功能，并包括对公式和演示示例问题的简要讨论。压裂井公式的数值实现也根据文献中可用的分析解决方案进行了验证，然后是各个包的输入说明。

3.2 裂缝-井方程的公式化

通过将代表有限半径井眼筛选部分的 1-D 线元素叠加到代表多孔介质的 3-D 单元上，为 FWL4 包制定压裂井，如图 3.1 所示。 Sudicky 等人（1995）详细介绍了压裂井单元的推导和使用，并提供了示范性示例。出于一般性的考虑，我们将他们的技术应用于 BCF4 包中实现的可变饱和流量公式。该公式进一步扩展到 FWL5 包，其中使用井函数在多孔基质块和井节点之间发送流量，而不是使用叠加假设。这将在本节稍后讨论。 FWL4 和 FWL5 套件中通过管状裂缝的流量均由一维可变饱和流量方程控制，如下所示：

式中，l为沿井筛的距离（L），K′l为井的导水率（LT-1），k′rw为井内的相对渗透率，S′w为井的无量纲饱和厚度。井单元（饱和厚度与单元厚度之比），S's 为井眼特定储存量（L–1），h' 为井内水头（L）。

伪土壤函数用于定义压裂井中的不饱和流动行为，如 BCF4 包文档（第 2 章）中所述，以便井单元内的地下水位条件由垂直平衡条件准确表示。井眼的饱和水力传导率，根据层流理论，应用 Hagen-Poiseuille 方程估算，可表示为 (Sudicky et al, 1995)：

其中D是水的密度(ML–3)，g是重力加速度，:是水的粘度(ML–1 T–1)，rw是井的半径(L)。通常，对于大直径井，井水力传导率值比多孔基体的块间水力传导率大几个数量级。因此，仿真结果对使用等式(25)可能产生的误差不敏感。

如前所述，FWL4 包使用井节点与多孔基体节点的叠加。另一方面，FWL5 软件包将它们保持为单独的节点，并使用 Theim 方程井函数提供矩阵块与其关联井节点之间的流动。将 Bennett 等人 (1982) 提出的概念应用到任意层 k 的裂缝井公式中，基质块与其井节点之间的通量由下式给出： 

其中 Qmwk 是 k 层井通量矩阵； krk 为相对渗透率，Kk 和 Thickk 为基质块的水平电导和厚度，hmk 为 k 层基质块中的水头，hwk 为 k 层井节点的水头，Sk 为表皮系数re 是矩阵块的有效半径，定义为（Halford 和 Hanson，2002；Neville 和 Tonkin，2004）。

其中 △x 和 △y 是包含孔的单元格的水平网格块尺寸。

然后，FWL5 程序包会同时求解井中的水头和多孔基质块节点，这些节点彼此不同；因此，多孔基质网格块不需要在井周围进行细化，以获得更好的井流量、水头和提取条件的准确性。

FWL4 和 FWL5 软件包要求用户提供的信息包括：井位、特定井筛选的顶层和最底层的标识、井底高程、井半径、抽取 (-ve) 或注入 (+ve) ）来自井的速率（L3T–1），以及井筒特定储存值（L–1）。此外，FWL5 软件包还需要输入表面摩擦系数和在任何应力周期激活任何井时的井内初始水头。 FWL4和FWL5模块组织几何信息，使用方程(25)计算Kl，并将离散裂缝井方程(24)应用到各自的井节点。然后，在每个应力周期，将抽取或注入应用于井中最底部的筛选层。该公式将应力适当地分布在井节点之间。最后，当需要时，井节点处的边界条件（即，从规定的通量到规定的水头）的切换自动执行，如下所述。

每口井按规定的速率继续抽取，直至井内水位达到井底。井的底部标高实际上位于被筛选的最低节点层中的任何位置。当最低井节点处的测压水头对应于井的底部高程时，水头值保持固定，并且允许井产量低于需求并达到零。如果在后期，含水层系统恢复并且井产量增加到大于需求，或者如果需求下降到产量以下，则井恢复以规定的抽水速率抽水。

输出文件中报告了总井产量和每层对井产量的贡献。在以下部分中，提供说明性示例来验证和演示 FWL4 和 FWL5 软件包的应用。

3.3 验证示例

通过将数值结果与 Theis (1935) 以及 Papadopulos 和 Cooper (1967) 获得的解析解进行比较，验证了压裂井包 (FWL4) 选项，该解析解针对具有恒定流量井的均质含水层中的水位下降。与 Theis 解不同，Papadopulos 和 Cooper 的解考虑了井筒储存。

因此，进行了两种模拟，一种有井眼存储，另一种没有井眼存储，以测试 FWL4 中使用的数值方案。

Sudicky 等人 (1995) 描述了验证问题。在此问题中，半径为 0.1 m 的全穿透井位于尺寸为 1,000 m x 1,000 m x 5 m 的均质各向同性含水层的中心。该井的抽水速度为 10 –3 m3 /s。

横向边界规定了 20 m 的恒定水头，并且假定含水层的顶部和底部边界是不渗透的。含水层的水力传导率为10–4 m/s，比容为10–4 m–1。使用 FWL4 软件包执行瞬态仿真。

建模方法和结果

该域被离散为 5 层，每层 1 英尺厚，具有 35 行和 35 列可变网格间距（从井附近的 111.1 英尺到 0.1 英尺）。含水层的顶部和底部高程分别设置为 0 和 –5 m。 MODFLOW 要求的 VCONT 值（垂直导水率除以任意两层之间的层间距离）计算为 10–4 s–1。含水层的主要蓄水系数（具体蓄水时间块厚度）指定为10 – 4。含水层侧面规定恒定水头20 m，并假设初始水头均匀为20 m。 0.1 m 半径井位于区域中心（第 18 行、第 18 列），允许以 10–3 m3 /s 的速率泵送 1000 秒。由于井已完全贯通，井顶和井底分别位于1层和5层，井底标高指定为–5 m。为了将井压降的时间变化与 Theis 和 Papadopulos 的解析解进行比较，在特定井筒储水量分别为零和 0.1 m-1 的情况下进行了两次瞬态模拟。对于这两种模拟，PCG4 软件包用于求解矩阵方程。图 3.2 绘制了井面的压降与时间曲线以及 Theis 和 Papadopulos 解析解。如图所示，使用FWL4 Package获得的MODFLOWSURFACT/MODHMS结果与各自的解析解非常吻合，从而验证了FWL4 Package中采用的压裂井公式。

3.4 应用示例

提供了两个示例来演示FWL4包的应用，第三个示例演示了FWL5包的验证和应用。第一个示例说明了完全 3D 情况，考虑了从无承压含水层抽水的井。第二个示例说明了准 3-D 情况，考虑了从由限制床分隔的三个含水层中抽水的井。在这两种情况下，井都被故意过度抽水，并在第二个压力期恢复含水层。作为对井边界条件实施严格性的全面检查，系统的响应也使用不调用 FWL4 包的独立模拟进行评估。对于全 3-D 情况，从压裂井模拟获得的井处水头值被强加为验证模拟的时间相关水头条件，并将结果与从原始模拟获得的结果进行比较，其中 FWL4 包被使用了。最后一个示例说明了使用带有准 3-D 设置的 FWL5 软件包从由约束床分隔的两个含水层进行井抽水。该问题的受限情况与解析解进行比较进行验证，该问题的非受限情况则展示了包装的干燥能力。

3.4.1 问题 1——无承压含水层抽水的全三维模拟

该问题考虑了图 3.3 所示的 300 英尺厚的均质、各向同性无承压含水层。建模域是一个 75,000 英尺 x 75,000 英尺的正方形，含水层的顶部和底部高程分别为 +50 英尺和 –250 英尺。整个区域含水层的初始地下水位为零。图中左（西）边界为恒水头边界，其余边界为无流边界。一口直径为 0.4 英尺的全穿透井和 10 –6 英尺 –1 的井眼特定存储位于该区域的中心。 10 ft 3 /s 容量的泵运行 2.5 年后，停止从井中抽水。 （见图 3.4。）然后允许含水层恢复 4 年。含水层参数包括：

建模方法

该域被离散为 3 层、15 行和 15 列（△x = △y = 5,000 ft 和 △z = 100 ft）。 VCONT（垂直导水率除以两个相邻层之间的层间距离）的值计算为 10 –6 s–1。为每层指定 10-4 的恒定主存储系数（特定存储时间块厚度）。

顶部和中间离散层的西侧规定恒定水头为零。

直径 0.4 英尺的井位于区域中心（第 8 行、第 8 列），并在所有 3 层中进行屏蔽直至含水层底部。前 2.5 年的提款需求为 10 ft 3 /s。请注意，当含水层无法提供指定的抽水速率（即，当井中的水位到达井底时）时，Fracture-Well (FWL4) 包会自动调整抽水速率。在第二个压力期，抽水会停止 4 年，以便含水层恢复。由于井已完全贯通，井的顶部和底部分别位于第 1 层和第 3 层，井底标高指定为 –250 英尺。

使用 BCF4（输入 LAYCON=43 的可变饱和流选项）和 FWL4 封装执行瞬态仿真。 PCG4 软件包用于求解矩阵方程。用于模拟的输入数据文件列在附录 C 中。

结果与讨论

井内水位和总井通量随时间的变化分别如图 3.5 和图 3.6 所示。如图 3.5 所示，井内水位在 2.15 年内从初始值（0）逐渐下降至井底（-250 英尺）。水位保持在井底直至停止抽水。在第二个（恢复）期间，井中水位从-250 英尺上升到-17.13 英尺。总井通量的相应变化可以在图3.6 中观察到。在第一个应力期间，从井中的总抽油量从 10 英尺 3 /秒减少到 7.63 英尺 3 /秒。当抽水停止时，总井通量减少至零。

采用图 3.5 中描述的瞬态水头分布，使用上述相同的水文条件执行规定的井口模拟，但不调用 FWL4 包。本次仿真的目的是验证之前使用 FWL4 封装进行仿真时获得的结果的准确性。从之前的模拟中获得的瞬态水头值是针对无承压含水层第 3 层、第 8 行和第 8 列（即底部井节点）的单元指定的（图 3.3）。

为了适应规定的瞬态水头分布，单元被视为具有大界面水力传导率（109 ft/s）的一般水头边界节点。

图 3.6 比较了验证井模拟和压裂井模拟运行获得的模拟结果。两组结果非常吻合，从而验证了 FWL4 套件对规定抽水量井进行严格处理的准确性。

3.4.2 问题2——从三含水层系统中抽取的准三维模拟

该问题涉及由两个 10 英尺厚的零渗透性限制层（含水层）分隔的三个含水层，如图 3.7 所示。请注意，使用非零泄漏可以轻松处理泄漏的弱透水层。建模域是一个尺寸为 75,000 英尺的正方形。三个含水层的顶部和底部高程如下：

整个域内所有含水层的初始测压水头均为零。图中左（西）边界为零水头常水头边界，其余边界为无流边界。一口直径为 0.4 英尺的全穿透井和 10 –6 英尺 –1 的井眼特定储存量位于该区域的中心。以 5 ft 3 /s 的速度运行三年后，将停止从井中抽水（图 3.8）。然后允许含水层恢复 9 年。三个含水层的水力参数如下：

建模方法

域中的每个含水层都被视为节点层，并且限制层（含水层）中的垂直渗漏为零。每个含水层均匀分为 15 行和 15 列 ()x = )y = 5,000 ft)。第 1 层、第 2 层和第 3 层分别厚 90 英尺、40 英尺和 30 英尺。通过将 VCONT（含水层之间的垂直渗漏）的值指定为零，可以纳入限制层对流量的影响。每个含水层指定恒定的初级存储系数（特定存储时间块厚度）0.001。在顶部和中间离散层的西侧规定了恒定的零水头。一口直径为 0.4 英尺的井位于域的中心（第 8 行、第 8 列），前 3 年允许以 5 英尺 3 /秒的速率抽水。请注意，裂缝井包会自动将总抽取量分配给各个含水层。如果总需求不可行（即井内水位达到井底时），则将井回抽率调整至含水层系统能够供给的值。在随后的 9 年内停止抽水以恢复含水层。由于该井已完全贯通，井的顶部和底部分别位于第 1 层和第 3 层，井底标高指定为 -130 ft。井筒具体储存量为 10–6 ft–1。

使用 BCF4（输入 LAYCON=43 的可变饱和流选项）和 FWL4 封装执行瞬态仿真。 PCG4 软件包用于求解矩阵方程。除了总井通量之外，FWL4 软件包还打印每个时间步长每个井节点的通量和水头值。

结果与讨论

井节点水头随时间的变化如图 3.9 所示。每个含水层对井的通量贡献如图 3.10 所示。最初，1 号含水层的贡献占井取水总量的 90% 以上（图 3.10）。当井节点的水头下降到含水层底部标高时，这种贡献就会下降。

然后，含水层 2 和含水层 3 开始对所施加的抽水压力产生更多影响。

当井中的水位到达井底时（t = 2.3 年），井的总产量开始降至低于 5 ft 3 /s 的指定需求，并在 3 年时降至 1.3 ft 3 /s。当在第二个胁迫期（3年后）停止抽水时，含水层1和含水层2继续向井供水，而含水层3从井取水。

图 3.11 至 3.14 分别显示了 2.0、3.0、4.5 和 7.5 年模拟的井节点处测压头的横截面剖面。三个含水层中的测压头显示为沿第 8 行穿过井的东西向垂直剖面（见图 3.7）。如图 3.11 所示，在 t = 2 年时，顶部、中部和底部含水层的井测压水位分别为 –39.98 英尺、–89.97 英尺和 –116.00 英尺。图 3.12 表示第一个应力期结束时含水层的状态。停止抽水一年半后（t=4.5 年），由于通过井眼从两个上覆含水层获得补给，观察到含水层 3 显着恢复（图 3.13）。图 3.14 显示了第二个应力期（t=7.5 年）中期的测压水头，表明井中水位从 –130 英尺的井底水平上升至 –58.25 英尺。随着测压水头的增加水平，含水层 3 转换为完全承压含水层。含水层 2 由来自含水层 1 的水通过井以及从恒水头边界向东补给。在第二个应力期结束时（t=12 年），井中水位上升至 –23.21 英尺，含水层 2 完全受到限制（另见图 3.9）。

3.4.3 问题 3——使用 FWL5 软件包从两含水层系统中抽取的准三维模拟

该模拟的问题设置改编自 Neville 和 Tonkin (2004)，并考虑了跨两个含水层的开放井。该区域空中呈正方形，边长为 47,000 英尺，中间有一个抽水井，每天抽水 62,840 英尺。沿域的横向边缘应用无流边界条件。该问题考虑了两种情况——第一种是受限模拟情况，与 Papadopulos (1966) 的解析解进行比较以进行验证；第二种是无限制模拟情况，演示了在良好干燥条件下的应用。

两种情况下含水层的顶部和底部高程如下。

顶部含水层的初始水头统一设置为 10 英尺。底部含水层的初始水头统一设置为 30 英尺。顶部含水层的水力传导率和比储存值为 K=100 ft/d，S=0.0001，并且对于底部含水层，K=400 ft/d，S=0.0001。该井完全穿透两个含水层，井眼半径为 0.5 英尺。

建模方法

使用 FWL4 和 FWL5 软件包解决了这个问题。域中的每个含水层都被视为单独的节点层，通过限制床的垂直泄漏为零。 FWL5 模拟使用 100 列和 100 行的均匀间隔网格，网格间距为 470 英尺。井中的初始水头与相邻层中的初始水头相同。

FWL4 模拟使用可变间距网格，邻近井的网格间距更细。第一种情况（受限）模拟持续 1.2 天，以排除稍后横向边界效应可能干扰解析解结果的比较。第二种情况（无侧限）模拟持续160天，以使井内的头达到井底条件。

结果与讨论

第一个（封闭）模拟情况下井内水位和通量的时间变化分别如图 3.15 和 3.16 所示。 FWL4 和 FWL5 案例的结果与解析解相匹配，从而验证了 FWL4 和 FWL5 包。

请注意，图中的比较是使用 MNW 包进行的模拟，它与解析解的效果非常好（Neville 和 Tonkin，2004 年）。此外，FWL4 包通过井周围的精细网格提供了两个含水层之间正确的水头和通量分配。图 3.17 和 3.18 显示了第二个（无侧限）模拟情况下井内水头和通量的时间变化。 FWL4 和 FWL5 包的结果非常匹配，表明具有不同井近似值的两个包提供了相似的结果。此外，两种包装都能成功地应对良好的干燥条件。

3.5 FWL4 和 FWL5 包的输入指令

3.5.1 概述

压裂井（FWL4 和 FWL5）包是井（WEL1）包的替代品，模拟高导管状压裂元件的井。 FWL4 和 FWL5 软件包通过将与井相关的所有节点与具有高电导率的一维有限直径圆柱形井单元相连来模拟井。从概念上讲，压裂井是通过将代表井眼筛选部分的这些井元素叠加到代表多孔基质的 3-D 单元上来在 FWL4 包中制定的。从概念上讲，压裂井在 FWL5 包中被公式化为代表井眼筛选部分的单独井元素，这些井元素根据井函数（蒂姆方程）从代表多孔基质的相关 3-D 单元接收流量。与 WEL1 封装相比，FWL4 和 FWL5 封装的优点包括：

• 代码自动计算与压裂井相关的每个单独节点的抽取量，以求和到用户指定的井的总抽取量。因此，无需事先估计每个节点对穿透两个或多个节点的物理井的贡献。

• 对于超抽含水层系统，井中的水位保持在井底水平，井总抽水量自动调整为井在该阶段可以供应的数量。对于多个含水层，每个含水层对油井的贡献会自动调整为压力过大的系统的产量。如果出现抽水需求小于含水层在井底压力条件下可产出的情况（即，在特定抽水周期内减少井抽水量或增加补给量），则井中的水位会因此而测压井节点处的水头将反弹以反映新的应力。 FWL4 和 FWL5 套件的这一功能甚至对于单层系统或完全从单层抽水的系统也很有用。

• FWL4 和FWL5 套件用于井眼存储。应该注意的是，在抽水区域仍然需要对模拟域进行更精细的面积离散，以准确地描述 FWL4 程序包井附近的水位下降，但 FWL5 程序包不需要这样做，因为连接网格的对数函数FWL4 配方中未合并孔细胞块，而是使用叠加。

FWL4 封装只能与 BCF4 封装的新 LAYCON 选项一起使用（即输入 LAYCON = 40 或 43）。

3.5.2FWL4输入

Fracture-Well (FWL4) 包的输入是从 IUNIT (21) 中指定的单元读取的。要调用 FWL4 包，模型网格中的所有层必须使用输入 LAYCON 值 40 或 43。

FWL4 包的输入数据主要包括井位、每口井被屏蔽的顶层和底层、井的底部高程和半径、抽取 (-ve) 或注入 (+ve) 速率以及井眼具体信息存储（L–1）。每口井的信息均在给定压力期间的一条记录中指定。应力期内的记录数量取决于域中的井数。

MXFWEL——模拟中的最大井数。其中包括未堵塞的井，这些井可能不会抽水，但仍充当层之间的垂直连接。

IFWLCB——是一个标志和一个单元号。

如果 IFWLCB > 0，则它是每当设置 ICBCFL（请参阅输出控制）时将记录孔的逐个单元流动项的单元编号。

如果 IFWLCB = 0，则不会打印或记录孔的逐个细胞流动项。

如果 IFWLCB < 0，则每当设置 ICBCFL 时都会打印井注入速率。

（如果 IFBLCB=-1，每个 FWL 单元的通量贡献也会被打印。）

KR = FCONST *（半径）2

其中 FCONST =ρg/8μ（参见压裂井公式的方程 2）

半径 = 井筛半径 [L]； 

ρ = 水的密度 [M/L3 ]； 

g = 重力加速度 [L/T 2 ]；和μ=水的动力粘度[M/LT]。

FCONST 的值取决于准备输入数据集时使用的长度和时间单位。下表列出了不同长度和时间单位的 FCONST 值：

FCONT 尺寸值[1/LT]

如果仿真中使用的长度和时间单位与上表中给出的不同，则必须使用上面给出的关系计算 FCONST 的等效值。

IFREWL4——是以自由格式读取FWL4文件的数据项2和3的标志。

如果 IFREWL4 = 0，则按照输入指令中所述的默认格式读取 FWL4 文件。

如果 IFREWL4 = 1，则以自由格式读取 FWL4 文件的数据项 2 和 3。

MUTWL4——是一个禁止在输出文件中打印 FWL4 输入的标志。

如果 MUTWL4 = 0，则将每个应力周期的 FWL4 输入读取打印到输出文件。

如果 MUTWL4 = 1，则抑制在每个应力周期将 FWL4 输入打印到输出文件，但仍会打印警告。

ITMP——是给定压力期内的井数（包括所有未堵塞的井）。

注：除非钻探新井或堵塞旧井，否则 ITMP 保持不变且等于 MXFWEL。换句话说，未放电/补给的现有井必须包含在每个应力周期的数据中，并将其 Q 值设置为零。

LSTART——是包含孔筛顶部的模型单元的层号。

LEND——是包含井底的模型单元的层号。

行——是包含孔的模型单元的行号。

列——是包含孔的模型单元的列号。

高程——是井底的高程。从物理上讲，该标高应位于井的底层 (LEND) 内。

半径——是井的半径，其单位与模拟中使用的任何其他长度单位相同。将井的半径设置为零可以有效地断开井裂缝单元，因为裂缝井电导率（方程（1）和（2）中的K R ）变为零。 （由于井体积为零，井的储存量也变为零）。

Q——是进入井中的体积通量（L3 T–1）。正值表示井中注入水，负值表示井中排出水。

储存量——井眼特定储存量（S′s），由以下公式给出：S′s=（1/Ls）/（rc2/rw2），假设含水层饱和模拟，其中Ls为筛网总长度，rc为直立水管（或套管）半径，rw为筛网半径。请注意，井眼储存的影响对于大直径井来说相当大，特别是在模拟的早期阶段。通过将存储期限设置为零，可以忽略井眼存储效应。对于无侧限模拟，因井去饱和而产生的井筒存储自动纳入方程（4）中的井。

3.5.3 输入指令中使用的字段说明

 3.5.4 断裂井 (FWL5) 封装输入

FWL5 包的输入是从 IUNIT (37) 中指定的单位读取的，按照流量模拟的要求（在文件中默认扩展名 WL5）。首先提供所有井的输入井属性，包括井位置、每口井被屏蔽的顶层和底层、底部高程、井半径和井眼特定存储。

在模拟开始时，所有孔均设置为非活动状态。然后，对于每个应激期，可以激活或随后停用任意数量的井。对于在任何压力时期激活的每口井，都会寻求井的抽水速率的附加信息（抽水为负，注入为正）、表示集肤效应电导的包含/读数的标志以及指示输入的标志井中的初始水头条件。如果这些标志处于活动状态，则需要额外的相应输入。请注意，每当将非活动状态激活时，都需要输入井中的初始水头，以防止非物理值 (HDRY) 表达井的压力状态。 FWL5 软件包执行这些检查，如果检测到在没有关联初始水头输入的情况下激活了井，则将中止模拟，并在输出列表文件中显示适当的消息。此外，为每个物种的传输模拟指定了注入流体的浓度。请注意，传输所需的额外输入变量显示在阴影框中。

请注意，LSTART 和 LEND 是井 LWELL 的起始层和结束层，并通过数据项 2 为所有井定义。

** 仅当 INHD = 1 时才输入第 6 项（该孔的 ITMP 和 IBOUNDF5 也必须大于零）。

6. 数据： 格式：HDFWL5 (K), K=LSTART, LEND 8F10.3 ** 仅当 INHD = 1 且 ITRAN 0 时才输入第 7 项。（对于 ≠ 此孔，ITMP 和 IBOUNDF5 也必须大于零）。

7. 数据： 格式：CONFWL5 (K, LSPEC), K=LSTART, LEND 8F10.3 注：对于模拟的 LSPEC = 1、NSPECI 污染物种类，应重复第 7 项（对于每个种类）。

3.5.5 FWL5 输入指令中使用的字段说明

MXFWEL5——模拟中的最大井数。其中包括可能不活跃（汇或源）但仍充当层之间垂直连接的井。

IFWL5CB——是流量模拟的标志和单元号。

如果 IFWL5CB > 0，则它是每当设置 ICBCFL（请参阅输出控制）时将记录孔的逐个单元流动项的单元编号。文件的默认扩展名为 CW5。

如果 IFWL5CB = 0，则不会打印或记录孔的逐个细胞流动项。

如果 IFWL5CB < 0，则每当设置 ICBCFL 时都会打印井注入率。

（如果 IFWL5CB=-1，每个 FWL 单元的通量贡献也会被打印。）

FCONST5——用于计算井的水力传导率的乘数，计算公式为： 

 KR = FCONST5 *（半径）2

其中 FCONST = ρg/8μ

半径 = 井筛半径 [L]；

 ρ = 水的密度 [M/L3 ]；

g = 重力加速度 [L/T 2 ]；和

μ= 水的动力粘度[M/LT] 。

FCONST5 的值取决于准备输入数据集时使用的长度和时间单位。下表列出了不同长度和时间单位的 FCONST5 值：

FCONT5 Value in dimensions[1/LT]

FCONT5 尺寸值[1/LT]

第 4 章：补给渗流面边界 (RSF4) 条件包

4.1 概述

MODFLOW 代码使用 Recharge (RCH1) 包适应地下水补给边界条件。然而，RCH1 包在表示无约束系统方面有一定的局限性。如果无承压含水层饱和到地表，含水层吸收补给的能力就会降低，剩余的水会作为地表径流流失。 RCH1 包无法处理这种情况。相反，它继续为含水层提供补给（如在封闭系统中），水头不断上升至远高于地面。 MODFLOW-SURFACT 的新 RSF4 软件包旨在消除这种物理上不切实际的情况。如果地下水位低于用户规定的水池（积水）高度，RSF4 套件允许将所提供的补给到地下水系统中。

水池标高对应于无积水时的地面标高。

如果在任何时候地下水位达到水池高度，模拟只允许进行尽可能多的补给以维持规定的水池条件。剩余充值不被接受到模拟域中。模拟的输出反映了由于系统饱和至积水高度而导致系统体积补给通量 (L3 /T) 的减少。请注意，RSF4 包与 Streamflow-Routing 包（STR1 或 SFR1）分离，因此，该剩余补给不会以任何方式路由到规定的流中。

此外，对于 MODHMS 模拟，无侧限渗流选项与该域处于活动状态的地表水流域相冲突，因此应在这些位置设置非常高的积水标高，以防止该包与地表水流的干扰，如果无侧限渗流渗流选项用于 RSF4 包中。

新 RSF4 软件包的一项附加功能可以方便地输入瞬态补给/降水，而与 MODFLOW 的“应力周期”概念无关。使用此选项时，可以通过单独的文件输入充值时间序列，该文件确定域上的充值如何随时间变化，而不管压力周期设置如何，以提供输入便利，并且通常在充值时提供较小的输入文件大小变化很快，并且经常在模拟周期内变化。

除了地表补给-排水条件外，RSF4 软件包还可用于通过规定渗流面边界节点的高程来模拟渗流面边界条件。因此，边界处出现零通量，直到地下水位达到渗流面高程。如果地下水位达到渗流面高程，则含水层排水以维持渗流面高程（即大气压力）条件。如果在模拟过程中地下水位下降到渗流面高程条件以下，则渗流面边界处将保持零通量。

RSF4 软件包的无侧限渗流选项仅适用于 BCF4 软件包的严格可变饱和模拟选项（即输入 LAYCON = 40 或 43）。下面介绍 RSF4 包的功能，包括对公式和演示示例问题的简要讨论。

4.2 RSF4 边界条件包的制定

4.3 快速、瞬态充电条件验证和应用示例

 4.4 验证和应用示例

 4.4.1 问题 1——无承压含水层中平行排水沟的流量

4.4.2 问题2——方形路堤渗漏

4.5 Rsf4 封装输入指令

4.5.1 RSF4 输入

4.5.2 瞬态充电时间序列文件

第 5 章：自适应时间步进和输出控制 (ATO4) 包

5.1 概述

在 MODFLOW 代码中，时域（用于瞬态流模拟）使用后向差分公式进行离散化，时间步值通过指定以下参数预先确定：时间步数 (NSTP)、应力周期长度（PERLEN），以及使用几何级数增加时间步长的乘法因子（TSMULT）。下面列出了使用这种方案的限制。

• 如果矩阵求解器或非线性迭代方案（对于无约束情况）未能在给定时间步内收敛，则计算将中止。

• 即使获得了收敛的数值结果，预选的时间步长对于整体解决方案来说可能效率低下。

• 必须保持时间步长的几何级数。结果，用户在分配特定时间值来检查流动系统方面缺乏灵活性。

因此，只能在时间级数上的时间值请求输出，这需要由用户提前确定。

为了克服上述困难并提高解决方案的效率，开发了自适应时间步进和输出控制包（ATO4）。自适应时间步长方案根据给定计算的系统的预期非线性来选择时间步长大小。如果预期的非线性不显着，则选择较大的时间步长以积极地推进模拟。如果预期的非线性很严重，则选择较小的时间步长以确保该时间步长的收敛。如果解在给定的时间步长内未能收敛，则进一步减小时间步长，并重复求解。

对于自动时间步进方案，模拟结果的时间值是未知的；因此，ATO4 软件包还包括用于输出控制的新模块，该模块确定所需输出的时间值。调用 ATO4 包后，不再需要使用 MODFLOW 之前的输出控制 (OC) 包。请注意，原始 MODFLOW 的输出控制实现起来很繁琐且麻烦。另一方面，新 (ATO4) 方案的输出控制则简单明了。用户只需输入任何给定压力期间所需输出的时间值（以先前使用的时间单位），以及所请求输出的详细信息（即，头、回撤、质量预算和预算中的一项或全部）。细胞间的流动项是否被输出）。 ATO4 包自动识别输出的选定时间值，并调整时间步长以在这些打印时间执行计算。因此，例如，用户可以在 15 个月的压力期模拟的 10、11 和 15 个月后检查结果。 ATO4 包不仅在整个模拟过程中提供高效的时间步进，而且还确保在 10、11 和 15 个月（请求输出的时间值）执行计算。

最后，ATO 包提供了从非零起始时间值重新启动仿真的选项。 MODFLOW 的 OC 包不允许重新启动，并且开发新的输入文件来执行重新启动可能会变得极其繁琐，特别是在涉及多个压力周期时。

5.2 ATO4 包中使用的自适应时间步进的制定

5.3 验证示例

5.4 ATO4 封装的输入指令

5.4.1 概述

5.4.2 ATO4 输入

5.4.3 输入指令中使用的字段说明

第 6 章：预条件共轭梯度 (PCG4) 包

6.1 概述

MODFLOW 代码生成一个以有限差分形式描述地下水流系统的方程组。这组代数方程通常通过使用切片连续过松弛 (SSOR) 或强隐式过程 (SIP) 包来求解。对于复杂的现场模拟，这些求解例程可能缺乏效率和鲁棒性。 Hill (1994) 提供了预条件共轭梯度 (PCG2) 包，为 MODFLOW 提供了替代求解器。

PCG2 中的预处理方案是 Saad (1985) 提出的最小二乘多项式预处理器或 Meyer 等人 (1989) 的最优切比雪夫多项式预处理器。选择这些方案主要是出于计算机存储方面的考虑，并且在大规模现场研究中通常表现不佳。认识到这一点，我们合并了 PCG4 包，其中包含基于部分 LU 分解的 PCG 求解器作为预处理器。计算机存储要求比以前使用 MODFLOW 的求解器要求的要多。然而，PCG4 求解器简单、稳健且高效。该求解器已用于许多不同规模和复杂程度的地下水建模项目。本章的其余部分致力于 PCG4 包的制定和输入其使用说明。

6.2 标准PCG方案的制定

6.3 PCG4 封装的输入指令

6.3.1 PCG4 输入

6.3.2 输入指令中使用的字段说明

第 7 章：带回溯 NRB1 包的牛顿-拉夫森线性化

7.1 概述 . 

以水或空气为活性相的严格非饱和区流动模拟通常是高度非线性的，需要小时间步长和多次迭代才能解决实际情况。对于某些情况，一步稳态模拟可能无法与无约束情况下的 MODFLOW 皮卡德迭代方案收敛，或者对于非饱和流（空气或水）情况下的 MODFLOW-SURFACT/MODHMS（第 2 章）的修改皮卡德方案可能无法收敛，或用伪土壤函数确定无承压含水层中的地下水位高程，其中不需要严格的非饱和流模拟。提供二次收敛的Newton-Raphson方案可以极大地缓解收敛困难。这里讨论的 NRB1 包将 Newton-Raphson 方案与回溯算法集成到 MODFLOW 中以控制步长，以用于非线性情况，尤其是当 Picard 迭代过多或完全失败时，这会受益。

NRB1 包包括 Newton-Raphson 线性化和回溯方案，以稳定牛顿迭代。该公式被设计为与当前在每次迭代中求解所有节点处的头的 MODFLOW 方案兼容。

此外，该公式不需要对 MODFLOW 当前模拟的任何边界条件进行特殊处理。最后，回溯方案限制了任何迭代中残差的增加，而欠松弛技术有助于迭代之间解决方案的振荡行为。这些问题将在本章的其余部分讨论

7.2 牛顿拉夫森方案的制定

7.3 NRB1 封装的输入指令

第 8 章：观察节点 (OBS1) 包

8.1 概述

MODFLOWSURFACT/MODHMS 中包含观察节点包，以简化域内任何节点处瞬态模拟突破曲线的后处理。流量模拟输出显示模拟所请求节点处的压力-时间关系，传输输出包括模拟的每种污染物的浓度-时间关系。创建一个包含突破输出的单独输出文件，可以轻松对其进行后处理以提供突破曲线

8.2 OBS1 包输入指令

8.2.1 OBS1 输入

8.2.2 输入指令中使用的字段说明

8.3 OBS1 包的输出

第 9 章：参考资料