李志強,王亞娟,溫子希,張秋蘭,崔亞莉,張琳琳,楊國軍,張莉娜,劉翠珠,董鶴

(1.中國地質(zhì)大學(北京)水資源與環(huán)境學院，北京 100083;2.北京市水文總站,北京 100089)

地下水資源是水資源的重要組成部分。對于北京而言，地下水資源既是首都發(fā)展的基礎性資源，更是北京市發(fā)展的戰(zhàn)略資源和重要保障[1-4]。潮白河沖洪積扇中上游地區(qū)是北京市地下水資源最豐富的地區(qū)，在城市供水中占據(jù)了重要地位。由于多年連續(xù)超采,地下水位持續(xù)下降[5-8]。2014年12月南水北調(diào)中線工程正式通水，截至2021年3月已累計向北京供水超過62億m3，在緩解北京市供水緊張局面的同時，也為地下水的涵養(yǎng)提供了契機[9-11]，2015年通過小中河向懷河及潮白河牛欄山地區(qū)補水，通過北上臺及大水峪水庫對潮白河地區(qū)進行試驗性補水，共引水3 700多萬m3；2016—2019年均通過南水北調(diào)中線工程引水及水庫放水對潮白河流域進行補水。為貫徹落實京津冀協(xié)同發(fā)展戰(zhàn)略，深入落實北京城市總體規(guī)劃，大力推進北京市委市政府“用生態(tài)的方法解決生態(tài)的問題”、積極推進“藏水于地”的部署要求，2021年春季對潮白河實施了1個月試驗性生態(tài)補水，本文即基于此次生態(tài)補水展開相關研究。

南水北調(diào)水進京后，在滿足居民生活基本用水的情況下，按照“喝、存、補”的原則，部分來水通過水庫放水、自然河道補水的方式對水源地進行回補。補水后地下水水位、水量、水質(zhì)各方面對生態(tài)補水動態(tài)響應的研究成為關注重點。比如：永定河生態(tài)補水工程中，李海軍等[12]和湯萬龍等[13]利用監(jiān)測數(shù)據(jù)分析補水期間水位、水質(zhì)的變化，運用水均衡法分析生態(tài)補水對地下水儲存量的影響；胡立堂等[14]對永定河補水后地下水動態(tài)變化和補水滯后性進行分析，劃分了觀測孔地下水動態(tài)類型，為補水方案提供參考；霍麗濤等[15]則以2015—2016年潮白河調(diào)水回補地下水為例，分析了觀測孔水位變化和水質(zhì)變化。除了利用監(jiān)測數(shù)據(jù)分析補水對地下水的影響外，學者們還通過建立模型的手段對補水動態(tài)響應進行研究，其中：朱冰等[16]、彭勃等[17]和劉波等[18]采用Visual Modflow軟件構(gòu)建黃河三角洲地下水流數(shù)值模型，模擬分析濕地生態(tài)恢復區(qū)生態(tài)補水對地下水的影響；Sun等[19]基于人工神經(jīng)網(wǎng)絡等機器學習建立地下水模型預測了永定河補水入滲條件地下水位變化?？梢姌?gòu)建地下水流模型能較好地實現(xiàn)對地下水動態(tài)的模擬，有利于分析生態(tài)補水效果。GMS是由多種模塊構(gòu)成的一款可視化地下水數(shù)值模擬軟件，處理功能強大，適用范圍廣泛，是使用最多的地下水模擬平臺之一[20-21]。本文采用GMS軟件構(gòu)建潮白河流域地下水流數(shù)值模型，模擬分析地下水對2021年春季潮白河生態(tài)補水的動態(tài)響應，以期為科學實施生態(tài)補水工程提供參考。

1 研究區(qū)概況

