地下水建模方法和步驟課件

地下水建模方法和步驟中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院2012.8地下水建模方法和步驟中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院1地下水建模方法和步驟1.求解地下水運(yùn)動(dòng)方程的數(shù)值方法2.地下水?dāng)?shù)值模型建模步驟3.建模所需要的基本資料地下水建模方法和步驟1.求解地下水運(yùn)動(dòng)方程的數(shù)值方法2絕大部分?jǐn)?shù)學(xué)模型是無法用解析法求解的，數(shù)值化就是將數(shù)學(xué)模型轉(zhuǎn)化為可解的數(shù)值模型。1.數(shù)值方法1.1有限差分法1.2有限單元法1.3分有限差分法1.4半解析半數(shù)值法1.5邊界元法絕大部分?jǐn)?shù)學(xué)模型是無法用解析法求解的，數(shù)值化就是3(1)有限差分法原理(2)兩種方法建立有限差分方程(3)求解有限差分方程(4)收斂性和穩(wěn)定性概念(5)算例1.1有限差分法(1)有限差分法原理1.1有限差分法4(1)有限差分法的基本原理將連續(xù)的問題離散后求解：方法一．以地下水流基本微分方程及其定解條件為基礎(chǔ)，在滲流區(qū)剖分基礎(chǔ)上，用差商代替微商，將地下水流微分方程的求解轉(zhuǎn)化為差分方程（代數(shù)方程）求解。方法二．在滲流區(qū)剖分的基礎(chǔ)上，直接由達(dá)西定律和水均衡原理，建立各個(gè)均衡區(qū)的水均衡方程，即差分方程。矩形網(wǎng)格多邊形網(wǎng)格(1)有限差分法的基本原理將連續(xù)的問題離散后求解：矩形網(wǎng)格多5網(wǎng)格劃分的基本類型（1）先劃格線，格點(diǎn)位于網(wǎng)格中心均衡網(wǎng)格節(jié)點(diǎn)網(wǎng)格（2）先規(guī)定格點(diǎn)位置，再垂直平分兩相鄰結(jié)點(diǎn)的連線作格線，形成的網(wǎng)格即為水均衡區(qū)網(wǎng)格劃分的基本類型（1）先劃格線，格點(diǎn)位于網(wǎng)格中心均衡網(wǎng)格節(jié)6MODFLOW網(wǎng)格系統(tǒng)MODFLOW網(wǎng)格系統(tǒng)7導(dǎo)數(shù)的有限差商近似導(dǎo)數(shù)的定義

當(dāng)非常小的時(shí)候，有

上式右端項(xiàng)即為f(x)在x0處的差商。這樣定義的差商很容易理解，但不知道用差商代替微商所產(chǎn)生的誤差。下面利用泰勒公式導(dǎo)出差商及其誤差。方法一：差商代替微商(2)有限差分方程建立導(dǎo)數(shù)的有限差商近似導(dǎo)數(shù)的定義當(dāng)非常小的時(shí)候，有上8已知泰勒公式①由A得：

AB②由B得：

稱為f(x)在x0處的一階前向差商，為截?cái)嗾`差。稱為f(x)在x0處的一階后向差商，為截?cái)嗾`差。方法一已知泰勒公式①由A得：AB②由B得：稱9

③由A-B可以得：

