波浪-海流-微地形耦合的沉積動力模式建立及應用

蘆靜，夏長水，滕涌，劉學海

(1.中國海洋大學 海洋與大氣學院，山東 青島 266100；2.國家海洋局第一海洋研究所 海洋環(huán)境科學和數(shù)值模擬國家海洋局重點實驗室，山東 青島 266061)

?

波浪-海流-微地形耦合的沉積動力模式建立及應用

蘆靜1，2，夏長水2，滕涌2，劉學海2

(1.中國海洋大學 海洋與大氣學院，山東 青島 266100；2.國家海洋局第一海洋研究所 海洋環(huán)境科學和數(shù)值模擬國家海洋局重點實驗室，山東 青島 266061)

沙紋微地形普遍存在于海底，沙紋的消長能改變底部應力進而影響泥沙的運移。以往研究較多側重于波致沙紋，并已應用于波浪模式的底摩擦計算，而較少考慮波流聯(lián)合效應產生的沙紋，也未將其應用于綜合的水動力模式和沉積物輸運模式。本文在POM水動力模式中嵌入新南威爾士大學泥沙模式，通過耦合波流共同作用的微地形模型與波流相互作用底邊界層模型，發(fā)展了波浪-海流-微地形(沙紋)耦合的沉積動力模式。本文將該模式應用于澳大利亞Jervis灣，針對波主導和波流聯(lián)合主導沙紋兩種類型，分別進行了沙紋發(fā)展狀態(tài)、幾何形態(tài)的分布及懸浮泥沙的模擬。結果表明：波致沙紋比波流聯(lián)合作用的沙波具有更大的波高和波長，因此當波主導時沙紋對懸浮泥沙起著關鍵作用。通過考慮隨沙紋變化的粗糙度，相比于以往模式設置均一的粗糙度，該模型能對懸浮物濃度的驟升過程進行更精細的預測。

沉積動力模式；波流耦合；微地形；沙紋；底邊界層；底粗糙度；懸浮泥沙

1 引言

沙紋微地形為在沉積層表面形成的類似于水波紋的微細形態(tài)，其波長通常在幾厘米到幾十厘米。沙紋普遍存在于近岸非黏性底質的海底，偶爾也會在黏性底質出現(xiàn)。底應力作用是沙紋變化的主要動力因素，底應力的產生主要為波浪和底層流兩部分。按照沙紋主導機制，可分為波主導沙紋、流主導沙紋和波流聯(lián)合主導沙紋。

沙紋微地形可增加海底湍動能，因此對泥沙輸運具有重要的意義。沙紋出現(xiàn)能顯著增大底粗糙度，底部應力又是水動力和海底粗糙度共同作用的結果，因此沙紋的存在可使底部應力驟升，而底應力是決定底邊界層中泥沙再懸浮和沉積的關鍵。當海底為平底時，粗糙度由海底沉積物的粒度特征決定。一旦沙紋產生，將使底粗糙度顯著增大。Babanin等[1]對澳大利亞Lake George粉砂的實驗室研究表明沙紋的出現(xiàn)可使Nikuradse粗糙度Kb增大60倍?？蓜由臣y底形發(fā)展過程本身屬于推移質運動的一種形式，但它同時又影響懸移質運動。由于推移質輸沙量與沙紋運移速率和沙紋的波高成正比，沙紋可作為推移質運動的指征用來評估推移質運移的方向和大小[2]。懸浮泥沙運移受底部混合的影響極大，沙紋可以通過擴散和對流渦旋過程影響懸浮泥沙的夾帶機制，其計算依賴于不同形態(tài)函數(shù)類型以及參考濃度[3—5]。

波浪和海流均能產生沙紋，但一直以來波浪主導的沙紋得到研究重視[1,6-8]，而海流主導的沙紋研究較少。波、流聯(lián)合主導沙紋則限于復雜的相互作用機制，相關的研究更是少之又少，但在綜合的水動力模型中只考慮波浪產生的沙紋是不足的。Li和Amos[9-10]將沙紋引入一維沉積物輸運模型中，提出一個計算波流共同作用下的沙紋的算法，但其對沙紋增強的剪切速度的計算尚不完善。

雖然沙紋的存在非常普遍，但以往很少將其影響引入到沉積動力模式中，其對懸沙及推移質運移影響的研究更少。Tolman[11]給出了沙紋引起的波浪耗散的參數(shù)化方案并引入到WAVEWATCH波浪模式中, Ardhuin等[12]對該方案進行了修正。Smith等[8]將波致沙紋引起的底摩擦引入SWAN模式，提高了近岸的波浪模擬精度。在當前的沉積模式中，往往采用均一的海底粗糙度設置，既不能體現(xiàn)物理要素在空間上的細微分布，也不能體現(xiàn)微地形隨時間的變化，從而影響了對底層和近底層的動力機制及各物理要素分布變化的精細模擬。

本文發(fā)展了建立波浪-海流-微地形耦合的沉積動力模式，該模式通過改進參數(shù)化方案，建立模擬波流聯(lián)合作用下的沙紋幾何形態(tài)的模型，并實現(xiàn)沙紋模型與波流相互作用底邊界層模型[13-15]的耦合。首先，我們參考并改進了Li和Amos[9-10]及Tolman[11]的參數(shù)化方案，并編寫了獨立模塊來模擬波流聯(lián)合作用下的沙紋幾何形態(tài)；第二，改進沙紋增強的剪切速度計算模型，考慮沙紋的時空變化及其對底層剪切速度的影響，代替通常使用的均一粗糙度設置；第三，體現(xiàn)沙紋對懸浮泥沙濃度的作用，建立的模式將能夠更精細的模擬懸沙運移狀況。本文還將該模式用于Jervis灣懸浮泥沙的模擬，給出該海域沙紋的幾何形態(tài)分布及沙紋影響下懸沙濃度隨底層流及波浪的演變規(guī)律。

