潮汐波動對潛水含水層海水入侵規(guī)律的影響研究*

武雅潔，楊自良，程從敏，徐 淳

(1.中國海洋大學(xué)工程學(xué)院，山東 青島 266100；2.中國海洋大學(xué)山東省海洋工程重點實驗室，山東 青島 266100)

在自然環(huán)境下，濱海含水層咸淡水處于一種動態(tài)平衡狀態(tài)。然而，海平面的上升以及許多沿海地區(qū)地下水的過渡開采，導(dǎo)致舊的平衡被打破，海水向內(nèi)陸含水層遷移，造成海水入侵，含水層水質(zhì)惡化[1-2]。深入研究海水入侵過程的鹽度分布變化規(guī)律，對于科學(xué)管理、利用濱海地區(qū)地下水資源、防范濱海地區(qū)海水入侵危害有著重要意義。

海水入侵最關(guān)鍵的問題是研究咸淡水界面形態(tài)、位置及發(fā)展規(guī)律。早期的研究通常忽略潮汐等引起的海水面波動的影響,并假設(shè)海相邊界水位為平均海平面的固定邊界，如Ghyben和Herzberg早在100多年前就推導(dǎo)了用以確定咸淡水突變界面位置的Ghyben-Herzberg公式[3]。Henry[4]根據(jù)咸淡水混溶原理，將變密度流與鹽度輸移方程耦合求解，提出了垂直海岸線剖面上鹽度分布的解析解。

而對于濱海含水層，地下水流模式及鹽度分布規(guī)律通常受到潮汐作用的影響[5]。潮汐作用會引起含水層中地下水位的波動，這種現(xiàn)象已經(jīng)得到了較為廣泛的關(guān)注[6-8]。地下水位的波動必然會對鹽度輸移過程造成影響，目前越來越多的國內(nèi)外研究者已經(jīng)開始研究并考慮海水入侵過程中的潮汐作用，Robinson等[9]采用有限元法模擬了地下河口地下水流動和鹽分運移過程，并將模擬結(jié)果與實驗數(shù)據(jù)進行對比分析。Boufadel等[10]通過實驗和數(shù)值模擬兩種方法研究了穩(wěn)態(tài)與瞬態(tài)下溶質(zhì)的運移機理。Li等[11]考慮了浸潤面及毛細壓力等多種因素，研究了地下水與海水的循環(huán)模式。Vandenbohede等[12]通過野外觀測和地下水流模擬，研究了斜坡海岸潮汐作用下地下水流模式以及鹽度運移規(guī)律。

潮汐波動動態(tài)地改變了海相邊界條件，尤其當(dāng)海灘為小坡度傾斜海岸時，在潮汐落潮階段海灘上會形成滲出面，非線性變化的邊界使海相邊界條件變得更為復(fù)雜。近年來，國內(nèi)外學(xué)者對于濱海潛水層海水入侵的數(shù)值模擬過程中常用的一些數(shù)值模擬軟件，已經(jīng)不足以準(zhǔn)確模擬上述情形。在這方面，前人已作了一些探索，如荷蘭學(xué)者Post[13]在模擬濱海含水層海水入侵問題時,編寫了周期性邊界條件程序包,來模擬潮汐周期性波動邊界的影響,并將其嵌入MODFLOW與SEAWAT中，模擬結(jié)果與前人數(shù)據(jù)吻合較好。

本文基于OpenGeoSys科學(xué)模擬軟件，研究潮汐對潛水含水層海水入侵機制的影響，海相邊界分別設(shè)定為固定邊界和移動邊界兩種不同情形，建立斜坡海灘潛水層海水入侵的二維剖面數(shù)值模型。模擬結(jié)果表明，移動邊界能夠很好的體現(xiàn)出海相邊界隨潮汐變化的動態(tài)變化特征。在此基礎(chǔ)上探討了潛水含水層受潮汐影響下的水位波動及鹽度運移規(guī)律。

1 移動邊界處理方法

1.1 移動邊界特征

