基于有限元法的傍河水源地地下水資源評價及預(yù)測

雷 磊，龍子華，閆曉君，徐曉青

(1.中國電力工程顧問集團(tuán)華北電力設(shè)計院有限公司，北京 100120；2.國網(wǎng)北京市電力公司，北京 100031)

模型計算指的是進(jìn)行地下水流或溶質(zhì)運(yùn)移正反演計算，常用的方法主要是有限差分法、有限元法、邊界元法等。FEFLOW軟件是一套專門用于孔隙介質(zhì)中的地下水流動的三維有限元數(shù)值模擬軟件，近年來，該軟件在地下水水資源評價得到了廣泛應(yīng)用。經(jīng)多年的開采實(shí)踐證明，傍河取水是保證長期穩(wěn)定供水的有效途徑，本文研究的某傍河電廠水源地已運(yùn)行了約30年，周邊地下水位下降情況比較嚴(yán)峻，已影響到電廠正常的生產(chǎn)用水。在對水源地水文地質(zhì)條件現(xiàn)狀再分析的基礎(chǔ)上建立水文地質(zhì)概念模型和數(shù)學(xué)模型，在此基礎(chǔ)上運(yùn)用地下水模擬軟件FEFLOW建立研究區(qū)地下水流模型，并基于所建模型評價了當(dāng)前地下水的水資源和進(jìn)行了趨勢分析。

1 研究區(qū)概況

1.1 自然地理

研究區(qū)位于宣化縣城西北，為一條帶狀山間盆地，見圖1。洋河自西向東從盆地流過(現(xiàn)基本已干涸)，將盆地地形形態(tài)進(jìn)一步分為三部分，即沿洋河為河床、河漫灘及一、二級沖積階地；洋河以北為沖洪積傾斜平原；以南主要為坡洪積臺地。洋河階地為上迭階地，沿河斷續(xù)分布，階面較窄，高出河床5～10 m，沖洪積傾斜平原(沖洪積裙)和坡洪積臺地前緣高出河床約5～20 m左右。電廠位于洋河北岸二級階地，其配套水源地分二期，一期水源地位置在沙嶺子鎮(zhèn)至樣臺一帶洋河北岸的河漫灘地段，二期水源地位于洋河南岸的西前所村至古樹營村一帶洋河河漫灘地中。

1.2 水文地質(zhì)條件

研究區(qū)地貌成因?qū)儆跇?gòu)造侵蝕堆積盆地區(qū)，含水層巖性為中細(xì)砂、圓礫和卵石與灰黃色和褐黃色含礫粘性土互層，最厚可達(dá)160 m，為地下水賦存了良好的條件，地下水類型為弱承壓水～潛水。沿洋河兩岸地區(qū)富水程度極強(qiáng)，單位涌水量在65～90(m3/h·m)；沿洋河向南、向北富水性逐漸變?nèi)酰鶐r山區(qū)由于基巖巖性及裂隙不發(fā)育，富水程度較差；地下水溶解性總固體含量基本上小于1 g/l，地下水水化學(xué)類型一般為HCO3－Ca－Mg型水，個別為HCO3－Ca－Mg－Na型水。水源地運(yùn)行初期，研究區(qū)地下水的主要補(bǔ)給來源是洋河地表水，河床一個水文年平均滲漏量為1.154 m3/s，區(qū)域內(nèi)潛水水位埋深一般0～3 m，經(jīng)過十幾年地下水的開采，研究區(qū)水文地質(zhì)條件已發(fā)生了變化?，F(xiàn)洋河基本斷流，洋河已改為景觀河，河底和兩側(cè)都鋪設(shè)水泥石板，阻斷了地表水滲透補(bǔ)給地下水，現(xiàn)區(qū)域內(nèi)潛水水位埋深一般大于20 m，地下水補(bǔ)給現(xiàn)在主要為地下水側(cè)向徑流，其次為田間灌溉及大氣降水入滲補(bǔ)給；地下水主要以側(cè)向徑流、人工開采等方式排泄。地下水流向自西北流向東南。

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

2.1 計算區(qū)范圍的確定

從地下水流動理論出發(fā)，第四系孔隙含水系統(tǒng)滲流場數(shù)值模擬的范圍應(yīng)取至流動系統(tǒng)的自然邊界，即基巖與第四系的交界處。本次地下水資源評價的模型范圍東北與西南邊界為地表自然分水嶺，西北邊界為清水河和洋河的交匯處，東南邊界為盤長河與洋河的交匯處，模型區(qū)計算范圍與研究區(qū)面積一致，見圖1。

