时步差分方程求解带来的误差不会在以后的时步中积累，则差分模型是稳定的。否则模型不稳定，随着时步的推进，模型误差将不断积累增大，越来越偏离精确解。
对地下水流的有限差分模型，人们已经证明（陈崇希等，1989）隐式差分格式是无条件收敛和稳定的，而显式差分格式虽然满足收敛性，但一般只有在时间步长较小、网格剖分尺度相对较大的时候才满足稳定性。因此，绝大多数地下水流有限差分模型采用隐式差分格式进行求解。
二、基于矩形网格的地下水流差分方程
用矩形网格对模型空间进行离散，是地下水流数值模拟的常用方法，因为规则网格便于图形显示和数据处理。在矩形网格上建立二维地下水流的有限差分模型，既可以直接套用偏导数的差分格式，也可以根据水均衡原理进行推导。
考虑给一个承压含水层建立矩形网格的差分模型，在矩形网格的每个单元形心处设置一个节点，如图4.2中的圆点，节点的离散网格编号为（i，j），其坐标为（xi，yj）。对于模型边界以内的节点，一般都有8个相邻节点，如图4．2中节点O（i，j）分别与节点1～8相邻。如果含水层的渗透性主轴方向与网格坐标轴的方向一致，则仅用节点l～4就可以建立节点O所在矩形单元的地下水流差分方程。
yj+1547i-1162yjjj-13xii8i+1ox 图4.2 矩形网格差分模型的局部特征
设含水层均质各向异性，导水系数为Tx和Ty，储水系数为S。在t=tk时刻，地下水的水头在节点上的分布为{Hk}，根据达西定律，从周围四个节点（节点l～4）流向节点O的流量用差分格式可分别表示为
kH1k?H0kH2?H0kQ1?Ty?xi；Q2?Tx?yj；yj?1?yjxi?1?xiQ3?Ty?xiH?HH?H；Q4?Tx?yjyj?yj?1xi?xi?1k3k0k4k0
而第k个时步内，节点O所在矩形单元内的地下水储存量增加值表示为
kk?1?Vw?S?xi?yi(H0?H0) （4.2.23）
这个方程中隐含的假定是：单元内的平均水头变化与节点0处的水头变化相等。根据水均衡原理，单元内水体积的增加量应与侧向流量在时间步长内的积累和其他外部补给量相互平衡，即
?Vw?(?Qn?Qa)?tk （4.2.24）
n?14式中：Qa为其他外部补给流量。上式的展开式为
kkkk(H2?H0k)Tx(H3?H0)Ty(H4?H0k)TxQaSkk?1?????(H0?H0)(yj?1?yj)?yj(xi?1?xi)?xi(yj?yj?1)?yj(xi?xi?1)?xi?xi?yj?tkk(H1k?H0)Ty（4.2.25）
这就是对应节点0（i，j）的地下水流差分方程，它可以进一步变形为
kkkkC0H0?C1H1k?C2H2?C3H3?C4H4?F （4.2.26）
其中：
C1??C3??Ty(yj?1?yj)?yjTy(yj?yj?1)?yj；C2??；C4??Tx(xi?1?xi)?xiTx(xi?xi?1)?xi
k?1QaSH0SC0??(C1?C2?C3?C4)；F???tk?xi?yj?tk对于任何一个内部节点都可以写出形如式（4.2.25）的差分方程，然后运用边界条件写
出边界节点的差分方程。所有节点的差分方程联立求解，即得到tk时刻水头分布的离散形式。
如果含水层是非均质的，那么可以对每个节点所在的矩形单元设置不同的参数（导水系数T与储水系数S），例如，用
T0x，T0y；T1x，T1y；???；T8x，T8y；分别表示节点O及相邻节点所在矩形单元的x方向和y方向的导水系数。那么节点之间的单位界面流量，可以用一个平均意义上的导水系数进行计算，以节点1流向节点O的流量为例，即
kH1k?H0 （4.2.27） Q1?T?xiyj?1?yjy0，1其中：平均导水系数通常采用调和平均值
Ty0，1y2T0yT1 （4.2.28） ?yyT0+T1依此类推其他节点与节点O之间的平均导水系数。考虑这种非均质参数分布之后，式
（4.2.26）中的各项系数更新为
yT0,2T0，1C1??；C2??(yj?1?yj)?yj(xi?1?xi)?xixyT0,4T0，3C3??；C4??(yj?yj?1)?yj(xi?xi?1)?xix
QaS0H0k?1SC0??(C1?C2?C3?C4)；F???tk?xi?yj?tk注意：其中的外部补给流量Qa可能既与节点有关，又随时间变化。
如果研究对象是潜水含水层，把上述差分方程中的储水系数S改为给水度Sy，至于导
水系数T，则可以根据饱和厚度水头计算出每个节点所在单元的等效导水系数，以节点0为例，有
kT0?K0(H0?zb0) （4.2.29）
式中：K0为含水层在节点0处的渗透系数；zb0为含水层底板在节点0处的高程。这样就导致了一个新问题，在潜水含水层的差分方程（4.2.26）的中，系数C0～C4等与待求的
水头序列｛Hk｝有关，从而使式（4.2.26）变成非线性方程。解决这个问题的办法有以下几种。
1）显式法：用前一个时刻的已知水头序列｛Hk-1｝计算潜水含水层的导水系数，在时步内不再更新。当时间步长较大，水头变化明显时，这种方法精度较低。
2）预测-校正法：首先用显式法求出tk-1/2时刻的水头序列{Hk-1/2｝，时间步长设置为△tk/2，然后计算出中间时刻的导水系数，重新计算系数C0～C4，返回去由前一时刻的已知水头序列{Hk-1｝求出待求水头序列号{Hk｝，时间步长恢复为△tk。
3）迭代法：首先用显式法求出水头序列｛Hk｝，并计算出新的导水系数，重新计算系数C0～C4，返回去再次求解水头序列｛Hk｝，如此反复，直到前后两次计算得到的水头几乎相同（由允许误差控制）为止。这种方法较为精确，但计算量大。
基于矩形网格的地下水流三维模型，一般是进行分层处理，每层都是同样的矩形网格。网格单元为长方体，每个长方体内设置节点，编号的形式为（i，j，l），其中的l表示层序号。在建立地下水流差分方程时，除了同层的4个相邻长方体单元之外，还需要考虑上下相邻的2个长方体单元，因此每个内部节点一般包含6个相邻节点。关下分层矩形网格的三维地下水流差分模型，在有关MODFLOW模型的章节中将具体阐述。
三、基于多边形网格的地下水流差分方程
基于矩形网格的有限差分法虽然简单，但是在处理不规则的模型边界和局部加密方面存在缺陷，会导致网格浪费。基于多边形网格的有限差分可以适应不规则的边界，并且能够在不显著增加整体节点数的情况下加密局部网格。多边形网格的形成，需要先建立模型空间的三角形剖分离散网格，作为辅助网格，再用三角形单元各边的中垂线连接为各个节点所对应的多边形单元（图4.3）。
图4.3 一个平面多边形网格的局部图（灰色线为辅助三角形网格）
辅助三角形网格的剖分应注意（陈崇希等，1989）：①三角形的任一内角应小于90°，尽可能接近正三角形；② 三角形的顶点不能落在另外某个三角形的边上，即节点之间的连线最多与2个三角形相连；③节点的分布需要考虑地下水流动的控制条件，如河流、井点等。辅助三角形生成之后，再组建多边形网格，其方法类似于气象水文学中降雨量插值的泰森多边形法（Bear，1979），对每条三角形的边作垂直平分线，直到与其他的垂直平分线相交，这些垂直平分线连接成一个一个的多边形，每个多边形包围一个节点，作为该节点的控制区。对具有多边形特征的节点控制区建立地下水流偏微分方程的差分格式，可以利用封闭域内体积分与边界积分之间的关系（高斯定律），其方式与有限体积法一致。因此，在某种程度上说，基于多边形网格的有限差分法属于有限体积法。
不过，节点控制区的差分方程也可以在水均衡原理的基础上建立。设模型研究的对象是一个均质各向同性的承压含水层，从多边形网格的内部取一个节点控制区进行分析（图4.4）。
5Lii=4wiN=6312 图4.4 平面多边形网格的一个节点控制区
节点0与N个节点相邻，对于第i个节点，节点0与节点i之间连线的长度为Li，而对应的垂直平分线的宽度为wi，则在节点i与节点0之间流向控制区的流量Qi可根据达西定律计算为
Qi?TwiHi?H0 （4.2.30） Li式中：T为含水层的导水系数；Hi和H0分别为节点i与节点0的水头。所有相邻节点与控制区之间都可以建立这样的关系，因此从控制区周围流向控制区内的地下水的总流量应为
?Q?T?kii?1i?1NNwik(Hik?H0) （4.2.31） Li式（4.2.31）中把时间节点k也放进去了。在时间步长?tk内，控制区内地下水的总体积的增加量为
kk?1?Vw?SA0(H0?H0) （4.2.32）
式中：S为含水层储水系数；A0为节点0的控制区面积，式（4.2.32）中隐含的假设是
节点0处的水头变化与控制区内的平均水头变化相等。根据水均衡原理，控制区水体积的增加量应该与流入控制区的总水量相互平衡，即
(Qa??Qi)?t??Vw （4.2.33）
i?1N式中：Qa为其他的外部补给流量。把式（4.2.33）展开得
Qa?t?T?t?i?1Nwikkk?1(Hik?H0)?SA0(H0?H0) （4.2.34） Li此式可以进一步改写为
Nwikwk?1 （4.2.35） (SA0?T?t?)H0?T?t?iHik?Qa?t?SA0H0i?1Lii?1LiN这就是节点控制区的地下水流差分方程。对于多边形网格中的每个内部节点，都可以建
立相同形式的差分方程，与边界节点的差分方程一起形成求解模型的方程组。
每个相邻节点对应的wi/li数值，以及节点控制区的面积A0，可以根据辅助三角形网格的几何信息进行计算，具体方法见有关文献。
利用多边形网格建立三维有限差分模型的方法，一般是建立分层的多边形网格以模拟含水层的构造特征［图4.5（a）］。每个节点除了和同一层的若干节点相邻之外，还与上层或下层相同平面位置的节点相邻，以图4.5（b）中的节点0为例，其同层相邻节点为节点1～N，而节点u和节点l分别为上层相邻节点和下层相邻节点，它们与节点0连线的中垂面和平面控制区共同形成节点0的三维控制体，这个控制体的水量均衡方程为
kk?1(Qa??Qi)?t?(Qu?Ql)?t??Vw?SsV0(H0?H0) （4.2.36）
i?1N式中：Qu和Qi分别为控制体顶面和底面流向控制体内的流量；Ss为含水层的储水率；V0为节点O控制体的体积。控制体顶底面的垂向流量根据达西定律计算为