河道硬化對傍河地下水源補給結構及范圍的影響

劉 昶， 鄭明霞， 孫源媛， 蘇 婧， 傅雪梅， 丁鴻羽， 席北斗*

1.南昌大學資源環(huán)境與化工學院, 鄱陽湖環(huán)境與資源利用教育部重點實驗室， 江西 南昌 330031

2.中國環(huán)境科學研究院, 國家環(huán)境保護地下水污染模擬與控制重點實驗室， 北京 100012

傍河地下水源因水量充足、水質較好等優(yōu)點成為我國常見的水源類型. 然而，國內河道整治經常采用護岸與河床固化的方式，形成漿砌三面光河道來滿足防洪需求或達到減少河道滲透損失量與河水污染物下滲的目的. 但這一治理方式往往缺乏對河流生態(tài)環(huán)境的綜合考慮[1-2]，河道硬化阻礙了地表水與地下水的相互聯(lián)系，削減了河水對地下水的補給，降低了傍河地下水的自凈能力[3]，使傍河地下水的水量水質受到較大影響，危及地下水安全與可持續(xù)發(fā)展[4-6].

近年來，河流與地下水之間的交互作用機理和機制逐漸成為研究熱點. 馮創(chuàng)業(yè)等[7-8]在河流流域內開展了基礎性大型入滲試驗，采用水量平衡原理計算出河流對地下水的入滲量與河道的入滲能力. 任長江等[9-10]通過流域內不同植被覆蓋條件下的地表入滲試驗，擬合不同入滲模型，分析總結了不同植被覆蓋條件下的土壤入滲規(guī)律及地表水對地下水的補給. Marie-Jose等[11]將熱紅外遙感與直接超孔隙度測量相結合，精確捕捉了一條40 km河段上河水與傍河地下水的交互作用，進一步揭示了地表水與地下水交互作用在水環(huán)境系統(tǒng)上的重要性. 研究地表入滲的方法越來越豐富且先進，數(shù)值模擬法成為目前最為精確且可靠的研究方法之一. Singh等[12-17]通過建立地下水模型可直觀而精確地揭示地表水入滲地下水的狀況. Frei等[18]通用耦合河流-含水層的模型研究其交互動力學和交互模式，模型中分析并量化了地表河流對地下水的補給關系，揭示了河水與地下水的交互規(guī)律. XIE等[19]在空間上將河水與地下水的滲流作用分為垂向滲流與橫向滲流，利用熱量作為定量交換示蹤劑，結合模型研究河流-含水層的滲流作用. 雖然地表水與地下水之間的交互對于自然界水循環(huán)具有非常重要的作用，但由于缺乏對河道硬化造成的地表水入滲量及其對地下水影響的系統(tǒng)、量化的研究，河道硬化在我國河道整治工作中仍然是常見的工程措施，該措施不僅會減少河水對地下水的入滲量，還會使得地下水位下降，改變地下水的補給范圍與補給結構，對地下水安全產生威脅[20-24].

以河北省張家口市Y傍河地下水源地為研究對象，建立地下水數(shù)值模型，模擬研究區(qū)地下水流場，對研究區(qū)地下水進行水均衡分析，并運用質點追蹤技術對水源地井群的補給區(qū)進行刻畫. 通過對比河道硬化前后模型計算結果，分析并量化了河道硬化對Y傍河水源地的水均衡狀況、補給結構及其補給范圍的影響，結合歷史與近期水質數(shù)據(jù)對比分析河道硬化對地下水質條件的影響，以期為相關地下水保護提供參考和依據(jù).

1 材料與方法

1.1 研究區(qū)概況

1.1.1自然地理

Y傍河地下水源地位于張家口盆地，其地理位置與水源井分布如圖1所示. 張家口盆地為典型的斷裂型山間河谷盆地，東西長約90 km，南北寬約21 km，洋河自西向東貫穿其中，研究區(qū)是盆地中心的洋河帶狀平原，周圍低山環(huán)境，整個盆地的地勢由西向東傾斜，坡降約4%，海拔在550～1 000 m范圍內. 研究區(qū)屬北寒溫帶大陸性半干旱氣候區(qū)，2004—2018年平均降水量為394 mm，最大年降水量為496 mm(2010年)，最小年降水量為323 mm (2006年)，多年平均氣溫為5～8 ℃，最低氣溫為-26 ℃，最高氣溫為42 ℃.