2 耦合模式的建立

2.1 沉積物輸運模式

我們以普林斯頓大學海洋模式(POM)[16]作為水動力基礎，并將新南威爾士大學泥沙模式(UNSW-Sed)嵌入到POM模式中。UNSW-Sed由POMSED發(fā)展而來[17-19]，它具有以下特點：第一，通過引入通量Richardson數(shù)Rf，考慮了底邊界層化效應[17]。拖曳系數(shù)Cd為：

(1)

式中，A=5.5，是經驗常數(shù)；zr是參考高度；kbc是表觀粗糙度，它由波浪邊界層和底物理粗糙度聯(lián)合得到；κ是卡曼常數(shù)；Rf為通量Richardson數(shù)。根據(jù)Blumbery和Mellor[20]，通量Richardson數(shù)Rf可由二階湍封閉算法中的梯度Richardson數(shù)Ri計算：

Rf=0.725[Ri+0.186-(Ri2-

0.316Ri+0.034 6)1/2]，

(2)

第二，根據(jù)霧狀層動力學，考慮了懸浮泥沙對海水密度的貢獻及泥沙密度對水動力的影響，適于模擬異重流等情況。本文在UNSW-SED基礎上做了如下兩點改進：一是將微地形模型嵌入波流相互作用的底邊界層模型[13-15]，并考慮了波致底沙液化[21]對底床的侵蝕作用；二是改進了懸浮泥沙濃度輻射邊界條件。

2.2 波流共同作用的底邊界層模型

2.2.1 波流聯(lián)合作用下的最大底切應力

基于波流耦合的底邊界層模型[13-15]，波流共同作用下的最大底切應力可表示為：

(3)

式中，τc為流致底切應力，它是一個平均分量；τw為波致底切應力，它是擾動分量的最大值。不考慮平均應力分量τc對擾動分量τw的作用。u*cw為波流共同作用的剪切速度，波流共同作用下的底切應力τb是決定沉積物在底邊界層中再懸浮和沉積通量的關鍵因子。

2.2.2 波浪邊界層和底物理粗糙度聯(lián)合引起的表觀底粗糙度

采用Signell等[15]的方法引入波浪對底流速的作用，根據(jù)波浪摩擦/剪切速度u*w利用迭代過程來確定參考高度zr處的粗糙度。底邊界模型中，定義波浪邊界層和表觀粗糙度kb聯(lián)合引起的Nikuradse底粗糙度kbc，

(4)

2.3 沙紋微地形計算模型

將上文的波流耦合底邊界層模型和微地形模型結合起來，以考慮沙紋存在產生的底邊界效應。

2.3.1 波/流主導的沙紋類型

要計算沙紋的幾何形態(tài)，首先須判斷沙紋的主導類型。波浪剪切速度與穩(wěn)定流的剪切速度之比U*ws/U*cs可用來定義波流共同作用下不同種類的沙紋:當U*ws/U*cs<0.75時流主導；當U*ws/U*cs>1.25時，波浪起主導作用；當0.75

2.3.2 波流聯(lián)合作用下的沙紋發(fā)展狀態(tài)

沙紋的5個發(fā)展階段[9-10]由剪切速度與不同臨界值的比較界定：(1)無傳輸狀態(tài)，U*cweU*cr，其中U*cws為僅基于底沙粒徑的波流聯(lián)合剪切速度；(3)平衡狀態(tài)，U*cws≥U*cr且U*cwb

表1 主要變量的定義

由于沙紋僅為幾厘米到幾十厘米量級，在數(shù)值計算中通常難以直接分辨其形態(tài)，因此粗糙沙紋底床上的底摩擦通常難以模擬。本文采用參數(shù)化方案，以沙紋的波高和波長為描述沙紋幾何形態(tài)兩個指標，模擬它們隨時間和空間的變化。參數(shù)化方案參考Li和Amos[9-10]的實驗室和現(xiàn)場觀測，我們將該方案與波流耦合沉積動力學模型銜接，并進行了對計算沙紋增強的剪切速度做了改進(見2.3.4節(jié))。

沙紋波高和波長的計算公式如下：

(1)無傳輸狀態(tài)，U*cws

(2)弱傳輸狀態(tài)，U*cwsU*cr。

η=(19.6·U*cws/U*cr+20.9)D,λ=η/0.12，

(5)

式中，D為底沙的中值粒徑。

(3)平衡狀態(tài)，U*cws≥U*cr且U*cwb

a．當波主導，即當U*ws/U*cs≥1.25時，

η=27.14·U*cws/U*cr+16.36D,λ=η/0.12，

(6)

b．當波流聯(lián)合主導，即當U*ws/U*cs<1.25時，

η=(22.15·U*cws/U*cr+6.38)D,λ=η/0.15.

(7)

(4)破碎狀態(tài)，U*cwb≥U*bf時，

λ=535D,η/λ=0.15U*up-U*cwb/

U*up-U*bf.

(8)

(5)層移狀態(tài)，U*cwb≥U*up時，理論上無沙紋。此時模式中仍然設置最小值η=0.01 m,λ=0.1 m，因此沙紋波陡η/λ=0.1。

2.3.3 沙紋增強的底粗糙度

Nikuradse粗糙度Kb由下式計算：

