李 巍，郭抒燕，薛麗娟，丁曉夢(mèng)

（1.水利部海河水利委員會(huì)，天津 300170；2.天津市龍網(wǎng)科技發(fā)展有限公司，天津 300170）

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

1 研究區(qū)概況

1.1 水文地質(zhì)

大清河淀西平原區(qū)主要由山前沖洪積傾斜平原、沖洪積扇和中部沖積平原構(gòu)成，山前沖洪積平原、沖洪積扇呈扇狀交錯(cuò)分布于山前，含水砂層主要由砂礫石、粗砂、中砂、中細(xì)砂等各類砂、砂礫石組成，含水層下部無(wú)連續(xù)隔水層，垂向水力聯(lián)系良好，常因下部含水層的開(kāi)采而疏干，因此在山前常將這一含水層組與下部承壓含水層組看作是統(tǒng)一含水層組。

平原區(qū)地下水系統(tǒng)的上游為大清河沖洪積扇子系統(tǒng)，分為4個(gè)亞級(jí)，包括拒馬河沖洪積扇系統(tǒng)、瀑河-漕河沖洪積扇系統(tǒng)、唐河-界河沖洪積扇系統(tǒng)、大沙河-磁河沖洪積扇系統(tǒng)；下游為大清河古河道帶子系統(tǒng)。

1.2 含水組劃分

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

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

（2）第二含水層組。底界埋深一般120～210 m，在山前平原，與第一含水層組之間缺乏穩(wěn)定的隔水層，二者之間具有較好的水力聯(lián)系。自西向東發(fā)育2～3 套中細(xì)砂～中粗砂～砂礫石韻律層，含水層透水性與導(dǎo)水性均比第一含水層組強(qiáng)。目前，由于2個(gè)含水層組混合開(kāi)采，人工加強(qiáng)了二者之間的水力聯(lián)系。

（3）第三含水層組。底界埋深250～310 m，幾乎全部為淡水。山前平原含水層呈扇狀、扇裙?fàn)钫共迹?～4套中細(xì)砂～中粗砂～礫石韻律層組成，下段含水層更受到不同程度的風(fēng)化。典型平原區(qū)南部扇體內(nèi)單井涌水量20～50 m3（/h·m），扇間地帶單井涌水量10～20 m3（/h·m）。在大型扇體內(nèi)部與上覆第二含水層組之間無(wú)連續(xù)分布的隔水層，二者水力聯(lián)系良好，在其他地段一般都有5～10 m 的黏土或砂質(zhì)黏土分布，水力聯(lián)系變?nèi)酢?/p>

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

1.3 地下水開(kāi)發(fā)利用

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

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

2 基于優(yōu)化MODFLOW的淀西平原地下水?dāng)?shù)值模型

以下通過(guò)建立大清河淀西平原地下水?dāng)?shù)值模型，并結(jié)合南水北調(diào)受水區(qū)調(diào)水方案，預(yù)測(cè)在不同地下水開(kāi)發(fā)利用條件下的地下水位、水量變化，為指導(dǎo)地下水開(kāi)發(fā)利用、地下水管理提供科學(xué)支撐。

2.1 MODFLOW 的優(yōu)化

針對(duì)大區(qū)域模型中源匯項(xiàng)數(shù)據(jù)處理的難題，韓忠等人在MODFLOW2005 源程序的基礎(chǔ)上進(jìn)行了優(yōu)化，開(kāi)發(fā)了MODFLOW-RAW，并利用模型降雨、蒸發(fā)和開(kāi)采數(shù)據(jù)處理的新方法，簡(jiǎn)化模型的數(shù)據(jù)處理過(guò)程，縮短模型運(yùn)行時(shí)間[8，9]。其中，包括降雨子程序包（RCH）、蒸發(fā)子程序包（EVT）的處理方法改進(jìn)和面狀源匯項(xiàng)程序包（RAW）的開(kāi)發(fā)，MODFLOW-RAW 系列程序包不僅通過(guò)改變數(shù)據(jù)存儲(chǔ)方式減少了人工處理數(shù)據(jù)的工作量、加快了模型數(shù)據(jù)的讀取速度，且壓縮后的源匯項(xiàng)文件更有利于數(shù)據(jù)的修改及模型的調(diào)試。

同時(shí)，借助RAW 程序包具可同時(shí)處理最多4 個(gè)面狀源匯項(xiàng)的平行處理能力，可將區(qū)域農(nóng)業(yè)開(kāi)采、灌區(qū)開(kāi)采、灌區(qū)回灌和白洋淀入滲分別計(jì)入，使模型的前、后處理更加清晰、便捷。

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

