磨刀門鹽水入侵特性及其動(dòng)力機(jī)制研究

呂紫君，馮佳佳，孔 俊，羅照陽，潘明婕

(1.河海大學(xué) 江蘇省海岸海洋資源開發(fā)與環(huán)境安全重點(diǎn)試驗(yàn)室，南京 210098;2.南京市水利規(guī)劃設(shè)計(jì)院股份有限公司，南京 210006)

鹽水入侵是沿海河口地區(qū)特有的一種水文現(xiàn)象，多發(fā)于枯水期，是河口研究中關(guān)鍵問題之一。國(guó)外關(guān)于鹽水入侵的研究從20世紀(jì)50年代開始，早期研究工作主要集中在河口環(huán)流、鹽水入侵長(zhǎng)度、鹽淡水混合類型、最大渾濁帶等[1-3]。河口區(qū)的鹽淡水混合與輸運(yùn)是一個(gè)極其復(fù)雜的過程，鹽水入侵受向陸和向海兩種鹽度輸送驅(qū)動(dòng)力共同控制；Hansen[2]較早開始采用機(jī)制分解方法研究河口鹽度輸運(yùn)，提出以河口流速及鹽度的垂、橫向變化來描述河口的鹽度通量；Lerczak[4]根據(jù)實(shí)測(cè)資料研究了哈德遜河口的鹽度輸運(yùn)，將縱向鹽度輸運(yùn)分成平流輸運(yùn)、穩(wěn)定剪切輸運(yùn)及潮汐震動(dòng)輸運(yùn)。國(guó)內(nèi)關(guān)于鹽水入侵的研究從20世紀(jì)80年代開始，主要集中在長(zhǎng)江口[5]和珠江口[6]。

珠江口磨刀門水道沿途分布著眾多取水口，是珠海、中山、澳門等城市的主要飲用水水源地。21世紀(jì)以來，受航道疏浚、口門挖沙以及上游徑流量變化的影響，磨刀門鹽水入侵日趨嚴(yán)重，眾多學(xué)者針對(duì)這一問題開展了大量研究。賈良文等[7]對(duì)磨刀門水流及鹽度的變化過程進(jìn)行分析，得出磨刀門枯季鹽淡水類型屬于緩混合型；韓志遠(yuǎn)等[8]基于不同年份的潮位和地形資料，分析了磨刀門鹽水入侵加劇的原因；盧陳等[9]采用物理模型實(shí)驗(yàn)方法對(duì)不同潮差驅(qū)動(dòng)下鹽水入侵距離進(jìn)行研究，結(jié)果表明存在潮差臨界值使得鹽水入侵距離最短。由于實(shí)地測(cè)量的時(shí)空覆蓋有限，而物理試驗(yàn)的研究成本高，所以數(shù)值模擬已廣泛應(yīng)用于河口鹽水入侵的研究，Gong和Shen[10]采用ELCIRC模型對(duì)磨刀門鹽水入侵與徑流和潮差的響應(yīng)法則進(jìn)行研究，表明徑流量的大小是影響鹽水入侵的關(guān)鍵因素；Wang等[11]、陳文龍等[12]利用FVCOM模型對(duì)磨刀門水道的鹽水入侵規(guī)律進(jìn)行了分析，發(fā)現(xiàn)鹽水入侵最遠(yuǎn)出現(xiàn)在小潮之后的中潮。

上述成果多為對(duì)磨刀門鹽水入侵規(guī)律與影響因素的研究，相應(yīng)的物質(zhì)輸運(yùn)通量及其分解的研究偏少，鹽水入侵半月周期變化的動(dòng)力學(xué)機(jī)制尚不明確，同時(shí)磨刀門鹽度數(shù)值模擬的精度及效率有待進(jìn)一步提高。鑒于此，本文基于一種最新發(fā)展的河口海洋模型SCHISM[13]，構(gòu)建磨刀門水道三維水流鹽度數(shù)值模型，利用該模型對(duì)磨刀門鹽水入侵特性、縱向鹽度輸運(yùn)進(jìn)行分析，探究磨刀門鹽水入侵的動(dòng)力機(jī)制。

1 數(shù)值模型建立

1.1 數(shù)值模型

