劉尚辰，孫昭晨，梁書秀

(大連理工大學 海岸和近海工程國家重點實驗室，大連 116024)

新村潟湖位于海南島東南部的陵水縣，是一個典型的半封閉沙壩-潮汐汊道-潟湖海岸體系，地貌格局復(fù)雜，地形變化較大[1]。新村潟湖通過單一口門與外海相連，潟湖口門處北岸為新村碼頭堤岸，南岸是基巖組成的南灣猴島，口門夾于二者之間，位置比較固定，穩(wěn)定性高，口門附近漲潮三角洲、落潮三角洲有所發(fā)育[2]。受上述地貌影響，新村潟湖內(nèi)部水動力條件溫和。同時潟湖地處熱帶，常年水溫、鹽度適中，因而擁有豐富的浮游生物，潟湖內(nèi)水產(chǎn)養(yǎng)殖活動較多[3-4]。因大量養(yǎng)殖塘向潟湖水域擴張、養(yǎng)殖廢水直接排放、養(yǎng)殖密度過大、生活垃圾污染等問題，新村潟湖水質(zhì)下降，局部水體富營養(yǎng)化程度增高，赤潮時有發(fā)生[5-6]，潟湖環(huán)境已受到較大影響。要實現(xiàn)對新村潟湖的合理開發(fā)利用，需要對潟湖內(nèi)水動力環(huán)境進行深入研究。目前，針對新村潟湖潮汐汊道演變[2,7-9]、潮汐汊道及潟湖二者水動力關(guān)系[1,10-12]的研究較多，關(guān)于新村潟湖水體交換方面的研究較為少見。

海洋對污染物本身具有巨大的自凈能力[13]，潟湖中水體的自凈能力主要體現(xiàn)在由漲、落潮過程產(chǎn)生的水體交換中，潟湖內(nèi)部污染物在該過程中得到稀釋和擴散后輸運出潟湖，因此水體交換能力本質(zhì)上反映了一定的水體自凈能力[14-16]。對于水體交換能力和與之相關(guān)的水體自凈能力，國內(nèi)外已開展了大量工作，常用的計算模型有箱式模型[17-18]、隨機游動模型[19-20]、對流擴散模型[21-22]、時間尺度模型[23]等。本文擬在新村潟湖已有自然資料基礎(chǔ)上，通過MIKE21數(shù)值模型建立新村潟湖及臨近海域的潮流數(shù)值模型，利用該模型對不同潮汐過程作用下潟湖內(nèi)保守物質(zhì)的輸運擴散過程以及粒子運動過程進行數(shù)值模擬，從水體自凈速率[24]、半交換時間和停留時間3個方面入手，對不同潮汐過程作用下新村潟湖的水體交換能力進行探討。

1 數(shù)值模型

1.1 控制方程

本文采用由丹麥水資源及水環(huán)境研究所(DHI)開發(fā)的水動力模型MIKE21中水動力模塊、輸運模塊以及粒子追蹤模塊進行計算?？刂品匠滩捎萌S不可壓縮雷諾平均下的Navier-Stokes方程，在服從Boussinesq和靜水壓力假定條件下，對水平方向上的動量方程和連續(xù)性方程沿水深方向進行積分平均處理。不同坐標系下，控制方程可呈不同表現(xiàn)形式，這里統(tǒng)一采用笛卡爾坐標系形式進行說明。水動力模塊控制方程如下

連續(xù)性方程

(1)

動量方程

(2)

(3)

(4)

(5)

(6)

式中：τsx和τsy為表面風應(yīng)力分量；τbx和τby為底部切應(yīng)力分量；sxx、sxy、syx、syy為輻射應(yīng)力分量；Txx、Txy、Tyx、Tyy為粘滯應(yīng)力分量。

輸運模塊控制方程如下

(7)

粒子追蹤模塊控制方程如下

dXt=a(t,Xt)dt+b(t,Xt)ξtdt

(8)

式中：a(t,Xt)為因流、風、底摩擦等影響產(chǎn)生的漂流移動項；b(t,Xt)為因分子運動、湍流等產(chǎn)生的擴散移動項；ξt為隨機數(shù)；Xt為粒子位置。在該模塊中采用Langevin方程作為計算基礎(chǔ)，離散后方程如下

Yn+1=Yn+a(t,Xt)YnΔn+b(t,Xt)YnΔWn

(9)

1-a 選自Google Earth 1-b 選自Landsat影像數(shù)據(jù)圖1 新村潟湖地理位置Fig.1 The location of Xincun Lagoon

式中：Y表示軌跡位置；ΔWn=Wt-Ws∈N(μ=0，σ2=Δn)為服從標準高斯分布的維納過程W的增量。