④由A+B可以得：AB稱為f(x)在x0處的一階中心差商，為截?cái)嗾`差。稱為f(x)在x0處的二階中心差商，為截?cái)嗾`差。方法一③由A-B可以得：④由A10對(duì)于偏導(dǎo)數(shù)（偏微商），類似可以得到相應(yīng)的差商：方法一(2)有限差分方程建立(續(xù))對(duì)于偏導(dǎo)數(shù)（偏微商），類似可以得到相應(yīng)的差商：方法一(2)有11一維控制方程差分格式顯式差分格式隱式差分格式方法一控制方程網(wǎng)格剖分nx個(gè)一維控制方程差分格式顯式差分格式隱式差分格式方法一控制方程網(wǎng)12取右圖所示得微小六面體。設(shè)與x,y,z,方向?qū)?yīng)得主滲透系數(shù)分別為Kx,Ky,Kz；建立均衡期t時(shí)段內(nèi)，微小均衡六面體的水量守恒方程。方法二：達(dá)西定律和水均衡原理(2)有限差分方程建立(續(xù))取右圖所示得微小六面體。設(shè)與x,y,z,方向?qū)?yīng)得主13基于達(dá)西定律，x，y，z方向流入—流出分別為:t時(shí)段內(nèi)，側(cè)向流入與源匯項(xiàng)導(dǎo)致六面體水量變化量為：ABC(2)有限差分方程建立(續(xù))方法二：達(dá)西定律和水均衡原理DA+B+C+D源匯項(xiàng)基于達(dá)西定律，x，y，z方向流入—流出分別為:t時(shí)段內(nèi)，側(cè)14六面體內(nèi)地下水儲(chǔ)存量的變化為由水均衡原理得三維地下水流動(dòng)方程的有限差分格式(2)有限差分方程建立(續(xù))方法二：達(dá)西定律和水均衡原理六面體內(nèi)地下水儲(chǔ)存量的變化為由水均衡原理得三維地下水流動(dòng)方程15有限差分法:三維（MODFLOW）差商代替微商有限差分法:三維（MODFLOW）差商代替微商16(3)差分方程求解一維顯式差分格式網(wǎng)格個(gè)數(shù)為ni直接求解(3)差分方程求解一維顯式差分格式網(wǎng)格個(gè)數(shù)為ni直接求解17(3)差分方程求解一維隱式差分格式網(wǎng)格個(gè)數(shù)為ni迭代求解方程組PCGSIPSORWHSSAMGGMGMODFLOW(3)差分方程求解一維隱式差分格式網(wǎng)格個(gè)數(shù)為ni迭代求解方程18(4)差分方程的收斂性和穩(wěn)定性截?cái)嗾`差：用差商代替微商時(shí)，地下水流動(dòng)方程產(chǎn)生的誤差為截?cái)嗾`差。收斂性：當(dāng)空間步長(zhǎng)和時(shí)間步長(zhǎng)趨于0時(shí)，有限差分方程的精確解趨于地下水流動(dòng)問題微分方程定解問題的精確解。則稱該差分格式是收斂的。穩(wěn)定性：如果在求解差分方程過程中，某時(shí)間步引入某個(gè)誤差，而在以后的各時(shí)段計(jì)算中，該誤差不再擴(kuò)大，則稱該差分格式是穩(wěn)定的。一維顯示格式的收斂條件和穩(wěn)定條件是：(4)差分方程的收斂性和穩(wěn)定性截?cái)嗾`差：用差商代替微商時(shí)，地19(6)算例：顯式有限差格式設(shè)兩條河流平行、完全切割含水層，含水層等厚、均質(zhì)各向同性。應(yīng)用實(shí)例：河間地塊承壓水流模型(6)算例：顯式有限差格式設(shè)兩條河流平行、完全切割含20步驟：（1）基礎(chǔ)資料的分析（2）概念模型（3）數(shù)學(xué)模型（4）數(shù)值方法及計(jì)算機(jī)程序（5）參數(shù)（6）結(jié)果分析步驟：（1）基礎(chǔ)資料的分析21建立數(shù)學(xué)模型（1）模型概化由所述水文地質(zhì)條件，可以概化為一維承壓水流問題。（2）建立坐標(biāo)系（如圖），將地下水流動(dòng)系統(tǒng)空間結(jié)構(gòu)放在坐標(biāo)系內(nèi)，從而量化各變量的取值范圍。本例，取x-軸原點(diǎn)位于左端河，右側(cè)為正向，設(shè)兩河流間距為L(zhǎng).（3）數(shù)學(xué)模型建立數(shù)學(xué)模型（1）模型概化22差分方程及其解法—顯式格式①將（0—L）分成N等份，1)網(wǎng)格剖分：②取時(shí)間步長(zhǎng)，記（n=0、1、2、3、4……）記，(i=0,1,2,3，4……N)2）建立差分方程：在網(wǎng)格系統(tǒng)中任意取一點(diǎn)設(shè)是問題的解，則在處有記為（i，n）差分方程及其解法—顯式格式①將（0—L）分成N等份，1)23用差商代替微商：

將上述兩式舍去余項(xiàng)，代入方程并記為顯然該式具有截?cái)嗾`差得到顯式格式（續(xù)1）用差商代替微商：將上述兩24引入無量綱變量：將該式子代入得到：

（i=1，2，3，.....N-1），（n=1，2，3，.....）

