第四章 模型维度和边界的设置

. if one advances confidently in the direction of his dreams. he will meet with success
unexpected in his common hours. He will put something behind, will pass an invisible
boundary.
Henry David Thoreau

"如果一个人充满信心地朝着他的梦想前进，他将在日常生活中遇到意想不到的成功。他将留下一些东西，越过一道无形的界限。" - 亨利·大卫·梭罗

Contents 

4.1 Spatial Dimensions 118 

4.1.1 Two-Dimensional Models 118 

4.1.1.1 Areal Models 118 

4.1.1.2 Profile Models 125 

4.1.2 Three-Dimensional Models 133 

4.2 Selecting Boundaries 134 

4.2.1 Physical Boundaries 136

 4.2.1.1 Contacts with Geologic Units of Low Hydraulic Conductivity 136 

4.2.1.2 Shear Zones and Faults 137

 4.2.1.3 Surface Water Features 137 

4.2.1.4 Freshwater–Seawater Interface 140 

4.2.2 Hydraulic Boundaries 144 

4.2.2.1 Streamlines (Groundwater Divides) 145 

4.2.2.2 Equipotential Lines (Constant Head/Constant Flow) 145 

4.3 Implementing Boundaries in a Numerical Model 145 

4.3.1 Setting Boundaries in the Grid/Mesh 145

 4.3.2 Specified Head Boundaries 147 

4.3.3 Specified Flow Boundaries 148 

4.3.4 Head-dependent Boundaries 152 

4.3.4.1 Surface Water Bodies 153 

4.3.4.2 Drains 155 

4.3.4.3 Evapotranspiration (ET) from the Water Table 157 

4.3.4.4 Lateral Boundary Flows and Distant Boundaries 158 

4.4 Extracting Local Boundary Conditions from a Regional Model 159 

4.5 Simulating the Water Table 162 

4.5.1 Fixed Nodes 164

4.5.2 Movable Nodes 165 

4.5.3 Variably Saturated Codes 166 

4.6 Common Modeling Errors 170 

4.7 Problems 170

 References 176

目录

4.1 空间维度 118

4.1.1 二维模型 118

4.1.1.1 面积模型 118

4.1.1.2 剖面模型 125

4.1.2 三维模型 133

4.2 选择边界 134

4.2.1 物理边界 136

4.2.1.1 与低水力导流性地质单元的接触 136

4.2.1.2 剪切带和断层 137

4.2.1.3 地表水特征 137

4.2.1.4 淡水-海水界面 140

4.2.2 水头边界 144

4.2.2.1 流线（地下水分水岭） 145

4.2.2.2 等势线（定水头/定流） 145

4.3 在数值模型中实施边界 145

4.3.1 在网格/网格中设置边界 145

4.3.2 指定水头边界 147

4.3.3 指定流边界 148

4.3.4 水头依赖的边界 152

4.3.4.1 地表水体 153

4.3.4.2 排水 155

4.3.4.3 蒸散蒸发（ET）从水位表 157

4.3.4.4 侧向边界流和远距离边界 158

4.4 从区域模型中提取本地边界条件 159

4.5 模拟水位表 162

4.5.1 固定水头节点 164 4.5.2 可移动水头节点 165

4.5.3 可变饱和度代码 166

4.6 常见建模错误 170

4.7 问题 170

参考文献 176

Boxes 

Box4.1 Two-Dimensional or Three-DimensionaldMore about the D-F Approximation 122 

Box 4.2 Profile Models 125 

Box 4.3 Spreadsheet Solution of a Finite-Difference Profile Model 128 

Box 4.4 The Freshwater–Seawater Interface 141 

Box 4.5 Large Water Budget Errors Arising from an HDB 153 

Box 4.6 What Controls the Water Table? 167

盒子

盒子4.1 二维或三维：关于D-F逼近的更多信息 122

盒子4.2 剖面模型 125

盒子4.3 有限差分剖面模型的电子表格解决方案 128

盒子4.4 淡水-海水界面 141

盒子4.5 由于HDB而产生的大型水预算错误 153

盒子4.6 什么控制着水位？ 167

4.1 空间维度

4.1.1 二维模型

4.1.1.1 平面模型

在二维（2D）平面模型中，流动在两个水平维度中进行模拟，以表示承压或非承压含水层中的流动（图4.1(a)和(b)）。

方程（3.13a）是承压含水层的控制方程，方程（3.13b）是非承压含水层的控制方程。模型中表示的所有要素的水力效应（例如，抽水井，内部水源，如地表水，以及周边边界条件）完全穿透含水层，即贯穿含水层的整个深度。二维平面模型的稳态解是地图视图中的一系列头部的二维数组；瞬态解计算每个时间步的头部数组。在承压含水层中，头部表示势面；在非承压含水层中，头部表示水位面。二维平面模型还可能在模型域的不同区域表示非承压和承压条件（图4.1(c)）。

平面模型仅由侧边界定义。模型的顶部和底部边界是在流动的二维表达中固有的，不由用户指定。水平流动的要求意味着二维平面模型的底部本质上是无流边界。承压含水层在底部由承压层或其他相对不透水的材料界定（图4.2）。承压含水层的顶部也由承压层界定。如果出于模型目的假定承压层的顶部或底部是不透水的，则在控制方程（3.13a）中没有穿越承压层的泄漏。如果承压层是渗透的，方程中的R表示泄漏，qL，即水通过承压层的垂直流动。

其中是源水位，可以是在封闭层上方的非承压含水层（图4.2）或另一个承压含水层； 是含水层中计算得到的水位；和分别是封闭层的垂直水力导度和厚度。泄漏可以发生在位于感兴趣的含水层上方和下方的含水层。参数 被定义为垂直渗透度。垂直渗透度的倒数是垂直阻力，即。尽管在概念上类似，垂直阻力具有直观的特性，即当厚度接近零时，阻力趋近于零，与渗透性相反，当 接近零时，渗透性趋近于无穷大。

非承压含水层在底部本质上由被视为不透水的材料界定，在顶部由水位面界定。在底部边界上的泄漏（在底部边界不完全不透水的情况下）可以通过方程（4.1）进行模拟，或者通过从现场数据中指定 进行模拟。水位面不被视为边界，因为它是通过计算得到的解决方案。

二维平面模型使用杜皮-福尔赫默（D-F）逼近法。简而言之，D-F逼近法假定流动主要是水平的（见盒子4.1中的图B4.1.1），并且垂直水力梯度可以忽略不计（尽管仍然可以表示垂直流动；参见盒子4.1）。当垂直方向上的水头变化为零（∂h/∂z=0）时，水平平面上任何点（x，y）处的水头在有限含水层中等于势面标高，在非有限含水层中等于水位面标高（从含水层底部测量）（图4.3）（即，水位面的标高等于饱和厚度，b）。使用D-F逼近法的解与三维（3D）流动的解之间的比较表明，在距离引起垂直流动的水力要素（例如，3D要素，如部分穿透的地表水体、抽水井和地下水分水岭）大于2.5d的距离处，D-F模型计算的水头几乎无法与3D模型中的水头区分开（Haitjema，2006，p. 788），其中：

在方程（4.2a）中，\(K_h\) 和 \(K_v\) 分别是水平和垂直水力导度，\(b\) 是饱和厚度。因此，在各向同性条件下（\(K_h/K_v = 1\)），在距离一个三维要素两到三个饱和厚度之外，垂直流动可以忽略不计。在区域地下水流系统中（例如，见盒子4.1中的图B4.1.1），导致垂直流动的水力要素位于系统边界，这使得对于单个水力要素所需的距离加倍，并且意味着当系统的长度 \(L\) 大于 \(5d\) 时，D-F流是一个很好的近似值（Haitjema，2006，第788页）：

在许多地下水系统中，距离 \(L\) 远大于饱和厚度 \(b\)，因此使用D-F逼近的平面流模型适用于各种各样的地下水流问题（盒子4.1）。

在D-F逼近法下，不存在渗流面。渗流面是非承压含水层中的一个放水面，其中压力等于零（大气压力）；沿着渗流面的水头等于面的标高。渗流面可能在山坡、河岸（图4.3）以及隧道、矿井和大口径井中形成（第6.2节）。D-F模型在靠近放水面的地方性能较差，水位面的坡度较大，有强烈的垂直流分量，但当距离渗流面足够大时，根据上述讨论的准则定义，解是准确的。

准确模拟水位表作为自由表面（移动边界）是复杂的，因为它需要对垂直流动进行严格的表示，并在水位表上施加非线性边界条件（Neuman和Witherspoon，1971；Diersch，2014，第216-218，405-406页；第4.5节）。正如Diersch（2014，第406页）所观察到的，实际上，“自由表面问题通常只以非严格的方式解决”。

D-F模型是一种非严格的方法（盒子4.1）。有关模拟水位表和相关渗流面的其他选项将在第4.5节中讨论。

图4.3 在D-F逼近法下非承压含水层中的水平流动（蓝色箭头）。在放水面附近和水位表附近，D-F逼近法不准确。真实流场中的等势线（红色线）在水位表和放水面附近偏转，存在垂直梯度和渗流面；D-F逼近法假设垂直等势线（显示为蓝色）。在真实流场中，水位表在放水面之上与地表水体的自由表面相交（红线），形成一个渗流面。

在D-F逼近法下，水位表是连续的，与地表水体的自由表面相遇，没有渗流面。

4.1.1.2 剖面模型

剖面模型表示地下水流系统的垂直切片（横截面）中的二维流动（见盒子4.2中的图B4.2.1）。其控制方程为方程（B4.2.1）（盒子4.2）。该模型由剖面的顶部、底部和侧边的边界条件来定义。

盒子4.2 剖面模型

该盒子介绍了使用标准有限差分（FD）和有限元（FE）地下水流代号建立剖面模型的说明。为了方便起见，参考使用了FD网格。

剖面模型的纵轴必须与地下水流路径平行，因为所有水必须在剖面内流动；不能有水垂直于剖面流动。剖面的厚度设置为一个单位，或者设置为与流动平行于剖面的含水层宽度（图B4.2.1）。通常，剖面是一个矩形，有四个边界（即顶部、底部和两侧边界），但当剖面不是矩形时，可以适应更复杂的边界几何形状。通常，顶部边界表示水位表，底部边界被模拟为无流边界，但也可以使用指定水头、指定流量或水头相关条件进行模拟。剖面两端的侧边界通常表示地下水分水岭，并被模拟为无流边界条件（盒子4.3）。剖面模型通常被构建为3D模型的多层切片，我们称之为切片取向，但也可以构建为单层模型，我们称之为层取向。

图B4.2.1 剖面模型，其纵轴与地下水流平行，由紫色箭头表示。

还显示了水位等值线（以米为单位）。在三维模型中模拟的切片取向是剖面建模的首选取向。

图B4.2.2 在层取向中（底部三个图），剖面被模拟为一个面积上的二维模型。层的厚度等于剖面的宽度。为了比较，顶部的图中显示了切片取向。

剖面在x-z平面内的流动的瞬态控制方程由方程（3.12）推导而来，其中 ∂h/∂y = 0：

其中（1/T）是应用于剖面顶部单元（）的体积充水速率（）。瞬态模拟很繁琐，因为流动路径随时间变化，剖面的取向必须调整以符合变化的流场（图B4.2.1）。因此，3D模型更适用于瞬态模型，而剖面模型通常模拟稳态流动，其中在方程（B4.2.1）中 。

切片取向 切片取向是剖面建模的自然取向。切片取向中的剖面简单地是3D模型的一个横截面；剖面模型具有多个层，剖面的厚度等于Δy或Δx，具体取决于取向（例如，图B4.2.2）。3D控制方程（方程3.12）适用，含水层参数（Kx = Ky，Kz，Ss）和层厚度输入方式与3D模型相同。充水量根据FD或FE代号中提供的充水选项输入。由于垂直切片包括多个模型层，模型输入包括逐层输入表示该切片的单行数据；因此，与层取向相比，数据输入略显繁琐，但在GUI中，数据的组装更加方便。与3D模型一样，水位节点的贮水率等于比流系数，而水位以下的节点的贮水率表示由于含水层的压缩和水的膨胀而释放的水（第5.5节；图5.27）。

