基于Boussinesq型方程的錢塘江涌潮數(shù)值模擬

黃婷，張懷,2*，石耀霖

1 中國科學院大學地球與行星科學學院，計算地球動力學重點實驗室，北京 100049 2 南方海洋科學與工程廣東省實驗室 (珠海)，廣東珠海 519080

0 引言

涌潮是一種發(fā)生在漲潮期間的潮汐現(xiàn)象，來潮沿著河流或狹窄海灣逆流而上，形成波長和頻率不斷變換的潮波.涌潮的形成基于來潮、河流條件以及河道地形之間的微妙平衡，并非所有的河口都會發(fā)生涌潮.據(jù)估計，全世界發(fā)生涌潮的河口有四百多個，分布在除南極洲以外的所有大陸，河口附近的環(huán)境、生態(tài)以及文化都將受到涌潮所帶來的影響(Chanson,2009).杭州灣是著名的涌潮地，而自澉浦以上的錢塘江出口段(圖1)是世界聞名的強涌潮河口(譚維炎等,1995)，因錢塘江河床泥沙易沖易淤，為治理錢塘江河口環(huán)境，20 世紀60 年代以來，治江縮窄工程不斷推進(尤愛菊等,2010)，河流邊界內(nèi)移，使得錢塘江河口處形成形態(tài)復雜多樣的涌潮.最具代表性的潮景是大缺口的“交叉潮”、鹽官的“一線潮”以及老鹽倉的“回頭潮”，每年都會吸引無數(shù)中外游客慕名而來，是我國獨特的自然旅游資源.對錢塘江涌潮的研究，不僅具有學術意義，也具有重要的實際應用價值.

圖1 錢塘江河口地形示意圖黑色輪廓線是20世紀初錢塘江河流邊界，藍色線是江底地形等高線.Fig.1 Schematic diagram of the Qiantang EstuaryThe black outline indicates the boundary of the Qiantang River in the early 20th century.The blue line indicates the contour line of the topography of the river bottom.

幾個世紀以來，國內(nèi)外的學者一直致力研究涌潮現(xiàn)象.Saint-Venant(1871)對一維渠道內(nèi)的水流進行觀測，提出了著名的圣維南方程，并將其應用于法國塞納河的涌潮，通過計算來潮的Froude數(shù)，判斷涌潮的發(fā)生.一維模型基于理想假設，適用于流速方向易于評估，并且在整個涌潮持續(xù)時間內(nèi)固定的河道(Clamond and Dutykh,2013).趙雪華(1985)應用一維圣維南方程，有機地結合特征線法和特征差分格式，計算了錢塘江自杭州閘口至下游澉浦之間河段涌潮時的水位和流量情況.于普兵和潘存鴻(2010)應用迎風格式有限體積法求解一維圣維南方程，建立了錢塘江涌潮一維數(shù)值模型，并模擬計算從上游富春江電站至下游澉浦河段的潮位情況.

錢塘江涌潮存在形狀良好的波狀涌潮和潮頭破碎的強涌潮等多種形態(tài)，且波狀涌潮的頻散作用較大，需要在控制方程中添加頻散項(潘存鴻等,2008)，以往的模擬大多針對的是涌潮的傳播過程并不考慮頻散效應.綜上所述，本文在前人研究的基礎上建立了新的二維高階Boussinesq型方程，采用有限體積法，考慮移動的干濕邊界，對錢塘江涌潮進行數(shù)值模擬研究，探討了在考慮頻散作用下錢塘江涌潮的形成與傳播特征.

1 模型與數(shù)值方法

1.1 控制方程

1.1.1 Boussinesq型方程

Madsen和S?rensen(1992)從經(jīng)典Boussinesq 方程出發(fā)，提出了緩坡地形上的一種Boussinesq型方程,并通過理論分析和應用驗證該方程比經(jīng)典Boussinesq方程適用的水深范圍更廣.經(jīng)典Boussinesq方程適用于水深與波長比值h0/λ0<1/5的水深區(qū)域，而Madsen和S?rensen的方程在調(diào)節(jié)適當?shù)念l散系數(shù)后可將應用水深擴展到h0/λ0<1/2的區(qū)域.Beji和Nadaoka(1996)同樣從經(jīng)典Boussinesq方程出發(fā)，通過對非線性項和頻散項的數(shù)學關系代換，修正了Madsen和S?rensen的方程，提出了一種考慮地形變化的二維高階Boussinesq型方程.我們在Beji和Nadaoka工作的基礎上，導出在考慮河道彎曲變化和強烈地形起伏情形下不可壓縮流體的二維高階Boussinesq型方程，表示為

(1)

(2)

(3)

其中n為曼寧系數(shù)，F(xiàn)=F(fx,fy)為科氏力，有

fx=fv,fy=-fu,

(4)

