第五章 地下水污染预测

第五章地下水污染预测

主要内容概述（了解）

近似解析法（了解）

数值法（熟悉）

数理统计法（熟悉）灰色预测法（了解）

国外地下水模型软件2/5/20231第一节概述

地下水污染预测预报包括：预报污染地下水分布边界的推进情况（即经过一定时间t后边界推进的距离x，或预测推进一定距离需要多长时间）；确定污染地下水中污染物质的浓度在空间的分布及随时间的变化情况；预测污染地下水能否侵入附近的水源地，污染物质达到水源地所需的时间；预测防治地下水污染的措施的效果等。2/5/20232第一节概述

预测评价的基本步骤如下：建立预测模型；模型识别（平衡性、稳定性、敏感性）；模型参数估计；模型方程的验证；地下水污染预测。2/5/20233第二节近似解析解前面介绍了一些解析解公式，但应用起来仍十分复杂，不便于实际应用。实际中可以进一步简化：完全忽略弥散作用，将污染地下水的运动视为“活塞式”的推挤淡水的运动，假设两者始终保持明显的铅直分界面；将地下水动力运移和弥散分开考虑，先以上述方法计算，然后计算由于弥散所形成的过渡混合带，修正按平均锋面运移的距离。2/5/20234第二节近似解析解

计算平均渗透锋面的运移距离；求出平均渗透锋面处污染物的浓度；确定分界面弥散带和变形带，在平均渗透锋面距离上加上这两个带就可确定污染水实际锋面的运移位置：

相对密度差岩层倾角2/5/20235第二节近似解析解承压（或潜水）含水层污染地下水在承压含水层的运动，如果补给来源稳定，可视为似稳定运动，属于一维平面平行流，则渗透速度为：而实际流速为：地下水的天然单宽流量2/5/20236第二节近似解析解承压（或潜水）含水层初始时刻污染锋面的位置某时刻污染锋面的位置天然地下水中污染物的浓度初始分界面以内，地下水中污染物的浓度若是潜水，则取平均厚度代替上述各式中的m。2/5/20237第三节数值解由于含水层以及污染物质迁移的复杂性，解析解不易求得，通常采用数值解法。数值法在剖分后认为每一单元内均质同性，将方程线性化。有限差分法有限单元法边界元法2/5/20238一、有限差分法有限差分法的基本思路是按照时间步长和空间步长将时间和空间区域剖分成若干网格，用未知函数在网格节点上的值所构成的差分商近似代替所用偏微分方程中出现的各阶导数，从而把表示变量连续变化关系的偏微分方程离散为有限个代数方程，然后解此线性代数方程，以求出污染物浓度在各网格节点上不同时刻的解。

2/5/20239导数的差分近似设f(x)为任一足够光滑的函数，把f(x)沿x的正向展开为泰勒级数：上式称为f(x)的一阶向前差分2/5/202310导数的差分近似如果把f(x)对x的负向展开为泰勒级数：上式称为f(x)的一阶向后差分2/5/202311导数的差分近似把上面两个差分方程相减，可得：上式称为f(x)的一阶中心差分同理可得f(x)的二阶导数的差分公式：2/5/202312导数的差分近似若把上述一元函数f(x)换为浓度函数，并按网格节点编号，则对x的一阶向前、向后和中心差分分别为：同样可以写出对时间t的向前、向后和中心差分公式2/5/202313导数的差分近似C(x,y,z,t)对x的二阶差分公式为：因此，在任何一个节点处，浓度C(x,y,z,t)的一阶和二阶偏导数均可用该节点及其相邻节点上浓度值的线性组合表示。2/5/202314一维弥散方程的差分格式设在无限含水层中，存在一维均匀流场，其渗透速度V=nu，流动方向为x轴正向。则一维弥散方程为：2/5/202315一维弥散方程的差分格式对时间区域[0,T]和空间区域[0,L]都作等距剖分，设时间步长为△t，空间步长为△x，把第i个节点xi处在tn时刻的浓度记为Cin。2/5/202316一维弥散方程的差分格式显示差分格式将及差分公式中的浓度取为tn时刻的值，便可得到其显示差分格式：2/5/202317一维弥散方程的差分格式显示差分格式的稳定条件在应用显示差分格式时，剖分的步长必须满足：2/5/202318一维弥散方程的差分格式隐示差分格式将及差分公式中的浓度取为tn＋1时刻的值，便可得到其隐示差分格式：2/5/202319一维弥散方程的差分格式隐示差分格式的稳定条件在应用隐示差分格式时，剖分的步长必须满足：2/5/202320一维弥散方程的差分格式隐示差分格式的矩阵形式2/5/202321一维弥散方程的差分格式Crank-Nicolson差分格式取显示差分格式和隐示差分格式的平均，便可得到C-N差分格式：2/5/202322一维弥散方程的差分格式Crank-Nicolson差分格式整理后可得：2/5/202323追赶法：方程的统一形式：2/5/202324追赶法：引入边界条件：当i=1时：当i=M时：2/5/202325追赶法：令A·C=δ2/5/202326追赶法：追赶法解方程组的方法：先把对角矩阵A分解为：A=B·W其中：2/5/202327追赶法：由于A=B·W2/5/202328追赶法：由此可得：

