大清河淀西平原地下水數值模型研究

李 巍，郭抒燕，薛麗娟，丁曉夢

（1.水利部海河水利委員會，天津 300170；2.天津市龍網科技發展有限公司，天津 300170）

我國華北地下水水資源短缺、地下水超采問題由來已久[1-4]，出現一系列環境、生態問題[5-7]。數值模型是研究地下水問題的有效手段之一，韓忠等人通過對MODFLOW 的二次開發，解決了面狀補給的問題，豐富了源匯項包，提高了模型運行效率，形成優化后的MODFLOW-RAW 主程序、RAW 面狀源匯項包和改進后的RCH、EVT 程序包[8，9]?；诖笄搴拥砦髌皆叵滤r業開采數據以面狀為主的特點，為滿足地區發展和國家工農業發展規劃的規定及地下水資源的可持續發展利用，有必要基于MODFLOW-RAW 對大清河淀西平原地下水資源做出新的評價。而建立大清河淀西平原地下水數值模型是研究、管理地下水的重要技術手段和有效方法。

1 研究區概況

1.1 水文地質

大清河淀西平原區主要由山前沖洪積傾斜平原、沖洪積扇和中部沖積平原構成，山前沖洪積平原、沖洪積扇呈扇狀交錯分布于山前，含水砂層主要由砂礫石、粗砂、中砂、中細砂等各類砂、砂礫石組成，含水層下部無連續隔水層，垂向水力聯系良好，常因下部含水層的開采而疏干，因此在山前常將這一含水層組與下部承壓含水層組看作是統一含水層組。

平原區地下水系統的上游為大清河沖洪積扇子系統，分為4個亞級，包括拒馬河沖洪積扇系統、瀑河-漕河沖洪積扇系統、唐河-界河沖洪積扇系統、大沙河-磁河沖洪積扇系統；下游為大清河古河道帶子系統。

1.2 含水組劃分

由區域水文地質剖面圖以及查閱歷史資料可知，大清河淀西平原分4個含水巖組。

（1）第一含水層組。底界面埋深一般小于50 m，在山前平原沖洪積扇和扇間、扇前地帶具有不同的水文地質特征。在沖洪積扇地區，含水層粒度大、厚度30～50 m、垂向連續性強，屬單層或雙層結構，透水性強，導水系數多大于5 000 m2/d，單位涌水量為20～30 m3（/h·m）。含水層直接裸露，或者被薄層砂質黏土覆蓋，具有強入滲補給和儲存條件，且山前地區與山區河谷含水體相連，具有側向徑流補給條件。

（2）第二含水層組。底界埋深一般120～210 m，在山前平原，與第一含水層組之間缺乏穩定的隔水層，二者之間具有較好的水力聯系。自西向東發育2～3 套中細砂～中粗砂～砂礫石韻律層，含水層透水性與導水性均比第一含水層組強。目前，由于2個含水層組混合開采，人工加強了二者之間的水力聯系。

（3）第三含水層組。底界埋深250～310 m，幾乎全部為淡水。山前平原含水層呈扇狀、扇裙狀展布，由3～4套中細砂～中粗砂～礫石韻律層組成，下段含水層更受到不同程度的風化。典型平原區南部扇體內單井涌水量20～50 m3（/h·m），扇間地帶單井涌水量10～20 m3（/h·m）。在大型扇體內部與上覆第二含水層組之間無連續分布的隔水層，二者水力聯系良好，在其他地段一般都有5～10 m 的黏土或砂質黏土分布，水力聯系變弱。

（4）第四含水層組。底界為第四系基底，由于少有勘探井揭穿本含水巖組底板，底板埋深及地下水厚度不清。該層地下水水力性質均為承壓水，TDS較第三含水層組略有增高。

1.3 地下水開發利用

由于多年的開采，形成了保定市區漏斗、一畝泉漏斗和高蠡清漏斗。保定市區漏斗形成于1974 年低水期，但季節性非常明顯，進入20 世紀80 年代以后，就轉為常年性漏斗，漏斗中心地下水位埋深及影響面積從1980 年的15.14 m 和136.3 km2發展到2000 年的32.0 m 和241.0 km2，20 a 間最大埋深增加了16.83 m，影響面積增大了104.7 km2。保定市區漏斗的發展主要原因是市區自備井的無序取水。

平原區內包括四大灌區，分別為房淶涿灌區、易水灌區、唐河灌區、沙河灌區。灌區多引水庫或地表河流水系水進行農業灌溉，但是近年來，由于降雨減小地表徑流被上游水庫攔截等因素影響，灌區采用地下水的比例有大幅度增大，這對地下水模型的建立與驗證產生一定影響。