其中f=2Ωsinφ為Coriolis參數(shù)，Ω為地球自轉角速度，φ為緯度值.A=A(Ax,Ay)表示水平動量擴散

(5)

其中Am為水平渦流黏度系數(shù)，可以簡單地采用渦流黏性去模擬湍流混合和波破碎時的能量耗散(Chen et al.,2000;Kennedy et al.,2000).B=(Bx,By)稱為Boussinesq項或是頻散項，是關于時間和空間的混合高階導數(shù)，是在靜水壓力平衡基礎上的一個附加項

(6)

其中h為靜水深，參數(shù)β調(diào)節(jié)著方程的頻散程度.

1.1.2 方程的線性頻散性

長波假設意味著，波的特征振幅η0與特征水深h0相比較小，即系數(shù)ε=η0/h0?1,波的傳播體現(xiàn)為波的相位傳播，且波的相速度與stokes一階理論(即Airy線性波理論)表述的一致，于是可以在長波假設下對控制方程進行線性頻散性分析.

線性頻散關系表示為：

(7)

式中ω為波的頻率，k=2π/λ為波數(shù)，h為靜水深.進而得到波速c與水深的線性頻散關系

(8)

將式(8)在kh=0處進行泰勒展開有

(9)

截取到O(k4h4)項，即方程的精度達到O(σ4)階.

如前所述的控制方程(1)的頻散關系表示為(Beji and Nadaoka,1996)

(10)

于是系統(tǒng)的相速度可以表示為

(11)

圖2 非靜水壓力梯度的影響虛線表示初始波的形態(tài)，實線表示當前波的形態(tài)，空心箭頭表示額外水平壓力梯度方向.Fig.2 The influence of non-hydrostatic pressure gradientThe dashed line represents the shape of the initial wave,the solid line represents the shape of the current wave,and the hollow arrow represents the direction of the different horizontal pressure gradients.

1.2 數(shù)值方法

1.2.1 網(wǎng)格劃分

在二維空間，采用節(jié)點中心格式設置空間交錯的均勻網(wǎng)格，變量η、u和v分別在三套網(wǎng)格中進行求解，如圖3所示.

圖3 交錯網(wǎng)格節(jié)點與控制體積示意圖三個物理量η，u和v所在計算節(jié)點分別表示為圓形、三角形和正方形，所在控制體積分別填充正交線、左斜杠和右斜杠表示.Fig.3 Schematic diagram of node and control volume of the computational gridThe nodes of three independent variables η,u,and v are expressed as circular,triangular,and square respectively.An orthogonal line,left slash,and right slash represents the corresponding control volume.

1.2.2 連續(xù)性方程

對連續(xù)性方程(1)，在時間上采用向前歐拉格式，在空間上使用有限體積法.于ηi,j節(jié)點處，在以節(jié)點為中心的控制體積ΔVi,j內(nèi)對連續(xù)性方程進行積分，引用高斯定理得

+Δx[(Dv)n-(Dv)s]=0,

(12)

式中S為控制體積ΔVi,j的外邊界，n為邊界S外法線的單位向量，下標e、w、n、s分別指示邊界S的東、西、北、南四條邊界.Δx和Δy分別為在x和y方向上的空間步長，控制體積ΔVi,j=Δx×Δy.邊界處的水深使用相鄰節(jié)點的水深平均值，得到

(13)

式中Δt為時間步長，k為當前時間步.

1.2.3 動量守恒方程

形如上述動量守恒方程的保守部分可以表示為關于變量φ的形式

(14)

式中ν為擴散系數(shù).于φi,j節(jié)點處，在以節(jié)點為中心的控制體積ΔVi,j內(nèi)對方程進行積分得到

(15)

再在時間上采用向前歐拉格式得到：

(16)

式中m為S的各個邊界組分(e、w、n、s)，Δm為邊界長度，φm為邊界處的變量值，整理后得

(17)

l表示邊界處的切線方向單位向量.

對流項中C表示Courant數(shù)，展開有

上標“+”表示邊界處速度方向與邊界法線方向相同，上標“-”表示邊界處速度方向與邊界法線方向相反，即可以概括為

(19)

對于擴散項展開則有

(20)

對于控制體積邊界處的值可以由其上下游處的節(jié)點值插值得到，如東邊界處的值φe可以表示為如下形式：

(21)

式中函數(shù)Ψ稱為通量限制器，與迎風側的梯度與逆風側的梯度比值r有關，r表示為

(22)

式中φu表示當前邊界處的上游節(jié)點值，φu u表示φu的上游節(jié)點值，φd表示當前邊界處的下游節(jié)點值，φd d表示φd的下游節(jié)點值，如圖4所示.