圖1 Y傍河地下水源地地理位置與水源井分布

1.1.2水文地質

研究區(qū)屬山間盆地水文地質區(qū)，主要分布山地基巖裂隙水和平原區(qū)第四系孔隙地下水. 由圖2可見，盆地腹部的河流沖積平原第四系厚度可達100～300 m，地下水埋深為3～60 m，含水層巖性以中粗砂、砂礫石為主，單井涌水量10～120 m3h，為HCO3-CaMg型水，歷次水文地質勘察得出的含水層滲透系數(shù)在1～150 md，近期通過抽水試驗得到滲透系數(shù)為7.36～163.87 md，孔隙度為0.25～0.32. Y傍河地下水源地處于盆地第四系厚度最厚區(qū)域，由3個含水層組和3個不連續(xù)的隔水層組成，下部略具承壓性，其開采層位主要接受來自上游沖洪積扇地下水的側向流入補給，其排泄方式以開采和向下游徑流為主，地下水流方向大體為自西向東.

圖2 研究區(qū)水文地質剖面

1.1.3水源地基本情況

Y傍河地下水源地為傍洋河沖積平原水源地，位于張家口盆地洋河帶狀沖積平原東段的西端，地勢開闊平坦. 水源井群共有生產井18眼，沿洋河北岸呈東西向分布兩列，距離河岸100 m左右(見圖1). 該水源地生產井深96～120 m，含水層由第四系全新統(tǒng)沖積砂礫石孔隙含水層和第四系全新統(tǒng)沖洪積砂礫卵石孔隙含水層構成，補給條件良好，主要接受洋河沖積平原、洋河河水的地表水系補給與萬全區(qū)城東河和城西河沖洪積扇的地下水補給. 根據(jù)建設綜合勘察研究設計院有限公司編寫的《河北省張家口盆地地下水詳查報告》，其最大供水能力為150×103m3d，根據(jù)水廠供水量記錄計算得到實際平均取水量為54×103m3d，用水高峰期取水量為90×103m3d. 2010年張家口市政府開展洋河綜合治理工程，為增強洋河防洪能力，降低河水對地下水的污染，對河道澆筑混凝土進行防滲處理，河道形成了“三面光”狀態(tài).

1.2 研究方法

該研究中主要研究方法為數(shù)值模擬法，選用美國鹽湖城大學環(huán)境研究所開發(fā)的地下水模擬軟件GMS (Groundwater Modeling System)[25-26]來建立研究區(qū)地下水數(shù)值模型，通過對模型參數(shù)和源匯項進行識別校正，使得數(shù)學模型能如實反映研究區(qū)水文地質條件，模擬的地下水流場能有效擬合實際水位情況，從而揭示研究區(qū)地下水的水均衡狀況. 在此基礎上，使用軟件模塊MODPATH，對水源井群進行反向質點追蹤，通過示蹤粒子在地下水中向前 1 000 d的運移軌跡[27]，得到水源井群抽水的捕獲區(qū)域，即地下水源的補給范圍.

地下水流動的數(shù)學模型是基于質量守恒定律而建立的，水流控制偏微分方程：

(1)

式中：h為地下水水頭，m；Kxx、Kyy和Kzz分別為含水層在x、y、z方向上的滲透系數(shù)，md；W為單位體積流量，代表流進匯或來自源的水量，d-1；Ss為含水層孔隙介質的貯水率，m-1；t為時間，d.

1.3 概念模型

1.3.1邊界概化

確定模型范圍和劃分邊界時盡可能以自然邊界或地貌單元分界線作為模型邊界. 該模型南部和北部邊界為山區(qū)和平原區(qū)的交界地帶，接受山前地下水側向徑流補給，定義為定流量邊界；東部和西部邊界分別為盆地地下水流入和流出邊界，均定義為定水頭邊界；地表河流洋河定義為河流邊界. 模型范圍及邊界條件如圖3所示.