圖1 模擬區(qū)及計算區(qū)邊界

2.2 邊界性質(zhì)的概化

根據(jù)研究區(qū)地形地貌、水文地質(zhì)及工程地質(zhì)條件等，可將研究區(qū)的邊界概化見圖2。

BC、DA邊界：概化為零通量邊界，此兩處邊界為地表自然分水嶺，與地下水流場的流線平行，為隔水邊界。

CD邊界：概化為流量邊界(第二類排泄邊界)，即在水源地的影響范圍之外，地下水的流動仍保持原來的流動狀態(tài)，向外流出，但由于開采量增大，會襲奪一部分水量流入井內(nèi)，從而使側(cè)向流出減少。

AB邊界：概化為流量邊界(第二類補(bǔ)給邊界)，補(bǔ)給來源主要為上游地下水的側(cè)向流入。

在垂向上模型的上邊界為潛水水面，接受大氣降水、田間灌溉滲漏入滲補(bǔ)給，以人工開采、大氣蒸發(fā)來排泄地下水。下邊界取第四系單一結(jié)構(gòu)潛水含水層的下界面為相對隔水層。

圖2 水文地質(zhì)概念模型示意圖

2.3 含水層的概化

研究區(qū)內(nèi)有供水意義的含水層主要為第四系全新統(tǒng)和上更新統(tǒng)松散巖類孔隙水，含水層巖性主要為卵礫石和中細(xì)砂，含水層之間并無連續(xù)穩(wěn)定的隔水層，相互間水力聯(lián)系密切，并且區(qū)域內(nèi)地下水皆為統(tǒng)一開采，故可將計算區(qū)含水層概化為非均質(zhì)各向異性的統(tǒng)一的潛水含水層。本次模擬模型的底部邊界根據(jù)巖性和第四系含水層的厚度以及下部滲透系數(shù)較差的粘土層厚度確定，模型的頂部根據(jù)實(shí)測地面標(biāo)高插值得出，三維建模見圖3。

圖3 含水層頂?shù)装鍢?biāo)高示意圖

3 地下水?dāng)?shù)值模擬

3.1 數(shù)學(xué)模型

根據(jù)計算區(qū)水文地質(zhì)概念模型，對應(yīng)的數(shù)學(xué)模型選用一層非均質(zhì)各向異性三維非穩(wěn)定流數(shù)學(xué)模型，所建立的數(shù)學(xué)模型可表示見式(1)。

式中：Ω 為模擬區(qū)域；L2；Kxx、Kyy和Kzz分別為X、Y和Z方向的滲透系數(shù)(L/T)，Kxx=Kyy；H為水頭值(L)；H0為初始地下水位；ε為源匯項(xiàng)(L/T)；S為給水度；n為邊界面的外法線方向；Г為AB、CD側(cè)邊界；B為底邊界和BC、DA側(cè)邊界。

3.2 網(wǎng)格剖分

對計算區(qū)在空間上的離散包括平面上的網(wǎng)格剖分和垂向上的分層。本次研究區(qū)的面積為230 km2，研究區(qū)剖分的單元格共有57005個，節(jié)點(diǎn)57240，剖分單元格為三角形，平均邊長約120 m，面積約6000m2。在電廠水源地和自來水廠水源地區(qū)域進(jìn)行了網(wǎng)格加密，加密三角形平均單元格邊長約0.2 m，如圖4 所示，垂向上概化為一層。

圖4 模型網(wǎng)格剖分圖

3.3 邊界條件

將底邊界處理為隔水邊界；上邊界作為開放邊界，考慮入滲及蒸發(fā)；將側(cè)邊界條件概化為流量邊界和零通量邊界。

3.4 參數(shù)性質(zhì)的概化

模型中每個單元格的含水層的厚度、地形高程及底板高程由ARCGIS軟件自動插值獲得。根據(jù)滲流場的特征，將研究區(qū)劃分為5個參數(shù)區(qū)，每個參數(shù)區(qū)根據(jù)利用已有抽水試驗(yàn)資料計算滲透系數(shù)和給水度，結(jié)合不同巖性的經(jīng)驗(yàn)值，確定不同巖性的水平、垂向滲透系數(shù)和給水度，模型擬合、核正后參數(shù)見表1。模型參數(shù)分區(qū)及滲透系數(shù)賦值見圖2(a)。降雨入滲系數(shù)依據(jù)地形地貌和出露巖性確定。