以往的海水入侵研究常將海相邊界設(shè)置為固定邊界條件，即設(shè)為定水頭和定濃度邊界。然而海相邊界常受到潮汐波動的影響，潮汐引起的海水位變化會對海水入侵過程造成影響，將海相邊界設(shè)置為固定邊界顯然是不合理的。為考慮潮汐波動作用，后來研究者們常將海相邊界概化為“L”型的垂直邊界，即假設(shè)海岸線水平(水陸交界面)位置保持不變，潮汐波動僅引起海水位在邊界位置上垂向變化，這種簡化對于傾角較大的巖質(zhì)海岸是合理的。

然而對于傾角較小的砂質(zhì)海岸(見圖1)，由于海水水位的波動不僅會造成水位垂向的變化，而且會造成海岸線位置沿水平方向移動，因此將海相邊界概化為“L”型邊界也不太合理。這種情況下，海相邊界會在垂直和水平方向都存在與潮汐相關(guān)的周期性變化，我們將這種邊界稱為移動邊界。移動邊界條件下的海水入侵問題求解更為復(fù)雜，解析法很難求出其解析解，通常采用數(shù)值模擬的方法。

圖1 移動邊界示意圖Fig.1 Schematic diagram of moving boundary

1.2 移動邊界判定

潮汐引起的海水水位變化可以分解為多個規(guī)則的簡諧振動，每個簡諧振動稱為一個分潮，實際運用中通常選取11個較大的分潮，就可以得到較為準(zhǔn)確的海面水位變化規(guī)律，即

(1)

式中：Hs(t)表示海水水位；H0表示平均海平面位置；Ai、Ti與φi分別表示各個分潮的振幅、周期和遲角。

現(xiàn)假設(shè)海平面水位為0，存在一個振幅為A、周期為T、遲角為0的潮汐波。

1.3 移動邊界條件模擬方法

本文的研究工作基于德國亥姆霍茲環(huán)境研究中心開發(fā)的OpenGeoSys科學(xué)模擬軟件，它是一個由C++實現(xiàn)、面向?qū)ο蟮慕Ｔ闯绦蚬ぞ甙?，?shù)值計算采用有限元方法，能夠模擬多孔介質(zhì)或裂隙介質(zhì)中單個或者耦合的傳熱-流動-力學(xué)-化學(xué)過程。

在模擬潮汐作用下的海水入侵過程時，采用以下方法來設(shè)定海相移動邊界條件：

(1) 每個時間步開始計算時，通過方程(1)動態(tài)地計算每個時間步的海水水位值Hs。

(2) 將海相邊界上所有離散點的初始水頭都設(shè)置為Hs，鹽度設(shè)為海水鹽度Cs。

(3) 運行OpenGeoSys科學(xué)模擬軟件，計算得到所有海相邊界離散點的壓強p與濃度C。

(4) 對于海相邊界上所有的離散點，若該點高程z

(5) 若該點高程z>Hs且壓強p>0，則該點處在滲出面上，即令該點壓強值p=0。

(6) 將與條件(2)、(3)都不滿足的海相邊界上其它點設(shè)置為零流量邊界。

(7) 重復(fù)步驟(3)～(6)，直至計算收斂，進入下一時間步。

(8) 重復(fù)步驟(2)～(7)，直至計算完成。

圖2給出了在OpenGeoSys中實現(xiàn)動邊界處理的每個時間步長的數(shù)值模擬流程圖。

圖2 動邊界模擬流程圖Fig.2 Flowing chart of moving boundary modelling

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

2.1 模型構(gòu)建

本文在Kuan[15]等物理模型試驗的基礎(chǔ)上，建立如圖3所示的實驗室尺度下的濱海潛水含水層海水入侵的數(shù)值模型。其中，潛水含水層長為3.4 m，厚度為0.7 m，假設(shè)潛水含水層介質(zhì)各向同性，底板水平隔水，不考慮降雨入滲等影響。左側(cè)為內(nèi)陸淡水邊界，有定流量的淡水補給。右側(cè)為坡度為1∶3.8的斜坡海岸，海岸水平長度為1.482 m，海岸最右端端點處高度為0.27 m。在不考慮潮汐作用時，右側(cè)海水位維持在平均海平面位置保持不變，在考慮潮汐作用時，右側(cè)海水位會隨潮汐波動而上下波動，則此時海相邊界應(yīng)設(shè)置為移動邊界。