研究區(qū)位于北京市東北部，包括密云、懷柔、順義3個區(qū)縣的部分平原區(qū)，研究范圍及邊界條件見圖1，總體地形西北高、東南低，地面由北向南傾斜。研究區(qū)屬溫帶大陸性半濕潤季風氣候，多年平均氣溫為11.7 ℃，日極端最高氣溫可達41.3 ℃，日極端最低氣溫可低至-23.3 ℃。年均(2000—2020年)降水量為542.58 mm，見圖2，全年降水最多的月份為6—9月，降水量占年降水總量的82%以上。

圖1 研究區(qū)范圍及邊界條件

圖2 2000年至2020年蘇莊站降雨量

研究區(qū)內(nèi)主要水系為潮白河水系。上游為潮河和白河。白河發(fā)源于河北省沽源縣，流經(jīng)赤城縣，進入北京境內(nèi)，由延慶、懷柔匯入密云水庫；潮河發(fā)源于河北省豐寧縣，經(jīng)欒平、密云注入密云水庫。潮河、白河出庫后在密云縣河槽村匯合為潮白河[22]。根據(jù)前人研究成果[23]，依據(jù)地貌單元、地層沉積環(huán)境、地下水的賦存條件以及水力聯(lián)系等特征，研究區(qū)屬于潮白河地下水子系統(tǒng)，在平原區(qū)以及山間溝谷地帶普遍分布有第四系松散沉積物，主要特征為北薄南厚，東薄西厚；顆粒從北向南由粗變細；層次也由單一到多層。根據(jù)地下水埋藏條件、水動力特征和水質(zhì)特征，以及該區(qū)水文地質(zhì)和開采利用條件，把第四系松散沉積巖類細分為單一結(jié)構(gòu)的砂卵礫石層、2～3層結(jié)構(gòu)的砂卵礫石層、多層結(jié)構(gòu)的砂礫石層夾少量砂層。

2 數(shù)據(jù)及研究方法

2.1 數(shù)據(jù)基礎

地下水模型構(gòu)建主要需要鉆孔資料、補給項、排泄項和模型校準數(shù)據(jù)等。為構(gòu)建長時間序列地下水模型，主要數(shù)據(jù)如下。

收集研究區(qū)DEM高程數(shù)據(jù)、100個鉆孔資料和5條剖面，結(jié)合開采層位選取相對穩(wěn)定厚度較大的隔水層頂板作為分層界線來確定主要含水層的頂?shù)装鍢烁摺?/p>

含水介質(zhì)參數(shù)和水文地質(zhì)參數(shù)初值根據(jù)抽水試驗、相關水文地質(zhì)報告及前人研究成果[23]確定。

密云、懷柔、順義3區(qū)縣的氣象站2000—2020年日尺度降雨量數(shù)據(jù)。

密云、懷柔、順義3區(qū)縣2000—2020年工業(yè)、生活、農(nóng)業(yè)年開采量；懷柔應急水源地、水源八廠、潮-懷應急水源地等集中水源地開采量。

密云、懷柔、順義3區(qū)縣內(nèi)長觀孔位置及2000—2020年逐月水位埋深。補水期間(2021年4月30日—5月27日)監(jiān)測井日尺度地下水位監(jiān)測數(shù)據(jù)。

南水北調(diào)2015—2021年補水路徑及補水具體情況數(shù)據(jù)。

數(shù)據(jù)主要來源為北京市水文總站，其他未說明數(shù)據(jù)均采用前人研究數(shù)據(jù)[23]。

2.2 研究方法

2.2.1模型原理