顯式格式（續(xù)2）引入無量綱變量：將該式子代入得到：（i=1，2253）顯示差分方程的求解計(jì)算各結(jié)點(diǎn)初始時(shí)刻水頭值利用差分方程計(jì)算各結(jié)點(diǎn)t1時(shí)刻水頭值利用邊界條件計(jì)算邊界結(jié)點(diǎn)水頭值重復(fù)2、3步，直到計(jì)算出擬計(jì)算的各個(gè)時(shí)刻的水頭值顯式格式（續(xù)3）3）顯示差分方程的求解計(jì)算各結(jié)點(diǎn)初始時(shí)刻水頭值顯式格式（續(xù)326算例（續(xù)4）在上述模型中，設(shè)L=1000米取空間步長(zhǎng)為200米，時(shí)間步長(zhǎng)為0.25天，分別計(jì)算各節(jié)點(diǎn)各時(shí)刻的水頭值。算例（續(xù)4）在上述模型中，設(shè)L=1000米取空間步長(zhǎng)為20027Time/dayx=0mx=200mx=400mx=600mx=800mx=1000m02010101010100.252012.5101010100.502013.7510.6251010100.752014.53111.25010.15610101.002015.07811.79710.39110.039101.252015.48812.26610.65410.11710算例（續(xù)5）Time/dayx=0mx=200mx=400mx=628Time/dayx=0mx=200mx=400mx=600mx=800mx=1000m02010101010100.252020101010100.502010201010100.75203002010101.0020-10101.252010算例（續(xù)6）如果⊿t=1,則Time/dayx=0mx=200mx=400mx=629(6)算例：隱式格式在上述模型中，設(shè)L=1000米取空間步長(zhǎng)為200米，時(shí)間步長(zhǎng)為0.25天，用隱式差分格式計(jì)算各節(jié)點(diǎn)個(gè)時(shí)刻的水頭值。(6)算例：隱式格式在上述模型中，設(shè)L=1000米取空間步長(zhǎng)30在這個(gè)例子中，解：隱式格式一般方程為于是有根據(jù)初始條件得根據(jù)邊界條件得在這個(gè)例子中，解：隱式格式一般方程為于是有根據(jù)初始條件得根據(jù)31由初始條件和邊界條件由此解得t1時(shí)刻的水頭值為在上述方程中取n=0,可以得到計(jì)算t1時(shí)刻水頭值的方程所以上述方程變成同理，可計(jì)算t2時(shí)刻的水頭值由初始條件和邊界條件由此解得t1時(shí)刻的水頭值為在上述方程中取322.地下水?dāng)?shù)值模型建模步驟2.地下水?dāng)?shù)值模型建模步驟33模擬步驟建立概念模型建立數(shù)學(xué)模型數(shù)值方法及軟件（編程）參數(shù)擬合模擬：模型校正與檢驗(yàn)參數(shù)敏感性分析預(yù)測(cè)模擬軟件模擬步驟建立概念模型軟件34一、概念模型（模型概化）

根據(jù)詳細(xì)的地形地貌、地質(zhì)、水文地質(zhì)、構(gòu)造地質(zhì)、水文地球化學(xué)、巖石礦物、水文、氣象、工農(nóng)業(yè)利用情況等模擬的區(qū)域:含水層類型:潛水（無壓）、承壓、混合、多層維數(shù):一維、二維、三維水流狀態(tài)：穩(wěn)定流/非穩(wěn)定流、飽和流/非飽和流介質(zhì)狀況:均質(zhì)和非均質(zhì)/各向同性和各向異性孔隙/裂隙/雙重介質(zhì)流體的密度差邊界條件和初始條件

必要時(shí)需進(jìn)行一系列的室內(nèi)試驗(yàn)與野外試驗(yàn),以獲取有關(guān)參數(shù),如滲透系數(shù)、彌散系數(shù)、分配系數(shù)、反應(yīng)速率常數(shù)等。

一、概念模型（模型概化）根據(jù)詳細(xì)的地形地貌、地質(zhì)、水35收集研究區(qū)已有水文地質(zhì)資料概念模型邊界概化內(nèi)部結(jié)構(gòu)概化完成模型概化圖地表水體水頭邊界流量邊界斷層接觸邊界隔水邊界水頭邊界流量邊界巖體巖層接觸邊界天然分水嶺據(jù)含水層組類型、結(jié)構(gòu)、巖性含水介質(zhì)地下水運(yùn)動(dòng)狀態(tài)水文地質(zhì)參數(shù)確定層組的均質(zhì)、非均質(zhì)，各向同行、異性，穩(wěn)定、非穩(wěn)定，潛水、承壓水空間分布導(dǎo)水、儲(chǔ)水系數(shù)、主滲透方向目的層與相鄰層關(guān)系均質(zhì)、非均質(zhì)，各向同、異性層、紊流，二、三維時(shí)間概化滲透、儲(chǔ)水系數(shù)，給水度及單位涌水量，含水層分布規(guī)律，地下水流場(chǎng)，水化學(xué)場(chǎng)、溫度場(chǎng)的空間概化平面圖剖面圖研究范圍、主要居民點(diǎn)、標(biāo)志性地形、地貌、地表河流、湖泊、開采井、地下水天然露頭、含水層的各類鉆孔地表地理要素、含水層結(jié)構(gòu)、地質(zhì)結(jié)構(gòu)、地下水水位、各類源匯項(xiàng)及其性質(zhì)收集研究區(qū)已有水文地質(zhì)資料概念模型邊界概化內(nèi)部結(jié)構(gòu)概化完成模36二、數(shù)學(xué)模型

