地下水數(shù)值模擬計算與預測

郭 銳

(遼寧省水文局, 沈陽 110003)

?

地下水數(shù)值模擬計算與預測

郭 銳

(遼寧省水文局, 沈陽 110003)

為了保障阜新地下水的開采量能夠滿足阜新供水的需求，針對阜新凌河地下水狀況進行了預測。以阜新凌河地下水源區(qū)為例，建立地下水流數(shù)學模型，運用該模型對地下水流系統(tǒng)進行預測。應用GMS軟件，用有限差分法的MODFLOW程序進行求解。在1方案總計6.42萬m3/d和2方案總計8.42萬m3/d的開采情景下，分別預測2015年阜新凌河地下水狀況。結果表明所建立的數(shù)學模型能夠較為真實地刻畫研究區(qū)地下水系統(tǒng)的特征，仿真性強。從而得出文中所提出的兩種方案均可行。

凌河; 地下水; 數(shù)值模擬計算; 開采預測

大凌河是遼寧省沿渤海西部的河流，發(fā)源于建昌縣，流經(jīng)朝陽、北票、義縣等市縣，流域面積23 837 km2,河長397 km。流域97%的面積處于遼西低山丘陵區(qū)，荒山禿嶺，水土流失嚴重，河流含沙量大，氣候干旱，是遼寧省內嚴重的資源性缺水地區(qū)，且暴雨集中，洪水來時陡漲陡落。通過查閱相關資料得知[1]，研究區(qū)內的含水層主要有：全新統(tǒng)沖洪積砂礫石孔隙潛水含水層、上更新洪沖積砂礫石承壓含水層和層狀巖類裂隙水。

為了滿足阜新市需水的要求，現(xiàn)需要對現(xiàn)有1號水源井滲渠進行改造，以增加1號水源井供水量；在義縣鐵路橋至錦阜高速公路橋之間新建2處水源井，其中一處擬在鐵路下游600 m新建2水源、另一處在錦阜高速公路橋上游550 m新建3水源，新建及改造的水源達到滿足向凌河水源一泵站供水條件要求，供水能力達到6萬m3/d。為了保障地下水開采量的可行性，需要對凌河地下水狀況進行預測。

1 水文地質概念模型

水文地質概念模型是根據(jù)建模的目的，簡化實際水文地質條件并組織相關數(shù)據(jù)，以便能夠分析地下水系統(tǒng)，為建立地下水流數(shù)值模型提供依據(jù)[2]。通過對研究區(qū)水文地質條件的概化，確定模型的范圍和邊界條件、水文地質結構、地下水流場、水文地質參數(shù)與源匯項，為建立地下水數(shù)值模型奠定基礎。本次數(shù)值模擬采用GMS軟件，用有限差分法的MODFLOW程序求解。

1.1 研究區(qū)范圍

研究區(qū)位于義縣大凌河河谷平原區(qū),北側西起萬佛堂以東、經(jīng)雙臺子—頭溝屯—紅墻子—平房子、東至星星屯；南側西起四方臺、經(jīng)前五里南—八家子、東至西地溝，東西長約18.5 km、南北長約8 km；為一個完整的河谷平原水文地質單元，面積122 km2。

1.2 含水層結構概化

根據(jù)區(qū)內水文地質調查情況，第四系含水層為研究對象，第四系厚度3～14 m，研究區(qū)地下水總的徑流方向是沿著大凌河方向自西向東。地下水在多孔介質中的流動符合質量守恒定律和達西定律；本次模擬只考慮上層第四系含水層，忽略垂向水量交換，故將該區(qū)地下水運動概化為平面二維流；由于介質的非均勻性造成水文地質參數(shù)隨空間變化，體現(xiàn)了系統(tǒng)的非均質性，因而可概化為非均質、各向同性含水層。因此將評價區(qū)含水系統(tǒng)概化為非均質各向同性、平面二維非穩(wěn)定流地下水模型。

1.3 研究區(qū)邊界條件概化

垂向邊界：上部邊界為潛水面，是位置不斷變化的水量交換邊界，有降水入滲和人工開采等；下部以白堊紀砂巖作為隔水底板，概化為隔水邊界。

側向邊界：AB，EF為大凌河河谷平原區(qū)邊界，概化為定流量邊界，AB邊界接受研究區(qū)外的地下潛流，為入流邊界，研究區(qū)范圍內地下水經(jīng)過EF邊界，以地下潛流形式排出研究區(qū)，EF為出流邊界；CD，GH，IJ為支流河谷平原與大凌河干流河谷平原交接處，概化為定流量邊界，均為入流邊界；BC，DE，F(xiàn)G，JA為第四系地層與山區(qū)基巖交界處，概化為零流量邊界(表1)。

表1 研究區(qū)邊界條件

2 地下水流數(shù)學模型及求解