层取向

在层取向中（图B4.2.2），剖面模型的网格/网格设置为一个面积上的二维模型，因此需要模型师在视角上进行改变。层取向依赖于方程（B4.2.1）和方程（3.13a）之间的相似性，方程（3.13a）是受限含水层中的二维水平流的控制方程。如果方程（B4.2.1）中的 z 等于方程（3.13a）中的 y，这两个方程是相同的。类似地，方程（B4.2.1）中的 Kx、Kz、Ss 和 W* 分别代替方程（3.13a）中的 Tx、Ty、S 和 R。层的厚度 b 等于剖面的宽度，使得 Tx = Kxb 和 Ty = Kzb；S = Ssb 和 R = W*b。参数 S 是贮水率，对于水位节点等于比流系数，在其他地方表示由于含水层的压缩和水的膨胀而释放的水（第5.5节；图5.27）。

层取向的主要优势在于，原始模型输出可以被视为一个2D头值数组，而不是对于多层逐行查看。在使用电子表格解决剖面模型时，层取向也更为方便（盒子4.3）。然而，在GUI中，切片取向有助于剖面建模，因为它们在横截面中显示模型的输入和输出。由于大多数应用建模都是通过GUI完成的，切片取向是大多数剖面模型的首选表示。

剖面模型的解是一个垂直平面中的2D头值数组，可以从中绘制等势线并推断流向（盒子4.3的图B4.3.3(c)；图4.4）。剖面模型假设沿着剖面没有垂直流动；它们沿着流线定向，所有流动发生在剖面的垂直平面内。

图4.4 剖面模型中的等势线（淡灰色线）和流路径（浓重蓝色线），带有示意的流动箭头。在上层底部，水力导数有两个数量级的对比有效地创建了无流边界。图是使用TopoDrive（Hsieh, 2001）创建的。

盒子4.3 有限差分剖面模型的电子表格解法

通过近似地简化地下水流方程和边界条件而导致的FD方程的迭代解可以借助电子表格获得。在电子表格建模中，电子表格中的每个单元格都是一个FD单元格。该技术适用于简单的问题，例如Toth（1962）描述的均匀、各向同性含水层横截面中的稳态区域地下水流问题。

Toth的概念模型如图B4.3.1所示，数学模型如图B4.3.2所示。

线性水位表形成了顶部边界，并由指定的水头边界表示（有关水位表上指定水头条件的注意事项，请参见盒子4.6中的警示讨论）。侧边界表示由无流边界模拟的区域地下水分水岭。系统底部的无流边界表示不透水材料。

控制方程是拉普拉斯方程：

其中 z 是垂直坐标。对于 Dx 和 Dz 都是常数且相等的均匀规则网格的方程（B4.3.1）的有限差分逼近是：

其中 m 代表当前迭代水平的头，m+1 代表新迭代水平的头。方程（B4.3.2）是以适合通过Gauss-Seidel点迭代解 的形式编写的，其中在计算时，头部会被更新到下一个迭代水平。请注意，每个FD单元格中的水头等于相邻四个单元格中水头的平均值，而 的解与水力导数和网格间距无关。在方程（B4.3.2）中表示的五个节点组成了2D FD计算器，称为五点有限差分星操作符（图3.4）。

图B4.3.1 展示了一个区域地下水系统横截面中的边界和示意流线的概念模型（参考 Toth, 1962）。

图B4.3.2 数学模型显示了图B4.3.1中概念模型的控制方程和边界条件（参考 Toth, 1962）。

在层取向的剖面模型（盒子4.2）中，使用电子表格建立（图B4.3.3）。对于我们版本的Toth的问题，我们假设在(x,z) = (0,100)的出水区头等于100米，并且水位线的坡度为0.05米/米。我们使用20米的节点间距（Dx = Dz = 20米），因此电子表格有七行和13列，其中第7行和第1列以及第13列（图B4.3.3(b)和(c)中的A列和M列）包含域外的虚拟（或“虚构”）节点。虚拟节点用于通过将虚拟节点上的头设置为边界另一侧的计算头来实现无流边界条件。这样，跨越边界的头梯度为零。电子表格设置为块中心的FD网格，以便无流边界位于边界单元的外侧边缘。通过为边界单元使用不同的FD方程，可以使用点中心的网格（问题P4.1）。

电子表格中每个单元格的FD方程都是根据方程（B4.3.2）（图B4.3.3(b)）编写的。例如，电子表格的第3行第4列（单元格D3）中的方程是：

单元格沿边界的方程被修改以包含边界条件。电子表格的顶部行包含从100米到110米的指定头值，这些是沿顶部（水位线）边界的头，使用图B4.3.2中的线性方程计算。

图B4.3.3 图B4.3.2问题的数值模型。 (a) 块中心网格显示边界位置和问题域内11列和6行节点的放置。 (b) 问题的Excel电子表格模型中每个单元格的FD方程。指定的头边界值（以米为单位）输入到第一行以表示水位线。

列A和M以及第七行的虚拟节点在问题域之外，用于在边界沿着FD单元的外侧边缘实施无流边界条件。 (c) 显示以米为单位的头。为了计算水位线上的通量，水力导数等于10米/天（单元格B10）。水位线处的流量（Q）按横截面的1米宽计算为m3/天。还计算了水位线上的总补给（RTotal）和总放水（DTotal）以及水预算中的误差（即RTotal和DTotal之间的差异）。

通过将虚拟节点的头设置为跨越边界的相邻单元格中的头（图B4.3.3(b)），模拟了其他三个边界的无流条件。这样，第一列和最后一列以及最后一行的单元格中的头反映了边界上的头（图B4.3.3(c)）。角节点未包括在五点星计算器（方程（B4.3.2））中。一旦制定了适当的方程，就可以设置电子表格以迭代地解决每个单元格的头问题（请参见问题P4.1）。对于我们的问题，假设水力导数为10米/天，剖面（层）厚度为1米，计算水位线上的通量。这个问题的水预算表述为：水位线上的流入（RTotal）等于水位线上的流出（DTotal）。水预算中的误差是RTotal和DTotal之间的差异。头和通量的解决方案（图B4.3.3(c)）可以与使用点边界的解决方案进行比较（问题P4.1）。

电子表格建模是一种有用的教学工具，可以说明如何为每个FD单元编写方程并通过点迭代解决FD模型。感兴趣的读者可以在Bair和Lahm（2006）中找到有关地下水模型电子表格解决方案的更多信息。

然而，从实际的角度来看，即使在建模项目的初期阶段，电子表格解决方案也没有任何优势。为电子表格中的每个单元格输入和检查方程是乏味而繁琐的。设置和测试复杂电子表格模型所需的时间可能等于或大于设置和运行标准流程代码所需的时间。

Therefore, profile models cannot simulate radial flow to a pumping well or 3D flow around point sources and sinks or other 3D flow effects (Box 4.2).

An axisymmetric profile model is designed to simulate radial flow in a pie-spaced problem domain. Linear or point sources and sinks can be simulated as long as the source or sink is located at the origin of the profile (Fig. 4.5). The thickness of the profile widens from the origin (r = 0) according to an angle q. Although special purpose codes are available to solve the governing equation in radial coordinates (e.g., Reilly, 1984), standard codes based on Cartesian coordinates can also simulate axisymmetric profiles by the appropriate adjustment of input parameters (e.g., see Langevin, 2008 for details of the procedure for the MODFLOW family of codes).

Profile models have been used to study regional groundwater flow in cross section (e.g., Cardenas, 2008; Winter, 1976; Freeze and Witherspoon, 1967) and hyporheic interchange beneath a stream (e.g., Woessner, 2000). Axisymmetric profile models are typically used to simulate aquifer (pumping) tests; Langevin (2008) used axisymmetric models to simulate radial flow to a pumping well, a push-pull test, and upconing of saline water as a result of pumping. However, in practice most applications of groundwater models to engineering and management problems use a 3D model to simulate 3D flow effects because it usually is not reasonable to restrict flow to one horizontal dimension or to assume radial flow throughout the problem domain.

图4.5 轴对称剖面。(a) 轴对称剖面模型的有限元网格，用于模拟向腔体（即，一个中央穿孔间隔的井，如左侧所示）注入的瞬时地下水流（改编自Keller等人，1989)。(b) 轴对称剖面，用于模拟到部分穿透抽水井的流动，该井位于右侧FD网格所示的含水层扇形部分的点上。剖面的厚度和分配给单元格的传导度随着离井的距离而增加。还需要调整蓄积来反映剖面厚度的变化。抽水速率根据含水层楔形部分的角度进行调整（此处等于20°）。三维模型的网格假设径向流动，并使用对称性仅对含水层的一个四分之一进行建模，显示在图的左侧（改编自Land，1977)。

4.1.2 Three-Dimensional Models

三维（3D）模型模拟了所有三个空间维度中的水流，并可以明确包括影响地下水流的所有水文地层单元。通常，每个水文地层单元都是一个模型层（图2.7和图2.9(b)），但在需要计算水文地层单元内部头的垂直分布时，会使用更多的层次（第5.3节;见Box 5.3的图B5.3.2）。3D模型具有顶部、底部以及侧边界，用户必须指定所有周边边界的条件（即，指定头、指定流或头依赖条件;第3.3节）。

上边界通常是水位面。下边界通常是与相对不透水材料的接触，并由不渗流边界条件表示，或者可以通过指定头、指定流或头依赖条件（方程（4.1））模拟下边界的渗漏。支配方程是方程（3.12），解是瞬时模拟的每个时间步的3D头数组。对于稳态问题，解是单一的3D头数组。

在过去，有时使用准3D模型来模拟区域尺度的地下水系统，但完全的3D模型在很大程度上已经取代了准3D模型。

在准3D模型中（图4.6），每个含水层都由一个模型层代表，但约束层未明确包含在网格/网格中。相反，通过泄漏项间接模拟由介质约束层的存在引起的垂直流动阻力。由约束单元引起的垂直流动阻力由垂直渗透率（K_z'/b'）表示，模型中的含水层之间的垂直流量为：

\[ \]

其中，是约束层上方的含水层中的水头，是约束层下方的含水层中的水头；是约束层的厚度，是其垂直水力导数。方程（4.3）假定垂直阻力主要由控制，因此含水层中的流动主要是水平的。该模型被认为是准3D，因为含水层（模型层）中的流动是水平的，但是模型层之间的垂直流动由方程（4.3）表示，而无需在网格/网格中明确表示约束层。

在准3D模型中，流动被假定严格垂直通过约束层，并且没有从约束层内部释放水的情况（例如，由于从约束层上方和/或下方的含水层抽水引起）。当含水层的水力导数比约束层小两个数量级以上时（图4.4），通过假定约束层中的垂直流动引入的误差通常小于5％；然而，在许多实际应用中，忽略约束层中的储存可能引起显著误差（Steltsova，1976）。可以使用包括约束层比释放的瞬时泄漏项来模拟约束层中的存储（有关详细信息，请参见Bredehoeft和Pinder，1970；Trescott等人，1976；Anderson和Woessner，1992，第198页）。准3D模型不应与粒子跟踪一起使用，因为在网格或网格中没有约束层会导致计算流路径和传播时间时的错误（第8.2节）。

在大多数应用建模问题中，准3D模型的局限性表明对于大多数应用，大多数情况下都应使用完全3D模型。然而，准3D方法可以有效地模拟由一系列含水层和约束层组成的特定多含水层系统（例如，见Box 4.4的图B4.4.1(b)）。

图4.6 示意图显示了美国南卡罗来纳州萨凡纳河站的七层准3D模型的水文地质和模型层。约束层未在FD网格中表示，但通过用泄漏项表示约束层的垂直阻力以及通过泄漏项表示层间流动，间接地包含在模型中（Clark和West，1998）。

4.2 SELECTING BOUNDARIES

Physical and hydraulic boundaries of the conceptual model were discussed in Section 2.3 and the three types of mathematical boundary conditions (specified head including constant head; specified flow including no flow; and head-dependent conditions) were introduced in Section 3.3. Conceptually, boundary conditions define flow into and out of the problem domain. Boundary conditions that refer to hydraulic conditions along the perimeter of the model domain are called perimeter boundary conditions (Fig. 4.7(a)).