在OpenGeoSys軟件中，建立變密度流情況下的滲流-彌散模型，即水流控制方程選取用來描述非飽和流動的Richards Flow方程，鹽度擴散過程選取溶質(zhì)輸移方程，不考慮含水層中的揮發(fā)、吸附、生物及化學(xué)反應(yīng)等作用。

水力特征函數(shù)采用van Genuchten水力特征函數(shù)，即

(2)

其中：S為飽和度；Sr表示孔隙水殘余飽和度；pc表示毛細壓強，與飽和度S度相關(guān)；a表示進氣壓強水頭比例因子，反應(yīng)巖土空隙的大??；m為指數(shù)參數(shù)，表示孔隙尺寸分布參數(shù)。

變密度流密度ρ與濃度的關(guān)系按下式計算，即

ρ(c)=ρ0[1+βc(c-c0)]。

(3)

其中：c表示溶質(zhì)質(zhì)量濃度；ρ0表示淡水密度；c0表示淡水濃度；βc表示溶質(zhì)膨脹系數(shù)。

(改自文獻[15]。Modifed from bibliography [15].)圖3 濱海潛水含水層海水入侵數(shù)值模型Fig.3 Schematic diagram of seawater intrution numerical model in unconfined coastal aquifers

本文數(shù)值模型采用有限元法進行模擬研究，將整個模擬區(qū)域通過三角形網(wǎng)格進行劃分，共有5 670個節(jié)點，包含三角形網(wǎng)格10 944個。為了更加精確體現(xiàn)海水侵入潛水含水層的整個過程，對海岸線潮間帶及潮間帶以下區(qū)域的網(wǎng)格進行局部加密，網(wǎng)格密度增大為其他區(qū)域網(wǎng)格密度的一倍。數(shù)值模型中采用的相關(guān)參數(shù)取值見表1。

表1 海水入侵模型參數(shù)取值Table 1 Parameters values adopted in seawater intrusion model

2.2 模型驗證

依據(jù)物理模型試驗分別模擬有、無潮汐兩種條件下海水入侵過程。平均海平面為0.398 m，潮汐周期T為62 s，振幅A為3.2 cm，海相邊界設(shè)為移動邊界，內(nèi)陸淡水邊界設(shè)為定流量邊界，流量為20 mL/min。并將數(shù)值模擬結(jié)果與試驗結(jié)果進行比較，以驗證數(shù)值模擬結(jié)果的有效性。為更清晰刻畫鹽度輸移結(jié)果，選取0.5鹽度等值線作為咸淡水的分界線，模擬結(jié)果如圖4所示。

模擬結(jié)果表明，不考慮潮汐作用時，海水從海相邊界處開始入侵潛水含水層，在含水層中形成單一的海水入侵鹽水楔。而考慮潮汐作用時，由于潮汐波動引起的潮間帶區(qū)域海水-地下水交換，因此會在潮間帶區(qū)域內(nèi)的潛水含水層中形成半橢圓形狀的高鹽度區(qū)，稱為上鹽水楔。上鹽水楔的形成改變了含水層中鹽度的分布規(guī)律，從圖4中可以看出，潮汐作用減小了鹽水楔向內(nèi)陸延伸的距離，這在一定程度上減緩了海水入侵的程度。數(shù)值模擬結(jié)果與Kuan等[15]的實驗結(jié)果吻合較好，表明所建模型以及移動邊界的處理方法能夠準(zhǔn)確模擬潮汐波動動態(tài)邊界條件下的濱海潛水含水層海水入侵過程。

3 潮汐對海水入侵的影響分析

為進一步探討潮汐波動對潛水含水層海水入侵規(guī)律的影響，引入不同的淡水邊界條件和潮汐條件，分析潛水含水層淡水邊界和潮汐條件對海水入侵過程的影響。不同模型條件參數(shù)取值見表2。

(黑色實線表示數(shù)值模擬0.5鹽度等值線，藍色空心點表示試驗數(shù)據(jù)。Black lines indicate 0.5 isochlors of the salt condition from the numerical simulations.Blue circles denote the laboratory experiments.)圖4 潛水含水層有無潮汐條件下的數(shù)值模擬結(jié)果Fig.4 Numerical simulation results of tide and non-tide in unconfined coastal aquifer