針對(duì)淀西平原區(qū)第四系松散含水層組建模。該區(qū)域地質(zhì)條件復(fù)雜，目前的水文地質(zhì)工作基礎(chǔ)尚不足以詳細(xì)刻畫(huà)控制全區(qū)的三維物理滲流結(jié)構(gòu)，能夠取得的較為詳細(xì)的地質(zhì)調(diào)查成果為前述的4層結(jié)構(gòu)的含水層系統(tǒng)，對(duì)區(qū)域地下水系統(tǒng)結(jié)構(gòu)已經(jīng)做了相當(dāng)?shù)母呕１敬蔚叵滤當(dāng)?shù)值模型對(duì)區(qū)域地下水含水層組進(jìn)行進(jìn)一步的概化，將第一含水層組和第二含水層組合并作為淺層地下水系統(tǒng)，主要模擬潛水循環(huán)運(yùn)動(dòng)規(guī)律，同時(shí)也包含與潛水含水層緊密相關(guān)的微承壓含水層的循環(huán)在內(nèi)；將第三層和第四層合并作為統(tǒng)一的深層地下水系統(tǒng)，主要模擬區(qū)域第三層承壓水的循環(huán)運(yùn)動(dòng)規(guī)律，同時(shí)也包含循環(huán)量較小的第四層承壓水的循環(huán)在內(nèi)。

含水層厚度大、分布廣、地下水流以水平方向?yàn)橹?，垂向運(yùn)動(dòng)較弱但不可忽視，地下水運(yùn)動(dòng)空間概化為三維流；源匯項(xiàng)和水位均隨時(shí)間發(fā)生變化，表現(xiàn)為明顯的非穩(wěn)定流特征；水文地質(zhì)參數(shù)隨空間變化，體現(xiàn)非均質(zhì)性，在同一點(diǎn)的水平方向上，參數(shù)無(wú)明顯方向性，可視為各向同性。

綜上，由于第一含水層組和第二含水層組水力聯(lián)系緊密，且實(shí)際開(kāi)采多以混采井居多，故將第一含水層組和第二含水層組合并概化為第一含水層，將第三含水層組作為第二含水層，第四含水層組作為第三含水層。其中，第一含水層為潛水含水層，第二、三含水層為承壓含水層。綜合考慮地下水開(kāi)發(fā)利用層位的特點(diǎn)，參照含水層的發(fā)育程度、含水層的滲透性、地下水水力性質(zhì)、地下水動(dòng)態(tài)特征等因素，將模型在垂向上概化為三層的非均質(zhì)各項(xiàng)同性、空間三維結(jié)構(gòu)的非穩(wěn)定地下水流動(dòng)系統(tǒng)。

2.3 地下水?dāng)?shù)學(xué)模型

根據(jù)水文地質(zhì)概念模型，對(duì)非均質(zhì)各向同性、空間三維結(jié)構(gòu)的非穩(wěn)定地下水流動(dòng)系統(tǒng)，可用如下方程的定解問(wèn)題來(lái)描述：

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

2.4 地下水?dāng)?shù)值模型

2.4.1 模型離散

（1）時(shí)間離散。模擬期為2000—2010 年，為反映年內(nèi)的周期變化，以月單位作為一個(gè)應(yīng)力期，每個(gè)應(yīng)力期為一個(gè)時(shí)間步長(zhǎng)。

（2）空間離散。根據(jù)區(qū)域的含水層結(jié)構(gòu)、邊界條件和地下水流場(chǎng)特征，將模擬區(qū)第一層剖分為220行、160 列的1 km×1 km 規(guī)則網(wǎng)格，第一層有效單元共12 554個(gè)，下部第二、三層7 983個(gè)。

2.4.2 定解條件

（1）初始水位。以2000 年1 月的地下水水位資料為基礎(chǔ)，參考長(zhǎng)序列數(shù)據(jù)，繪制初始流場(chǎng)，采用內(nèi)插法和外推法獲取各含水層的初始水位。

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

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

2.4.3 源匯項(xiàng)處理

以區(qū)縣為單位，蒸發(fā)量采用改進(jìn)的MODFLOWRAW 中的EVT 蒸發(fā)程序包處理，根據(jù)蒸發(fā)極限埋深、蒸發(fā)速率及模擬水位確定。降雨入滲采用改進(jìn)的MODFLOW-RAW 中的RCH 程序包處理。灌區(qū)開(kāi)采和回灌以及其他開(kāi)采包括生活工業(yè)用水均采用面狀補(bǔ)排處理方式。開(kāi)采山前側(cè)向流量采用GMS中的WELL井流程序包處理。具體應(yīng)用中，在每個(gè)網(wǎng)格中心點(diǎn)上設(shè)置一個(gè)注水井，然后將線狀的流量數(shù)據(jù)都分配到網(wǎng)格中心點(diǎn)上，將側(cè)流量以點(diǎn)井的形式導(dǎo)入模型中。這種源匯處理方式效率高，便于調(diào)參過(guò)程中源匯數(shù)據(jù)的修改。補(bǔ)給排泄項(xiàng)對(duì)應(yīng)的程序包，詳見(jiàn)表1。