圖4 插值示意圖黑色實心圓和空心圓分別表示節(jié)點位置和相應節(jié)點處物理量的值，黑色實心三角和空心三角分別表示邊界位置和相對應的物理量的值.Fig.4 Schematic diagram of interpolationThe black solid circle and the hollow circle represent the position of the node and the value of the corresponding physical variables.The black solid triangle and the hollow triangle represent the position of the boundary and the value of the corresponding physical variables,respectively.

根據(jù)不同的使用需要，通量限制器函數(shù)Ψ具有不同的形式，最常被使用的格式見表1(Versteeg and Malalasekera,2007).

表1 常見的通量限制器函數(shù)表達形式(Versteeg and Malalasekera,2007)Table 1 Most popular limiter functions found in the literature (Versteeg and Malalasekera,2007)

Sweby(1984)根據(jù)r-Ψ關系給出了格式符合TVD方案的充要條件，還提出了對達到二階以上高階精度格式的要求.本文采用的QUICK格式可以保證TVD特性，且至少具有二階精度.在均勻正交網(wǎng)格下，對于擴散項，使用TVD-QUICK格式有

(23)

于是得到方程保守部分的最終形式

(24)

對于控制方程中的Boussinesq項在空間上采用有限差分法，我們將在算法測試中對一維方程的Boussinesq項進行具體離散.

1.2.4 穩(wěn)定性分析

(25)

1.3 涌潮的生長

初始速度場分布與初始波形為

(26)

式中u0和a均為常數(shù).

為簡便操作，取初始靜水深h=1，得到控制方程

(27)

η，u，x，t均為無量綱變換之后的變量，其中頻散項有

(28)

對方程(27)使用數(shù)值格式中介紹的方法進行數(shù)值離散，將Boussinesq項離散之后得到

(29)

我們分別計算頻散系數(shù)β=0和β=1/5的情況，并參照不考慮頻散項的靜水壓力模型，測試中選擇u0=0.1，dx=0.6，dt=0.25，a=6(Peregrine,1966)，計算結果見圖5.

圖5 波狀涌潮的生長過程.從下往上的曲線簇分別表示時間從0開始，以5增長直至20的計算結果灰色實線、紅色實線分別表示使用Peregrine 1966年淺水模型和Boussinesq模型的數(shù)值計算結果，黑色虛線、黑色實線分別表示使用文中所述方程及算法在β=0和β=1/5時所得的計算結果.Fig.5 Growth of an undular bore.The curves from bottom to top represent the calculation results of time from 0 to 20 with 5 as the interval The gray solid lines and the red solid lines represent the computation results by Peregrine 1966 using the shallow water model and the Boussinesq model respectively,while the black dashed lines and the black solid lines represent the computation results using the equation and algorithm described above when β=0 and β=1/5 respectively.

從圖5中可以看到，在模擬時間內(nèi)，僅考慮靜水壓力的數(shù)值模型計算結果(灰色實線)，波前緣逐漸變陡，但波后方仍舊保持平滑，波前緣漲幅保持在0.1；使用Boussinesq型方程的計算結果在初始一段時間內(nèi)和靜水壓力理論的結果重合，而后波前緣逐漸隆起形成波峰，在波后方出現(xiàn)波谷，再現(xiàn)了波的頻散(二次自由面起伏)特征，形成了波狀涌潮.當頻散系數(shù)β=0時，使用如前所述的算法計算結果(黑色虛線)與Peregrine 1966年的測試結果(紅色實線)相當吻合，波前緣漲幅在0.15左右.而當頻散系數(shù)β=1/5時考慮了更高階的垂直壓力分布改變對水平壓力梯度變化的影響，因此計算得出的漲幅也會更大(黑色實線)，波前緣漲幅在0.17左右.由測試結果可見，在涌潮模擬中對靜水壓力的校正是非常重要的，方程需要考慮頻散項的影響.

1.4 模擬區(qū)域與初邊值設置

海面的潮汐變化可以分解為許多簡單的、規(guī)則的簡諧振動，每個振動稱為一個分潮，大多數(shù)地區(qū)的潮汐主要成分是由月球引起的主太陰半日潮.錢塘江河口潮汐現(xiàn)象即為一種半日潮.2000年,浙江省錢塘江管理局與浙江省河口海岸研究所對錢塘江河口進行了一次大規(guī)?？疾煨杂砍庇^測，錢塘江河口在地形上有兩個顯著的特點(林炳堯等,2002),一是杭州灣河口呈喇叭形；二是自乍浦往上至倉前、七堡之間,存在宏大的沙坎,其頂部比乍浦河床高出10 m左右.錢塘江潮波的變化、涌潮的形成和發(fā)展與這兩種地形特點密切相關.澉浦是我國潮差最大的海域之一,但是當潮波通過杭州灣時,變形的程度并不明顯，澉浦下游的潮位變化過程仍然接近簡諧曲線.一般情況下涌潮自大缺口至鹽官之間形成、發(fā)展，最終在七堡附近湮滅，大潮時期涌潮最遠可到達杭州以上的聞家堰至富陽河段(Pan and Lu,2011).