However, boundary conditions in groundwater flow codes can also represent internal sources and sinks (Fig. 4.7(b)), which in effect are internal boundary conditions. Perimeter boundary conditions are required in finite-difference (FD) and finite-element (FE) models. Only internal boundaries are specified in analytic element (AE) models (Section 3.4) because the perimeter boundary in an AE model is automatically set at infinity.

The head along a specified head boundary is fixed, which anchors the hydraulic gradient and control flux at the boundary. A specified head can have a large influence on a simulation because it has an infinite capacity to supply and remove water from the model. A specified flow boundary fixes the flux at the boundary so that boundary flows are not influenced by changes in head at the boundary. A head-dependent boundary (HDB) uses an assigned boundary head and the simulated head at the boundary to control flows across the boundary. In this section, we re-visit these three types of mathematical boundary conditions and relate them to the physical and hydraulic boundaries of the conceptual model (Section 2.3). Both physical and hydraulic boundaries can be represented by all three types of mathematical boundary conditions. In Section 4.3 we discuss how to implement the three types of mathematical boundaries in a numerical model.

Of course, boundaries in the numerical model should be consistent with boundaries in the conceptual model. Whenever possible, physical boundaries should be used instead of hydraulic boundaries. Usually physical boundaries can be found by increasing the size of the model domain to extend out to physical boundaries present in the field (Section 4.4).

After boundaries are selected, the modeler should visualize the general flow pattern that will be induced by the conditions assigned to the boundaries and confirm that the flow pattern makes hydrogeologic sense and that inflow and outflow locations are consistent with what is known about the hydrogeology of the modeled area.

Figure 4.7 Physical boundaries. (a) Two-dimensional areal FD grid showing perimeter boundaries defined by physical boundaries. Relatively impermeable bedrock outcrops are no-flow boundaries; specified head boundaries represent wetlands, Lake Wausau, and the Eau Claire River, WI, USA (modified from Kendy and Bradbury, 1988). (b) Relatively impermeable bedrock across a fault creates a physical boundary for the alluvial aquifer. The intermittent stream is separated from the water table by a thick unsaturated zone but contributes water to the aquifer via percolating conditions (Fig. 4.16(d)) (Niswonger and Prudic, 2005).

4.2.1 Physical Boundaries

Physical boundaries are the most robust and defensible type of perimeter or internal boundary as they represent physical features that are easily identified in the field (Fig. 4.7). Physical boundaries include contacts with units of low hydraulic conductivity (e.g., between valley fill and bedrock or between permeable sandstone and less permeable shale), surface water features such as rivers, lakes, wetlands, and in some coastal settings the interface between freshwater and seawater.

4.2.1.1 Contacts with Geologic Units of Low Hydraulic Conductivity

Low hydraulic conductivity units may form no-flow boundaries. Lateral no-flow boundaries may be placed at outcrops of relatively impermeable bedrock (Fig. 4.7(a) and (b)). The bottom boundary of most models is specified as a no flow condition at the contact with relatively impermeable material (e.g., Figs. 4.6, 2.3, 2.9(b), 2.12).

When there is at least two orders of magnitude contrast in hydraulic conductivity 136 Applied Groundwater Modeling between hydrogeologic units, flow along the contact is dominantly horizontal and the contact could be represented as a no-flow boundary (Fig. 4.4; Freeze and Witherspoon, 1967). The justification for placing a boundary at such a contact follows from Darcy’s law: given horizontal flow in two hydrogeologic units under the same horizontal hydraulic gradient, the unit with a transmissivity two or more orders of magnitude lower will convey less than 1% of the flow, which in most cases is sufficiently small that it can be neglected. In practice, the contrast in hydraulic conductivity at a geologic boundary is usually much greater than two orders of magnitude.

4.2.1.2 Shear Zones and Faults 

Faults and shear zones are fractures along which movement has displaced geologic materials (Fig. 4.7(b)). Usually hydraulic conductivities of the geological materials in and near the fault (shear) zone are altered. Geologic materials may be crushed to create higher conductivity breccia or lower hydraulic conductivity fine-grained gouge in the fault; zones adjacent to the fault have either enhanced or reduced hydraulic conductivity and/or, in some settings, the hydraulic conductivity is not appreciably changed (Scholz and Anders, 1994; Caine et al, 1996). Geological characterization of faults (e.g., Caine and Minor, 2009; Rawling et al, 2001) and analysis of groundwater flow patterns (e.g., Jeanne et al, 2013; Bense and Person, 2006; Seaton and Burbey, 2006) show that groundwater flow is often disrupted and that flow patterns through the fault zone can be complicated. In some cases a fault may be a conduit or a barrier to flow (Fig. 4.8) but in geologically complicated aquifers the same fault may act as both barrier and conduit (Bense and Person, 2006).

4.2.1.3 Surface Water Features 

Major surface water features, including oceans, large lakes, and streams, are physical boundaries that are perennial sources or sinks for groundwater. The stage in a major surface water body exerts hydraulic control and is not appreciably changed by stresses such as groundwater pumping (Figs. 4.7(a), 2.5(a), 2.9(b)). For example, a very large and deep lake is a physical boundary that robustly forms a constant head boundary for regional groundwater systems. A major river (Fig. 4.7(a)) can also form a robust physical boundary.

A major surface water feature that exerts strong hydraulic control may not physically penetrate the full thickness of the aquifer but nevertheless can hydraulically form a fully penetrating boundary if there is a groundwater divide beneath the feature (e.g., the rivers in Fig. 4.9(a), 2.1 and 2.6(a) and the Delaware Bay in Fig. 2.9(b)). When hydraulic control of the surface water feature is relatively weak, the divide does not fully penetrate the aquifer and some groundwater flows under the feature (e.g., the ditch in Fig. 4.9(b)). Zheng et al (1988) found that the formation of a fully penetrating divide beneath drainage ditches was inversely proportional to the 

regional slope of the water table and directly proportional to the head gradient across the streambed sediments. A combination of a specified head or HDB (to represent surface water) and no-flow boundary (beneath the surface water feature) could be used to represent a strong partially penetrating surface water body. However, partially penetrating surface water bodies away from the model perimeter are best represented by head-dependent conditions so that underflow beneath the surface water feature can change as hydrologic conditions change (Section 4.3).

Surface water bodies that are not directly connected to the groundwater system are not strong physical boundaries. For example, surface water may exert little or no control over groundwater flow in deep confined aquifers and surface water in arid settings may be hydraulically disconnected from the groundwater system owing to separation by a thick unsaturated zone (Fig. 4.7(b)).

Figure 4.8 Flow across fault zones shown in schematic cross sections of an unconfined aquifer. Simulations of the profile were done in TopoDrive (Hsieh, 2001). Equipotential lines (faint gray lines) and flowpaths (heavy blue lines) are shown. There is a two order of magnitude contrast between the hydraulic conductivity (K) of the aquifer and the fault zone. Both aquifer and fault are isotropic. (a) Fault as a conduit. The fault zone is shaded in blue; K of the fault is larger than K of the aquifer. Groundwater flows up the fault and discharges in the valley bottom. (b) Fault as a barrier (dam). The fault zone is shaded in pink; K of the fault is smaller than K of the aquifer. There is an abrupt drop in the water table across the fault and water is dammed against the fault.

When represented with specified head or head-dependent conditions, surface water features potentially supply an infinite amount of water to the model, which may be appropriate for large perennial bodies of water. However, ephemeral and headwater streams may dry up seasonally and failure to account for their seasonal disappearance may introduce significant error into the model (Mitchell-Bruker and Haitjema, 1996). Advanced options for simulating surface water features that allow for adjustment of surface water flow and stream or lake level are discussed in Chapter 6.

Figure 4.9 Fully and partially penetrating surface water bodies. (a) Schematic cross section through an unconfined aquifer showing groundwater divides beneath a topographic high and beneath a stream. The stream partially penetrates the aquifer physically but is hydraulically fully penetrating. The lake (at right) is both physically and hydraulically fully penetrating (Granneman et al, 2000). (b) Cross section through an unconfined aquifer showing streamlines in the vicinity of a shallow ditch. Streamlines flow beneath the ditch indicating underflow. The ditch is both physically and hydraulically partially penetrating (modified from Zheng et al, 1988).

4.2.1.4 FreshwatereSeawater

Interface In coastal aquifers, fresh groundwater discharges to the ocean along the freshwatere seawater interface (Fig. 4.10). Under field conditions, the interface moves in response to tidal fluctuations, groundwater pumping, and changes in recharge, creating a transition zone, or zone of dispersion. When the interface is relatively stable, the zone of dispersion is narrow and a no-flow boundary could be specified at a representative average position of the interface (Fig. 2.3). However, if the interface is not stable (or may not be stable under conditions represented in a forecasting simulation), a static boundary is likely not appropriate. When it is not appropriate to simulate the freshwatereseawater interface as a no-flow boundary, special purpose codes may be used to simulate a sharp interface without mixing between freshwater and seawater (Box 4.4) or to simulate mixing and flow in the interface by including density effects and dispersion (Sections 12.2 and 12.3).

Figure 4.10 Freshwatereseawater interface in a coastal aquifer showing the transition from freshwater to seawater in the zone of dispersion. The interface acts as a barrier to groundwater flow; freshwater flows upward along the interface and discharges to the ocean (Barlow, 2003).

Box 4.4 The Freshwater–Seawater Interface 

Groundwater is an important source of water for coastal communities throughout the world; hence, seawater intrusion is a critical concern, especially because its effects are largely irreversible. Groundwater in coastal aquifers flows seaward, mixes with subsurface saline water along the interface between freshwater and seawater, and discharges at the ocean floor (Fig. 4.10).

Groundwater pumping causes the interface to move inland resulting in seawater intrusion.

Mixing at the interface occurs by diffusion and dispersion creating a transition zone. Diffusion is caused by the difference in solute concentration between freshwater and seawater. Dispersion results from mixing of freshwater and seawater in response to transience in groundwater flow caused, for example, by intermittent pumping and tidal fluctuations. Such transience causes the interface to move back and forth.

Figure B4.4.1 The freshwatereseawater interface: (a) under hydrostatic conditions as assumed by the GhybeneHerzberg relation (Eqn (B4.4.1b)) (Barlow, 2003); (b) in a multiaquifer system simulated using a quasi-three-dimensional model (Section 4.1). The offset in the interface between aquifers (along EF in the figure) is small when vertical resistance between layers is small (i.e., leakance is large). The offset is relatively large when there is a confining bed between aquifer layers (Fitts et al, 2015).