1.2 模型計算域及邊界條件

新村潟湖地處海南陵水縣東南側(cè)，面向南海，地理位置見圖1。計算區(qū)域包含新村潟湖及鄰近外海，范圍大致在109.89°E～110.04°E、18.33°N～18.46°N。數(shù)值模型采用單元中心有限體積法進行離散，利用不重疊的三角形單元進行計算，在口門、潮汐汊道附近[11]及形狀復(fù)雜固邊界處對網(wǎng)格進行加密，網(wǎng)格總數(shù)為28 316，節(jié)點數(shù)為14 532，模型計算網(wǎng)格及驗證點分布見圖2。

圖2 模型計算網(wǎng)格、地形及驗證點分布Fig.2 The mesh grids, terrain and verification points of model

模型邊界包含固邊界和開邊界，即陸水邊界和水水邊界。固邊界規(guī)定垂直于海岸邊界流速為零，開邊界采用水位邊界進行驅(qū)動，由于新村潟湖無明顯徑流流入[2]，模型不考慮徑流影響。開邊界水位數(shù)據(jù)由中國近海潮汐潮流數(shù)學模型Chinatide計算得到[25]，Chinatide模型可預(yù)測給定位置下基于Q1、P1、O1、K1、N2、M2、S2、K2、Sa等9個主要分潮調(diào)和常數(shù)的潮位變化，模型精度5′×5′，模型范圍可覆蓋中國各海域，并且在模型開發(fā)時已通過各海域長期驗潮站資料對模型進行計算精度檢測。

1.3 模型驗證

由于缺乏近期新村潟湖實測資料，根據(jù)潟湖水動力方面已有研究資料，對小潮期間、大潮期間以及時長半個月潮汐作用下的模型進行驗證[10,12,26]，驗證點分布見圖2，其中S1、S2、S3驗證點用于小潮和大潮期間模型驗證，S4驗證點用于為期半個月潮汐作用下的模型驗證。計算發(fā)現(xiàn)，小潮和大潮期間驗證點位置處模擬值與實測值存在0～2 h左右相位差，量值模擬較好，相位調(diào)整后的小潮、大潮期間驗證結(jié)果見圖3；為期半個月潮汐作用下的水位模擬結(jié)果與測量值吻合良好，驗證結(jié)果見圖4。小潮和大潮期間的漲落潮過程基本符合水體運動情況，由于潟湖地形復(fù)雜，水上水下建筑物、 植被均會對流場產(chǎn)生一定影響，因此流速模擬有些偏差?？傮w而言，模型能夠準確反映水動力特性。

圖3 潮汐驗證結(jié)果(25 h)Fig.3 Model verification during tide period(25 h)

圖4 潮汐驗證結(jié)果(15 d)Fig.4 Model verification during tide period(15 d)

1.4 水體交換計算方案

新村潟湖水體交換能力與水動力條件關(guān)系密切，該區(qū)域水動力主要受潮汐因素控制[2]。潮汐漲落現(xiàn)象主要由天體的周期性運動引起，由于太陽、月球、地球運轉(zhuǎn)狀況不同，潮汐會出現(xiàn)如日不等、月不等、年不等等現(xiàn)象，因而與之相關(guān)的新村潟湖水體交換過程也會有所不同。潟湖整體水體交換是一個時間較長的過程，有時長達數(shù)月，故本文擬對潟湖長期水體交換過程進行計算，并由此展開討論。

已有研究表明該海域的潮汐類型為不規(guī)則的全日混合潮[27]，由于缺乏長期觀測資料，并且 Chinatide數(shù)值模型計算考慮分潮數(shù)較多，模擬得到的潮位應(yīng)用在潟湖驗證計算中效果較好，故采用Chinatide數(shù)值模型對該區(qū)域長期潮位進行計算。圖5為2017～2019年的潮位計算數(shù)據(jù)和潮位包絡(luò)線；圖6包含了季度平均水位過程線和潮差變化曲線，季度平均水位過程線表示包含該天在內(nèi)的前后共90 d的平均歷時水位過程線。受天文因素影響，潮差和平均水位具有長周期性、季節(jié)性等特征，由圖5、圖6可見，該地潮差在春秋季整體變化范圍較小，而在冬夏季整體變化較大。結(jié)合新村潟湖當?shù)爻毕^程、潮差和季度平均水位的長期變化特征，選用為期40 d的具有一定對稱性和相似性的4種潮汐過程對該區(qū)域進行研究，以考慮不同潮汐過程作用對水體交換的影響，所選過程所在時間段見圖6，詳細過程曲線見圖7。過程1-1和過程1-2取自秋春季，二者分別與季度平均高水位和季度平均低水位相對應(yīng)，期間各階段潮差相差較小、潮差的歷時變化范圍較??；過程2-1和過程2-2分別取自冬夏季，二者分別與季度平均高水位和季度平均低水位相對應(yīng)，期間各階段潮差相差較大、潮差歷時變化范圍較大，通過以上所選4種潮汐過程，可考慮季度平均水位和長期潮差變化對潟湖水體交換的影響。