表2 不同條件海水入侵模擬參數(shù)值Table 2 Cases for parameters values performed for seawater intrution models

3.1 潮汐對海水入侵過程的影響

潮汐作用使得咸淡水交換過程更為頻繁，在漲潮時，海水由高潮線附近的海灘滲入潛水含水層，落潮時，由低潮線附近的海灘出流。與無潮條件相比，這種海水-地下水循環(huán)模式以及地下水位的波動必然會對海水入侵過程產(chǎn)生影響，從而影響潛水含水層中的海水入侵過程。

3.1.1 潛水含水層中的潮汐效應(yīng) 選取Q1A2 為研究對象，在模型中距離淡水邊界1.6、1.9和2.2 m位置處選取三點，來獲取地下水位的變化情況，不同位置處水頭變化幅度隨時間的關(guān)系曲線如圖5所示。

圖5 濱海潛水含水層不同位置水位波動圖Fig.5 Water table fluctuations of different locations in unconfined coastal aquifer

模擬結(jié)果表明，周期性的潮汐波動會引起潛水含水層中水位的周期性波動，且潛水含水層中水位波動周期與潮汐周期一致。位置離海岸線最近，則水位波動振幅最大，水位波動振幅隨著與海岸線距離的增大而依次減小，表明潮汐的影響在向內(nèi)陸的傳播過程中不斷衰減，越靠近內(nèi)陸，水位受潮汐的影響越小，水位波動最終會衰減為零。同時，地下水波動的相位與潮汐波動相比存在滯后性，滯后時間隨著離岸的水平距離的增加而不斷增大。

另外，地下水波動的平均水位要高于初始水位值，即潮汐作用會引起地下水位超高。Li等[7]總結(jié)得到水位超高是由于沙灘傾角導(dǎo)致的邊界條件的不對稱、浸潤面的形成以及控制方程的非線性共同作用下的結(jié)果。

3.1.2 地下淡水排放通道變化 地下淡水排放通道為地下淡水通過含水層邊界流入海洋的區(qū)域。選取Q1A2、Q1A0 為研究對象，模擬潮汐波動對地下淡水排放通道的影響，模擬結(jié)果如圖6所示。

(0表示淡水、1表示海水。0 for freshwater,1 for seawater.)圖6 潛水含水層有無潮汐條件下的淡水排放通道Fig.6 Fresh groundwater discharging zone under the condition of seawater intrusion of non tide and tide in unconfined coastal aquifer

在不考慮潮汐作用時，海平面和海岸的交點與鹽水楔和海岸的交點兩點之間的海岸為潛水含水層中地下淡水排放通道。潮汐作用不僅會改變含水層中鹽度的分布，同時也會改變含水層中地下水淡水的排放通道。

從圖6中可以看出，在考慮潮汐作用時，由于上鹽水楔的存在，上鹽水楔與下鹽水楔之間區(qū)域成為淡水排放的通道，與無潮條件下相比，潮汐作用使下鹽水楔和海岸的交點向低潮線方向移動，淡水排放通道也同時向低潮線方向移動。

3.2 潮汐振幅對海水入侵的影響

分別選取潮汐振幅不同的三種工況Q1A1、Q1A2、Q1A3為研究對象，探討潮汐振幅對海水入侵的影響，不同潮汐振幅條件下海水入侵數(shù)值模擬結(jié)果如圖7所示?？梢钥闯?，潮汐振幅大小對潛水含水層海水入侵過程有著顯著的影響。潮汐振幅越小，上鹽水楔的分布范圍越小，潛水含水層中地下水淡水排放通道向低潮線偏移距離越短。這是因為潮汐是引起潮間帶海岸區(qū)域水體循環(huán)的主要動力因素，潮汐振幅的減小減弱了潮間帶水體循環(huán)的動力因素，從而使得入侵淡水形成的上鹽水楔范圍隨之減小。同時，由于上鹽水楔的形成會限制下鹽水楔的入侵距離，上鹽水楔的減小會令下鹽水楔向內(nèi)陸移動，入侵潛水含水層的距離增大。

圖7 潛水含水層不同潮汐振幅海水入侵數(shù)值模擬結(jié)果Fig.7 Numerical simulation results of seawater intrusion with different tidal amplitude in unconfined coastal aquifer