There is widespread interest in the general topic of saltwater intrusion, which includes
seawater intrusion as well as movement of brines into inland freshwater aquifers (e.g.,
http://www.swim-site.org/). The literature on modeling saltwater intrusion encompasses a
variety of methods from relatively simple to complex (Werner et al., 2013). The simplest
expression for the freshwater–seawater interface, which derives from the well-known
Ghyben–Herzberg relation, states that the depth of the interface below sea level, z, is approximately 40 times the head above sea level, h (Fig. B4.4.1(a)). This simple mathematical statement follows from the observation that under hydrostatic conditions (i.e., both freshwater
and seawater are static) the head on the interface is the same from both the freshwater
and the seawater sides so that:

\rho_s gz=\rho_s g\left (  z+h\right )

where \rho_s is the density of seawater, \rho_f is the density of freshwater, and g is gravity. Rearranging
Eqn (B4.4.1a) gives the Ghyben–Herzberg relation:

For \rho_f=1.0 gm/cm3 (=gm/ml) and \rho_s=1.025 gm/cm3 , α= 40.

Of course, in reality neither freshwater nor seawater is static; both move upward and mix along the interface. Freshwater discharges to the ocean floor while seawater cycles from the ocean through the subsurface and discharges back to the ocean (as shown in Fig. 4.10).

Although there is some exchange of water between freshwater and seawater in the transition zone, the interface effectively acts as a no-flow boundary. When the interface is relatively stable and the transition zone is narrow it may be appropriate (depending on the modeling objective) to represent the interface as a static no-flow boundary where freshwater discharge to the ocean floor is simulated by a head-dependent or specified head boundary condition (Fig. 2.3).

However, for some problems it is necessary to forecast the transient movement (or the final long term average, equilibrium or steady-state, position) of the interface in response to a hydrologic stress such as freshwater pumping or rising sea level. There are two categories of models for forecasting the position of the interface: sharp interface (also called interface flow) models and variable density models. Both approaches allow for the contrast in density between fresh and seawater but the sharp interface approach assumes there is no mixing along the interface and no transition zone. A variable density model couples a variable density groundwater flow model to a model for solute transport based on the advection–dispersion equation (Sections 12.2 and 12.3). Simulating the freshwater–seawater interface with a coupled model requires small nodal spacing in the vicinity of the interface and typically involves long run times. Sharp interface numerical models, which typically run much faster than coupled models, are described below. There are also a number of analytical solutions that solve for the steady-state sharp interface; the reader is referred to Werner et al (2013; their section 4.1) and Bear and Cheng (2010, pp. 613–620) for summaries.

In the sharp interface approach, the governing equation for groundwater flow (Eqn (3.12)) is written in terms of freshwater head in the freshwater zone and native (seawater) head in the seawater zone (Fig. B4.4.2). The interface is simulated as a moving boundary between the two regimes where continuity in both head and flux are imposed along the interface. The mathematics are relatively complicated and the interested reader is referred to Bear and Cheng (2010, pp. 601–605) for details. The approach can be simplified by assuming horizontal flow following D-F theory (Section 4.1; Box 4.1). Then numerical codes for single density flow can be modified 

to solve for the interface, or interfaces in a multiaquifer system. Bakker et al (2013) used this approach to develop the SWI2 Package for MODFLOW, which solves for the transient movement of the interface in single aquifers (Fig. B4.4.2) as well as multiaquifer systems (Fig. B4.4.1(b)).

Fitts (2013, p. 224) points out that for 2D, transient simulations, storativity (Section 5.4) at the freshwater–seawater interface for both confined and unconfined aquifers should be represented as:

Sy=S = n_e\frac{ \rho_f }{ \rho_s- \rho_f}

where Sy is specific yield (i.e., storativity of an unconfined aquifer), S is confined storativity, and ne is effective porosity (Box 8.1). This definition of storativity is necessary because the landward movement of the interface replaces freshwater with seawater and effectively removes freshwater from storage.

Steady-state sharp interface models are also available. For example, Bakker and Schaars (2013) presented a D-F based approach to solve for the steady-state position of the interface; however, a disadvantage of their approach is that the top aquifer in a multiaquifer system cannot be unconfined. An approach similar to Bakker and Schaars (2013) is used in some AE codes where interface discharge potentials (Section 3.4) developed by Strack (1976; also see Bear and Cheng, 2010, pp. 616–619) are used to compute the steady-state position of the interface. GFLOW (Haitjema, 2007) uses Strack’s formulation to simulate interface flow in a single aquifer; AnAqSim (see review by McLane, 2012) can simulate both single and multiaquifer interface flow, and also allows the top aquifer in a multiaquifer system to be unconfined (Fitts et al, 2015).

 

Figure B4.4.2 Cross section of a sharp interface model as simulated with MODFLOWs SWI2 Package. (a) Conceptual model showing native (seawater) head and equivalent freshwater head at the interface. The freshwatereseawater interface (dotted line) separates freshwater (zone 1) from subsurface seawater (zone 2). (b) One-layer model of a coastal aquifer. The thickness of the layer varies in space; vertical variations in density within the layer represent freshwater and subsurface seawater. The code solves for the transient movement of the interface (Bakker et al, 2013).

4.2.2 Hydraulic Boundaries

When it is not possible or convenient to specify physical boundaries around the perimeter of the problem domain, the modeler might consider using hydraulic boundaries. Hydraulic boundaries are delineated by streamlines (groundwater divides) and equipotential lines (lines of constant head).

Ideally, hydraulic boundaries, like physical boundaries, define hydraulic conditions that are invariant in time and stable in space. However, hydraulic conditions in the field may change in response to stresses so that hydraulic boundary locations and definitions in the model require reevaluation. Hydraulic boundaries are more likely than physical boundaries to be specified inappropriately. Although hydraulic boundaries are never ideal, they are more suited for steady-state problems than for transient problems because conditions along hydraulic boundaries are likely to change as a result of transient stresses.

Given these issues, alternatives for setting boundary conditions (discussed in Section 4.4) should be considered before selecting a hydraulic boundary.

4.2.2.1 Streamlines (Groundwater Divides) 

Groundwater flows parallel to a streamline so that, by definition, no water crosses a streamline. Thus, streamlines can serve as groundwater divides that form hydraulic noflow boundaries (Figs. 4.11(a) and 4.12). Groundwater divides that coincide with regional topographic highs (Figs. 4.9(a) and 4.11(a)) are often relatively stable features and may form perimeter boundaries. But even regional groundwater divides can shift in response to changes in pumping and to a lesser extent to changes in recharge and changes in stage at regional sinks (e.g., Sheets et al, 2005). One way to determine if a hydraulic streamline boundary is affecting modeling results is to replace the no-flow boundary with specified heads and re-run the model. If flow to or from the specified boundary nodes is insignificant, the assignment of a hydraulic no-flow boundary at that location is appropriate (Problem P4.6).

4.2.2.2 Equipotential Lines (Constant Head/Constant Flow) 

An equipotential line, a line of constant head, may be used to form a constant head hydraulic boundary (Fig. 4.12), or specified flow rates may be calculated across the equipotential line and used to specify boundary flows. In practice, a specified head hydraulic boundary is the least desirable type of boundary condition because heads that are not tied to a physical feature are seldom stable in space or time. Moreover, a specified head boundary potentially provides infinite amounts of water to and/or removes infinite amounts of water from the model.

4.3 IMPLEMENTING BOUNDARIES IN A NUMERICAL MODEL 

This section presents methods for representing specified head, specified flow, and headdependent conditions in a numerical model both as perimeter boundary conditions and internal boundaries. We refer to FD grids and FE meshes as introduced in Section 3.5 and discussed in detail in Sections 5.1, 5.2, 5.3, but for convenience we use FD notation when presenting equations for implementing boundary conditions. Generally the same equations, with slightly different notation, apply to FE models. Although the boundary options discussed below are frequently used to represent streams, lakes, springs, and wetlands, more sophisticated options are available and are discussed in Sections 6.4 through 6.7. The presentation of methods in this section is general in nature; therefore, the modeler should always consult the code’s user’s manual for information about codespecific boundary implementation.

4.3.1 Setting Boundaries in the Grid/Mesh 

In a block-centered FD grid (e.g., MODFLOW and MODFLOW-USG) nodes are located in the center of the FD cell/block (Fig. 4.11(b)). Specified head boundaries are located on the nodes but specified flow boundaries are located on the edge of the block surrounding the node (Fig. 4.11(b); Fig. B4.3.3(a) in Box 4.3). In a point-centered FD grid (e.g., HST3D; Kipp, 1987) nodes are located at the intersection of grid lines and both specified head and specified flow boundaries are placed on the nodes (Fig. 4.11(c)). (The difference between block-centered and point-centered FD models is explored in Problem P4.1.) In an FE mesh, both specified head and specified flow boundaries are placed on the nodes (Fig. 4.11(d) and (e)). Head-dependent boundaries relate flow to a user-specified boundary head and are located on the node in a point-centered FD grid and in an FE mesh and at the edge of the block in a block-centered FD grid.

Figure 4.11 Boundary representation in FD grids and FE meshes of a two-dimensional areal model.

(a) For purposes of illustration, relatively impermeable rock at the mountain front forms a physical no-flow boundary. Streamlines, defined from a water-table map, form hydraulic no-flow boundaries.

The fully penetrating river is a physical specified head boundary. (b) Block-centered FD grid showing that no-flow boundaries are located at the edges of FD cells and specified heads are located on the nodes. The grid is larger than the problem domain. (c) Point-centered FD grid showing that both no flow and specified head boundaries are placed directly on the nodes. The grid coincides with the problem domain. (d) Triangular FE mesh. Node numbers are shown; element numbers are circled.

Both no flow and specified head boundaries are located directly on the nodes. (e) Quadrilateral FE mesh. Node numbers are shown; element numbers are circled. Both no flow and specified head boundaries are located directly on the nodes.

Figure 4.12 Hydraulic boundaries. Schematic water-table contour maps for a regional problem domain (on the left) bounded by physical features and a local problem domain (on the right) with three hydraulic boundaries taken from the solution of the regional problem; the circled dot represents a pumping well (modified from Townley and Wilson, 1980).

4.3.2 Specified Head Boundaries 

A specified head boundary is implemented by fixing head values at the nodes along the boundary; hence, specified boundary heads do not change in response to hydrologic stresses. In steady-state models, specified boundary heads are invariant but most codes allow the user to input a time series of heads to update boundary values during transient simulations. For example, this is done with the Time-Variant Specified-Head Package (CHD Package) in MODFLOW. Under field conditions, the head at the location of a designated specified head node may in fact change with time or the flux across the boundary may be limited by physical or hydraulic constraints. Then, a specified head boundary may cause erroneous simulated heads. Specified head boundaries are best used to represent large bodies of water (major rivers, lakes, reservoirs, and oceans) that are not affected by stresses in the system such as pumping and changes in recharge rate. Less commonly, a specified head condition can be used to solve for associated discharge in a pumping well. For example, a specified head condition allows the head in the well to be fixed at the desired elevation for problems involving dewatering without the need to know the pumping rate that will produce that head. More commonly, however, pumping wells are simulated as specified flow conditions by specifying the pumping rate (Section 4.3.3 and Section 6.2).

In 2D areal models specified head boundaries necessarily fully penetrate the aquifer, physically (Fig. 4.3) or hydraulically. The specified boundary head represents the vertically averaged head in the aquifer at the boundary and is enforced throughout the full thickness of the aquifer. Surface water bodies in the interior of the model can be simulated using specified head conditions but the modeler should be aware that this assumes the surface water body fully penetrates the aquifer. For example, a stream represented by specified head nodes bisects a 2D areal problem domain into two independent regions each of which could be modeled separately. Similarly, placing a no-flow boundary beneath a topographic divide (Fig. 4.9(a)) effectively divides the problem domain in two. (This issue is further explored in Problems P4.3b and P4.4.) Partially penetrating surface water bodies in the model interior (e.g., Fig. 4.9(b)) are better represented using head-dependent conditions (see below).

4.3.3 Specified Flow Boundaries

A specified flow boundary is implemented by setting the flow at the boundary. Then heads at boundary nodes are calculated by the code and can change as the simulation progresses in time. Most codes allow the modeler to input time-varying boundary flows in transient simulations. Specified flow can be implemented along the lateral boundaries of a 2D or 3D model to represent horizontal flow across the boundary, such as groundwater entering or leaving the problem domain as underflow (e.g., Figs. 4.13(a) and (b), 2.4 and 2.15). A specified flow boundary can also be placed at the top of a 3D model to represent the vertical flow across the water table representing recharge (Fig. 4.13(a); Section 4.5), or at the bottom of the model to represent leakage to or from an underlying unit that is not included in the model. Specified flow conditions are also used to represent point sources and sinks such as pumping and injection wells (Section 6.2).

Conceptually, the same hydraulic effect can be achieved using either specified head or specified flow conditions. However, it is usually preferable to specify flow rather than head. A specified head boundary fixes the head at the boundary so that flow across the boundary is dependent on the head gradient, which is calculated by the model; hence, the calculated flow may not be representative of field conditions. A specified flow boundary, on the other hand, maintains a constant realistic flow of water into or out of the model, as designated by the modeler. Moreover, simulated heads are more sensitive to parameter changes during model calibration when a specified flow boundary is used, which helps with parameter estimation (Chapter 9).

A no-flow condition (i.e., flow is specified to be zero) is the default perimeter boundary condition in both FD and FE codes; the user must activate another option in the code to over-ride the default. Within the computer code itself, the implementation 148 Applied Groundwater Modeling of no flow conditions in FD codes involves ghost nodes (Fig. B4.3.3, Box 4.3; Problem P4.1). For example, MODFLOW internally assigns default zero values of hydraulic conductivity to ghost nodes located outside the boundaries of the model in order to create no-flow conditions at the edge of the boundary cells (Fig. 4.14). In an FE code, the boundary flux is automatically set to zero when the FE equations are assembled into the global matrix equation. Details are given by Wang and Anderson (1982, pp. 117, 126e127). In both FD and FE codes, the user may implement a nonzero specified flow perimeter boundary by introducing flow into the boundary node via the code’s recharge or pumping/injection well boundary conditions (Figs. 4.13(a), 4.14 and 4.15). In that way, water is introduced into the boundary node as a sink or source thereby inducing flow at the boundary.

 

Figure 4.13 Implementation of specified flow conditions. (a) In an FD grid, a volume of water is placed into an FD cell/block (or extracted from the cell/block) using wells (Q) or areal recharge (R). Lateral flows (i.e., underflow) can be introduced using a code’s well or recharge option or using a headdependent boundary. For underflow input via a well, Q= UΔxΔz where U is the lateral flux (L/T).

When underflow (U) is input as recharge, R = UΔz/Δy. (b) In an FE mesh, diffuse flow (Qs) along the boundary is discretized along the sides of triangular elements and then assigned to nodes (modified from Townley and Wilson, 1980).

Flow across the boundary is typically estimated from field data or information taken from the literature. In FD models, flow is input to an FD cell as a flux (L/T) or a volumetric flow rate (L3 /T) (Figs. 4.13(a) and 4.14). A common application of a specified flow boundary is to represent a pumping or injection well (Section 6.2). Conceptually, water enters or leaves the cell through one of its faces and then enters the faces of adjacent active cells. For example, underflow conceptually enters or leaves through a side face; in 3D models recharge enters the top face of the cell and leakage enters or leaves through side and bottom faces. In an FE model, the flux at the boundary is proportionately distributed to each boundary node (Fig. 4.13(b)).

Figure 4.14 Default no-flow boundaries in a two-dimensional block-centered FD grid are implemented by setting transmissivity (T) equal to zero in inactive cells/ghost cells (shaded) outside the problem domain. The ghost cells with T = 0 that are used to implement no flow conditions along the groundwater divides and the ghost cells to the right of the underflow boundary are not shown.

Constant head and specified flow (injection wells) conditions are imposed in boundary cells to cancel the effect of the default no-flow boundaries (adapted from McDonald and Harbaugh, 1988).

Although specified flow boundaries are preferred over specified head boundaries, the exclusive use of specified flow boundaries generally is not recommended even when hydrogeologically defensible. Recall that the governing equation (Eqn (3.12)) is written in terms of derivatives or differences in head. The solution is satisfied when the gradients of head are consistent with the specified boundary flux. However, when a steady-state problem is set up only in terms of gradients, the problem is nonunique because gradients are relative. In that case, many different arrays of head values can produce the same gradients. Nonuniqueness is avoided and a good solution attained by specifying head for at least one perimeter or internal boundary node to give the model a reference elevation from which to calculate heads. A head-dependent boundary condition can also serve this purpose provided the associated resistance is not extremely large (Section 4.3.4).

In a transient simulation, the heads specified in the initial conditions (Section 7.4) provide reference elevations for the solution and a unique solution can be obtained without specified head or head-dependent conditions. However, because errors in the initial heads can cause errors in the solution, the output should be reported as a difference from the initial condition (e.g., drawdown) when all boundary conditions in a transient model are set as specified flow. Hence, even in transient modeling it is usually advisable to specify a head value as an internal boundary in a model that uses all specified flow boundaries (Problem P4.5).

Figure 4.15 FD grid showing implementation of head-dependent boundary (HDB) conditions for boundary flows (designated as general head boundary cells), drains, and evapotranspiration. Stream cells and lake cells use more sophisticated representation of HDB conditions as discussed in Sections 6.5 and 6.6, respectively (Gannett et al, 2012).

4.3.4 Head-dependent Boundaries

When a head-dependent boundary (HDB) is implemented, the code calculates flow across the boundary using the hydraulic gradient between a user-specified boundary head and the model-calculated head at a boundary node.

Mathematically, the volumetric flow rate, Q (L3 /T), across an HDB is computed using an equation of the general form:

Q = C Δ h = C（h_B h_{i,j}^{k}）

C = KA/L

where Δh is the difference between the user-specified boundary head, hB, and the model-calculated head near the boundary, h_{i,j}^{k}; C is conductance (L2 /T), which is computed using a representative hydraulic conductivity (K) times a representative area (A), divided by the distance (L) between the locations of hB and hi,j,k. The hydraulic gradient Δh/L can represent either horizontal or vertical flow. An advantage of an HDB is that it gives the modeler flexibility in choosing the location of the specified boundary head, hB, which does not have to be located directly on the boundary.

Another advantage is that in transient simulations hi,j,k changes with time as the simulation progresses and the simulated boundary flow (Q) also is automatically updated.

The flexibility in specifying both the location of the boundary head and the direction of the associated hydraulic gradient means that HDB conditions can represent a wide variety of hydrogeological situations including vertical flow to and from streams, lakes, wetlands, and other surface water bodies; flow to drains; evapotranspiration from the water table; lateral and bottom boundary flows; and boundaries outside the model domain (Fig. 4.15). MODFLOW, for example, has several options for implementing HDBs including basic options provided by the River (RIV) Package, the Drain (DRN) Package, the Evapotranspiration (EVT) Package, and the General Head Boundary (GHB) Package and more advanced options provided by the Lake (LAK) Package and the Wetlands Package. The way in which conductance, C, is specified varies with the specific application of an HDB, as discussed in the sections below. Under some conditions, an HDB can generate anomalous large water budget errors (Box 4.5).

Box 4.5 Large Water Budget Errors Arising from an HDB 

A large error in the model’s water budget (i.e., greater than 1.0%) is generally indicative of conceptual errors in model design, nonconvergence of the solution, and/or data entry errors (Sections 3.6 and 3.7). However, large water budget errors (up to 200%) may be generated for some formulations of HDB conditions (Sections 4.3; 6.4–6.7) even when the head solution is acceptable (i.e., the head closure criterion is met; Section 3.7). This type of water budget error does not necessarily invalidate the results of the model as long as the head closure criterion is met and flow to or from the HDB is identified as the cause of the water budget error. The way in which the error occurs is explained below.

Large errors in the computed water budget can occur as an artifact of representing the transfer of water between a perimeter or internal HDB and the aquifer. Such anomalous conditions may occur (if they occur at all) when HDB conditions are used to represent distant boundaries (Section 4.4) or large lakes (Sections 4.3 and 6.6). Recall that when implementing an HDB, the modeler assigns a boundary head and a value for conductance (Eqn (4.4)). Anomalously large volumes of water can be transferred when the assigned conductance is set very large and the gradient at the boundary is not equal to zero. The flow across the boundary is calculated using the simulated head gradient at the HDB and the large assigned conductance, which may produce an anomalously large flux of water. Even a small head gradient at the boundary could produce large fluxes of water if conductance is large. The head closure criterion is still met because the head residual at the HDB node is small.

Whether such conditions have a significant impact on the water budget depends on the relative magnitude of the sum of such anomalous computed flows compared to other components of the water budget. The water budget error can often be eliminated by reducing the value of the conductance until the model produces an acceptable water budget error while still meeting the head closure criterion and simulating appropriate responses at the HDB.

4.3.4.1 Surface Water Bodies 

Head-dependent conditions are frequently used in both 2D areal and 3D models to represent exchange of water at surface water bodies including streams, lakes, and wetlands (Fig. 4.15). The surface water body is not represented explicitly in the model (i.e., it does not occupy space within the grid/mesh); rather, an HDB is specified for the node influenced by the surface water feature (Fig. 4.16(a)). Recall from Section 3.3 that vertical flux, qz, across an HDB is calculated as follows:

q_z =-K'_ z \frac{h_{i,j,k} - h_s}}{b'}