圖7 典型潮汐過程Fig.7 Typical tidal process

由于潟湖內(nèi)部地形較復(fù)雜，水動力條件各區(qū)域差別較大，故在空間上水體交換能力可能會有一定差異，因此采用以歐拉法為基礎(chǔ)的物質(zhì)輸運方法和以拉格朗日法為基礎(chǔ)的粒子追蹤方法對潟湖整體水體交換能力及其空間分布開展研究。兩種方法初始條件設(shè)置如下：將計算區(qū)域劃分為兩部分，以口門為界，邊界內(nèi)部潟湖區(qū)域內(nèi)設(shè)置保守物質(zhì)濃度為1，邊界以外濃度為0，同時在潟湖區(qū)域內(nèi)隨機均勻撒入保守物質(zhì)粒子。初始濃度設(shè)置及粒子布置見圖8。

圖8 初始濃度設(shè)置及粒子布置Fig.8 Initial distribution of concentration and particles

2 結(jié)果與討論

2.1 水體自凈速率

水體交換能力與水體自凈能力密切相關(guān)，通過水體自凈速率可以反映出潟湖整體水體交換能力。水體自凈速率代表在區(qū)域自然、物理因素作用下，水體所含某物質(zhì)總量對時間的變化率，該指標包含了水體交換和物質(zhì)自身擴散性的影響，水體自凈速率表達式如下[24]

(10)

式中：M表示水體中某保守物質(zhì)總量；Ml(t)表示水體自凈速率，Ml(t)>0時代表該保守物質(zhì)由研究區(qū)域輸運至域外，為方便比較，對M做歸一化處理，即M=M(t)/M(0)。圖9-a給出了4種潮汐過程作用下的水體自凈速率；為比較不同季度平均水位的影響，圖9-b分別將潮差變化相近的過程1-1和過程1-2、過程2-1和過程2-2的水體自凈速率曲線置于一幅圖中。從圖9-a可見，潟湖水體自凈速率曲線在一定程度上與潮汐過程曲線形狀相似，同時其峰值變化與高潮位及低潮位的變化趨勢相近，可見該區(qū)水體交換過程與潮汐運動密切相關(guān)；在潮差較大階段，自凈速率峰值相對較大，潮差較小階段，自凈速率峰值相對較小，說明自凈速率與潮差具有一定的正相關(guān)關(guān)系，并且自凈速率正向峰值大多數(shù)遠大于臨近的負向峰值，可見區(qū)域內(nèi)物質(zhì)具有不斷向外輸運和交換的趨勢。從圖9-b中可見，同一圖中的兩個潮汐過程自凈速率曲線變化差別不大，而不同圖中的兩條曲線則呈現(xiàn)較大差別，說明不同季度平均水位對自凈速率的影響有限，而潮汐歷時變化對自凈速率有較大影響；潮汐過程2-1、2-2瞬時自凈速率峰值在潮汐作用前300 h明顯大于潮汐過程1-1、1-2，而在潮汐作用后期，由于在潮汐過程1-1、1-2作用下保守物質(zhì)仍有很多沒有輸運出潟湖，二者的瞬時自凈速率峰值仍相對較大。

因潮汐往復(fù)運動影響，由圖9可見，水體自凈速率曲線也在正負向不斷交替變化，取潮汐過程曲線中每次高潮位對應(yīng)的時間，計算相應(yīng)時間潟湖內(nèi)剩余保守物質(zhì)占初始總量的比例P，如圖10-a所示，在潮汐過程2-1、2-2作用下，P在前期迅速下降，并且作用時間相同下，P明顯小于潮汐過程1-1、1-2作用下的值，可見在大潮差潮汐階段物質(zhì)更易向外輸運；由于潮汐過程2-1、2-2與1-1、1-2相比，小潮階段潮差小、歷時短、潮差變化相對小，故P曲線呈現(xiàn)出更為明顯的兩段平臺段，此時水體交換速率較慢。按照水體自凈速率表達式，規(guī)定此時自凈速率計算式如下

(11)