表1 模擬模型核正參數(shù)

3.5 源匯項(xiàng)的處理

田間入滲補(bǔ)給量則是分不同灌區(qū)，根據(jù)各灌區(qū)的地下水位埋深、包氣帶巖性和灌溉方式及時間等因素影響，確定回滲系數(shù)，然后根據(jù)收集到各鄉(xiāng)鎮(zhèn)的相關(guān)農(nóng)業(yè)用水量計算各單元的農(nóng)田灌溉回滲量，按各月的用水比例進(jìn)行分?jǐn)?，將其概化為單位面狀補(bǔ)給量。

地下水開采用well程序包模擬。區(qū)內(nèi)現(xiàn)狀地下水開采以電廠水源地生產(chǎn)用水、農(nóng)業(yè)開采用水和自來水廠開采用水為主，本次模擬根據(jù)實(shí)際調(diào)查的結(jié)果進(jìn)行分配。其中，電廠水源地生產(chǎn)用水和自來水廠開采用水均按井邊界條件賦值。農(nóng)業(yè)開采用水由于開采井?dāng)?shù)量大，且開采具有季節(jié)性，在每年4～5月、7～8月為灌溉期，為簡化農(nóng)業(yè)開采賦值，又不失準(zhǔn)確性，在本次模型中，我們將農(nóng)業(yè)開采量總額比農(nóng)業(yè)開采井分布面積值，在模型中以面源來進(jìn)行賦值，其值為已扣去灌溉入滲補(bǔ)給量的大小。

潛水蒸發(fā)量根據(jù)多年平均蒸發(fā)量以及模擬區(qū)內(nèi)的水位埋深分布情況將其分配到各個單元格上。

洋河地表水對地下水的補(bǔ)給量在本模型中可不考慮。由于評價區(qū)內(nèi)地表水體(洋河)河道已改為景觀河道，河槽下及兩側(cè)均鋪設(shè)水泥石板，造成地表水與地下水聯(lián)系不密切，且洋河上游建設(shè)有橡皮壩，造成洋河常年斷流，地表水入滲量極小，可忽略。

3.6 初始水位

采用2014年1月實(shí)測地下水等水位線為初始流場，對模擬期后期2016年10月的流場進(jìn)行擬合。模型按非穩(wěn)定流運(yùn)行，通過初始流場的擬合，在調(diào)整參數(shù)的同時，可以判斷邊界處理的合理性，經(jīng)反復(fù)調(diào)試計算，流場擬合結(jié)果見圖5。

圖5 流場擬合圖

4 模型識別與校正

4.1 水位擬合檢驗(yàn)

為了驗(yàn)證對地下水系統(tǒng)模擬是正確的，研究區(qū)內(nèi)設(shè)有3個水位觀測點(diǎn)(圖5)，分別位于洋河北岸和南岸，分布較均勻，具有較好的代表性。通過不斷地改變水文地質(zhì)參數(shù)，重復(fù)計算，直到整個模擬期內(nèi)觀測水位與計算水位的水位差在精度要求的范圍內(nèi)。

從圖5可以看出，觀測水位和計算水位在計算區(qū)范圍內(nèi)變化趨勢總體保持一致，即上升和下降趨勢一致。依據(jù)計算區(qū)2014年1月至2016年10月地下水觀測水位資料，時間步長以天為單位。研究區(qū)內(nèi)三個觀測孔的擬合水位過程如圖6。從圖6可以看出，整體曲線擬合很好，通過以上模型的識別和校驗(yàn)，所建立的地下水?dāng)?shù)學(xué)模型能夠較為真實(shí)地反映研究區(qū)地下水系統(tǒng)的實(shí)際情況。

圖6 水位擬合圖

4.2 基于模型的地下水均衡

FEFLOW軟件提供了記錄計算過程的輸出文件，本次模擬的總補(bǔ)給量為1.35×104m3/d，總排泄量為1.85×104m3/d，均衡差為－0.50×104m3/d，呈現(xiàn)負(fù)均衡，這與該地區(qū)地下水位逐年下降的趨勢吻合。

4.3 研究區(qū)地下水水資源的趨勢分析

