水力屏障和截滲墻在海水入侵防治中的數值模擬研究

呂盼盼，宋 健，吳劍鋒，吳吉春

（南京大學地球科學與工程學院/表生地球化學教育部重點實驗室，江蘇 南京 210023）

海水入侵是由淡水和海水之間的密度差異驅動的咸-淡水界面向陸地方向遷移并達到平衡的自然現象[1?2]。但是，隨著社會經濟的發展以及氣候變化等外部因素，濱海地區的地下水超采、土地利用類型的改變以及海平面上升等現象破壞了原有的咸淡水界面的動態平衡，進一步加劇了海水入侵的程度[3?7]。海水入侵已在全球范圍內發展，造成了濱海地區可利用淡水資源的減少以及開采井的報廢等不利影響，極大地制約了濱海地區的社會經濟發展[8?9]。在以往研究中，依據濱海含水層中過渡帶的類型可將海水入侵研究分為突變界面模型和過渡帶模型，其中突變界面模型假設海水和淡水為非混溶流體，分界面固定，但是在實際場地條件下海水和淡水具有可混溶性，因此過渡帶模型更符合海水入侵的規律[10?11]。在人類活動引發海水入侵的濱海地區，需要耗費大量經濟成本構建復雜的工程措施防止濱海含水層中可利用淡水資源的進一步污染，并且治理的周期十分漫長。

當前，國內外學者已提出多種方法應對海水入侵問題，如控制地下水開采量、人工補給地下水、建立地下水截滲墻等[12?16]。如Ebeling等[17]利用SEAWAT和FloPy建立了一個二維變密度海水入侵模型，研究混合水力屏障方案的可行性和最優管理策略，并與單個抽注水力屏障進行對比分析。Kaleris等[18]利用數值模擬方法研究了截滲墻的深度、距海岸線的距離、截滲墻滲透率以及含水介質的各向異性對海水入侵的影響，并評估了其對海岸附近地下水開采的保護作用。Abdoulhalik等[19]提出了結合不透水截滲墻和半透水地下壩的混合物理屏障作為控制海水入侵的新屏障系統。該系統可引起淡水推動咸水朝海岸線方向向上抬升，可顯著減小濱海含水層咸水體的入侵范圍。

本文從建立水力屏障和物理屏障（截滲墻）2個角度，構建基于SEAWAT-2000程序的室內二維砂箱試驗數值模型，并從物理屏障（截滲墻）的角度構建基于山東龍口地區某典型剖面的海水入侵數值模型。結果揭示了不同水力或物理屏障布設方式對咸淡水界面運移規律的影響，可為場地條件下濱海含水層海水入侵防治中的水力屏障和截滲墻布設提出合理化建議與決策。

1 研究方法

1.1 海水入侵模型

基于SEAWAT-2000程序構建數值模型模擬二維砂箱試驗以及山東龍口地區典型剖面在不同注水井、截滲墻布設方案下的咸水入侵過程[20]。SEAWAT是由美國地質調查局開發的用于模擬變密度地下水水流和鹽分遷移的模型。水流方程與溶質運移方程為：

式中：ρ—流體密度；

Kfi—主軸方向的等效淡水滲透系數；

ρf—淡水密度；

Z—計算點高程；

hf—等效淡水水頭；

Sf—等效淡水單位儲水系數；

t—時間；

θ—有效孔隙度；

C—溶質濃度；

ρs—源（匯）流體密度；

qs—源（匯）單位體積流量；

Ck—溶質k的溶解濃度；

Dij—水動力彌散張量；

vi—滲流速度；

—源（匯）流中溶質k的濃度；

Rn—化學反應項。

該程序設計的基本原理是：以MODFLOW模擬地下水流過程，MT3DMS模擬溶質運移過程，并將兩者耦合求解變密度條件下地下水流和溶質運移的耦合模型。SEAWAT已成為目前研究變密度地下水問題最常用的模型之一，并表現出較好的穩定性與有效性[21?25]。

