《地下水建模方法和步骤》由会员分享，可在线阅读，更多相关《地下水建模方法和步骤（50页珍藏版）》请在金锄头文库上搜索。
1、地下水建模方法和步骤地下水建模方法和步骤中国地质大学中国地质大学(武汉武汉)环境学院环境学院2012.82021/6/161地下水建模方法和步骤地下水建模方法和步骤1.求解地下水运动方程的数值方法求解地下水运动方程的数值方法2.地下水数值模型建模步骤地下水数值模型建模步骤3.建模所需要的基本资料建模所需要的基本资料 2021/6/162 绝大部分数学模型是无法用解析法求解的，数绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。值化就是将数学模型转化为可解的数值模型。1.数值方法数值方法1.1有限差分法有限差分法1.2有限单元法有限单元法1.3分有限差分法分有限差分
2、法1.4半解析半数值法半解析半数值法1.5边界元法边界元法2021/6/163o(1)有限差分法原理有限差分法原理o(2)两种方法建立有限差分方程两种方法建立有限差分方程o(3)求解有限差分方程求解有限差分方程o(4)收敛性和稳定性概念收敛性和稳定性概念o(5)算例算例1.1有限差分法有限差分法2021/6/164(1)有限差分法的基本原理有限差分法的基本原理将连续的问题离散后求解：将连续的问题离散后求解：方法一方法一以地下水流基本微分以地下水流基本微分方程及其定解条件为基础，方程及其定解条件为基础， 在在渗流区剖分基础上，用差商代渗流区剖分基础上，用差商代替微商，将地下水流微分方程替微商，将
3、地下水流微分方程的求解转化为差分方程（代数的求解转化为差分方程（代数方程）求解。方程）求解。方法二方法二在渗流区剖分的基础在渗流区剖分的基础上，直接由达西定律和水均衡上，直接由达西定律和水均衡原理，建立各个均衡区的水均原理，建立各个均衡区的水均衡方程，即差分方程。衡方程，即差分方程。矩形网格矩形网格多边形网格多边形网格2021/6/165网格划分的基本类型网格划分的基本类型o（1）先划格线，格）先划格线，格点位于网格中心点位于网格中心均均衡衡网网格格节节点点网网格格o（2）先规定格点位）先规定格点位置，再垂直平分两相置，再垂直平分两相邻结点的连线作格线，邻结点的连线作格线，形成的网格即为水均形
4、成的网格即为水均衡区衡区2021/6/166MODFLOW网格系统网格系统2021/6/167导导数数的的有有限限差差商商近近似似导数的定义导数的定义 当当非常小的时候，有非常小的时候，有 上式右端项即为上式右端项即为f(x)(x)在在x0 0处的差商。处的差商。 这样定义的差商很容易理解，但不知道用差商代替微商所产生的误差。下面利用泰勒公式导出差商及其误差。方法一：方法一：差商代替微商差商代替微商(2)有限差分方程建立有限差分方程建立2021/6/168已知泰勒公式 由A得： AB 由B 得： 称称 为为f(x)在在x0处的处的一阶前向差商，一阶前向差商， 为为截断误差截断误差。称称 为为f
5、(x)在在x0处的处的一阶后向差商，一阶后向差商， 为为截断误差截断误差。方法一方法一2021/6/169 由A-B可以得： 由A+B可以得：AB称称 为为f(x)在在x0处的处的一阶中心差商，一阶中心差商， 为为截断误差截断误差。称称 为为f(x)在在x0处的处的二阶二阶中心中心差商，差商， 为为截断误差截断误差。方法一方法一2021/6/1610l对于偏导数（偏微商），类似可以得到相应的差商：对于偏导数（偏微商），类似可以得到相应的差商：方法一方法一(2)有限差分方程建立有限差分方程建立(续续)2021/6/1611一维控制方程差分格式一维控制方程差分格式显式差分格式显式差分格式隐式差分格
6、式隐式差分格式方法一方法一控制方程控制方程网格剖分网格剖分nx个个2021/6/1612 取右图所示得微小六面取右图所示得微小六面体。设与体。设与x,y,z,方向对应得方向对应得主渗透系数分别为主渗透系数分别为Kx, Ky,Kz；建立均衡期；建立均衡期 t时时段内，微小均衡六面体的段内，微小均衡六面体的水量守恒方程。水量守恒方程。方法二：方法二：达西定律和水均衡原理达西定律和水均衡原理(2)有限差分方程建立有限差分方程建立(续续)2021/6/1613基于基于达西定律达西定律，x，y，z方向流入方向流入流出分别为流出分别为:l t时段内，侧向流入与源汇项导致六面体水量变化量为：时段内，侧向流入
7、与源汇项导致六面体水量变化量为：ABC(2)有限差分方程建立有限差分方程建立(续续)方法二：方法二：达西定律和水均衡原理达西定律和水均衡原理DA+B+C+D源汇项源汇项2021/6/1614六面体内地下水储存量的变化为六面体内地下水储存量的变化为由水均衡原理得三维地下水流动方程的有限差分格式由水均衡原理得三维地下水流动方程的有限差分格式(2)有限差分方程建立有限差分方程建立(续续)方法二：方法二：达西定律和水均衡原理达西定律和水均衡原理2021/6/1615有限差分法有限差分法:三维（三维（MODFLOW）差商代替微商差商代替微商2021/6/1616(3)差分方程求解差分方程求解l一维显式差
8、分格式网格个数为网格个数为ni直接求解直接求解2021/6/1617(3)差分方程求解差分方程求解l一维隐式差分格式一维隐式差分格式网格个数为网格个数为ni迭代求解迭代求解方方程程组组PCGSIPSORWHSSAMGGMGMODFLOW2021/6/1618(4)差分方程的收敛性和稳定性差分方程的收敛性和稳定性l截断误差截断误差：用差商代替微商时，地下水流动方程产生用差商代替微商时，地下水流动方程产生的误差为截断误差。的误差为截断误差。l收敛性收敛性：当空间步长和时间步长趋于当空间步长和时间步长趋于0 0时，有限差分方时，有限差分方程的精确解趋于地下水流动问题微分方程定解问题的程的精确解趋于地
9、下水流动问题微分方程定解问题的精确解。则称该差分格式是收敛的。精确解。则称该差分格式是收敛的。l稳定性稳定性：如果在求解差分方程过程中，某时间步引入如果在求解差分方程过程中，某时间步引入某个误差，而在以后的各时段计算中，该误差不再扩某个误差，而在以后的各时段计算中，该误差不再扩大，则称该差分格式是稳定的。大，则称该差分格式是稳定的。一维一维显示格式显示格式的收敛条件和稳定条件是：的收敛条件和稳定条件是：2021/6/1619(6)算例：算例：显式有限差格式显式有限差格式 设两条河流平行、完全切割含水层，含水层等厚、均质设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。各向同性。应用实例
10、：河间地块承压水流模型应用实例：河间地块承压水流模型2021/6/1620步骤：步骤：（1 1）基础资料的分析）基础资料的分析（2 2）概念模型）概念模型（3 3）数学模型）数学模型（4 4）数值方法及计算机程序）数值方法及计算机程序（5 5）参数）参数（6 6）结果分析）结果分析 2021/6/1621建立数学模型建立数学模型（1 1）模型概化）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。由所述水文地质条件，可以概化为一维承压水流问题。（2 2）建立坐标系）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内，（如图），将地下水流动系统空间结构放在坐标系内，从而量化各变量的
11、取值范围。本例，取从而量化各变量的取值范围。本例，取x-x-轴原点位于左端河，右侧轴原点位于左端河，右侧为正向，设两河流间距为为正向，设两河流间距为L. L. （3 3）数学模型）数学模型2021/6/1622差分方程及其解法差分方程及其解法显式格式显式格式将（0L）分成 N 等份，1)1)网格剖分：网格剖分：取时间步长取时间步长 ，记，记 （n n=0=0、1 1、2 2、3 3、44）记 ，(i=0,1,2,3，4N)2 2）建立差分方程）建立差分方程：在网格系统中任意取一点：在网格系统中任意取一点设设是问题的解，则在是问题的解，则在处有处有记为（记为（i i，n n）2021/6/162
12、3用差商代替微商用差商代替微商： 将上述两式舍去余项，代入方程并记将上述两式舍去余项，代入方程并记为为显然该式具有截断误差得到得到显式格式（续显式格式（续1）2021/6/1624引入无量纲变量引入无量纲变量：将该式子代入得到：将该式子代入得到： （i=1，2，3，.N-1），（n=1，2，3，.） 显式格式（续显式格式（续2）2021/6/16253）显示差分方程的求解）显示差分方程的求解l计算各结点初始时刻水头值计算各结点初始时刻水头值l利用差分方程计算各结点利用差分方程计算各结点t1t1时刻水头值时刻水头值l利用边界条件计算边界结点水头值利用边界条件计算边界结点水头值l重复重复2 2、3
13、 3步，直到计算出拟计算的各个时刻的水头步，直到计算出拟计算的各个时刻的水头值值显式格式（续显式格式（续3）2021/6/1626算例（续算例（续4）在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200200米，时间步长为米，时间步长为0.250.25天，分别计算各节点天，分别计算各节点各各时刻的时刻的水头值。水头值。2021/6/1627Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 m02010101010100.252012.5101010100.502013.7510.6251010100.752014.531
14、11.250 10.156 10101.002015.078 11.797 10.391 10.039 101.252015.488 12.266 10.654 10.117 10算例（续算例（续5）2021/6/1628Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 m02010101010100.252020101010100.502010201010100.7520300 2010101.0020-10 101.252010算例（续算例（续6）如果如果 t=1, t=1,则则2021/6/1629(6)算例：算例：隐式格式隐式格式在上述模型
15、中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200米，时间步长为米，时间步长为0.25天，用隐式差分格式计算各节天，用隐式差分格式计算各节点个时刻的水头值。点个时刻的水头值。2021/6/1630在这个例子中，在这个例子中，解：隐式格式一般方程为解：隐式格式一般方程为于是有于是有根据初始条件得根据初始条件得根据边界条件得根据边界条件得2021/6/1631由初始条件和边界条件由初始条件和边界条件由此解得由此解得t1t1时刻的水头值为时刻的水头值为在上述方程中取在上述方程中取 n=0, n=0,可以得到计算可以得到计算t1t1时刻水头值的方程时刻水头值的方程所以上述方程变成所以
16、上述方程变成同理，可计算同理，可计算t2t2时刻的水头值时刻的水头值2021/6/16322.地下水数值模型建模步骤地下水数值模型建模步骤2021/6/1633模拟步骤模拟步骤o建立概念模型建立概念模型 o建立数学模型建立数学模型 o数值方法及软件数值方法及软件（编程）（编程） o参数参数 o拟合模拟：模型校拟合模拟：模型校正与检验正与检验 o参数敏感性分析参数敏感性分析 o预测模拟预测模拟软件软件2021/6/1634一、概念模型（模型概化）一、概念模型（模型概化） 根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、岩石矿
17、物、水文、气象、工农业利用情况等岩石矿物、水文、气象、工农业利用情况等模拟的区域模拟的区域: :含水层类型含水层类型: : 潜水（无压）、承压、混合、多层潜水（无压）、承压、混合、多层维数维数: :一维、二维、三维一维、二维、三维水流状态：稳定流水流状态：稳定流/ /非稳定流、饱和流非稳定流、饱和流/ /非饱和流非饱和流介质状况介质状况: :均质和非均质均质和非均质/ /各向同性和各向异性各向同性和各向异性孔隙孔隙/ /裂隙裂隙/ /双重介质双重介质流体的密度差流体的密度差边界条件和初始条件边界条件和初始条件 必要时需进行一系列的室内试验与野外试验必要时需进行一系列的室内试验与野外试验, ,
18、以获取有关参数以获取有关参数, , 如渗透系如渗透系数、弥散系数、分配系数、反应速率常数等。数、弥散系数、分配系数、反应速率常数等。 2021/6/1635收集研究区已有收集研究区已有水文地质资料水文地质资料概念模型概念模型边界概化边界概化内部结构概化内部结构概化完成模型概化图完成模型概化图地表水体地表水体水头边界水头边界流量边界流量边界断层接触边界断层接触边界隔水边界隔水边界水头边界水头边界流量边界流量边界岩体岩层接触边界岩体岩层接触边界天然分水岭天然分水岭据含水层组类型、据含水层组类型、结构、岩性结构、岩性含水介质含水介质地下水运动状态地下水运动状态水文地质参数水文地质参数确定层组的均质、
19、非均质，确定层组的均质、非均质，各向同行、异性，稳定、非各向同行、异性，稳定、非稳定，潜水、承压水稳定，潜水、承压水空间分布空间分布导水、储水系数、主渗透方向导水、储水系数、主渗透方向目的层与相邻层关系目的层与相邻层关系均质、非均质，各向同、异性均质、非均质，各向同、异性层、紊流，二、三维层、紊流，二、三维时间概化时间概化渗透、储水系数，给水度及单位涌渗透、储水系数，给水度及单位涌水量，含水层分布规律，地下水流水量，含水层分布规律，地下水流场，水化学场、温度场的空间概化场，水化学场、温度场的空间概化平面图平面图剖面图剖面图研究范围、主要居民点、标志性地形、地貌、地表河流、研究范围、主要居民点、
20、标志性地形、地貌、地表河流、湖泊、开采井、地下水天然露头、含水层的各类钻孔湖泊、开采井、地下水天然露头、含水层的各类钻孔地表地理要素、含水层结构、地质结构、地下水水位、地表地理要素、含水层结构、地质结构、地下水水位、各类源汇项及其性质各类源汇项及其性质2021/6/1636二、数学模型二、数学模型 三维地下水流动问题控制方程三维地下水流动问题控制方程第二类边界条件第二类边界条件第一类边界条件第一类边界条件初始条件初始条件2021/6/1637 绝大部分数学模型是无法用解析法求解的，数绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。值化就是将数学模型转化为可解的数
21、值模型。三、数值方法及软件（或编程）三、数值方法及软件（或编程）有限差分法有限差分法有限单元法有限单元法积分有限差分法积分有限差分法半解析半数值法半解析半数值法边界元法边界元法2021/6/1638有限差分法有限差分法: MOFLOW系列系列GMS 中中MODFLOWVisual MODFLOWProcessing MODFLOW有限单元法有限单元法: FEFLOW积分有限差分法积分有限差分法: TOUGH2，TOUGH REACT软件软件2021/6/1639四、模型参数四、模型参数o含水层参数：渗透系数，弹性释水系数（重力给水含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等度），
22、孔隙度等o源汇项：源汇项：n大气降水入渗系数（分区、数值）大气降水入渗系数（分区、数值）n蒸发排泄系数蒸发排泄系数n地表水体（河流、湖泊、水库等）水位、底面高程、底地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等）面岩性特征（厚度、渗透系数等）n渠系灌溉入渗系数渠系灌溉入渗系数n人工开采（点状、面状）人工开采（点状、面状）o边界条件边界条件o初始条件初始条件2021/6/1640参数的不确定性参数的不确定性钻孔太少，钻孔太少，地层资料少，地层资料少，钻孔多，含钻孔多，含水层结构会水层结构会发生变化发生变化2021/6/1641五、模拟：模型校正（参数识别五、模拟：模
23、型校正（参数识别） 将模拟结果与实测结果比较将模拟结果与实测结果比较,进行参数调整进行参数调整, 使模使模拟结果在给定的误差范围内与实测结果吻合。拟结果在给定的误差范围内与实测结果吻合。 调参过程是一个复杂而辛苦的工作调参过程是一个复杂而辛苦的工作, 所调整的参数所调整的参数必须符合模拟区的具体情况。必须符合模拟区的具体情况。 人机交互与自动调参相结合。尽管自动调参程序人机交互与自动调参相结合。尽管自动调参程序( (如如PEST ) , PEST ) , 也不能代替人的工作。也不能代替人的工作。2021/6/1642五、模拟：模型检验五、模拟：模型检验 模型验证是在模型校正的基础上模型验证是在
24、模型校正的基础上, 进一步进一步调整参数调整参数, 使模拟结果与第二次实测结果吻合使模拟结果与第二次实测结果吻合, 以进一步提高模型的置信度。以进一步提高模型的置信度。2021/6/1643六、灵敏度分析六、灵敏度分析 校正后的模型受参数值的时空分布、边界校正后的模型受参数值的时空分布、边界条件、水流状态等条件、水流状态等不确定度不确定度的影响。灵敏度的影响。灵敏度分析就是为了确定不确定度对校正模型的影分析就是为了确定不确定度对校正模型的影响程度。响程度。2021/6/1644七、预测七、预测 用校正的参数值进行预测用校正的参数值进行预测, 预测时需估算未预测时需估算未来的水流状态。来的水流状
25、态。2021/6/1645后续检查与模型的重建（完善）后续检查与模型的重建（完善） 后续检查在模拟研究结束数年后进行。收集新的野后续检查在模拟研究结束数年后进行。收集新的野外数据以确定预测结果是否正确。如果模拟结果精确外数据以确定预测结果是否正确。如果模拟结果精确, , 则该模型对该模拟区来说是有效的。由于场址的唯一则该模型对该模拟区来说是有效的。由于场址的唯一性性, , 故模型只对该模拟区有效。故模型只对该模拟区有效。 后续检查应后续检查应在预测结束足够长的时间后进行在预测结束足够长的时间后进行, , 以以便有足够的时间发生明显的变化。便有足够的时间发生明显的变化。 2021/6/1646模
26、型的再设计模型的再设计 一般来说一般来说, 后续检查会发现系统性能的变后续检查会发现系统性能的变化化, 从而导致概念模型和模型参数的修改。一从而导致概念模型和模型参数的修改。一般来说般来说, 所有模拟研究都应该进行到第五步所有模拟研究都应该进行到第五步, 即校正灵敏度分析。即校正灵敏度分析。 2021/6/1647地地裂裂缝缝自然地理及水文地质条件自然地理及水文地质条件边边界界敏敏感感性性分分析析数学模型数学模型网网格格化化概念模型概念模型介质类型、结构特征介质类型、结构特征地下水补、径、排特征地下水补、径、排特征数值模型数值模型数值模型的建立数值模型的建立边边界界条条件件初初始始条条件件源源
27、汇汇项项降降雨雨入入渗渗地地表表水水体体入入渗渗灌灌溉溉、渠渠系系蒸蒸发发排排泄泄人人工工开开采采介介质质参参数数观观测测孔孔动动态态拟拟合合流流场场拟拟合合水水均均衡衡对对比比分分析析观观测测孔孔动动态态检检验验流流场场检检验验水水均均衡衡对对比比检检验验拟合调参拟合调参识识别别结结果果不不符符合合要要求求修修正正概概念念模模型型小结：地下水流动小结：地下水流动模型构建过程模型构建过程识识别别结结果果不不符符合合要要求求修修正正概概念念模模型型模型应用模型应用水水资资源源量量评评价价预预测测开开采采方方案案研研究究环环境境地地质质问问题题地地面面沉沉降降地地下下水水污污染染研研究究大大型型工
28、工程程对对环环境境的的影影响响模型识别模型识别模型检验模型检验2021/6/16483.地下水流数值模拟资料需求地下水流数值模拟资料需求o空间展布相关资料空间展布相关资料n地表高程、分层数据地表高程、分层数据n边界位置边界位置o含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等o源汇项：源汇项：n大气降水入渗系数（分区、数值）大气降水入渗系数（分区、数值）n蒸发排泄系数蒸发排泄系数n地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等）数等）n渠系灌溉入渗系数渠系灌溉入渗系数n人工开采（点状、面状）：开采井位置、井结构、开采量动态人工开采（点状、面状）：开采井位置、井结构、开采量动态o边界条件：边界类型、水头或流量边界条件：边界类型、水头或流量o初始条件：统测水位初始条件：统测水位o水位动态观测资料：观测孔位置、结构、水位时间变化水位动态观测资料：观测孔位置、结构、水位时间变化2021/6/1649 结束语结束语若有不当之处，请指正，谢谢！若有不当之处，请指正，谢谢！