圖3 模型邊界條件示意

根據(jù)研究區(qū)實地勘探結果并結合張家口盆地水文地質資料對模型邊界進行確定. 該模型定義西部邊界水頭值為750 m，東部邊界水頭值為590 m；根據(jù)地下水位數(shù)據(jù)，應用達西定律計算得到，南北邊界側向流量數(shù)值，北部定流量邊界劃分為三部分，取值分別為5.00×104、10.00×104、15.00×104m3d，南部邊界劃分為兩部分，取值分別為0.45×104、0.30×104m3d.

1.3.2模型結構概化

研究區(qū)地下水以水平運動為主，流速緩慢，滲流符合達西定律，其介質是非均質的. 在模擬過程中，巖性相似的相鄰含水層通常合并為一個含水層[28]. 據(jù)勘探資料及水文地質剖面圖，模型最終概化為2個含水層組的非均質各向異性三維穩(wěn)定流模型. 根據(jù)場地高程DEM圖和鉆孔資料獲取各層頂?shù)装鍢烁? 對研究區(qū)滲透系數(shù)進行分區(qū)，如圖4所示，其中a區(qū)、b區(qū)、c區(qū)、d區(qū)、e區(qū)、f區(qū)、g區(qū)的滲透系數(shù)分別賦值為60、6、55、40、120、140、8 md. 根據(jù)巖層屬性，垂向滲透系數(shù)設為水平滲透系數(shù)的14，給水度、儲水率、有效孔隙度等參數(shù)值取對應地層屬性的經驗值[29-30].

注： a、b、c、d、e、f、g為對應的滲透系數(shù)分區(qū)編號.

1.3.3源匯項

降雨入滲量和灌溉入滲量結合蒸發(fā)量通過面狀補排程序包RCH處理計算. 研究區(qū)極限蒸發(fā)深度為4 m，根據(jù)2019年地下水水位統(tǒng)測結果，計算出地下水位埋深小于4 m的面積為9.49 km2. 模型將降雨補給與入滲補給系數(shù)整合，依據(jù)張家口市氣象資料與研究區(qū)水文地質條件，對研究區(qū)蒸散發(fā)及補給參數(shù)進行分區(qū)，如圖5所示，其中A區(qū)、B區(qū)、C區(qū)、D區(qū)、E區(qū)、F區(qū)的蒸散發(fā)及補給參數(shù)分別賦值為-0.005 0、0.000 2、0.001 0、0.002 0、0.003 0、0.000 2 md.

注： A、B、C、D、E、F為對應的蒸散發(fā)及補給系數(shù)分區(qū)編號.

河流入滲量通過河流程序包RIV處理計算. 考慮到洋河河床河堤進行過硬化處理，而河流參數(shù)主要為河床的水力傳導系數(shù)，由底積層的巖性、厚度及河床寬度、計算單元的河段長度決定[31]，計算公式：

C=KLWM

(2)

式中：K為河床土壤水力傳導系數(shù)，md；L為河流映射長度，代表河流在模型中運算的有效長度，m；W為河流寬度，m；M為河床厚度，m.

根據(jù)現(xiàn)場地質試驗，將洋河分為4個河段賦值，各河段分別賦值為0.5、0.2、0.3、0.1 m2(d·m)，該賦值范圍遠小于天然條件下的河流水力傳導系數(shù)，以期刻畫河道硬化的實際情況.

水源井開采量通過井流程序包WEL處理計算. 根據(jù)張家口市市政水務2014—2018年Y水廠供水量記錄，計算得到18口水源井實際的平均單井取水量為3.0×103td，用水高峰期的平均單井取水量為4.7×103td. 考慮后續(xù)城市人口增長造成用水量增加，模型內每口井抽水量取值為5.0×103td.

1.4 數(shù)學模型求解與參數(shù)靈敏度分析