Kb=27.7η2/λr+170D(θs-0.05)0.5+2.5D.

(9)

Kb分為3部分：(1)底沙粒徑產生的粗糙度為2.5D，D為底沙中值粒徑；(2)170D(θs-0.05)0.5為推移質運動引起的粗糙度；(3)沙紋產生的粗糙度為27.7η2/λr。θs為Sheilds參數(shù)，由波流聯(lián)合作用下的底切應力決定。

2.3.4 沙紋增強的剪切速度及其改進

以往常用Nielsen[22]的結果計算沙紋增強的剪切速度U*cwe，但其僅基于粒徑的剪切速度U*cws和沙紋波陡，并只是簡單線性關系的表達式U*cwe_NL=U*cws/(1-πη/λ)，而U*cws由不考慮沙紋的且均一的Kb計算得到。而本文采用改進的沙紋增強的波流聯(lián)合剪切速度U*cwe_Kb，它由沙紋增強的Nikuradse粗糙度Kb(公式(9))導出。

我們改進的沙紋增強剪切速度U*cwe_Kb所使的沙紋增強粗糙度經波浪修正得到。通過考慮波流相互作用修正粗糙度Kb(波)以獲得Kbc(流)，Kb和Kbc用于計算波致底切應力的波浪摩擦因子fw及拖曳系數(shù)Cd，然后分別計算出流致剪切速度U*c和波致剪切速度U*w，再由公式(3)疊加得到U*cwe_Kb。也就是說，一方面,我們通過Kb引入了沙紋的效應，而在另一方面，我們通過修正Kb(沙紋增強的粗糙度)來引入表面波的效應。因此，這種方法比使用Nielsen[22]的表達式更加符合物理過程。

3 Jervis灣沉積模型的建立

3.1 模式區(qū)域地理位置、概況及網格分布

Jervis灣位于澳大利亞東南海岸、新南威爾士州南部，面積136 km2，是一個半封閉港[23]。Jervis灣南部陸地岬角設有海軍基地。Jervis灣的平均深度為15 m，北半部的淺水陸架廣闊，通過寬3.75 km的灣口與外陸架相連。灣內潮流較弱且主要由M2分潮控制[23]。灣口潮流流速約為0.07 m/s，而在灣北部潮流流速則小于0.01 m/s[24]。

Jervis灣沉積物以粗沉積物(CSM)為主，CSM的平均含量約占77%；但在Jervis灣北部的沉積物中，細顆粒泥沙物質(FSM,粒徑介于2～20 μm之間)的組分含量較高(達47%)。在Jervis灣也可以觀測到最細的黏土和淤泥沉積物(小于2 μm)，但細顆粒泥沙在Jervis灣中并不占主導地位(組分含量小于7%)[25]。

模式區(qū)域為31.45°～37.95°S，150.0°～156.95°E(圖1b)。網格設置為正交曲線網格系統(tǒng)，并在Jervis灣局部加密(圖1b)，網格數(shù)為234×258。x方向的水平分辨率為438～7 403 m，y方向的水平分辨率為556～6 922 m。垂向分為21個均勻的sigma層。

3.2 邊界條件

關于潮汐邊界的設置，我們考慮了4個主要分潮，即M2、S2、K1和O1。風場強迫和外部波浪邊界條件由ECMWF(Centre for Medium-Range Weather Forecasts)全球模式獲得。我們使用的ECMWF數(shù)據(jù)的分辨率為0.125°×0.125°，時間間隔為6 h。SWAN模式用10 m風速的U和V分量來驅動。ECMWF波浪資料被應用在3個開邊界上，這些開邊界條件由波浪和涌浪疊加的有效波高、平均波周期和平均波向組成。

3.3 泥沙參數(shù)及沙紋臨界值的選取

表2給出了泥沙輸運模式中系數(shù)的選取。根據(jù)澳大利亞地球科學院于2007、2008和2009年的調查結果，粒徑覆蓋范圍200～400 μm。在本文最關注的E站和F站，砂質組分所占比例超過85%，E站附近的平均粒徑大約介于300～350 μm之間，而F站附近的粒徑大于300 μm。因此，本文設置粒徑D為中間值340 μm，以確保應用Li和Amos[9-10]的方法和參數(shù)的合理性，他們的方法基于實驗室觀測和現(xiàn)場測量，也設置了0.34這個粒徑值。根據(jù)Miller等[26]的Yalin方法用來獲得臨界Shields參數(shù)θcr, 然后由此得出推移質運動的臨界剪切速度U*cr。懸移質啟動的臨界應力τcrs由Bagnold[27]得出。Jervis灣主要為非黏性底質，根據(jù)Gibbs等[28]的表達式，使用中值粒徑計算的沉降速度Ws具有10-3m/s的量級。動力學黏性系數(shù)μ的變化范圍根據(jù)不同水溫等條件變化范圍為0.001～0.001 2 N·s/m2。層移(推移質運動的最高級別)臨界Shields參數(shù)θ*up根據(jù)Li和Amos[29]的一個通用化的標準導出。

圖1 模型區(qū)和網格系統(tǒng)
Fig.1 Model domain and grid distribution
a為Jervis灣網格, E站(35.113°S, 150.708°E)位于Jervis灣西南水深11 m處，F(xiàn)站(35.075°S,150.765°E)位于灣口東北部，Tg站為潮位計的布放位置；b為模型區(qū)水深和曲線網格系統(tǒng)，其中Jervis灣局部加密
a.Jervis Bay grid Setting, site E(35.113°S, 150.708°E) located off the southwest coast in Jervis Bay with depth of 11 m, and Site F
(35.075°S, 150.765°E) northeast of the bay entrance, the diamond (Tg) denotes the location of the tidal gauge; b. bathymetry and curvilinear grid system with a locally refined grids in Jervis Bay