where here q_z is vertical flux of water to or from a surface water body, h_{i,j,k} is the groundwater head computed by the model beneath the surface water body; hs is the surface water head specified by the modeler (i.e., the specified boundary head); K'_ z is the vertical hydraulic conductivity at the interface between the groundwater system and the surface water body (usually the vertical hydraulic conductivity of streambed,

lakebed, or wetland sediments) and b0 is the thickness of the interface (usually the thickness of streambed, lakebed, or wetland sediments).

Equation (4.5) is conceptually identical to Eqn (4.1) for leakage through a confining bed. In Eqn (4.5), the connection between the surface water feature and the groundwater system is characterized by leakance, K0 z/b0 , or its inverse vertical resistance, b0 /K0 z, of the streambed, lakebed, or wetland sediments. Then, K/L in Eqn (4.4b) equals K0 z/b0 in Eqn (4.5). To compute conductance, A in Eqn (4.4b) is the horizontal area perpendicular to flow, calculated from the width, W, and length, LR, of the surface water segment or reach (A ¼ WLR) (Fig. 4.16(e)). When the surface water feature covers the entire surface of the cell or element, A is equal to the surface area of the relevant cell or element, but the width of the surface water feature can be smaller than the width of the cell or element (Fig. 4.16(a)) and similarly the length of the surface water feature may not coincide with the length of the cell or element (Fig. 4.16(e)).

HDBs can represent both gaining (flow out of the groundwater system; Fig. 4.16(b)) and losing (flow into the groundwater system, Fig. 4.16(c)) conditions. Some specialized HDBs (e.g., the River Package for MODFLOW) disconnect the surface water feature from the aquifer when the head in the aquifer falls below the bottom of the interface sediments (SBOT in Fig. 4.16(d)). Then, percolating conditions occur and flow out of the surface water body is independent of groundwater head and equal to C(hsSBOT) (Fig. 4.16(d)). When conditions are near but not at percolation, the surface water feature is said to be steeply mounded.

Most codes provide basic HDB conditions to represent vertical flux between the aquifer and surface water bodies. However, such simple HDB conditions, which use fixed surface water levels, only approximately represent the exchange of water between groundwater and surface water bodies. More advanced methods for simulating streams, lakes, and wetlands in a groundwater model are discussed in Sections 6.5, 6.6, and 6.7, respectively.

Figure 4.16 Implementation of HDB conditions for representing surface water bodies using a stream in an FD cell for illustration. (a) Representation of the stream in an FD cell. The stream is conceptualized to be embedded in the cell and to exchange water with the aquifer but the stream does not occupy space within the grid. (Representation in an element of an FE mesh is similar.) As shown, the stage of the stream is lower than the head in the cell and the width of the stream is less than the width of the cell. (b) When the stream is gaining, the head in the aquifer, h_{i,j,k}, is higher than the head in the stream, hs. The elevation of the bottom of the streambed sediments is SBOT; the thickness of the sediments is b0 . QGW is the volumetric rate (L3 /T) of groundwater discharge to the stream. (c) For a losing stream h_{i,j,k} < hs and QGW is the volumetric rate (L3 /T) of induced recharge from the stream to the aquifer. (d) Under percolating conditions, the stream is separated from the aquifer and QGW is constant. (e) Discretization of a stream into 12 reaches. The width, W, of the stream is much less than the grid spacing (Dx); the length of the stream reach, LR, is not equal to the length of the cell (Δy). Each reach can have different values for hs, SBOT, K'z/b' , as well as LR and W (modified from McDonald and Harbaugh, 1988).

4.3.4.2 Drains 

Head-dependent conditions can also simulate drains including closed drains, tunnels, and mines, and springs and seeps (Figs. 4.15 and 4.17). The head in the drain is the specified boundary head, hB, in Eqn (4.4a). In most standard representations of drains, the drain is active only when the head in the aquifer, hi,j,k, is higher than the elevation of the drain, hB. That is, drains only remove water from the problem domain; a drain does not contribute water to the model. Discharge to the drain increases as the head rises above the specified drain elevation. If h_{i,j,k}≤  h_B there is no discharge from the drain cell and Q in Eqn (4.4a) is zero by definition. For springs and seeps, the head in the drain is equal to the elevation of the land surface at the location of the spring or seep. Diffuse seepage (e.g., to wetlands) is simulated by placing drain nodes in the general area where seepage is likely to occur (Fig. 4.15) (e.g., Batelaan and De Smedt, 2004).