GMS(groundwater modeling simulation)軟件是國內(nèi)外常用的地下水流模擬軟件，其包含的MODFLOW程序，是基于有限差分法刻畫地下水流運動規(guī)律的計算機程序。此程序具體求解方法是通過把研究區(qū)域進行空間和時間上的離散，建立每個網(wǎng)格單元的水均衡方程式，所有網(wǎng)格單元的方程聯(lián)立成大型的線性方程組，耦合初始條件及邊界條件后，迭代求解方程組即可求得每個網(wǎng)格單元的水頭值。GMS 10.0版本以后，嵌入了MODFLOW-USG(MODFLOW-Unstructured Grid)程序，該程序基于體積有限差分法對網(wǎng)格進行剖分[24-25]。區(qū)別于以往剖分方式，MODFLOW-USG程序?qū)卧竦亩x不再局限于行、列、層，剖分方式更加靈活多變，可以使用任何類型的網(wǎng)格，即剖分時具有任意的面和節(jié)點的單元格，這種特性可以使得尖滅等地質(zhì)特征進行更真實建模，這也促使在進行網(wǎng)格剖分時可以對重點關注區(qū)域(如河道、集中水源地)進行局部加密處理。

2.2.2模型建立

根據(jù)對研究區(qū)水文地質(zhì)等條件的分析，確定模擬范圍為北起密云水庫山前一線，南至順義區(qū)蘇莊；東西界限為順義行政區(qū)劃界限各外推約1 km2。模型總面積約1 440 km2。

研究區(qū)北、西北、東北3個方向接受山區(qū)側(cè)向補給，為自然邊界，邊界流入量根據(jù)資料確定；南部與平原區(qū)接壤，與平原區(qū)的地下水發(fā)生水量交換，為通用水頭邊界，邊界流入流出量根據(jù)模型通用水頭包計算得到。潛水面為該模型的上邊界，地下水通過此邊界得到降水入滲補給、山前側(cè)向補給、河水滲漏補給、灌溉回歸入滲補給、渠系管網(wǎng)等入滲補給、生態(tài)補水等。模型的底邊界為基巖，故將其看作不透水邊界。

此次模型網(wǎng)格加密方法選擇的四叉樹網(wǎng)格剖分，對補水河道、重點關注的水源八廠、懷柔應急水源地等地進行了加密。其中,基礎網(wǎng)格為200 m×200 m，細化網(wǎng)格100 m×100 m?？臻g上將模型分為四層:第一層為潛水含水層;其他層為承壓水含水層。模型期為2000年1月至2020年12月，以1個月為1個應力期，共252個應力期；選擇2000年1月至2011年12月為模型識別期，2012年1月至2020年12月作為模型驗證期，以2000年1月流場作為初始流場，建立長時間序列地下水流模型。

模型中面狀補給項包括降雨入滲補給和灌溉回歸滲漏補給，面狀排泄項包括蒸發(fā)和部分人工開采(如工業(yè)、生活、農(nóng)業(yè)開采)。降雨入滲補給根據(jù)每個降雨入滲分區(qū)把降雨量處理成為降雨強度，降雨入滲分區(qū)依據(jù)降雨入滲系數(shù)不同劃分，利用RCH(recharge)程序包處理輸入模型。蒸發(fā)項則將蒸發(fā)量處理為蒸發(fā)強度，用EVT(evapotranspiration)程序包處理輸入模型。滲漏補給和部分人工開采則以注水井或抽水井的形式將數(shù)據(jù)處理成相應格式用井(Well)程序包輸入模型。線狀的補排項主要有山前側(cè)向補給量、南水北調(diào)補給水量和部分通用水頭邊界流入流出量，如山前側(cè)向補給量主要應用達西(Darcy)定律計算得到每個網(wǎng)格的補給量，應用ArcGIS將數(shù)據(jù)拓撲到相應位置網(wǎng)格中心點上，并計算轉(zhuǎn)換成單個網(wǎng)格上的井流量，利用Well模塊處理。通用水頭邊界的流入流出量根據(jù)Darcy定律和流場形態(tài)，通過GHB(general head boundary)程序包計算得出后代入模型。

2.2.3識別及驗證

