截滲墻滲漏對黃水河下游地區地下水水質的影響研究?

劉貫群，張淑琪，李 熊

(1.中國海洋大學海洋環境與生態教育部重點實驗室，山東 青島 266100；2.中國海洋大學環境科學與工程學院，山東 青島 266100；3.中國海洋大學山東省海洋環境地質工程重點實驗室，山東 青島 266100)

地下水資源是中國許多城市生活、生產、生態用水的重要水源，對當地經濟社會發展起著重要作用。濱海地區是人口高度集中、經濟快速發展的地區，因對地下淡水的需求量過高導致過度開采，破壞了咸水與內陸淡水的動態平衡，最終出現海水侵入淡水含水層的現象[1]。目前，全球已有幾十個國家和地區出現海水入侵現象，導致地下水水質惡化，引發了一系列生態環境問題，嚴重制約經濟社會的可持續發展[2-4]。黃水河流域是龍口市工農業生產和居民生活的重要供水水源地，是該市經濟發展的基礎保障。20世紀70年代后期黃水河流域出現海水入侵現象，到20世紀80年代海水入侵現象加重[5]，為了控制黃水河區域海水倒灌，1995年在黃水河下游距海岸線 1.2 km處修建一座地下水庫截滲墻，截滲墻具有攔截入海地下潛流和阻擋海水入侵的雙重作用。截滲墻的修建有效阻止了截滲墻兩側的水力聯系，同時也使內陸地區殘留了一定區域的咸水體。

了解濱海含水層海水入侵的影響因素對于地下水資源的合理開發利用具有重要意義，多位學者研究發現海平面上升、氣候變化、地下水開采對海水入侵程度產生影響[6-9]。針對截滲墻這一影響因素，國內外多位學者研究了截滲墻類型、高度、位置等因素變化對海水入侵的影響[10-12]，少有學者考慮截滲墻墻體出現滲漏的可能性，而水利工程施工工序較多，施工人員、外部環境及材料等因素將導致在施工過程中或建成后出現滲漏的情況[13]，海水仍可侵入淡水含水層，故有必要研究其對地下水水質的影響程度。本文運用SEAWAT模塊對黃水河下游地區建立三維變密度地下水數值模擬模型，針對不同截滲墻滲漏情況設計多種預測方案，分析未來研究區地下水氯離子濃度的演變趨勢。

1 研究區概況

龍口市是山東省沿海地區重要開放城市，位于膠東半島的西北部，渤海灣南畔。龍口市南部地勢較高，北部地勢平坦，其地表水主要集中在南部山區，地下水主要分布在北部平原區。黃水河為龍口市境內最大的河流，發源于棲霞縣北部山區，從東南向西北流經龍口市境內流入渤海，干流全長55.43 km，流域面積 1 034.47 km2。研究區域位于黃水河下游，地處龍口市東北部的平原區，地理坐標為120°30′20″E—120°38′16″E, 37°37′16″N—37°44′36″N，總面積約85 km2，具體位置如圖1所示。

圖1 研究區地理位置示意圖

研究區屬于暖溫帶半濕潤季風型大陸性氣候[14]，四季分明，年平均氣溫11.8 ℃，年平均蒸發量1 150～1 250 mm。根據研究區內的諸由觀水文站1966—2020年降水量資料顯示，年平均降水量為557 mm，降水多集中在6—9月份，約占全年總降水量的72%，7—8月份降水量最多(見圖2)。

圖2 諸由觀水文站降水量圖

2 海水入侵數值模擬研究

2.1 水文地質概念模型

黃水河下游地區主要含水層為第四系含水層，含水層厚度在20～55 m之間。潛水含水層可大致分為兩層，上層含水層巖性以礫質粗砂為主，下層含水層主要由粗砂、中砂組成，局部為細砂，上層含水層滲透性比下層含水層強。研究區水文地質剖面如圖3所示，剖面位置參見圖1。依據含水層巖性、埋藏條件及開采條件，將研究區潛水含水層在垂直方向上劃分為兩層，概化為非均質各向同性的三維潛水含水層。

圖3 水文地質剖面圖

