郭 銳

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

?

地下水?dāng)?shù)值模擬計(jì)算與預(yù)測(cè)

郭 銳

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

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

凌河; 地下水; 數(shù)值模擬計(jì)算; 開(kāi)采預(yù)測(cè)

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

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

1 水文地質(zhì)概念模型

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

1.1 研究區(qū)范圍

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

1.2 含水層結(jié)構(gòu)概化

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

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

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

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

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

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

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

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

(1)

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

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

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

運(yùn)用基于有限差分法的modflow程序求解以上微分方程的定解問(wèn)題。根據(jù)評(píng)價(jià)區(qū)水文地質(zhì)結(jié)構(gòu)和源匯項(xiàng)的特征，選擇MODFLOW中相應(yīng)的子程序包來(lái)實(shí)現(xiàn)地下水流的模擬。

2.3 空間離散

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

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

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

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

降水入滲系數(shù)的大小與包氣帶巖性、地下水埋深、降水量、降水強(qiáng)度、降水前土壤含水量、地形地貌、地面植被等因素有關(guān)[6]，前兩者起主導(dǎo)作用。通過(guò)對(duì)這些因素分析，將研究區(qū)分為4個(gè)降水入滲分區(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ù)與給水度是通過(guò)抽水試驗(yàn)資料計(jì)算得到，參考已有評(píng)價(jià)結(jié)果[7]，結(jié)合各種巖性的經(jīng)驗(yàn)滲透系數(shù)值，綜合確定各計(jì)算分區(qū)滲透系數(shù)及給水度的大致取值范圍。

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

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

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

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

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

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

(2)

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

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

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

pr=p·α·F×10-1

(3)

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

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

4) 地下水開(kāi)采量：扣除豐安用水80萬(wàn)m3/a,合計(jì)取水1 525.1萬(wàn)m3/a，折成日量為41 784 m3/d。

5) 地下水側(cè)向排泄量：計(jì)算方法同地下水側(cè)向補(bǔ)給量。

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

2.6 數(shù)學(xué)模型的識(shí)別和驗(yàn)證

由于對(duì)目的層的水文地質(zhì)條件有比較清楚的認(rèn)識(shí)，各水文地質(zhì)參數(shù)初值范圍也比較客觀，通過(guò)細(xì)致的調(diào)參擬合，模型識(shí)別取得了較好的結(jié)果，參數(shù)識(shí)別結(jié)果見(jiàn)下表。計(jì)算流場(chǎng)與識(shí)別流場(chǎng)基本相同。

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

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

2.7 模擬計(jì)算區(qū)地下水可動(dòng)用的儲(chǔ)存量的計(jì)算

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

W=1/3×h×μ×F

(4)

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

含水層厚度采用4 m，給水度采用0.16，計(jì)算得模擬計(jì)算區(qū)地下水可動(dòng)用的儲(chǔ)存量為2 603萬(wàn)m3。

2.8 地下水狀況預(yù)測(cè)

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

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

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

在方案1條件下，汛期河道滲漏量78 441 m3/d，占總開(kāi)采量的74.0%，需襲奪河水0.908 m3/s。非汛期河道滲漏量91 011 m3/d；占總開(kāi)采量的85.9%，需襲奪河水1.053 m3/d。年內(nèi)水均衡量為-1 979 m3/d，為負(fù)均衡，年均衡量為72.2萬(wàn)m3，可見(jiàn)區(qū)內(nèi)補(bǔ)給量小于消耗量，需動(dòng)用儲(chǔ)存量(按可動(dòng)用儲(chǔ)存量的60%計(jì)算)，年均衡量?jī)H占可動(dòng)用儲(chǔ)量的4.6%?？梢?jiàn)區(qū)內(nèi)水源地取水6.42萬(wàn)m3/d是可行的。

2.8.2 地下水狀況預(yù)測(cè)開(kāi)采2方案 新建2，3號(hào)井與改1號(hào)井以設(shè)計(jì)流量開(kāi)采，即，改1井開(kāi)采量1.4萬(wàn)m3/d、新2井1.9萬(wàn)m3/d、新3井2.7萬(wàn)m3/d、現(xiàn)狀2號(hào)井開(kāi)采量0.42萬(wàn)m3/d，現(xiàn)狀3，4號(hào)井以原有開(kāi)采量1半開(kāi)采，即：現(xiàn)3井開(kāi)采1.2萬(wàn)m3/d、現(xiàn)4井0.8萬(wàn)m3/d。總計(jì)8.42萬(wàn)m3/d的開(kāi)采情景下，2015年研究區(qū)地下水狀況。

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

汛期河道滲漏量84 561 m3/d，占地下水實(shí)際開(kāi)采量67.1%,需襲奪河水0.979 m3/s；非汛期河道滲漏量99 313 m3/d，占地下水實(shí)際開(kāi)采量78.8%，需襲奪河水1.149 m3/s。年內(nèi)水均衡量為-2 981 m3/d，為負(fù)均衡，年均衡量為108.8萬(wàn)m3，可見(jiàn)區(qū)內(nèi)補(bǔ)給量小于消耗量，需動(dòng)用儲(chǔ)存量，年均衡量占可動(dòng)用儲(chǔ)量的4.2%?？梢?jiàn)區(qū)內(nèi)水源地取水8.42萬(wàn)m3/d是可行的。

3 結(jié) 論

通過(guò)對(duì)地下水計(jì)算流場(chǎng)和流場(chǎng)的擬合統(tǒng)計(jì)分析，說(shuō)明含水層結(jié)構(gòu)、邊界條件概化、水文地質(zhì)參數(shù)的選取及源匯項(xiàng)的選取的是合理的，所建立的數(shù)學(xué)模型較為真實(shí)地刻畫(huà)了研究區(qū)地下水系統(tǒng)的特征，仿真性強(qiáng)，可以運(yùn)用該模型進(jìn)行地下水流系統(tǒng)的預(yù)測(cè)。

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

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

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

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

[4] 王金生,王長(zhǎng)申,滕彥國(guó).地下水可持續(xù)開(kāi)采量評(píng)價(jià)方法綜述[J].水利學(xué)報(bào),2006,37(5):525-533.

[5] 李勤奮,方正.上海市地下水可開(kāi)采量模型計(jì)算及預(yù)測(cè)[J].上海地質(zhì),2000(2):36-43.

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

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

[8] 師幸生.太原市地下水開(kāi)采量控制目標(biāo)探討[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

國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目(50979012)

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

TV211.1

1005-3409(2015)02-0295-04