第五章 地下水污染預(yù)測

第五章地下水污染預(yù)測

主要內(nèi)容概述（了解）

近似解析法（了解）

數(shù)值法（熟悉）

數(shù)理統(tǒng)計(jì)法（熟悉）灰色預(yù)測法（了解）

國外地下水模型軟件2/5/20231第一節(jié)概述

地下水污染預(yù)測預(yù)報(bào)包括：預(yù)報(bào)污染地下水分布邊界的推進(jìn)情況（即經(jīng)過一定時(shí)間t后邊界推進(jìn)的距離x，或預(yù)測推進(jìn)一定距離需要多長時(shí)間）；確定污染地下水中污染物質(zhì)的濃度在空間的分布及隨時(shí)間的變化情況；預(yù)測污染地下水能否侵入附近的水源地，污染物質(zhì)達(dá)到水源地所需的時(shí)間；預(yù)測防治地下水污染的措施的效果等。2/5/20232第一節(jié)概述

預(yù)測評價(jià)的基本步驟如下：建立預(yù)測模型；模型識別（平衡性、穩(wěn)定性、敏感性）；模型參數(shù)估計(jì)；模型方程的驗(yàn)證；地下水污染預(yù)測。2/5/20233第二節(jié)近似解析解前面介紹了一些解析解公式，但應(yīng)用起來仍十分復(fù)雜，不便于實(shí)際應(yīng)用。實(shí)際中可以進(jìn)一步簡化：完全忽略彌散作用，將污染地下水的運(yùn)動(dòng)視為“活塞式”的推擠淡水的運(yùn)動(dòng)，假設(shè)兩者始終保持明顯的鉛直分界面；將地下水動(dòng)力運(yùn)移和彌散分開考慮，先以上述方法計(jì)算，然后計(jì)算由于彌散所形成的過渡混合帶，修正按平均鋒面運(yùn)移的距離。2/5/20234第二節(jié)近似解析解

計(jì)算平均滲透鋒面的運(yùn)移距離；求出平均滲透鋒面處污染物的濃度；確定分界面彌散帶和變形帶，在平均滲透鋒面距離上加上這兩個(gè)帶就可確定污染水實(shí)際鋒面的運(yùn)移位置：

相對密度差巖層傾角2/5/20235第二節(jié)近似解析解承壓（或潛水）含水層污染地下水在承壓含水層的運(yùn)動(dòng)，如果補(bǔ)給來源穩(wěn)定，可視為似穩(wěn)定運(yùn)動(dòng)，屬于一維平面平行流，則滲透速度為：而實(shí)際流速為：地下水的天然單寬流量2/5/20236第二節(jié)近似解析解承壓（或潛水）含水層初始時(shí)刻污染鋒面的位置某時(shí)刻污染鋒面的位置天然地下水中污染物的濃度初始分界面以內(nèi)，地下水中污染物的濃度若是潛水，則取平均厚度代替上述各式中的m。2/5/20237第三節(jié)數(shù)值解由于含水層以及污染物質(zhì)遷移的復(fù)雜性，解析解不易求得，通常采用數(shù)值解法。數(shù)值法在剖分后認(rèn)為每一單元內(nèi)均質(zhì)同性，將方程線性化。有限差分法有限單元法邊界元法2/5/20238一、有限差分法有限差分法的基本思路是按照時(shí)間步長和空間步長將時(shí)間和空間區(qū)域剖分成若干網(wǎng)格，用未知函數(shù)在網(wǎng)格節(jié)點(diǎn)上的值所構(gòu)成的差分商近似代替所用偏微分方程中出現(xiàn)的各階導(dǎo)數(shù)，從而把表示變量連續(xù)變化關(guān)系的偏微分方程離散為有限個(gè)代數(shù)方程，然后解此線性代數(shù)方程，以求出污染物濃度在各網(wǎng)格節(jié)點(diǎn)上不同時(shí)刻的解。