該研究利用GMS里的MODFLOW程序建立了地下水流動模型，采用有限差分的離散方法剖分網格，模擬范圍東西方向最長為49.4 km，南北方向最長為25.1 km. 研究區(qū)剖分為132行、260列、2層，共 32 014 個活動單元，網格大小為190 m×190 m. 模擬期為2019年2月. 解算器選擇穩(wěn)定性高、收斂性好的共軛梯度解算器PCG2. 完成地下水數(shù)值模擬的模型構建和數(shù)據(jù)輸入后，將概念模型映射到MODFLOW數(shù)值模型中. 檢查無誤后保存模型并運行.

通過模型參數(shù)的靈敏度分析明確哪些參數(shù)對模型輸出結果和動態(tài)影響較大，進而減小因模型參數(shù)的不確定性產生的誤差. 對模型中含水層滲透系數(shù)、降雨入滲補給系數(shù)、河床水力傳導系數(shù)、定流量邊界水力傳導系數(shù)、河水水位和含水層孔隙度6個主要參數(shù)進行靈敏度分析，模擬水源地捕獲區(qū)最遠距離的變幅作為靈敏度分析的依據(jù). 選取各參數(shù)分別調整±5%和±10% 對數(shù)學模型進行靈敏度分析. 由圖6可見，含水層滲透系數(shù)與河床水力傳導系數(shù)是對捕獲區(qū)最遠距離影響的主要敏感因子，其他4個參數(shù)影響則較弱，表明Y傍河地下水源地的補給范圍主要受到含水層屬性與河流入滲的影響. 由于河床水力傳導系數(shù)是模型研究河道硬化的調整因素，因此含水層滲透系數(shù)的獲取與賦值大幅影響了模型結果的精度. 通過增加研究區(qū)野外試驗，獲取更密的含水層滲透系數(shù)分布，再通過克里金插值法對含水層滲透系數(shù)進行更細致合理的賦值，以減小因參數(shù)不確定性對模型運行結果可靠性產生的影響.

圖6 靈敏度分析結果

2 結果與討論

2.1 模型識別與驗證

水源地周邊共有25個水位觀測井，其水位數(shù)據(jù)可較為清晰地描繪水源地的地下水水位狀況. 先將水位數(shù)據(jù)繪制等水位線，將模擬流場等水位線與實測等水位線比對，結果如圖7所示.

圖7 模擬流場與實測等水位線對比

將25個水位觀測數(shù)據(jù)導入模型進行校準. 根據(jù)校核目標的填充顏色查看模型校正情況，當模型校正目標顏色條變小且為綠色時，可認為模型各項參數(shù)調至最優(yōu). 由圖8可見：在擬合效果圖中校核目標的條形有17個點為綠色，8個點為黃色，未出現(xiàn)紅色點位，說明觀測值與計算值的差值在校核置信范圍區(qū)間內；計算水位擬合誤差平均值為0.48 m，模型擬合效果達到要求，水流模擬結果可靠，能較好地反映地下水情況.

注： 圖(a)中綠色或黃色點位為觀測井，不同顏色表示水位擬合校核結果，其左側誤差條位于點上方時表示計算水位高于實測水位，反之則表示計算水位低于實測水位. 當觀測值與計算值差值在校核置信范圍(0～100%)時，用綠色充填；當誤差超出置信區(qū)間范圍(100%)但小于200%時，則用黃色充填；當誤差大于200%時，則用紅色充填.

2.2 模型運行結果

通過對模型進行合理調參，參數(shù)修正并不斷迭代[32]，最終使擬合效果達到要求，水流模擬結果可靠，得到能較好地反映地下水情況的模型. 由圖9可見，運行結果可以觀察到研究區(qū)地下水流場情況，其主要受地形與地層構造控制，地下水流向主要為由西向東，流向在盆地東側轉向東南方向，地下水水頭在580～780 m之間.

注： X、Y、Z分別表示空間坐標軸正東向、正北向和垂直水平向上方向.

2.3 河道硬化對傍河水源地補給區(qū)范圍影響分析

為了研究水源井群取水的補給范圍，應用MODPATH模塊模擬計算出水源井群抽水的補給區(qū). 基于已建好的MODFLOW水流模型進行質點追蹤.