因此，陸地邊界流量相同時，潮汐振幅的減小使得上鹽水楔收縮，下鹽水楔向陸地擴張。

3.3 陸源流量對海水入侵的影響

從圖8可以看出，在相同潮汐邊界條件下，內(nèi)陸淡水邊界的流量越小，則潛水含水層中海水入侵淡水的距離越遠。同時，內(nèi)陸邊界流量為20 mL/min時，潮汐對鹽水楔成長的影響明顯小于內(nèi)陸邊界流量為15 mL/min時，這表明在潮汐振幅相同時，陸地邊界流量越小受潮汐影響越大。這是因為當(dāng)?shù)髁繙p小時，在潮汐作用下形成的上鹽水楔會進一步在潛水含水層中擴散，而分布范圍更大的上鹽水楔會再次限制下鹽水楔的入侵距離。

圖8 潛水含水層不同淡水流量海水入侵數(shù)值模擬結(jié)果Fig.8 Numerical simulation results of seawater intrusion with different freshwater flux in unconfined coastal aquifer

因此，潮汐條件相同時，內(nèi)陸淡水邊界流量越小，海水入侵距離越遠，并且受潮汐的影響越大。

3.4 多個分潮對海水入侵的影響

潮汐引起的海水位波動是多個分潮共同作用下的結(jié)果，為進一步研究多個分潮的影響，假設(shè)向海一側(cè)同時存在兩個分潮，探討兩個分潮共同作用下的海水入侵作用，兩個分潮其中一個潮周期為62 s，振幅為4.2 cm，另一個潮周期為124 s，振幅為2.1 cm，參數(shù)取值如表2所示。

圖9給出了分別在Q1A2、Q1A4以及兩個分潮共同作用下的潛水含水層海水入侵模擬結(jié)果。當(dāng)海側(cè)潮汐振幅為0.021 m時，由于潮汐動力因素太弱，因此在模擬結(jié)果中并未在潛水含水層中形成上鹽水楔，此時的鹽度分布基本上與穩(wěn)態(tài)時保持一致。但在兩個分潮共同作用下，上鹽水楔的擴散范圍比只考慮單一分潮的二者都要大，這是由于多個分潮的存在增強了海洋中的水動力作用，促進了潮間帶海水向潛水含水層的擴散，使得上鹽水楔進一步成長，下鹽水楔入侵距離減小。

4 結(jié)論

本文建立了潮汐影響下復(fù)雜海相邊界海水入侵的數(shù)值模型，實現(xiàn)了OpenGeoSys模擬海水入侵時移動邊界的動態(tài)處理，得到了與以固定邊界為海相邊界的傳統(tǒng)海水入侵模型不同的鹽度分布規(guī)律，重點探討了潮汐波動對潛水含水層海水入侵規(guī)律的影響，達到了深入研究海水入侵機理的目的。本文得到的主要結(jié)論有以下幾個方面：

(1) 潮汐波動會引起潛水含水層地下水水位超高；潮汐引起的水位波幅隨著離岸的水平距離的增加而不斷衰減，相位與潮汐波相比存在滯后，滯后時間隨著離岸的水平距離的增加而不斷增大。

(2) 斜坡海灘上的潮汐波動會引起潮間帶區(qū)域相對較快的水體循環(huán)，這種海水-淡水循環(huán)促使了上鹽水楔的形成，上鹽水楔的存在改變了海水入侵過程中鹽度的分布及淡水排放通道，并且潮汐的存在會在一定程度上減緩海水入侵程度。

(3) 陸地邊界流量相同時，潮汐振幅的減小使得上鹽水楔收縮，下鹽水楔向陸地擴張。

(4) 潮汐條件相同時，內(nèi)陸淡水邊界流量越小，海水入侵距離越遠，并且受潮汐的影響越大。

(5) 多個分潮的存在使得上鹽水楔加劇擴散，下鹽水楔向海收縮。

本文只進行了實驗室尺度下物理模型的數(shù)值模擬分析，對于實際尺度下的應(yīng)用，需要考慮更為復(fù)雜的降雨入滲、地下水開采、各向異性介質(zhì)等問題，還有待進一步研究。