三維地下水流動(dòng)問題控制方程第二類邊界條件第一類邊界條件初始條件二、數(shù)學(xué)模型三維地下水流動(dòng)問題控制方程第二類邊界條件第一類37絕大部分?jǐn)?shù)學(xué)模型是無法用解析法求解的，數(shù)值化就是將數(shù)學(xué)模型轉(zhuǎn)化為可解的數(shù)值模型。三、數(shù)值方法及軟件（或編程）有限差分法有限單元法積分有限差分法半解析半數(shù)值法邊界元法絕大部分?jǐn)?shù)學(xué)模型是無法用解析法求解的，數(shù)值化就是38有限差分法:MOFLOW系列GMS中MODFLOWVisualMODFLOWProcessingMODFLOW有限單元法:FEFLOW積分有限差分法:TOUGH2，TOUGHREACT軟件有限差分法:MOFLOW系列軟件39四、模型參數(shù)含水層參數(shù)：滲透系數(shù)，彈性釋水系數(shù)（重力給水度），孔隙度等源匯項(xiàng)：大氣降水入滲系數(shù)（分區(qū)、數(shù)值）蒸發(fā)排泄系數(shù)地表水體（河流、湖泊、水庫(kù)等）水位、底面高程、底面巖性特征（厚度、滲透系數(shù)等）渠系灌溉入滲系數(shù)人工開采（點(diǎn)狀、面狀）邊界條件初始條件四、模型參數(shù)含水層參數(shù)：滲透系數(shù)，彈性釋水系數(shù)（重力給水度）40參數(shù)的不確定性鉆孔太少，地層資料少，鉆孔多，含水層結(jié)構(gòu)會(huì)發(fā)生變化參數(shù)的不確定性鉆孔太少，地層資料少，鉆孔多，含水層結(jié)構(gòu)會(huì)發(fā)生41五、模擬：模型校正（參數(shù)識(shí)別）將模擬結(jié)果與實(shí)測(cè)結(jié)果比較,進(jìn)行參數(shù)調(diào)整,使模擬結(jié)果在給定的誤差范圍內(nèi)與實(shí)測(cè)結(jié)果吻合。調(diào)參過程是一個(gè)復(fù)雜而辛苦的工作,所調(diào)整的參數(shù)必須符合模擬區(qū)的具體情況。人機(jī)交互與自動(dòng)調(diào)參相結(jié)合。盡管自動(dòng)調(diào)參程序(如PEST),也不能代替人的工作。五、模擬：模型校正（參數(shù)識(shí)別）將模擬結(jié)果與實(shí)測(cè)結(jié)果42五、模擬：模型檢驗(yàn)

模型驗(yàn)證是在模型校正的基礎(chǔ)上,進(jìn)一步調(diào)整參數(shù),使模擬結(jié)果與第二次實(shí)測(cè)結(jié)果吻合,以進(jìn)一步提高模型的置信度。五、模擬：模型檢驗(yàn)?zāi)Ｐ万?yàn)證是在模型校正的基礎(chǔ)上,43六、靈敏度分析

校正后的模型受參數(shù)值的時(shí)空分布、邊界條件、水流狀態(tài)等不確定度的影響。靈敏度分析就是為了確定不確定度對(duì)校正模型的影響程度。六、靈敏度分析校正后的模型受參數(shù)值的時(shí)空分布、邊界44七、預(yù)測(cè)

用校正的參數(shù)值進(jìn)行預(yù)測(cè),預(yù)測(cè)時(shí)需估算未來的水流狀態(tài)。七、預(yù)測(cè)用校正的參數(shù)值進(jìn)行預(yù)測(cè),預(yù)測(cè)時(shí)需估算未來45后續(xù)檢查與模型的重建（完善）

后續(xù)檢查在模擬研究結(jié)束數(shù)年后進(jìn)行。收集新的野外數(shù)據(jù)以確定預(yù)測(cè)結(jié)果是否正確。如果模擬結(jié)果精確,則該模型對(duì)該模擬區(qū)來說是有效的。由于場(chǎng)址的唯一性,故模型只對(duì)該模擬區(qū)有效。后續(xù)檢查應(yīng)在預(yù)測(cè)結(jié)束足夠長(zhǎng)的時(shí)間后進(jìn)行,以便有足夠的時(shí)間發(fā)生明顯的變化。

后續(xù)檢查與模型的重建（完善）后續(xù)檢查在模擬研究結(jié)束數(shù)46模型的再設(shè)計(jì)