2/5/20239導(dǎo)數(shù)的差分近似設(shè)f(x)為任一足夠光滑的函數(shù)，把f(x)沿x的正向展開為泰勒級數(shù)：上式稱為f(x)的一階向前差分2/5/202310導(dǎo)數(shù)的差分近似如果把f(x)對x的負(fù)向展開為泰勒級數(shù)：上式稱為f(x)的一階向后差分2/5/202311導(dǎo)數(shù)的差分近似把上面兩個(gè)差分方程相減，可得：上式稱為f(x)的一階中心差分同理可得f(x)的二階導(dǎo)數(shù)的差分公式：2/5/202312導(dǎo)數(shù)的差分近似若把上述一元函數(shù)f(x)換為濃度函數(shù)，并按網(wǎng)格節(jié)點(diǎn)編號，則對x的一階向前、向后和中心差分分別為：同樣可以寫出對時(shí)間t的向前、向后和中心差分公式2/5/202313導(dǎo)數(shù)的差分近似C(x,y,z,t)對x的二階差分公式為：因此，在任何一個(gè)節(jié)點(diǎn)處，濃度C(x,y,z,t)的一階和二階偏導(dǎo)數(shù)均可用該節(jié)點(diǎn)及其相鄰節(jié)點(diǎn)上濃度值的線性組合表示。2/5/202314一維彌散方程的差分格式設(shè)在無限含水層中，存在一維均勻流場，其滲透速度V=nu，流動(dòng)方向?yàn)閤軸正向。則一維彌散方程為：2/5/202315一維彌散方程的差分格式對時(shí)間區(qū)域[0,T]和空間區(qū)域[0,L]都作等距剖分，設(shè)時(shí)間步長為△t，空間步長為△x，把第i個(gè)節(jié)點(diǎn)xi處在tn時(shí)刻的濃度記為Cin。2/5/202316一維彌散方程的差分格式顯示差分格式將及差分公式中的濃度取為tn時(shí)刻的值，便可得到其顯示差分格式：2/5/202317一維彌散方程的差分格式顯示差分格式的穩(wěn)定條件在應(yīng)用顯示差分格式時(shí)，剖分的步長必須滿足：2/5/202318一維彌散方程的差分格式隱示差分格式將及差分公式中的濃度取為tn＋1時(shí)刻的值，便可得到其隱示差分格式：2/5/202319一維彌散方程的差分格式隱示差分格式的穩(wěn)定條件在應(yīng)用隱示差分格式時(shí)，剖分的步長必須滿足：2/5/202320一維彌散方程的差分格式隱示差分格式的矩陣形式2/5/202321一維彌散方程的差分格式Crank-Nicolson差分格式取顯示差分格式和隱示差分格式的平均，便可得到C-N差分格式：2/5/202322一維彌散方程的差分格式Crank-Nicolson差分格式整理后可得：2/5/202323追趕法：方程的統(tǒng)一形式：2/5/202324追趕法：引入邊界條件：當(dāng)i=1時(shí)：當(dāng)i=M時(shí)：2/5/202325追趕法：令A(yù)·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追趕法：又因?yàn)閃·C=Y，即：2/5/202332追趕法：根據(jù)上述等式，可得：2/5/202333追趕法：追：根據(jù)B·Y=δ，求出bi和wi，并求出yi；趕：根據(jù)W·C=Y，求出Cin+1。2/5/202334二維彌散方程的差分格式若多孔介質(zhì)是各向同性的，流場為均勻的一維流場，滲透速度V=nu，流動(dòng)方向取x軸的正向，與流速垂直為y方向，對流擴(kuò)散方程可簡化為：2/5/202335二維彌散方程的差分格式將研究區(qū)域網(wǎng)格化2/5/202336二維彌散方程的差分格式建立矩形網(wǎng)差分格式隱式差分格式2/5/202337二維彌散方程的差分格式于是，內(nèi)節(jié)點(diǎn)的差分方程為2/5/202338二維彌散方程的差分格式整理后可得2/5/202339二維彌散方程的差分格式邊界條件初始條件2/5/202340二維彌散方程的差分格式將網(wǎng)格取為正方形，即△x=△y=h，方程就可簡化為：2/5/202341二維彌散方程的差分格式交替方向差分格式（ADI）ADI法的實(shí)質(zhì)是將二維隱式差分格式化為一維的隱式差分格式。在第n層與第n+1層兩層正中間引入一個(gè)過渡層n+1/2。2/5/202342二維彌散方程的差分格式交替方向差分格式（ADI）