根據(jù)實際涌潮的特點，我們的模型區(qū)域上游邊界設置在閘口處，下游邊界設置在乍浦下游，計算區(qū)域范圍取值為120.15°E—121.25°E/30.0°N—30.7°N(圖1)，地形數(shù)據(jù)水平分辨率精度約90 m.河流東邊界給定水面高度變化，使用與實際記錄相近的正弦波形函數(shù)來模擬潮汐引起的強迫場(林炳堯等,2002)，根據(jù)杭州灣的潮差分布(Fan et al.,2015)以及實測資料中的錢塘江涌潮的特性(林炳堯,2008)，選擇振幅amp=2 m，周期period=12 h.河流西邊界在模擬的時間內(nèi)假定為水深恒定速度自由，即水位邊界條件為給定的恒定水深，速度邊界設置為自由傳輸邊界.模型考慮水波的上岸淹沒情況，因此網(wǎng)格的干濕情況是時刻變化的，形成移動的干濕邊界，網(wǎng)格干濕判斷條件如下，取最小水深Dmin=0.01 m，即有

(30)

經(jīng)各式模型的測試與實測校對表明，錢塘江河口底部的曼寧系數(shù)較小，約為0.01(潘存鴻等,2008)，或者更小，取曼寧系數(shù)n=0.01，頻散系數(shù)β=1/5，初始時間步長Δt=0.5 s.

2 模擬結果與分析

我們模擬了潮波進入錢塘江河口，上行直至上游閘口處，形成涌潮的整個發(fā)展過程.結果主要展示涌潮經(jīng)過不同河段時的波高場以及速度場的變化，旨在從涌潮高度、涌潮傳播速度、潮景等涌潮傳播表征進行分析，探究在不同區(qū)域河段涌潮的特征.將陸地區(qū)域的名稱，作為其所對應相鄰的錢塘江河段的暫代名稱，用于對涌潮經(jīng)過此河段區(qū)域的定性分析.對錢塘江涌潮的觀測資料顯示，涌潮的形成基本在尖山一帶，自大缺口至鹽官河段發(fā)展達到最大，再往上行,涌潮高度逐漸減小,直至消失,其上界一般在七堡附近，最遠也可到達杭州以上(譚維炎等,1995;林炳堯等,2002;潘存鴻等,2008).

因關注的重點在涌潮本身，我們將濕單元的數(shù)據(jù)進行提取并繪制成波高場圖和速度場圖.在我們的模型計算結果中，涌潮進入形成、發(fā)展和消亡三個階段的位置和實際記錄結果吻合.涌潮形成于尖山河段，該河段下游的水深在40 m左右，相對較深，此時的涌潮高度較小，潮頭并不明顯，潮波因江底沙洲和沿岸地形的影響分成兩股波，在尖山和大缺口河段間形成“交叉潮”(圖6a).隨后涌潮從大缺口進入了較為平直狹窄的鹽官甬道，大量水體涌入，表現(xiàn)出“一線潮”(圖6b)，且強度隨著深入不斷增大(圖6c).在老鹽倉拐角處，因轉折角度較大，涌潮被反射形成“回頭潮”(圖6d).鹽官上游之后的河段水深較淺，基本處于10 m以下，江底摩擦產(chǎn)生的能量耗散作用明顯，再加上彎曲地形與涌潮之間的碰撞，涌潮強度開始減弱.涌潮后續(xù)經(jīng)過幾個拐角上溯至七堡河段(圖6e)，之后進入閘口河段走向消亡階段(圖6f).2.2小節(jié)將會對三處典型潮景的涌潮特征進行具體分析.

2.1 涌潮特征的沿程變化

涌潮高度、涌潮速度以及Froude數(shù)是表征涌潮特性的三個重要指標.涌潮高度ΔH一般指涌潮引起的平均水面上漲高度，ΔH=hd-hu，hu、hd分別表示涌潮潮頭前方(即波前緣的上游方向)水深和潮頭后方(即波前緣的下游方向)水深，如圖7所示，相對漲幅表示為ΔH/hu.

使用涌潮潮頭傳播的速度描述涌潮速度，涌潮速度c的流量守恒方程和動量守恒方程可以表示為(以指向上游方向為流速的正方向)

c(hd-hu)=udhd-uuhu,

(31)

(32)

其中uu、ud分別表示潮前流速和潮后流速(圖7).

圖7 涌潮傳播示意圖uu、ud分別表示潮前流速和潮后流速，hu、hd表示潮前水深和潮后水深，c表示涌潮流速.Fig.7 Definition sketch of a bore propagationuu and ud are respectively the flow velocity before and after the tidal bore,hu and hd indicate the flow depths before and after the tidal bore respectively,c is the bore celerity.

