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

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

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

近似解析法（了解）

數(shù)值法（熟悉）

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

國外地下水模型軟件12/9/20241第一節(jié)概述

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

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

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

相對密度差巖層傾角12/9/20245第二節(jié)近似解析解承壓（或潛水）含水層污染地下水在承壓含水層的運(yùn)動，如果補(bǔ)給來源穩(wěn)定，可視為似穩(wěn)定運(yùn)動，屬于一維平面平行流，則滲透速度為：而實(shí)際流速為：地下水的天然單寬流量12/9/20246第二節(jié)近似解析解承壓（或潛水）含水層初始時刻污染鋒面的位置某時刻污染鋒面的位置天然地下水中污染物的濃度初始分界面以內(nèi)，地下水中污染物的濃度若是潛水，則取平均厚度代替上述各式中的m。12/9/20247第三節(jié)數(shù)值解由于含水層以及污染物質(zhì)遷移的復(fù)雜性，解析解不易求得，通常采用數(shù)值解法。數(shù)值法在剖分后認(rèn)為每一單元內(nèi)均質(zhì)同性，將方程線性化。有限差分法有限單元法邊界元法12/9/20248一、

有限差分法有限差分法的基本思路是按照時間步長和空間步長將時間和空間區(qū)域剖分成若干網(wǎng)格，用未知函數(shù)在網(wǎng)格節(jié)點(diǎn)上的值所構(gòu)成的差分商近似代替所用偏微分方程中出現(xiàn)的各階導(dǎo)數(shù)，從而把表示變量連續(xù)變化關(guān)系的偏微分方程離散為有限個代數(shù)方程，然后解此線性代數(shù)方程，以求出污染物濃度在各網(wǎng)格節(jié)點(diǎn)上不同時刻的解。

12/9/20249導(dǎo)數(shù)的差分近似設(shè)f(x)為任一足夠光滑的函數(shù)，把f(x)沿x的正向展開為泰勒級數(shù)：上式稱為f(x)的一階向前差分12/9/202410導(dǎo)數(shù)的差分近似如果把f(x)對x的負(fù)向展開為泰勒級數(shù)：上式稱為f(x)的一階向后差分12/9/202411導(dǎo)數(shù)的差分近似把上面兩個差分方程相減，可得：上式稱為f(x)的一階中心差分同理可得f(x)的二階導(dǎo)數(shù)的差分公式：12/9/202412導(dǎo)數(shù)的差分近似若把上述一元函數(shù)f(x)換為濃度函數(shù)，并按網(wǎng)格節(jié)點(diǎn)編號，則對x的一階向前、向后和中心差分分別為：同樣可以寫出對時間t的向前、向后和中心差分公式12/9/202413導(dǎo)數(shù)的差分近似C(x,y,z,t)對x的二階差分公式為：因此，在任何一個節(jié)點(diǎn)處，濃度C(x,y,z,t)的一階和二階偏導(dǎo)數(shù)均可用該節(jié)點(diǎn)及其相鄰節(jié)點(diǎn)上濃度值的線性組合表示。12/9/202414一維彌散方程的差分格式設(shè)在無限含水層中，存在一維均勻流場，其滲透速度V=nu，流動方向?yàn)閤軸正向。則一維彌散方程為：12/9/202415一維彌散方程的差分格式對時間區(qū)域[0,T]和空間區(qū)域[0,L]都作等距剖分，設(shè)時間步長為△t，空間步長為△x，把第i個節(jié)點(diǎn)xi處在tn時刻的濃度記為Cin。12/9/202416一維彌散方程的差分格式顯示差分格式將及差分公式中的濃度取為tn時刻的值，便可得到其顯示差分格式：12/9/202417一維彌散方程的差分格式顯示差分格式的穩(wěn)定條件在應(yīng)用顯示差分格式時，剖分的步長必須滿足：12/9/202418一維彌散方程的差分格式隱示差分格式將及差分公式中的濃度取為tn＋1時刻的值，便可得到其隱示差分格式：12/9/202419一維彌散方程的差分格式隱示差分格式的穩(wěn)定條件在應(yīng)用隱示差分格式時，剖分的步長必須滿足：12/9/202420一維彌散方程的差分格式隱示差分格式的矩陣形式12/9/202421一維彌散方程的差分格式Crank-Nicolson差分格式取顯示差分格式和隱示差分格式的平均，便可得到C-N差分格式：12/9/202422一維彌散方程的差分格式Crank-Nicolson差分格式整理后可得：12/9/202423追趕法：方程的統(tǒng)一形式：12/9/202424追趕法：引入邊界條件：當(dāng)i=1時：當(dāng)i=M時：12/9/202425追趕法：令A(yù)·C=δ12/9/202426追趕法：追趕法解方程組的方法：先把對角矩陣A分解為：A=B·W其中：12/9/202427追趕法：由于A=B·W12/9/202428追趕法：由此可得：

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-112/9/202429追趕法：由于A=B·W

→

B·(W·C)=δ又令：W·C=Y其中：Y=(y1,y2,…,yM)T→B·Y=δ12/9/202430追趕法：由以上等式可知：

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值。12/9/202431追趕法：又因?yàn)閃·C=Y，即：12/9/202432追趕法：根據(jù)上述等式，可得：12/9/202433追趕法：追：根據(jù)B·Y=δ，求出bi和wi，并求出yi；趕：根據(jù)W·C=Y，求出Cin+1。12/9/202434二維彌散方程的差分格式若多孔介質(zhì)是各向同性的，流場為均勻的一維流場，滲透速度V=nu，流動方向取x軸的正向，與流速垂直為y方向，對流擴(kuò)散方程可簡化為：12/9/202435二維彌散方程的差分格式將研究區(qū)域網(wǎng)格化12/9/202436二維彌散方程的差分格式建立矩形網(wǎng)差分格式隱式差分格式12/9/202437二維彌散方程的差分格式于是，內(nèi)節(jié)點(diǎn)的差分方程為12/9/202438二維彌散方程的差分格式整理后可得12/9/202439二維彌散方程的差分格式邊界條件初始條件12/9/202440二維彌散方程的差分格式將網(wǎng)格取為正方形，即△x=△y=h，方程就可簡化為：12/9/202441二維彌散方程的差分格式交替方向差分格式（ADI）

ADI法的實(shí)質(zhì)是將二維隱式差分格式化為一維的隱式差分格式。在第n層與第n+1層兩層正中間引入一個過渡層n+1/2。12/9/202442二維彌散方程的差分格式交替方向差分格式（ADI）

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

12/9/202444二維彌散方程的差分格式交替方向差分格式（ADI）后半時間步長：

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

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

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

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

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

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

MT3D99是美國S.S.Papodopulos&Associates公司的Zheng博士設(shè)計的模擬三維地下水溶質(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方法）。

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

12/9/202469FEFLOW

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

12/9/202471VisualModflowPumpingWells

12/9/202472VisualModflowFlowandTransportProperties

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

12/9/202473VisualModflowBo