2 基于優化MODFLOW的淀西平原地下水數值模型

以下通過建立大清河淀西平原地下水數值模型，并結合南水北調受水區調水方案，預測在不同地下水開發利用條件下的地下水位、水量變化，為指導地下水開發利用、地下水管理提供科學支撐。

2.1 MODFLOW 的優化

針對大區域模型中源匯項數據處理的難題，韓忠等人在MODFLOW2005 源程序的基礎上進行了優化，開發了MODFLOW-RAW，并利用模型降雨、蒸發和開采數據處理的新方法，簡化模型的數據處理過程，縮短模型運行時間[8，9]。其中，包括降雨子程序包（RCH）、蒸發子程序包（EVT）的處理方法改進和面狀源匯項程序包（RAW）的開發，MODFLOW-RAW 系列程序包不僅通過改變數據存儲方式減少了人工處理數據的工作量、加快了模型數據的讀取速度，且壓縮后的源匯項文件更有利于數據的修改及模型的調試。

同時，借助RAW 程序包具可同時處理最多4 個面狀源匯項的平行處理能力，可將區域農業開采、灌區開采、灌區回灌和白洋淀入滲分別計入，使模型的前、后處理更加清晰、便捷。

2.2 水文地質概念模型

針對淀西平原區第四系松散含水層組建模。該區域地質條件復雜，目前的水文地質工作基礎尚不足以詳細刻畫控制全區的三維物理滲流結構，能夠取得的較為詳細的地質調查成果為前述的4層結構的含水層系統，對區域地下水系統結構已經做了相當的概化。本次地下水數值模型對區域地下水含水層組進行進一步的概化，將第一含水層組和第二含水層組合并作為淺層地下水系統，主要模擬潛水循環運動規律，同時也包含與潛水含水層緊密相關的微承壓含水層的循環在內；將第三層和第四層合并作為統一的深層地下水系統，主要模擬區域第三層承壓水的循環運動規律，同時也包含循環量較小的第四層承壓水的循環在內。

含水層厚度大、分布廣、地下水流以水平方向為主，垂向運動較弱但不可忽視，地下水運動空間概化為三維流；源匯項和水位均隨時間發生變化，表現為明顯的非穩定流特征；水文地質參數隨空間變化，體現非均質性，在同一點的水平方向上，參數無明顯方向性，可視為各向同性。

綜上，由于第一含水層組和第二含水層組水力聯系緊密，且實際開采多以混采井居多，故將第一含水層組和第二含水層組合并概化為第一含水層，將第三含水層組作為第二含水層，第四含水層組作為第三含水層。其中，第一含水層為潛水含水層，第二、三含水層為承壓含水層。綜合考慮地下水開發利用層位的特點，參照含水層的發育程度、含水層的滲透性、地下水水力性質、地下水動態特征等因素，將模型在垂向上概化為三層的非均質各項同性、空間三維結構的非穩定地下水流動系統。

2.3 地下水數學模型

根據水文地質概念模型，對非均質各向同性、空間三維結構的非穩定地下水流動系統，可用如下方程的定解問題來描述：

式中：Ω為滲流區域；h為含水層的水位標高（m）；Kx、Ky、Kz分別為水平和垂向滲透系數（m/d）；ε為含水層的源匯項（1/d）；h0為含水層的初始水位分布（m）；S為含水介質的儲水率（1/m）；Kn為邊界面法向方向的滲透系數（m/d）；τ1為滲流區域的側向和下邊界；q為τ1邊界的流量（m/d），流入為正，流出為負，隔水邊界為0；n為邊界面的法線方向；Q為流量（m3/s）；K為滲透系數（m/d）；H為外擴邊界水頭（m）；H0為模型邊界水頭（m）；a為單個單元格（網格）長度（m）；b為單個單元格（網格）寬度（m）；L為滲流路徑長度（m）；t為時間（d）。

2.4 地下水數值模型

2.4.1 模型離散

（1）時間離散。模擬期為2000—2010 年，為反映年內的周期變化，以月單位作為一個應力期，每個應力期為一個時間步長。

（2）空間離散。根據區域的含水層結構、邊界條件和地下水流場特征，將模擬區第一層剖分為220行、160 列的1 km×1 km 規則網格，第一層有效單元共12 554個，下部第二、三層7 983個。

2.4.2 定解條件

（1）初始水位。以2000 年1 月的地下水水位資料為基礎，參考長序列數據，繪制初始流場，采用內插法和外推法獲取各含水層的初始水位。