表2 泥沙模塊中的主要參數(shù)

3.4 數(shù)值實驗設置

我們進行了不同的數(shù)值實驗(表3)，包括考慮或不考慮沙紋的比較，并比較了使用U*cwe_Kb或U*cwe_NL之間的沙紋狀態(tài)及沙紋引起的懸浮泥沙濃度(suspended sediment concentration,SSC)差異。對于沒有沙紋的實驗，Kb設定為均一的0.05 m，這相當于有或者沒有沙紋的平均值。

表3 沙紋數(shù)值實驗

3.5 觀測儀器布放及數(shù)據(jù)分析方法

水位的觀測在Tg站(圖1a)使用Manly液壓實驗室(MHL)開發(fā)的電磁測潮儀(EG)，間隔15 min采集水位數(shù)據(jù)。有效波高的驗證選擇了的分別位于Batemans灣和悉尼外海兩個波浪騎士浮標。浮標施放地點位置見圖1b由白色菱形標注。Batemans灣和悉尼外海浮標觀測點的水深分別是55.7 m和93.9 m。波浪騎士浮標適于在較深或中等水深處投放，浮標施放位置的水深及鄰近條件可以使得波浪折射、衍射及淺水效應和摩擦衰減的影響微小而忽略。我們在Jervis灣布放了ADP實施高頻流速的測量，坐底式ADP每小時觀測1次，每次持續(xù)流速采樣時間間隔為1 s。布放于2009年11月16日在E站(35.113°S，150.708°E，水深11 m)進行，并回收于2009年12月18日。E站位于Jervis灣西南，站點E的位置見圖1a白色圓圈。

由信噪比(SNR)分析獲得SSC作為一種替代性的測量已被廣泛應用[30-31]。校準方程通常遵循一種對數(shù)的形式。在本文中，我們使用Xavier等[31]的表達式，實驗室結果顯示相關系數(shù)達0.99。

4 模擬結果

4.1 水動力模式驗證

4.1.1 大風過程中的浮標觀測的有效波高驗證

強風過程可以用來檢驗測試模型預測風浪的能力，因此我們選擇一個風速大于10 m/s的時間段用來驗證波浪模式，涵蓋從澳大利亞東部標準時間2009年10月10日10:00之后的457 h(圖2)。圖2表示風速大于10 m/s的條件下，在Batemans灣和悉尼外海的有效波高浮標觀測與模擬的比較。

圖2 2009年10月Batemans灣(a)和悉尼外海(b)近岸有效波高浮標觀測與模擬比較Fig.2 Validation rom offshore buoys in October 2009 in Batemans Bay (a) and Sydney (b)

圖2的波浪模擬結果顯示，在波浪近岸復雜因素可以忽略的情況下，SWAN模擬結果與實測符合較好，且與風速的發(fā)展趨向及數(shù)量變化基本符合。在此大風過程中，有效波高的決定因素基本來自本地風場。

4.1.2 水位的驗證

圖3給出了Jervis灣西南Tg站2008年和2011年兩次水位模擬結果與基于電磁驗潮計的資料的對比。

我們使用了潮汐站Tg(圖1a)的水位來驗證模式的各分潮(圖3)。模式與觀測結果的均方根(RMS)誤差計算公式如下：

(10)

式中，A和φ分別是潮汐振幅和位相；下標m和o分別表示模擬和觀測。表4給出了2008年和2011年用于模型強迫的4個主要分潮的觀測和模擬的結果的均方根誤差統(tǒng)計。

表4 2008年和2011年各分潮觀測和模擬結果之間的均方根誤差對比

圖3 2008年(a)和2011年(b)Tg站的觀測和模型潮位對比Fig.3 Comparison of tidal level from observation and model results at the Site Tg for 2008 (a) and 2011 (b)

圖4 模擬的2009年11月Jervis灣局部沙紋類型(a)和發(fā)展階段(b)Fig.4 Distributions of ripple types, ranges and geometry with Jervis Bay locally amplifieda.流/波主導類型：色標值-1代表當U*ws/U*cs<0.75時流主導, 1代表當U*ws/U*cs>1.25時波主導, 0代表當0.751.25, and 0 represents combined flow range while0.75

從圖3及表4來看，模式計算的各分潮振幅和相位的時間序列及誤差均達到使用要求。

4.2 沙紋的水平分布

4.2.1 沙紋的波流主導類型以及發(fā)展階段

不同的海底動力條件下的沙紋形態(tài)各異。波致海底表面摩擦(剪切)速度與穩(wěn)定流產生的摩擦(剪切)速度的比值，U*ws/U*cs被用來計算波流聯(lián)合流作用下的沙紋類型(浪/流主導)(2.3.1節(jié))。圖4和圖5的平均分布，是指本文模擬的2009年11月期間的平均狀態(tài)，這期間涵蓋波浪增大和減小的狀態(tài)。我們將11月每小時的輸出結果進行月平均。

圖4給出了模式模擬的Jervis灣內各種沙紋的分布及變化。針對沙紋類型，模式輸出-1代表流主導，輸出0代表波流共同主導，輸出1代表波主導。平均得到的值就會位于-1到1這個區(qū)間上，由此得到圖4a。平均值越接近1，波浪的成分就越大，越接近-1，流的成分就越大。針對沙紋發(fā)展階段，模式輸出0代表無傳輸階段，模式輸出1代表弱傳輸階段，模式輸出2代表平衡狀態(tài)，模式輸出3代表破碎狀態(tài)，模式輸出4代表層移狀態(tài)。平均值位于0到4這個區(qū)間上，由此得到圖4b。平均值由0增加到4的過程中，水動力作用逐漸增強，沙紋開始出現(xiàn)，波高波長由小變大，最后沙紋發(fā)生破碎，沙紋消失。