一般來說,后續(xù)檢查會(huì)發(fā)現(xiàn)系統(tǒng)性能的變化,從而導(dǎo)致概念模型和模型參數(shù)的修改。一般來說,所有模擬研究都應(yīng)該進(jìn)行到第五步,即校正靈敏度分析。

模型的再設(shè)計(jì)一般來說,后續(xù)檢查會(huì)發(fā)現(xiàn)系統(tǒng)性能的47地裂縫自然地理及水文地質(zhì)條件邊界敏感性分析數(shù)學(xué)模型網(wǎng)格化概念模型介質(zhì)類型、結(jié)構(gòu)特征地下水補(bǔ)、徑、排特征數(shù)值模型數(shù)值模型的建立邊界條件初始條件源匯項(xiàng)降雨入滲地表水體入滲灌溉、渠系蒸發(fā)排泄人工開采介質(zhì)參數(shù)觀測(cè)孔動(dòng)態(tài)擬合流場(chǎng)擬合水均衡對(duì)比分析觀測(cè)孔動(dòng)態(tài)檢驗(yàn)流場(chǎng)檢驗(yàn)水均衡對(duì)比檢驗(yàn)擬合調(diào)參識(shí)別結(jié)果不符合要求修正概念模型小結(jié)：地下水流動(dòng)模型構(gòu)建過程識(shí)別結(jié)果不符合要求修正概念模型模型應(yīng)用水資源量評(píng)價(jià)預(yù)測(cè)開采方案研究環(huán)境地質(zhì)問題地面沉降地下水污染研究大型工程對(duì)環(huán)境的影響模型識(shí)別模型檢驗(yàn)地裂縫自然地理及水文地質(zhì)條件邊界敏感性分析數(shù)學(xué)模型網(wǎng)格化概念483.地下水流數(shù)值模擬資料需求空間展布相關(guān)資料地表高程、分層數(shù)據(jù)邊界位置含水層參數(shù)：滲透系數(shù)，彈性釋水系數(shù)（重力給水度），孔隙度等源匯項(xiàng)：大氣降水入滲系數(shù)（分區(qū)、數(shù)值）蒸發(fā)排泄系數(shù)地表水體（河流、湖泊、水庫(kù)等）水位、底面高程、底面巖性特征（厚度、滲透系數(shù)等）渠系灌溉入滲系數(shù)人工開采（點(diǎn)狀、面狀）：開采井位置、井結(jié)構(gòu)、開采量動(dòng)態(tài)邊界條件：邊界類型、水頭或流量初始條件：統(tǒng)測(cè)水位水位動(dòng)態(tài)觀測(cè)資料：觀測(cè)孔位置、結(jié)構(gòu)、水位時(shí)間變化3.地下水流數(shù)值模擬資料需求空間展布相關(guān)資料49地下水建模方法和步驟中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院2012.8地下水建模方法和步驟中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院50地下水建模方法和步驟1.求解地下水運(yùn)動(dòng)方程的數(shù)值方法2.地下水?dāng)?shù)值模型建模步驟3.建模所需要的基本資料地下水建模方法和步驟1.求解地下水運(yùn)動(dòng)方程的數(shù)值方法51絕大部分?jǐn)?shù)學(xué)模型是無法用解析法求解的，數(shù)值化就是將數(shù)學(xué)模型轉(zhuǎn)化為可解的數(shù)值模型。1.數(shù)值方法1.1有限差分法1.2有限單元法1.3分有限差分法1.4半解析半數(shù)值法1.5邊界元法絕大部分?jǐn)?shù)學(xué)模型是無法用解析法求解的，數(shù)值化就是52(1)有限差分法原理(2)兩種方法建立有限差分方程(3)求解有限差分方程(4)收斂性和穩(wěn)定性概念(5)算例1.1有限差分法(1)有限差分法原理1.1有限差分法53(1)有限差分法的基本原理將連續(xù)的問題離散后求解：方法一．以地下水流基本微分方程及其定解條件為基礎(chǔ)，在滲流區(qū)剖分基礎(chǔ)上，用差商代替微商，將地下水流微分方程的求解轉(zhuǎn)化為差分方程（代數(shù)方程）求解。方法二．在滲流區(qū)剖分的基礎(chǔ)上，直接由達(dá)西定律和水均衡原理，建立各個(gè)均衡區(qū)的水均衡方程，即差分方程。矩形網(wǎng)格多邊形網(wǎng)格(1)有限差分法的基本原理將連續(xù)的問題離散后求解：矩形網(wǎng)格多54網(wǎng)格劃分的基本類型（1）先劃格線，格點(diǎn)位于網(wǎng)格中心均衡網(wǎng)格節(jié)點(diǎn)網(wǎng)格（2）先規(guī)定格點(diǎn)位置，再垂直平分兩相鄰結(jié)點(diǎn)的連線作格線，形成的網(wǎng)格即為水均衡區(qū)網(wǎng)格劃分的基本類型（1）先劃格線，格點(diǎn)位于網(wǎng)格中心均衡網(wǎng)格節(jié)55MODFLOW網(wǎng)格系統(tǒng)MODFLOW網(wǎng)格系統(tǒng)56導(dǎo)數(shù)的有限差商近似導(dǎo)數(shù)的定義