1.2 海水入侵管理措施

為探究注水井的井位、流量及截滲墻的位置、貫穿深度等參數對咸淡水界面運移規律的影響，本文采用SEAWAT模型，模擬咸淡水界面達到自然條件下的穩定狀態，設置注水井或截滲墻的不同工況。注水井布設于自然穩定狀態下鹽水楔的外部，為僅底部滲水非完整井的點源補給。截滲墻沿地表向地下貫穿且假設為不透水介質，模擬均為單井（或單個截滲墻）。注水井系列模擬方案的變量僅為注水井位置或補給強度，截滲墻系列模擬中變量僅為墻體水平位置或貫穿深度。在加入注水井或截滲墻后，咸淡水界面達到新的穩定狀態，此時記錄鹽水楔前鋒到達的位置，并對比初始狀態下鹽水楔前鋒位置，計算表征鹽水楔驅退程度的海水入侵回退系數（R）[26]，即為注水井和截滲墻管理措施的有效性。回退系數（R）為：

式中：L0—初始穩定條件下鹽水楔前鋒到海岸線的 距離；

L—設置注水井或截滲墻措施后再次達到穩定時鹽水楔前鋒到海岸線的距離。

2 二維砂箱數值模型

2.1 試驗模型

模擬區為二維砂箱，可概化為二維、均質、各向同性的潛水含水層（圖1）。模擬范圍長90.0 cm，寬1 cm，高41.6 cm，為典型的二維剖面。將無外界干預下咸淡水界面達到自然平衡的狀態設為初始條件。模擬區域頂部為零通量邊界；底部為隔水邊界；左側設置定水頭、定濃度邊界，邊界濃度以TDS表征設為33.6 g/L；右側設置定水頭邊界，初始濃度為0 g/L。理想砂箱實驗的概念模型見圖1。采用SEAWAT-2000程序對模型求解，利用有限差分方法將其在空間上剖分成41列和91行的正方形網格，差分網格間隔為Δx=Δz=1 cm。在天然咸淡水界面達到穩定狀態后，設定不同工程措施的管理期為0.5 h。該理想算例主要參數取自文獻[26]，見表1。

圖1 理想砂箱概念模型示意圖Fig.1 Conceptual model of the ideal sandbox

表1 理想砂箱試驗地下水數值模型主要參數[26]Table 1 Input parameters for the numerical model of the ideal sandbox

2.2 模擬結果與分析

2.2.1 注水井工程措施評價

（1）注水井流量

為探究注水井流量對咸淡水界面的影響，在初始咸水楔前鋒附近設置注水井（圖1，假設僅在注水井底部進水），并改變流量進行不同補給情景下的數值模擬。穩定后記錄鹽水楔前鋒位置L，根據式（3）計算不同補給情景下的海水入侵回退系數。海水入侵回退系數（R）與補給流量（Q）的關系見圖2。

圖2 海水入侵回退系數R與補給流量Q的關系Fig.2 Repulsion ratio with respect to recharge quantity

模擬結果表明，在相同補給點位下，海水入侵回退程度隨補給強度的增大而增大，回退系數R與補給流量Q呈非線性正相關關系。當Q大于0.005 m3/d時，隨著Q的增大R產生的增量明顯減小，且緩慢趨近于1。研究結果表明隨著注水井流量的增大，淡水體水頭不斷抬升，造成濱海含水層中地下水向海洋的排泄量增加，引起咸淡水界面向海岸線方向移動，實現了驅退海水入侵前鋒的目的。

（2）注水井位置

為探究不同注水井位置對咸淡水界面的影響，識別最大化海水入侵回退距離的最優注水井位置，在初始咸水楔外部不同位置以20 cm×10 cm的矩形網格節點形式均勻設置一系列相同補給強度（4.2×10?5m3/d）注水井，并在初始鹽水楔前鋒附近局部加密了注水井分布。圖3表示海水入侵回退系數R與注水井（僅底部進水）位置的關系圖。

