潮白河春季生態(tài)補(bǔ)水及地下水響應(yīng)

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

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

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

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

1 研究區(qū)概況

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 研究方法

2.2.1模型原理

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

2.2.2模型建立

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

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

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

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

2.2.3識(shí)別及驗(yàn)證

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

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

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

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

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

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

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

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

4 結(jié)果與討論

4.1 地下水水位響應(yīng)

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

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

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

4.2 補(bǔ)水影響范圍

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

圖6 補(bǔ)水期末水位變幅

4.3 補(bǔ)水對儲(chǔ)存量的影響

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

(1)

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

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

4.4 地下水位的滯后效應(yīng)

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

圖7 觀測井水位滯后響應(yīng)

5 結(jié) 論

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

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

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

致謝：

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