2.1 地下水流數(shù)學模型

根據(jù)水文地質概念模型，將研究區(qū)地下水流系統(tǒng)概化為非均質、各向同性、二維、潛水非穩(wěn)定地下水流系統(tǒng)[3]，運用式(1) 進行計算：

(1)

式中：H0——水頭初始值(m)；?——第一類邊界Γ1上的已知函數(shù)；q——第二類邊界Γ2上的單位寬度側向徑流量(m2/d)；n——邊界Γ2的外法線方向；T——導水系數(shù)(m/d)；S——貯水系數(shù)；w——單位時間在垂向從單位體積含水層中流入或流出的水量(m3)，包括，降水入滲補給地下水量和位于相應單元的抽、水量。

2.2 數(shù)學模型的求解

本次地下水數(shù)值模擬采用了當前國際上利用較廣的GMS軟件。GMS(Groundwater Modeling System)是美國Brigham Young大學環(huán)境模擬研究實驗室開發(fā)的先進的、基于概念模型的地下水環(huán)境模擬軟件。GMS全面包括了模擬地下水流每一個階段所需的工具，如邊界概化、建模、后處理、調參、可視化[4]。

運用基于有限差分法的modflow程序求解以上微分方程的定解問題。根據(jù)評價區(qū)水文地質結構和源匯項的特征，選擇MODFLOW中相應的子程序包來實現(xiàn)地下水流的模擬。

2.3 空間離散

計算區(qū)面積為122 km2，利用GMS的3 D GRID模塊實現(xiàn)網(wǎng)格剖分。共剖分矩形網(wǎng)格單元2 684個，每個單元個面積195 m×235 m。模擬期為2011年1月1日—2012年4月30日，以1天為1個應力期，1天1個時間步長。2012年5月1日—2015年12月31日為預測期，1個月為1個應力期，1天1個時間步長。初始水位流場采用1月1日水文觀測數(shù)據(jù)，對比流場采用4月20日實測值。

2.4 水文地質參數(shù)分區(qū)及初值的確定

用于地下水流模型的水文地質參數(shù)主要有兩類，一類是用于計算各種地下水補排量的參數(shù)和經(jīng)驗系數(shù)，如大氣降水入滲系數(shù)；另一類是含水層的水文地質參數(shù)，主要包括潛水含水層的滲透系數(shù)和給水度[5]。

本次計算所采用的參數(shù)有滲透系數(shù)、給水度、降水入滲系數(shù)等參數(shù)，主要是根據(jù)區(qū)內勘察和現(xiàn)有資料確定，根據(jù)參數(shù)分布規(guī)律和滲流場特征，對各項參數(shù)進行分區(qū)。

降水入滲系數(shù)的大小與包氣帶巖性、地下水埋深、降水量、降水強度、降水前土壤含水量、地形地貌、地面植被等因素有關[6]，前兩者起主導作用。通過對這些因素分析，將研究區(qū)分為4個降水入滲分區(qū)。分區(qū)2為大凌河河漫灘地段，地下水埋深淺、地表巖性粗降水入滲系數(shù)稍大，取值0.25～0.3；分區(qū)1與分區(qū)3為河流階地降水入滲系數(shù)較分區(qū)1稍小，分區(qū)1地面坡度較分區(qū)3大，因此分區(qū)3降水入滲系數(shù)取值較分區(qū)1稍大，分區(qū)1系數(shù)取值0.18～0.22，分區(qū)3取值0.2～0.25；分區(qū)4為義縣縣城范圍，為硬覆蓋區(qū)，降水入滲系數(shù)最小，取值0.15～0.18。

滲透系數(shù)與給水度是通過抽水試驗資料計算得到，參考已有評價結果[7]，結合各種巖性的經(jīng)驗滲透系數(shù)值，綜合確定各計算分區(qū)滲透系數(shù)及給水度的大致取值范圍。

表2 滲透系數(shù)與給水度初值

2.5 研究區(qū)源匯項的確定

研究區(qū)孔隙含水層主要接受河道滲漏補給、地下徑流側向補給、降水入滲補給[8]；地下水主要消耗于人工開采和地下水側向徑流排泄，由于研究區(qū)潛水蒸發(fā)微弱，忽略潛水蒸發(fā)量計算。

1) 河流滲漏補給量。天然條件下，河流側向補給量隨年份和豐枯季節(jié)變化。在模型處理時將其概化為RⅣER邊界，水位根據(jù)義縣水位站監(jiān)測數(shù)據(jù)差值確定，模型RⅣER模塊自動計算地表水與地下水的水力聯(lián)系強度，河流滲漏補給量采用達西定律由模型自動識別計算。

2) 地下水側向補給量。根據(jù)達西定律，各個流量邊界段側向徑流量按如下公式計算：