b1=β1

b1w1=γ1→

w1=γ1/b1

b2+a2w1=β2

biwi=γi

→

wi=γi/bi

(i=1……M-1)

bi+aiwi-1=βi

→

bi=βi-aiwi-1

bM=βM’-aM’wM-12/5/202329追赶法：由于A=B·W

→

B·(W·C)=δ又令：W·C=Y其中：Y=(y1,y2,…,yM)T→B·Y=δ2/5/202330追赶法：由以上等式可知：

b1y1=δ1’→y1=δ1’/b1

a2y1+b2y2=δ2→

y2=(δ2-a2y1)/b2aiyi-1+biyi=δi→

yi=(δi-aiyi-1)/biaM’yM-1+bMyM=δM→

yM=(δM-aM’yM-1)/bM至此可以求出所有的Y(y1,y2,…,yM)T值。2/5/202331追赶法：又因为W·C=Y，即：2/5/202332追赶法：根据上述等式，可得：2/5/202333追赶法：追：根据B·Y=δ，求出bi和wi，并求出yi；赶：根据W·C=Y，求出Cin+1。2/5/202334二维弥散方程的差分格式若多孔介质是各向同性的，流场为均匀的一维流场，渗透速度V=nu，流动方向取x轴的正向，与流速垂直为y方向，对流扩散方程可简化为：2/5/202335二维弥散方程的差分格式将研究区域网格化2/5/202336二维弥散方程的差分格式建立矩形网差分格式隐式差分格式2/5/202337二维弥散方程的差分格式于是，内节点的差分方程为2/5/202338二维弥散方程的差分格式整理后可得2/5/202339二维弥散方程的差分格式边界条件初始条件2/5/202340二维弥散方程的差分格式将网格取为正方形，即△x=△y=h，方程就可简化为：2/5/202341二维弥散方程的差分格式交替方向差分格式（ADI）ADI法的实质是将二维隐式差分格式化为一维的隐式差分格式。在第n层与第n+1层两层正中间引入一个过渡层n+1/2。2/5/202342二维弥散方程的差分格式交替方向差分格式（ADI）

分两步走：在前半时间步长内，对x偏导用隐式表示，对y的偏导用显示表示；在后半时间步长内，对x偏导用显示表示，对y的偏导用隐式表示。2/5/202343二维弥散方程的差分格式交替方向差分格式（ADI）前半时间步长：

2/5/202344二维弥散方程的差分格式交替方向差分格式（ADI）后半时间步长：

2/5/202345二维弥散方程的差分格式计算步骤：（1）根据初始条件或前一时段计算结果，可知tn时刻的浓度值，那么第一个方程中只有三个未知数（n+1/2时刻）。如果固定j值，对i的联立方程组就是三对角线的（追赶法），即沿x方向追赶，可求出过渡层（n+1/2）层上的全部浓度值。

2/5/202346二维弥散方程的差分格式计算步骤：（2）根据前一步求出的过渡层（n+1/2）层上的浓度值，再按同样的方法求出第二个方程中的（n+1）时刻的浓度值。即对于一个固定的i值，方程也是一个三对角方程组（追赶法），沿y方向追赶。

2/5/202347有限单元法基本原理在有限差分法中，研究区域被一系列的网格节点所代替，在各个节点上以差分代替导数而得到求解相应偏微分方程的线性代数方程组。有限单元法是把研究区域剖分为有限个子区域，在每个子区域上用某种插值函数来近似待求解的未知函数，从而得到求解相应偏微分方程的线性代数方程组。