當(dāng)非常小的時(shí)候，有

上式右端項(xiàng)即為f(x)在x0處的差商。這樣定義的差商很容易理解，但不知道用差商代替微商所產(chǎn)生的誤差。下面利用泰勒公式導(dǎo)出差商及其誤差。方法一：差商代替微商(2)有限差分方程建立導(dǎo)數(shù)的有限差商近似導(dǎo)數(shù)的定義當(dāng)非常小的時(shí)候，有上57已知泰勒公式①由A得：

AB②由B得：

稱為f(x)在x0處的一階前向差商，為截?cái)嗾`差。稱為f(x)在x0處的一階后向差商，為截?cái)嗾`差。方法一已知泰勒公式①由A得：AB②由B得：稱58

③由A-B可以得：

④由A+B可以得：AB稱為f(x)在x0處的一階中心差商，為截?cái)嗾`差。稱為f(x)在x0處的二階中心差商，為截?cái)嗾`差。方法一③由A-B可以得：④由A59對(duì)于偏導(dǎo)數(shù)（偏微商），類似可以得到相應(yīng)的差商：方法一(2)有限差分方程建立(續(xù))對(duì)于偏導(dǎo)數(shù)（偏微商），類似可以得到相應(yīng)的差商：方法一(2)有60一維控制方程差分格式顯式差分格式隱式差分格式方法一控制方程網(wǎng)格剖分nx個(gè)一維控制方程差分格式顯式差分格式隱式差分格式方法一控制方程網(wǎng)61取右圖所示得微小六面體。設(shè)與x,y,z,方向?qū)?yīng)得主滲透系數(shù)分別為Kx,Ky,Kz；建立均衡期t時(shí)段內(nèi)，微小均衡六面體的水量守恒方程。方法二：達(dá)西定律和水均衡原理(2)有限差分方程建立(續(xù))取右圖所示得微小六面體。設(shè)與x,y,z,方向?qū)?yīng)得主62基于達(dá)西定律，x，y，z方向流入—流出分別為:t時(shí)段內(nèi)，側(cè)向流入與源匯項(xiàng)導(dǎo)致六面體水量變化量為：ABC(2)有限差分方程建立(續(xù))方法二：達(dá)西定律和水均衡原理DA+B+C+D源匯項(xiàng)基于達(dá)西定律，x，y，z方向流入—流出分別為:t時(shí)段內(nèi)，側(cè)63六面體內(nèi)地下水儲(chǔ)存量的變化為由水均衡原理得三維地下水流動(dòng)方程的有限差分格式(2)有限差分方程建立(續(xù))方法二：達(dá)西定律和水均衡原理六面體內(nèi)地下水儲(chǔ)存量的變化為由水均衡原理得三維地下水流動(dòng)方程64有限差分法:三維（MODFLOW）差商代替微商有限差分法:三維（MODFLOW）差商代替微商65(3)差分方程求解一維顯式差分格式網(wǎng)格個(gè)數(shù)為ni直接求解(3)差分方程求解一維顯式差分格式網(wǎng)格個(gè)數(shù)為ni直接求解66(3)差分方程求解一維隱式差分格式網(wǎng)格個(gè)數(shù)為ni迭代求解方程組PCGSIPSORWHSSAMGGMGMODFLOW(3)差分方程求解一維隱式差分格式網(wǎng)格個(gè)數(shù)為ni迭代求解方程67(4)差分方程的收斂性和穩(wěn)定性截?cái)嗾`差：用差商代替微商時(shí)，地下水流動(dòng)方程產(chǎn)生的誤差為截?cái)嗾`差。收斂性：當(dāng)空間步長(zhǎng)和時(shí)間步長(zhǎng)趨于0時(shí)，有限差分方程的精確解趨于地下水流動(dòng)問題微分方程定解問題的精確解。則稱該差分格式是收斂的。穩(wěn)定性：如果在求解差分方程過程中，某時(shí)間步引入某個(gè)誤差，而在以后的各時(shí)段計(jì)算中，該誤差不再擴(kuò)大，則稱該差分格式是穩(wěn)定的。一維顯示格式的收斂條件和穩(wěn)定條件是：(4)差分方程的收斂性和穩(wěn)定性截?cái)嗾`差：用差商代替微商時(shí)，地68(6)算例：顯式有限差格式設(shè)兩條河流平行、完全切割含水層，含水層等厚、均質(zhì)各向同性。應(yīng)用實(shí)例：河間地塊承壓水流模型(6)算例：顯式有限差格式設(shè)兩條河流平行、完全切割含69步驟：（1）基礎(chǔ)資料的分析（2）概念模型（3）數(shù)學(xué)模型（4）數(shù)值方法及計(jì)算機(jī)程序（5）參數(shù)（6）結(jié)果分析步驟：（1）基礎(chǔ)資料的分析70建立數(shù)學(xué)模型（1）模型概化由所述水文地質(zhì)條件，可以概化為一維承壓水流問題。（2）建立坐標(biāo)系（如圖），將地下水流動(dòng)系統(tǒng)空間結(jié)構(gòu)放在坐標(biāo)系內(nèi)，從而量化各變量的取值范圍。本例，取x-軸原點(diǎn)位于左端河，右側(cè)為正向，設(shè)兩河流間距為L(zhǎng).（3）數(shù)學(xué)模型建立數(shù)學(xué)模型（1）模型概化71差分方程及其解法—顯式格式①將（0—L）分成N等份，1)網(wǎng)格剖分：②取時(shí)間步長(zhǎng)，記（n=0、1、2、3、4……）記，(i=0,1,2,3，4……N)2）建立差分方程：在網(wǎng)格系統(tǒng)中任意取一點(diǎn)設(shè)是問題的解，則在處有記為（i，n）差分方程及其解法—顯式格式①將（0—L）分成N等份，1)72用差商代替微商：

將上述兩式舍去余項(xiàng)，代入方程并記為顯然該式具有截?cái)嗾`差得到顯式格式（續(xù)1）用差商代替微商：將上述兩73引入無量綱變量：將該式子代入得到：