It is important to remember that water discharged to drains is removed from the model. This may not be appropriate for simulating seepage to springs and wetlands where water may flow overland and re-enter the groundwater system downgradient. Seepage to wetlands can be simulated using other approaches as discussed in Section 6.7. A drain package for MODFLOW (DRT1 by Banta, 2000) allows water to flow from drain nodes to the groundwater system but only at user designated locations. In general, drains are not appropriate for any surface water features that have losing reaches because water cannot leave a drain node and enter the groundwater system.

Conductance of the drain, C in Eqn (4.4b), is affected by the material surrounding the drain. For closed drains (Fig. 4.17) conductance depends on the size and density of openings in the drain, the presence of chemical precipitates around the drain, and the hydraulic conductivity and thickness of backfill around the drain. Conductance is often estimated during calibration by adjusting conductance values until simulated flows to the drain match measured flows. For mines and tunnels (e.g., Zaidel et al, 2010), drain nodes can represent a seepage face (Fig. 4.17) and conductance is related to the geological material that forms the sidewalls.

Figure 4.17 Examples of drains with associated cross-sectional area used to compute conductance (composite of images modified from Yager, 1987; Fipps et al, 1986 and McDonald and Harbaugh, 1988).

4.3.4.3 Evapotranspiration (ET) from the Water Table 

In many hydrogeologic settings, ET occurs mainly or exclusively from the unsaturated zone rather than from the water table. In that case, ET from the saturated zone is zero and the loss of infiltration via ET in the unsaturated zone is accounted for by inputting a net recharge rate to the groundwater model equal to precipitation minus ET (Box 5.4).

However, if the water table is close to the land surface, there may be direct evaporation from the water table and/or phreatophytes (plants whose roots extend into the water table) may extract groundwater through transpiration. ET ceases when the water table drops to or below the root zone (often called the extinction depth, d in Fig. 4.18). In temperate climates, ET varies seasonally with high rates during the growing season and little or no ET during plant senescence.

In groundwater models, ET is usually expressed as a flux across the water table. In 2D areal models, ET is represented by internal HDBs; in 3D models, the water table typically is the upper boundary and ET is represented by HDB conditions at water-table nodes (Fig. 4.15). In a basic representation of ET, the modeler assigns the extinction depth (d); a boundary head (hs), which is usually equal to the land surface elevation; and a maximum ET rate (RETM). The maximum ET rate occurs when the water table equals or exceeds hs. In between the location of the boundary head (land surface) and the extinction depth, the volumetric rate of water loss via ET varies linearly, where ΔET = C Δh (Eqn (4.4a)) and C and Δh are:

C=\frac{R_{ETM}}A}{d}=\frac{Q_{ETM}}}{d}

\Delta h=h_{i,j,k}-(h_s-d))

Figure 4.18 Representation of ET as a head-dependent boundary (Eqn (4.6)) showing extinction depth (d), land surface elevation (hs), hs-d, and calculated head (h).

In Eqn (4.6a), A is the surface area of the cell (e.g., DxDy) or element; in Eqn (4.6b), hi,j,k is the head at the water table calculated by the model (Fig. 4.18). Here, C is not a conductance but is defined so that there is a linear increase in ET between the extinction depth where ET is zero, and the land surface where ET is at its maximum. Other functions can also be used to represent the relation between ET and Dh (e.g., Banta, 2000).

The field information for estimating the extinction depth is most often obtained from estimates of plant rooting depth. The maximum ET rate (RETM) is commonly estimated from remote sensing and climatological information. As is true for many point measurements made in the field, local factors confound point measurements of RETM (e.g., Lott and Hunt, 2001) and in practice accurate estimation of RETM is complicated. Moreover, scaling up point measurements to represent ET in a cell or element of a groundwater model is a challenging problem. ET estimation techniques and issues are continually being researched. The discussion of these issues is beyond the scope of our text; more detail is provided by Abtew and Melesse (2012), Goyal and Harmsen (2013), and Moene and van Dam (2014), among others.

4.3.4.4 Lateral Boundary Flows and Distant Boundaries

An HDB condition is often used to simulate lateral boundary flows (Fig. 4.15), including underflow (Figs. 4.13(a) and 2.15), and flows to and from a distant boundary located outside of the model domain (Fig. 4.19). In MODFLOW, lateral boundary flows and distant boundaries are represented using the General Head Boundary (GHB) Package.

Figure 4.19 Head-dependent boundary used to represent flow between the modeled area and a distant physical boundary, shown here as a large lake. The boundary flow (Q) is controlled by the head at the distant boundary, shown here as the head in a large lake, hB. C is conductance (Eqn (4.4b)).

For lateral flows at perimeter boundaries, the user-specified boundary head, hB in Eqn (4.4a), is the head at or near the perimeter boundary and conductance (Eqn (4.4b)) reflects conditions at the boundary. An HDB can also represent lateral flows to or from a distant boundary (e.g., the lake in Fig. 4.19). In that case, the HDB effectively extends the model to a distant physical feature without expanding the grid/mesh.

The head at the distant physical boundary (e.g., lake level in Fig. 4.19) is the userspecified boundary head. Conductance (Eqn (4.4b)) is computed with K equal to the average horizontal hydraulic conductivity between the perimeter boundary of the model and the distant physical boundary (e.g., the shoreline of the lake in Fig. 4.19); L is the distance between the model boundary and the distant boundary; A is the area of the face of the cell or element perpendicular to flow. In this way, the perimeter boundary of the model is tied to a physical feature without extending the grid/mesh to its physical location.

Although an HDB is usually more defensible than an inferred hydraulic boundary, professional judgment is required to decide how far a physical boundary can reasonably be located from the perimeter of the model. That is, over what distance can the average hydraulic conductivity between the distant physical feature and the perimeter of the model be reasonably estimated? Some applications of this type of HDB condition have involved large distances. For example, Handman and Kilroy (1997) used a distant HDB to represent springs located 16 miles from the perimeter of the model, where the specified boundary head was set at the elevation of the springs.

4.4 EXTRACTING LOCAL BOUNDARY CONDITIONS FROM A REGIONAL MODEL

Selecting perimeter boundaries is often difficult because physical boundaries and suitable hydraulic boundaries (e.g., stable groundwater divides) may not be located near the problem domain. An HDB can be used to tie the perimeter boundary of the model to a head at a distance physical boundary (Section 4.3) or the problem domain can be expanded so that the boundaries are located on distance physical features. Expanding the problem domain requires expanding the grid/mesh and increasing the number of nodes, which increases the computational burden. In an FE mesh and an unstructured Control Volume Finite Difference (CVFD) grid, a fine nodal network may be embedded within a coarser regional network of nodes (Section 5.1) so that perimeter boundaries are placed at distant physical features while minimizing the number of nodes in the far-field area of coarse nodal spacing.

Perimeter boundaries for local scale models can also be defined by extracting heads and flows from large-scale groundwater flow models (Fig. 4.20). These techniques work well with standard FD and FE models. Two of these methods are telescopic mesh refinement (TMR) (Ward et al, 1987; Leake and Claar, 1999), which applies to both FD and FE models, and local grid refinement (LGR), which was developed for MODFLOW (Mehl and Hill, 2006). In TMR, boundary conditions for FD or FE models covering successively smaller geographic areas are assigned based on the heads and flows computed by a larger scale FD or FE model. Each model is run independently and boundary conditions are extracted successively after each run. Some graphical user interfaces (GUIs) (e.g., Groundwater Vistas) have a TMR option that: extracts heads and flow along a designated boundary location for the smaller scale model; calculates specified heads or flows along the boundaries of the smaller scale model; imports that information into the smaller scale model as boundary conditions and runs the smaller scale model. Alternatively, results from a regional AE model can be used to specify boundary conditions for a local FD or FE TMR model (e.g., Hunt et al, 1998; Feinstein et al, 2003). The GUI for the AE code GFLOW, for example, exports heads or flows from the solution of the AE model directly to an FD or FE model as specified head or specified flow boundaries (Fig. 4.21).

LGR is conceptually similar to TMR except there is an iterative feedback between the regional and local models to update the calculated boundary conditions (e.g., see Feinstein et al, 2010; Hoard, 2010). In practice, TMR is usually preferred over LGR because LGR can be cumbersome to use, can substantially increase runtimes, and may be less efficient than using a single globally refined grid (Vilhelmsen et al, 2012). TMR is appropriate for most problems provided the perimeter of the local scale model is sufficiently far from the area of interest so that results are not dominated by conditions along the boundaries. However, TMR procedures are more complicated and cumbersome when conditions along the extracted boundaries change with time. In transient TMR models, extracted boundary conditions must be updated as the rate and direction of groundwater flow across the boundaries changes with time (e.g., see Buxton and Reilly, 1986). LGR models can automate the transient updating of inset perimeter boundaries.

Figure 4.21 Hydraulic boundary conditions for a 3D FD model (basin model) are extracted from the solution of a 2D regional analytic element (AE) model. Lake and stream analytic elements are outlined in blue and pink. Heads (dashed lines) calculated by the AE model were used to compute fluxes along the perimeter boundaries (outlined in red) of the FD model. Fluxes extracted from the AE model were uniformly distributed vertically along the perimeter of the five layer FD model (modified from Hunt et al, 1998).

4.5 SIMULATING THE WATER TABLE 

In a rigorous formulation of groundwater flow under unconfined conditions, the water table is a moving free surface boundary (Fig. 4.22(a)) whose location depends on the solution of the model. Solving the free surface boundary problem is challenging but some analytical solutions are available for special cases (Bruggeman, 1999). Simulating the water table as a moving boundary in a numerical model (e.g., Diersch, 2014, pp. 416e421) requires a code that allows for movable nodes. More often, in practice, the boundary condition on the water table is simplified to make the solution tractable using a grid/mesh with nodes that are fixed in space. In 2D areal models of unconfined groundwater flow using the D-F approximation (Section 4.1; Box 4.1), head at the water table is the dependent variable calculated as the solution (Fig. 4.3). In 2D profile models and 3D models, however, the water table is usually the upper boundary of the model. The modeler can specify the head at the water table (Box 4.3) or, more commonly, flux across the water table (Box 4.6). Generally speaking, it is preferable to specify flux (recharge) across the water table rather than head (Box 4.6).

In either case, at least the approximate elevation of the water table is needed in order to locate the boundary correctly. In this section, we present some options for simulating the water table in numerical models, but first we review a few basic principles.

Groundwater flows in response to gradients in total head (h). Total head (or head) is the sum of pressure head (Fig. 4.2) and elevation head (z). At the water table, pressure is equal to atmospheric (zero) pressure so that head at the water table equals elevation head. The fundamental dilemma in problems with a water-table boundary is that we need to know the head at the water table to set the location of the boundary when the water table head is what we would like the model to calculate. Another complication is the associated seepage face (Figs. 4.3, 4.22(b) and (c)), whose position is also unknown. Similar to conditions along the water table, the pressure along the seepage face is zero (atmospheric), so that the head at the seepage face equals elevation of the seepage 162 Applied Groundwater Modeling 

face. The removal of water through a seepage face can be simulated using an HDB with drain nodes (Section 4.3; Fig. 4.17) and there are also options for simulating a seepage face along the well casing of a pumping well (Section 6.2). However, standard groundwater flow codes typically do not solve for the location of a seepage face, which is usually not a concern because for most groundwater problems the seepage face is a small feature relative to the size of the flow system. Solving for the location of the seepage face, however, can be important to engineering problems involving slope and tunnel stability, dewatering, and flow through a dam.

Below we discuss three options for simulating the water table using: (1) standard FD and FE codes with fixed nodes; (2) FE codes with movable nodes; (3) variably saturated flow codes, which simulate the unsaturated and saturated zones as a continuum. The first method is the most widely used in practice; however, if finding the location of the seepage face is a primary modeling objective, the water table and seepage face should be simulated using options 2 or 3.