由式(31)、(32)可以分別得到

(33)

(34)

聯(lián)合整理式(33)、(34)之后得到由潮前、潮后流速表示的涌潮傳播速度c為

(35)

(36)

由式(36)可以看到，F(xiàn)roude數(shù)與潮前水深和潮后水深的比值相關，F(xiàn)r越大，涌潮引起的前后水面漲幅落差也就越大.而c同時可以由Fr表示為

(37)

可見涌潮速度c與Froude數(shù)Fr、流速以及水深有關.我們選取靠近河道中心位置處的涌潮剖面，計算涌潮上述三個特征以及漲幅的沿程變化情況，具體計算結果見表 2.

表2 涌潮特征沿程變化Table 2 Variation of bore characteristics along the river

在不考慮老鹽倉處的“回頭潮”的情況下，我們將表中的數(shù)據(jù)按照錢塘江涌潮上溯的傳播方向繪制了涌潮特征沿程變化，見圖8.涌潮的Fr隨ΔH自涌潮的形成階段、發(fā)展階段直至消亡階段，整體呈現(xiàn)出先增大再減小的規(guī)律.自鹽官以上河段的水深相差不大，在10 m左右，因此c也整體呈現(xiàn)先增大再減小的規(guī)律.涌潮于鹽官河段處上溯，強度一直增大，于外五工段區(qū)域河道處規(guī)模達到最大，ΔH約為3.5 m，F(xiàn)r約為1.54，c約為9 m·s-1，此時段的涌潮潮頭陡立，自由表面可能失穩(wěn)，波開始破碎.

圖8 涌潮特征沿程變化曲線方塊、三角、圓分別表示涌潮速度c、涌潮高度ΔH和涌潮的Froud數(shù)Fr.Fig.8 Curves of variation of bore characteristics along the riverThe celerity c,bore height ΔH,and Froude number Fr are expressed as square,triangular,and circular.

涌潮經(jīng)過老鹽倉處的急轉彎，在反射和推擠作用下形成“回頭潮”，“回頭潮”的強度并不大，ΔH在0.7 m左右，F(xiàn)r接近1.繼續(xù)上溯的涌潮強度開始減弱但規(guī)模依舊較大，三工段至倉前的整個河段涌潮變化不大，ΔH約為2.5 m，F(xiàn)r接近1.4，c約為8 m·s-1.在經(jīng)過七堡至閘口河段，涌潮規(guī)模和強度都在減弱，F(xiàn)r減小到近1.2，ΔH也僅有1 m左右，涌潮處在消亡階段.閘口河段處的水深較七堡處要大，于是出現(xiàn)在閘口河段處的涌潮速度要稍微大于七堡河段處的涌潮速度，約為6 m·s-1.

從實測資料分析中發(fā)現(xiàn)，涌潮高度一般于鹽官河段達到極值，最大可達3 m以上，涌潮傳播速度在3～10 m·s-1之間，此后涌潮開始衰減，涌潮高度也逐漸減小，沿程的涌潮速度也呈現(xiàn)從小到大，再從大到小的變化規(guī)律，最大流速出現(xiàn)在鹽官河段(林炳堯,2008;潘存鴻等,2008)，可見我們的模擬結果和實際觀測現(xiàn)象所得的特點是吻合的.

2.2 涌潮的發(fā)展與潮景的形成

潮景的形成是涌潮不同發(fā)展階段的標志，我們復演并選取了最具代表性的三個河段位置處的潮景，分別是大缺口處的“交叉潮”，鹽官處的“一線潮”以及老鹽倉處的“回頭潮”.

2.2.1 交叉潮

尖山河段是錢塘江河口的一個轉折處，河口近似呈“之”字形，從此處潮波由較為深廣的水域經(jīng)大缺口進入淺窄且平直的鹽官河段.大缺口附近河口呈平躺的“U”字形，因泥沙淤積，錢塘江河口下游發(fā)育一巨大江中沙洲(圖1)，在大缺口下方尖山河段分汊河勢以及江中沙洲的影響下，此處常常發(fā)育“交叉潮”(潘存鴻等,2008).

潮波在通過沙洲時被分成了兩股波，一股沿著北岸向西北傳播，另一股沿著南岸向西南方向傳播(圖9a)，潮波帶來整個尖山河段水面的上漲，并未形成明顯的涌潮現(xiàn)象.兩股暗潮相互作用，使得西北向潮波得到加強，于尖山河段下游形成較為清晰整齊的波，涌潮形成(圖9b)，涌潮的方向指向河流的上游.事實上，涌潮的形成并不是瞬時的，而是一個連續(xù)的過程，因此無法具體準確地判斷涌潮形成的精確位置與時間，當肉眼可分辨涌潮形態(tài)時，我們即認為涌潮形成.在模擬時間的80 min左右，從圖9c的涌潮剖面中可以看到，尖山河段潮頭漲幅可以達到1.5 m，涌潮高度僅約0.3 m，計算出Froude數(shù)在1左右，涌潮形成初期波形平緩，波列整齊.