2/5/202348边界元法基本原理通过把求解区域的边界剖分为若干单元，将求函数的解简化为求节点的函数值，求解积分方程就化为求解一组线性代数方程，在求出边界上的物理量以后，计算区域内部的任一点的未知量可通过边界上的已知量求出。与有限元法相比，具有降低维数、输入数据准备简单、计算工作量小、精度比有限元法高等优点。在处理非均质、非线性、非稳定等复杂问题时计算效率低，不能将对流项化成沿边界的积分。2/5/202349第四节数理统计法地下水水质受多种因素综合影响，具有不确定性，可用数理统计法预测水质变化；主要包括：时间序列分析法、相关分析法、概率统计法。2/5/202350（1）时间序列分析法基本思路：分析时间序列的变化特征，选择适当的模式建立预测模型，利用模型进行趋势外推预测，对模型预测值进行评价和修正，得到预测结果。2/5/202351（1）时间序列分析法-平均数平均数预测：对动态数列中若干个时期的数值，求其算术平均值或加权平均值，作为下一个时期的预测值。2/5/202352（1）时间序列分析法-指数平滑一次指数平滑预测：对动态数列中若干个时期的数值，求其算术平均值或加权平均值，作为下一个时期的预测值。2/5/202353（1）时间序列分析法-指数平滑一次指数平滑法是以最近周期的一次指数平滑值作为下一周期的预测值，即：其中平滑系数a的选择直接影响预测效果。如果时间序列长期趋势稳定，应取较小的a值；若时间序列不稳定，应取较大的a值。2/5/202354（1）时间序列分析法-指数平滑a越小，对数据的平滑能力越强，但对数据变化的敏感性越差；a越大，对数据的平滑能力越差，但对数据变化的敏感性越强。2/5/202355（1）时间序列分析法-指数平滑初始值S0(1)对预测结果的影响很大，一般可取第一个数据（x1）作为初始值；若总体数据量较少，则可取最初几个数据的平均值作为初始值。2/5/202356（1）时间序列分析法-指数平滑二次指数平滑预测：对一次指数平滑值再做一次指数平滑：2/5/202357（1）时间序列分析法-指数平滑二次指数平滑值不直接用于预测，而是根据滞后偏差建立线性预测模型：其中：2/5/202358（1）时间序列分析法-指数平滑若一、二次指数平滑值不满足预测要求，应采用三次指数平滑法建立非线性预测模型：2/5/202359（1）时间序列分析法-指数平滑三次指数平滑法建立的非线性预测模型：2/5/202360（2）相关分析法线性相关预测：建立直角坐标下数据散点图，用一元线性回归方程近似描述数据点的分布规律：2/5/202361（2）相关分析法根据最小二乘法(要求残差RSS最小)可求得：用相关系数r来判断回归方程实用性：2/5/202362（2）相关分析法曲性相关预测：将直角坐标下散点图拟合形成的有规律的曲线进行函数转换，变为前述的一元线性回归方程进行计算。如：2/5/202363国外地下水模拟软件的发展现状与趋势

利用数值模型对地下水流和溶质运移问题进行模拟的方法以其有效性、灵活性和相对廉价性逐渐成为地下水研究领域的一种不可或缺的重要方法，并受到越来越大的重视和广泛的应用。一个完整的地下水模拟过程包含3个部分：前处理、模型计算和后处理。2/5/202364国外地下水模拟软件的发展现状与趋势前处理是指在进行模拟计算之前对计算过程中所需数据的整理、组织、输入及计算网格的编号与生成。模型计算是进行地下水流动或水质运移正反演计算，常用的方法主要有：有限差分法、有限元法、边界元法等。后处理是将计算所产生的结果数据，用图形或表格显示或存放起来，以供研究人员方便地进行分析和使用。2/5/202365国外地下水模拟软件的发展现状与趋势通过近二十年的研究与发展，国际上已经形成了一批非常有影响的地下水模拟的软件，它们今天在国际地下水模拟研究领域依旧非常活跃，如MODFLOW、FEFLOW、MT3DMS、MT3D99、PEST、MODPATH、UCODE等。2/5/202366MODFLOW

MODFLOW是由美国地质调查局的McDonald和Harbaugh于80年代开发出来的一套专门用于孔隙介质中三维有限差分地下水流数值模拟的软件。程序结构的模块化（模拟抽水引起地面沉降的子程序包、模拟水平流动障碍的子程序包等）。离散方法的简单化（有限差分法360×360×18个网格单元）。MODFLOW引进了应力期(StressPeriod)概念。求解方法的多样化（强隐式法、逐次超松弛迭代法、预调共轭梯度法）。2/5/202367MT3D99

MT3D99是美国S.S.Papodopulos&Associates公司的Zheng博士设计的模拟三维地下水溶质运移程序MT3D(1990)的增强版，成为目前世界上首屈一指的溶质运移模拟软件。MT3D99能够模拟地下水系统中的平流、扩散、衰减、溶质化学反应、线性与非线性吸附作用等现象，能够对承压含水层，不承压含水层，承压与不承压交替的含水层以及倾斜的和单元厚度变化的含水层进行空间离散。

MT3D99提供了丰富的求解方法（标准有限差分法、基于Eulerian-Lagrangian的粒子跟踪方法和高阶有限体积TVD方法）。

2/5/202368PEST2000PEST2000是由澳大利亚WatermarkComputing公司开发的PEST新版本，是功能强大的模型独立的参数估计程序。PEST2000利用一个强有力的数值反演算法来“控制”运行中的模型，程序在每次模拟之后自动调整所选择的模型参数，直到将校正的目标最小化为止。PEST2000引进了参数估计的最新技术──预测分析。

2/5/202369FEFLOW

FEFLOW是加拿大Waterloo水文地质公司开发的基于三维有限元的地下水模拟可视化软件包。它能够解决下列地下水模拟问题：完全瞬时、半瞬时、稳态地下水流动与溶质运移；随时间变化的实体属性和约束边界条件；饱和与不饱和流动；包含潜水面的承压与不承压含水层；带有非线性吸附作用、衰变、对流、弥散的化学质量运移；考虑贮存、对流、热散失、热运移的流体和固体热量运移；密度变化的流动(海水入侵等)。2/5/202370VisualModflowGraphicalUserInterface

2/5/202371VisualModflowPumpingWells

2/5/202372VisualModflowFlowandTransportProperties

Kx,Ky,Kz,Ss,Sy,n,h0Dx,Dy,Dz,ρ,Kd,K,andC0

2/5/202373VisualModflowBoundary

－s