1、地下水建模方法和步骤地下水建模方法和步骤中国地质大学中国地质大学(武汉武汉)环境学院环境学院2012.8地下水建模方法和步骤地下水建模方法和步骤1.求解地下水运动方程的数值方法求解地下水运动方程的数值方法2.地下水数值模型建模步骤地下水数值模型建模步骤3.建模所需要的基本资料建模所需要的基本资料 绝大部分数学模型是无法用解析法求解的，数绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。值化就是将数学模型转化为可解的数值模型。1.数值方法数值方法1.1有限差分法有限差分法1.2有限单元法有限单元法1.3分有限差分法分有限差分法1.4半解析半数值法半解析半数值法1.5
2、边界元法边界元法o(1)有限差分法原理有限差分法原理o(2)两种方法建立有限差分方程两种方法建立有限差分方程o(3)求解有限差分方程求解有限差分方程o(4)收敛性和稳定性概念收敛性和稳定性概念o(5)算例算例1.1有限差分法有限差分法(1)有限差分法的基本原理有限差分法的基本原理将连续的问题离散后求解：将连续的问题离散后求解：方法一方法一以地下水流基本微分以地下水流基本微分方程及其定解条件为基础，方程及其定解条件为基础，在在渗流区剖分基础上，用差商代渗流区剖分基础上，用差商代替微商，将地下水流微分方程替微商，将地下水流微分方程的求解转化为差分方程（代数的求解转化为差分方程（代数方程）求解。方程
3、）求解。方法二方法二在渗流区剖分的基础在渗流区剖分的基础上，直接由达西定律和水均衡上，直接由达西定律和水均衡原理，建立各个均衡区的水均原理，建立各个均衡区的水均衡方程，即差分方程。衡方程，即差分方程。矩形网格矩形网格多边形网格多边形网格网格划分的基本类型网格划分的基本类型o（1）先划格线，格）先划格线，格点位于网格中心点位于网格中心均衡网格均衡网格节点网格节点网格o（2）先规定格点位）先规定格点位置，再垂直平分两相置，再垂直平分两相邻结点的连线作格线，邻结点的连线作格线，形成的网格即为水均形成的网格即为水均衡区衡区MODFLOW网格系统网格系统x000()()()f xxf xfxx 导导数数
4、的的有有限限差差商商近近似似导数的定义导数的定义 当当非常小的时候，有非常小的时候，有 上式右端项即为上式右端项即为f(x)(x)在在x0 0处的差商。处的差商。这样定义的差商很容易理解，但不知道用差商代替微商所产生的误差。下面利用泰勒公式导出差商及其误差。xxfxxfxfx)()(lim)(0000方法一：方法一：差商代替微商差商代替微商(2)有限差分方程建立有限差分方程建立2000()()()()()2!ff xxf xfxxx 2000()()()()()2!ff xxf xfxxx 000()()()()f xxf xfxOxx 已知泰勒公式 由A得：AB 由B 得：000()()()
5、()f xf xxfxOxx 称称 为为f(x)在在x0处的处的一阶前向差商，一阶前向差商，为为截断误差截断误差。xxfxxf)()(00)(xO 称称 为为f(x)在在x0处的处的一阶后向差商，一阶后向差商，为为截断误差截断误差。xxxfxf)()(00)(xO 方法一方法一200002()2()()()()()f xxf xf xxfxOxx 由A-B可以得：由A+B可以得：2000()()()()2f xxf xxfxOxx(4)23400000()()()()()()()()()2!3!4!fxfxff xxf xfxxxxx (4)23400000()()()()()()()()()
6、2!3!4!fxfxff xxf xfxxxxx AB称称 为为f(x)在在x0处的处的一阶中心差商，一阶中心差商，为为截断误差截断误差。xxxfxxf2)()(002)(xO 称称 为为f(x)在在x0处的处的二阶二阶中心中心差商，差商，为为截断误差截断误差。2000)()()(2)(xxxfxfxxf2)(xO 方法一方法一l对于偏导数（偏微商），类似可以得到相应的差商：对于偏导数（偏微商），类似可以得到相应的差商：ttxHttxHttxH),(),(),(000000 xtxHtxxHxtxH),(),(),(00000020000002002)(),(),(2),(),(xtxxHtx
7、HtxxHxtxH方法一方法一(2)有限差分方程建立有限差分方程建立(续续)一维控制方程差分格式一维控制方程差分格式ttxhxtxhT),(),(22ttxHttxHxtxxHtxHtxxH),(),()(),(),(2),(00002000000ttxHttxHxttxxHttxHttxxH),(),()(),(),(2),(00002000000显式差分格式显式差分格式隐式差分格式隐式差分格式方法一方法一1,3,2 )(21211 nxitHHxHHHninininini控制方程控制方程1,3,2 )(21211111 nxitHHxHHHninininini网格剖分网格剖分nx个个),(
8、|tzyxxxv),(|tzyxxvtzyvxtzyxx),(|)(方向流入tzyvxtzyxxx),(|)(方向流出 取右图所示得微小六面取右图所示得微小六面体。设与体。设与x,y,z,方向对应得方向对应得主渗透系数分别为主渗透系数分别为Kx,Ky,Kz；建立均衡期；建立均衡期 t时时段内，微小均衡六面体的段内，微小均衡六面体的水量守恒方程。水量守恒方程。方法二：方法二：达西定律和水均衡原理达西定律和水均衡原理(2)有限差分方程建立有限差分方程建立(续续)基于基于达西定律达西定律，x，y，z方向流入方向流入流出分别为流出分别为:l t时段内，侧向流入与源汇项导致六面体水量变化量为：时段内，侧
9、向流入与源汇项导致六面体水量变化量为：tyxzyHHyHHtzxvvkjikjikjikjiyyyyy|,1,1,tyxzxHHxHHtzyvvkjikjikjikjixxxxx|,1,1tyxzzHHzHHtyxvvkjikjikjikjizzzzz|,1,1,ABC(2)有限差分方程建立有限差分方程建立(续续)方法二：方法二：达西定律和水均衡原理达西定律和水均衡原理DtyxzA+B+C+D源汇项源汇项六面体内地下水储存量的变化为六面体内地下水储存量的变化为由水均衡原理得三维地下水流动方程的有限差分格式由水均衡原理得三维地下水流动方程的有限差分格式),(),(tzyxHttzyxHzyxHz
10、yxssnkjinkjiskjikjikjikjikjikjikjikjikjikjikjikjiHHzyxtyxztyxzzHHzHHtyxzyHHyHHtyxzxHHxHH,1,1,1,1,1,1,1(2)有限差分方程建立有限差分方程建立(续续)方法二：方法二：达西定律和水均衡原理达西定律和水均衡原理有限差分法有限差分法:三维（三维（MODFLOW）差商代替微商差商代替微商(3)差分方程求解差分方程求解l一维显式差分格式 )(2)(2 )(2)(2111212212212331324322122321tHHxHHHtHHxHHHtHHxHHHtHHxHHHnninninninninninn
11、inninninninninnnnnnnnnn网格个数为网格个数为ni直接求解直接求解(3)差分方程求解差分方程求解l一维隐式差分格式一维隐式差分格式 )(2)(2 )(2)(2111211112212211121331321413122122131211tHHxHHHtHHxHHHtHHxHHHtHHxHHHnninninninninninninninninninninnnnnnnnnn网格个数为网格个数为ni迭代求解迭代求解方方程程组组PCGSIPSORWHSSAMGGMGMODFLOW(4)差分方程的收敛性和稳定性差分方程的收敛性和稳定性l截断误差截断误差：用差商代替微商时，地下水流动方
12、程产生用差商代替微商时，地下水流动方程产生的误差为截断误差。的误差为截断误差。l收敛性收敛性：当空间步长和时间步长趋于当空间步长和时间步长趋于0 0时，有限差分方时，有限差分方程的精确解趋于地下水流动问题微分方程定解问题的程的精确解趋于地下水流动问题微分方程定解问题的精确解。则称该差分格式是收敛的。精确解。则称该差分格式是收敛的。l稳定性稳定性：如果在求解差分方程过程中，某时间步引入如果在求解差分方程过程中，某时间步引入某个误差，而在以后的各时段计算中，该误差不再扩某个误差，而在以后的各时段计算中，该误差不再扩大，则称该差分格式是稳定的。大，则称该差分格式是稳定的。102一维一维显示格式显示格
13、式的收敛条件和稳定条件是：的收敛条件和稳定条件是：2()()OtOx(6)算例：算例：显式有限差格式显式有限差格式 设两条河流平行、完全切割含水层，含水层等厚、均质设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。各向同性。应用实例：河间地块承压水流模型应用实例：河间地块承压水流模型步骤：步骤：（1 1）基础资料的分析）基础资料的分析（2 2）概念模型）概念模型（3 3）数学模型）数学模型（4 4）数值方法及计算机程序）数值方法及计算机程序（5 5）参数）参数（6 6）结果分析）结果分析 建立数学模型建立数学模型0,0tLxLx 022eHHKMtx)(00 xHHt)(),(210t
14、HtHLXx（1 1）模型概化）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。由所述水文地质条件，可以概化为一维承压水流问题。（2 2）建立坐标系）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内，（如图），将地下水流动系统空间结构放在坐标系内，从而量化各变量的取值范围。本例，取从而量化各变量的取值范围。本例，取x-x-轴原点位于左端河，右侧轴原点位于左端河，右侧为正向，设两河流间距为为正向，设两河流间距为L.L.（3 3）数学模型）数学模型0t差分方程及其解法差分方程及其解法显式格式显式格式将（0L）分成 N 等份，Nlx 1)1)网格剖分：网格剖分：xixittntn取
15、时间步长取时间步长 ，记，记 （n n=0=0、1 1、2 2、3 3、44）记 ，(i=0,1,2,3，4N)2 2）建立差分方程）建立差分方程：在网格系统中任意取一点：在网格系统中任意取一点),(nitx),(txH),(nitx设设是问题的解，则在是问题的解，则在处有处有记为（记为（i i，n n）22(,)nneiniiHHKMx ttx22(,)nneiniiHHKMx ttx1()nnniiiHHHOttt2211222()()nnnniiiiHHHHOxxx用差商代替微商用差商代替微商：将上述两式舍去余项，代入方程并记将上述两式舍去余项，代入方程并记niHnih为为11122(,
16、)()nnnnniiiiiinehhhhhKMx txt2)()(xOtO显然该式具有截断误差得到得到显式格式（续显式格式（续1）111(12)(,)nnnniiiiinethhhhx t引入无量纲变量引入无量纲变量：将该式子代入得到：将该式子代入得到：（i=1，2，3，.N-1），（n=1，2，3，.）显式格式（续显式格式（续2）2)(xtKMe11122(,)()nnnnniiiiiinehhhhhKMx txt3）显示差分方程的求解）显示差分方程的求解l计算各结点初始时刻水头值计算各结点初始时刻水头值l利用差分方程计算各结点利用差分方程计算各结点t1t1时刻水头值时刻水头值l利用边界条件
17、计算边界结点水头值利用边界条件计算边界结点水头值l重复重复2 2、3 3步，直到计算出拟计算的各个时刻的水头步，直到计算出拟计算的各个时刻的水头值值)(00 xHHt显式格式（续显式格式（续3）22eHHKMtx)(00 xHHt)(),(210tHtHLXx算例（续算例（续4）004.0em10)(m,20)(m,20 m/d,821ttMK在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200200米，时间步长为米，时间步长为0.250.25天，分别计算各节点天，分别计算各节点各各时刻的时刻的水头值。水头值。0 x20m,0 x10m,)(0 xH4/1)200(*
18、004.025.0*20*8)(22xtKMeTime/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 m02010101010100.252012.5101010100.502013.7510.6251010100.752014.531 11.250 10.156 10101.002015.078 11.797 10.391 10.039 101.252015.488 12.266 10.654 10.117 10nieninininithhhh,111/)21()2(4141)4121(4111111nininininininihhhhhhh算例（续算
19、例（续5）Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 m02010101010100.252020101010100.502010201010100.7520300 2010101.0020-10 101.252010算例（续算例（续6）1)200(*004.01*20*8)(22xtKMenieninininithhhh,111/)21(nininininininihhhhhhh11111)121(如果如果 t=1,t=1,则则(6)算例：算例：隐式格式隐式格式22eHHKMtx)(00 xHHt)(),(210tHtHLXx004.0e在
20、上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200米，时间步长为米，时间步长为0.25天，用隐式差分格式计算各节天，用隐式差分格式计算各节点个时刻的水头值。点个时刻的水头值。m10)(m,20)(10m,)(m,20 m/d,821ttxHMK11122332210111112131211/./.2121.212121nNnNenNnNenNnennennnennNnNnnnhththththhthhhhhh4/1)200(004.025.0208)(22xtKMe在这个例子中，在这个例子中，.0,1,2,3,4,m,10,m2050nhhnnm10,m10,m10,
21、m1004030201hhhh154321011413121125.025.0225.00025.0225.00025.0225.00025.02nnnnnnnnnnhhhhhhhhhh解：隐式格式一般方程为解：隐式格式一般方程为于是有于是有根据初始条件得根据初始条件得根据边界条件得根据边界条件得.0,1,2,3,4,m,10,m2050nhhnnm10,m10,m10,m1004030201hhhh1504030210011413121125.025.0225.00025.0225.00025.0225.00025.02hhhhhhhhhh由初始条件和边界条件由初始条件和边界条件5.1210
22、10151025.01010102025.010225.00025.0225.00025.0225.00025.0214131211hhhh由此解得由此解得t1t1时刻的水头值为时刻的水头值为m10,m10,m10,m1014131211hhhh在上述方程中取在上述方程中取 n=0,n=0,可以得到计算可以得到计算t1t1时刻水头值的方程时刻水头值的方程所以上述方程变成所以上述方程变成同理，可计算同理，可计算t2t2时刻的水头值时刻的水头值2.地下水数值模型建模步骤地下水数值模型建模步骤模拟步骤模拟步骤o建立概念模型建立概念模型 o建立数学模型建立数学模型 o数值方法及软件数值方法及软件（编程
23、）（编程）o参数参数 o拟合模拟：模型校拟合模拟：模型校正与检验正与检验 o参数敏感性分析参数敏感性分析 o预测模拟预测模拟软件软件一、概念模型（模型概化）一、概念模型（模型概化）根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、岩石矿物、水文、气象、工农业利用情况等岩石矿物、水文、气象、工农业利用情况等模拟的区域模拟的区域:含水层类型含水层类型:潜水（无压）、承压、混合、多层潜水（无压）、承压、混合、多层维数维数:一维、二维、三维一维、二维、三维水流状态：稳定流水流状态：稳定流/非稳定流、饱和流非稳定流、饱和流/非饱和流非
24、饱和流介质状况介质状况:均质和非均质均质和非均质/各向同性和各向异性各向同性和各向异性孔隙孔隙/裂隙裂隙/双重介质双重介质流体的密度差流体的密度差边界条件和初始条件边界条件和初始条件 必要时需进行一系列的室内试验与野外试验必要时需进行一系列的室内试验与野外试验,以获取有关参数以获取有关参数,如渗透系如渗透系数、弥散系数、分配系数、反应速率常数等。数、弥散系数、分配系数、反应速率常数等。收集研究区已有收集研究区已有水文地质资料水文地质资料概念模型概念模型边界概化边界概化内部结构概化内部结构概化完成模型概化图完成模型概化图地表水体地表水体水头边界水头边界流量边界流量边界断层接触边界断层接触边界隔水
25、边界隔水边界水头边界水头边界流量边界流量边界岩体岩层接触边界岩体岩层接触边界天然分水岭天然分水岭据含水层组类型、据含水层组类型、结构、岩性结构、岩性含水介质含水介质地下水运动状态地下水运动状态水文地质参数水文地质参数确定层组的均质、非均质，确定层组的均质、非均质，各向同行、异性，稳定、非各向同行、异性，稳定、非稳定，潜水、承压水稳定，潜水、承压水空间分布空间分布导水、储水系数、主渗透方向导水、储水系数、主渗透方向目的层与相邻层关系目的层与相邻层关系均质、非均质，各向同、异性均质、非均质，各向同、异性层、紊流，二、三维层、紊流，二、三维时间概化时间概化渗透、储水系数，给水度及单位涌渗透、储水系数
26、，给水度及单位涌水量，含水层分布规律，地下水流水量，含水层分布规律，地下水流场，水化学场、温度场的空间概化场，水化学场、温度场的空间概化平面图平面图剖面图剖面图研究范围、主要居民点、标志性地形、地貌、地表河流、研究范围、主要居民点、标志性地形、地貌、地表河流、湖泊、开采井、地下水天然露头、含水层的各类钻孔湖泊、开采井、地下水天然露头、含水层的各类钻孔地表地理要素、含水层结构、地质结构、地下水水位、地表地理要素、含水层结构、地质结构、地下水水位、各类源汇项及其性质各类源汇项及其性质二、数学模型二、数学模型 三维地下水流动问题控制方程三维地下水流动问题控制方程1|(,)BHH x y t00|(,
27、)tHHx y2|(,)BHTq x y tntHqzHKzyHKyxHKxszyx)()()(第二类边界条件第二类边界条件第一类边界条件第一类边界条件初始条件初始条件 绝大部分数学模型是无法用解析法求解的，绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。数值化就是将数学模型转化为可解的数值模型。三、数值方法及软件（或编程）三、数值方法及软件（或编程）有限差分法有限差分法有限单元法有限单元法积分有限差分法积分有限差分法半解析半数值法半解析半数值法边界元法边界元法有限差分法有限差分法:MOFLOW系列系列GMS 中中MODFLOWVisual MODFLOWPro
28、cessing MODFLOW有限单元法有限单元法:FEFLOW积分有限差分法积分有限差分法:TOUGH2，TOUGH REACT软件软件四、模型参数四、模型参数o含水层参数：渗透系数，弹性释水系数（重力给含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等水度），孔隙度等o源汇项：源汇项：n大气降水入渗系数（分区、数值）大气降水入渗系数（分区、数值）n蒸发排泄系数蒸发排泄系数n地表水体（河流、湖泊、水库等）水位、底面高程、底地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等）面岩性特征（厚度、渗透系数等）n渠系灌溉入渗系数渠系灌溉入渗系数n人工开采（点状、面状
29、）人工开采（点状、面状）o边界条件边界条件o初始条件初始条件参数的不确定性参数的不确定性钻孔太少，钻孔太少，地层资料少，地层资料少，钻孔多，含钻孔多，含水层结构会水层结构会发生变化发生变化五、模拟：模型校正（参数识别五、模拟：模型校正（参数识别）将模拟结果与实测结果比较将模拟结果与实测结果比较,进行参数调整进行参数调整,使模使模拟结果在给定的误差范围内与实测结果吻合。拟结果在给定的误差范围内与实测结果吻合。调参过程是一个复杂而辛苦的工作调参过程是一个复杂而辛苦的工作,所调整的参数所调整的参数必须符合模拟区的具体情况。必须符合模拟区的具体情况。人机交互与自动调参相结合。尽管自动调参程序人机交互与
30、自动调参相结合。尽管自动调参程序(如如PEST),PEST),也不能代替人的工作。也不能代替人的工作。五、模拟：模型检验五、模拟：模型检验 模型验证是在模型校正的基础上模型验证是在模型校正的基础上,进一步进一步调整参数调整参数,使模拟结果与第二次实测结果吻合使模拟结果与第二次实测结果吻合,以进一步提高模型的置信度。以进一步提高模型的置信度。六、灵敏度分析六、灵敏度分析 校正后的模型受参数值的时空分布、边界校正后的模型受参数值的时空分布、边界条件、水流状态等条件、水流状态等不确定度不确定度的影响。灵敏度的影响。灵敏度分析就是为了确定不确定度对校正模型的影分析就是为了确定不确定度对校正模型的影响程
31、度。响程度。七、预测七、预测 用校正的参数值进行预测用校正的参数值进行预测,预测时需估算未预测时需估算未来的水流状态。来的水流状态。后续检查与模型的重建（完善）后续检查与模型的重建（完善）后续检查在模拟研究结束数年后进行。收集新的野后续检查在模拟研究结束数年后进行。收集新的野外数据以确定预测结果是否正确。如果模拟结果精确外数据以确定预测结果是否正确。如果模拟结果精确,则该模型对该模拟区来说是有效的。由于场址的唯一则该模型对该模拟区来说是有效的。由于场址的唯一性性,故模型只对该模拟区有效。故模型只对该模拟区有效。后续检查应后续检查应在预测结束足够长的时间后进行在预测结束足够长的时间后进行,以以便
32、有足够的时间发生明显的变化。便有足够的时间发生明显的变化。模型的再设计模型的再设计 一般来说一般来说,后续检查会发现系统性能的变后续检查会发现系统性能的变化化,从而导致概念模型和模型参数的修改。一从而导致概念模型和模型参数的修改。一般来说般来说,所有模拟研究都应该进行到第五步所有模拟研究都应该进行到第五步,即校正灵敏度分析。即校正灵敏度分析。地裂缝地裂缝自然地理及水文地质条件自然地理及水文地质条件边界敏感性分析边界敏感性分析数学模型数学模型网格化网格化概念模型概念模型介质类型、结构特征介质类型、结构特征地下水补、径、排特征地下水补、径、排特征数值模型数值模型数值模型的建立数值模型的建立边界条件
33、边界条件初始条件初始条件源汇项源汇项降雨入渗降雨入渗地表水体入渗地表水体入渗灌溉、渠系灌溉、渠系蒸发排泄蒸发排泄人工开采人工开采介质参数介质参数观测孔动观测孔动态拟合态拟合流场拟合流场拟合水均衡对比分析水均衡对比分析观测孔动态检验观测孔动态检验流场检验流场检验水均衡对比检验水均衡对比检验拟合调参拟合调参识别结果识别结果不符合要不符合要求修正概求修正概念模型念模型小结：地下水流动小结：地下水流动模型构建过程模型构建过程识别结果识别结果不符合要不符合要求修正概求修正概念模型念模型模型应用模型应用水资源量评价水资源量评价预测开采方案预测开采方案研究环境地质问题研究环境地质问题地面沉降地面沉降地下水污
34、地下水污染研究染研究大型工程大型工程对环境的对环境的影响影响模型识别模型识别模型检验模型检验3.地下水流数值模拟资料需求地下水流数值模拟资料需求o空间展布相关资料空间展布相关资料n地表高程、分层数据地表高程、分层数据n边界位置边界位置o含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等o源汇项：源汇项：n大气降水入渗系数（分区、数值）大气降水入渗系数（分区、数值）n蒸发排泄系数蒸发排泄系数n地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等）数等）n渠系灌溉入渗系数渠系灌溉入渗系数n人工开采（点状、面状）：开采井位置、井结构、开采量动态人工开采（点状、面状）：开采井位置、井结构、开采量动态o边界条件：边界类型、水头或流量边界条件：边界类型、水头或流量o初始条件：统测水位初始条件：统测水位o水位动态观测资料：观测孔位置、结构、水位时间变化水位动态观测资料：观测孔位置、结构、水位时间变化
