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

(1.中國海洋大學(xué)工程學(xué)院，山東 青島 266100；2.中國海洋大學(xué)山東省海洋工程重點(diǎn)實(shí)驗(yàn)室，山東 青島 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]采用有限元法模擬了地下河口地下水流動和鹽分運(yùn)移過程，并將模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比分析。Boufadel等[10]通過實(shí)驗(yàn)和數(shù)值模擬兩種方法研究了穩(wěn)態(tài)與瞬態(tài)下溶質(zhì)的運(yùn)移機(jī)理。Li等[11]考慮了浸潤面及毛細(xì)壓力等多種因素，研究了地下水與海水的循環(huán)模式。Vandenbohede等[12]通過野外觀測和地下水流模擬，研究了斜坡海岸潮汐作用下地下水流模式以及鹽度運(yùn)移規(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é)模擬軟件，研究潮汐對潛水含水層海水入侵機(jī)制的影響，海相邊界分別設(shè)定為固定邊界和移動邊界兩種不同情形，建立斜坡海灘潛水層海水入侵的二維剖面數(shù)值模型。模擬結(jié)果表明，移動邊界能夠很好的體現(xiàn)出海相邊界隨潮汐變化的動態(tài)變化特征。在此基礎(chǔ)上探討了潛水含水層受潮汐影響下的水位波動及鹽度運(yùn)移規(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ī)則的簡諧振動，每個簡諧振動稱為一個分潮，實(shí)際運(yùn)用中通常選取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++實(shí)現(xiàn)、面向?qū)ο蟮慕Ｔ闯绦蚬ぞ甙?，?shù)值計(jì)算采用有限元方法，能夠模擬多孔介質(zhì)或裂隙介質(zhì)中單個或者耦合的傳熱-流動-力學(xué)-化學(xué)過程。

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

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

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

(3) 運(yùn)行OpenGeoSys科學(xué)模擬軟件，計(jì)算得到所有海相邊界離散點(diǎn)的壓強(qiáng)p與濃度C。

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

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

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

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

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

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

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

2 海水入侵?jǐn)?shù)值模型

2.1 模型構(gòu)建

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

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

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

(2)

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

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

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

(3)

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

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

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

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

2.2 模型驗(yàn)證

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

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

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

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

(黑色實(shí)線表示數(shù)值模擬0.5鹽度等值線，藍(lán)色空心點(diǎn)表示試驗(yàn)數(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位置處選取三點(diǎn)，來獲取地下水位的變化情況，不同位置處水頭變化幅度隨時間的關(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

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

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

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

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

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

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

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

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

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

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

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

潮汐引起的海水位波動是多個分潮共同作用下的結(jié)果，為進(jìn)一步研究多個分潮的影響，假設(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)時保持一致。但在兩個分潮共同作用下，上鹽水楔的擴(kuò)散范圍比只考慮單一分潮的二者都要大，這是由于多個分潮的存在增強(qiáng)了海洋中的水動力作用，促進(jìn)了潮間帶海水向潛水含水層的擴(kuò)散，使得上鹽水楔進(jìn)一步成長，下鹽水楔入侵距離減小。

4 結(jié)論

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

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

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

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

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

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

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