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

5.1 General

在MODFLOW代码中，时间域（用于瞬态流动模拟）使用反向差分公式进行离散化，通过指定以下参数预先确定时间步长值：时间步数（NSTP），应力周期长度（PERLEN）以及一个乘法因子（TSMULT），该因子使用几何级数递增时间步长大小。使用这种方案的局限性如下。

• 如果矩阵求解器或非线性迭代方案（用于非均质情况）在给定的时间步长内无法收敛，则计算将中止。

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

• 必须保持时间步长大小的几何级数。因此，用户在分配要检查的流系统的特定时间值时缺乏灵活性。

因此，只能在位于时间推进系列上的时间值处请求输出，这需要用户提前确定。

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

模拟结果的时间值对于自动时间步进方案是未知的；因此，ATO4 Package还包括用于输出控制的新模块，该模块确定所需输出的时间值。调用ATO4 Package后，MODFLOW先前的输出控制（OC）Package的使用变得不再必要。请注意，原始MODFLOW的输出控制实施起来繁琐且麻烦。另一方面，新方案（ATO4）的输出控制简单而直接。用户只需输入希望在任何给定应力周期的哪些时间值（以前使用的时间单位）进行输出以及所请求输出的详细信息（即是否输出任何或全部水头、降水、质量预算和单元间流动术语），ATO4 Package将自动识别所选输出的时间值并调整时间步长大小以在这些打印时间进行计算。因此，例如，用户可以在15个月的应力周期模拟后的10、11和15个月时检查结果。ATO4 Package不仅将在整个模拟中提供高效的时间步进，而且还将确保在请求输出的10、11和15个月进行计算的时间值。

最后，ATO包提供了从非零起始时间值重新启动模拟的选项。MODFLOW的OC包不允许重新启动，并且在涉及多个应力周期时，为执行重新启动而开发新的输入文件可能会变得非常繁琐。

5.2 Formulation of Adaptive Time-Stepping Used in the ATO4 Package

ATO4 Package中实施的自适应时间步长选择方案简单而有效。用户提供了每个长度为PERLEN的应力周期的初始时间步长（DELT），最大时间步长（TMAX），最小时间步长（TMIN），时间步长放大因子（TSMULT）和时间步长减小因子（TSDIV）。初始时间步长DELT用于应力周期的第一个时间步。如果在最大迭代次数（MXITER）的35%内实现了收敛，则将时间步长增加TSMULT的乘法因子。对于MXITER的35%到65%之间的收敛率，时间步长与先前的DELT保持不变，并且系统的收敛行为被认为对于DELT的此值是最佳的。然而，如果实现收敛所需的迭代次数超过了MXITER的65%，则将时间步长减小TSDIV的除法因子。如果无法在特定时间步中实现收敛，则将时间步长减小5.0倍，并使用较小的时间步长大小重新进行计算。此外，如果达到打印时间（即需要打印/保存输出的目标时间值），则选择时间步长以包括打印时间值。如果将DELT的值减小到小于TMIN，则模拟将中止。如果遇到严重困难，这种技术提供了退出程序模拟的手段。最后，DELT的值受到TMAX的限制，这对于时间离散化的准确性是重要的，即使非线性程度较轻或不显著。用户可以选择使用TMVEC数组在指定的时间打印，或在指定的时间步数（NPSTP）之后，或两者都可以。

5.3 Verification Example

为验证ATO4 Package并展示其实施，进行了一次非均质含水层中的抽水瞬变模拟。所选问题在BCF4 Package的文档中有描述，模拟使用了原始MODFLOW的时间步进和输出控制。在计算效率和鲁棒性方面进行了这两个模拟运行的比较。

在这个例子中，进行的瞬变分析考虑了一个厚度为300英尺的非均质含水层，如图5.1所示。建模域对应于一个边长为75,000英尺的正方形。含水层的顶部和底部分别位于海拔50英尺和-250英尺。整个域内的含水层初始地下水位为零。域内受到3.0 x 10^–9 ft/s的均匀连续补给。图中左侧（西侧）边界是一个恒定（零）水头边界，其余边界为不流动边界。在含水层底部的100英尺处设置了15口井，并且位置如图5.1所示。每口井以3.85 ft^3/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 

该区域被离散成3个层，每个层厚100英尺，由均匀方块的15行15列组成，面积为5,000平方英尺。因此，MODFLOW所需的VCONT值（垂直水力导数除以两个层之间的间隔）计算为10^(-7) s^(-1)。每个层的恒定主要储存系数（比流特定存储乘以块厚度）为10^(-4)。在顶部和中间离散层的西侧规定了零头。最初，假定域中的水头为零。在第一次4年内允许位于第3层的15口井以3.85立方英尺/秒的速率抽水。在随后的12年中停止抽水，以使含水层恢复。

进行了两次独立的瞬态模拟运行，分别使用（1）现有的MODFLOW时间步长和输出控制，以及（2）使用新的自动时间步长和输出控制（ATO4）包。这两个模拟分别被称为Case（1）和Case（2）。对于两种情况，选择了输入LAYCON = 43的非约束流选项。在两种情况下都使用SIP包来解决矩阵方程。Case（1）的时间步控制参数包括：MXITER = 50，DELT = 7.875 x 10^6 s，