第五章 地下水污染预测

1、13/7/20221第五章第五章 地下水污染预测地下水污染预测 主要内容主要内容l 概述（了解）概述（了解）l 近似解析法近似解析法（了解）（了解）l 数值法数值法（熟悉）（熟悉）l 数理统计法数理统计法（熟悉）（熟悉）l 灰色预测法（了解）灰色预测法（了解）l 国外地下水模型软件国外地下水模型软件23/7/20222第一节第一节 概述概述 地下水污染预测预报包括：地下水污染预测预报包括：l 预报污染地下水分布预报污染地下水分布边界的推进情况边界的推进情况（即经过（即经过一定时间一定时间t后边界推进的距离后边界推进的距离x，或预测推进一定，或预测推进一定距离需要多长时间）；距离需要多长时间）；

2、l 确定污染地下水中污染物质的确定污染地下水中污染物质的浓度浓度在空间的分在空间的分布及随时间的变化情况；布及随时间的变化情况；l 预测污染地下水能否侵入附近的预测污染地下水能否侵入附近的水源地水源地，污染，污染物质达到水源地所需的时间；物质达到水源地所需的时间；l 预测防治地下水污染的预测防治地下水污染的措施措施的效果等。的效果等。33/7/20223第一节第一节 概述概述 预测评价的基本步骤如下：预测评价的基本步骤如下：l 建立预测模型；建立预测模型；l 模型识别（平衡性、稳定性、敏感性）；模型识别（平衡性、稳定性、敏感性）；l 模型参数估计；模型参数估计；l 模型方程的验证；模型方程的验

3、证；l 地下水污染预测。地下水污染预测。43/7/20224第二节第二节 近似解析解近似解析解前面介绍了一些解析解公式，但应用起来仍十前面介绍了一些解析解公式，但应用起来仍十分复杂，不便于实际应用。分复杂，不便于实际应用。实际中可以进一步简化：实际中可以进一步简化：l完全忽略弥散作用，将污染地下水的运动视为完全忽略弥散作用，将污染地下水的运动视为“活活塞式塞式”的推挤淡水的运动，假设两者始终保持明显的推挤淡水的运动，假设两者始终保持明显的铅直分界面；的铅直分界面；l将地下水动力运移和弥散分开考虑，先以上述方法将地下水动力运移和弥散分开考虑，先以上述方法计算，然后计算由于弥散所形成的过渡混合带，