（i=1，2，3，.....N-1），（n=1，2，3，.....）

顯式格式（續(xù)2）引入無量綱變量：將該式子代入得到：（i=1，2743）顯示差分方程的求解計(jì)算各結(jié)點(diǎn)初始時(shí)刻水頭值利用差分方程計(jì)算各結(jié)點(diǎn)t1時(shí)刻水頭值利用邊界條件計(jì)算邊界結(jié)點(diǎn)水頭值重復(fù)2、3步，直到計(jì)算出擬計(jì)算的各個(gè)時(shí)刻的水頭值顯式格式（續(xù)3）3）顯示差分方程的求解計(jì)算各結(jié)點(diǎn)初始時(shí)刻水頭值顯式格式（續(xù)375算例（續(xù)4）在上述模型中，設(shè)L=1000米取空間步長(zhǎng)為200米，時(shí)間步長(zhǎng)為0.25天，分別計(jì)算各節(jié)點(diǎn)各時(shí)刻的水頭值。算例（續(xù)4）在上述模型中，設(shè)L=1000米取空間步長(zhǎng)為20076Time/dayx=0mx=200mx=400mx=600mx=800mx=1000m02010101010100.252012.5101010100.502013.7510.6251010100.752014.53111.25010.15610101.002015.07811.79710.39110.039101.252015.48812.26610.65410.11710算例（續(xù)5）Time/dayx=0mx=200mx=400mx=677Time/dayx=0mx=200mx=400mx=600mx=800mx=1000m02010101010100.252020101010100.502010201010100.75203002010101.0020-10101.252010算例（續(xù)6）如果⊿t=1,則Time/dayx=0mx=200mx=400mx=678(6)算例：隱式格式在上述模型中，設(shè)L=1000米取空間步長(zhǎng)為200米，時(shí)間步長(zhǎng)為0.25天，用隱式差分格式計(jì)算各節(jié)點(diǎn)個(gè)時(shí)刻的水頭值。(6)算例：隱式格式在上述模型中，設(shè)L=1000米取空間步長(zhǎng)79在這個(gè)例子中，解：隱式格式一般方程為于是有根據(jù)初始條件得根據(jù)邊界條件得在這個(gè)例子中，解：隱式格式一般方程為于是有根據(jù)初始條件得根據(jù)80由初始條件和邊界條件由此解得t1時(shí)刻的水頭值為在上述方程中取n=0,可以得到計(jì)算t1時(shí)刻水頭值的方程所以上述方程變成同理，可計(jì)算t2時(shí)刻的水頭值由初始條件和邊界條件由此解得t1時(shí)刻的水頭值為在上述方程中取812.地下水?dāng)?shù)值模型建模步驟2.地下水?dāng)?shù)值模型建模步驟82模擬步驟建立概念模型建立數(shù)學(xué)模型數(shù)值方法及軟件（編程）參數(shù)擬合模擬：模型校正與檢驗(yàn)參數(shù)敏感性分析預(yù)測(cè)模擬軟件模擬步驟建立概念模型軟件83一、概念模型（模型概化）