圖9 尖山河段涌潮圖(a—b)展示了不同時刻潮波引起的尖山河段的水面變化.紅色虛線表示涌潮的波前緣，紅色實線表示截取的涌潮剖面.圖(c)展示了圖(b)中的涌潮西北-東南方向的剖面形態(tài).Fig.9 Tidal bore in JianshanFigures (a—b)show the variation of water surface at Jianshan at different times respectively.The red dashed line indicates the wavefront of the tidal bore wave,and the red solid line indicates the selected tidal bore profile.Figure (c)shows the bore profile in figure (b)from NW to SE.

在波狀涌潮上溯的過程中于尖山至大缺口河段形成“交叉潮”(圖10).尖山至大缺口河段是一個近似由西向北的轉彎，在潮波與南岸邊發(fā)生反射推擠作用的影響下，南岸形成了一股近似由南向北傳播的涌潮(圖10a)，另一股涌潮沿著東南向西北行進，兩股涌潮的潮頭波高皆約為1.5 m，兩股波前形成一個鈍角的“人”字形.此時涌潮高度較小，波形完好，顯示出良好的波的特征.兩股涌潮繼續(xù)上行，進入大缺口河段入口，兩組波列相互干涉，波前形成了一個近似交叉的“十”字，即所謂的“交叉潮”(圖10b)，波前緣的夾角依然是鈍角.在“交叉潮”后方清晰可見由干涉產(chǎn)生交錯排列的明暗紋，波峰處漲幅在1.5 m左右，波谷處漲幅在1～1.5 m之間，波后的水面上漲約為2 m.“交叉潮”向大缺口行進，波前緣夾角逐漸減小(圖10c)，向北傳播的涌潮接近北岸，涌潮的形態(tài)也重新呈現(xiàn)“人”字形.

圖10 “交叉潮”波高場圖(a—d)展示了不同時刻尖山河段的水面變化，紅色虛線表示涌潮的波前緣，涌潮從尖山河段上溯至大缺口河段的過程中形成“交叉潮”.Fig.10 Wave height field of cross-shaped tidal boreFigures (a—d)show the variation of water surface at Jianshan at different times respectively.The red dashed line indicates the wavefront of the tidal bore wave.The cross-shape tidal bore is formed when the bore goes upstream from Jianshan to Daquekou.

“交叉潮”形成時，涌潮波的傳播速度雖然接近20 m·s-1，但流場的速度較小，約為1 m·s-1，波的傳播速度比流體速度大得多.速度方向的分布和兩股潮的行進方向一致，波峰的位置處流場速度較大，而波谷位置處的流場速度較小，流場也呈現(xiàn)出與波場類似的明暗條紋分布，如圖11所示.

圖11 “交叉潮”速度場箭頭方向表示速度方向，箭頭長度表示速度大小，紅色虛線表示涌潮的波前緣位置，速度場因速度方向和大小的變化呈現(xiàn)出與波紋類似的明暗條紋.Fig.11 Velocity field of crossed tidal boreThe direction of the arrow indicates the direction of velocity while the magnitude indicates the strength,and the red dashed line indicates thewavefront of the tidal bore wave.Due to the variation of strength and direction,the velocity field presents light and dark fringes.

2.2.2 一線潮

鹽官以上的上游河段水深較淺，在10 m以下，因此涌潮在進入淺窄且平直的鹽官河段后，水面高度迅速上漲，潮頭走向幾乎垂直于岸邊，且波列以近乎平行的姿態(tài)上溯，于鹽官河段附近形成較為理想的“一線潮”潮景(圖12).“一線潮”是涌潮的發(fā)展階段，涌潮強度開始增大.鹽官河段處形成強度較大的“一線潮”，涌潮高度約為3 m.之后“一線潮”保持著相似的形態(tài)，強度繼續(xù)變大，通過整個鹽官河段到達外五工段河道區(qū)域.

圖12 “一線潮”波高場圖(a—c)展示了不同時刻潮波引起的鹽官河段的水面變化，涌潮通過鹽官河段時，潮頭處波高等值線變密，潮頭急速抬升，形成“一線潮”，紅線表示截取的涌潮剖面.Fig.12 Wave height field of thread-shape tidal boreFigures (a—c)show the variation of water surface at Yanguan at different times respectively.When the tidal bore reaches the Yanguan outlet,the wave height contour at the bore head becomes denser and the bore head rises steeply,then the thread-shape tidal bore forms.The red line indicates the selected tidal bore profile.