（2）邊界條件。模型根據水文地質條件，可將區域劃分為2 類邊界，即流量邊界和通用水頭邊界。區域西部山前側向流入邊界，月平均徑流量為2.11×106m3；東部的邊界以行政區劃分，定為通用水頭邊界。模擬邊界，如圖1所示。

圖1 淀西平原模擬邊界示意

2.4.3 源匯項處理

以區縣為單位，蒸發量采用改進的MODFLOWRAW 中的EVT 蒸發程序包處理，根據蒸發極限埋深、蒸發速率及模擬水位確定。降雨入滲采用改進的MODFLOW-RAW 中的RCH 程序包處理。灌區開采和回灌以及其他開采包括生活工業用水均采用面狀補排處理方式。開采山前側向流量采用GMS中的WELL井流程序包處理。具體應用中，在每個網格中心點上設置一個注水井，然后將線狀的流量數據都分配到網格中心點上，將側流量以點井的形式導入模型中。這種源匯處理方式效率高，便于調參過程中源匯數據的修改。補給排泄項對應的程序包，詳見表1。

表1 模型中源匯項的數據對照

2.4.4 模型識別與校驗

通過模擬水位與實測水位的對比來檢驗長序列模型的可靠性。對研究區內19 個觀測點共計2 267個數據進行了誤差分析，其平均絕對誤差為2 m。其中，11 個觀測孔的平均絕對誤差小于1.5 m，4 個觀測孔的平均絕對誤差介于1.5~2 m，4個觀測孔的平均絕對誤差大于2 m。綜上，模擬區內誤差2 m 之內的觀測孔約占80%，整體擬合效果較好。擬合誤差分布，詳見表2。

表2 擬合誤差分布

2.5 模型應用

基于擬合驗證后的大清河淀西平原地下水數值模型，分別建立現狀和南水北調實施后2 種開采條件下的10 a預測模型，通過方案對比，定量分析地下水調控對地下水資源的影響。南水北調壓采條件與現狀開采條件下各市地下水儲變量統計，詳見表3。

表3 南水北調壓采條件與現狀開采條件下各市地下水儲變量統計 108 m3/a

從儲變量變化上對2 個方案進行比較，可以看出，南水北調壓采實施后，淺層地下水儲變量虧損減少4.19×108m3/a，約下降32%，淺層地下水虧損的狀況得到了較好緩解。淺層地下水-30 m 水位等值線的范圍有了明顯縮小，面積約從0.62×104km2縮減到0.41×104km2，下降34%，降落漏斗程度減輕。深層地下水的儲變量也有所減少，但由于區域內對深層地下水的開采相比淺層要小很多，因而減少幅度不大?？傮w上而言，在南水北調壓采實施后由于減少了人工開采量，有效地控制了現有降落漏斗的擴大和加深，對受水區的地下水恢復起到了很大作用，同時對緩解受水區的地下水虧損狀況和平原區的地下水緊張局面、涵養和保護地下水資源起到了很好的積極作用。

3 結論

（1）對MODFLOW 的源程序進行了優化應用。使用改進后的MODFLOW 程序模塊，實現了基于面狀處理方式的各源匯項的分類處理，在保證水均衡與模擬精度要求合理的前提下，可以較好地提高地下水流數值模型的運算速率。通過模型的外部建模方式更好地實現了各個運行文件參與平臺的集成以及后期數據的分項調整功能，彌補了MODFLOW 源匯項處理過于簡單的問題。

（2）構建了水文地質概念模型。通過水文地質條件分析確定了淀西平原為非均質各向同性、空間三維結構、非穩定地下水流系統的水文地質概念模型，確定了模型的范圍結構、邊界條件、水文地質參數及各源匯項。

（3）建立了長時間序列地下水流數值模型。在概念模型的基礎上建立了大清河淀西平原的地下水流數值模型，運用模擬期地下水長期觀測孔水位、各年各層地下水水位等值線對模型進行識別驗證，通過這些實測數據對區域內地下水開采量進行估算和評價。總體看來，淺層含水層擬合效果良好，較好地反映了實際水文地質條件。

（4）建立了不同情境下的預測模型，預測分析了地下水流場變化趨勢。從總體的預測結果來看，在南水北調壓采方案實施后地下水環境得到了較好的改善，有效地控制了現有降落漏斗的擴大和加深，對緩解受水區地下水的虧損狀況和平原區的地下水緊張局面、涵養和保護地下水資源起到了很好的積極作用，對受水區的地下水恢復起到了很大作用。