根據(jù)詳細(xì)的地形地貌、地質(zhì)、水文地質(zhì)、構(gòu)造地質(zhì)、水文地球化學(xué)、巖石礦物、水文、氣象、工農(nóng)業(yè)利用情況等模擬的區(qū)域:含水層類型:潛水（無壓）、承壓、混合、多層維數(shù):一維、二維、三維水流狀態(tài)：穩(wěn)定流/非穩(wěn)定流、飽和流/非飽和流介質(zhì)狀況:均質(zhì)和非均質(zhì)/各向同性和各向異性孔隙/裂隙/雙重介質(zhì)流體的密度差邊界條件和初始條件

必要時(shí)需進(jìn)行一系列的室內(nèi)試驗(yàn)與野外試驗(yàn),以獲取有關(guān)參數(shù),如滲透系數(shù)、彌散系數(shù)、分配系數(shù)、反應(yīng)速率常數(shù)等。

一、概念模型（模型概化）根據(jù)詳細(xì)的地形地貌、地質(zhì)、水84收集研究區(qū)已有水文地質(zhì)資料概念模型邊界概化內(nèi)部結(jié)構(gòu)概化完成模型概化圖地表水體水頭邊界流量邊界斷層接觸邊界隔水邊界水頭邊界流量邊界巖體巖層接觸邊界天然分水嶺據(jù)含水層組類型、結(jié)構(gòu)、巖性含水介質(zhì)地下水運(yùn)動(dòng)狀態(tài)水文地質(zhì)參數(shù)確定層組的均質(zhì)、非均質(zhì)，各向同行、異性，穩(wěn)定、非穩(wěn)定，潛水、承壓水空間分布導(dǎo)水、儲(chǔ)水系數(shù)、主滲透方向目的層與相鄰層關(guān)系均質(zhì)、非均質(zhì)，各向同、異性層、紊流，二、三維時(shí)間概化滲透、儲(chǔ)水系數(shù)，給水度及單位涌水量，含水層分布規(guī)律，地下水流場(chǎng)，水化學(xué)場(chǎng)、溫度場(chǎng)的空間概化平面圖剖面圖研究范圍、主要居民點(diǎn)、標(biāo)志性地形、地貌、地表河流、湖泊、開采井、地下水天然露頭、含水層的各類鉆孔地表地理要素、含水層結(jié)構(gòu)、地質(zhì)結(jié)構(gòu)、地下水水位、各類源匯項(xiàng)及其性質(zhì)收集研究區(qū)已有水文地質(zhì)資料概念模型邊界概化內(nèi)部結(jié)構(gòu)概化完成模85二、數(shù)學(xué)模型

三維地下水流動(dòng)問題控制方程第二類邊界條件第一類邊界條件初始條件二、數(shù)學(xué)模型三維地下水流動(dòng)問題控制方程第二類邊界條件第一類86絕大部分?jǐn)?shù)學(xué)模型是無法用解析法求解的，數(shù)值化就是將數(shù)學(xué)模型轉(zhuǎn)化為可解的數(shù)值模型。三、數(shù)值方法及軟件（或編程）有限差分法有限單元法積分有限差分法半解析半數(shù)值法邊界元法絕大部分?jǐn)?shù)學(xué)模型是無法用解析法求解的，數(shù)值化就是87有限差分法:MOFLOW系列GMS中MODFLOWVisualMODFLOWProcessingMODFLOW有限單元法:FEFLOW積分有限差分法:TOUGH2，TOUGHREACT軟件有限差分法:MOFLOW系列軟件88四、模型參數(shù)含水層參數(shù)：滲透系數(shù)，彈性釋水系數(shù)（重力給水度），孔隙度等源匯項(xiàng)：大氣降水入滲系數(shù)（分區(qū)、數(shù)值）蒸發(fā)排泄系數(shù)地表水體（河流、湖泊、水庫(kù)等）水位、底面高程、底面巖性特征（厚度、滲透系數(shù)等）渠系灌溉入滲系數(shù)人工開采（點(diǎn)狀、面狀）邊界條件初始條件四、模型參數(shù)含水層參數(shù)：滲透系數(shù)，彈性釋水系數(shù)（重力給水度）89參數(shù)的不確定性鉆孔太少，地層資料少，鉆孔多，含水層結(jié)構(gòu)會(huì)發(fā)生變化參數(shù)的不確定性鉆孔太