MODFLOW-SURFACT卷一：流动模块

摘要

MODFLOW-SURFACT是由HydroGeoLogic,      Inc.开发的完全集成的地下水流和溶质运移代码，基于广受欢迎的美国地质调查局模块化三维（3-D）地下水流模拟代码MODFLOW。本卷（代码文档的第I卷）描述了添加到MODFLOW中的新流动模块，以增强其地下水流模拟能力和计算鲁棒性。详细说明了新增增强功能的物理和数学概念，并讨论了这些概念如何实现到 MODFLOW的模块化结构中。提供示例问题以验证代码，并使用户熟悉其应用。预期读者熟悉原始MODFLOW文档（McDonald和Harbaugh，1988）。

MODFLOW使用块心有限差分方法来模拟地下水流。可以执行受限和无限制层的全或准三维模拟。通过其新的流动包，MODFLOW-SURFACT增强了执行无限制模拟的方案，以严格模拟含水层的除湿/再湿润，并克服了先前版本的MODFLOW遇到的数值困难。MODFLOW- SURFACT提供了使用轴对称几何离散化域的选项，以便有效地模拟抽水试验、水库抽水/恢复试验等。通常允许的外部应力包括恒定头、恒定通量、面积补给、蒸发蒸腾、排水和河流。此外，MODFLOW-SURFACT提供了一个严格的井提取包、无限制的补给边界条件和渗透面边界条件。最后，MODFLOW-SURFACT包括自适应时间步长和输出控制程序的选项，以及额外的预处理共轭梯度（PCG）解包。MODFLOW-SURFACT的第二版包含额外的功能，包括严格的饱和-非饱和湿度运动模拟能力、空气流动模拟能力和用于改进鲁棒性的牛顿-拉弗逊线性化包。MODFLOW-SURFACT的新流动包是用FORTRAN 77编写的，并使用CompaqVisual Fortran Professional Edition 6.6®编译器编译。建议最低配置为带有64Mb RAM的任何基于Pentium的个人电脑上无需修改即可运行。

第一章：介绍

1.1  通用

美国地质调查局模块化流模型（MODFLOW）可能是全球最流行的地下水代码，被用作支持地下水调查的仿真工具。MODFLOW已经通过HydroGeoLogic, Inc.大大增强，以允许用户以更强大和高效的方式处理复杂的现场问题。增强的代码被命名为MODFLOW- SURFACT（即具有新的流动和污染物传输包）。给新流包分配的SURF首字母缩写代表以下内容：

•S-健壮而高效（矩阵解法、非线性迭代与自适应时间步长、输出控制以及轴对称）方案

•U-Unconfined flow 的严格处理

•R-"非积水补给"和"渗漏面边界条件"

•F-将裂缝井表示法提升到可以严格处理井条件。

HydroGeoLogic, Inc.开发的七个新模块包已经添加到原始的MODFLOW代码中（McDonald和 Harbaugh，1988年）。这些包括：

(1)  BCF4-这个包括轴对称分析选项的块中心流动包，采用可变饱和度表述和伪土壤函数，对无限制流动进行严格处理。该包还包括 Version 2的严格水流和气流模拟功能。

(2)  FWL4/FWL5-Fracture-Well Packages（FWL4和FWL5）使用一维（1-D）断裂管元素对井采条件进行严格处理，以模拟井。

(3)  RSF4-"可容纳非塘水充注和渗流边界条件的补给-渗流面包"

(4)  ATO4-自适应时间步长和输出控制包是一个使用积极的时间步进方案进行瞬态解决方案的软件包，具有自动生成和控制时间步长以及更好的输出控制功能。

(5)  PCG4-Preconditioned Conjugate Gradient Package包含了一种新的PCG矩阵求解器，作为 MODFLOW中现有的Strongly Implicit Procedure（SIP）和Slice Successive Over Relaxation（SSOR）求解器的稳健替代方案。

(6)  NRB1-Newton-Raphson线性化包结合回溯策略，提高了对非受限或非饱和流（空气或水）模拟的鲁棒性。

(7)  OBS1-观测节点软件包，用于记录选定节点处的水文曲线。

这些附加的流动包在MODFLOW-SURFACT中得到了记录。以下章节提供了它们的简要描述、公式、验证和应用示例、输入说明以及样本数据文件。本文档的用户应当具有对原始 MODFLOW代码及其文档的合理了解。对于不熟悉MODFLOW的用户，应参考以下参考文献：

McDonald, M.G., and A.W. Harbaugh (1988), Techniques of Water-Resources Investigations of the United States Geological Survey, Chapter A1: A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model, Book 6, Modeling Techniques.

1.2  MODFLOW-SURFACT中新流包的属性

新的流程套餐和增强功能提供以下关键好处，如下所述。

网格块的完全脱水和重新饱和处理-原始的美国地质调查局MODFLOW的BCF包在瞬态模拟中存在网格块干燥/重新湿润的困难。

准确地勾勒和跟踪地下水位位置，考虑到非饱和区域的流动、延迟产流以及垂直流组成部分-原始的美国地质调查局MODFLOW的BCF包忽略了这些特征。

轴对称流动模拟选项-原始的美国地质调查局（USGS）MODFLOW的BCF包忽略了这些特性。

多层井的总流量自动且正确地分配到井节点——原始的美国地质调查局（USGS）MODFLOW软件中的WEL模块要求对来自多个含水层或多个模型网格层的井抽水进行先验分配。

当井筒储存和超抽水井的需求大于井的供水时，原始的USGSMODFLOW中的WEL包会失败。

处理非集水或规定的集水重新充足条件-原始的美国地质调查局MODFLOW软件的RCH模块忽略了这些特征。

处理渗流边界条件-美国地质调查局原始的MODFLOW RCH模块忽略了这些特征。

A自动生成和控制时间步长的自适应时间步长方案——原始的美国地质调查局（USGS）MODFLOW的OC包忽略了这些功能。

B更好地组织和控制模拟输出-原始的USGS MODFLOW的OC程序忽略了这些功能。强大高效的PCG矩阵解决方案选项。

对非饱和状态下的水分运移进行建模的能力。

使用双峰水力导度函数实现非平衡优先不饱和流的能力。模拟非饱和空气流动的能力。

稳健高效的牛顿-拉弗逊线性化选项。

可变各向异性选项-此选项允许按单元格基础包括水平各向异性。请注意，MODFLOW要求各向异性在每个模型层内保持均匀。

垂直导水性输入选项-使用此选项，输入网格块的垂直导水性，并在内部计算层间饱和渗透系数。请注意，MODFLOW要求输入层间饱和渗透系数。

美国地质调查局（USGS）为MODFLOW的后续版本（包括MODFLOW-96和MODFLOW-2000）开发的其他模块也已纳入此版本的MODHMS中。从MODFLOW-96中纳入的模块包括瞬态渗漏包（TLK）、直接求解器包（DE4）、河流路由包（STR1）、水力流动障碍包（HFB）、层间储存包（IBS）、恒定（时变）水头包（CHD）、一般有限差分包（GFD）以及流量和水头边界（时变）包（FHB）。

从MODFLOW 2000中纳入的模块包括水库包（RES1）、河流流量路由包（SFR1）、湖泊包（LAK3）、测量点包（GAGE）和蒸发蒸腾片段包（ETS1）。请注意，为了使LAK3包与本文中的BCF包正确配合工作，LAYCON应设置为非围限选项（LAYCON=3, 13, 23, 33, 43）。如果未满足此条件，则输出清单文件中将反映此警告。

•    用于将MODPATH与FWL4、FWL5和ATO4模块配合使用的修改。此改进将时间步信息提取到二进制*.mph文件中（当使用ATO4模块时），以便与MODPATH的修改版本一起使用。原始版本的MODPATH仅利用OC模块的时间步信息，并不识别MODHMS中FWL4和FWL5模块的单元格输出：请注意，由于MODPATH不包括表面水模块，因此只能将其用于MODHMS的地下部分。

•曲线网格选项 - 此选项允许在平面上使用非矩形网格。因此，对于每个层中的每个节点，需要输入DELR和DELC，并且在第一层下方垂直堆叠多个层。请注意，原始的MODFLOW网格在平面上是矩形的。

•LAK2模块，由GeoTrans的Greg Council在1996年开发，已获得作者许可并已实施。有关该模块的输入说明和实施细节可从GeoTrans获取。

•已包含观测节点模块，可列出在任何规定的观测位置随时间的水头或污染物浓度的突破情况。

MODFLOW-SURFACT的所有模块均完全兼容，并且可以在单个模拟中混合使用来自原始MODFLOW和MODFLOW-SURFACT的模块。MODFLOW的模块化结构得到了保留，并且增强功能的模块从主程序中调用。新的公式和计算方案已经通过解析解和其他数值模型进行验证，并且在附加包中实施。MODFLOW-SURFACT的附加包与原始MODFLOW代码之间保持完全兼容性。MODFLOW-SURFACT的输入准备简单明了，并遵循MODFLOW格式结构。流动建模的输入数据基本上是原始MODFLOW所需的数据，增强功能需要的额外数据很少。传输模拟利用所有MODFLOW数据集，并且消除了提供的输入信息的重复。传输参数通过附加数据集提供给代码，并且边界条件通过MODFLOW边界数据集实施。

1.3  MODFLOW-SURFACT操作和输入选项