4、修计算，然后计算由于弥散所形成的过渡混合带，修正按平均锋面运移的距离。正按平均锋面运移的距离。53/7/20225第二节第二节 近似解析解近似解析解 计算平均渗透锋面的运移距离计算平均渗透锋面的运移距离 ； 求出平均渗透锋面处污染物的浓度求出平均渗透锋面处污染物的浓度 ； 确定分界面弥散带确定分界面弥散带 和变形带和变形带 ，在平，在平均渗透锋面距离上加上这两个带就可确定污染均渗透锋面距离上加上这两个带就可确定污染水实际锋面的运移位置：水实际锋面的运移位置： nKmtLnDtLLLLLpdpd/cos/)(5 . 0LCdLpL相对密度差相对密度差岩层倾角岩层倾角63/7/20226第二节第二

5、节 近似解析解近似解析解l 承压（或潜水）含水层承压（或潜水）含水层 污染地下水在承压含水层的运动，如果补给来源污染地下水在承压含水层的运动，如果补给来源稳定，可视为似稳定运动，属于一维平面平行流，则稳定，可视为似稳定运动，属于一维平面平行流，则渗透速度为：渗透速度为：而实际流速为：而实际流速为：常数mqLHHKVex)(12nVuxx地下水的天然单宽流量地下水的天然单宽流量L73/7/20227第二节第二节 近似解析解近似解析解l承压（或潜水）含水层承压（或潜水）含水层eeeeCCCtxCCqmnLttmnqxxL00),(初始时刻污染锋面的位置初始时刻污染锋面的位置某时刻污染锋面的位置某时

6、刻污染锋面的位置天然地下水中污染物的浓度天然地下水中污染物的浓度初始分界面以内，地下水中污染物初始分界面以内，地下水中污染物的浓度的浓度若是潜水，则取平均厚度代替上述各式中的若是潜水，则取平均厚度代替上述各式中的m。L83/7/20228第三节第三节 数值解数值解由于含水层以及污染物质迁移的复杂性，由于含水层以及污染物质迁移的复杂性，解析解不易求得，通常采用数值解法。数解析解不易求得，通常采用数值解法。数值法在剖分后认为每一单元内均质同性，值法在剖分后认为每一单元内均质同性，将方程线性化。将方程线性化。 有限差分法有限差分法 有限单元法有限单元法 边界元法边界元法93/7/20229一、一、

7、有限差分法有限差分法有限差分法的有限差分法的基本思路基本思路是按照时间步长和空间步是按照时间步长和空间步长将时间和空间区域剖分成若干网格，用未知函长将时间和空间区域剖分成若干网格，用未知函数在网格节点上的值所构成的差分商数在网格节点上的值所构成的差分商近似代替近似代替所所用偏微分方程中出现的各阶导数，从而把表示变用偏微分方程中出现的各阶导数，从而把表示变量连续变化关系的偏微分方程离散为有限个代数量连续变化关系的偏微分方程离散为有限个代数方程，然后解此方程，然后解此线性代数方程线性代数方程，以求出污染物浓，以求出污染物浓度在各度在各网格节点网格节点上不同时刻的解。上不同时刻的解。 103/7/2

8、02210导数的差分近似导数的差分近似设设f(x)为任一足够光滑的函数，把为任一足够光滑的函数，把f(x)沿沿x的正向的正向展开为泰勒级数：展开为泰勒级数：上式称为上式称为f(x)的的一阶向前差分一阶向前差分 333222! 3! 2)()(dxfdxdxfdxdxdfxxfxxf)()()(xRxxfxxfdxdfxxfxxfdxdf)()(113/7/202211导数的差分近似导数的差分近似如果把如果把f(x)对对x的负向展开为泰勒级数：的负向展开为泰勒级数：上式称为上式称为f(x)的的一阶向后差分一阶向后差分 333222! 3! 2)()(dxfdxdxfdxdxdfxxfxxf)()