研究中采用的模型為SCHISM[13]，該模型是基于SELFE[14]開發(fā)出的一個(gè)跨尺度湖泊-河流-河口-海洋數(shù)值模型，計(jì)算區(qū)域水平方向上采用三角形和四邊形混合網(wǎng)格，可以靈活地?cái)M合復(fù)雜的地形和岸線；垂向上采用傳統(tǒng)的混合SZ坐標(biāo)或者新增的LSC2坐標(biāo)[15]，混合SZ坐標(biāo)不僅能適應(yīng)復(fù)雜不規(guī)則地形，而且能有效改善模型垂向的空間分辨率，克服河口區(qū)因地形變化大而引起的模擬水平梯度失真的難題；而LSC2坐標(biāo)在保證深海區(qū)垂向網(wǎng)格精度的同時(shí)，避免淺水區(qū)垂向?qū)訑?shù)過多。模型采用高階ELM方法處理動(dòng)量方程中的對(duì)流項(xiàng)，采用迎風(fēng)格式或隱式二階TVD格式處理物質(zhì)輸運(yùn)方程中的對(duì)流項(xiàng)，既可提高計(jì)算的時(shí)間步長(zhǎng)，又可保證模型計(jì)算穩(wěn)定性和效率。采用GLS模型求解摻混系數(shù)，并提供了與通用的GOTM模塊的接口。關(guān)于SCHISM模型的其他詳細(xì)介紹，可見參考文獻(xiàn)[13]，模型在水動(dòng)力模型基礎(chǔ)上耦合了波浪、生態(tài)、水質(zhì)、風(fēng)暴潮等模塊，已經(jīng)用于全球的水動(dòng)力相關(guān)問題的研究與預(yù)報(bào)，如北海和波羅海[16]等。

1.2 計(jì)算區(qū)域與網(wǎng)格

磨刀門三維數(shù)值模型共有2個(gè)開邊界，上游開邊界位于百頃頭以上的外海大橋附近，給流量過程，下游至磨刀門口外-15 m等深線附近，給水位過程；邊界流量、水位及鹽度均從已建的珠江口二維模型中提取，穩(wěn)定的初始鹽度場(chǎng)通過循環(huán)運(yùn)行模型得出。模型水平方向采用非結(jié)構(gòu)三角形網(wǎng)格，網(wǎng)格共23 283個(gè)，節(jié)點(diǎn)共13 168個(gè)；模型垂向采用混合SZ坐標(biāo)，平均海平面45 m以下采用Z坐標(biāo)，分1層；45 m以上采用S坐標(biāo)，分11層，每層厚度隨地形變化。網(wǎng)格分辨率的大小直接影響著模型的計(jì)算效率和精度，故分辨率分區(qū)域給定，口門區(qū)域分辨率取為1 500 m，而為確保鹽度模擬精度，河道內(nèi)分辨率取在40 m以下。

1.3 模型計(jì)算參數(shù)

模型計(jì)算時(shí)間為2009年1月2日～1月31日，共30 d，時(shí)間步長(zhǎng)為150 s；模型全域給恒定曼寧系數(shù)，大小取0.015；模型為斜壓模型，采用冷啟動(dòng)，不考慮底床變形和風(fēng)應(yīng)力，采用GLS湍流混合模型求解摻混系數(shù)，采用隱式TVD2格式處理輸運(yùn)方程中的對(duì)流項(xiàng)，采用高階ELM方法處理動(dòng)量方程中的對(duì)流項(xiàng)，干濕網(wǎng)格判斷最小水深均為0.1 m。

圖1 磨刀門水深地形及測(cè)點(diǎn)位置圖Fig.1 Topography of the Modaomen estuary and the locations of measurement stations

2 模型驗(yàn)證

采用2009年1月的實(shí)測(cè)資料對(duì)建立的數(shù)值模型進(jìn)行驗(yàn)證，驗(yàn)證點(diǎn)位置分布見圖1，其中三灶站、燈籠山站、竹銀站為水位測(cè)站，M1 ～ M8點(diǎn)為流速和鹽度測(cè)點(diǎn)。為定量評(píng)估驗(yàn)證精度，引入均方根誤差(RMSE)，其表達(dá)方式如下

(1)

式中：Xmod為模型計(jì)算值；Xobs為實(shí)際測(cè)量值；N為觀測(cè)次數(shù)。鑒于篇幅所限，本文只列出了部分驗(yàn)證成果，如圖2～圖4所示，圓點(diǎn)為實(shí)測(cè)值，實(shí)線為模擬值。驗(yàn)證表明，計(jì)算值與實(shí)測(cè)值吻合較好，誤差均在允許范圍之內(nèi)?？傮w上，模型能較好地反演磨刀門水流及鹽度的變化過程，因此建立的三維水流鹽度數(shù)值模型可用于鹽水入侵機(jī)制的研究。