從波形剖面可以清楚地看到涌潮強度不斷增大的過程.涌潮在進入鹽官河道入口時，依然可見波前緣后方一系列形態(tài)良好的波(圖13a).較為平緩的波前緣成長為陡立的潮頭，且潮頭前的水面會出現(xiàn)略微下降，在鹽官處形成明顯的間斷面(圖13b)，計算得到涌潮的Froude數(shù)約為1.45.直至在鹽官河段的出口處，外五工段的附近，涌潮的潮頭一直在上漲，達到將近4 m(圖13c)，F(xiàn)roude數(shù)達到1.54.

圖13 “一線潮”涌潮剖面展示了圖12中選取的涌潮剖面，涌潮自東向西通過鹽官河段，潮頭逐漸陡立，涌潮強度增強，二次面起伏現(xiàn)象明顯，具有強頻散作用.Fig.13 Profile of thread-shape tidal boreThis figure shows the tidal bore profile selected in Fig.12.The tidal bore passes through the Yanguan outlet from east to west,the tidal head becomes steep,the intensity of the tidal bore increases,secondary free-surface undulations are reproduced,and the bore has a strong dispersion effect

涌潮經(jīng)過鹽官河段的速度場分布近乎與岸邊平行，涌潮的波前緣存在明顯的速度分界.在鹽官河段入口附近，兩股流場交匯處的速度相近，速度大小約為1 m·s-1，波前緣處流體速度接近為零(圖14a).隨著涌潮深入，潮波進入甬道，潮后速度逐漸增大，于鹽官處形成“一線潮”時潮后的下游流場速度約為3 m·s-1(圖14b).“一線潮”上溯速度繼續(xù)增大，鹽官河段的出口處連接老鹽倉拐角，因此在外五工段附近，南岸的速度略微大于北岸(圖14c)，靠近南岸的流場速度在3 m·s-1左右.經(jīng)計算的鹽官河段涌潮速度在6～9 m·s-1，與流體速度僅相差幾倍，可見涌潮已不再是簡單的淺水模式，具有很強的頻散效應.

圖14 “一線潮”速度場圖(a—c)展示了不同時刻鹽官河段的速度分布，箭頭方向表示速度方向，箭頭長度表示速度大小，紅色虛線表示涌潮的波前緣位置.Fig.14 Velocity field of thread-shape boreFigures (a—c)show the velocity distribution at Yanguan at different times respectively.The direction of the arrow indicates the direction of velocity while the magnitude indicates the strength. The red dashed line indicates the wavefront of the tidal bore wave.

2.2.3 回頭潮

老鹽倉拐角處夾角接近130°，“一線潮”于拐角處發(fā)生反射，在老鹽倉形成了“回頭潮”的潮景.彎道處造成的碰撞和反射使得涌潮波與岸邊作用加強，且較淺的水深使得江底摩擦作用顯著，涌潮能量被耗散，于是在“回頭潮”出現(xiàn)之后，涌潮強度開始減弱.涌潮發(fā)生反射時，整個流場由下游來潮、反射波以及上游河流三個部分的流場組成(圖15b)，外拐角處水面上漲最大可達4 m.沿著河道上行的涌潮疊加部分反射波繼續(xù)上溯，進入三工段區(qū)域河道.三工段河道與鹽官河段相似，狹長且河面寬度較為均勻，涌潮進入河道后，以“一線潮”的形態(tài)前進，其規(guī)模較鹽官河段的“一線潮”要小，涌潮強度開始減弱，涌潮高度減小到約為2 m(圖15c).反向往下游傳播的 “回頭潮”，潮頭漲幅在3 m左右，涌潮高度約為0.7 m(圖15d).涌潮經(jīng)過老鹽倉所帶來的水面上漲在3 m左右，老鹽倉內(nèi)拐角內(nèi)側岸邊處的一些洼地被水填滿.“回頭潮”發(fā)生后，上溯涌潮的強度開始減弱，在經(jīng)過七堡河段的拐角之后，涌潮進入消亡階段.

圖15 “回頭潮”波高場圖(a—d)展示了不同時刻潮波引起的老鹽倉河段的水面變化，紅色虛線表示涌潮的波前緣，涌潮在老鹽倉處發(fā)生反射，形成自西向東傳播的“回頭潮”.Fig.15 Wave height field of returned tidal boreFigures (a—d)show the variation of water surface at Laoyancang at different times respectively.The red dashed line indicates the wavefront of the tidal bore wave.The tidal bore is reflected at Laoyancang,returned tidal bore forms and spreads from west to east.