由圖4a可見，波浪在水深小于10 m的近岸海區(qū)起主導作用(U*ws/U*cs>1.25)，而流主導(U*ws/U*cs<0.75)的區(qū)域在Jervis灣中部水深大于20 m的海區(qū)，波流聯(lián)合主導0.75

4.2.2 模擬沙紋幾何形態(tài)

圖5a和5b分別為11月每小時的輸出的沙紋波高和波長的月平均。圖5給出了模擬的Jervis灣內沙紋幾何形態(tài)(即平均波高和波長)的水平分布及變化。

圖5 模擬的Jervis灣局部沙紋幾何形態(tài)的平均分布Fig.5 Averaged distributions of ripple geometry with Jervis Bay locally amplified

圖5a和5b顯示出沙紋波高和波長的分布及大小與圖4a和4b中的沙紋類型和狀態(tài)有相似的分布趨勢。沙紋波高η(10～14 cm)和波長λ(9～12 cm)通常在波主導的區(qū)域較高；此外，灣口也呈現(xiàn)出較大的沙紋波高和波長，其原因可能是潮流較大。這意味著，波浪主導的沙紋的波高和波長通常高于流主導的沙紋波高和波長。在Jervis灣外側東北部的近岸區(qū)域，沙紋從破碎狀態(tài)發(fā)展到層移狀態(tài)，沙紋波高和波長在平衡狀態(tài)達到最高，從破碎狀態(tài)發(fā)展到層移狀態(tài)的過程中逐漸減小。

4.3 沙紋隨底層流及波浪的演變規(guī)律

4.3.1 懸沙濃度驗證以及沙紋對其作用

E站位于東南沿岸的外海11.0 m水深處，ADP于2009年11月在此布放。E站大多數(shù)時候的沙紋類型為波浪主導(圖4a和圖7b)。F站選在灣口的東北處，這里的潮流較強。因此F站具有典型的流主導的沙紋特征(圖4a和圖8b)。F站沒有ADP流觀測，選取F站的原因是為了確定我們在波占主導地位的區(qū)域得出的結論是否只適用于波主導的沙紋還是同樣適用于流占主導的沙紋區(qū)域。

圖6顯示了2009年11月E站和F站處沙紋相關的各項模擬結果。如圖6所示，考慮沙紋的SSC變化更為劇烈與也與觀測符合更好。從圖6a和6b來看，考慮沙紋的SSC都比不考慮沙紋時變化更為劇烈。考慮沙紋的情況下，底切應力(bottom shear stress, BSS)能達到5.0 N/cm2，SSC達到60 mg/L，而在無沙紋(設置平均Kb=0.05 m)的實驗中，模擬的最大BSS僅為2.0 N/cm2，SSC為20 mg/L。然而，在BSS低于0.8 N/cm2的部分，考慮沙紋情況下BSS和SSC低于不考慮沙紋情況。

圖6 2009年11月26日00:00至28日17:00 E站(a)和F站(b)的SSC以及E站(c)和F站(d)的底層流速和波軌道速度，其中E站有ADP觀測Fig.6 a and b are SSC from 00:00 on 26th to 17:00 on 28th November, 2009, c and d are bottom tidal velocity and wave orbital velocity, a and c are for Site E with ADP measurements

圖7 E站(a)和F站(b)分別使用U*cwe_Kb和U*cwe_NL計算的沙紋狀態(tài)時間序列Fig.7 Time series of ripple stages at Site E (a) and Site F (b) under U*cwe_Kb and U*cwe_NL respectively最上方框為波主導，中間灰色方框代表波流聯(lián)合主導，最下層方框為主導；3種不同波或流主導的狀態(tài)下又分別包括無傳輸狀態(tài)、弱傳輸狀態(tài)、平衡狀態(tài)、沙紋的破碎狀態(tài)和層移狀態(tài)The top box represents wave dominant, the middle box represents combined flow dominant, the bottom box represent current dominant. These 3 wave/current-dominant stages include no transport ranges, weak transport ranges, equilibrium ranges, break off ranges and sheet flow statuses

圖8 E站(a)和F站(b)沙紋幾何形態(tài)的時間序列及使用U*cwe_Kb和U*cwe_NL的比較，虛線中間灰色區(qū)域與ADP觀測時間對應Fig.8 Time series at Site E(a)and Site F(b)ripple dimensions and their comparison between using U*cwe_Kb and U*cwe_NL，time period between dashed lines corresponds to measurement time

圖9 2009年11月E站(a，c)和F站(b，d)拖曳系數(shù)Cd(a，b)和波浪摩擦因子fw(c，d)時間序列Fig.9 Time series of drag coefficient Cd(a，b)and fw(c，d) in November 2009 at Site E(a，c) and Site F(b，d)

