地下水建模方法和步骤.ppt

地下水建模方法和步骤,中国地质大学(武汉)环境学院 2012.8,地下水建模方法和步骤,1.求解地下水运动方程的数值方法 2.地下水数值模型建模步骤 3.建模所需要的基本资料,绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。,1.数值方法,1.1有限差分法 1.2有限单元法 1.3分有限差分法 1.4半解析半数值法 1.5边界元法,(1)有限差分法原理 (2)两种方法建立有限差分方程 (3)求解有限差分方程 (4)收敛性和稳定性概念 (5)算例,1.1有限差分法,(1)有限差分法的基本原理,将连续的问题离散后求解： 方法一以地下水流基本微分方程及其定解条件为基础， 在渗流区剖分基础上，用差商代替微商，将地下水流微分方程的求解转化为差分方程（代数方程）求解。 方法二在渗流区剖分的基础上，直接由达西定律和水均衡原理，建立各个均衡区的水均衡方程，即差分方程。,矩形网格,多边形网格,网格划分的基本类型,（1）先划格线，格点位于网格中心,均衡网格,节点网格,（2）先规定格点位置，再垂直平分两相邻结点的连线作格线，形成的网格即为水均衡区,MODFLOW网格系统,导数的有限差商近似,导数的定义,当,非常小的时候，有,上式右端项即为f(x)在x0处的差商。,这样定义的差商很容易理解，但不知道用差商代替微商所产生的误差。下面利用泰勒公式导出差商及其误差。,方法一：差商代替微商,(2)有限差分方程建立,已知泰勒公式, 由A得：,A,B, 由B 得：,称 为f(x)在x0处的一阶前向差商， 为截断误差。,称 为f(x)在x0处的一阶后向差商， 为截断误差。,方法一,由A-B可以得：,由A+B可以得：,A,B,称 为f(x)在x0处的一阶中心差商， 为截断误差。,称 为f(x)在x0处的二阶中心差商， 为截断误差。,方法一,对于偏导数（偏微商），类似可以得到相应的差商：,方法一,(2)有限差分方程建立(续),一维控制方程差分格式,显式差分格式,隐式差分格式,方法一,控制方程,网格剖分nx个,取右图所示得微小六面体。设与x,y,z,方向对应得主渗透系数分别为Kx, Ky,Kz；建立均衡期t时段内，微小均衡六面体的水量守恒方程。,方法二：达西定律和水均衡原理,(2)有限差分方程建立(续),基于达西定律，x，y，z方向流入流出分别为:,t时段内，侧向流入与源汇项导致六面体水量变化量为：,A,B,C,(2)有限差分方程建立(续),方法二：达西定律和水均衡原理,D,A+B+C+D,源汇项,六面体内地下水储存量的变化为,由水均衡原理得三维地下水流动方程的有限差分格式,(2)有限差分方程建立(续),方法二：达西定律和水均衡原理,有限差分法:三维（MODFLOW）,差商代替微商,(3)差分方程求解,一维显式差分格式,网格个数为ni,直接求解,(3)差分方程求解,一维隐式差分格式,网格个数为ni,迭代求解,方程组,PCG SIP SOR WHS SAMG GMG,MODFLOW,(4)差分方程的收敛性和稳定性,截断误差：用差商代替微商时，地下水流动方程产生的误差为截断误差。 收敛性：当空间步长和时间步长趋于0时，有限差分方程的精确解趋于地下水流动问题微分方程定解问题的精确解。则称该差分格式是收敛的。 稳定性：如果在求解差分方程过程中，某时间步引入某个误差，而在以后的各时段计算中，该误差不再扩大，则称该差分格式是稳定的。,一维显示格式的收敛条件和稳定条件是：,(6)算例：显式有限差格式,设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。,应用实例：河间地块承压水流模型,步骤：,（1）基础资料的分析 （2）概念模型 （3）数学模型 （4）数值方法及计算机程序 （5）参数 （6）结果分析,建立数学模型,（1）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。 （2）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内，从而量化各变量的取值范围。本例，取x-轴原点位于左端河，右侧为正向，设两河流间距为L. （3）数学模型,差分方程及其解法显式格式,将（0L）分成 N 等份，,1)网格剖分：,取时间步长 ，记 （n=0、1、2、3、4）,记 ，(i=0,1,2,3，4N),2）建立差分方程：在网格系统中任意取一点,设,是问题的解，则在,处有,记为（i，n）,用差商代替微商：,将上述两式舍去余项，代入方程并记,为,显然该式具有截断误差,得到,显式格式（续1）,引入无量纲变量：,将该式子代入得到：,（i=1，2，3，.N-1），（n=1，2，3，.）,显式格式（续2）,3）显示差分方程的求解,计算各结点初始时刻水头值 利用差分方程计算各结点t1时刻水头值 利用边界条件计算边界结点水头值 重复2、3步，直到计算出拟计算的各个时刻的水头值,显式格式（续3）,算例（续4）,在上述模型中，设L=1000米,取空间步长为200米，时间步长为0.25天，分别计算各节点各时刻的水头值。,算例（续5）,算例（续6）,如果 t=1,则,(6)算例：隐式格式,在上述模型中，设L=1000米,取空间步长为200米，时间步长为0.25天，用隐式差分格式计算各节点个时刻的水头值。,在这个例子中，,解：隐式格式一般方程为,于是有,根据初始条件得,根据边界条件得,由初始条件和边界条件,由此解得t1时刻的水头值为,在上述方程中取 n=0,可以得到计算t1时刻水头值的方程,所以上述方程变成,同理，可计算t2时刻的水头值,2.地下水数值模型建模步骤,模拟步骤,建立概念模型 建立数学模型 数值方法及软件（编程） 参数 拟合模拟：模型校正与检验 参数敏感性分析 预测模拟,软件,一、概念模型（模型概化）,根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、岩石矿物、水文、气象、工农业利用情况等 模拟的区域: 含水层类型: 潜水（无压）、承压、混合、多层 维数:一维、二维、三维 水流状态：稳定流/非稳定流、饱和流/非饱和流 介质状况: 均质和非均质/各向同性和各向异性 孔隙/裂隙/双重介质 流体的密度差 边界条件和初始条件 必要时需进行一系列的室内试验与野外试验, 以获取有关参数, 如渗透系数、弥散系数、分配系数、反应速率常数等。,二、数学模型,三维地下水流动问题控制方程,第二类边界条件,第一类边界条件,初始条件,绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。,三、数值方法及软件（或编程）,有限差分法 有限单元法 积分有限差分法 半解析半数值法 边界元法,有限差分法: MOFLOW系列 GMS 中MODFLOW Visual MODFLOW Processing MODFLOW 有限单元法: FEFLOW 积分有限差分法: TOUGH2，TOUGH REACT,软件,四、模型参数,含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等 源汇项： 大气降水入渗系数（分区、数值） 蒸发排泄系数 地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等） 渠系灌溉入渗系数 人工开采（点状、面状） 边界条件 初始条件,参数的不确定性,钻孔太少，地层资料少，钻孔多，含水层结构会发生变化,五、模拟：模型校正（参数识别）,将模拟结果与实测结果比较,进行参数调整, 使模拟结果在给定的误差范围内与实测结果吻合。 调参过程是一个复杂而辛苦的工作, 所调整的参数必须符合模拟区的具体情况。 人机交互与自动调参相结合。尽管自动调参程序(如PEST ) , 也不能代替人的工作。,五、模拟：模型检验,模型验证是在模型校正的基础上, 进一步调整参数, 使模拟结果与第二次实测结果吻合, 以进一步提高模型的置信度。,六、灵敏度分析,校正后的模型受参数值的时空分布、边界条件、水流状态等不确定度的影响。灵敏度分析就是为了确定不确定度对校正模型的影响程度。,七、预测,用校正的参数值进行预测, 预测时需估算未来的水流状态。,后续检查与模型的重建（完善）,后续检查在模拟研究结束数年后进行。收集新的野外数据以确定预测结果是否正确。如果模拟结果精确, 则该模型对该模拟区来说是有效的。由于场址的唯一性, 故模型只对该模拟区有效。 后续检查应在预测结束足够长的时间后进行, 以便有足够的时间发生明显的变化。,模型的再设计,一般来说, 后续检查会发现系统性能的变化, 从而导致概念模型和模型参数的修改。一般来说, 所有模拟研究都应该进行到第五步, 即校正灵敏度分析。,地裂缝,自然地理及水文地质条件,边界敏感性分析,数学模型,网格化,概念模型,介质类型、结构特征,地下水补、径、排特征,数值模型,数值模型的建立,边界条件,初始条件,源汇项,降雨入渗,地表水体入渗,灌溉、渠系,蒸发排泄,人工开采,介质参数,观测孔动态拟合,流场拟合,水均衡对比分析,观测孔动态检验,流场检验,水均衡对比检验,拟合调参,识别结果不符合要求修正概念模型,小结：地下水流动模型构建过程,识别结果不符合要求修正概念模型,模型应用,水资源量评价,预测开采方案,研究环境地质问题,地面沉降,地下水污染研究,大型工程对环境的影响,模型识别,模型检验,3.地下水流数值模拟资料需求,空间展布相关资料 地表高程、分层数据 边界位置 含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等 源汇项： 大气降水入渗系数（分区、数值） 蒸发排泄系数 地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等） 渠系灌溉入渗系数 人工开采（点状、面状）：开采井位置、井结构、开采量动态 边界条件：边界类型、水头或流量 初始条件：统测水位 水位动态观测资料：观测孔位置、结构、水位时间变化,