MODFLOW-SURFACT可以以独立模式运行，也可以作为一个包含图形用户界面的完整套件的一部分。如果用户仅获取MODFLOW-SURFACT本身，或希望在用户界面系统之外执行模拟，则第一种选项是相关的。

本节重点介绍第一种选项（即以独立模式运行MODFLOW-SURFACT）。在运行代码之前，用户需要为给定的模拟问题创建输入数据文件。为了构建地下水流模拟的输入数据文件，应采取两个关键步骤：

1.   R参考原始的MODFLOW用户手册（McDonald和Harbaugh，1988年），并按照手册中的说明为相关的流动模块创建输入文件。

2.如果希望使用新的SURF模块，则继续阅读本文档的相关章节，并根据所选模块的附加输入说明编辑相应的数据文件。

如果用户已经有现有的MODFLOW输入文件，则只需进行第2步。

1.4  文档组织和使用指南

该文档分为九个章节。以下是其余章节及其目的的概述。

第2章介绍了BCF4模块，其中包含新的可变饱和和轴对称流动公式。该章节包括公式描述、验证和应用示例，以及BCF4模块的输入说明。由于BCF4替代了以前的BCF模块，第2章为用户提供了在IUNIT(1)中设置块中心流输入文件的最一般的指导。

第3章介绍了FWL4和FWL5模块，这些模块包含了处理多层井和特殊井条件（例如，井筒储存和超泵水）的新方案。对于希望严格处理多层井并考虑特殊井条件的用户，应查阅本章。

第4章介绍了RSF4模块，该模块包含了处理充水边界条件（带有规定或零蓄水高程） 和渗流边界条件的新方案。对于处理具有复杂边界条件的非围限流问题的用户，应查阅本章。

第5章介绍了ATO4模块，该模块包含自适应时间步长方案和时间步长尺寸的自动控制，以及模拟输出控制。这些方案使得模拟能够更有效地进行。因此，对于对此类选项感兴趣的用户，应阅读本章。

第6章介绍了PCG4模块，该模块包含预处理共轭梯度（PCG）和Orthomin矩阵解法方案。对于有兴趣尝试PCG求解器的用户，应查阅本章。

第7章介绍了NRB1模块，用于牛顿-拉夫逊线性化。请注意，PCG4模块需要与NRB1模块一起使用，因为它是MODFLOW中唯一能处理不对称矩阵的求解器。

第8章介绍了OBS1模块，用于记录指定观测节点的时间变化水头。

第9章提供了前面章节引用的参考文献列表。文档还包括五个附录（A至E），其中包含了描述第2至5章中所述示例问题的附加信息和输入数据文件。用户可以选择使用这些文件进行测试运行，并深入了解MODFLOW-SURFACT提供的新功能。另外四个附录（F至I）提供了MODFLOW-SURFACT/MODHMS扩展模块的输入说明。
 

第二章：块中心流（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.2  A 利用可变饱和流动公式来建立非围限流动方程

水在可变饱和系统中的三维运动可以表示为（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.2  B 变饱和水流方程的建立

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

和（van Genuchten，1977）：

式中，n、γ为经验参数，Se为有效含水饱和度，定义为Se=(Sw-Swr)/(1-Swr)，其中Swr为残余水饱和度。当选择参数n使得n=1+2/γ时，Brooks-Corey表达式会生成与vanGenuchten函数类似的曲线（vanGenuchten，1980）。

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

其中S*是临界点或断点饱和值，其中流量从毛细管发生变化

主导到非毛细管主导，反之亦然，kr*是对应于S*的相对渗透率，δ是代表有效大孔隙度或其他有助于非毛细管主导流动的结构特征的拟合参数。还提供了表格输入选项，可以适应其他功能形式。

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

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

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

2.3  A 非承压流动验证示例

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

所考虑的问题涉及以3658ft3/d的速率从均质且各向同性的无承压含水层中抽水的完全穿透井。含水层的初始饱和厚度为40ft，井在含水层厚度较低的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英尺x2,000英尺的正方形，井位于中心。为了利用对称性，我们仅离散流域的一个象限并考虑井流量Q的四分之一。

3D网格由20层组成，每层厚2英尺，有21行和21列。水平网格间距（Δx和Δy）从0.5英尺到100英尺不等。VCONT（垂直导水率除以两个相邻节点层之间的层间距离）的值计算为10.5d-1。每层的恒定主存储系数（SF1，定义为特定存储时间块厚度）被指定为2.4  x10–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定义为：

sD = 4BKbs/Q

and

ty = Kbt/r2Sy

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

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

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

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

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

Huang等人(1996)详细介绍了这个问题，涉及非饱和垂直土柱中的一维渗透，160cm长，离散成80层2cm厚。整个域的初始水头为-104cm，代表极度干燥的条件，并且在t=0时在顶部应用50cm/d的恒定渗透通量。对典型的粗质地土壤进行了长达0.8天的分析，其水力特性如下：

Saturated hydraulic conductivity, K=100 cm/d

Specific yield (i.e., porosity), Sy= 0.45

Specific storage, Ss=1×10-5

van Genuchten α=0. 2cm-1

van Genuchten β=5

Residual water saturation, Swr=0.11

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

建模方法和结果

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

图2.1 无量纲缩减与无量纲时间关系显示了数值解和解析解的比较

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

该问题涉及截面尺寸为水平15厘米和垂直10厘米的矩形土板中的二维非饱和流动。使用Δx=1cm和Δz=1cm的恒定节点间距来离散化域。整个域的初始水头为-90厘米。沿着域的左边缘的顶部4cm被规定为6cm的边界水头，并且域的右边缘被维持在-90cm的规定水

头。其余边界（即顶部边界、底部边界和左边缘下6厘米）为无流动条件。土壤的水力特性如下：

Saturated hydraulic conductivity, K=1 cm/d

Specific yield (i.e., porosity),

Sy= 0.45

Specific storage, Ss=1×10-14

van Genuchten α=0.005 cm-1

van Genuchten β=2

Brooks Corey parameter, n=2

Residual water saturation, Swr=0.3333

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

建模方法和结果

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

2.4  应用实例

提供了两个示例问题来演示新BCF4包的可变饱和流选项的应用。在这两个示例中，都允许发生足够的抽水以使无承压含水层的部分去饱和。第一个示例是稳态泵送模拟。第二个例子是瞬态模拟，其中在四年的压力期后关闭油井以恢复含水层。该场景还使用MODFLOW中现有的润湿选项（McDonald等人，1991年）和STAFF3D有限元代码

（HydroGeoLogic，Inc.，1992年）进行建模。下面对不同方法获得的结果进行比较和讨论。

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

该测试问题考虑了图2.2所示的300英尺厚的无承压含水层。建模域是一个尺寸为75,000英尺x75,000英尺的正方形。含水层的顶部和底部海拔分别为50英尺和–250英尺。该区域受到3.0x10–9英尺/秒的均匀、连续的垂直补给。图中左（西）边界为恒水头边界，其余边界为无流边界。含水层底部100英尺处有15口井，每口井的抽水率为0.95立方英尺/天。井位如图2.2所示。含水层参数包括：

Horizontal hydraulic conductivity, Kh (Kxx and Kyy) =10-4 ft/s

Vertical hydraulic conductivity, Kv (Kzz) =10-5ft/s

建模方法

该域被均匀离散为 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.C1 用于气相流模拟的喷射井附近气压的稳态和瞬态分布

图 2.2 无承压含水层系统示意图

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

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

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

Horizontal hydraulic conductivity, Kh (Kxx, Kyy) =10–4 ft/s

Vertical hydraulic conductivity, Kv (Kzz) =10–5 ft/s

Specific Yield, Sy=0.03

Specific Storage, Ss=10–6 ft–1

Modeling Approach

与问题1一样，域被均匀离散为3层、15行和15列。MODFLOW所需的VCONT（垂直导水率除以相邻节点层之间的层间距离）的值计算为10–7s–1。每层指定恒定的主存储系数（SF1定义为特定存储时间块厚度）10–4。顶部和中间离散层的西侧规定恒定水头为零。最初，假设域中的头为零。位于第3层的15口井在前4年内被允许以3.85ft3/s的速率进行抽水。在随后的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.3 第1层中的稳态水头分布-情况(a)和情况(b)之间的比较。

图2.4 第2层中的稳态水头分布-情况(a)和情况(b)之间的比较。

图2.5 第3层中的稳态头分布-情况(a)和情况(b)之间的比较。

图 2.6 用于无承压含水层系统瞬态分析的每口井所施加的抽水应力。

图 2.7 抽水（第一）周期（t=4年）结束时第1层的水头分布显示了先前和新模拟选项之间的比较

 

图2.8a 使用先前的MODFLOW选项获得泵送（第一）周期（t=4年）结束时第 3 层的水头分布。

图2.8b 使用MODFLOW-SURFACT/MODHMS中新的可变饱和流量选项获得泵送（第一）周期（t=4年）结束时第3层的水头分布。

 

图2.9 恢复（第二）期结束时（t=16年）第1层的水头分布显示了先前和新模拟选项之间的比较。

图2.10 恢复（第二）期结束时（t=16年）第3层的水头分布显示了先前和新模拟选项之间的比较。

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

Observation Point #	Layer	Row	Column
1	1	7	11
2	2	11	11
3	3	11	10

三个观察点处水头的时间变化（图2.11）表明，两种情况[情况(1)和(2)]非常一致，直到情况

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

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