為了評價研究區(qū)水資源今后趨勢變化，我們選取當(dāng)?shù)?996～2005年這10年降雨量作為本次模型降水量的輸入值。因?yàn)檫@10年不僅包括了豐水期、平水期及枯水期不同年份的組合，而且年平均降水量偏小，較能體現(xiàn)今后降水量減少的趨勢。本文從兩種情況來對比分析，一種是保持現(xiàn)開采規(guī)模不變的情況，另一種是電廠水源地停采地下水采用城市中水作為生產(chǎn)用水或部分采用中水部分采用地下水的情況。10年后，其水位預(yù)測見圖7、圖8。

圖7 開采規(guī)模不變情況10年后水位預(yù)測圖

圖8 電廠水源地停止開采情況10年后水位預(yù)測圖

從水位預(yù)測圖上可以看出，如果研究區(qū)開采規(guī)模不變，10年后，研究區(qū)水位下降了約12 m，相當(dāng)于水位平均下降1.2 m/a。如果電廠停采地下水，采用城市中水作為生產(chǎn)用水，10年后，水位下降了約8 m，相當(dāng)于水位平均下降0.8 m/a。說明，電廠采用中水，緩解了水位下降的速率，但并沒有改變研究區(qū)水位持續(xù)下降。這是因?yàn)檠芯繀^(qū)內(nèi)人工開采為地下水的主要排泄方式，其中，電廠生產(chǎn)用水只占用了其中之一，約17%；農(nóng)業(yè)灌溉用水比例較大，灌溉持續(xù)用水依然會導(dǎo)致區(qū)域水位下降。但電廠采用中水后，評價區(qū)內(nèi)水源地周邊漏斗有所恢復(fù)，水位下降速率變緩，對整個研究區(qū)的地下水環(huán)境也有所改善。

在評價區(qū)內(nèi)，以AS－56觀測孔為監(jiān)測點(diǎn)，對比兩種情況下，水位變幅情況，見圖9。從圖中，我們可以看出電廠采用城市中水作為生產(chǎn)用水后，水位下降較正常情況下開采，下降幅度小很多。在0～900天范圍內(nèi)，年水位變幅基本持平，主要是因?yàn)檫@兩年為豐水年，年降雨量在466.7 mm，給當(dāng)?shù)厮Y源有了很大補(bǔ)給，之后年份為平水年和枯水年，降雨量較豐水年份有所減少，水位又開始下降。

圖9 兩種情況下AS-56觀測孔內(nèi)水位變幅圖

5 結(jié)論

通過對電廠傍河水源地的水文地質(zhì)條件進(jìn)行再分析，建立概念模型，依據(jù)水文地質(zhì)概念模型對應(yīng)的數(shù)學(xué)模型選用一層非均質(zhì)各向異性三維非穩(wěn)定流數(shù)學(xué)模型。經(jīng)過模型的識別和校驗(yàn)，所建立的地下水?dāng)?shù)學(xué)模型能夠較為真實(shí)地反映研究區(qū)地下水系統(tǒng)的實(shí)際情況。

本次模擬的總補(bǔ)給量為1.35×104m3/d，總排泄量為1.85×104m3/d，均衡差為－0.50×104m3/d，呈現(xiàn)負(fù)均衡。通過對兩種開采情況地下水水資源趨勢變化可以看出，如果保持現(xiàn)開采量不變，10年后地下水水位持續(xù)下降，由于將來研究區(qū)內(nèi)農(nóng)業(yè)的持續(xù)發(fā)展，水位的下降速率可能比預(yù)測更加嚴(yán)峻，電廠生產(chǎn)用水更加得不到保障。因此，從電廠安全用水考慮，電廠將來可采用城市中水，滿足生產(chǎn)用水。并且，根據(jù)電廠采用城市中水預(yù)測的結(jié)果，電廠減少開采地下水，水位下降的速率變緩，遇到豐水期時，地下水水位還能恢復(fù)，可有效地對區(qū)域內(nèi)地下水環(huán)境進(jìn)行改善。

[1]丁繼紅，周德亮，馬生忠.國外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢[J].勘查科學(xué)技術(shù)，2002，(1).

[2]劉國義，等.傍河取水水源地地下水資源評價及實(shí)例研究[J].安全與環(huán)境工程，2008，15(2).

[3]郇心善，等.沙嶺子電廠水源地供水水文地質(zhì)勘測報告(詳勘階段)[R].北京：水利電力部電力規(guī)劃設(shè)計院，1987.

[4]馬小波.基于有限差分法的康灘水源地地下水資源評價[J].地下水，2016，38(6).