Figure 4.22 Representation of hydraulic conditions at the water table and seepage face. (a) The water table is a streamline when there is no recharge (left hand side figure) but is not a streamline when recharge is present (right hand side figure). In both cases, the pressure head at the water table is zero. (b) A seepage face (DC) along a streambank (left hand side figure). Schematic flowlines and arrows are shown. The location of the water table (DE) and the point of intersection of the water table with the streambank (D) are unknown. Right hand side figure shows detailed schematic depiction of flow near the seepage face. The pressure head at the seepage face is zero so that head at the seepage face is equal to elevation head(modified from Fitts, 2013). (c) The water table computed as the surface of zero pressure in a variably saturated model. The aquifer is shown in cross section with vertical exaggeration ¼ 10. Equipotential lines are computed in the unsaturatedesaturated continuum and are closely spaced near the discharge face at the ocean (shaded in green). The ocean level and seepage face are also shown (modified from Ataie-Ashtiani, 2001).

4.5.1 Fixed Nodes

In standard groundwater flow codes nodes are fixed in space, which means that the calculated head at the water table, (hi,j,k)WT, is not necessarily equal to the elevation of the water table, zi,j,k. Instead, (hi,j,k)WT, is expected to fall within the range (zi,j,k  Δz/2) where Δz is the nodal spacing around zi,j,k (Fig. 4.23). If the calculated head is higher than the top of the water-table cell (usually set at the elevation of the land surface), the water-table cell converts to confined conditions and the cell is said to be flooded (Section 5.3). If the calculated head is lower than the bottom of the water table cell, the cell is dewatered (see below).

Usually water-table nodes are assigned specified flow boundary conditions (Box 4.6) and then the code solves for head in the water-table cells. When (hi,j,k)WT is not equal to the elevation of the node (zi,j,k), which is usually the case, the water-table boundary is not represented exactly. Clemo (2005) presented a modification to MODFLOW-2000 to adjust the placement of the node within a cell during the simulation, but for most applications the error in the unadjusted water-table head is acceptable without such adjustments.

In the most straightforward case of using fixed nodes, the water table remains within an unconfined layer (see Section 5.3). Spatially, the water table might also be present in parts of several layers, each of which is the upper active layer in the model (e.g., Fig. 4.6). If the water table drops below the bottom of a water-table layer, nodes in the layer are de-watered. Modelers struggled with dry nodes for over two decades.

In one simple approach used in older codes, the dry nodes were removed from the array of active nodes for which head is calculated. Procedures for removing dry nodes and adding them back when rewet were developed but were inherently unstable (Doherty, 2001; Painter et al, 2008) and led to incorrect solutions. A more successful approach, which is appropriate for problems where the saturated thickness of the water-table layer is not expected to change appreciably, is to simulate unconfined water-table layers as confined layers (Section 5.3). This approach linearizes the problem and produces acceptable solutions for many problems (Sheets et al, 2015; Juckem et al, 2006).

The best solution is realized in modern versions of MODFLOW (e.g., MODFLOW-NWT Niswonger et al, 2011) where the problem of dry nodes is resolved more effectively and robustly by using an improved method for formulating and solving the FD equations (see review by Hunt and Feinstein, 2012; test case presented by Bedeker et al, 2012).

Figure 4.23 Water table in a three-dimensional FD grid showing that head calculated at the watertable node (h) is higher than the bottom elevation of the top layer of the model but is not necessarily equal to the elevation of the node. (A similar situation occurs in a fixed node FE mesh.)

4.5.2 Movable Nodes 

Movable nodes are more easily incorporated into FE codes because the FE method references the location of the node within 3D space to construct even a fixed node mesh (Section 3.5; Fig. 3.5). However, additional programming (not included in all FE codes) is required to update the location of a movable node as the solution progresses. The FE code FEFLOW (Diersch, 2014) and some special purpose FE codes (e.g., FREESURF by Neuman, 1976, Neuman and Witherspoon (1971); AQUIFEM-N by Townley, 1990) allow for movable nodes that track the water table (Fig. 4.24(a)). Movable nodes are placed exactly on the water table and the head at the node is exactly equal to the elevation of the node (i.e. (hi,j,k)WT=zi,j,k). Movable nodes also make it possible to calculate the position of the seepage face (Fig. 4.24(b)). Of course when nodes move, the shape of the affected elements change (deform); codes that include movable nodes also account for deformable elements.

Figure 4.24 Movable nodes and deformable elements (shaded) in FE meshes. (a) Movable nodes are placed along the water-table boundary (modified from Mitten et al, 1988). (b) Movable nodes are placed along the water-table boundary and along the exit face and in the interior of a permeable earthen dam. The model solves for the location of the water table and associated seepage face; nodes 25 and 30 are on the seepage face (modified from Neuman, 1976).

4.5.3 Variably Saturated Codes 

The most realistic and theoretically rigorous way to simulate the water table uses a variably saturated code that represents the unsaturated and saturated zones as a continuum (Box 6.2; Section 12.2). The infiltration rate is specified across the upper boundary of the model, which is typically the land surface. The code solves for water pressure or pressure head in the subsurface continuum and the water table is determined as the surface of zero (atmospheric) pressure (Fig. 4.22(c)). The location of the seepage face is determined iteratively (Cooley, 1983; Neuman, 1973). Of course, a model that includes the unsaturated zone is more complicated than one representing only the saturated zone, requires a more complicated governing equation, and has much longer runtimes than the other two approaches.

Box 4.6 What Controls the Water Table? 

The water table is an important boundary in many groundwater models but the water-table configuration is usually poorly known. It is a common misconception that the water table is always a subdued replica of the land surface, i.e., that the water table is controlled by topography. In many hydrogeologic settings the water table is controlled by recharge rather than topography (Haitjema and Mitchell-Bruker, 2005; Shahbazi et al, 1968). The distinction is important because a water table controlled by topography is often represented by a specified head boundary condition whereas a water table controlled by recharge is usually represented by a specified flow boundary.

Groundwater modelers are tempted to specify heads at the water table by extrapolating field measurements of the water table to mimic the land surface. However, even small extrapolation errors in defining the water-table surface (for a 3D model) or the water-table profile (for a 2D profile model) can cause significant errors in the flow of water supplied to the model as recharge and removed as discharge (e.g., Stoertz and Bradbury, 1989). Recall from Section 4.3, that a specified head boundary potentially supplies infinite amounts of water to the problem domain and potentially accepts infinite amounts from the problem domain. Therefore, an incorrectly specified water-table surface potentially could provide unconstrained flows into and out of the specified boundary heads. Such flows are unlikely to represent actual flows under field conditions. Alternatively, the modeler may specify recharge as flux across the water table so that flow to or from the water table is constrained to reasonable values. Although recharge rates are difficult to measure (Section 5.4), approximate values of recharge can be estimated with some degree of confidence (Box 5.4). Moreover, initial values of recharge rates input to a model almost always are adjusted and improved during calibration.

Specified head water-table boundaries are still occasionally used to simulate topography driven groundwater flow (e.g., Zlotnik et al, 2011, 2015), which is also called Tothian flow in reference to groundbreaking work by Toth (1962, 1963) . Toth ’s work was pioneering in that it shifted the focus of hydrogeology from well hydraulics to flow system analysis. He developed a profile model for regional groundwater flow with specified heads along a linear water table (Toth, 1962 ; Box 4.3); later he used a sinusoidal function (Toth, 1963 ) to specify heads along the water table. The latter model produced a distinctive pattern of nested local, intermediate and regional flow cells (Fig. B4.6.1). Diagrams like Fig. B4.6.1 became icons for regional groundwater flow (Toth, 2005 ).

Nested flow cells, however, are not present when the water table is controlled by recharge (Fig. B4.6.2; also see Haitjema and Mitchell-Bruker, 2005). Moreover, most practical problems in applied groundwater modeling involve recharge driven flow as can be demonstrated using a dimensionless parameter, h_D/d

\frac{h_D}{d}=\frac{R}{k} \frac{L^2}{8bd}

developed by Haitjema and Mitchell-Bruker (2005) from the parameter grouping in Eqn (B3.2.4) (Box 3.2). Equation (B4.6.1) is used to help determine conditions under which the water table is controlled by topography versus recharge. In Eqn (B4.6.1), h_D is the head at the groundwater divide for 1D flow in an unconfined aquifer bounded by two streams separated by a distance L; d represents the terrain rise above a horizontal datum (Fig. B4.6.3); R is recharge rate; K is hydraulic conductivity; and b is thickness as shown in Fig. B4.6.3. When hD is approximately equal to d (i.e., hD/d ≈1) the water table is not completely controlled by either topography or recharge. When h_D/d > 1, the water table intersects the land surface and is controlled by topography. Under such topography-driven flow, the water table is controlled by recharge under topographic highs and by discharge to topographic lows. When h_D/d < 1, the water table is controlled by recharge.

Figure B4.6.1 Schematic diagram of a regional flow system when the water table is controlled by topography, based on Toth ’s (1963) profile model. A sinusoidal specified head condition, intended to mimic topography, was imposed at the water table. The model simulates nested local, intermediate, and regional flow cells (Winter et al, 1998).

Figure B4.6.2 Water table controlled by recharge in a 2D profile model. The water table was specified using heads determined from a Hele-Shaw analog model (Section 1.2) in which uniform recharge was infiltrated at a sinusoidal land surface. The water table does not follow the sinusoidal function of the land surface and nested flow cells are not present (modified from Shahbazi et al, 1968).

Figure B4.6.3 Conceptual model of one-dimensional flow under the D-F approximation in an unconfined aquifer under uniform recharge, R. The maximum terrain rise, d, is the largest vertical distance between the datum (defined by the heads at the boundaries) and the land surface. The vertical scale is greatly exaggerated for purposes of illustration.

For Toth  ’s (1963) profile model h_D/d = 40 (calculated from Eqn (B4.6.1) with L/b= 4, L/d = 400, R/K = 0.2; Haitjema and Mitchell-Bruker, 2005), which means that, as expected, the water table is controlled by topography. In the Toth model, where R/K = 0.2, low hydraulic conductivity material is subjected to high recharge. In hydrogeologic settings with more permeable hydrostratigraphic units, hD/d is usually less than 1 and the groundwater system is recharge controlled. This follows from the observation that R/K is commonly in the range 106 to 103 in most field settings. For example, h_D/d is much less than 1 (=0.0088) when R/K= 1.4*10^-6 (e.g., K =50 m/d; R= 25.6 mm/year = 10 inches/year); L/b = 100 and L/d= 500 (Haitjema and Mitchell-Bruker, 2005). The analysis of Haitjema and MitchellBruker (2005) suggests that topography driven flow as shown in Fig. B4.6.1 is not commonly encountered in groundwater problems focused on aquifers rather than aquitards.

In most groundwater flow problems, water tables are better represented by specified flow boundary conditions rather than specified heads. Specified flow conditions allow the model to calculate the heads at the water table, are more likely to approximate the true flow to and from the water table, and are conceptually more appropriate for most aquifer systems. Topographic control of the water table occurs when a groundwater system is characterized by low hydraulic conductivity and high recharge, high vertical anisotropy, and/or deep flow systems where the ratio of system length, L, to saturated thickness, b, is small (e.g., L/b =4 in Toth ’s problem; also see Lemieux et al, 2008; Haitjema and Mitchell-Bruker, 2005). However, streambed morphology driven flow, a concept similar to topographically driven flow, can be effectively used to simulate groundwater flow in a hyporheic zone (the region below a streambed where surface water mixes with groundwater). In models of hyporheic flow, heads are specified along the top of the streambed, which forms the upper boundary of the groundwater model (e.g., Zlotnik et al, 2011).

4.6 COMMON MODELING ERRORS

• The modeler uses specified head and head dependent boundaries (HDBs) to represent surface water features but neglects to check whether the amount of water exchanged with the groundwater system is reasonable and consistent with the conceptual model.

Specified head boundaries and HDBs (Box 4.5) may transfer unrealistic amounts of water into or out of the problem domain.

• A drain is used to simulate a surface water feature that has both gaining and losing reaches. Basic representation of a drain using an HDB does not allow water to enter the groundwater system and therefore the drain HDB cannot simulate loss of water from the drain.

