第五章 地下水污染預測

第五章地下水污染預測

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

近似解析法（了解）

數(shù)值法（熟悉）

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2/5/202369FEFLOW

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