9、()(xOxxxfxfdxdfxxxfxfdxdf)()(123/7/202212导数的差分近似导数的差分近似把上面两个差分方程相减，可得：把上面两个差分方程相减，可得：上式称为上式称为f(x)的的一阶中心差分一阶中心差分同理可得同理可得f(x)的二阶导数的差分公式：的二阶导数的差分公式：xxxfxxfdxdf2)()(222)()(2)(xxxfxfxxfdxfd133/7/202213导数的差分近似导数的差分近似若把上述一元函数若把上述一元函数f(x)换为浓度函数，并按网格换为浓度函数，并按网格节点编号，则对节点编号，则对x的一阶向前、向后和中心差分的一阶向前、向后和中心差分分别为：分别为

10、：xCCxcxCCxcxCCxcnkjinkjinkjinkjinkjinkjinkjinkjinkji2|,1,1),(,1,),(,1),(同样可以写出同样可以写出对时间对时间t的向的向前、向后和中前、向后和中心差分公式心差分公式143/7/202214导数的差分近似导数的差分近似C(x,y,z,t)对对x的二阶差分公式为：的二阶差分公式为：2,1,1),(222|xCCCxcnkjinkjinkjinkji因此，在任何一个节点处，浓度因此，在任何一个节点处，浓度C(x,y,z,t)的一阶和二阶偏导数均可用该节点及其的一阶和二阶偏导数均可用该节点及其相邻节点上浓度值的线性组合表示。相邻节点

11、上浓度值的线性组合表示。153/7/202215一维弥散方程的差分格式一维弥散方程的差分格式设在无限含水层中，存在一维均匀流场，其渗透设在无限含水层中，存在一维均匀流场，其渗透速度速度V=nu，流动方向为，流动方向为x轴正向。则一维弥散方轴正向。则一维弥散方程为：程为： Tt,0 Lx0 22xcuxcDtcL163/7/202216一维弥散方程的差分格式一维弥散方程的差分格式对时间区域对时间区域0,T和空间区域和空间区域0,L都作等距剖分，都作等距剖分，设时间步长为设时间步长为t，空间步长为，空间步长为x，把第，把第i个节点个节点xi处在处在tn时刻的浓度记为时刻的浓度记为Cin。173/7

12、/202217一维弥散方程的差分格式一维弥散方程的差分格式显示差分格式显示差分格式 将将 及及 差分公式中的浓度取为差分公式中的浓度取为tn时刻的值，时刻的值，便可得到其显示差分格式：便可得到其显示差分格式：xC22xCxCCuxCCCDtCCnininininiLnini22112111niLniLniLniCxtuxtDCxtDCxtuxtDC122121)2()21()2(183/7/202218一维弥散方程的差分格式一维弥散方程的差分格式显示差分格式的稳定条件显示差分格式的稳定条件 在应用显示差分格式时，剖分的步长必须满足：在应用显示差分格式时，剖分的步长必须满足：1212xtuxtD

13、L及193/7/202219一维弥散方程的差分格式一维弥散方程的差分格式隐示差分格式隐示差分格式 将将 及及 差分公式中的浓度取为差分公式中的浓度取为tn1时刻的时刻的值，便可得到其隐示差分格式：值，便可得到其隐示差分格式：xC22xCxCCuxCCCDtCCnininininiLnini2211112111111niniLniLniLCCxtuxtDCxtDCxtuxtD11212112)2()21()2(203/7/202220一维弥散方程的差分格式一维弥散方程的差分格式隐示差分格式的稳定条件隐示差分格式的稳定条件 在应用隐示差分格式时，剖分的步长必须满足：在应用隐示差分格式时，剖分的步长

14、必须满足：xDuL2|213/7/202221一维弥散方程的差分格式一维弥散方程的差分格式隐示差分格式的矩阵形式隐示差分格式的矩阵形式1nCA223/7/202222一维弥散方程的差分格式一维弥散方程的差分格式Crank-Nicolson差分格式差分格式 取显示差分格式和隐示差分格式的平均，便可得取显示差分格式和隐示差分格式的平均，便可得到到C-N差分格式：差分格式：)22(21)22(211111211111112111xCCuxCCCDxCCuxCCCDtCCnininininiLnininininiLnini233/7/202223一维弥散方程的差分格式一维弥散方程的差分格式Crank-

15、Nicolson差分格式差分格式 整理后可得：整理后可得：niLniLniLniLniLniLCxtuxtDCxtDCxtuxtDCxtuxtDCxtDCxtuxtD1221211212112)2()1 ( 2)2()2()1 ( 2)2(1nAC243/7/202224追赶法：追赶法：iniiniiniiCCC11111方程的统一形式：方程的统一形式：253/7/202225追赶法：追赶法：102111121111nnnCxDCC引入边界条件：引入边界条件：当当i=1时：时：当当i=M时：时：MMMMMMMnMMnMMnMnMnMCCCCC2211111111其中则令263/7/202226

16、追赶法：追赶法：令令AC=MMnMMMMMMMCCCCCA121112111133322211273/7/202227追赶法：追赶法：追赶法解方程组的方法：追赶法解方程组的方法：先把对角矩阵先把对角矩阵A分解为：分解为：A=BW其中：其中：MMMMbabababB112211111121MwwwW283/7/202228追赶法：追赶法：由于由于A=BW1112111332333221222111MMMMMMMMMMwabawbwabawbwabawbwabawbbBWA293/7/202229追赶法：追赶法：由此可得：由此可得： b1=1 b1w1=1 w1=1/b1 b2+a2w1=2 bi

17、wi=i wi=i/bi (i=1M-1) bi+aiwi-1=i bi=i-aiwi-1 bM=M-aMwM-1303/7/202230追赶法：追赶法：由于由于A=BW B(WC)=又令：又令：WC=Y 其中：其中：Y=(y1, y2, , yM)T BY=MMMMMMMMyyyybabababB12112111221313/7/202231追赶法：追赶法：由以上等式可知：由以上等式可知： b1y1=1 y1=1 /b1 a2y1+b2y2=2 y2=(2-a2y1)/ b2 aiyi-1+biyi=i yi=(i-aiyi-1)/ bi aMyM-1+bMyM=M yM=(M-aMyM-1

18、)/ bM至此可以求出所有的至此可以求出所有的Y(y1, y2, , yM)T值。值。323/7/202232追赶法：追赶法：又因为又因为WC=Y，即：，即：MMnMMMyyyyCCCCwwwCW12111211211111333/7/202233追赶法：追赶法：根据上述等式，可得：根据上述等式，可得：MnMMnMMnMniiininnnnyCyCwCCwyCyCwCyCwC111111111213212112111343/7/202234追赶法：追赶法：追：追：根据根据BY=，求出，求出bi和和wi，并求出，并求出yi；赶：赶：根据根据WC=Y，求出，求出Cin+1。353/7/202235

19、二维弥散方程的差分格式二维弥散方程的差分格式若多孔介质是各向同性的，流场为均匀的若多孔介质是各向同性的，流场为均匀的一维流场，渗透速度一维流场，渗透速度V=nu，流动方向取，流动方向取x轴轴的正向，与流速垂直为的正向，与流速垂直为y方向，对流扩散方方向，对流扩散方程可简化为：程可简化为： xcuycDxcDtcTL2222363/7/202236二维弥散方程的差分格式二维弥散方程的差分格式将研究区域网格化将研究区域网格化tntyjcyxiaxNttMcdyMabxnji,/,/ )(,/ )(373/7/202237二维弥散方程的差分格式二维弥散方程的差分格式建立矩形网差分格式建立矩形网差分格

20、式l隐式差分格式隐式差分格式tCCtcxCCxcnjinjinjinji,1,1,11,12211,1,11,2221,11,1,12222yCCCycxCCCxcnjinjinjinjinjinji383/7/202238二维弥散方程的差分格式二维弥散方程的差分格式于是，内节点的差分方程为于是，内节点的差分方程为xCCVyCCCDxCCCDtCCnjinjinjinjinjiTnjinjinjiLnjinji2221,11,1211,1,11,21,11,1,1,1,393/7/202239二维弥散方程的差分格式二维弥散方程的差分格式整理后可得整理后可得njinjiTnjiTnjiLnjiT

21、LnjiLCCyTDCyTDCxtVxtDCytDxtDCxtVxtD,11,211,21,121,221,12)2()221()2(403/7/202240二维弥散方程的差分格式二维弥散方程的差分格式边界条件边界条件初始条件初始条件),(),(),(),(1,10,1,1,0ninMinininjnjMnjnjtdxCCtcxCCtybCCtyaCC),(00,jijiyxCC413/7/202241二维弥散方程的差分格式二维弥散方程的差分格式将网格取为正方形，即将网格取为正方形，即x=y=h，方程就，方程就可简化为：可简化为：njinjiTnjiTnjiLnjiTLnjiLCTCDhTCD

22、hCtVhtDhCtDtDhhCtVhtDh,11,211,21,121,221,1211)(21)22(1)(21423/7/202242二维弥散方程的差分格式二维弥散方程的差分格式交替方向差分格式（交替方向差分格式（ADI） ADI法的实质是法的实质是将二维隐式差分格式化为一将二维隐式差分格式化为一维的隐式差分格式维的隐式差分格式。 在第在第n层与第层与第n+1层两层正中间引入一个过层两层正中间引入一个过渡层渡层n+1/2。433/7/202243二维弥散方程的差分格式二维弥散方程的差分格式交替方向差分格式（交替方向差分格式（ADI） 分两步走：分两步走：在前半时间步长内，对在前半时间步长

23、内，对x偏导用隐式偏导用隐式表示，对表示，对y的偏导用显示表示；在后半时间步长的偏导用显示表示；在后半时间步长内，对内，对x偏导用显示表示，对偏导用显示表示，对y的偏导用隐式表的偏导用隐式表示。示。443/7/202244二维弥散方程的差分格式二维弥散方程的差分格式交替方向差分格式（交替方向差分格式（ADI） 前半时间步长：前半时间步长： 21,1,1,12,1,1,2222/212121212121yCCCDxCCVxCCCDtCCnjinjinjiTnjinjinjinjinjiLnjinji453/7/202245二维弥散方程的差分格式二维弥散方程的差分格式交替方向差分格式（交替方向差分

24、格式（ADI） 后半时间步长：后半时间步长： 211,1,11,1,12,1,1,1,2222/212121212121yCCCDxCCVxCCCDtCCnjinjinjiTnjinjinjinjinjiLnjinji463/7/202246二维弥散方程的差分格式二维弥散方程的差分格式计算步骤：计算步骤：（1）根据初始条件或前一时段计算结果，可知）根据初始条件或前一时段计算结果，可知tn时刻的浓度值，那么第一个方程中只有三个未时刻的浓度值，那么第一个方程中只有三个未知数（知数（n+1/2时刻）。如果固定时刻）。如果固定j值，对值，对i的联立的联立方程组就是三对角线的（追赶法），即方程组就是三对

25、角线的（追赶法），即沿沿x方向方向追赶追赶，可求出过渡层（，可求出过渡层（n+1/2）层上的全部浓度）层上的全部浓度值。值。 473/7/202247二维弥散方程的差分格式二维弥散方程的差分格式计算步骤：计算步骤：（2）根据前一步求出的过渡层（）根据前一步求出的过渡层（n+1/2）层上的）层上的浓度值，再按同样的方法求出第二个方程中的浓度值，再按同样的方法求出第二个方程中的（n+1）时刻的浓度值。即对于一个固定的）时刻的浓度值。即对于一个固定的i值，值，方程也是一个三对角方程组（追赶法），方程也是一个三对角方程组（追赶法），沿沿y方方向追赶向追赶。 483/7/202248有限单元法基本原理有

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

27、的边界剖分为若干单元，将求函数通过把求解区域的边界剖分为若干单元，将求函数的解简化为求节点的函数值，求解积分方程就化为的解简化为求节点的函数值，求解积分方程就化为求解一组线性代数方程，在求出边界上的物理量以求解一组线性代数方程，在求出边界上的物理量以后，计算区域内部的任一点的未知量可通过边界上后，计算区域内部的任一点的未知量可通过边界上的已知量求出。的已知量求出。与有限元法相比，具有降低维数、输入数据准备简与有限元法相比，具有降低维数、输入数据准备简单、计算工作量小、精度比有限元法高等优点。单、计算工作量小、精度比有限元法高等优点。在处理非均质、非线性、非稳定等复杂问题时计算在处理非均质、非线

28、性、非稳定等复杂问题时计算效率低，不能将对流项化成沿边界的积分。效率低，不能将对流项化成沿边界的积分。 503/7/202250第四节第四节 数理统计法数理统计法地下水水质受多种因素综合影响，具有不地下水水质受多种因素综合影响，具有不确定性，可用数理统计法预测水质变化；确定性，可用数理统计法预测水质变化；主要包括：时间序列分析法、相关分析法、主要包括：时间序列分析法、相关分析法、概率统计法。概率统计法。513/7/202251（1）时间序列分析法）时间序列分析法基本思路基本思路：分析时间序列的变化特征，选：分析时间序列的变化特征，选择适当的模式建立预测模型，利用模型进择适当的模式建立预测模型，

29、利用模型进行趋势外推预测，对模型预测值进行评价行趋势外推预测，对模型预测值进行评价和修正，得到预测结果。和修正，得到预测结果。523/7/202252（1）时间序列分析法）时间序列分析法-平均数平均数平均数预测平均数预测：对动态数列中若干个时期的数值，：对动态数列中若干个时期的数值，求其算术平均值或加权平均值，作为下一个时求其算术平均值或加权平均值，作为下一个时期的预测值。期的预测值。niinCnC111niiinCWnC11)(1533/7/202253（1）时间序列分析法）时间序列分析法-指数平滑指数平滑一次指数平滑预测一次指数平滑预测：对动态数列中若干个时期：对动态数列中若干个时期的数值

30、，求其算术平均值或加权平均值，作为的数值，求其算术平均值或加权平均值，作为下一个时期的预测值。下一个时期的预测值。 221)1()1 ()1 (ttttCaaCaaaCS)1(121)1()1 ()1 ()1 ( ttttttSaaCCaaaCaaCS543/7/202254（1）时间序列分析法）时间序列分析法-指数平滑指数平滑一次指数平滑法是以最近周期的一次指数平滑一次指数平滑法是以最近周期的一次指数平滑值作为下一周期的预测值，即：值作为下一周期的预测值，即：其中平滑系数其中平滑系数a的选择直接影响预测效果。的选择直接影响预测效果。如如果时间序列长期趋势稳定，应取较小的果时间序列长期趋势稳定

31、，应取较小的a值；值；若时间序列不稳定，应取较大的若时间序列不稳定，应取较大的a值。值。)1(1)1(1)1 (ttttSaaCSC553/7/202255（1）时间序列分析法）时间序列分析法-指数平滑指数平滑a越小，对数据的平滑能力越强，但对数据变化的敏感越小，对数据的平滑能力越强，但对数据变化的敏感性越差；性越差；a越大，对数据的平滑能力越差，但对数据变越大，对数据的平滑能力越差，但对数据变化的敏感性越强。化的敏感性越强。563/7/202256（1）时间序列分析法）时间序列分析法-指数平滑指数平滑初始值初始值S0(1)对预测结果的影响很大，一般可取第一个数对预测结果的影响很大，一般可取第

32、一个数据（据（x1）作为初始值；若总体数据量较少，则可取最）作为初始值；若总体数据量较少，则可取最初几个数据的平均值作为初始值。初几个数据的平均值作为初始值。573/7/202257（1）时间序列分析法）时间序列分析法-指数平滑指数平滑二次指数平滑预测二次指数平滑预测：对一次指数平滑值再做一：对一次指数平滑值再做一次指数平滑：次指数平滑：)2(1)1()2()1 (tttSaaSS583/7/202258（1）时间序列分析法）时间序列分析法-指数平滑指数平滑二次指数平滑值不直接用于预测，而是根据滞二次指数平滑值不直接用于预测，而是根据滞后偏差建立线性预测模型：后偏差建立线性预测模型：其中：其中

33、：TbaCttTt)2()1(2tttSSa1)2()1(tttSSaab593/7/202259（1）时间序列分析法）时间序列分析法-指数平滑指数平滑若一、二次指数平滑值不满足预测要求，应采若一、二次指数平滑值不满足预测要求，应采用三次指数平滑法建立非线性预测模型：用三次指数平滑法建立非线性预测模型：(3)1)2()3()1 (tttSaaSS603/7/202260（1）时间序列分析法）时间序列分析法-指数平滑指数平滑三次指数平滑法建立的非线性预测模型：三次指数平滑法建立的非线性预测模型：)3()2()1(333ttttSSSa)34()810(2)56()12()3()2()1(2ttt

34、tSaSaSaaab2TtttTtcTbaC2)12()3()2()1(2ttttSSSaac613/7/202261（2）相关分析法相关分析法线性相关预测线性相关预测：建立直角坐标下数据散点图，：建立直角坐标下数据散点图，用一元线性回归方程近似描述数据点的分布规用一元线性回归方程近似描述数据点的分布规律：律：bxay623/7/202262（2）相关分析法相关分析法根据最小二乘法根据最小二乘法(要求残差要求残差RSS最小最小)可求得：可求得：用相关系数用相关系数r来判断回归方程实用性：来判断回归方程实用性：22222)()(iiiiiiiiiiiiixnxyxnyxbxnxxyxyxa22)