• The modeler uses a 2D profile model to simulate pumping wells. A profile model assumes that there is no flow through the sides of the profile and therefore cannot simulate radial flow to a well. Pumping wells must be simulated using an axisymmetric profile model, a 2D areal model, or 3D model.

• A 2D profile model is not aligned along a groundwater flowpath. Profile models simulate flow only within the thickness of the profile and must be aligned with groundwater flow.

• The modeler selects equipotential lines to define hydraulic boundary conditions for a model designed to determine the long term impacts of pumping. Under field conditions, pumping may affect heads at the locations of the selected equipotential lines used to specify boundary conditions thereby invalidating the model’s boundary conditions. Furthermore, specified head conditions based on equipotential lines provide the model with an unlimited supply of water and thereby may underestimate the impact of pumping by incorrectly keeping simulated drawdowns low.

• A model to determine the long-term impacts of pumping uses hydraulic boundary conditions defined by streamlines to set lateral no-flow boundaries. In the field, the effects of pumping may reach the hydraulic no-flow boundaries and inappropriately affect the expansion of the cone of depression, thereby causing simulated drawdowns that are too large.

• The water table is simulated using specified heads (Box 4.6). The model simulates unrealistic nested flow systems owing to unrealistic flows into and out of the water-table nodes.

4.7 PROBLEMS

Chapter 4 problems introduce boundary concepts and examine the effect of boundary assignment on modeling results. Some of the problems require a spreadsheet and/or a groundwater flow code (either FD or FE); some suggested codes are listed on the Web site for our book (http://appliedgwmodeling.elsevier.com/). Consult the code’s 

user’s manual for specific instructions on implementing boundary conditions. Differences in boundary implementation in FD and FE approaches may cause small discrepancies in the solutions. In this chapter, we introduced pumping and injection wells to simulate boundary flows (Figs. 4.13(a), 4.14 and 4.15) and also mentioned that an internal specified flow boundary condition can represent a pumping or injection well (Section 4.3) as in Problem P4.4. Section 6.2 provides more information on wells as internal sinks.

P4.1 Representation of no-flow and specified head boundaries in an FD grid is illustrated using a spreadsheet to simulate 2D steady-state groundwater flow in the isotropic and homogeneous aquifer shown in Figs. B4.3.1 and B4.3.2 in Box 4.3. We will solve different versions of this problem using block-centered and point-centered FD grids. The top boundary of the model is the water table, which is specified using the equation given in Fig. B4.3.2 (Box 4.3). No-flow boundaries are implemented using Eqn (B4.3.2) in Box 4.3 by incorporating ghost nodes as explained in Box 4.3. For example, along the left hand side no-flow boundary of the block-centered grid in Fig. B4.3.3 (Box 4.3) the head at ghost node A5 is set equal to the calculated head at B5; the no-flow boundary is located half-way between nodes A5 and B5. The FD equations for the block-centered grid, as implemented in the spreadsheet, are shown in Fig. B4.3.3(b) (Box 4.3).

The five-point star operator (Section 3.5), shown in Fig. P4.1(a), is modified for nodes along the left hand side boundary (Fig. P4.1(b), (c)).

a. Use a spreadsheet model to solve the FD equations for heads in a blockcentered grid as given in Fig. B4.3.3(a), (b) (Box 4.3). To solve the spreadsheet in Excel, go to: tools>options>calculations and check manual; also check the iteration box. Press F9 to begin the calculation. Use at least 1000 iterations and an error tolerance of 1E-4m. Also include a water budget calculation. Check your solution against the solution in Fig. B4.3.3(c) (Box 4.3). Draw equipotential lines (using a contour interval of 1 m) and construct flowlines on a properly scaled figure.

Figure P4.1 (a) The five-point star computational module for an interior node (filled circle) in a twodimensional FD grid. The numbers refer to the weighting of heads in the FD equation (Eqn (B4.3.2) in Box 4.3). (b) The computational module for a boundary node (shaded) in a block-centered FD grid; the no-flow boundary is to the left of the node. The head at the ghost node at i-1,j (not shown) equals the head at i,j. (c) The computational module for a boundary node (shaded) in a point-centered FD grid; the no-flow boundary is directly on the node. The head at the ghost node at i-1,j (not shown) equals the head at i+1,j. (d) Point-centered FD grid for the profile model in Box 4.3.

b. Use a point-centered grid with 11 columns and 6 rows (Fig. P4.1(d)). Because no-flow boundaries are placed differently in a point-centered grid, the geometry of the problem domain is different (Fig. P4.1(d)) from the block-centered grid (Fig. B4.3.3(a) in Box 4.3) and the FD equations along the boundaries are also different. (Hint: Write the FD equations at the left hand side no-flow boundary of the point-centered grid using the computational module shown in Fig. P4.1(c). Construct analogous computational modules for the right hand side and bottom boundaries.) c. Solve the problem again using a standard FD code (e.g., MODFLOW) and an FE code (e.g., FEFLOW). Since this is the first problem requiring the use of an FD or FE code, you should test the assignment of values for solver closure criteria and assess solution convergence prior to reporting results (Section 3.7). Compare your solutions with those in parts a and b.

P4.2 In this problem, we will construct a profile model to compute groundwater flow under a dam and solve the model using either an FD or FE code. An impermeable concrete dam 60 m long and 40 m wide is constructed over an isotropic and homogeneous silty sand (K ¼ 10 m/d) that is 26.5 m thick and fills a river valley underlain by impermeable bedrock (Fig. P4.2). The depth of the reservoir pool is 30 m and the tail water stream elevation below the dam is 10 m. The crosssectional area of interest is 200 m long centered on the dam. As you formulate your numerical model think carefully about how to represent boundary conditions. Test values for solver closure criteria and assess solution convergence prior to reporting results (Section 3.7).

a. Whenever possible, simulated results should be checked against some independent calculation. In this case, you can solve the problem graphically using a flow net (Boxes 5.2 and 8.2). Take the dimensions from the Figure and produce a scaled drawing of the profile without vertical exaggeration. Create a 

flow net with curvilinear squares (see Fig. B5.2.1(a) in Box 5.2) and compute groundwater seepage per unit width under the dam (see Box 8.2).

b. Construct a numerical model; assign side boundaries of the profile model at x ¼ 0 and x ¼ 200 m. Use the model to calculate heads and the groundwater flow rate per unit width of the dam (from the calculated water budget).

c. Compare and contrast the graphical solution in part (a) with the numerical solution in part (b) for both heads and flow. If the solutions are significantly different, there is an error in one or both solutions. Correct the error(s) and recomputed heads and seepage per unit width of the dam.

d. Suppose the porous material below the dam is anisotropic so that Kx ¼ 10Kz.

Run the numerical model again but this time include anisotropy. Compute groundwater seepage per unit width of the dam. Why is the seepage rate different from the solution for isotropic and homogeneous conditions? (You can also construct a graphical solution for this problem using a transformed section. See Box 5.2.) P4.3 The town of Hubbertville is planning to expand its water supply by constructing a pumping well in an unconfined gravel aquifer (Fig. P4.3). The well is designed to pump constantly at a rate of 20,000 m3 /day. Well construction was halted by the State Fish and Game Service who manage the Green Swamp Conservation area.

The agency claimed that pumping would “significantly reduce” groundwater discharge to the swamp and damage waterfowl habitat. The town claimed the fully penetrating river boundary to the north and the groundwater divide located near the center of the valley would prevent any change in flow to the swamp.

a. Construct a 2D areal steady-state model of the aquifer between the river and swamp for conditions prior to pumping using the information in Fig. P4.3.

Represent the river and swamp boundaries as constant head boundaries with head set at 1000 m. The side boundaries are no-flow boundaries. Justify this assignment of boundary conditions. Use a constant nodal spacing of 500 m.

Run the model and produce a contour map of heads. Draw the water-table profile in a north-south cross section and label the simulated groundwater divide between the river and the swamp. Compute the discharge to Green Swamp.

b. Using the steady-state heads derived in (a), locate the groundwater divide in the central portion of the valley. Run the model using first a no-flow boundary and then a specified head boundary at the location of the groundwater divide.

Compare results with those in part (a). Compute the discharge to Green Swamp under each representation. What is the effect of assigning an internal boundary? c. Run the model in part (a) again but this time use a HDB to represent the river.

The stage of the river is 1000 m and the width is 500 m. The vertical hydraulic 

conductivity of the riverbed sediments is 5 m/day and the thickness of the sediments is 1 m. The elevation of the bottom of the sediments is 995 m. Compare results with those in part (a).

d. Run the model in part (c) again but this time assume the width of the river is 5 m. What is the effect of reducing the width of the river? P4.4 Return to the model you designed for Problem P4.3 and place a well at the location indicated in Fig. P4.3. Pump the well at a constant rate of 20,000 m3 /day.

Run the model under steady-state pumping conditions three times under each of the following three representations of model boundaries: (1) physical boundaries shown in Fig. P4.3; (2) an internal no-flow boundary at the groundwater divide between the river and the swamp; the location of the divide is determined from the solution of Problem P4.3; (3) an internal specified head boundary at the groundwater divide between the river and the swamp; the location of the divide is determined from the solution of Problem P4.3.

a. Discuss the effects of the internal boundary conditions imposed on the groundwater divide on the resulting head distributions. Compare northesouth watertable profiles drawn through the node representing the pumping well for the three pumping simulations. Compute the discharge to the Green Swamp under each set of boundary conditions.

b. Compare groundwater discharge to the swamp under the prepumping scenario in Problem P4.3 with the results under the pumping scenarios. In light of the modeling results, consider what might be meant by “significantly reduced” as used by the state agency (see discussion in Problem P4.3). Make a list of physical, geochemical and ecological conditions that potentially could be affected by a change in groundwater flow to the Green Swamp.

P4.5 Replace the constant head boundaries at both the river and swamp in Problem P4.3a with specified flow boundaries. Use the water balance results from Problem P4.3a to calculate the boundary fluxes at the river and swamp. Note that all the boundaries of the model are now specified flow boundaries, including the zero flow lateral boundary conditions.

a. Run the model with specified flow boundaries using starting heads of 1000 m and then a second time with starting heads of 2000 m. (Note: Some GUIs will warn that using all specified flow boundaries can create errors, or will not permit the model to execute under these conditions.) Compare the results with those in Problem P4.3a and explain the differences.

b. Take one of the models you designed in part (a) and replace one constant flux node on either the river or swamp boundary with a specified head node equal to 1000 m. Run the model to steady state. Compare the results with those in part (a) and with results from Problem P4.3a. Explain the differences.

P4.6 It is often tempting to define boundaries of an areally extensive aquifer using hydraulic no-flow boundary conditions placed at some distance from the area of interest.

The two no-flow boundaries (streamlines) define a flow tube (or under special conditions, a streamtube; see Fig. B8.2.1 in Box 8.2). One way to determine if this boundary assignment will adversely affect the modeling results is to replace the no-flow boundaries with specified heads. If flow occurs to or from these specified head nodes, the assignment of no-flow boundaries is not appropriate. This method can be illustrated using the Green Swamp problem. Let us assume that the Slate 

Mountains do not exist. Instead the no-flow boundaries define a 4500 m wide flow tube within a larger regional flow system.

a. Replace the no-flow boundaries that delineate the flow tube in Fig. P4.3(a) with specified heads taken from the solution of Problem P4.3a. Run the model to check that it reproduces the steady-state heads of Problem P4.3a. Calculate the eastewest fluxes at the flow tube boundaries and the discharge to the Green Swamp. Be sure your modeled area is 4500 m wide.

b. Simulate the steady-state flow field with the pumping well using the specified head boundaries that define the flow tube. Examine flow to and from the specified head boundaries and compare them to the prepumping rates calculated in part (a). What is the discharge to Green Swamp? Are the hydraulic boundaries that define the flow tube sufficiently removed from the pumping well so that they act as no-flow boundaries? In other words is the flow to and from the specified head nodes along the flow tube boundaries sufficiently small so that there is effectively no flow from the side boundaries into the model domain under pumping conditions?