在波通常占主導地位E站(圖4a)，使用U*cwe_NL計算的SSC與觀測相比明顯被低估了。波浪在近岸區(qū)域起著關鍵作用，因此該低估不能忽視。再懸浮的臨界底切應力為1.2 N/m2(剪切速度0.034 9 m/s)，推移質運動的臨界剪切應力為0.16 N/m2(剪切速度0.012 7 m/s) (表2)。如果我們人為設置再懸浮的臨界剪切應力，給定一個較低值，如0.50 N/m2, 使用U*cwe_NL時SSC的峰值也有會被模擬出來，然而0.50 N/m2(相當于剪切速度0.022 m/s)遠小于之前Li和Amos[9-10]的臨界值標準(1.2 N/m2)，這個臨界值針對0.34 mm粒徑由Shields參數(shù)計算得到。對于波浪并不主導的F站,使用U*cwe_Kb或者U*cwe_NL模擬的BSS在波浪較大時仍然呈現(xiàn)出較大的差異。

4.3.2 沙紋發(fā)展狀態(tài)及幾何形態(tài)的演變規(guī)律

圖7給出了E站和F站不同類型沙紋(波/流主導)的發(fā)展狀態(tài)。由圖7可見，U*cwe_Kb和U*cwe_NL之間的差異主要出現(xiàn)在從“無傳輸狀態(tài)”到“弱傳輸狀態(tài)”的過渡過程。當考慮沙紋時，低估U*cwe導致了底床停留在無傳輸狀態(tài)，而不是進展到弱傳輸狀態(tài)。同時，圖7c也說明E站沙紋是主要由波主導的。平衡態(tài)沙紋相對出現(xiàn)較少，E站以無傳輸和弱傳輸狀態(tài)為主。

F站屬于波流聯(lián)合主導的海區(qū)(圖4a)。F站接近Jervis灣灣口，這里潮流流速是E站流速的2倍。在波浪較小時，F(xiàn)站的沙紋類型可以達到流主導的無傳輸狀態(tài)(圖7b)以及波流聯(lián)合占主導的弱傳輸狀態(tài)。

圖8給出了E站和F站沙紋幾何形態(tài)的時間序列。如圖8所示，流或波流聯(lián)合主導的沙紋(圖8b)比波浪主導的沙紋(圖8a)在沙紋波高η和波長λ的變化上更加平緩。沙紋波高和波長在沙紋產生后迅速上升，到達頂峰后逐漸減少。沙紋波高變化范圍為1～2 cm，而沙紋波長變化范圍為8～13 cm。沙紋波高η和波長λ都在平衡狀態(tài)達到最大值。使用U*cwe_Kb或U*cwe_NL，沙紋幾何形態(tài)(波高和波長)并未出現(xiàn)明顯變化。究其原因，U*cwe僅僅是用以通過和臨界值比較，決定底邊界狀態(tài)是否從無傳輸升級到弱傳輸狀態(tài)；而U*cws和U*cwb則通過和臨界值比較來判斷沙紋處于平衡、破碎還是層移狀態(tài)。

如圖9所示，使用Nielson(1986)方法計算的拖曳系數(shù)Cd_NL略大于使用我們Kb方法計算的Cd_Kb。拖曳系數(shù)Cd經過波浪修正后得到Cd_new。然而，在波浪存在時，fw_Kb明顯高于fw_NL。因此fw是造成U*cwe_Kb甚至比U*cwe_NL大的關鍵。U*cwe_Kb與U*cwe_NL的差別主要體現(xiàn)在波致剪切速度U*w上。

5 討論

5.1 隨沙紋變化的粗糙度

粗糙度是影響的底切應力BSS和SSC的關鍵因素。當不考慮沙紋時，以前的模式往往假設了均一的Nikurade粗糙度Kb或者Z0(Kb/30)。例如，在POM模式的水動力模擬中Z0被設置為0.01 m，以及在Ecomsed泥沙模式中Z0被設置為0.001 4 m作為默認值。這個默認的Kb通常是一個平均狀態(tài)，這個平均態(tài)比僅僅基于粒徑的Kb(gs)高，但往往比沙紋增強的Kb(gs+bl+rp)小。導致的結果是，平均的Kb不能體現(xiàn)沙紋發(fā)展的動態(tài)過程，因而不能更好地模擬SSC的變化。

5.2 波浪效應的數(shù)值穩(wěn)定性

5.3 平均或最大波浪軌道速度的選取

當考慮波浪效應對粗糙度的影響時，我們使用平均波浪軌道速度。因為當用波浪修正表粗糙度Kb而得到表觀粗糙度Kbc(Kbc為計算拖曳系數(shù)Cd所需)，是波浪的平均狀態(tài)影響波浪在波流共同作用中所占的比例(公式(4)，系數(shù)β)而不是最大狀態(tài)。類似地, Smith等[8]通過積分波浪譜來獲得平均近底軌道速度Ub并用于計算波致沙紋。

然而，當確定底部泥沙是否會發(fā)生再懸浮時我們使用波浪軌道速度最大值來計算波致底切應力，因為是波浪擾動分量的最大值決定了底沙是否能夠進入水體[14]。這時，流致底切應力為平均量，而不是一個周期性變化的值。因此波流共同作用下的最大底切應力τb,max(公式(3))是平均流致切應力與最大的波浪致切應力的疊加。

6 結論

6.1 沙紋微地形對懸浮泥沙的影響

當沙紋出現(xiàn)時，波浪可以大幅提高懸浮泥沙濃度SSC?？紤]沙紋之后，SSC比不考慮沙紋時變化更為劇烈。因此，與以往采用均一的底粗糙度相比，考慮沙紋之后實際觀測的SSC驟升過程被更準確地模擬出來。此外，我們改進了計算沙紋增強的剪切速度的方案，記為U*cwe_Kb，當波浪主導時，使用我們提出的剪切速度U*cwe_Kb可明顯提高SSC的模擬量值，與觀測也更為符合。在無傳輸和弱傳輸狀態(tài)，使用U*cwe_NL與U*cwe_Kb,造成的SSC差異也不甚明顯，因為此時沙紋波高和波長都較小。