利用同時期的地下水流場或者長觀孔歷時曲線來擬合模擬水位，識別地下水流模型中水文地質(zhì)參數(shù)、邊界值和其他源匯項，使模型愈加符合研究區(qū)的水文地質(zhì)條件，以便更精確地定量研究模擬區(qū)的補給與排泄，預報給定水資源開發(fā)利用方案下的地下水位變化。利用區(qū)內(nèi)4個長時間序列觀測孔對地下水流模型進行參數(shù)率定，模型主要參數(shù)率定結(jié)果見表1。選取國內(nèi)外研究常用的評價因子決定系數(shù)(R2)和均方根誤差(RMSE)作為評價模型參數(shù)優(yōu)化的評價標準。R2越接近1，表明模擬值越接近實測值，模擬效果越好;RMSE則反映了模型的模擬值與真實觀測值的接近程度[26]。模型驗證結(jié)果(圖3)顯示：4個觀測孔決定系數(shù)范圍為0.93～0.97；NH1、NH2、NH3三孔的RMSE范圍為3.19～3.78，另一個孔的RMSE為6.97，受開采和邊界流入量影響，數(shù)值偏大。評價因子R2和RMSE值表明模型源匯項、參數(shù)基本合理，驗證了模型的可靠性，模型可以用于生態(tài)補水過程分析。

圖3 長時間序列模型水位擬合結(jié)果

表1 模型主要參數(shù)率定結(jié)果

3 潮白河流域生態(tài)補水的模型應用

3.1 潮白河春季生態(tài)補水概況

潮白河流域試驗性補水在2014年底開始實施，利用“密云水庫調(diào)蓄工程”，在汛期通過小中河向懷河及潮白河牛欄山地區(qū)補水，通過北上臺水庫和大水峪水庫向潮白河地區(qū)補水。2018年9月，南水引水及水庫放水累計引水量達到4.56億m3；截至2019年9月，南水引水及水庫放水累計引水量超過5.7億m3。2021年4月30日起，潮白河流域?qū)嵤?2年來的首次較大規(guī)模試驗性補水，主要通過白河、潮河、雁棲河、懷河、小中河向白河干流集中補水，計劃補水總量達到2.4億m3。截至2021年5月27日，累計出庫及南水北調(diào)放水總量達2.06億m3。

3.2 潮白河春季生態(tài)補水模擬

在長時間序列地下水流模型的基礎上，將模型應用于潮白河流域2021年春季試驗性生態(tài)補水中，模型邊界條件及參數(shù)不變，預測期為2021年4月30日—2021年7月31日，1 d為1個應力期，共93個應力期；初始流場為2021年4月末流場。補水期內(nèi)(2021年4月30日—5月27日)無降水，為了充分體現(xiàn)地下水對生態(tài)補水的響應過程，補水期結(jié)束后的模型應力期內(nèi)也沒有考慮降水的影響，故模型補給項主要包括灌溉入滲、側(cè)向流入、水庫滲漏及潮白河對河道補水；研究區(qū)內(nèi)長期超采至今，導致地下水埋深大于蒸發(fā)極限埋深(4 m)，故排泄項未考慮蒸發(fā)，主要為人工開采。將數(shù)據(jù)處理成相應的格式輸入到GMS中，利用模型分析評價潮白河2021年春季生態(tài)補水效果和地下水響應。

4 結(jié)果與討論

4.1 地下水水位響應

本次補水沿潮白河流域10 km范圍內(nèi)共布置了65眼觀測井作為補水監(jiān)測井，對所搜集的監(jiān)測數(shù)據(jù)進行水位擬合，誤差基本在0～1 m的可接受范圍內(nèi)，其中4眼觀測井擬合情況見圖4。模型模擬水位和實際水位數(shù)據(jù)較為接近，模型參數(shù)和源匯項基本符合實際情況。分析擬合孔水位變化可知：7號孔反應迅速，補水后不久水位即開始上升，這與該孔位于上游且離河道較近有關；其他擬合孔水位在補水一段時間后，水位上升明顯。截取剖面A-A′，見圖5，隨著時間的推移，到7月31日距離主河道兩岸的水位開始消散，而距離河道較遠的地方水位逐漸回升。24號和27號擬合水位出現(xiàn)0～1 m波動，擬合效果一般，可能是在補水初期受到懷柔應急水源地和水源八廠開采的影響。

