第四章 地下水运动的数值模型
解析解虽然具有精确可靠的特点，但采用解析解反映自然状态和复杂人类活动干扰下的地下水运动是相当困难的。因此，当含水层的条件严重偏离现有解析模型的简化假设时，人们通过数值模型来获得近似的地下水流场及演变趋势。
第一节 地下水流数值方法概述
地下水流的数学模型采用偏微分方程描述地下水流的时间和空间连续状态，而数值模型则是采用离散（非连续）时空模型中水头的分布与演变对数学模型进行近似描述。从精确数学模型到近似数值模型的转化，虽然会损失一些精度，但使复杂地下水流问题的分析得以通过机械计算实现，而且误差也是可控的。
把偏微分方程求解的数值方法引入到地下水流问题的求解始于20世纪70年代，主要方法包括有限差分法、有限元法和边界元法，此后又发展了有限分析法、多重网格法和无网格法等。这些方法的共同特点是将模型空间及边界离散为由一系列的节点以及联系这些节点的单元（无网格法除外），含水层的水头在这些节点上定义，从而实现了水头分布空间连续函数向离散变量的转化，表示为
?2fk?fkfk?fk?1a?b?c?e，0?x?2L2?x?x?tk?tk?tk?tk?1fk(x)?2fk?1?fk?1fk?fk?1a?b?c?e，0?x?2L2?x?x?tk?tk?f?f(a?b?c)，0?x?2L2e?x?xd2fkdfkekefk?1a?b?f?(c?)?0，0?x?2L2dxdx?tk?tkfk?fk?1?fk(x?0)?uH(x,y,z)??Hp;p?1,2,3,???,M? （4.1.1）
式中；H为含水层的水头；x、y、z为空间坐标；p为数值模型的节点；M为节点的数目。
与此同时，时间也被离散化为一系列具有先后顺序的时刻，不同时刻节点的水头分布方式也不同，每个时刻节点的水头分布结果总是由前一个时刻演变过来，即
2k?1k?1
?Hp;p?1,2,3,???,M;t?tk;?k??Hk?1p;p?1,2,3,???,M;t?tk??tk; （4.1.2）
?式中：k为时刻；t为时间；?tk为时间步长。数值模型的核心部分，是在地下水流数学模型的基础上确定节点在流场中的相互关系，并形成如下形式的方程组
CHp1k?11?Cp2H2?????CppHp?????CpMHM?Fp,p?1,2???,M（4.1.3）
k?1k?1k?1k式中：系数
C反映了节点P与节点L在流场中的关系，可称为关联系数；pLFkp通常是由
前一个时刻决定的常量。式（4.1.3）除包括地下水流控制方程的离散近似以外，还包括了边界条件的离散化处理。当k=0时，
Fk中包含了初始条件。如果模拟稳定流，则pFkp与时
间无关，上标k可以略去。使用矩阵符号，式（4.1.3）可缩写为
?C?M?M?H???F? （4.1.4）
实际在处理模型时，一个节点通常只与相隔最近的若干节点有直接的关系，在式（4.1.3）中大多数关联系数为零，而且一般有对称性，即
CPL?CLP。因此，矩阵?C?通常是对称的
稀疏矩阵，降低了方程组求解的困难。尽管如此，区域地下水流数值模型的节点数往往很大，而且在某些非线性条件下，?C?本身是待求水头的函数，方程组的求解需要使用迭代法而耗费相当多的时间。
不同类型的地下水流数值模型，最核心的差异在于形成式（4.1.3）的方式，关联系数具有不同的计算公式。原则上，对于完全相同的地下水流数学模型，以及完全相同的离散节点序列和离散时间序列，不同的数值模型应给出近乎同样的结果，只是计算精度略有差异。已有研究表明，某种类型的有限差分法与某种类型的有限元法存在一定的等价性，但更多的对比分析尚待进一步探讨。实际上，研究者在检验数值模型的可靠性时，总是与解析解的模型对比。当不同方法的数值模型用相同的解析解作验证时，他们之间的等价性就可以得到一定程度的证明。然而，人们对数值模型的偏好并不取决于数值模型之间的等价性，而往往取决于某种数值方法的熟悉程度，或者基于某种数值方法专业软件的方便快捷。
专业软件在地下水流数值模型的推广应用中发挥了重要作用，如有限差分模拟程序MODFLOW（McDonald and Harbaugh，1988）和HST3D（Kipp，1987）、有限元模拟程序FEMWATER（Lin et al．，1997）和FEFLOW（Diersch and Kolditz，1998）等，国内有陈崇希等开发的有限差分模拟程序PGMS（陈崇希等，2007）。其中，以美国地质调查局发布的有限差分模拟程序MODFLOW最为著名，本章最后将介绍其特点。
第二节 地下水流有限差分模型
有限差分法是地下水流数值模拟最早使用的一种方法。这种方法古典而简洁、物理意义明了、编程容易，因而被广泛流传。地下水流的有限差分模型已经从最初的一维均质模型发展为目前的三维非均质模型，并有专业软件加以实现。
一、有限差分法的基本原理
有限差分法是求解偏微分方程边值问题和初值问题的一种数值方法，其实质是把连续的模型空间离散化为规则或不规则的网格点，利用导数的差分近似形式代替偏微分方程形成差分方程组，通过求解方程组得到离散点的待求变量作为连续场的一种近似结果。
模型空间的离散化从形态上主要分为规则网格与不规则网格。在二维平面空间，规则网格由一系列与坐标轴平行的直线组成，直线交叉形成的单元（cell）为矩形；在三维空间，规则网格由一系列与坐标轴正交的平面组成，平面相互切割形成的块体（block）为长方体。不规则网格在平面上由一系列互不重叠的多边形填充模型空间而形成，一般通过分层建立三维空间网格。一维问题空间离散化很简单，就是在坐标轴上生成一系列具有一定间隔的点，或称节点（node），相邻两个节点之间为线段，待求变量在这些节点上的分布代表了在坐标轴上连续分布的近似状态。有限差分法中的“有限”是指网格中的节点、单元或块体数目是有限的，每个线段、单元或块体的尺度也是有限的。
首先用一维问题来阐述导数的差分格式。设待求变量f沿x轴的分布为二阶连续函数
f(x)，现将x轴离散化为节点(x0,x1,x2,???,xi,???,xN)，它们的待求变量为
(f0,f1,f2,???,fi,???,fN)，则f(x)在节点xi处的一阶导数表示为
?f?x?x?xifi?1?fi????xi? （4.2.1）
xi?1?xi式中：?xi?xi?1?xi，为节点的间距。去掉式（4.2.1）中的Taylor展开式余项，就得到一阶导数前向差分格式。一阶导数还可以表示成后向差分格式为
?f?x?x?xifi?fi?1 （4.2.2）
xi?xi?1式中：?xi?1?xi?xi?1。显然，第一个节点x0处的一阶导数只能表示为前向差分格式，而最后一个节点xN处的一阶导数只能表示为后向差分格式。
有了一阶导数的差分格式之后，再次使用中心差分格式，就可以建立二阶导数的差分格式为
?f?x22?x?xi(?f?x?xxi?1?xi?1)i?1?(?f)i?1 （4.2.3）
其中：一阶导数在节点xi?1采用后向差分格式，在节点xi?1采用前向差分格式，则可以得到二阶导数的差分格式为
?2f?x2?x?xi(xi?xi?1)fi?1?(xi?1?xi?1)fi?(xi?1?xi)fi?1 （4.2.4）
(xi?1?xi?1)(xi?1?xi)(xi?xi?1)如果相邻节点的间距都相等，即?xi??xi?1??x，则上式可简化为
?2f?x2?x?xifi?1?2fi?fi?1 （4.2.5）
2?x2建立了一阶和二阶导数的差分格式，就可以用它们来改写所研究问题的偏微分方程。
下面用一个简单的例子来说明有限差分模型的建立方法。设有一边值问题的数学模型为
d2fdfa2?b?c?0,0?x?2L （4.2.6） dxdxf(x?0)?u （4.2.7）
dfdx?v （4.2.8）
x?2L在建立有限差分模型时，把模型空间[0，2L]离散化为三个节点。对第二个节点(x0?0,x1?L,xL)x1建立控制方程（4.2.6）的差分格式为 3?2af2?2f1?f0f2?f1?b?c?0 （4.2.9）
2L2L第一个节点x0为边界节点，直接利用边界条件式（4.2.7）得
f0?u （4.2.10）
第三个节点x2也是边界节点，把边界条件式（4.2.8）表示成差分格式为
f2?f1?v （4.2.11） L这样，由式（4.2.6）至式（4.2.8）描述的数学模型就转化为了由方程组[式（4.2.9）至式（4.2.11）]构成的差分模型。求解这个差分方程组，得
2bL22cL22bL22cL2f0?u，f1?(L?)v?u?，f2?(2L?)v?u? （4.2.12）
aaaa待求变量的差分模型结果(f0，f1，f2)就反映了数学模型精确解f(x)得特征。
如果我们遇到的是一个初边值问题，例如，把上述关于f(x)的边值问题改写为如下关于f(x，t)的数学模型
?2f?f?fa2?b?c?e，0?x?2L （4.2.13） ?x?x?tf(x?0,t?0)?u （4.2.14）
?f(x，t?0)?v （4.2.15）
?xx?2Lf(x，t?0)?0，0?x?L （4.2.16）
其中，式（4.2.16）是初始条件。对于这种初边值问题，我们需要先解决关于时间的偏导数。把时间坐标轴离散化为时步(t0，t1，t2，???，tk，???，tM)，则第k时刻的控制方程可以写为
?2fk?fkfk?fk?1a2?b?c?e，0?x?2L （4.2.17） ?x?x?tk式中：?tk?tk?tk?1为时间步长。时间偏导数取后向差分格式，方程左侧偏导数的自变量全部取第k时刻的数值，这种时间上的处理称为隐式差分格式。控制方程也可以处理为显式
差分格式，即
?2fk?1?fk?1fk?fk?1a?b?c?e，0?x?2L （4.2.18） 2?x?x?tk这样可以直接根据前一个时刻的变量分布fk?1计算当前时刻的fk，为
f?fkk?1?tk?2fk?1?fk?1?(a?b?c)，0?x?2L （4.2.19） e?x2?x关于时间的偏导数问题解决以后，就可以把上述初始问题的求解转化为一系列离散时刻关于
fk(x)的边值问题的求解。例如，采用隐式差分格式，则fk(x)的数学模型为
d2fkdfkekefk?1a?b?f?(c?)?0，0?x?2L （4.2.20） 2dxdx?tk?tkfk(x?0)?u （4.2.21）
dfkdx?v （4.2.22）
x?2L利用边值问题的有限差分法，可以把上述数学模型转化为有限差分模型，得到各个节点上的待求变量，只要把前一时刻求解出来的结果fk-1代入到式（4.2.20）即可。数值模型求解的流程（图4.1）。
k=0 fk?f(t?0) k?k?1,?t?tk?tk?1 k?1?fik 差分模式：fi????否 k = M 是 结束
图4.1 初边值问题有限差分模型的求解流程示意图
有限差分模型必须满足收敛性和稳定性才是有价值的。收敛性的含义是当离散网格的△x、△y、△z、△t趋于零时，有限差分模型的离散解应趋于数学模型的精确解，即随着网格剖分尺度和时间步长的缩小，模型的误差也越来越小。稳定性则与初边值问题在时间上的推进有关，初边值问题在有限差分模型中是转化为各个时步的边值问题求解的，如果前一个