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

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

(1.中國電力工程顧問集團華北電力設計院有限公司，北京 100120；2.國網北京市電力公司，北京 100031)

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

1 研究區概況

1.1 自然地理

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

1.2 水文地質條件

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

2 水文地質概念模型

2.1 計算區范圍的確定

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

圖1 模擬區及計算區邊界

2.2 邊界性質的概化

根據研究區地形地貌、水文地質及工程地質條件等，可將研究區的邊界概化見圖2。

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

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

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

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

圖2 水文地質概念模型示意圖

2.3 含水層的概化

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

圖3 含水層頂底板標高示意圖

3 地下水數值模擬

3.1 數學模型

根據計算區水文地質概念模型，對應的數學模型選用一層非均質各向異性三維非穩定流數學模型，所建立的數學模型可表示見式(1)。

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

3.2 網格剖分

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

圖4 模型網格剖分圖

3.3 邊界條件

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

3.4 參數性質的概化

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

表1 模擬模型核正參數

3.5 源匯項的處理

田間入滲補給量則是分不同灌區，根據各灌區的地下水位埋深、包氣帶巖性和灌溉方式及時間等因素影響，確定回滲系數，然后根據收集到各鄉鎮的相關農業用水量計算各單元的農田灌溉回滲量，按各月的用水比例進行分攤，將其概化為單位面狀補給量。

地下水開采用well程序包模擬。區內現狀地下水開采以電廠水源地生產用水、農業開采用水和自來水廠開采用水為主，本次模擬根據實際調查的結果進行分配。其中，電廠水源地生產用水和自來水廠開采用水均按井邊界條件賦值。農業開采用水由于開采井數量大，且開采具有季節性，在每年4～5月、7～8月為灌溉期，為簡化農業開采賦值，又不失準確性，在本次模型中，我們將農業開采量總額比農業開采井分布面積值，在模型中以面源來進行賦值，其值為已扣去灌溉入滲補給量的大小。

潛水蒸發量根據多年平均蒸發量以及模擬區內的水位埋深分布情況將其分配到各個單元格上。

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

3.6 初始水位

采用2014年1月實測地下水等水位線為初始流場，對模擬期后期2016年10月的流場進行擬合。模型按非穩定流運行，通過初始流場的擬合，在調整參數的同時，可以判斷邊界處理的合理性，經反復調試計算，流場擬合結果見圖5。

圖5 流場擬合圖

4 模型識別與校正

4.1 水位擬合檢驗

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

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

圖6 水位擬合圖

4.2 基于模型的地下水均衡

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

4.3 研究區地下水水資源的趨勢分析

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

圖7 開采規模不變情況10年后水位預測圖

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

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

在評價區內，以AS－56觀測孔為監測點，對比兩種情況下，水位變幅情況，見圖9。從圖中，我們可以看出電廠采用城市中水作為生產用水后，水位下降較正常情況下開采，下降幅度小很多。在0～900天范圍內，年水位變幅基本持平，主要是因為這兩年為豐水年，年降雨量在466.7 mm，給當地水資源有了很大補給，之后年份為平水年和枯水年，降雨量較豐水年份有所減少，水位又開始下降。

圖9 兩種情況下AS-56觀測孔內水位變幅圖

5 結論

通過對電廠傍河水源地的水文地質條件進行再分析，建立概念模型，依據水文地質概念模型對應的數學模型選用一層非均質各向異性三維非穩定流數學模型。經過模型的識別和校驗，所建立的地下水數學模型能夠較為真實地反映研究區地下水系統的實際情況。

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

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

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

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

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