圖3 M2、M6測(cè)點(diǎn)流速驗(yàn)證Fig.3 Calibration of flow velocity at M2 and M6 station

圖4 M2、M6測(cè)點(diǎn)鹽度驗(yàn)證Fig.4 Calibration of salinity at M2 and M6 stations

3 結(jié)果分析

3.1 磨刀門鹽水入侵特征

圖5 M6點(diǎn)表、底層鹽度時(shí)間過程線Fig.5 Temporal variation of surface and bottom salinity at M6 station

珠江河口的潮汐類型為不規(guī)則半日潮，一天中出現(xiàn)2次高潮和2次低潮，潮差和潮歷時(shí)明顯不等。圖5給出了M6點(diǎn)表、底層的鹽度時(shí)間變化過程，為說明鹽度隨大小潮潮型的變化，同時(shí)給出了三灶站的潮位過程。由圖可知，M6點(diǎn)的底層鹽度變化在0～9 ppt之間，表層鹽度變化在0～6 ppt之間；且鹽水入侵過程有著明顯的半月周期變化規(guī)律，表底層鹽度在小潮后期驟然升高，并在中潮時(shí)達(dá)到峰值，之后隨著潮差增大，鹽度慢慢變小。

圖6和圖7分別給出了小潮、大潮期間磨刀門水道潮平均的表、底層的流場(chǎng)和鹽度場(chǎng)分布情況，小潮期間(圖6)，表層余流在整個(gè)磨刀門水道包括洪灣水道及鶴州水道都是向海運(yùn)動(dòng)，磨刀門水道流速有明顯的橫向差異，東側(cè)流速明顯較西側(cè)大，這是因?yàn)闁|側(cè)為深槽，是重要潮流通道，流速相對(duì)較大，而西側(cè)為淺灘，底摩阻相對(duì)較大，削弱了潮流動(dòng)力，流速相對(duì)較??；底層余流在磨刀門水道主要向陸運(yùn)動(dòng)，在西岸淺灘區(qū)域則向海運(yùn)動(dòng)，在洪灣水道上段向海運(yùn)動(dòng)，下段則向陸運(yùn)動(dòng)；鹽水沿磨刀門水道東側(cè)深槽入侵，底部10 ppt鹽水可入侵到掛定角以上4 km附近，表層10 ppt等鹽度線則處于口門大橫琴附近，表底層鹽度差異大，鹽水分層明顯。大潮期間(圖7)，表層余流在磨刀門口外海區(qū)、口內(nèi)河道和洪灣水道仍是向海運(yùn)動(dòng)，底層余流在磨刀門水道深槽區(qū)向陸運(yùn)動(dòng)，在淺灘則向海運(yùn)動(dòng)；與小潮相比，磨刀門水道內(nèi)鹽水入侵減弱，底部10 ppt鹽水回落到掛定角以上1.5 km處，鹽水分層較弱。

從潮周期平均的流場(chǎng)和鹽度場(chǎng)分布可以看出，磨刀門水道東側(cè)深槽是高鹽水入侵磨刀門的重要通道，小潮期間磨刀門水道攔門沙至掛定角段底部積聚較高的鹽水，隨著潮差增大，積聚的鹽水被輸送至上游。

圖6 小潮期間表層(左)及底層(右)流場(chǎng)、鹽度分布Fig.6Surface(left)andbottom(right)meancurrentandtidally-averagedsalinityatneaptide圖7 大潮期間表層(左)及底層(右)流場(chǎng)、鹽度分布Fig.7Surface(left)andbottom(right)meancurrentandtidally-averagedsalinityatspringtide

3.2 磨刀門鹽水入侵動(dòng)力機(jī)制分析

3.2.1 鹽淡水混合特征

為了進(jìn)一步分析磨刀門鹽水入侵的動(dòng)力機(jī)制，利用數(shù)值模型計(jì)算結(jié)果，繪制了漲急、落急時(shí)刻A-A縱斷面的瞬時(shí)鹽度場(chǎng)，縱斷面位置如圖1所示。選取分層系數(shù)(n)來描述不同時(shí)刻鹽度的分層特征，n為測(cè)點(diǎn)表、底層的鹽度差與測(cè)點(diǎn)垂線平均鹽度的比值，定義為