如圖3所示，位于鹽水楔前鋒左上角的W3井（x=0.40 m，z= 0.05 m）的回退系數達到最大值，約為21%，表明W3井補給的淡水形成了一個從含水層下部產生的水力屏障，達到了最優的驅退海水入侵前鋒的管理效果。鹽水楔前鋒處W1井的位置達到19%的回退系數，此時產生的補給強度仍然有效，但沒有W3井的效果明顯，距離W3井右邊10 cm的W4井也能達到約14%的回退系數。總體上，注水點距W3井越遠，海水入侵回退系數越小。

圖3 二維砂箱試驗模型海水入侵回退系數R與注水井位置關系Fig.3 Repulsion ratio with respect to recharge well position

模擬結果表明：注水井注入淡水最有效的位置是靠近W3井的位置。如果將淡水補給直接施加在W3井附近，則產生的效果最佳，但當施加的淡水補給遠離W3井的位置時效果逐漸減弱；表面補給（如W15、W16、W17、W18）在控制海水入侵方面不如鹽水楔前鋒附近的注入措施更有效。對于完整井的補給井，同樣可以確定為W3、W8、W12、W16構成最優組合，最有利于驅替海水入侵的前鋒。由于注水井工程成本隨著深度的增加而增大，為節約成本，可適當縮減井的貫穿深度，將注水井布設在靠近W3處。對于表面補給，最有效的位置應在W3井的正上方。

2.2.2 截滲墻工程措施評價

為研究截滲墻對防治海水入侵的影響，模型達到自然條件下的穩定狀態后，加入不透水截滲墻并改變其水平位置和貫穿深度完成系列模擬，并計算不同情景模式下的海水入侵回退系數R。貫穿深度以無量綱的貫穿率h/H表示，截滲墻的水平位置以無量綱的1-Xb/L0表示。其中：h表示截滲墻貫穿深度，H表示含水層厚度，Xb表示截滲墻的水平位置，L0為初始狀態下鹽水楔前鋒位置。

（1）截滲墻水平位置的影響

基于不同水平位置截滲墻條件下的模擬結果，計算海水入侵回退系數R。從圖4（a）可以看出，當1-Xb/L0<0時，即截滲墻水平位置位于初始鹽水楔前鋒靠岸一側時，海水入侵回退系數R值為負數，鹽水楔進一步向淡水區推進。當1-Xb/L0=0時，即截滲墻水平位置恰位于初始鹽水楔前鋒位置時，海水入侵回退系數R值接近0。當1-Xb/L0>0時，即截滲墻水平位置位于初始鹽水楔前鋒靠海一側時，海水入侵回退系數R值為正數，鹽水楔形體向海岸線方向回退。因此，在不同貫穿率下，截滲墻水平位置越靠近海岸線，海水入侵回退系數越大，截滲墻工程效果越顯著。在不同貫穿率下回退系數R與1-Xb/L0呈線性關系。

圖4 二維砂箱試驗模型回退系數R與截滲墻水平位置1-Xb/L0和截滲墻貫穿率h/H的關系Fig.4 Repulsion ratio（R）with respect to horizontal barrier location expressed in（1-Xb/L0）and respect to penetration ratio（h/H）

（2）截滲墻貫穿深度的影響

基于不同貫穿深度截滲墻條件下的模擬結果，計算海水入侵回退系數R。從圖4（b）可以看出，沒有加入截滲墻時，鹽水楔前鋒位置沒有變化，初始值R=0。當1-Xb/L0>0時，隨著h/H增加，即截滲墻向含水層底部貫穿時，回退系數R增加；當貫穿率h/H接近1時，鹽水楔前鋒被阻擋在截滲墻位置之外，實現鹽水楔最大回退值。當1-Xb/L0<0時，即截滲墻位于鹽水楔靠岸一側時，回退系數R為負值，表明隨著截滲墻貫穿深度增加，截滲墻阻止了淡水水流對海水入侵的排斥作用，導致鹽水楔前鋒進一步向淡水體推進，加劇了海水入侵程度。因此，截滲墻貫穿深度對控制海水入侵是否有效取決于截滲墻水平位置。在不同的截滲墻水平位置下，R與h/H呈非線性關系。