对以块为中心的流(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.12 三个观察点处头部的时间变化显示了MODFLOW-MODHMS和STAFF3D中新模拟选项之间的比较。

2.6  BCF4 封装的输入指令

2.6.1     General

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 Input

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

1a. Data:

ISS

IBCFCB HDRY IWDFLG WETFCT IWETIT IHDWET

Format:

I10

I10

F10.0

I10

F10.0

I10

I10

*** 仅当使用表格保留和渗透函数时才输入项目 1b（IREALSL = 5 或 6）

1b.  Data:   N_HSK

Format: I10

2.   Data: Format:

LAYCON (NLAY) (Maximum of 200 layers) 40I2

（如果层数为  40  或更少，则使用一条记录；否则，根据需要使用更多记录。）

**仅当需要空气相流模拟时才输入第 3 项（IREALSL=3 或 4）。

**仅当模拟不是轴对称 (IAXIS # 0) 时才输入第 4 项。

**仅当 IANIXY=0 时才输入项目 4a

4a.  Data:

Utility Module:

TRPY (NLAY)

U1DREL.   请注意，U1DREL 是读取真实一维 (1-D) 数组

（即第 3 项的 TRPY）的实用程序模块。

**Enter item 4b only if IANIXY=1

4b. Data:

Utility Module:

TRPYNU (NCOL, NROW)

U2DREL.   请注意，项目 4b 为模型网格的每一层输入一次。

**Enter item 5a only if ICURVL=0

5a.  Data:

Utility Module:

DELR (NCOL) U1DREL

**仅当模拟不是轴对称 (IAXIS # 0) 时才输入第 6 项。

**仅当 ICURVL=0 时才输入项目 6a

6a. Data:

Utility Module:

DELC (NROW) U1DREL

以下二维数组的子集用于描述模型网格每一层的属性。 每层所需的阵列取决于每层的输入 LAYCON 值（参见第 2 项）、模拟是瞬态（ISS = 0）还是稳态（ISS 不为 0），以及润湿功能是否处于活动状态（IWDFLG 非 0) 对于 BCF2 和 BCF3 的润湿方案。 请注意，对于 BCF4 封装中引入的新的可变饱和流量方案，IWDFLG 应设置为零（即输入 LAYCON=43）。 如果不需要数组，则必须将其省略。 不存在需要所有数组的情况。 首先读取第 1 层所需的数组（第 8-15 项）； 然后是第2层的数组，依此类推，直到所有层的属性数组都输入完毕。

**仅当模拟是瞬态时（ISS = 0），才输入第 8 项。

8.   Data:

Utility Module:

SF1 (NCOL, NROW) U2DREL

**仅当该层的输入 LAYCON 值为以下之一时，才输入项目 9： 0、2、10、12、20、22

9.   Data:

Utility Module:

TRAN (NCOL, NROW) U2DREL

请注意，对于输入 LAYCON=40 或 43，水力传导率、底部和顶部高程分别在第 10、11 和 14 项中输入。

*仅当该层的输入 LAYCON 值为以下之一时，才输入项目 10、11 和 12： 1, 3, 11, 13, 21, 23, 31, 33, 40, 43

10. Data:

Utility Module:

HY (NCOL, NROW) U2DREL

12. Data:

Utility Module:

BOT (NCOL, NROW) U2DREL

**仅当该层不是网格的最底层且 IVHYC = 0 时，才输入第 13 项。

13. Data:

Utility Module:

VCONT (NCOL, NROW) U2DREL

**仅当模拟是瞬态的 (ISS = 0) 并且该层的输入 LAYCON 值是以下之一时，才输入项目 14：2, 3, 12, 13, 22, 23, 33, 43

14. Data:

Utility Module:

SF2 (NCOL, NROW) U2DREL

**仅当该层的输入 LAYCON 值为以下之一时，才输入项目 15：2, 3, 12, 13, 22, 23,

33, 40, 43

15. Data:

Utility Module:

TOP (NCOL, NROW) U2DREL

**仅当润湿能力处于活动状态（IWDFLG 非零）并且该层的输入 LAYCON 值为以下之一时，才输入项目 16： 1, 3, 11, 13, 21, 23, 31, 33

16. Data:

Utility Module:

WETDRY (NCOL, NROW) U2DREL

注意：WETDRY 数组与输入 LAYCON = 40 或 43 无关，因此，新的可变饱和流方案不需要第 14 项。

***仅当 IREALSL=7 且输入 LAYCON 为 43 时才输入项目 21。

21. Data:

Utility Module:

S_STAR (NCOL, NROW) U2DREL

***Enter item 22 only if IREALSL=7 and input LAYCON is 43.

22. Data:

Utility Module:

D_STAR (NCOL, NROW) U2DREL

**仅当需要表格保留输入（IREALSL= 5 或 6）并且输入 LAYCON 为 43 时，才输入项目 23、24 和 25。请注意，项目 23、24 和 25 对每个插值点输入一次，然后对每个插值点重复输入一次。 表格输入的剩余 N_HSK 插值点。

23. Data:

Utility Module:

H_VAL (NCOL, NROW) U2DREL

24. Data:

Utility Module:

25. Data:

Utility Module:

S_VAL (NCOL, NROW) U2DREL

K_VAL (NCOL, NROW) U2DREL

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

ISS – 是稳态标志

如果 ISS = 0，则模拟是瞬态的。

如果 ISS 不为 0，则模拟处于稳态。

IBCFCB      –是一个标志和一个单元号。

If IBCFCB > 0,  它是每当设置 ICBCFL 时将记录逐个单元流量项的单元号。

（参见 McDonald 和 Harbaugh [1988] 的输出控制。）

If IBCFCB = 0,  不会打印或记录逐个单元的流动术语。

If IBCFCB < 0,  每当设置  ICBCFL  时，就会打印每个常压头单元的流量。

HDRY  –   是分配给在模拟过程中转换为干燥的细胞的头部。

尽管该值在模型计算中不起任何作用，但在检查模型输出的模拟水头时，它可作为一个有用的指标。 因此，HDRY 类似于基本 (BAS) 包中的 HNOFLO，它是在模型模拟开始时分配给非流动单元的值。

注意：当层的输入 LAYCON 值为 40 或 43 时，HDRY 值不会分配给模拟期间干燥的单元，如其他      LAYCON    选项 0、1、2、3、10、11、12、13、20 中那样 、21、22、23、31 和 33。

IWDFLG    –   是确定润湿功能是否处于活动状态的标志。

If IWDFLG = 0, 润湿能力无效。

如果 IWDFLG 不为 0，则润湿功能处于活动状态。

注意：如果第 2 项中的任何输入LAYCON 值为 40 或 43，则不使用IWDFLG。

WETFCT  –   是在单元从干转换为湿时最初建立的水头计算中包含的一个因子。

（请参阅 IHDWET。）润湿功能未激活时留空（即 IWDFLG = 0）。

注：该标志与输入 LAYCON = 40 或 43 条件无关，并且对于这些选项可以留空。

IWETIT  – 是尝试润湿细胞的迭代。 在每次 IWETIT 迭代期间都会尝试细胞润湿。 如果 IWETIT 为 0，则更改为 1。润湿功能未激活时留空（即 IWDFLG = 0）。

注：该标志与输入 LAYCON = 40 或 43 条件无关，并且对于这些选项可以留空。

IHDWET  –   是一个标志，用于确定使用哪个方程来计算变湿单元的初始水头。

如果 IHDWET = 0，则 McDonald 等人的方程 (3a)。 （1991）使用：

h = BOT + WETFCT (hn - BOT).

如果 IHDWET 不为 0，则 McDonald 等人的方程

（3b）。 （1991）使用：

h = BOT + WETFCT（阈值）。

当润湿功能处于非活动状态时留空（即 IWDFLG = 0）。 注：该标志与输入 LAYCON = 40 或 43 条件无关，并且对于这些选项可以留空。

注意：对于轴对称仿真方法 (IAXIS >0)，基本封装输入中的 NROW 必须设置为 1。

IREALSL -- 是一个标志，表示实际土壤湿度函数用于定义地下水位以上非饱和区域的流量.

如果IREALSL=0,   拟土壤关系用于定义地下水位。如果 IREALSL > 0，则严格处理非饱和区的流动。

如果IREALSL=1,   van Genuchten 函数用于非饱和层的滞留和相对渗透特性，并模拟水流。

If IREALSL=2, 采用van  Genuchten函数表征保湿性，采用

Brooks-Corey函数表征非饱和层的相对渗透率特征并模拟水流。

If IREALSL=3, 利用van  Genuchten函数计算非饱和层的滞留

特性和相对渗透特性并模拟气流。模拟非饱和层的相对渗透特性和气流。

If IREALSL=4, van  Genuchten函数用于表征保湿性，Brooks-

Corey函数用于非饱和层的相对渗透特性，并模拟气流。

If IREALSL=5, 使用保水性和相对渗透率曲线的表格输入来表征非饱和层，并模拟水流。

If IREALSL=6, 使用保水性和相对渗透性曲线的表格输入来

表征非饱和层，并模拟空气流动。

If IREALSL=7,  Mohanty (1997) 的双峰水力传导率曲线用于表征具有优先路径的非饱和区水流，如第 2.2B 节中所述。

Note:    This flag is required only for input LAYCON=43, 53, 44 or 54.

ICNTRL     --     is a flag determining selection of mid-point or upstream weighting of the relative permeability term.

If ICNTRL = 0,  upstream weighting is used (suggested option). If ICNTRL = 1,    midpoint weighting is used.

Note: This flag is required only for input LAYCON=43, 53, 44 or 54.

IVHYC --     is the vertical hydraulic conductivity input flag.

If IVHYC = 0,    vertical hydraulic conductivities are not read.

Leakances are read for multi-layer simulations.

If IVHYC = 1 ,   vertical hydraulic conductivities are read.

Leakances are computed from vertical conductivity of each layer as per a volume weighted harmonic mean.

IANIXY      --     is the index for variable anisotropy on a cell-by-cell basis.

If IANIXY = 0,  horizontal anisotropy is uniform within each layer as in the original MODFLOW formulation.

If IANIXY=1,    horizontal anisotropy is non-uniform within a layer and is input on a cell-by-cell basis for all layers.

ICURVL     --     is the index for using a curvilinear grid in the areal plane.

If ICURVL =1,   a curvilinear grid is used with input of DELR and DELC required for every node within a layer.

Note that layers are stacked vertically below each other, hence, DELR and DELC for a node in lower layers are the same as for the node vertically above it.

If ICURVL =0,   a rectangular grid is used, as in the original

MODFLOW.

IOLDHDR --     is a flag indicating if old header formats (MODFLOW-style) are to be used for saving cell-by-cell flow or budget terms in the binary file.

If IOLDHDR =1, the old MODFLOW style headers are used for cell- by-cell flow term output in the binary file.

If IOLDHDR =0, a new header format is used that includes the total time value (TOTIM) and the stress-period time value PERTIM in the header for cell-by-cell budget term outputs.

N_HSK --     is the number of interpolation points used for tabulated input of moisture retention curves and relative permeability functions for unsaturated zone layers, when tabular input of these functions is requested (via IREALSL=5 or 6).

LAYCON -- is the layer-type index array.    Each element of LAYCON holds the identification for the respective layer.  Read one value for each layer.

There is a limit of 200 layers. If there are 40 or fewer layers, use one record; otherwise, use two or more records. Leave unused elements in a record blank. LAYCON is the parameter through which the new options of the BCF4 Package may be implemented, as discussed below.

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

0    - Strictly confined

1    - Strictly unconfined

-     - 在整个模拟过程中，层的透射率和存储系数保持恒定。

-     该层的透射率各不相同。 它是根据饱和厚度和导水率计算得出的。 存储系数随时间恒定。 此代码仅对第 1 层有效。

2    - Confined/Unconfined     -      该层的透射率保持恒定。  存储系数可以在受限值

和非受限值之间交替。

3    - Confined/Unconfined     -      该层的透射率允许变化。     它是根据饱和厚度和导

水率计算得出的。  存储系数可以在受限值和非受限值之间交替。

上述描述可以总结如下（Goode and Appel，1992）：第二位数字为

the input LAYCON

value

      Layer Characteristics 

0    Strictly confined; T and S constant in time

1    Strictly unconfined; T depends on head, S constant in time

2    Confined/Unconfined; T constant in time, S depends on head because of storage conversion

3    Confined/Unconfined; T and S depend on head

请注意，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 的存储值：

Input LAYCON Value	Stored LAYCON Value	Stored LAYAVG Value
0	0	0
1	1	0
2	2	0
3	3	0
10	0	10
11	1	10
12	2	10
13	3	10
20	0	20
21	1	20
22	2	20
23	3	20
31	1	30
33	3	30
40	0	40
43	3	40
44	4	40
50	0	50
53	3	50
54	4	50

**请注意，MODHMS 不允许输入 LAYCON 值 30、32、41、42、51、52。 LAYAVG（因子为十的十位数字）的含义如下：

Stored

LAYAVG Value         Meaning      

0    Harmonic mean interblock transmissivity (McDonald and Harbaugh, 1988)

10  Arithmetic mean interblock transmissivity (Goode and Appel, 1992)

20  Logarithmic mean interblock transmissivity (Goode and Appel, 1992)

30 Arithmetic mean saturated thickness times logarithmic mean hydraulic conductivity (Goode and Appel, 1992)

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

RHOWP	–	Density of water.
RHOAP	–	Density of air at reference pressure.
VISW	–	Viscosity of water.
VISG	–	Viscosity of air.
COMPWAT	–	Compressibility of water.
COMPAIR	–	Compressibility of air.
ATMGP	–	Standard atmospheric pressure.
GRAV	–	Gravitational acceleration.
INITW	--	Index for ambient water head input
	=	0 when water is at hydrostatic equilibrium with heads computed from a prescribed water table elevation
	=1	when ambient head in water is input for all nodes of the domain.

RHOWP、RHOAP、VISW、VISG、COMPWAT、COMPAIR、ATMGP 和

GRAV 的值取决于准备输入数据集时使用的单位。 下表列出了这些参数在标准条件下、采用不同常用单位组的值。

Parameter	SI Units	cgs Units	fps Units
RHOWP	1000 kg/m3	1 g/cm3	1.939 slugs/ft3
RHOAP	1.777 kg/m3	1.777×10-3 g/cm3	.003445 slugs/ft3
VISW	1×10-3 N.s/m2	1.0 cP	6.72×10-4 Rb/ft.s (Rbf.s/ft2)
VISG	1.983×10-5 N.s/m2	1.983×10-2 cP	1.3327×10-5 Rb/ft.s (Rbf.s/ft2)
COMPWAT	1-8 Pa-1	10-9 cm2/dyn	4.788×10-7 ft2/Rbf
COMPAIR	1.777×10-5 Pa-1	1.777×10-6 cm2/dyn	8.5083×10-4 ft2/Rbf
ATMGP	1.01×105 Pa	1.01×106 dyn/cm2	2.1094×103 Rbf/ft2
GRAV	9.8066 m/s2	980.66 cm/s2	32.2 ft/s2

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

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

LOCAT –     is an index for the location of values to fill-in the array of concern.

<0, read an unformatted record containing the array values

= 0, set all the array values equal to constant (CNSTNT)

>0, read the formatted records containing the array values.

CNSTNT  –  is a constant to which all array values are set if LOCAT is equal to zero or by which all array values are multiplied if LOCAT is not equal to zero.

DELR  – is the cell width along rows. Read one value for each of the NCOL columns. This is a single array with one value for each column. For an axi- symmetric simulation (IAXIS>0), DELR represents radial width ()r) of a cylindrical cell.

DELC  –     is the cell width along columns. Read one value for each of the NROW rows. This is a single array with one value for each row.

Note:    DELC array is NOT read for an axi-symmetric simulation (IAXIS>0).

SF1       -- is the primary storage coefficient, the storability. Read only for a transient simulation (steady-state flag, ISS, is 0). For input LAYCON equal to 1, 11, 21, or 31, SF1 will always be specific yield, while for input LAYCON equal to 2, 12, 22, 3, 13, 23, 33, 40, or 43, SF1 will always be confined storage coefficient. For input LAYCON equal to 0, 10, or 20, SF1 would normally be confined storage coefficient. The primary storage coefficient divided by the block thickness is equal to the specific storage (Ss) of equations (1) and (2).

Note (1): An input LAYCON value of 0, 10, or 20 can also be used to simulate water-table conditions where drawdowns are expected to remain everywhere a small fraction of the saturated thickness and where there is no layer above, or flow from above is negligible. In this case, specific yield values would be entered for SF1.

TRAN  –  is the transmissivity along grid rows. TRAN is multiplied by TRPY to obtain transmissivity along columns. TRAN is read only for a layer represented by one of the following input LAYCON values: 0, 2, 10, 12,

20, 22.

HY –  is the hydraulic conductivity along grid rows. HY is multiplied by TRPY to obtain hydraulic conductivity along columns. HY is read only for a layer represented by one of the following input LAYCON values: 1, 3, 11, 13, 21,

23, 31, 33, 40, 43.

BOT     -- is the elevation of the aquifer bottom. Read only for a layer represented by one of the following input LAYCON values: 1, 3, 11, 13, 21, 23, 31, 33,

40, 43.

VCONT  -- is the vertical hydraulic conductivity divided by the thickness, from one layer to the layer below. The value for a cell is the hydraulic conductivity divided by thickness for the material between the node in that cell and the node in the cell below. Because there is not a layer beneath the bottom layer, VCONT is not specified for the bottom layer.

SF2       -- is the secondary storage coefficient, the specific yield. SF2 is read only if the simulation is transient (steady-state flag, ISS, is 0) and input LAYCON value for the layer is one of the following: 2, 3, 12, 13, 22, 23, 33, 43.

TOP      -- is the elevation of the aquifer top. TOP is read only for a layer represented by one of the following input LAYCON values: 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     --     is the van Genuchten parameter $ for unsaturated soils. VANBT is read only when the input LAYCON value for the layer is 43, and the real soil moisture functions are requested with IREALSL…0.

 VANSR     --     is the residual saturation level for unsaturated soil. VANSR is read only when the input LAYCON value for the layer is 43, and the real soil moisture functions are requested with IREALSL…0.

BROOK      --     is the Brooks-Corey exponent for computing relative permeability of unsaturated soil. BROOK is read only when the input LAYCON value for the layer is 43, and the real soil moisture functions, with Brooks- Corey permeability relations are requested with IREALSL=2 or 4.

S_STAR      –     is the critical saturation parameter S* (see equation 5b of Section 2.2B) for unsaturated soils with the bimodal conductivity function of Mohanty et al (1997). S_STAR is read only when the input LAYCON value for the layer is 43, and the bimodal conductivity function is requested with IREALSL=7.

D_STAR     –     is the effective macroporosity fitting parameter (see equation 5b of Section 2.2B) for unsaturated soils when using the biomodal conductivity function of Mohanty et al (1997). D_STAR is read only when the input LAYCON value for the layer is 43, and the bimodal conductivity function is requested with IREALSL=7.

H_VAL –     is the pressure head value for an interpolation point of tabular input for moisture retention and relative permeability functions. H_VAL is read only when the input LAYCON value for the layer is 43, and tabular input of functions is requested with IREALSL=5 or 6.

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

K_VAL–是与 H_VAL 相关值的插值点对应的活动流动相的相对渗透率值。当层的输入 LAYCON 值为 43 且 IREALSL=5 或 6 时请求函数的表格输入时，K_VAL 为只读。

第四章：补给-渗流面边界（RSF4）条件包

4.1  概述

MODFLOW代码使用recharge（RCH1）软件包适应地下水补给边界条件。但是，RCH1包在表示未受限制的系统方面有一定的局限性。如果无侧限含水层饱和到地表，含水层吸收补给的能力就会降低，剩余的水会作为地表径流排

出。RCH1软件包无法处理这种情况。相反，它继续向含水层提供补给（如在密闭系统中），水头不断上升，远远高于地表。MODFLOW-SURFACT的新RSF4软件包旨在消除这种物理上不现实的条件。如果地下水位低于用户规定的水池（积水）高程，RSF4包允许向地下水系统补给水。水池高程与无积水时的地面高程相对应。如果在任何时候，地下水位达到水池高程，模拟只允许进行尽可能多的补给，以保持规定的水池条件。剩余的充值不被接受进入模拟域。模拟的输出反映了由于系统饱和到积水高度而导致的系统体积补给通量（L3/T）的减少。请注意，RSF4包与流路由包（STR1或SFR1）断开连接，因此，该剩余充值不会以任何方式路由到规定的流中。

此外，对于MODHMS模拟，无侧限渗流选项与陆上流动域相冲突，因为陆上流动区域处于活动状态，因此，如果RSF4软件包中使用无侧限渗漏选项，则应在这些位置设置非常高的积水高程，以防该软件包干扰陆上流动。

新RSF4软件包的额外功能允许方便地输入瞬态补给/降水，而不依赖于MODFLOW的“应力期”概念。当使用此选项时，可以通过单独的文件输入充值时间序列，该文件确定域上的充值如何随时间变化，而不考虑压力期设置，以在充值快速且经常在模拟期内变化时提供输入便利性和通常较小的输入文件大小。

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

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

4.2  RSF4边界条件包的公式

RSF4软件包可应用于任何涉及补给和大气渗流面边界条件的情况。该包所需的输入数据包括再充电率，该再充电率存储在二维（2-D）阵列（RECH）中，每个水平单元位置有一个元素，如之前的MODFLOW再充电包中所述。此外，为了允许无限制的补给条件，RSF4包需要水池高程

（PNDEL）作为输入。

只要单元的地下水位高程（HNEW）保持在其水池高程以下，整个补给率（RECH）就会应用于边界。如果在非线性迭代过程中，地下水位超过水池水位，则维持水池水位，并将系统的补给降至水池高程条件下系统可接受的最大补给率。如果在模拟过程中的任何时候，地下水位低于规定的水池高程，则对系统应用完全补给率。此外，如果水池高程高于地表，则使用改良的 Picard方法（Celia等人，1990），通过额外的积水量来增加蓄水期。

除非标志IPNDPOR打开，否则假设池塘会发生在单元的整个区域。在这种情况下，网格块内的小湖和池塘的体积（在瞬态情况下）可能会被视为网格块大小的一小部分，并在阵列 PNDPOR中额外读取。再充电模块的逐小区流文件还包括用于瞬态模拟的ponded存储项的预算值。

还应注意的是，逐单元流量文件或质量预算项中列出的“补给”负值表示在渗流面条件下从“补给”边界节点发生的流量。

RSF4软件包还模拟侧向渗流面边界条件。对于划定为可能的渗流面节点的域内的任何节点，该代码提供到（或从）系统的零流量，直到地下水位高程高于渗流面高程，在渗流面高程条件下，水从系统排出。

4.3  快速、瞬态充电条件

在MODFLOW的应力期内，补给和降水被认为是恒定的。对于快速降雨事件，在模拟过程中将短时间段分解为几个应力段会变得乏味。因此，提供了一种选择，在模拟中改变补给或 降水，与应力期无关。瞬态再充电是分区域施加的，RSF4包读取一个索引MXZRCH，指示施加瞬态再充电的最大区域数量。如果MXZRCH为零，则模拟不会调用瞬态充电包。然而，如果存在瞬态再充电，RSF4包会读取一个额外的整数阵列，识别每个区域的区域位置。分区值为零表示该位置上没有出现瞬态再充电区，并且默认的再充电值应用于该位置，如RECH阵列中应力期的读数。在分区值不为零的位置，RECH阵列中应力期的默认再充电将添加到该分区的适当值 中，该值来自单独文件中提供的再充电时间序列数据。此文件，具有默认扩展名。RTS，包含所有区域的充电时间序列表。文件的每一行都包含事件的开始时间（TSTART）、事件的结束时间

（TEND）、乘数（用于单位转换），然后是每个区域的一列补给/降水率。如果后续TSTART值与前一行的TEND值不匹配，则在过渡期内对每个区域进行零充电。因此，不需要将长的中断期作为单独的一行包含在中。RTS文件。

事件之间或事件内的时间段不需要固定，因为此模块设计用于自适应时间步进方案。自

适应时间行进方案确保每个TSTART和TEND都能被模拟准确命中，以确保事件得到准确表示。在一个事件内或事件之间，自适应时间步进方案试图最大化计算效率。

应当注意，RECH阵列中提供的再充电被添加到RTS文件中提供的值，以在任何节点处提供总的再充电。因此，可以通过RSF包的RECH阵列在应力周期的基础上提供施加的灌溉或其他人类补给源，同时可以通过RTS文件提供快速变化的降水事件。

4.4  验证和应用实例

提供了两个示例问题来演示 RSF4 包的验证和应用。 第一个问题考虑的是一个无承压含水层，它从地表补给中接收水，然后排入长的平行排水沟。 第二个问题涉及沿垂直渗水面排水的路堤。 将这些问题与适用于各自情况的解析和数值解决方案进行比较。

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

该问题涉及各向同性、均质无承压含水层中流向平行排水沟的水流，该含水层受到降水的垂直补给，如图 4.1 所示。 由于问题的二维性质，只能分析一个垂直切片。 排水沟之间的距离为 1,000 英尺，含水层深度为 75 英尺。含水层的初始水头为 50 英尺。该区域受到 0.01 英尺/天的恒定降水补给。 排水沟保持 50 英尺的恒定水头。含水层参数包括：

Hydraulic conductivity (K)= 1 ft/d

Specific yield, Sy=0.2

建模方法

如前所述，问题是二维的，并且仅选择单位宽度 (DELC = 1 ft) 的一个切片 (NROW = 1)进行分析。 此外，排水沟之间存在一条对称线（x = 500 英尺），因此只有域的前半部分被离散化，如图 4.2 所示。 该域被离散为 10 层（每层厚度 7.5 英尺）和 20 列（每列长度 25 英尺）。 含水层的底部位于基准面 (0 英尺)，顶部位于 75 英尺。每层的 VCONT 值（垂直电导率除以两个相邻层之间的层间距离）为 0.1333 d–1。 排水节点处的初始水头值和边界水头位于 50 英尺处。所有层均使用输入 LAYCON 值 43，以利用新的可变饱和流公式，并选择 SIP 包来求解方程 。 模拟了两种情况：案例 1 向域应用了 0.01 ft/d 的补给，案例 2 向域应用了 0.02 ft/d的补给。

结果与讨论

案例 1 和案例 2 不同时间值的模拟结果分别如图 4.3 和 4.4 所示。 对于案例 1 模拟，地下水位保持在地表以下，并且整个补给量在整个模拟过程中都应用于该域。 这种情况下的稳态解与图 4.3 中 Bear (1979) 的解析解进行了比较，显示出极好的匹配。 对于情况 2，地下水位在流域的某些部分到达地表。 其中，地下水位保持在水池高程（地表），地下水系统的补给量减少，如模拟输出所示。

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

该问题涉及方形路堤的稳定渗漏。 选择它是为了测试 RSF4 软件包模拟渗流面的能力。 Cooley (1983) 提供了该问题的严格（可变饱和流）数值解。 几何形状、边界条件和各向同性土壤参数如图 4.5 所示。 一个 10 m x 10 m 的各向同性路堤，其各向同性电导率特性为 K = 0.0l m/d，左侧可容纳 10  m  的水，右侧可保持 2  m  的水头。 然而，渗流面可能沿着该边界形成。由于该问题的二维性质，仅分析该问题的垂直横截面。

建模方法

问题配置是二维的。 因此，仅选择单位宽度（DELC=1 m）的一个切片（NROW=1）进行分析。在垂直方向上，域被离散为 12 列和 10 层，厚度为 1 m。 在左右边界处，使用0.1 m的柱宽

（DELR）来准确捕捉边界效应。 沿左侧边界应用 10 m 的水头，并在右侧边界的底部两个单元格应用 2 m 的水头。 沿右侧边界的其余单元可能是渗流面，并通过 RSF4 软件包纳入模拟中。 MODFLOW 所需的 VCONT 值为 0.01 d–1。 每层使用输入 LAYCON 值 43，以使用 BCF4包的可变饱和流动公式，并使用 SIP 包求解矩阵方程。 用于模拟的输入数据文件列在附录 D中。

结果与讨论

从初始水头 10 m 开始对问题进行稳态模拟。 图 4.6 显示了沿域底部的模拟地下水位高程和测压头。 这些结果与使用 Cooley (1983) 和 Huyakorn 等人 (1986) 的有限差分和有限元模型获得的结果非常吻合。

图  4.1 示意图描绘了无承压含水层中平行排水沟的流动。

图 4.2 用于模拟平行排水管流动的离散化方案

图 4.3 并行排水问题的数值解和解析解的比较（补给= 0.01 英尺/天）

图 4.4 平行排水问题的地下水位瞬态移动（补给 = 0.02 英尺/天）。

图 4.5 方形路堤的稳定渗流

4.5  RSF4 包的输入指令

4.5.1     RSF4 输入

Input to the Recharge Seepage Face (RSF4) Package is read from the unit specified

in IUNIT(8). The additional input variables required by RSF4 are shown in a shaded box. Note that the new RSF4 Package encompasses the Recharge (RCH1) Package of the original MODFLOW (McDonald and Harbaugh, 1988).

Input for transient recharge conditions occurs via two files. The first is the recharge-seepage face package input file on IUNIT (8), which indicates the presence of transient recharge to the model and number of zones of recharge. The second file contains the recharge time series (default extension .RTS) for each of the zones. This file is format free and opened internally when transient recharge is applied.

1.   Data:      NRCHOP    IRCHCB      IUNCNF      MXSEEP     MXZRCH    IPNDPOR

      Format:  I10  I10  I10  I10  I10  I10

**Enter items 2 through 6 for each stress period.

2.   Data: Format:

INRECH I10

INIRCH I10

**Enter item 3a only if INRECH $ 0.

3a. Data:

Utility Module:

RECH(NCOL, NROW)

U2DREL. Note that U2DREL is the utility module that reads two-dimensional array (i.e., RECH for item 3a).

**Enter item 3b only if INRECH > 0 and MXZRCH > 0. The tabulated, zonal recharge values are added to the default recharge values read above in item 3a.

3b. Data:

Utility Module:

*Enter item 4 only if the recharge option (NRCHOP) is equal to 2 and INIRCH $ 0.

4.   Data:

Utility Module:

IRCH(NCOL, NROW) U2DINT

**Enter item 5 (both 5a and 5b) only if the unconfined option is used (IUNCNF> 0) and INPNDEL $ 0.

5a. Data:      PNDEL(NCOL, NROW)

Utility Module:  U2DREL

**Enter item 5b only if ponding porosity flag is on (IPNDPOR=1)

5b. Data:

Utility Module:

PNDOR(NCOL, NROW) U2DREL

**Enter item 6 only if INSEEP > 0.

6.   Data: Format:

(Input item 6 normally consists of one record for each seepage face boundary cell. If INSEEP is negative or zero, item 6 is not read.)

Explanation of Fields Used in Input Instructions

NRCHOP -- is the recharge option code (flag). Recharge rates are defined in a two- dimensional array, RECH, with one value for each vertical column. Accordingly, recharge is applied to one cell in each vertical column, and the option code determines which cell in the column is selected for recharge.

1    -  Recharge is only to the top grid layer.

2    -  Vertical distribution of recharge is specified in array IRCH.

3    - Recharge is applied to the highest active cell in each vertical column. A constant-head node intercepts recharge and prevents deeper infiltration.

IRCHCB   -- is a flag and a unit number for flow simulation.

If IRCHCB >     0, it is the unit number on which cell-by-cell flow

terms will be recorded whenever ICBCFL (see Output Control) is set.

If IRCHCB #     0, cell-by-cell flow terms will not be printed or

recorded.

IUNCNF     -- is the unconfined option code (flag).

If IUNCNF > 0, recharge seepage face boundary condition is employed and recharge (RECH) will be applied until water table reaches the ponding elevations described in PNDEL array.

If IUNCNF # 0,  the specified recharge (RECH) is applied always like in a confined case.

MXSEEP    -- is the maximum number of seepage boundary cells active at one time.

MXZRCH  -- is a zonal time-series flag for input of recharge

If MXZRCH = n, recharge varies as a time series, and n is the total number of recharge zones used in the simulation.

If MXZRCH = 0, MODFLOW’s input structure for recharge is used with variations only at stress periods.

Note 1:

Note 2:

Note 3:

MXZRCH > 0 indicates that recharge varies at a different (possibly much smaller) time scale than the stress period of MODFLOW, and that an independent time series for recharge will be provided (in file *.RTS) for each zone as detailed in Section 4.5.2.

A zonal index array identifies the grid blocks which lie within a particular zone. An index value of zero allows for the default application of MODFLOW’s recharge for that stress period (as read into the RECH array).

Zonal locations for the recharges in the time series may be varied between stress periods. Thus, the zones themselves may be varied on a stress period basis. This is controlled by the flag INRECH, which also controls reading of recharge rates of each stress period.

IPNDOR   -- is a flag for reading the ponding porosity

If IPNDPOR = 0,  ponding porosity is not read and the full cell area is

assumed to be ponded.

If IPNDPOR = 1,  an array of ponding porosity is read. The ponding

porosity is the fraction of the total cell area that is wetted by lakes/ponds, and is used to determine the total ponded storage, if ponding occurs.

INRECH   -- is the RECH read flag.

If INRECH $ 0,  an array of recharge rates, (RECH) is read.

If INRECH < 0, recharge rates from the preceding stress period are used.

Note:

If MXZRCH > 0, INRECH also controls input array for zonal indexes (IZNRCH). Hence, zonal locations may be changed at every stress period.

INIRCH      -- is the IRCH read flag. When NRCHOP is two,

If INIRCH $ 0,   an array of layer numbers (IRCH) is read.

If INIRCH < 0,  the array (IRCH) used in the preceding stress period is reused.

Note: When NRCHOP is one or three, INIRCH is ignored.

INPNDEL   -- is the PNDEL read flag, when the unconfined option is used (IUNCNF

> 0).

If INPNDEL $ 0,       an array of ponding elevations (PNDEL) is read.

If INPNDEL < 0,      ponding elevations from the preceding stress period

are used.

INSEEP      -- is a flag and a counter of seepage boundary cells.

If INSEEP < 0,   seepage face boundary data from the preceding stress period will be reused.

If INSEEP $ 0,   INSEEP is the number of seepage face boundaries during the current stress period.

RECH    -- is the recharge rate (L/T). Read only if INRECH is greater than or equal to zero.

IZNRCH   -- is an integrater array identifying the zone value of each node.

Note: The default value of recharge as read into the array RECH at each stress period will be applied to the cells with a zone value of zero.

IRCH   -- is the layer number array that defines the layer in each vertical column where recharge is applied. Read only if NRCHOP is two and if INIRCH is greater than or equal to zero.

PNDEL  -- is the ponding elevation array that defines the ponding elevation in each vertical column. Typically, the ponding elevation coincides with ground- surface elevation if there is no water pooling above the ground surface. Read only if IUNCNF is greater than zero and INPNDEL is greater than or equal to zero.

PNDPOR -- is the ponding porosity array that defines the fraction of cell area that is occupied by ponds/lakes, when water level is above the top of the domain, up to the ponding elevation PNDEL. This fraction is used in transient simulations to track volumes/rates of ponded water.

Layer Row Column

-- is the layer number of the cell representing seepage face boundary.

-- is the row number of the cell representing seepage face boundary.

-- is the column number of the cell representing seepage face boundary.

Elevation     -- is the elevation of the seepage face.

4.5.2     Transient Recharge Time-Series File

The transient recharge time-series file is opened internally on FORTRAN Unit 96, when transient recharge conditions are applied. The file is in free format and contains one row for each event within the simulation. Data entries on each row are as follows:

Data:     TSTARTR    TENDR TMULTR     RZ(1)  RZ(2).    RZ(MXZRCH)

Format: Free             

Explanation of fields used in input.

TSTARTR - Starting time of a rainfall/precipitation event TENDR    - Ending time of a rainfall/precipitation event

Note that TSTARTR and TENDR are referenced to the total simulation time and not to the stress period time, when multiple stress periods occur in a simulation.

TMULTR    - Multiplying factor applied to recharge rates of all the zones for the given event. This scaling is convenient for unit conversions.

RZ(1)

RZ(2)

•

•

•

-     Recharge rate in Zone 1 for this event.

-     Recharge rate in Zone 2 for this event.

RZ(MXZRCH) - Recharge rate in Zone MXZRCH for this event.

A sample RTS file is shown in Table 4.1. In this example, the domain is divided into 8 zones of recharge. No rainfall or precipitation occurs for the first 2.375 days of simulation, so the simulation would proceed rapidly till the first recharge event, adapting its time-step size depending on convergence behavior. At 2.375 days, a rainfall event occurs over some of the zones, for a one hour period, till 2.417 days. At 2.417 days, the intensity (and/or zone) of the recharge event changes, with these entries occurring in the second row of the file, and so on. Note that no rainfall event occurs between 2.5 and 4 days, and an entry is not required in the file for this situation.

Table 4.1 Example input file for recharge time series

TSTARTR(d)	TENDR(d)	TMULTR	RZ(1)ft/d	RZ (2) ft/d	RZ (3)	RZ (4)	RZ (5)	RZ (6)	RZ (7)	RZ (8)
2.375	2.417	1.0	0	0	0.225	0.177	.033	0	0	0
2.417	2.458	1.0	0	0	0	0.13	0.13	0.022	0	0
2.458	2.500	1.0	0	0	0	0	0.088	0	0	0
4.000	4.042	1.0	0.169	0.551	0.551	0.337	0.178	0	0	0
4.042	4.083	1.0	0.151	0.447	0.333	0.106	0.08	0	0.691	0.633
4.083	4.125	1.0	0.100	0.050	0	0	0	0.150	0.379	0.300
4.125	4.167	1.0	0	0	0	0			0.083	0.022

Note: MXZRCH = maximum number of recharge zones = 8

第 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  包中使用的自适应时间步长的制定

ATO4包中实现的自适应时间步长选择方案简单而有效。用户提供初始时间步长(DELT)、最大时间步长(TMAX)、最小时间步长(TMIN)、时间步长放大因子(TSMULT)和时间步长缩小因子(TSDIV)对于长度为PERLEN的每个应力周期。初始时间步长DELT用于应力周期的第一个时间步长。如果在最大迭代次数(MXITER)的35%以内实现收敛，则时间步长将增加乘数TSMULT。对于MXITER35%到65%之间的收敛率，时间步长与之前的DELT相比没有变化，并且系统的收敛行为被认为对于该DELT值是最佳的。然而，如果实现收敛所需的迭代次数超过MXITER的65%，则时间步长将减小除法因子TSDIV。如果特定时间步长未实现收敛，则时间步长将减小5.0倍，并对较小的时间步长重新进行计算。另外，如果达到打印时间(即，需要输出打印/保存的目标时间值)，则选择时间步长以包括打印时间值。如果DELT的值减小到小于TMIN，则模拟中止。如果遇到严重困难，该技术提供了退出程序模拟的方法。最后，DELT的值被限制为小于或等于TMAX。即使非线性轻微或不显着，此约束对于时间离散化的准确性也很重要。用户可以灵活地使用TMVEC阵列在指定时间进行打印，或在指定数量的时间步长(NPSTP)后进行打印，或两者兼而有之。

5.3  验证示例

对无承压含水层中的抽水进行瞬态模拟，以验证ATO4包并演示其实施。

BCF4包文章中描述了所选问题，其中使用原始MODFLOW的时间步进和输出控制来执行模拟。在计算效率和鲁棒性方面对这两次模拟运行进行了比较。

本例中的瞬态分析考虑了300英尺厚的无承压含水层，如图5.1所示。建模域对应于尺寸为75,000英尺x75,000英尺的正方形。含水层的顶部和底部分别位于50英尺和-250英尺的海拔处。整个区域含水层的初始地下水位为零。该域受到3.0x10–9ft/s的均匀连续补给。图中的左（西）边界是恒定（零）水头边界，其余边界是无流边界。在含水层底部100英尺范围内筛选了15口井，其位置如图5.1所示。每口井的抽水速度高达3.85ft3/s。四年后，井中的抽水就会停止。（见图5.2。）然后允许含水层恢复12年。含水层参数为：

Horizontal hydraulic conductivity, K xx. Kyy = 10–4 ft/s Vertical hydraulic conductivity, K zz  = 10–5 ft/s Specific Yield, S y      = 0.03

Specific Storage, S s =  10–6 ft–1

图 5.1 无承压含水层系统示意图。

图 5.2 用于无承压含水层系统瞬态分析的每口井所施加的抽水应力

Modeling Approach

The domain is discretized into 3 layers, each 100 ft thick with 15 rows and 15 columns of a uniform block of 5,000 ft areally. Thus, the value of VCONT required by MODFLOW (vertical hydraulic conductivity divided by inter-layer distance between any two layers) is calculated as 10–7 s–1. A constant primary storage coefficient (specific storage times block thickness) of 10–4 is specified for each layer. A constant head of zero is prescribed on the west side of top and middle discretized layers. Initially, the head in the domain is assumed to be zero. Fifteen wells located in Layer 3 are allowed to pump at a rate of 3.85 ft 3/s over the first 4 years. Pumping is stopped during the subsequent period of 12 years, to allow aquifer recuperation.

Two separate transient simulation runs are performed, using (1) the existing MODFLOW time-stepping and output control, and (2) using the new automatic time- stepping and output control (ATO4) package. These two simulations are referred as Case

(1) and Case (2), respectively. The unconfined flow option of input LAYCON = 43 was selected for both cases. The SIP Package is used to solve the matrix equations for both cases. Time-step control parameters for Case (1) include: MXITER = 50, DELT = 7.875 x 10 6 s, and TSMULT = 1. Time-step control parameters for Case (2) include:

MXITER = 50, DELT = 7.875 x 10 6 s, TSMULT = 1.7, TSDIV = 2.0, TMAX = 2

x 10 7 s, and TMIN = 8 x 10 5 s. The input data files using the ATO4 Package [Case (2)] are listed in Appendix E.

Results and Discussions

This problem has been examined previously, in the context of the BCF4 package. This section focuses on comparison of the old and new time-stepping schemes, as applied to the example problem. Figures 5.3 and 5.4 compare head distributions in layers 1 and 3 respectively, at the end of the pumping period (t=4 yr) for both cases, and Figures 5.5 and 5.6 compare head distributions in layers 1 and 3 respectively, at the end of the simulation period (t=16 yr). Comparison of the results is excellent, and the ATO4 package performance is verified.

The number of nonlinear iterations required to complete the simulation were 544 iterations for Case 1, and 412 iterations for Case 2. The CPU times for the Case 1 and Case 2 simulations were 33 seconds and 19 seconds, respectively. The ATO4 package automatically selects the optimal time-step size for Case 2, thereby making the simulation more efficient than for Case 1. Note that the simulation for Case 1 was aborted (due to failure of convergence) if a TSMULT value of 1.7 was used to increment the time-step size, as was done in Case 2. Further, if the maximum iteration count were more restrictive (say MXITER=25 instead of 50), Case 1 simulations would abort, while Case 2 simulations would cut the time-step size for unsuccessful iterations, and proceed until the simulation is completed.

Finally, note that the ATO4 Package used in Case 2 allows for flexibility in output, and the user may demand output at 3.2 years and 4 years of simulation as is noted in the output file. The simulation will perform optimally (using large time-step values possibly exceeding 0.8 years) until it comes to 3.2 years, upon which the simulation solves the system at both 3.2 and 4 years. The previous MODFLOW versions do not allow this, and predetermining the time-step at which one wishes to have output, can be tedious.

图 5.3 抽水（第一）周期（t=4 年）结束时第 1 层的水头分布显示了先前和新的时间步进方案之间的比较。

图 5.4 （第一个）应力期（t=4 年）结束时第 3 层的水头分布显示了先前和新时间步进方案之间的比较。

图5.5恢复（第二）期结束时（t=16年）第1层的水头分布显示了之前和新的时间步进方案之间的比较。

                                                                                          

图5.6恢复（第二）期结束时（t=16年）第3层的水头分布显示了先前和新的时间步进方案之间的比较。

5.4 INPUT INSTRUCTIONS FOR THE ATO4 PACKAGE

5.4.1 General

Input to the ATO4 Output Control is read from the unit specified in IUNIT (22). The previous output option, which is specified in IUNIT (12), is ignored if the new ATO4 Package is activated. If both IUNIT (12) and IUNIT (22) are zero, no output control data are read, and default output control is used. Under the default, head and total budget are printed at the end of every stress period (McDonald and Harbaugh, 1988). Similar to the old output control package, all printer output goes to unit 6 as specified in the main program. If necessary, the unit number for printer output can be changed to meet the requirements of a particular computer.

The ATO package further includes flags and parameters to restart (or continue) a simulation from a non-zero time value. This capability is active when the variable TCNTNU is non-zero and represents the time value (in simulation units of time) from which the simulation is restart. The first time-step size value of the restart simulation is also then input in variable DTCNTNU (in simulation units of time). When this is the case, the simulation will wind itself up to time TCNTNU (appropriately winding stress period information on boundary conditions and time values of the recharge or ET time-series files) and proceed from that point onwards with the first time-step size being DTCNTNU. Initial heads (and concentrations for transport) as well as values of initial storage in interception in the IPT package and initial values of the ON/OFF switch of regulated structures in the CHF package (for MODHMS) are obtained from the respective initial conditions as input for the simulation. Further flags are available for direct input of starting heads from the binary HDS file (and/or starting concentrations from the CON file) of a previous MODHMS simulation for the same problem. When the flag (and unit number) ICNTNU is set to a positive number, the file containing starting heads (all binary head outputs of a previous run) is opened on the specified unit number (and with default extension of HDC—which may be obtained by renaming the previous simulation’s HDS file to extension HDC) and the file is wound up till the specified time (TCNTNU) is reached. The head values (in all simulated domains) available in the file at this time value (or greater, if data is not available for time TCNTNU) are used for starting the current restart simulation. Similarly, when the flag and unit number, ICNTNUC is set to a positive value, and when a transport simulation is performed, the file containing starting concentrations is opened on the specified unit number (and containing default extension of COC—which may be obtained by renaming the previous simulation’s CON file to extension COC). This file is then wound up to the time value of TCNTNU (or greater) to obtain initial concentrations for all simulation domains (which include immobile and mobile domains in the subsurface for dual porosity cases, and, the OLF surface and the CHF domain as appropriate) and all species of the restart simulation.

5.4.2 ATO4 Input

1.                                                                                        Data:    IHEDFM      IDDNFM IHEDUN IDDNUN ICONFM                                                                   ICONUN Format:            I10         I10  I10  I10                                                                                                   I10        I10

TCNTNU                                                                           DTCNTNU ICNTNU      ICNTNUC                                                                                                   IUNITMB F10.3              F10.3            I10                                                                                                     I10        I10

**Enter items 2 through 5 for each stress period.

2.                                                                                        Data:    INCODE IHDDFL IBUDFL ICBCFL ICONFL ICONBD ICCCFL Format: I10                                  I10 I10  I10  I10  I10  I10

3.                                                                                        Data:

Format:                                                                               DELT TMIN TMAX TSMULT

F10.0  F10.0  F10.0                                                        F10.0    TSDIV

F10.0                                                                                  NPRTS

I10                                                                                      NPSTP

I10

4.                                                                                        **Enter item 4 only if NPRTS > 0. Data:                                                                                                   TMVEC(NPRTS)

Format:                                                                               8F10.0 (number of records = NPRTS/8 + 0 or 1)

**Enter item 5 only if INCODE $ 0.

5.                                                                                        Data: Format:     Hdpr I10       Ddpr Hdsv I10                                                                                                   I10 Ddsv I10      Conpr Consv I10 I10

(Record 5 is read 0, 1, or NLAY times, depending on the value of INCODE.

5.4.3 Explanation of Fields Used in Input Instructions

IHEDFM  –  is a code for the format in which heads will be printed.

IDDNFM –  is a code for the format in which drawdowns will be printed. Format codes have the same meaning for both head and drawdown. A positive format code indicates that each row of data is printed completely before starting the next row. This means that when there are more columns in a row than will fit on one line, additional lines are used as required to complete the row. This format is called the wrap format. A negative format code indicates that the printout is broken into strips where only that number of columns that will fit across one line are printed in a strip. As many strips are used as are required to print the entire model width. This format is called the strip format. The absolute value of the format code specifies the printout format as follows.

0 - (10G11.4)                                                                      3 - (15F7.1) 6 - (15F7.4)  9 - (20F5.2)  12 - (10G11.4)

1 - (11G10.3)                                                                      4 - (15F7.2) 7 - (20F5.0)  10 - (20F5.3)      

2 - (9G13.6)                                                                        5 - (15F7.3) 8 - (20F5.1)  11 - (20F5.4)      

Note if IREALSL…0, the drawdown is a meaningless parameter and this switch is associated with printout of saturations of the active phase.

IHEDUN                                                                            --          is the unit number to which heads will be written if they are saved on disk.

IDDNUN --                                                                        is the unit number to which drawdowns will be written if they are saved on disk.

Note if IREALSL…0, the drawdown is a meaningless parameter and this switch is associated with printout of saturations of the active phase. Thus, active phase saturations are printed to the DDN files when IREALSL…0.

ICONFM                                                                            -

-

is a code for the format in which concentrations are to be printed. The format codes are listed above (see the explanation of IDDNFM).

 TCNTNU – Flag and time value for continuation of a previous simulation. If TCNTNU = 0, the simulation is NOT a continuation of a previous simulation. If TCNTNU > 0, the simulation continues from a previous simulation, from the time value provided for this input.

DTCNTNU -  time-step size value used for the continuation simulation from time value of TCNTNU. If the simulation is not a continuation simulation, this variable is not required/used.

ICNTNU – Flag and unit number for reading restart heads from an external file created by a previous simulation.

If ICNTNU = 0, do not look for restart heads in a file created by a previous simulation.

If ICNTNU > 0, it is the Fortram unit number from which to read heads for the restart simulation (on file with extension HDC). The file will be wound to the appropriate time value (TCNTNU) to obtain the required starting heads.

ICNTNUC –  Flag and unit number for reading restart concentrations from an external file created by a previous simulation (when a transport simulation is performed).

If ICNTNUC = 0, do not look for restart concentrations in a file created by a previous simulation.

If ICNTNUC > 0, it is the Fortram unit number from which to read concentrations for the restart simulation (on file with extension COC). The file will be wound to the appropriate time value (TCNTNU) to obtain the required starting concentrations.

IUNITMB --  Flag and unit number for saving nodal mass balance terms to a binary file.

If IUNITMB # 0 Node-by-node flux balances are not saved.

If IUNITMB > 0 Node-by-node flux balances are saved to a binary file with default extension NFB and the value of IUNITMB is the Fortran Unit number on which the file is written.

INCODE  --  is the head/drawdown output code. It determines the number of records in input item 5.

If INCODE < 0,                                                                 layer-by-layer specifications from the last time steps are used. Input item 5 is not read

If INCODE = 0,                                                                 all layers are treated the same way. Input item 5 will consist of one record.

If INCODE > 0,                                                                 input item 5 will consist of one record for each layer.

IHDDFL                                                                            –    is a head and drawdown output flag.

If IHDDFL = 0,  neither heads nor drawdowns will be printed or saved on disk.

If IHDDFL … 0, head and drawdowns (or saturations) will be printed or saved according to the flags for each layer specified in input item 5.

IBUDFL  --  is a budget print flag for flow simulation.

If IBUDFL < 0, overall volumetric budget will be printed every “n” timesteps where n=abs (IBUDFL)

If IBUDFL = 0,  overall volumetric budget will not be printed.

If IBUDFL > 0, overall volumetric budget will be printed as specified by either NPRT or NPSTP in Item 5 input.

(Note that the overall volumetric budget will always be printed at the end of a stress period, even if the value of IBUDFL is zero.)

ICBCFL  --  is a cell-by-cell flow-term flag.

If ICBCFL = 0,  cell-by-cell flow terms are not printed or saved.

If ICBCFL … 0,  cell-by cell flow terms are printed or recorded on disk depending on flags set in the component of flow packages, i.e., IWELCB, IRCHCB, etc.

DELT                                                                                 --   is the first time step size of the stress period.

TMIN                                                                                 --   is the minimum value of time step size allowed in the stress period.

TMAX                                                                               --   is the maximum value of time step size allowed in the stress period.

TSMULT                                                                            –    is the time step multiplier. TSMULT should NOT be less than 1.0. This value overwrites the TSMULT value read in item 9 of Basic Package Input.

TSDIV                                                                               --   is the reduction factor of time step size. TSDIV must be greater than 1.0.

NPRTS                                                                               –    is the number of print times specified in item 4. If no print time is specified (i.e., NPRTS=0), overall volumetric budget and heads/ drawdowns will be printed or saved at the end of the stress period                                           according to the flags specified in input item 5. If NPRTS=0, item 4 must be skipped.

NPSTP                                                                               –    is an output printout control parameter. NPSTP=n means that the overall volumetric budget, heads/drawdowns, and cell-by-cell flow terms will be printed or saved at every nth time step of the stress period according to the flags specified in input items 2 and 5. Set NPSTP equal to zero if the step wise printing/saving is NOT required.

TMVEC                                                                             --   is a 1-D array containing NPRTS number of print time values at which printing or saving of volumetric budget, heads/drawdowns, or cell-by-cell flow terms is desired in a particular stress period. The printing/saving is performed according to the flags specified in input items 2 and 5. Time is assumed to be zero at the beginning of a stress period. The print times must be specified in ascending order of magnitude and the difference between two successive values must be at least TMIN.

                                                                                                Note: This item must be skipped if NPRTS=0.

Hdpr                                                                                   --   is the output flag for head printout.

                                                                                                 If Hdpr = 0, head is not printed for the corresponding layer.

                                                                                                 If Hdpr … 0, head is printed for the corresponding layer.

Ddpr                                                                                   --   is the output flag for drawdown printout.

                                                                                                 If Ddpr = 0, drawdown is not printed for the corresponding layer.

                                                                                                 If Ddpr … 0, drawdown is printed for the corresponding layer.

Hdsv                                                                                   --   is the output flag for head save.

                                                                                                 If Hdsv = 0, head is not saved for the corresponding layer.

                                                                                                 If Hdsv … 0, head is saved for the corresponding layer.

Ddsv                                                                                   --   is the output flag for drawdown save.

                                                                                                 If Ddsv = 0, drawdown is not saved for the corresponding layer.

                                                                                                 If Ddsv … 0, drawdown is saved for the corresponding layer.

Conpr                                                                                 --   is the output flag for concentration printout.

                                                                                                 If Conpr = 0, concentration is not printed for the corresponding layer.

                                                                                                 If Conpr … 0,  concentration is printed for the corresponding layer.

Consv                                                                                 --   is the output flag for concentration save.

                                                                                                 If Consv = 0, concentration is not saved for the corresponding layer.

                                                                                                 If Consv … 0,  concentration is saved for the corresponding layer.