(2)

圖8-a為小潮期間漲急、落急時(shí)刻磨刀門縱斷面鹽度分布圖，圖上橫軸為與外海A點(diǎn)的距離，流速箭頭代表沿河道的流速(下同)。小潮期間，三灶站的漲潮潮差為0.93 m，落潮潮差為0.68 m，潮汐動(dòng)力較小，縱斷面的鹽度等值線走勢(shì)比較平緩，沿程鹽度分層明顯。漲急時(shí)刻，口外底部水體鹽度變大，最高鹽度可達(dá)30 ppt，而表層鹽度值較低，隨著漲潮動(dòng)力增強(qiáng)，口外高鹽水團(tuán)沿河道底部向上游推進(jìn)，表層水體受其擠壓作用，形成明顯的高、低鹽水分界面；此時(shí)，n2= 0.77，n4= 1.65，n6= 0.01，磨刀門水道的鹽水幾乎處于高度分層狀態(tài)，且越往上游，分層越明顯，直到上游燈籠山附近，鹽水濃度很低，鹽度均勻分布，分層現(xiàn)象消失。落急時(shí)刻，中層及表層鹽水隨落潮流退至外海，但口門至洪灣水道段底部仍積蓄著高鹽水；此時(shí)，n2= 0.47，n4= 0.91，n6= 1.28，水道下游的鹽水處于緩混合狀態(tài)，而洪灣水道至燈籠山的鹽水仍處于高度分層狀態(tài)。

中潮期間(圖8-b)，三灶站的漲潮潮差為1.01 m，落潮潮差為1.73 m，與小潮期間相比，縱斷面的鹽度等值線趨勢(shì)有所傾斜，鹽水分層減弱。漲急時(shí)刻，口外鹽水楔進(jìn)入河道，與表層低鹽水混合后整體向上游入侵，水道內(nèi)鹽水濃度明顯變高，鹽水入侵距離增加；此時(shí)，n2= 0.51，n4= 0.98，n6= 1.46，在鹽水楔鋒面所在范圍內(nèi)，鹽水處于高度分層狀態(tài)。落急時(shí)刻，口外高鹽水稍有退卻，鹽水入侵距離達(dá)到峰值；此時(shí)，n2= 0.34，n4= 0.55，n6= 0.35，鹽淡水混合充分，河道內(nèi)的鹽水基本處于緩混合狀態(tài)。

大潮期間(圖8-c)，三灶站的漲潮潮差為1.31 m，落潮潮差為2.01 m。漲急時(shí)刻，口外高鹽水團(tuán)進(jìn)入河道，相比中潮漲急時(shí)刻，磨刀門水道內(nèi)鹽水濃度變低，鹽水入侵距離減??；此時(shí)，n2= 0.62，n4= 1.24，n6= 0.70，磨刀門水道的鹽水分層加強(qiáng)。落急時(shí)刻，口外高鹽水隨落潮流退出河道，此時(shí)，n2= 0.41，n4= 0.89，n6= 0.67，鹽水基本處于緩混合狀態(tài)，大潮潮汐動(dòng)力強(qiáng)，一天中鹽水濃度變化劇烈，鹽水“大進(jìn)大出”。

3.2.2 鹽通量分解

本文采用Lerczak等[4]提出的計(jì)算鹽度通量的方法，將河口的鹽度輸運(yùn)分為平流輸運(yùn)(advection transport)、穩(wěn)定剪切輸運(yùn)(steady shear transport)、潮汐震蕩輸運(yùn)(tidal oscillatory transport)三部分。斷面總鹽通量可以表示如下

FS==<(u0+uE+uT)(S0+SE+ST)dA>

(3)

式中：u為軸向流速，m/s；s為鹽度，ppt；A為橫斷面面積，m2；下標(biāo)0為潮平均的余流；下標(biāo)E為剪切流；下標(biāo)T為潮流波動(dòng)。MacCready[17]認(rèn)為公式中的交叉項(xiàng)可以近似為零，得出鹽度通量的分解表達(dá)式如下

FS=<(u0s0uESE+uTST)dA>=Qfs0+FE+FT

(4)