對于無硬化狀態(tài)下的模擬，通過調整河流參數(shù)來模擬分析，主要調整河流傳導系數(shù)(河床土壤飽和導水率)與河水水位. 在無硬化處理狀態(tài)下，河水水位取值參考硬化處理前的歷史資料，河流傳導系數(shù)取值參考實地測量參數(shù)結果并代入式(2)計算獲得，使河流傳導系數(shù)取值更接近未硬化處理的自然河流，將河流傳導系數(shù)設定為20～30 m2(d·m).

研究區(qū)含水層介質主要為砂卵礫石，孔隙度取值為0.3. 該模型對每口井設置示蹤粒子16個，分別對硬化和無硬化的情況進行 1 000 d時長捕獲模擬，結果如圖10所示.

圖10 硬化與無硬化狀態(tài)下數(shù)值模擬水源地1 000 d補給區(qū)對比

由圖10可見：無硬化狀態(tài)下 1 000 d補給區(qū)上游部分偏向河流，為沿河流分布的狹長型區(qū)域，最長距離為 5 062 m，最寬距離為 1 002 m，面積為4.74 km2；硬化狀態(tài)下 1 000 d補給區(qū)較無硬化狀態(tài)下整體向遠離河流方向偏離，其寬度更寬，呈上游窄下游寬的“喇叭形”，沿河流方向上最長距離減少了271 m，垂直河流方向的最寬距離增加了210 m，面積增加了0.77 km2. 根據(jù)不同狀態(tài)下等水位線對比發(fā)現(xiàn)，硬化狀態(tài)時，補給區(qū)地下水水位下降了2～6 m.

河道經硬化處理后，水源井群補給范圍發(fā)生了改變，河水與地下水水力聯(lián)系變弱造成了其補給范圍向遠離河流的方向偏移，來自北部上游地下水的補給比例增大，其補給范圍的面積也進一步擴大，給水源地的用水安全帶來新的潛在風險. 觀察等水位線變化可知，硬化處理后河水對河流附近地下水補給作用減弱；河流北岸的等水位線方向也受到改變，北岸地下水流向向東南方向偏移；水源井群上游等水位線變密，下游等水位線變疏，說明井群上游地下水流速變快，下游地下水流速變慢，在地下水水位下降的情況下，水源地取水對當?shù)氐叵滤鲌龅臄_動影響變得更為強烈.

2.4 河道硬化對傍河水源地水均衡影響分析

河道硬化條件下，對井群 1 000 d補給區(qū)域進行水均衡分析. 由表1可見：井群 1 000 d補給區(qū)地下水系統(tǒng)的總補給量為251.00×103m3d，其中，上游邊界流入補給量為主要補給來源，占總補給量的82.88%；降雨入滲補給量與河流入滲補給量則分別占總補給量的14.31%和2.81%. 井群 1 000 d補給區(qū)地下水系統(tǒng)的總排泄量為251.00×103m3d，下游邊界流出排泄量與地下水開采量為主要排泄項，共計占總排泄量的98.61%；蒸散發(fā)排泄量占總排泄量的1.39%；河流排泄量為0. 井群 1 000 d補給區(qū)地下水系統(tǒng)深層對淺層含水層的越流水量為3.48×103m3d，淺層對深層含水層越流水量為3.90×103m3d.

水均衡分析(見表1)表明，水源地地下水的主要補給來源于上游地下水的流入，而河流入滲補給量與降雨入滲補給較少；地下水的去向主要為下游邊界流出和地下水開采，地下水未對河流進行排泄. 相比于無硬化處理的情況，硬化處理狀態(tài)下的河水對地下水的河流入滲流量由53.85×103m3d減至7.05×103m3d，減幅為86.91%；地下水對河水的河流排泄量由13.31×103m3d減至0，減幅為100.00%，地下水停止向河水排泄. 此外，上游邊界流入補給量減少了2.07×103m3d，減幅為0.99%；下游流出排泄量減少了35.56×103m3d，減幅為18.95%.