老鹽倉拐角處，上下兩股流場交匯，速度場的分布情況和波高場分布相似(圖16)，由來潮、上游河流以及反射波三個速度場組成，涌潮的波前緣存在明顯的速度分界.反射波的速度方向仍然是總體指向上游的，拐角處的速度最小，僅在1 m·s-1左右(圖16b).指向上游的老鹽倉反射波和上溯的涌潮相疊加，以約為3 m·s-1的速度進入三工段區(qū)域河道，與上游的速度流場的交界處呈“一”字形態(tài)(圖16c).部分反射波形成“回頭潮”反向向東傳播，減慢了上溯流場的速度，“回頭潮”的強度也不斷減弱，涌潮同樣具有很強的頻散效應.

圖16 “回頭潮”速度場圖(a—d)展示了不同時刻老鹽倉河段的水流速度分布，箭頭方向表示速度方向，箭頭長度表示速度大小，紅色虛線表示涌潮的波前緣位置.Fig.16 Velocity field of the returned tidal boreFigures (a—d)show the velocity distribution at Laoyancang at different times respectively.The direction of the arrow indicates the direction of velocity while the magnitude indicates the strength.The red dashed line indicates the wavefront of the tidal bore wave.

3 結論與展望

通過對錢塘江河口涌潮的數(shù)值模擬研究，我們得到以下結論：

(1)錢塘江河口涌潮的形成和發(fā)展以及潮景的產(chǎn)生和錢塘江的地形特征密切相關，模擬結果與實際觀測中得到的結果可以很好吻合.潮波自錢塘江河口經(jīng)過江底沙洲的作用以及與江岸間的相互作用，在尖山至大缺口河段發(fā)育形成“交叉潮”；在鹽官河段，地形較為平坦，兩岸長直且近似平行，涌潮由較為深廣的水域進入淺窄的甬道時，波列逐漸均勻，形成“一線潮”；在老鹽倉拐角處，涌潮發(fā)生反射，形成“回頭潮”；“一線潮”潮景不僅發(fā)生在鹽官河段，三工段附近也出現(xiàn)了波動形態(tài)良好的波狀涌潮，同樣呈“一”字形.

(2)本次模擬復演了最具代表性的“交叉潮”、“一線潮”以及“回頭潮”三大潮景，鹽官河段“一線潮”表現(xiàn)的頻散效應最強.潮景的形成標志著涌潮進入不同發(fā)展階段，“交叉潮”出現(xiàn)在涌潮初期，F(xiàn)roude數(shù)接近1，涌潮高度在0.5 m左右；涌潮由大缺口進入鹽官河段之后形成“一線潮”，進入涌潮發(fā)展階段，并達到最大規(guī)模，“一線潮”的Froude保持在1.45以上，涌潮高度最大可達3 m，且涌潮的頻散作用較強；隨后涌潮到達老鹽倉形成“回頭潮”，同時涌潮強度被減弱，繼續(xù)上溯的涌潮進入三工段河道區(qū)域形成規(guī)模較小的“一線潮”，F(xiàn)roude數(shù)在1.3～1.4之間，涌潮高度約為2 m；七堡之后，上溯的涌潮進入消亡階段，F(xiàn)roude數(shù)減小至1.3以下，涌潮高度約為1 m.

(3)錢塘江涌潮過程中，F(xiàn)roude數(shù)和涌潮高度呈正相關，涌潮形成至進入消亡階段，F(xiàn)roude數(shù)隨涌潮高度呈現(xiàn)逐漸增大再逐漸減小的規(guī)律.尖山河段處涌潮Froude數(shù)接近1，鹽官河段涌潮的Froude數(shù)最大可達1.54，老鹽倉之后，涌潮的Froude數(shù)逐漸減小.

(4)Boussinesq型方程可以很好地再現(xiàn)錢塘江河口涌潮的形成、發(fā)展和消亡階段，對于鹽官河段處這類具有較強頻散效應的涌潮，展現(xiàn)了涌潮中的二次自由面起伏現(xiàn)象.

我們在前人工作的基礎上，進一步發(fā)展了錢塘江涌潮數(shù)值模型，使用考慮非線性、頻散以及耗散作用的Boussinesq型方程對錢塘江河口涌潮進行模擬，模擬結果與實際觀測中體現(xiàn)的特征相吻合.

在后續(xù)的工作中，考慮使用Boussinesq型方程結合更高精度的地形數(shù)據(jù)以及實時的涌潮觀測資料，可以運用在鞏固堤防、維護河道、保護涌潮等數(shù)值模擬方面的工作，推進日后對錢塘江涌潮的保護和深入研究探索.同時在強烈風暴作用下，風浪過程對水位的影響也十分重要(張舒羽和潘存鴻,2013;Luettich and Westerink,2015)，進一步工作中可以在本文使用的數(shù)值模式基礎上考慮風浪與潮流的相互作用，增強涌潮防災減災的實用性價值.

致謝感謝中國科學院大學計算地球動力學重點實驗室提供的計算平臺.