35、()()(yyxxyyxxriiii2)(iibxayRSS633/7/202263（2）相关分析法相关分析法曲性相关预测曲性相关预测：将直角坐标下散点图拟合形成：将直角坐标下散点图拟合形成的有规律的曲线进行函数转换，变为前述的一的有规律的曲线进行函数转换，变为前述的一元线性回归方程进行计算。元线性回归方程进行计算。如：如：baxy bXaYyYxXlglg,lg则令643/7/202264国外地下水模拟软件的发展现状与趋势国外地下水模拟软件的发展现状与趋势 利用数值模型对地下水流和溶质运移问题进行利用数值模型对地下水流和溶质运移问题进行模拟的方法以其有效性、灵活性和相对廉价性模拟的方法以其有

36、效性、灵活性和相对廉价性逐渐成为地下水研究领域的一种不可或缺的重逐渐成为地下水研究领域的一种不可或缺的重要方法，并受到越来越大的重视和广泛的应用。要方法，并受到越来越大的重视和广泛的应用。一个完整的地下水模拟过程包含一个完整的地下水模拟过程包含3个部分：个部分：前处前处理、模型计算和后处理理、模型计算和后处理。653/7/202265国外地下水模拟软件的发展现状与趋势国外地下水模拟软件的发展现状与趋势前处理前处理是指在进行模拟计算之前对计算过程中是指在进行模拟计算之前对计算过程中所需数据的整理、组织、输入及计算网格的编所需数据的整理、组织、输入及计算网格的编号与生成。号与生成。模型计算模型计算

37、是进行地下水流动或水质运移正反演是进行地下水流动或水质运移正反演计算，常用的方法主要有：有限差分法、有限计算，常用的方法主要有：有限差分法、有限元法、边界元法等。元法、边界元法等。后处理后处理是将计算所产生的结果数据，用图形或是将计算所产生的结果数据，用图形或表格显示或存放起来，以供研究人员方便地进表格显示或存放起来，以供研究人员方便地进行分析和使用。行分析和使用。663/7/202266国外地下水模拟软件的发展现状与趋势国外地下水模拟软件的发展现状与趋势通过近二十年的研究与发展，国际上已经形成通过近二十年的研究与发展，国际上已经形成了一批非常有影响的地下水模拟的软件，它们了一批非常有影响的地

38、下水模拟的软件，它们今天在国际地下水模拟研究领域依旧非常活跃，今天在国际地下水模拟研究领域依旧非常活跃，如如MODFLOW、FEFLOW、MT3DMS、MT3D99、PEST、MODPATH、UCODE等。等。 673/7/202267MODFLOW MODFLOW是由美国地质调查局的是由美国地质调查局的McDonald和和Harbaugh于于80年代开发出来的一套专门用年代开发出来的一套专门用于孔隙介质中于孔隙介质中三维有限差分三维有限差分地下水流数值模拟地下水流数值模拟的软件。的软件。l程序结构的程序结构的模块化模块化（模拟抽水引起地面沉降的子（模拟抽水引起地面沉降的子程序包、模拟水平流动

39、障碍的子程序包等）。程序包、模拟水平流动障碍的子程序包等）。l离散方法的离散方法的简单化简单化（有限差分法（有限差分法36036018个个网格单元）。网格单元）。lMODFLOW引进了引进了应力期应力期(Stress Period)概念。概念。 l求解方法的多样化（强隐式法、逐次超松弛迭代求解方法的多样化（强隐式法、逐次超松弛迭代法、预调共轭梯度法）。法、预调共轭梯度法）。 683/7/202268MT3D99 MT3D99是美国是美国S.S.Papodopulos & Associates公司的公司的Zheng博士设计的模拟博士设计的模拟三维地下水溶质三维地下水溶质运移运移程序程序M

40、T3D(1990)的增强版，成为目前世的增强版，成为目前世界上首屈一指的溶质运移模拟软件。界上首屈一指的溶质运移模拟软件。lMT3D99能够模拟地下水系统中的平流、扩散、衰能够模拟地下水系统中的平流、扩散、衰减、溶质化学反应、线性与非线性吸附作用等现象，减、溶质化学反应、线性与非线性吸附作用等现象，能够对承压含水层，不承压含水层，承压与不承压能够对承压含水层，不承压含水层，承压与不承压交替的含水层以及倾斜的和单元厚度变化的含水层交替的含水层以及倾斜的和单元厚度变化的含水层进行空间离散。进行空间离散。 lMT3D99提供了丰富的求解方法（标准有限差分法、提供了丰富的求解方法（标准有限差分法、基于

41、基于Eulerian-Lagrangian的粒子跟踪方法和高阶有的粒子跟踪方法和高阶有限体积限体积TVD方法）。方法）。 693/7/202269PEST2000 PEST2000是由澳大利亚是由澳大利亚Watermark Computing公公司开发的司开发的PEST新版本，是功能强大的模型独立的新版本，是功能强大的模型独立的参数估计参数估计程序程序。lPEST2000利用一个强有力的数值反演算法来利用一个强有力的数值反演算法来“控制控制”运行中的模型，程序在每次模拟之后自运行中的模型，程序在每次模拟之后自动调整所选择的模型参数，直到将校正的目标动调整所选择的模型参数，直到将校正的目标最小化为止。最小化为止。lPEST2000引进了参数估计的最新技术引进了参数估计的最新技术预测预测分析。分析。 703/7/202270FEFLOW FEFLOW是加拿大是加拿大Wate