6.2 不同類型沙紋的分布

波浪主導區(qū)域分布在近岸水深小于10 m的區(qū)域，而流主導的區(qū)域覆蓋灣的中部以及Jervis灣灣口水深大于20 m的海區(qū)。Jervis灣波浪主導的沙紋波高(10～14 cm)和沙紋波長(9～12 cm)通常高于流主導的沙紋，且流主導的沙紋在沙紋波高和波長上比波浪主導的沙紋隨時間的變化更為平緩。使用U*cwe_NL與U*cwe_Kb在沙紋幾何形態(tài)上沒有明顯的差異。

無傳輸狀態(tài)主要出現(xiàn)在水動力條件較弱的Jervis灣中部，而弱傳輸狀態(tài)覆蓋近岸區(qū)域。平衡狀態(tài)主要出現(xiàn)在西南海岸的狹窄的帶狀海區(qū)以及灣口兩側的近岸海區(qū)，平衡狀態(tài)具有較大的沙紋波高和波長。Jervis灣少有破碎狀態(tài)而沒有層移狀態(tài)。U*cwe僅僅是用于和臨界值比較作為判別，以確定底邊界狀態(tài)是否從無傳輸?shù)饺鮽鬏敔顟B(tài)，因此低估U*cwe導致底部形態(tài)仍然停留在無傳輸狀態(tài)，而不是進一步發(fā)展到弱傳輸狀態(tài)，因此使用U*cwe_NL并無明顯不同。

6.3 展望

波狀海底微地形根據(jù)尺度不同，可分為沙丘等不同類型。本文研究的沙紋尺度較小，波長為幾厘米至十幾厘米。今后可推廣到大尺度沙丘或者沙脊。如南海波長為10 m量級的沙紋和中國蘇北淺灘著名的沙脊等。

目前沙紋底邊界層模型中底粗糙度、拖曳系數(shù)等的計算主要采用參數(shù)化方法，這些參數(shù)化方案雖然已經較為準確地反映出真實的運動過程，但對起物理機制的把握尚不完備，今后可進一步發(fā)展理論模型。

[1] Babanin A V, Young I R, Mirfenderesk H. Field and laboratory measurements of wave-bottom interaction[M]//Presented at the Coasts and Ports: Coastal Living-Living Coast. Adelaide: Institution of Engineers, 2005.

[2] Holmedal L E, Myrhaug D. Bed load transport under irregular waves plus current from Monte Carlo simulations of parameterized models with application to ripple migration rates observed in the field [J]. Coastal Engineering, 2004, 51(2): 155-72.

[3] Ribberink J S, Al-Salem A A. Sediment transport in oscillatory boundary layers in cases of rippled beds and sheet flow[J]. Journal of Geophysical Research Atmospheres, 1994, 1994(C6): 12707-12728.

[4] Werf J J V D, Ribberink J S, O'Donoghue T, et al. Modelling and measurement of sand transport processes over full-scale ripples in oscillatory flow[J]. Coastal Engineering, 2006, 53(8): 657-673.

[5] Thorne P D, Davies A G, Bell P S. Observations and analysis of sediment diffusivity profiles over sandy rippled beds under waves[J]. Journal of Geophysical Research Oceans, 2009, 114(C2): 309-321.

[6] Nielsen P. Dynamics & geometry of wave-generated ripples[J]. Journal of Geophysical Research, 1981, 86(C7): 6467-6472.

[7] Grant W D, Madsen S O. Movable bed roughness in unsteady oscillatroy flow[J]. Journal of Geophysical Research, 1982, 87(C1): 469-481.

[8] Smith G A, Babanin A V, Riedel P, et al. Introduction of a new friction routine into the SWAN model that evaluates roughness due to bedform and sediment size changes[J]. Coastal Engineering, 2011, 58(4): 317-326.

[9] Li M Z, Amos C L. Predicting ripple geometry and bed roughness under combined waves and currents in a continental shelf environment[J]. Continental Shelf Research, 1998, 18(9): 941-970.

[10] Li M Z, Amos C L. SEDTRANS96: the upgraded and better calibrated sediment-transport model for continental shelves[J]. Computers & Geosciences, 2001, 27(6): 619-645.

[11] Tolman H L. Subgrid modeling of moveable-bed bottom friction in wind-wave models[J]. Coastal Engineering, 1995, 26(1/2): 57-75.

[12] Ardhuin F, Drake T G, Herbers T H C. Observations of wave-generated vortex ripples on the north carolina continental shelf[J]. Journal of Geophysical Research, 2002, 107(10): 7-1-7-14.

[13] Grant W D, Madsen O S. Combined wave and current interaction with a rough bottom[J]. Journal of Geophysical Research, 1979, 84(C4): 1797-1808.

[14] Grant W D, Madsen O S. The continental-shelf bottom boundary layer[J]. Fluid Mechanics Annual Reviews, 1986, 18(1): 265-305.

[15] Signell R P, Beardsley R C, Graber H C, et al. Effect of wave-current interaction on steady wind-driven circulation in narrow, shallow embayments[J]. Journal of Geophisical Research, 1990, 95(C6): 9671-9678.

[16] Mellor G L, Donelan M A, Oey L Y. A surface wave model for coupling with numerical ocean circulation models[J]. Journal of Atmospheric & Oceanic Technology, 2008, 25(10): 1785-1807.

[17] Wang X H. Tide-induced sediment resuspension and the bottom boundary layer in an idealized estuary[J]. Journal of Physical Oceanography, 2002, 32(4): 3113-3131.