表1 不同條件下井群1 000 d補給區(qū)地下水均衡表

研究顯示，河道硬化不僅直接導致河流入滲量大幅減少，還間接使地下水的上游邊界流入量與下游邊界流出量隨之減少；地下水位的下降使地下水靜水位低于河床高程，導致地下水不再對河流進行排泄.

2.5 河道硬化對傍河水源地水質影響分析

研究區(qū)地下水水化學類型以HCO3-Ca型為主，河水pH、硬度及氨氮濃度均明顯高于地下水. 地下水水質例行監(jiān)測數(shù)據(jù)表明，河道硬化前研究區(qū)傍河地下水pH為7.8～8.1，總硬度為465～588 mgL(以CaCO3計)，氨氮濃度為0.02～0.30 mgL. 實地采樣測試發(fā)現(xiàn)，河道硬化后近年來地下水pH為6.73～7.82，屬中性水；礦化度為300～790 mgL，屬低、中礦化度水；總硬度為295～325 mgL(以CaCO3計)，屬軟水或稍硬水；氨氮濃度在河道硬化后均未檢出，這可能由于研究區(qū)地下水在河道硬化后切斷了河水向地下水的入滲途徑，使得河水中的污染物質，如氨氮等不能進入地下水，說明河道硬化對于防止河水對地下水的污染起到了一定的阻斷作用，使地下水水質與理化環(huán)境發(fā)生部分改變；另一方面，河道硬化后地下水位不斷下降，溶解氧與氧化還原電位監(jiān)測數(shù)據(jù)均出現(xiàn)下降趨勢，表明研究區(qū)地下水化學環(huán)境向還原性方向演化.

因此，河道硬化雖然能有效減少地表水污染物的入滲，但同時也使地下水位下降，地下水環(huán)境發(fā)生改變，進而可能導致地下水自凈能力降低[33-34]，地下水脆弱性隨之增加，導致地下水環(huán)境生態(tài)系統(tǒng)的抵抗力穩(wěn)定性變弱[2,24]，對地下水源的用水安全帶來諸多不利因素與潛在風險.

3 結論

a) 通過數(shù)值模型法模擬分析河道硬化條件下的傍河地下水源，結果表明，水源井群 1 000 d的補給范圍呈上游窄下游寬的“喇叭形”，其最長距離 4 791 m，最寬距離 1 212 m，面積為5.51 km2；對補給范圍進行水均衡計算，穩(wěn)態(tài)下地下水總補給量與總排泄量均為251.00×103m3d (不包括層間越流量)，主要的補給項為邊界流入與降雨入滲的補給，分別為208.04×103與35.91×103m3d，占比分別為82.88%與14.31%；主要的排泄項為邊界流出與地下水開采的排泄，分別為152.12×103與95.40×103m3d，占比分別為60.60%與38.01%.

b)通過模型模擬河床無硬化的自然狀態(tài)，分析河道硬化造成的影響. 對比發(fā)現(xiàn)，河道硬化使得其地下水總補給量與總排泄量減少了48.87×103m3d，河水對地下水的入滲流量減少了46.79×103m3d，入滲量減幅達86.91%，并且地下水停止了對河水的排泄；河道硬化同時使得水源井群 1 000 d補給范圍在沿河流方向上的最長距離減少了271 m，垂直河流方向上的最寬距離增加了210 m，面積增加了0.77 km2，補給區(qū)向遠離河流方向上偏移，相比自然河道條件下，補給區(qū)地下水水位下降了2～6 m.

c) 河道的硬化處理會改變河水與地下水的交互關系，使得傍河地下水的水力結構與補排條件發(fā)生改變，其傍河地下水源井群的補給范圍發(fā)生偏移，地下水水位下降，井群上游地下水流速變快，下游地下水流速變慢，水源地取水對當?shù)氐叵滤鲌龅挠绊懽兊酶鼮閺娏?；同時其地下水水質與理化環(huán)境均發(fā)生改變，傍河地下水自凈能力降低，地下水脆弱性隨之增加，給傍河地下水帶來新的潛在風險.