式中：FS為總鹽通量；Qfs0為平流輸運(yùn)項(xiàng)；FE為穩(wěn)定剪切項(xiàng)；FT為潮汐震蕩項(xiàng)。平流輸運(yùn)項(xiàng)是由上游流量以及外海stokes漂流引起的，主要是將上游鹽水輸運(yùn)至外海；穩(wěn)定剪切項(xiàng)主要由垂向剪切流引起，而潮汐震蕩項(xiàng)是由流速與鹽度的潮內(nèi)變化相互作用產(chǎn)生，它們則主要是將外海的高鹽水輸運(yùn)至河口上游。

以口門處的M2點(diǎn)為起始點(diǎn)，向上游河道每隔5 km取一個(gè)橫斷面，共選取6個(gè)橫斷面，計(jì)算并分析各斷面的鹽通量在一個(gè)潮周期過程中的變化情況。圖9給出了斷面2及斷面5的鹽通量分量Qfs0、FE和FT以及總鹽通量FS在一個(gè)完整潮周期過程中的變化情況，圖中正值表示向海輸運(yùn)；負(fù)值表示向陸輸運(yùn)，斷面具體位置如圖8所示。

圖9-c為斷面2的鹽通量分量，在斷面2處，Qfs0始終為正值，即平流輸送作用向海輸送鹽水，且在小潮期間平流輸送作用弱，在大潮期間平流輸送作用強(qiáng)；FE始終為負(fù)值，即穩(wěn)定剪切作用向陸輸送鹽水，且在小潮期間穩(wěn)定剪切作用強(qiáng)，在大潮期間穩(wěn)定剪切作用明顯減弱；潮汐震蕩作用FT在斷面2處向陸輸送鹽水，且在大潮期間最強(qiáng)，在小潮期間最弱；總鹽通量FS在小潮期間通過斷面2向河道內(nèi)輸送鹽水，而在大潮期間向外海輸送鹽水。

圖9-d為斷面5的鹽通量分量，在斷面5處，穩(wěn)定剪切輸送、平流輸送和潮汐震蕩的規(guī)律與斷面2基本相同，但輸送強(qiáng)度變??；FS在后半月周期的小潮期間向陸輸送鹽水，隨后在中潮轉(zhuǎn)大潮、大潮期間向海輸送鹽水，而在前半月周期(1月6日～1月17日)，F(xiàn)S基本為0，表明在Qfs0、FE和FT三者作用下，此期間磨刀門水道內(nèi)的日潮平均總鹽量處于平衡狀態(tài)。

9-a 三灶站水位 9-b 鹽水入侵長(zhǎng)度

9-c 斷面2的鹽通量分量 9-d 斷面5的鹽通量分量圖9 鹽通量分量過程圖Fig.9 Decomposition of salt transport flux

4 結(jié)論

基于SCHISM模型，對(duì)磨刀門水道鹽水入侵特性、縱向鹽度輸運(yùn)進(jìn)行了分析，得出以下結(jié)論：

(1)枯季磨刀門鹽水入侵一般特性為：小潮期間，磨刀門水道鹽水分層明顯，攔門沙至掛定角段底部積聚高鹽水，鹽水入侵距離較??；中潮期間，鹽水混合充分，小潮時(shí)積聚的底部高鹽水逐漸摻混至表層，并被輸送至上游，鹽水入侵距離最遠(yuǎn)；大潮期間，一天中鹽水濃度變化劇烈，鹽水“大進(jìn)大出”，鹽水入侵距離比中潮期間有所減小，但大于小潮期間的入侵距離。

(2)磨刀門日平均鹽度受平流輸送和穩(wěn)定剪切輸送共同作用，鹽量輸運(yùn)在大部分時(shí)間內(nèi)處于不平衡狀態(tài)。小潮期間，總鹽通量向陸，穩(wěn)定剪切作用強(qiáng)，平流輸運(yùn)作用弱，鹽量在河道里累積；中潮期間，總鹽通量由向陸轉(zhuǎn)為向海，河道內(nèi)的日平均鹽度最大，鹽水入侵距離最遠(yuǎn)；大潮期間，總鹽度通量向海，驅(qū)動(dòng)鹽水向陸輸送的穩(wěn)定剪切和潮汐震蕩作用小于驅(qū)動(dòng)鹽水向海輸送的平流輸送作用，鹽量被沖出河口。

參考文獻(xiàn)：

[1]PRITCHARD D W. Salinity distribution and circulation in the Chesapeake Bay estuarine system [J]. Journal of Marine Research, 1952, 11(2): 106-123.