分兩步走：在前半時(shí)間步長內(nèi)，對x偏導(dǎo)用隱式表示，對y的偏導(dǎo)用顯示表示；在后半時(shí)間步長內(nèi)，對x偏導(dǎo)用顯示表示，對y的偏導(dǎo)用隱式表示。2/5/202343二維彌散方程的差分格式交替方向差分格式（ADI）前半時(shí)間步長：

2/5/202344二維彌散方程的差分格式交替方向差分格式（ADI）后半時(shí)間步長：

2/5/202345二維彌散方程的差分格式計(jì)算步驟：（1）根據(jù)初始條件或前一時(shí)段計(jì)算結(jié)果，可知tn時(shí)刻的濃度值，那么第一個(gè)方程中只有三個(gè)未知數(shù)（n+1/2時(shí)刻）。如果固定j值，對i的聯(lián)立方程組就是三對角線的（追趕法），即沿x方向追趕，可求出過渡層（n+1/2）層上的全部濃度值。

2/5/202346二維彌散方程的差分格式計(jì)算步驟：（2）根據(jù)前一步求出的過渡層（n+1/2）層上的濃度值，再按同樣的方法求出第二個(gè)方程中的（n+1）時(shí)刻的濃度值。即對于一個(gè)固定的i值，方程也是一個(gè)三對角方程組（追趕法），沿y方向追趕。

2/5/202347有限單元法基本原理在有限差分法中，研究區(qū)域被一系列的網(wǎng)格節(jié)點(diǎn)所代替，在各個(gè)節(jié)點(diǎn)上以差分代替導(dǎo)數(shù)而得到求解相應(yīng)偏微分方程的線性代數(shù)方程組。有限單元法是把研究區(qū)域剖分為有限個(gè)子區(qū)域，在每個(gè)子區(qū)域上用某種插值函數(shù)來近似待求解的未知函數(shù)，從而得到求解相應(yīng)偏微分方程的線性代數(shù)方程組。