表1 模型中源匯項(xiàng)的數(shù)據(jù)對(duì)照

2.4.4 模型識(shí)別與校驗(yàn)

通過(guò)模擬水位與實(shí)測(cè)水位的對(duì)比來(lái)檢驗(yàn)長(zhǎng)序列模型的可靠性。對(duì)研究區(qū)內(nèi)19 個(gè)觀測(cè)點(diǎn)共計(jì)2 267個(gè)數(shù)據(jù)進(jìn)行了誤差分析，其平均絕對(duì)誤差為2 m。其中，11 個(gè)觀測(cè)孔的平均絕對(duì)誤差小于1.5 m，4 個(gè)觀測(cè)孔的平均絕對(duì)誤差介于1.5~2 m，4個(gè)觀測(cè)孔的平均絕對(duì)誤差大于2 m。綜上，模擬區(qū)內(nèi)誤差2 m 之內(nèi)的觀測(cè)孔約占80%，整體擬合效果較好。擬合誤差分布，詳見(jiàn)表2。

表2 擬合誤差分布

2.5 模型應(yīng)用

基于擬合驗(yàn)證后的大清河淀西平原地下水?dāng)?shù)值模型，分別建立現(xiàn)狀和南水北調(diào)實(shí)施后2 種開(kāi)采條件下的10 a預(yù)測(cè)模型，通過(guò)方案對(duì)比，定量分析地下水調(diào)控對(duì)地下水資源的影響。南水北調(diào)壓采條件與現(xiàn)狀開(kāi)采條件下各市地下水儲(chǔ)變量統(tǒng)計(jì)，詳見(jiàn)表3。

表3 南水北調(diào)壓采條件與現(xiàn)狀開(kāi)采條件下各市地下水儲(chǔ)變量統(tǒng)計(jì) 108 m3/a

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

3 結(jié)論

（1）對(duì)MODFLOW 的源程序進(jìn)行了優(yōu)化應(yīng)用。使用改進(jìn)后的MODFLOW 程序模塊，實(shí)現(xiàn)了基于面狀處理方式的各源匯項(xiàng)的分類處理，在保證水均衡與模擬精度要求合理的前提下，可以較好地提高地下水流數(shù)值模型的運(yùn)算速率。通過(guò)模型的外部建模方式更好地實(shí)現(xiàn)了各個(gè)運(yùn)行文件參與平臺(tái)的集成以及后期數(shù)據(jù)的分項(xiàng)調(diào)整功能，彌補(bǔ)了MODFLOW 源匯項(xiàng)處理過(guò)于簡(jiǎn)單的問(wèn)題。

（2）構(gòu)建了水文地質(zhì)概念模型。通過(guò)水文地質(zhì)條件分析確定了淀西平原為非均質(zhì)各向同性、空間三維結(jié)構(gòu)、非穩(wěn)定地下水流系統(tǒng)的水文地質(zhì)概念模型，確定了模型的范圍結(jié)構(gòu)、邊界條件、水文地質(zhì)參數(shù)及各源匯項(xiàng)。

（3）建立了長(zhǎng)時(shí)間序列地下水流數(shù)值模型。在概念模型的基礎(chǔ)上建立了大清河淀西平原的地下水流數(shù)值模型，運(yùn)用模擬期地下水長(zhǎng)期觀測(cè)孔水位、各年各層地下水水位等值線對(duì)模型進(jìn)行識(shí)別驗(yàn)證，通過(guò)這些實(shí)測(cè)數(shù)據(jù)對(duì)區(qū)域內(nèi)地下水開(kāi)采量進(jìn)行估算和評(píng)價(jià)?？傮w看來(lái)，淺層含水層擬合效果良好，較好地反映了實(shí)際水文地質(zhì)條件。

（4）建立了不同情境下的預(yù)測(cè)模型，預(yù)測(cè)分析了地下水流場(chǎng)變化趨勢(shì)。從總體的預(yù)測(cè)結(jié)果來(lái)看，在南水北調(diào)壓采方案實(shí)施后地下水環(huán)境得到了較好的改善，有效地控制了現(xiàn)有降落漏斗的擴(kuò)大和加深，對(duì)緩解受水區(qū)地下水的虧損狀況和平原區(qū)的地下水緊張局面、涵養(yǎng)和保護(hù)地下水資源起到了很好的積極作用，對(duì)受水區(qū)的地下水恢復(fù)起到了很大作用。