（3）截滲墻水平位置與貫穿深度綜合影響

根據多組模擬結果繪制在不同水平位置及貫穿深度下回退系數（R）變化的等值線圖，見圖5。

從圖5可以看出，截滲墻水平位置越靠近海岸線一側，即1-Xb/L0越大，海水入侵回退系數越大。且當1-Xb/L0大于0時，截滲墻貫穿深度越大，海水入侵回退系數越大。當1-Xb/L0小于0時，截滲墻貫穿深度越大，海水入侵回退系數越小，且為負值。因此，當截滲墻水平位置靠近海岸線、貫穿深度較大時達到抵擋海水入侵最佳效果。

圖5 二維砂箱試驗模型不同水平位置與貫穿深度截滲墻下的回退系數R的分布Fig.5 Repulsion ratio distribution based on location and penetration of the barrier

3 龍口濱海含水層截滲墻剖面模型

3.1 海水入侵數值模型

選取位于山東龍口市黃水河附近靠近渤海灣的某典型剖面AB為研究對象（圖6）。該研究區受海水入侵影響且在距海岸線1 km處建有地下截滲墻，實際條件與前期研究的理想算例較為符合。選取剖面長度為1 500 m，厚度為20 m。含水層巖性主要為第四系松散沉積物，由黃水河河相沉積物組成。第四系下伏古近系含煤地層，巖性為砂礫巖、泥灰巖、油頁巖、黏土巖及褐煤，透水性及富水性差，構成研究區相對隔水層。研究區農業灌溉以開采地下水為主要水資源利用模式之一，同時也是影響該區域地下水動態的主要因素。

圖6 山東龍口地區研究區域示意圖Fig.6 Location of the study area in Longkou of Shandong Province

將研究區含水層概化為垂向二維、均質、各向同性的潛水含水層（圖7）。將未構建截滲墻時濱海含水層中咸淡水界面達到的穩定狀態作為初始條件。含水層上部邊界為潛水面；下部邊界為隔水邊界并且相對隔水底板存在一定的地形坡度；北部海岸線邊界設為定水頭、定濃度邊界；南部邊界設置定水頭邊界。在研究區內距海1.1 km處設置一口底部進水部位長5 m的非完整抽水井，模擬農業灌溉抽水對濱海地下水系統的影響，抽水量設為50 m3/d。采用SEAWAT-2000程序對模型進行求解，利用有限差分方法將其在空間上剖分成151列、1行、20層，差分網格的間隔為Δx=10 m、Δy=Δz=1 m，模擬期為40 d。模型主要參數見表2。

圖7 研究區典型剖面概念模型示意圖Fig.7 Conceptual model of the typical profile in Longkou

表2 山東龍口地區典型剖面地下水數值模型參數Table 2 Input parameters used in the numerical model for the typical profile in Longkou

3.2 模擬結果與分析

3.2.1 截滲墻工程措施評價

（1）截滲墻水平位置的影響

基于不同水平位置截滲墻條件下的模型結果計算海水入侵回退系數R。回退系數R與截滲墻水平位置1-X/Xw的關系圖，見圖8（a）。截滲墻的水平位置以無量綱的1-X/Xw表示，貫穿深度以無量綱的貫穿率h/H表示。其中：X表示截滲墻的水平位置，Xw表示抽水井的水平位置，h表示截滲墻貫穿深度，H表示含水層厚度。

模擬結果表明：海水入侵的回退系數與截滲墻的水平位置呈正相關關系，截滲墻水平位置越靠近海岸線一側效果越好。當截滲墻水平位置位于抽水井靠海岸線一側時（1-X/Xw>0），海水入侵回退系數為正，截滲墻對鹽水楔有驅退效果。當截滲墻水平位置位于抽水井靠岸的一側時（1-X/Xw>0），海水入侵回退系數為負數，截滲墻起反作用，海水入侵進一步加劇。