2/5/202348邊界元法基本原理通過把求解區(qū)域的邊界剖分為若干單元，將求函數(shù)的解簡化為求節(jié)點(diǎn)的函數(shù)值，求解積分方程就化為求解一組線性代數(shù)方程，在求出邊界上的物理量以后，計(jì)算區(qū)域內(nèi)部的任一點(diǎn)的未知量可通過邊界上的已知量求出。與有限元法相比，具有降低維數(shù)、輸入數(shù)據(jù)準(zhǔn)備簡單、計(jì)算工作量小、精度比有限元法高等優(yōu)點(diǎn)。在處理非均質(zhì)、非線性、非穩(wěn)定等復(fù)雜問題時(shí)計(jì)算效率低，不能將對流項(xiàng)化成沿邊界的積分。2/5/202349第四節(jié)數(shù)理統(tǒng)計(jì)法地下水水質(zhì)受多種因素綜合影響，具有不確定性，可用數(shù)理統(tǒng)計(jì)法預(yù)測水質(zhì)變化；主要包括：時(shí)間序列分析法、相關(guān)分析法、概率統(tǒng)計(jì)法。2/5/202350（1）時(shí)間序列分析法基本思路：分析時(shí)間序列的變化特征，選擇適當(dāng)?shù)哪Ｊ浇㈩A(yù)測模型，利用模型進(jìn)行趨勢外推預(yù)測，對模型預(yù)測值進(jìn)行評價(jià)和修正，得到預(yù)測結(jié)果。2/5/202351（1）時(shí)間序列分析法-平均數(shù)平均數(shù)預(yù)測：對動(dòng)態(tài)數(shù)列中若干個(gè)時(shí)期的數(shù)值，求其算術(shù)平均值或加權(quán)平均值，作為下一個(gè)時(shí)期的預(yù)測值。2/5/202352（1）時(shí)間序列分析法-指數(shù)平滑一次指數(shù)平滑預(yù)測：對動(dòng)態(tài)數(shù)列中若干個(gè)時(shí)期的數(shù)值，求其算術(shù)平均值或加權(quán)平均值，作為下一個(gè)時(shí)期的預(yù)測值。2/5/202353（1）時(shí)間序列分析法-指數(shù)平滑一次指數(shù)平滑法是以最近周期的一次指數(shù)平滑值作為下一周期的預(yù)測值，即：其中平滑系數(shù)a的選擇直接影響預(yù)測效果。如果時(shí)間序列長期趨勢穩(wěn)定，應(yīng)取較小的a值；若時(shí)間序列不穩(wěn)定，應(yīng)取較大的a值。2/5/202354（1）時(shí)間序列分析法-指數(shù)平滑a越小，對數(shù)據(jù)的平滑能力越強(qiáng)，但對數(shù)據(jù)變化的敏感性越差；a越大，對數(shù)據(jù)的平滑能力越差，但對數(shù)據(jù)變化的敏感性越強(qiáng)。2/5/202355（1）時(shí)間序列分析法-指數(shù)平滑初始值S0(1)對預(yù)測結(jié)果的影響很大，一般可取第一個(gè)數(shù)據(jù)（x1）作為初始值；若總體數(shù)據(jù)量較少，則可取最初幾個(gè)數(shù)據(jù)的平均值作為初始值。2/5/202356（1）時(shí)間序列分析法-指數(shù)平滑二次指數(shù)平滑預(yù)測：對一次指數(shù)平滑值再做一次指數(shù)平滑：2/5/202357（1）時(shí)間序列分析法-指數(shù)平滑二次指數(shù)平滑值不直接用于預(yù)測，而是根據(jù)滯后偏差建立線性預(yù)測模型：其中：2/5/202358（1）時(shí)間序列分析法-指數(shù)平滑若一、二次指數(shù)平滑值不滿足預(yù)測要求，應(yīng)采用三次指數(shù)平滑法建立非線性預(yù)測模型：2/5/202359（1）時(shí)間序列分析法-指數(shù)平滑三次指數(shù)平滑法建立的非線性預(yù)測模型：2/5/202360（2）相關(guān)分析法線性相關(guān)預(yù)測：建立直角坐標(biāo)下數(shù)據(jù)散點(diǎn)圖，用一元線性回歸方程近似描述數(shù)據(jù)點(diǎn)的分布規(guī)律：2/5/202361（2）相關(guān)分析法根據(jù)最小二乘法(要求殘差RSS最小)可求得：用相關(guān)系數(shù)r來判斷回歸方程實(shí)用性：2/5/202362（2）相關(guān)分析法曲性相關(guān)預(yù)測：將直角坐標(biāo)下散點(diǎn)圖擬合形成的有規(guī)律的曲線進(jìn)行函數(shù)轉(zhuǎn)換，變?yōu)榍笆龅囊辉€性回歸方程進(jìn)行計(jì)算。如：2/5/202363國外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢

利用數(shù)值模型對地下水流和溶質(zhì)運(yùn)移問題進(jìn)行模擬的方法以其有效性、靈活性和相對廉價(jià)性逐漸成為地下水研究領(lǐng)域的一種不可或缺的重要方法，并受到越來越大的重視和廣泛的應(yīng)用。一個(gè)完整的地下水模擬過程包含3個(gè)部分：前處理、模型計(jì)算和后處理。2/5/202364國外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢前處理是指在進(jìn)行模擬計(jì)算之前對計(jì)算過程中所需數(shù)據(jù)的整理、組織、輸入及計(jì)算網(wǎng)格的編號與生成。模型計(jì)算是進(jìn)行地下水流動(dòng)或水質(zhì)運(yùn)移正反演計(jì)算，常用的方法主要有：有限差分法、有限元法、邊界元法等。后處理是將計(jì)算所產(chǎn)生的結(jié)果數(shù)據(jù)，用圖形或表格顯示或存放起來，以供研究人員方便地進(jìn)行分析和使用。2/5/202365國外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢通過近二十年的研究與發(fā)展，國際上已經(jīng)形成了一批非常有影響的地下水模擬的軟件，它們今天在國際地下水模擬研究領(lǐng)域依舊非?；钴S，如MODFLOW、FEFLOW、MT3DMS、MT3D99、PEST、MODPATH、UCODE等。2/5/202366MODFLOW

MODFLOW是由美國地質(zhì)調(diào)查局的McDonald和Harbaugh于80年代開發(fā)出來的一套專門用于孔隙介質(zhì)中三維有限差分地下水流數(shù)值模擬的軟件。程序結(jié)構(gòu)的模塊化（模擬抽水引起地面沉降的子程序包、模擬水平流動(dòng)障礙的子程序包等）。離散方法的簡單化（有限差分法360×360×18個(gè)網(wǎng)格單元）。MODFLOW引進(jìn)了應(yīng)力期(StressPeriod)概念。求解方法的多樣化（強(qiáng)隱式法、逐次超松弛迭代法、預(yù)調(diào)共軛梯度法）。2/5/202367MT3D99

MT3D99是美國S.S.Papodopulos&Associates公司的Zheng博士設(shè)計(jì)的模擬三維地下水溶質(zhì)運(yùn)移程序MT3D(1990)的增強(qiáng)版，成為目前世界上首屈一指的溶質(zhì)運(yùn)移模擬軟件。MT3D99能夠模擬地下水系統(tǒng)中的平流、擴(kuò)散、衰減、溶質(zhì)化學(xué)反應(yīng)、線性與非線性吸附作用等現(xiàn)象，能夠?qū)Τ袎汉畬樱怀袎汉畬?，承壓與不承壓交替的含水層以及傾斜的和單元厚度變化的含水層進(jìn)行空間離散。

MT3D99提供了豐富的求解方法（標(biāo)準(zhǔn)有限差分法、基于Eulerian-Lagrangian的粒子跟蹤方法和高階有限體積TVD方法）。

2/5/202368PEST2000PEST2000是由澳大利亞WatermarkComputing公司開發(fā)的PEST新版本，是功能強(qiáng)大的模型獨(dú)立的參數(shù)估計(jì)程序。PEST2000利用一個(gè)強(qiáng)有力的數(shù)值反演算法來“控制”運(yùn)行中的模型，程序在每次模擬之后自動(dòng)調(diào)整所選擇的模型參數(shù)，直到將校正的目標(biāo)最小化為止。PEST2000引進(jìn)了參數(shù)估計(jì)的最新技術(shù)──預(yù)測分析。

2/5/202369FEFLOW

FEFLOW是加拿大Waterloo水文地質(zhì)公司開發(fā)的基于三維有限元的地下水模擬可視化軟件包。它能夠解決下列地下水模擬問題：完全瞬時(shí)、半瞬時(shí)、穩(wěn)態(tài)地下水流動(dòng)與溶質(zhì)運(yùn)移；隨時(shí)間變化的實(shí)體屬性和約束邊界條件；飽和與不飽和流動(dòng)；包含潛水面的承壓與不承壓含水層；帶有非線性吸附作用、衰變、對流、彌散的化學(xué)質(zhì)量運(yùn)移；考慮貯存、對流、熱散失、熱運(yùn)移的流體和固體熱量運(yùn)移；密度變化的流動(dòng)(海水入侵等)。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