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

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

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

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

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

1 研究方法

1.1 海水入侵模型

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

式中：ρ—流體密度；

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

ρf—淡水密度；

Z—計算點高程；

hf—等效淡水水頭；

Sf—等效淡水單位儲水系數(shù)；

t—時間；

θ—有效孔隙度；

C—溶質(zhì)濃度；

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

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

Ck—溶質(zhì)k的溶解濃度；

Dij—水動力彌散張量；

vi—滲流速度；

—源（匯）流中溶質(zhì)k的濃度；

Rn—化學反應項。

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

1.2 海水入侵管理措施

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

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

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

2 二維砂箱數(shù)值模型

2.1 試驗模型

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

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

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

2.2 模擬結(jié)果與分析

2.2.1 注水井工程措施評價

（1）注水井流量

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

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

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

（2）注水井位置

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

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

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

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

2.2.2 截滲墻工程措施評價

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

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

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

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

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

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

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

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

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

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

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

3.1 海水入侵數(shù)值模型

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

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

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

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

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

3.2 模擬結(jié)果與分析

3.2.1 截滲墻工程措施評價

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

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

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

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

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

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

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

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

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

3.2.2 實際管理措施分析

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

4 結(jié)論

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

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

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

（4）山東龍口黃水河地區(qū)的截滲墻工程管理措施對防治當?shù)睾Ｋ肭志哂兄匾淖饔?。但是?shù)值模擬結(jié)果分析表明，在后續(xù)濱海地下水管理中應進一步優(yōu)化截滲墻空間布局。

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