Qc=K·I·B·M·ΔT

(2)

式中：Qc——地下水側向徑流量(104m3/a)；K——含水層滲透系數(shù)(m/d)；I——水力坡度；B——斷面寬度(m)；M——含水層厚度(m)；ΔT——計算時間(d)。

根據(jù)公式(2)，依據(jù)動態(tài)觀測資料確定邊界流量隨時間段的變化趨勢，根據(jù)不同時期流場選取各段的水力梯度，根據(jù)給定的水文地質參數(shù)初值，計算給定隨時間變化的側向補給(排泄)量。

3) 降水入滲補給量：降水入滲補給是本區(qū)地下水的重要補給源，其入滲量與降水量、潛水水位埋深和包氣帶巖性有關。降水入滲補給量的計算采用公式：

pr=p·α·F×10-1

(3)

式中：pr——降水入滲補給量(萬m3)；p——降水量(mm)；α——降水入滲補給系數(shù)(無因次)；F——計算區(qū)計算面積(km2)。

模型中，以日降水量按照不同分區(qū)，不同降水入滲系數(shù)以Recharge的形式面狀補給地下水。

4) 地下水開采量：扣除豐安用水80萬m3/a,合計取水1 525.1萬m3/a，折成日量為41 784 m3/d。

5) 地下水側向排泄量：計算方法同地下水側向補給量。

6) 水位、降水量采用義縣水文站數(shù)據(jù)。

2.6 數(shù)學模型的識別和驗證

由于對目的層的水文地質條件有比較清楚的認識，各水文地質參數(shù)初值范圍也比較客觀，通過細致的調參擬合，模型識別取得了較好的結果，參數(shù)識別結果見下表。計算流場與識別流場基本相同。

表3 降水入滲系數(shù)識別結果

表4 滲透系數(shù)與給水度識別結果

2.7 模擬計算區(qū)地下水可動用的儲存量的計算

按照區(qū)域水文地質條件，模擬計算區(qū)含水層厚度一般為3～5 m，計算分區(qū)給水度在0.08～0.18之間，模擬計算區(qū)面積為122 km2，用以下公式進行計算：

W=1/3×h×μ×F

(4)

式中：W——地下水可動用的儲存量(m3)；h——地下水面以下含水層厚度(m)；μ——區(qū)域給水度；F——計算面積(m2)。

含水層厚度采用4 m，給水度采用0.16，計算得模擬計算區(qū)地下水可動用的儲存量為2 603萬m3。

2.8 地下水狀況預測

由于研究區(qū)滲透系數(shù)較大，水源區(qū)開采主要是襲奪地表徑流量，河道水是控制研究區(qū)地下水可開采量的主要因素，在義縣站95%水位(59.42 m)保證率下，按連續(xù)枯水年降雨情況對水源井不同開采條件下，對研究區(qū)地下水狀況進行預測。

2.8.1 地下水狀況預測開采1方案 現(xiàn)狀3，4號井關閉,現(xiàn)狀2號井開采0.42萬m3/d、改1井開采1.4萬m3/d、新2井1.9萬m3/d、新3井2.7萬m3/d，總計6.42萬m3/d的開采情景下，預測2015年研究區(qū)地下水狀況。

按照預報方案采用非穩(wěn)定流模型進行模型模擬，可以得出研究區(qū)2015年汛期與非汛期地下水流場演變趨勢。在方案1條件下，研究區(qū)范圍內，在新3井與現(xiàn)2井附近形成了漏斗區(qū)，但是漏斗區(qū)面積很小，汛期(6—9月)0.2 km2，漏斗最大直徑585 m，非汛期(10—4月)0.24 km2漏斗最大直徑780 m，降深小于1.0 m，說明新井開采對研究區(qū)地下水動態(tài)的影響不大。整體上枯水期地下水水位較豐水期略有下降，但是下降不明顯，說明降水不是研究區(qū)地下水的主要補給源，側向入流與河道滲漏控制著研究區(qū)的地下水補給。

在方案1條件下，汛期河道滲漏量78 441 m3/d，占總開采量的74.0%，需襲奪河水0.908 m3/s。非汛期河道滲漏量91 011 m3/d；占總開采量的85.9%，需襲奪河水1.053 m3/d。年內水均衡量為-1 979 m3/d，為負均衡，年均衡量為72.2萬m3，可見區(qū)內補給量小于消耗量，需動用儲存量(按可動用儲存量的60%計算)，年均衡量僅占可動用儲量的4.6%?？梢妳^(qū)內水源地取水6.42萬m3/d是可行的。

