2. 地下水流基本方程与MODFLOW离散原理

要理解MODFLOW为什么能算出地下水位，必须先把三件事串起来：水在孔隙介质中怎样流动（达西定律）、水在流动中如何守恒（连续性方程）、以及这三者如何被写成计算机能求解的代数方程（有限差分离散）。本节按这条逻辑展开，并给出一个可手算核对的离散示例。

一、达西定律：从水头差到流量

达西定律描述的是层流条件下通过多孔介质的水流通量。以三维形式表达，x、y、z三个方向的比流量（单位面积流量）为：

qx = -Kx · ∂h/∂x
qy = -Ky · ∂h/∂y
qz = -Kz · ∂h/∂z

其中 h 为水头（m），K 为渗透系数（m/d），负号表示水流方向与水头降低方向一致，即从高水头流向低水头。这里必须明确符号约定：本文及后续图示统一采用“水头 h 向上为正、z 轴向上为正”的约定，因此 z 方向重力项已包含在水头定义中，不再单独出现。若改用 z 轴向下为正的约定，则 qz 表达式需相应调整，全文必须保持一致，否则边界条件与源汇项的正负号会出错。

达西定律只回答“流多少”，不回答“水从哪来、到哪去”。要保证每个计算单元里的水不凭空产生也不凭空消失，就需要连续性方程。

二、连续性方程与三维水流控制方程

取一个微小控制体，进入的水量减去流出的水量，加上源汇项，应等于该控制体内储水量的变化率。对不可压缩介质与常密度水，连续性方程可写为：

∂qx/∂x + ∂qy/∂y + ∂qz/∂z = Ss·∂h/∂t - W

其中 Ss 为储水率（1/m），W 为单位体积源汇项（1/d，抽水为负、注水或补给为正）。把达西定律代入连续性方程，得到三维地下水流控制方程：

∂/∂x(Kx·∂h/∂x) + ∂/∂y(Ky·∂h/∂y) + ∂/∂z(Kz·∂h/∂z)
    = Ss·∂h/∂t - W

稳态条件下 ∂h/∂t = 0，方程退化为只含空间项的椭圆型方程；非稳态则保留时间导数项，属于抛物型方程。MODFLOW 求解的正是这个方程在离散网格上的近似形式。

三、有限差分在MODFLOW中的实现逻辑

MODFLOW 采用单元中心有限差分法：把研究区剖分为若干层、行、列组成的长方体单元，水头定义在单元中心（节点），渗透系数按单元赋值。核心做法是用相邻节点水头差近似水头梯度，用单元间距近似导数。

以 x 方向为例，节点 (i,j,k) 与其右侧节点 (i,j+1,k) 之间的流量可写成：

Q(i,j+1/2,k) = CR(i,j,k) · [h(i,j+1,k) - h(i,j,k)]
CR = Kx · Δy · Δz / Δx

CR 称为水力传导系数（conductance），它把渗透系数、单元尺寸和网格间距合并成一个系数。对每个单元写出六个方向的流量之和等于储水变化加源汇项，就得到一条线性方程。所有单元方程联立，形成大型稀疏线性方程组：

A · h = b

MODFLOW 用迭代求解器（如 PCG、SOR）求解该方程组，得到每个单元的水头。边界条件通过修改方程实现：定水头边界把对应节点水头固定；隔水边界令该方向传导系数为零；河流、排水沟等则作为与单元水头相关的源汇项进入右端项 b。

四、一个可手算核对的离散示例

设一维稳态承压含水层，Δx = 100 m，Δy = Δz = 10 m，Kx = 5 m/d，左端定水头 h₁ = 30 m，右端定水头 h₂ = 20 m，中间仅一个活动单元。则：

CR = 5 × 10 × 10 / 100 = 5 (m²/d)
左边界流入：5 × (30 - h) 
右边界流出：5 × (h - 20)
稳态：5(30 - h) = 5(h - 20)  →  h = 25 m

手算结果 h = 25 m 与线性水头分布一致，可作为模型建立后的第一个校核点。若模型输出偏离该值，应优先检查传导系数计算、边界水头赋值和单位一致性，而不是急于调整参数。

五、离散实现中必须注意的几个细节

网格间距与传导系数方向性。CR 中分子是垂直于流向的截面积（Δy·Δz），分母是沿流向的间距（Δx），三个方向的传导系数不能混用。层状介质中 Kz 常远小于 Kx、Ky，若垂向网格过粗，会人为放大或缩小垂向交换量。

储水率与时间离散。非稳态问题中，∂h/∂t 用前后时间步的水头差除以时间步长近似。MODFLOW 默认采用向后差分（隐式），对时间步长较宽容，但时间步过大仍会抹平水位快速变化过程，建议在抽水初期使用较小步长并逐步放大。

边界条件的方程化处理。定水头边界不是“删掉方程”，而是把该节点方程替换为 h = 给定值，使相邻单元通过传导系数获得流量；隔水边界则令对应方向传导系数为零。理解这一点，才能在结果异常时判断是边界设置问题还是参数问题。