[18] Wang X H, Pinardi N, Malacic V. Sediment transport and resuspension due to combined motion of wave and current in the northern Adriatic Sea during a Bora event in January 2001: A numerical modeling study[J]. Continental Shelf Research, 2007, 27(5): 613-633.

[19] Song D H, Wang X H, Cao Z Y, et al. Suspended sediment transport in the Deepwater Navigation Channel, Yangtze River Estuary, China, in the dry season 2009: 2. Numerical simulations[J]. Journal of Geophysical Research, 2013, 118(10): 5568-5590.

[20] Blumberg, A F, Mellor, G L. A description of a three-dimensional coastal ocean circulation model[M]//Three-Dimensional Coastal Ocean Models. Washington, D.C.:American Geophysical Union, 1987.

[21] Lambrechts J, Humphrey C, McKinna L, et al. Importance of wave-induced bed liquefaction in the fine sediment budget of Cleveland Bay, Great Barrier Reef[J]. Estuarine, Coastal and Shelf Science, 2010, 89(2): 154-162.

[22] Nielsen P. Suspended sediment concentrations under waves[J]. Coastal Engineering, 1986, 10(1): 23-31.

[23] Wang X H. A numerical study of sediment transport in a coastal embayment during winter storms[J]. Journal of Coastal Research, 2001(34): 414-427.

[24] Holloway P E, Symonds G, Nunes V R. Observations of circulation and exchange processes in Jervis Bay, New South Wales[J]. Australian Journal of Marine and Freshwater Research, 1992, 43(6): 1487-515.

[25] CSIRO. Jervis Bay Baseline Studies Final Report[R]. CSIRO Division of Fisheries, Marmion Research Laboratories, 1994.

[26] Miller M C, McCave I N, Komar P D. Threshold of sediment motion under unidirectional currents[J]. Sedimentology, 1977, 24(4): 507-527.

[27] Bagnold R A. An approach to the sediment transport problem from general physics[R]. Washington:U. S. Govt. Print. Off.,1966.

[28] Gibbs R J, Matthews M D, Link D A. The relationship between sphere size and settling velocity[J]. Journal of Sedimentary Research, 1971, 41(1): 7-18.

[29] Li M Z, Amos C L. Sheet flow and large wave ripples under combined waves and currents: their field observation, model prediction and effects on boundary layer dynamics[J]. Continental Shelf Research, 1999, 19(5): 637-663.

[30] Salehi M, Strom K. Using velocimeter signal to noise ratio as a surrogate measure of suspended mud concentration[J]. Continental Shelf Research, 2011, 31(9): 1020-1032.

[31] Xavier B C, Silva I O, Guimar?es L G, et al. Estimation of suspended sediment concentration by acoustic scattering: an experimental and theoretical analysis for spherical particles[J]. Journal of Soils & Sediments, 2014, 14(7): 1325-1333.

Developing the wave-current-microtopography coupled model of sediment dynamics and its applications

Lu Jing1，2，Xia Changshui2，Teng Yong2，Liu Xuehai2

(1.CollegeofOceanicandAtmosphericSciences,OceanUniversityofChina,Qingdao266100,China;2.KeyLabofMarineScienceandNumericalModeling,TheFirstInstituteofOceanography,StateOceanicAdministration,Qingdao266061,China)

Ripple microtopography prevalently exist on coastal beds, which significantly change the bottom stress and then influence the sediment transport. Previous researchers mainly study the wave-dominant ripples and have applied them to the wave modelling. Wave-current combined flow is rarely discussed to generate ripples, and the combined flow-dominant ripples are rarely implemented to hydrodynamic and sediment models. We embedded the University of New South Wales sediment model into the POM model, and connected them with a wave-current interaction bottom boundary model coupled with a microtopography module under combined flow. We developed a wave-current-microtopography coupled model of sedimentology dynamics, and applied this coupled model to Jervis Bay, Australia. Developing stages and types are modeled, and the ripple height and length are simulated. The suspended sediment transport was analyzed under wave-dominant and combined flow separately. Simulated results show that the wave-dominant ripples have longer height and length. Therefore, ripples place an important role on suspended sediment when waves dominate. Through considering roughness that varies with microtopography, this model can predict the abrupt rising of suspended sediment concentration rather than setting an average uniform roughness.

sediment dynamical model;wave-current coupling; microtopography; ripples; bottom boundary layer; roughness; suspended sediments

2016-10-14；

2016-12-13。

國家自然科學基金項目(41072176)；中央級公益性科研院所基本科研業(yè)務費專項(2014T01,2015P03)；國家重點研發(fā)計劃課題(2016YFC0503602，2016YFB0201103，2017YFA0604101，2017YFA0604104)；留學回國人員科技活動項目擇優(yōu)資助。

蘆靜(1982—)，女，四川省成都市人，助理研究員，博士，主要從事物理海洋研究、沉積動力學數(shù)值模擬。E-mail:lujing@fio.org.cn

10.3969/j.issn.0253-4193.2017.07.002

P731.2

A

0253-4193(2017)07-0012-14

蘆靜, 夏長水, 滕涌, 等. 波浪-海流-微地形耦合的沉積動力模式建立及應用[J]. 海洋學報, 2017, 39(7): 12-25,

Lu Jing, Xia Changshui, Teng Yong, et al. Developing the wave-current-microtopography coupled model of sediment dynamics and its applications[J]. Haiyang Xuebao, 2017, 39(7): 12-25, doi:10.3969/j.issn.0253-4193.2017.07.002