計算得到的自凈速率見圖10-b，4種過程的自凈速率均逐漸趨向于0。在550～800 h階段，潮汐過程1-2由于平均水位較低，潮差和潮差歷時變化均較小，因而高潮位和低潮位的納潮量變化較小，故自凈速率在正負向往復(fù)變化并且幅值相差不大，此時物質(zhì)凈輸運情況不佳，該階段整體表現(xiàn)出的水體交換能力較弱，反映在P值上，可見該階段P值下降緩慢。

整體來看，在4種不同潮汐過程作用下，新村潟湖內(nèi)剩余保守物質(zhì)總量可在1 000 h內(nèi)下降到原總量40%以下，潮汐過程2-1、2-2作用下前期自凈速率變化較快，在500 h以內(nèi)便可達到該值，表現(xiàn)出很強的水體交換能力，1-2過程作用下保守物質(zhì)輸運出潟湖的速度最為緩慢，水體交換能力較弱。潮汐過程對水體交換能力有較強影響，在相同作用時間下，潮差較大過程，由于水體高低潮納潮量相差較大，物質(zhì)向外輸運能力較強，水體交換能力較強；與高低潮位歷時變化不同的是自凈速率峰值變化前期呈現(xiàn)出較大的非對稱性，該非對稱性與潮差歷時變化和潮汐非對稱性有關(guān)，由潮汐過程可見，該區(qū)漲潮歷時略大于落潮歷時，且所選過程前期均為潮差逐漸增大階段，落潮潮差大于漲潮潮差，物質(zhì)向外輸運趨勢更加明顯，故初期瞬時自凈速率峰值快速增長。

2.2 水體半交換時間

為研究潟湖水體交換能力在空間上的分布特征，采用時間尺度模型進行后續(xù)研究，各種時間尺度可以從不同角度描述海水交換的快慢程度。Luff等[28]提出了半交換時間的概念，即海域內(nèi)某種保守物質(zhì)通過對流擴散作用將濃度稀釋到初始濃度一半所需的時間。受潮汐影響，水體交換過程呈現(xiàn)出往復(fù)性，各位置濃度也具有一定往復(fù)性特征，為減少該影響，規(guī)定第一次在高潮位時濃度達到初始濃度一半所用的時間為潟湖該位置處的半交換時間。采用該水體半交換時間的概念，對4種潮汐過程作用下潟湖各區(qū)域的半交換時間進行了計算，結(jié)果見圖11。從圖11可見，半交換時間大致可分為0～200 h、200～400 h、400～800 h以及800 h以上4部分，其中半交換時間為0～400 h和800 h以上區(qū)域的面積相當，400～800 h的區(qū)域面積最小，且在潟湖北部半交換時間過渡較弱，潟湖南部過渡較明顯；靠近口門處的區(qū)域整體半交換時間很短，大致在400 h以下，水體交換能力強，潟湖中部北岸淺灘較多，地形復(fù)雜，水動力較弱，半交換時間明顯高于其他部分；在4種潮汐過程作用下，遠離口門的次級潟湖及附近區(qū)域濃度均未達到0.5以下，水體交換能力較弱。從不同潮汐過程的作用結(jié)果來看，相同作用時間長度，潮汐過程2-1、2-2作用下濃度降到0.5以下的區(qū)域明顯大于潮汐過程1-1、1-2作用下的區(qū)域；同時半交換時間空間分布在潟湖中部呈現(xiàn)渦旋形狀，一定程度上體現(xiàn)了物質(zhì)凈輸運的歷時過程。整體來看，從半交換時間的空間分布反映出，口門附近、潟湖中部和南部深水區(qū)水體交換能力較強，其中口門附近交換能力最強，而受地形和距離口門遠近的影響，潟湖中部北岸、最遠端次級潟湖及其與中部深水區(qū)連接處水體交換能力較弱。

圖11 水體半交換時間Fig.11 Half-life time of substance in lagoon

2.3 停留時間

水體交換過程在微觀層面上表現(xiàn)的即是物質(zhì)粒子或粒子團的輸移交換過程，本文采用停留時間對該角度體現(xiàn)出的水體交換能力進行了分析。為簡化計算，在初始時刻，于潟湖內(nèi)部隨機均勻布置保守型物質(zhì)粒子作為示蹤劑，追蹤每個粒子直到離開潟湖為止，當標識粒子首次到達研究區(qū)外時，認為在該粒子初始位置處的水質(zhì)點或物質(zhì)粒子與外界進行了交換，所用的時間即為停留時間。計算中發(fā)現(xiàn)部分粒子運動到了淺灘、復(fù)雜邊界等位置，受邊界條件影響無法繼續(xù)移動，因此對其進行了剔除，處理后的停留時間結(jié)果見圖12。

