第二章：块中心流（BCF4）模块

2.1 通用

MODFLOW代码采用块中心有限差分方法来求解地下水流方程。流域被离散化为行、列和层，以使每个节点代表一个多孔材料的矩形块，被称为一个单元格。结果有限差分网格中的节点表示无流、可变水头或恒定水头单元格，与单元格相关联的任何水力特性都是相对于相应的节点指定的。
对于非围限流模拟，原始MODFLOW代码的BCF1模块（McDonald和Harbaugh，1988年）将变水头单元格转换为无流单元格。这导致了流域的部分永久排除在模型模拟之外。如果MODFLOW不能在电位水平恢复时将这些单元格重新转换为可变水头单元格，则可能导致代码产生误导性或错误的模拟结果。
McDonald等人（1991年）尝试将无流单元格重新转换为可变水头单元格，以允许需要时重新饱和单元格。他们的再润湿方案在BCF2模块（块中心流模块，第2版）中实施。该方案利用相邻单元格的水头来确定是否将无流单元格重新转换为可变水头单元格。不幸的是，BCF2的重新润湿选项在重新湿润时或者当提取使相应单元格干涸时容易出现收敛和稳定性问题（McDonald等人，1991年）。这是因为用于重新激活干燥单元格的特设程序可能严重违反流动和质量守恒原理。此外，当相邻单元格的头依赖性渗透率降为零时，代码中使用的调和透过率平均方案在物理上变得不一致。Goode和Appel（1992年）记录了与BCF2重新湿润包相关的一系列问题以及一系列缓解方案，并向MODFLOW中添加了BCF3模块，以提供用于平均透过率的替代程序。鉴于上述重新湿润程序遇到的严重困难，需要一种严格的方法，以满足流动连续性要求，并允许非围限
层中水位的自由移动，而无需任何形式的强制转换。HydroGeoLogic公司通过使用3-D可变饱和流动公式并自动生成伪土壤水保留函数来解决这一创新问题，从而为此提供了一种创新方法，以将不饱和流动问题减少到在一个单元格中寻找水位（即压力头为零或大气压的高程）。该公式已被设计为提供水位的准确划分，并捕获非围限系统对抽水和补给的延迟响应。模拟的数据需求保持不变，因为使用伪土壤关系而不是真实土壤的本构关系。我们方法的关键优势包括：

一种严格处理三维流动（即，不需要Dupuit假设）的方法，不需要打开和关闭单元格。

一种具有良好收敛性和稳定性特性的强大数值方案，以及处理单元格完全脱水和重新饱和的所需能力。利用有限差分方法，将可变饱和公式与伪土壤函数编程到现有的MODFLOW的BCF3模块中，并将新版本称为BCF4模块。BCF4模块的新选项将块间导水性计算为块水力导水性、相对渗透率和平均流动面积的加权调和平均值的乘积。BCF4模块中使用可变饱和公式与伪土壤函数以及其应用将在接下来的章节中讨论。
BCF4模块还提供了使用轴对称圆柱坐标系进行轴对称分析的选项。该选项避免了通过笛卡尔（x，y，z）坐标进行完全三维公式的使用，从而在计算工作中节省了大量资源。轴对称分析选项适用于某些情况（即，单井流），其中轴对称几何形状是可接受的。
在BCF4模块的开发中，保留了MODFLOW的模块化性质和输入数据准备的结构。新功能的数据需求与原始MODFLOW相同，不需要额外的参数值。为了保持向后兼容性并允许用户比较各种计算方案，BCF4模块还包含了以前BCF3和BCF2模块的所有选项。在本章中，用户通过对公式、验证和示范性示例问题的简要讨论，系统地介绍了BCF4模块的新功能。

2.2A 利用可变饱和流动公式来建立非围限流动方程

水在可变饱和系统中的三维运动可以表示为（Huyakorn等人，1986年）

式中: x、y和z是笛卡尔坐标（长度单位）；
Kxx、Kyy和Kzz分别是沿x、y和z轴的水力导度的主要分量（长度除以时间）；krw是相对渗透率，它是水饱和度的函数；
h是水力头（长度单位）；
W是单位体积的体积通量，表示水的源和/或汇（时间的倒数）。
φ 是可排水孔隙度，取为比流量 Sy；
Sw是水的饱和度，是压力头的函数；
Ss 是多孔材料的比存储（长度的倒数）；
t 是时间（时间单位）。
对于完全饱和介质（即，Sw = 1.0），相对渗透率为单位，方程（1）简化为：

方程（2）是 MODFLOW 开发中使用的基本地下水流方程（McDonald 和 Harbaugh，1988）。 
因此，可变饱和流量方程简化为地下水位以下和受限系统中的传统地下水流量方程。 拟土关系用于定义函数关系：Sw = Sw (ψ) 和 krw = krw (Sw)，其中 ψ 是压头，定义为 ψ = 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为残余水饱和度。 当选择参数 n 使得 n = 1 + 2/γ 时，Brooks-Corey 表达式会生成与 van Genuchten 函数类似的曲线（van Genuchten，1980）。
按照 Mohanty 等人 (1997) 的建议，通过包含分段连续水力传导率函数，可以进一步
增强模拟能力，从而包含系统中优先流动路径的双峰描述。 因此，如果系统中的饱和度小
于临界值S*，则使用等式(4b)的van Genuchten函数。 然而，当饱和度高于 S* 时，将应用指数电导率曲线，该曲线会在接近饱和时迅速增加电导率，如下所示：

其中 S* 是临界点或断点饱和值，其中流量从毛细管发生变化
主导到非毛细管主导，反之亦然，kr* 是对应于 S* 的相对渗透率，δ 是代表有效大孔隙度或其他有助于非毛细管主导流动的结构特征的拟合参数。 还提供了表格输入选项，可以适应其他功能形式。
压头（ψ = 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.3A 非承压流动验证示例

BCF4 包通过将数值结果与 Neuman (1972) 的解析解进行了验证，该解析解针对无承压含水层中的井流，该含水层对抽水表现出显着延迟的产量响应。
所考虑的问题涉及以 3658 ft3/d 的速率从均质且各向同性的无承压含水层中抽水的完全穿透井。 含水层的初始饱和厚度为 40 ft，井在含水层厚度较低的 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 英尺的正方形，
井位于中心。 为了利用对称性，我们仅离散流域的一个象限并考虑井流量 Q 的四分之一。
3D 网格由 20 层组成，每层厚 2 英尺，有 21 行和 21 列。 水平网格间距（Δx 和 Δy）从 0.5 英尺到 100 英尺不等。VCONT（垂直导水率除以两个相邻节点层之间的层间距离）的值计算为 10.5 d-1。 每层的恒定主存储系数（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 定义为：