研究區含水層在東西兩側逐漸尖滅，東、西兩側視為隔水邊界。北部邊界視為定水頭邊界，采用渤海常年海平面平均值作為定水頭邊界的值，取值為0 m。研究區地下水位自南向北呈遞減趨勢，南部邊界受到上游地下水補給，地下水側向流入研究區域，因此將南部邊界視為給定流量邊界。對于溶質運移模型，東、西部邊界設置為零質量通量邊界。北部邊界定為定濃度邊界，濃度值采用研究區實測海水樣氯離子值，取值為19 000 mg/L。南部邊界定義為定濃度邊界條件。研究區域的上邊界為潛水含水層自由水面，含水層的下邊界為透水性差的粘土層，組成區域性隔水底板。

2.2 數學模型

變密度過渡帶模型考慮密度對水頭、流速和濃度的影響，能夠有效刻畫復雜水文地質條件、人類活動等因素影響下的咸淡水界面運移規律[15-17]，是當前海水入侵數值模擬的重要研究方向。過渡帶海水入侵的數學模型中須用兩個偏微分方程來描述, 方程(1)描述了密度不斷改變的混合液體(咸淡水混合物)的流動。

(1)

式中：hf為等效淡水水頭值(m)；Kfx、Kfy、Kfz為沿x、y、z方向的滲透系數(m/d)；Z為計算點高程(m)；t為時間(d)；ρ為實際液體密度(mg·L-1)；ρf為等效淡水的密度(mg·L-1)；Sf為等效淡水單位貯水率(1/m)；θ為有效孔隙度；C為溶解物質的濃度(mg·L-1)；ρs為源、匯項中溶解物質的密度(mg·L-1)。qs為單位時間進入單位體積含水層源、匯項的體積(d-1);Ω為計算區；h為(x,y,z)在t時刻的水頭值(m)；h0為點(x,y,z)處的初始水位；q(x,y,z)為第二類邊界上單位面積補給量，q=0代表隔水邊界,為第二類邊界條件。

方程(2)描述了混合液體中鹽分的運移。在考慮溶質運移過程，即鹽分運移過程時，選用化學性質上較為穩定的氯離子作為模擬組分。在海水入侵過程中,氯離子的運移方式以對流和彌散為主，暫不考慮運移過程中可能進行的化學反應。

(2)

2.3 變密度地下水數值模型

2.3.1 網格剖分與時間離散 利用Visual MODFLOW中的SEAWAT模塊建立研究區變密度地下水數值模擬模型，SEAWAT模塊廣泛應用于例如海水入侵這樣的變密度流問題，其模擬結果具有較高的準確性與可信度[18]。根據模型區域內地層、流場、邊界條件等特征，以及考慮到模型計算的精度和耗時長短，恰當地設計有限差分單元網格，經綜合考慮，平面上將研究區剖分為70行、65列邊長為200 m的正方形網格結構，垂向上將模型劃分為2層。

整個模擬期包括兩個階段:2019年4月18日—10月31日為模型識別期，2019年11月1日—2020年9月1日為模型驗證期。選擇研究區內24眼地下水觀測井的實地觀測資料進行模型識別與驗證，以研究區2019年 4 月 18 日的實測數據作為初始地下水流場和初始濃度場。初始地下水位等值線圖及初始氯離子濃度等值線圖如圖4所示。

圖4 初始地下水位等值線圖和氯離子等值線圖

2.3.2 水文地質參數及源匯項 水文地質參數關系到地下水模型的準確性，是進行精確科學計算地下水資源和地下水均衡的基礎。水文地質參數分區及初值根據研究區水文地質條件和前人抽水試驗資料確定[19-20]，針對研究成果較少的水文地質參數，則根據研究區不同巖性的經驗值進行確定[21]。最終將研究區劃分11個水文地質分區，具體水文地質參數分區見圖5。

圖5 含水層水文地質參數分區圖