（2）截滲墻貫穿深度的影響

基于不同貫穿深度截滲墻條件下的模型結果計算了鹽水侵入楔的回退系數R，如圖8（b）所示。

模擬結果表明：當截滲墻水平位置位于抽水井靠海一側時（1-X/Xw>0），截滲墻貫穿深度越大，海水入侵回退系數越大，截滲墻效果越好。此時截滲墻有效阻擋了抽水井抽水引起的淡水水頭降低所帶來的影響。當截滲墻水平位置位于抽水井的靠岸的一側時（1-X/Xw>0），截滲墻阻擋了內陸淡水區的補給，因此截滲墻貫穿深度越大，海水入侵回退系數越小，截滲墻效果越差，甚至起到反作用。

（3）截滲墻水平位置與貫穿深度綜合影響

根據多組模擬數據，得到截滲墻不同水平位置（X）及不同貫穿深度（Z）對海水入侵回退系數（R）的等值線圖，見圖9。當截滲墻位于注水井靠海一側，且離海較近、貫穿深度較大時，驅退海水入侵的效果最佳。

圖9 龍口濱海含水層截滲墻情況下R值等值線圖Fig.9 Repulsion ratio distribution based on location and penetration of the barrier

3.2.2 實際管理措施分析

山東龍口地區所建截滲墻位于距海岸線1.0 km處，貫穿深度為18 m。依據實際截滲墻工程情況進行海水入侵的數值模擬，并與未設置截滲墻的初始狀態進行對比，模擬結果如圖7。研究表明龍口地區現有截滲墻相對初始狀態可以有效降低海水入侵程度，海水入侵前鋒向海岸線方向后退了約230 m，回退系數達23.7%。因此，現有截滲墻工程措施對保護黃水河濱海地區的淡水資源具有重要的作用。但是，對該地區截滲墻工程措施的評價分析表明，為進一步提升截滲墻工程措施的治理效果，在后續濱海地下水管理中應重新考慮優化截滲墻的水平位置與貫穿深度，適當向海岸線方向擴展以達到更好的海水入侵回退效果。

4 結論

（1）注水井工程管理措施對海水入侵的影響與補給流量及井位有關。海水入侵的回退系數與注水井流量在一定范圍內成正相關關系。注入井位在靠近鹽水楔的位置為最佳淡水注入點，并形成以該點為中心的最優注入井布局。表面補給最有效的位置是最佳注水點正上方，濱海含水層的表面補給在控制海水入侵時，不如注入井更有效。

（2）截滲墻工程管理措施對海水入侵的影響與水平位置及貫穿深度有關。截滲墻水平位置越靠近海岸線一側，海水入侵回退系數越大，截滲墻效果越好。一般情況下截滲墻貫穿深度越大，海水入侵回退系數越大，海水入侵防治效果最優。

（3）當截滲墻水平位置超過某一界限后（理想算例中該界限為鹽水楔前鋒，山東龍口實例中該界限為抽水井位置），截滲墻貫穿深度越大，海水入侵回退系數越小且為負值，截滲墻對海水入侵防治起反作用。當截滲墻水平位置靠近海岸線、貫穿深度較大時，該工程措施可以達到治理海水入侵最優效果。

（4）山東龍口黃水河地區的截滲墻工程管理措施對防治當地海水入侵具有重要的作用。但是數值模擬結果分析表明，在后續濱海地下水管理中應進一步優化截滲墻空間布局。

研究結果可為場地條件下濱海地下水管理提供實際的管理措施建議。但本文實例模型未考慮含水介質非均質性以及降雨入滲等邊界條件的動態變化對管理措施的影響。未來應構建流域尺度下龍口地區的海水入侵數值模型，基于數值模型并采用優化技術對海水入侵嚴重區進行注水井、截滲墻等海水入侵防治工程的優化管理，提出經濟、高效的濱海地下水管理建議。