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

(1.中國(guó)海洋大學(xué)工程學(xué)院，山東 青島 266100；2.中國(guó)海洋大學(xué)山東省海洋工程重點(diǎn)實(shí)驗(yàn)室，山東 青島 266100)

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

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

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

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

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

1 移動(dòng)邊界處理方法

1.1 移動(dòng)邊界特征

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

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

圖1 移動(dòng)邊界示意圖Fig.1 Schematic diagram of moving boundary

1.2 移動(dòng)邊界判定

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

(1)

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

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

1.3 移動(dòng)邊界條件模擬方法

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

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

(1) 每個(gè)時(shí)間步開(kāi)始計(jì)算時(shí)，通過(guò)方程(1)動(dòng)態(tài)地計(jì)算每個(gè)時(shí)間步的海水水位值Hs。

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

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

(4) 對(duì)于海相邊界上所有的離散點(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)都不滿(mǎn)足的海相邊界上其它點(diǎn)設(shè)置為零流量邊界。

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

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

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

圖2 動(dòng)邊界模擬流程圖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ù)值模型。其中，潛水含水層長(zhǎng)為3.4 m，厚度為0.7 m，假設(shè)潛水含水層介質(zhì)各向同性，底板水平隔水，不考慮降雨入滲等影響。左側(cè)為內(nèi)陸淡水邊界，有定流量的淡水補(bǔ)給。右側(cè)為坡度為1∶3.8的斜坡海岸，海岸水平長(zhǎng)度為1.482 m，海岸最右端端點(diǎn)處高度為0.27 m。在不考慮潮汐作用時(shí)，右側(cè)海水位維持在平均海平面位置保持不變，在考慮潮汐作用時(shí)，右側(cè)海水位會(huì)隨潮汐波動(dòng)而上下波動(dòng)，則此時(shí)海相邊界應(yīng)設(shè)置為移動(dòng)邊界。

在OpenGeoSys軟件中，建立變密度流情況下的滲流-彌散模型，即水流控制方程選取用來(lái)描述非飽和流動(dòng)的Richards Flow方程，鹽度擴(kuò)散過(guò)程選取溶質(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)行模擬研究，將整個(gè)模擬區(qū)域通過(guò)三角形網(wǎng)格進(jìn)行劃分，共有5 670個(gè)節(jié)點(diǎn)，包含三角形網(wǎng)格10 944個(gè)。為了更加精確體現(xiàn)海水侵入潛水含水層的整個(gè)過(guò)程，對(duì)海岸線潮間帶及潮間帶以下區(qū)域的網(wǎng)格進(jìn)行局部加密，網(wǎng)格密度增大為其他區(qū)域網(wǎng)格密度的一倍。數(shù)值模型中采用的相關(guān)參數(shù)取值見(jiàn)表1。

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

2.2 模型驗(yàn)證

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

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

3 潮汐對(duì)海水入侵的影響分析

為進(jìn)一步探討潮汐波動(dòng)對(duì)潛水含水層海水入侵規(guī)律的影響，引入不同的淡水邊界條件和潮汐條件，分析潛水含水層淡水邊界和潮汐條件對(duì)海水入侵過(guò)程的影響。不同模型條件參數(shù)取值見(jiàn)表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 潛水含水層有無(wú)潮汐條件下的數(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 潮汐對(duì)海水入侵過(guò)程的影響

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

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

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

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

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

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

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

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

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

3.2 潮汐振幅對(duì)海水入侵的影響

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

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

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

3.3 陸源流量對(duì)海水入侵的影響

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

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

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

3.4 多個(gè)分潮對(duì)海水入侵的影響

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

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

4 結(jié)論

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

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

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

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

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

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

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