研究區潛水含水層主要受大氣降水補給和側向補給,近幾十年研究區降水量較少，黃水河上游王屋水庫無下泄流量，導致河水基本被水庫截留，黃水河一直處于斷流狀態，不與地下水進行水量交換，故在模型建立時暫不考慮河流作用。根據研究區內黃河營、諸由觀、蘭高3個雨量站的降水量數據,采用泰森多邊形法進行大氣降水補給分區(見圖6)，再根據降水入滲系數方法計算降水入滲補給量，3個雨量站的降水量數據均來自煙臺水文信息網的實際統測數據。南部邊界以側向補給的形式給研究區提供補給量，采用達西公式計算其側向徑流補給量。

圖6 降水分區圖

研究區地下水的主要排泄方式包括城市生活、工業、農業等人工開采；研究區地下水埋深4～18 m，均大于3 m，因此不考慮蒸發排泄。根據煙臺市水資源公報和統計年鑒，確定識別期和驗證期的地下水開采總量，農業開采量按照山東水科院提供的農業開采月比例系數分配到每月，工業開采和生活開采量按月平均分配。

2.4 模型識別與驗證

采用試估-校正法對模型進行校正，反復調整水文地質參數以達到較為理想的擬合效果，識別驗證后水文地質參數值見表1。圖7(a)和(b)分別為地下水位和氯離子濃度值的擬合情況，識別與驗證期水位擬合絕對誤差小于0.5 m的分別達到82%和78%，均高于70%；識別與驗證期氯離子濃度絕對誤差小于15%所占比分別為90%和87%，均高于50%，滿足精度要求，說明模型已達到較為理想的擬合效果。圖7(c)和(d)分別為研究區從北到南三眼觀測井識別期地下水位和氯離子濃度計算值與實測值的擬合情況，圖7(e)和(f)為驗證期擬合情況，可以看出模型計算結果與實測值的趨勢一致，說明通過識別和驗證，得到了能夠比較真實反映研究區地下水流場和氯離子運移實際情況的數值模型。

圖7 模型識別期與驗證期擬合情況

表1 各區水文地質參數值

3 截滲墻滲漏對地下水的影響研究

3.1 預測方案構建

對1966—2020年黃水河流域降水量進行降水頻率分析，得出研究區實際降水量變化規律為豐-枯-平-枯-枯；豐水年(P=25%)降水量為710.7 mm，平水年(P=50%)降水量為517.5 mm，枯水年(P=75%)降水量為400 mm。以2020年9月1日為基準年，預測未來5年地下水氯離子濃度變化趨勢,根據降水量年內變化規律，將一年劃分為4個應力期(9—10月、11—4月、5—6月、7—8月)，預測期共有20個應力期。

在豐-枯-平-枯-枯水年條件下，探究地下水庫截滲墻不同滲漏位置對研究區地下水水質的影響，具體預測方案如表2所示。首先，考慮截滲墻不同方位滲漏對研究區地下水的影響，將截滲墻邊界按地理位置和水文地質條件，劃分西部、中部和東部分區，并在3個分區內各選取1處滲漏點，每一處裂縫寬度取為1 m，對滲漏位置所在的單元格進行加密處理，將邊長為200 m的單元格剖分為寬度為1 m的單元格；截滲墻用Visual MODFLOW中的水平阻滯(HFB)程序包來模擬，考慮滲漏時移除滲漏位置處的墻體。其次，選取截滲墻西部滲漏的情況，分析截滲墻位于上層含水層和下層含水層的部分分別滲漏時對研究區地下水水質的影響，滲漏點的寬度為1 m且上部和下部滲漏點的平面位置相同。

表2 海水入侵預測方案

3.2 預測結果分析

以氯離子高于250 mg/L的面積作為海水入侵總面積，表3顯示了各預測方案下研究區海水入侵程度，S0代表現狀年(2020年9月1日)海水入侵程度。從表3和圖8中可以看出，在豐-枯-平-枯-枯水年條件下，考慮地下水庫截滲墻滲漏的情況時，研究區海水入侵范圍擴大。截滲墻西部和中部出現滲漏時的海水入侵范圍增加幅度相近，與現狀年相比分別增加9.52%和9.11%。當截滲墻東部位置的墻體出現滲漏時，對研究區地下水的影響最大，氯離子濃度高于250 mg/L面積高達23.38 km2，與現狀年相比增加了3.49 km2，增長了17.58%。這是由于現狀年時研究區東北部地下水為淡水(<250 mg/L)，當海水侵入時，地下水逐漸咸化，導致研究區高于250 mg/L的面積相比現狀年增加幅度最大。