[2]HANSEN D V, RATTRAY M. Gravitational circulation in straits and estuaries [J].J.mar.res,1965, 11(3): 319-326.

[3]SIMMONS H B, BROWN F R. Salinity effects on estuarine hydraulics and sedimentation [C]//Tanaka H. Process of 13th IAHR (3).Kyoto, Japan: International Association for HydroEnvironment Engineering and Research, 1969:192-216.

[4]LERCZAK J A, GEYER W R, CHANT R J. Mechanisms Driving the Time-Dependent Salt Flux in a Partially Stratified Estuary [J]. Journal of Physical Oceanography, 2006, 36(12): 2 296-2 311.

[5]CHEN W, CHEN K, KUANG C, et al. Influence of sea level rise on saline water intrusion in the Yangtze River Estuary, China [J]. Applied Ocean Research, 2016, 54:12-25.

[6]歐素英. 珠江三角洲咸潮活動(dòng)的空間差異性分析 [J]. 地理科學(xué), 2009, 29(1): 89-92.

OU S Y. Spatial difference about activity of saline water intrusion in the Pearl River Delta [J]. Scientia Geographica Sinica, 2009,29(1): 89-92.

[7]賈良文, 吳超羽, 任杰, 等. 珠江口磨刀門枯季水文特征及河口動(dòng)力過程 [J]. 水科學(xué)進(jìn)展, 2006, 17(1): 82-88.

JIA L W, WU C Y, REN J, et al. Hydrologic characteristics and estuarine dynamic process during the dry season in Modaomen Estuary of the Pearl River [J]. Advances in water science, 2006, 17(1): 82-88.

[8]韓志遠(yuǎn), 田向平, 劉峰. 珠江磨刀門水道咸潮上溯加劇的原因 [J]. 海洋學(xué)研究, 2010, 28(2): 52-59.

HAN Z Y, TIAN X P, LIU F. Study on the causes of intensified saline water intrusion into Modaomen estuary of Pearl River in recent years [J]. Journal of Marine Science, 2010, 28(2): 52-59.

[9]GONG W, SHEN J. The response of salt intrusion to changes in river discharge and tidal mixing during the dry season in the Modaomen Estuary, China [J]. Continental Shelf Research, 2011, 31(7): 769-788.

[10]盧陳, 袁麗蓉, 高時(shí)友, 等. 潮汐強(qiáng)度與咸潮上溯距離試驗(yàn) [J]. 水科學(xué)進(jìn)展, 2013, 24(2): 251-257.

LU C,YUAN L R,GAO S Y, et al. Experimental study on the relationship between tide strength and salt intrusion length [J]. Advances in water science, 2013, 24(2): 251-257.

[11]WANG B, ZHU J, WU H, et al. Dynamics of saltwater intrusion in the Modaomen Waterway of the Pearl River Estuary [J]. Science China Earth Sciences, 2012, 55(11): 1 901-1 918.

[12]陳文龍, 鄒華志, 董延軍. 磨刀門水道咸潮上溯動(dòng)力特性分析 [J]. 水科學(xué)進(jìn)展, 2014, 25(5): 713-723.

CHEN W L, ZOU H Z, DONG Y J. Hydrodynamic of saltwater intrusion in the Modaomen waterway [J]. Advances in water science, 2014, 25(5): 713-723.

[13]ZHANG Y J, YE F, STANEV E V, et al. Seamless cross-scale modeling with SCHISM [J]. Ocean Modelling, 2016, 102:64-81.

[14]ZHANG Y, BAPTISTA A M. SELFE: a semi-implicit Eulerian-Lagrangian finite-element model for cross-scale ocean circulation [J]. Ocean modelling, 2008, 21(3): 71-96.

[15]ZHANG Y J, ATELJEVICH E, YU H C, et al. A new vertical coordinate system for a 3D unstructured-grid model [J]. Ocean Modelling, 2015, 85:16-31.

[16]ZHANG Y J, STANEV E V, GRASHORN S. Unstructured-grid model for the North Sea and Baltic Sea: validation against observations [J]. Ocean Modelling, 2016, 97: 91-108.

[17]MACCREADY P, GEYER W R. Estuarine salt flux through an isohaline surface [J]. Journal of Geophysical Research Atmospheres, 2001, 106(C6): 11 629-11 638.