圖12 停留時間Fig.12 Residence time

從結(jié)果可見，靠近口門的區(qū)域停留時間在200 h以內(nèi)，粒子運動出潟湖時間最短，水體交換速度最快，并且口門內(nèi)漲潮三角洲對交換過程影響較小；距離口門較近的潟湖中部停留時間在200～400 h，水體交換速度略低于口門區(qū)域；在遠端次級潟湖及其與潟湖中部的連接處，停留時間整體在800 h以上，部分該區(qū)域粒子仍可輸運出潟湖；由于次級潟湖末端距離口門最遠，并且注入潟湖徑流過小，單純潮汐作用下，在40 d內(nèi)粒子無法輸運出潟湖，水體交換能力最弱；潟湖南部由于存在較大范圍淺灘區(qū)，且距離口門較遠，水動力條件很弱，在4種潮汐過程作用下均有大量粒子停留時間在800 h以上，水體交換速率很低；從潮汐過程1-1、1-2和2-1、2-2作用下的結(jié)果來看，潟湖停留時間存在明顯分區(qū)，如圖13所示，A區(qū)和B區(qū)平均停留時間結(jié)果見表1。

圖13 潟湖分區(qū)Fig.13 The division of lagoon

表1 各分區(qū)停留時間Tab.1 The resistance time of different divisions h

A區(qū)停留時間明顯低于B區(qū)，表明A區(qū)水體交換速率高于B區(qū)，并且潮汐過程2-1、2-2與潮汐過程1-1、1-2相比，在其作用下A區(qū)平均停留時間更短，水體交換速率更快。從停留時間和半交換時間二者的結(jié)果來看，二者在反映水體交換能力上有一定差別，停留時間以拉格朗日法為基礎(chǔ)，在粒子輸運過程中含有運動隨機性，但不包含與域外水體的混合過程，因此在水交換能力較弱區(qū)域，停留時間結(jié)果反映出的水體交換能力比半交換時間反映出的水體交換能力要弱。

3 結(jié)論

(1)經(jīng)計算可見，自凈速率可以反映出潟湖內(nèi)物質(zhì)的輸運趨勢，一個漲落潮階段，自凈速率為正時，物質(zhì)隨水體流出，自凈速率為負時，物質(zhì)隨水體流入，當該階段自凈速率積分為正時，代表潟湖內(nèi)物質(zhì)在該階段結(jié)束時被交換到了域外。從新村潟湖水體自凈速率來看，季度平均水位對新村潟湖水體交換能力影響較小，而潮差大小和潮差歷時變化對其影響較大；潮汐過程2-1、2-2下的瞬時自凈速率在作用前期峰值明顯大于潮汐過程1-1、1-2，潟湖內(nèi)保守物質(zhì)總量在潮汐過程2-1、2-2作用下也下降較快，可見在潮汐過程2-1、2-2對應(yīng)的冬夏季新村潟湖表現(xiàn)出了很強的水體交換能力。

(2)4種潮汐過程作用下，潟湖內(nèi)半交換時間可分為0～400 h、400～800 h和800 h以上三類區(qū)域，其中0～400 h和800 h以上區(qū)域面積相當，400～800 h區(qū)域面積最小，潮汐過程2-1、2-2作用下，可進行較強水體交換區(qū)域的面積大于其余兩個潮汐過程?？梢娦麓鍧暫w交換能力在空間上分布不均勻，大體上以口門為中心，半交換時間長短呈扇形分布，距離口門越近，受口門束窄影響，區(qū)域流速較快，潮流漲落通暢，半交換時間很短，水體交換能力很強，而遠離口門及水深較淺區(qū)域，水動力條件很弱，半交換時間很長，表現(xiàn)出較差的水體交換能力。

(3)停留時間的計算結(jié)果可將新村潟湖分為兩個水體交換能力明顯相異的區(qū)域，其與半交換時間得到的0～400 h和800 h以上兩個分區(qū)有一定相似性。同時，基于拉格朗日法計算的停留時間和基于歐拉法計算的半交換時間在表示潟湖水體交換能力上有所不同，由于停留時間的計算基于粒子運動，沒有表現(xiàn)出很強的水體混合的過程，故在水深較淺但離口門較近區(qū)域如潟湖中部南岸，受地形及水動力條件影響，該位置處停留時間表現(xiàn)出的水體交換能力要弱于半交換時間表現(xiàn)出的能力；從不同潮汐過程角度來看，潮汐過程2-1、2-2作用下，停留時間大于800 h的區(qū)域更偏潟湖末端，此時水體交換能力較強。