圖4 觀測孔水位響應及擬合結(jié)果

圖5 A-A′剖面水位變化

4.2 補水影響范圍

根據(jù)模型運行結(jié)果，考慮生態(tài)補水作為單一變量影響下地下水水位動態(tài)的變化。在補水預測模型基礎上，考慮研究區(qū)無補水情形下的地下水水位。將兩者補水期末的水位相減得到補水期間動態(tài)水位抬升圖，見圖6。分析可知，補水期末，沿著補水路徑周圍水位均有所上升，最大上升幅度接近18 m。以水位抬升值大于0劃定此次補水的影響范圍，計算得到研究區(qū)內(nèi)的影響面積約為842.8 km2，占研究區(qū)總面積的58.6%?？拷卑缀痈闪餮匕兜叵滤簧仙却笥?0 m，影響面積約為44.2 km2，占比約為3%。

圖6 補水期末水位變幅

4.3 補水對儲存量的影響

補水期間，補水對地下水的入滲是研究區(qū)的主要補給來源。根據(jù)補水期間地下水觀測資料，確定觀測孔初始水位值和補水期末水位差值，圈定各觀測孔所控制的面積，取用相應的給水度，算出整個研究區(qū)的儲存量變化值。

(1)

式中：ΔQ儲為儲存量的變化量，m3；μi為第i個觀測孔所對應的給水度；Δhi儲為第i個觀測孔的水位差值，m；Fi為第i個觀測孔所圈定的面積大小，km2。

根據(jù)式1計算可知,在總水量為2.06億m3的前提下，補水滲漏補給地下水量為1.57億m3(約占總水量的76%)，研究區(qū)第一含水層系統(tǒng)儲存量增加了1.49億m3。由此可見，生態(tài)補水對于潮白河沖洪積扇地區(qū)地下水的恢復具有重要的意義。

4.4 地下水位的滯后效應

地下水的抬升和生態(tài)補水有著密切的聯(lián)系，但是地下水水位的回升需要一定的時間，即存在一定的滯后性[26]。根據(jù)現(xiàn)有監(jiān)測井分布特點，選取距離潮白河河道兩岸的4個觀測井分析地下水水位對補水的滯后響應，見圖7,南岸為13號(3 km)、17號(6 km)，北岸為11號(3 km)、19號(6 km)。對比同側(cè)觀測井可以發(fā)現(xiàn)：13號、11號觀測井離河道較近，對補水的響應及時，水位抬升迅速，響應滯后時間較快(t0=3 d,t2=6 d);17號、9號距離河道較遠，地下水響應時間較慢(t1=30 d,t3=18 d)，水位上升趨勢較平緩。13號觀測井響應滯后時間最短，水位上升最多，分析其原因是從補水初期開始最先受到補水入滲持續(xù)的影響，17號觀測井響應滯后時間最長，水位上升緩慢，其原因是距河道較遠且僅受到潮白河干流補水的影響。

圖7 觀測井水位滯后響應

5 結(jié) 論

補水對研究區(qū)范圍內(nèi)水位的抬升有明顯效果，考慮有補水和無補水情形，補水后最大抬升值接近18 m。以水位抬升0為分界線，劃分此次補水影響范圍約為842.8 km2，占研究區(qū)總面積的58.6%，補水效果較好。

截至2021年5月27日，潮白河春季試驗性生態(tài)補水總量達2.06億m3，補水下滲補給地下水量為1.57億m3，約占補水總量的76%。補水1個月后，第一含水層系統(tǒng)儲存量增加了1.49億m3。生態(tài)補水對于潮白河沖洪積扇地區(qū)地下水的恢復具有重要的意義。

地下水水位對補水的響應存在一定的滯后性。離河道距離越近，響應越迅速，反之越慢。

致謝：

感謝北京市水文總站對此次研究提供的數(shù)據(jù)資料。