2.8.2 地下水狀況預測開采2方案 新建2，3號井與改1號井以設計流量開采，即，改1井開采量1.4萬m3/d、新2井1.9萬m3/d、新3井2.7萬m3/d、現(xiàn)狀2號井開采量0.42萬m3/d，現(xiàn)狀3，4號井以原有開采量1半開采，即：現(xiàn)3井開采1.2萬m3/d、現(xiàn)4井0.8萬m3/d?？傆?.42萬m3/d的開采情景下，2015年研究區(qū)地下水狀況。

按照預報方案采用非穩(wěn)定流模型進行模型模擬，可以得出研究區(qū)2015年汛期與非汛期地下水流場演變趨勢。在方案2條件下，新3號井與現(xiàn)2井范圍形成漏斗，漏斗面積汛期(6—9月)0.18 km2，漏斗最大直徑580 m，非汛期(10—4月)為0.33 km2，漏斗最大直徑877 m，降深小于1.0 m；現(xiàn)3井附近雖然沒有形成閉合漏斗，但是開采區(qū)降深加大，有形成漏斗的趨勢。整體上枯水期地下水水位較汛期略有下降，但是下降不明顯。

汛期河道滲漏量84 561 m3/d，占地下水實際開采量67.1%,需襲奪河水0.979 m3/s；非汛期河道滲漏量99 313 m3/d，占地下水實際開采量78.8%，需襲奪河水1.149 m3/s。年內水均衡量為-2 981 m3/d，為負均衡，年均衡量為108.8萬m3，可見區(qū)內補給量小于消耗量，需動用儲存量，年均衡量占可動用儲量的4.2%?？梢妳^(qū)內水源地取水8.42萬m3/d是可行的。

3 結 論

通過對地下水計算流場和流場的擬合統(tǒng)計分析，說明含水層結構、邊界條件概化、水文地質參數(shù)的選取及源匯項的選取的是合理的，所建立的數(shù)學模型較為真實地刻畫了研究區(qū)地下水系統(tǒng)的特征，仿真性強，可以運用該模型進行地下水流系統(tǒng)的預測。

本文為了保障阜新地下水的開采量能夠滿足阜新供水的需求，針對阜新凌河地下水狀況進行了預測。由于研究區(qū)滲透系數(shù)較大，水源區(qū)開采主要是襲奪地表徑流量，河道水是控制研究區(qū)地下水可開采量的主要因素，在義縣站95%水位(59.42 m)保證率下，按連續(xù)枯水年降雨情況對水源井不同開采條件下，對阜新凌河2015年地下水狀況進行預測，從而得出兩方案可見區(qū)內水源地取水6.42萬m3/d和可見區(qū)內水源地取水8.42萬m3/d均是可行的。

[1] 賈金生,田冰,劉昌明.Visual MODFLOW在地下水模擬中的應用[J].河北農(nóng)業(yè)大學學報,2003,26(2):71-78.

[2] 馬馳,石輝,盧玉東.MODFLOW在西北地區(qū)地下水資源評價中的應用[J].干旱區(qū)資源與環(huán)境,2006(2):89-93.

[3] 丁元芳,遲寶明,易樹平,等.Visual MODFLOW在李官堡水源地水流模擬中的應用[J].水土保持研究,2006,13(5):99-102.

[4] 王金生,王長申,滕彥國.地下水可持續(xù)開采量評價方法綜述[J].水利學報,2006,37(5):525-533.

[5] 李勤奮,方正.上海市地下水可開采量模型計算及預測[J].上海地質,2000(2):36-43.

[6] 魏鴻.開封市地下水資源供求預測研究[J].地下水,2001,23(4):168-169.

[7] 趙鳳雙.地下水開采量與可開采量探析[J].農(nóng)業(yè)與技術,2008,27(5):97-97.

[8] 師幸生.太原市地下水開采量控制目標探討[J].山西水利,2008,24(4):21-23.

Numerical Simulation and Prediction of Groundwater Sources

GUO Rui

(HydrologicalBureauofLiaoningProvince,Shenyang110003,China)

In order to guarantee the quantity of groundwater exploitation to meet the demand for the water rsupply of Fuxin, Ling River groundwater condition in Fuxin was predicted. Groundwater source area of Ling River was set as an example, a groundwater flow model was used to predict groundwater flow system. It was solved using the finite difference method modflow program of applying GMS software. Under a programs of 64 200 m3/d and the program of 84 200 m3/d mining scenarios groundwater conditions of Ling River in Fuxin were predicted in 2015. The results show that the mathematical model can more realistically describe the characteristics of groundwater systems, and the simulation is strong. We can come to a conclusion that the proposed two schemes are feasible based on the prediction results.

Ling River; groundwater; numerical simulation; mining prediction

2014-04-21

2014-06-06

國家自然科學基金重點項目(50979012)

郭銳(1983—),女(滿族),遼寧沈陽人,碩士,中級工程師,研究方向:主要從事水資源方向的研究。E-mail:65893242@qq.com

TV211.1

1005-3409(2015)02-0295-04