表3 不同預測方案下海水入侵程度對比

圖8 模型運行5年后截滲墻不同滲漏位置氯離子濃度變化圖

14號井靠近截滲墻的西部地區，距離西部、中部、東部滲漏點分別約190、3 000和5 800 m。圖9(a)為14號地下水井不同滲漏位置氯離子濃度變化圖，可以看出14號井氯離子濃度值受截滲墻西部滲漏的影響最大，受截滲墻中部和東部滲漏的影響較小，預計2025年9月1日的Cl-濃度比中部、東部滲漏時的Cl-濃度高1 100 mg/L左右。17號井靠近截滲墻中部地區，距離西部、中部、東部滲漏點分別約2 430、1 080和2 800 m。

圖9 不同位置單獨滲漏時各地下水井Cl-濃度變化圖

從圖9(b)中可以看出17號井受截滲墻中部滲漏的影響較大，截滲墻西部滲漏和東部滲漏對其影響相近，預計5年后中部滲漏時Cl-濃度比西部、東部滲漏時大約高80 mg/L。靠近截滲墻東部的3號井距離西部、中部、東部滲漏點分別為5 150、1 350和3 080 m，模型運行5年后，當東部出現滲漏時，氯離子濃度與西部和中部滲漏相比僅增加30 mg/L左右(見圖9(c))。因此，可以得出當截滲墻出現滲漏時，距離截滲墻越近的地下水井氯離子濃度受其滲漏的影響越大，且靠近一處滲漏位置的地下水井，受其他兩處滲漏時的影響較小。

表3和圖10反映了截滲墻上部和下部分別滲漏時對研究區地下水水質的影響，可以看出，模型運行5年后，與截滲墻下部滲漏相比，當上部出現滲漏時海水入侵總面積更大，地下水咸化程度更嚴重；上部滲漏時 14號井的氯離子濃度比下部滲漏高300 mg/L左右，這是由于上部含水層滲透性較強，當截滲墻出現滲漏，海水涌入時，向內陸入侵的速度更快，導致內陸地區地下水咸化程度更重。

圖10 截滲墻位于不同含水層位置滲漏時14號井中Cl-濃度變化圖

4 結論

基于SEAWAT模塊建立變密度地下水流模型和溶質運移模型，利用研究區域水文地質資料及地下水動態監測數據，通過反復調整水文地質參數，對模型進行了識別和驗證，得到了可以較好的反映研究區地下水流場和氯離子運移過程的地下水數值模型。為了探究地下水庫截滲墻滲漏對研究區地下水的影響，在豐-枯-平-枯-枯水年，針對截滲墻不同位置出現滲漏設計多種預測方案，預測結果如下：

(1)當截滲墻不同位置出現滲漏時，研究區地下水氯離子濃度與現狀年相比均升高。以地下水氯離子濃度高于250 mg/L的面積作為海水入侵總面積，西部滲漏和中部滲漏時，海水入侵面積增速相近。當截滲墻東部出現滲漏時，對研究區地下水的影響最大，海水入侵面積與現狀年相比增長了17.58%。

(2)當截滲墻出現滲漏時，距離截滲墻越近的地下水井受其滲漏的影響越大，且靠近一處滲漏位置的地下水井受其他兩處滲漏時的影響較小。

(3)針對截滲墻位于不同含水層部分的滲漏情況進行預測分析，得出截滲墻上部出現滲漏時地下水咸化程度更嚴重。

(4)截滲墻滲漏將會加重黃水河下游地區海水入侵程度，影響研究區地下水資源的開發利用，因此應高度重視黃水河地下水庫截滲墻滲漏的探測工作，可適當采取防滲加固措施，以確保截滲墻發揮良好的截滲效果。此外，應減少研究區地下水開采量，抬高截滲墻南側區域的地下水位，縮小截滲墻南北兩側的水力梯度，以減少截滲墻出現滲漏時所造成的影響。模擬及預測結果可為未來地下水管理提供理論依據，對經濟社會的可持續發展具有重要的現實意義。