蘇文亮，鄒志利*，張慶民

( 1. 大連理工大學 海岸及近海工程國家重點實驗室，遼寧 大連 116024)

1 引言

海灣地貌對自然環(huán)境變化（河流來水來沙減少、海面上升等）和人類活動的影響（灘涂圍墾、跨海大橋建設等）十分敏感，認識海灣地貌的形態(tài)特征和演化規(guī)律、科學地利用海灣對自然環(huán)境保護和社會發(fā)展具有重要意義。世界上多數的海灣（或者潮流占優(yōu)的河口灣）是寬度呈收縮型的，如我國的黃茅海和杭州灣都是喇叭形的河口灣，國外這樣的海灣如阿曼灣、泰晤士河口等。其中，黃茅海河口灣的寬度在灣口大約為24 km，在灣頂大約為2 km，灣頂收縮為灣口的1/12。與此對應的另一類數量較少的海灣是寬度呈擴張型的海灣，如我國的湛江灣。

海灣地貌形態(tài)除了取決于海洋地質結構的塑造外還受到海灣內泥沙運動輸移所導致的淤積和沖刷的影響，后者是形成海灣地形形態(tài)的主要動力因素。這一動力因素包括水動力（潮流、波浪和徑流）、泥沙輸移和地形演變三者耦合相互作用。目前世界上的許多海灣的水底形態(tài)已經趨近于一種平衡狀態(tài)，即已經形成了海灣平衡剖面[1－2]，但人們對于這些海灣地形平衡形態(tài)的研究（特別是理論和計算模擬的研究）仍然處于發(fā)展階段。早期的研究表明[3]，由于漲潮流和落潮流的不對稱性，海灣的凈輸沙是沿海灣向陸地端輸移的，這樣如果海灣水底初始時是平底，泥沙將沿海灣淤積，在靠近陸地的海灣區(qū)域逐漸形成一個斜坡水底，最終形成向岸方向逐漸抬升的平衡剖面。Schuttelaars和de Swart[4]采用水深平均水流運動和泥沙輸移的一維數學模型來分析短尺度（長度遠小于潮流波長）矩形海灣地形演變，發(fā)現矩形海灣平衡剖面是一向岸方向水深均勻減小的平面斜坡。Schuttelaars和de Swart[5]于2000年又對長尺度（長度大于潮流波長）的矩形海灣進行了數值模擬，所得到的平衡剖面呈現出靠近入口處沖刷，然后沿著向岸方向呈現上凸型變化的形態(tài)，這與實際海灣地形剖面不是很符合，他們把這歸因于離岸邊界水深取為固定值所導致。Hibma等[6]在上述模型的基礎上對邊界條件進行改進，得到了與實際較為符合的海灣地形形態(tài)。這些研究都是針對寬度不變的理想情況的海灣，Lanzoni和Seminara[7]考慮寬度逐漸變小的海灣，采用了寬度積分的一維淺水方程來計算地形平衡剖面，所得結果為下凹的形態(tài)。Todeschini等[8]采用數值模擬綜合分析不同海灣收縮程度、潮位幅值和水底摩擦條件下的海灣平衡剖面，發(fā)現海灣收縮程度越大、潮位幅值越大，平衡剖面的長度將越短。與以上研究都集中在寬度不變和寬度變小的海灣形態(tài)不同，Meerman等[9]考慮了寬度變大的海灣，預測其平衡剖面為上凸形。在理論解方面，Prandle[10]研究了三角形截面海灣的平衡剖面解析表達式，給出的海灣水深h與離開灣口距離x的關系為h～。除針對三角形斷面這一特殊形式的海灣平衡剖面的理論解，目前還沒有見到針對一般形式海灣的有關解。

本文給出了實際可能存在的收縮型、擴張型和矩形3種海灣形態(tài)的平衡剖面的理論結果和數值模擬結果。推導了海灣平衡剖面和時均懸沙濃度的解析表達式，應用水深平均水平二維計算模型數值模擬了3種海灣平衡剖面和時均懸沙濃度，二者結果進行了對比，驗證了所得解析解的適用性。

2 海灣平衡剖面近似表達式

海灣形成平衡剖面的條件是：海灣泥沙輸移的凈輸沙率在潮灣分布均勻或者在各個位置都為0。這一條件是指潮周期平均的結果，其并不意味著在一個潮周期內的任何時刻的輸沙率都為0，因此，海灣平衡形態(tài)是一個相對的動態(tài)平衡。本節(jié)將通過理論分析來給出海灣平衡剖面的近似解析表達式。

2.1 理論分析

圖1 坐標系（以擴張型海灣為例）Fig. 1 Coordinate system (divergent embayment)

式中，μ為泥沙擴散系數?？紤]海灣長度相對于潮波波長非常小的情況，這樣的海灣內潮流的水面升高和速度的幅值變化很小，所以可以取u的幅值空間分布是近似均勻[4]，即

式中，u0為速度幅值；為潮流圓頻率。進一步假定海灣寬度遠小于長度，有＜＜1，從而有以及（可由式（20）驗證）。對方程(5)取潮周期平均，在以上假定下可以在方程中忽略項和項，由于地形平衡時，時均（潮周期平均）含沙量時間變化率為0（），可得該方程時均后的首階近似為

2.2 矩形海灣

如前面所述，這里考慮的海灣為短尺度的，灣內潮位空間變化很小，于是可假設潮位沿空間分布均勻，整個海灣內水面有如下形式：

式中，h0為海灣入口處水深。上式事實上即是海灣平衡水深，因為其對應的速度u（式（6））和懸沙濃度C（式（8））滿足地形平衡方程（4）。由此可知，矩形海灣平衡時剖面呈平面斜坡形狀。Schuttelaars和de Swart[4]也推導出了以上矩形海灣的形如式（11）的平衡剖面表達式，但其所采用的懸沙濃度輸移方程與常用的式（5）不同，方程右端沉積和落淤項采用了水深積分的懸沙濃度，而這里采用的式（5）的濃度是水深平均體積含沙量。

以上推導僅應用了連續(xù)方程，沒有應用水流動量方程，這是因為對短尺度海灣，水面變化非常小，其空間梯度接近于0，動量方程可自動滿足[4]。

2.3 收縮型和擴張型海灣

對收縮型和擴張型海灣，考慮水深和寬度積分的流動連續(xù)方程（線性化的）

將速度u表達式（6）代入式（12），可得水深h控制方程

與矩形海灣均勻水面不同，收縮型或擴張型海灣由于寬度的變化可能使得水面也存在沿海灣長度方向的相應變化。因為海灣寬度變化是呈指數型的（式（1）），這里假定水面升高沿海灣方向變化（x方向）也具有指數型的，并且隨時間呈現周期性的正弦變化（對應于潮流流速的表達式（6），參見式（10）），即具有以下形式的表達式：

式中，η0為海灣入口處潮位的幅值；β為表達水面沿海灣變化的待定參數，這一參數的值對不同的海灣是不同的，在本研究中它們的取值將通過擬合數值解來確定。將以上水面的表達式代入方程(13)，可得水深的解為

式中，h0為海灣入口處水深。該解包括齊次解式（17）右端第一項和特解式（17）右端第二項。齊次解系數D是由海灣入口水深為h0確定的。式中潮位對平衡剖面的影響是由參數β來反映的，海灣寬度變化對平衡剖面的影響是由參數α反映的。上式給出的水深變化整體上呈現指數變化趨勢，既可表達上凸地形，又可表達下凹地形，具體形式可從第3節(jié)的有關結果中看出。

3 數值模型及驗證

數值計算采用水深平均水平二維計算模型。模型包括水流連續(xù)方程、動量方程、懸沙輸移方程和地形演變方程。各方程具體形式如下：

連續(xù)方程

式中，H為總水深，H=h+η；g為重力加速度；Cz為謝才系數，Cz=h1/6/nr，nr為糙率，這里取 0.1；np為孔隙率，取值為0.4；zb為底床高程；為飽和懸沙濃度；，。后面兩個參數的取值對應泥沙中值粒徑為0.2 mm[11]。

計算的邊界條件為：離岸邊界為給定潮位，海岸邊界為考慮水面爬坡的干濕邊界，兩側邊界為固壁邊界。

水動力方程和泥沙輸移方程采用交替方向隱式差分格式(ADI格式)求解，地形演變方程采用Liu等[12]在ENO格式基礎上提出的加權基本無振蕩格式(WENO格式)，該格式不僅具有ENO格式在間斷區(qū)分辨率高、在光滑區(qū)計算精度高的優(yōu)點，還具有截斷誤差階數較高的特點，可精確模擬地形方程可能存在的間斷解。

以上模型計算結果的驗證是通過與Todeschini等[8]的寬度積分的水平一維模型的計算結果進行對比來實現的。Todeschini等[8]應用該模型計算了一收縮型海灣的平衡剖面，海灣長L=40 km，寬度參數的值為0.88（式（1）），入射潮波為規(guī)則半日潮，潮位幅值為2.0 m。計算中，入口處水深h0=12.5 m，初始地形為平面斜坡，坡度為h/L=1∶3 333。與Todeschini等[8]的計算一致，岸線處邊界沒有考慮水面爬坡，所以關閉了干濕邊界條件，而是取海岸邊界的靜水水深為0.5 m。在達到平衡狀態(tài)時，潮位幅值、速度幅值和水底剖面的對比結果如圖2所示。圖中結果表明，本文水平二維模型的結果與Todeschini等[8]的水平一維模型的結果符合良好。

4 數值計算結果

下面將應用以上介紹的計算模型對收縮型、擴張型和矩形3種形態(tài)的海灣的平衡剖面分別進行數值模擬，同時利用所得結果驗證第2節(jié)中海灣平衡剖面和時均懸沙濃度的解析解，式（15）和式（20）。

4.1 收縮型海灣

圖2 水面幅值、 速度幅值、水底剖面h計算結果和水平一維模型的結果(Todeschini等[8])的對比Fig. 2 Comparison of the calculation results of surface elevation, velocity amplitude , bottom profile with horizontal one-dimensional model (Todeschini et al.[8])

海灣平衡剖面的數值模擬一般都是以百年為單位，屬于長期模擬。因此在數值計算中采用Roelvink[13]所述的加速因子。其原理是在一個潮周期計算完成后，以這個周期的地形變化的數倍代替數個周期的地形變化，以這個周期的潮位和速度代替數倍潮周期的潮位和速度。計算表明，該海灣的計算模擬在95年后可達到平衡狀態(tài)，對應的水面、速度和水底剖面結果如圖4所示。圖中水面幅值沿空間變化很小，速度幅值沿空間變化整體上也是很小，基本與第2節(jié)中的速度假設式(6)相符合，只是在靠近岸線（x=20 km）區(qū)域出現幅值降低，這是由于此處水面升高η遠小于水深h的假定已經不是很滿足，存在非線性影響。圖中水底剖面水深較大的前半段可以看出明顯的下凹形態(tài)，后端靠近岸線區(qū)域存在地形隆起，其形成原因是由于岸線附近是潮間帶，實際上存在泥沙凈淤積。計算中將岸線處理為干濕邊界，所以模擬出了在潮間帶區(qū)域的泥沙向岸凈輸運產生的泥沙沉積。

圖3 收縮型海灣地形幾何形態(tài)(a)和初始地形(b)Fig. 3 Geometry(a) and initial terrain (b) of convergent embayment

圖4 收縮型海灣速度幅值(a)、水面幅值空 間分布和地形剖面h演變(b)Fig. 4 Distribution of velocity amplitude(a), surface elevation and terrain profile evolution(b) in convergent embayment

圖5 收縮型海灣時均濃度差和懸沙濃度的空間分布Fig. 5 Distribution of suspended sediment concentrationand difference in convergent embayment

由圖5和圖6也可看出收縮型海灣水面、時均懸沙濃度和水深空間變化特征：水面幅值和時均懸沙濃度整體是略有變大，中間呈輕微下凹的形態(tài)；對應的水深h空間變化特征是呈下凹趨勢的。

4.2 擴張型海灣

取擴張型海灣的長度與上面收縮型海灣的算例一樣為L=40 km，入口處寬度為=5 km，寬度參數為α=-0.88。入射潮波也與收縮型海灣的算例一致：幅值為2.2 m的規(guī)則半日潮。海岸邊界為考慮水面爬坡的干濕邊界。初始水底剖面和上述收縮型海灣的相同，海灣平面幾何形態(tài)如圖1所示。

計算中發(fā)現擴張型海灣地形達到穩(wěn)定需要130年的時間，比收縮型海灣的95年要長很多。達到平衡后的水面幅值、速度幅值和水底剖面如圖7所示，圖中也給出了在實現平衡剖面過程中不同時刻（20～100年）的水底剖面。水底地形的空間變化特征是：與收縮型海灣不同，平衡剖面整體上呈上凸形態(tài)，在海灣入口7.5 km長的區(qū)域出現沖刷，之后出現淤積。水面幅值沿空間變化很小，速度幅值在前半部分變化不大，基本符合第2節(jié)中速度具有常數幅值的假定，但岸線附近部分存在有一段向上的突起，這是由于岸線處泥沙淤積導致水體向岸爬坡所致。初始時刻和地形平衡時懸沙濃度差值在海灣內的分布如圖8所示，可見，在達到平衡狀態(tài)時已經變得非常小。

圖6 收縮型海灣水面幅值、平衡剖面h數值結果和解析解結果對比Fig. 6 Comparison between numerical result and analytical result of surface elevation and equilibrium profile in convergent embayment

圖7 擴張型海灣速度幅值 (a)、水面幅值空間分布和地形剖面h演變(b)Fig. 7 Distribution of velocity amplitude (a), surface elevation and terrain profile evolution (b) in divergent embayment

將水面解析解（16）與圖中水面的前半段進行擬合（圖9），可得參數取值為β=-0.7。將該參數值代入式（14）和式（20）可計算出理論結果的平衡剖面和時均懸沙濃度，分別如圖9和圖8所示。通過與圖中數值結果對比可以看出，二者相互吻合，只是水深結果在靠近岸線（h=0）的淺水區(qū)域存在一些差別，這是由于不能滿足建立解析解所采用的水面升高遠小于水深h的線性假定所致。

下面對圖7a中x=16 km附近區(qū)域出現的速度突起給出進一步的說明。事實上，當取水底地形為理論解（17）時，該突起將消失。為了說明這一點，將計算水域中20 km區(qū)域水深固定取為理論水底，然后重新計算速度結果。圖10給出了相應的速度值，圖中也給出了數值計算水深及其對應的速度曲線。可見，對應于理論水深的速度在海岸附近呈現漸進緩慢變化、沒有出現突起，而對應于數值計算水深的速度是存在這一突起的，這表明，后者的速度突起是由水底隆起所引起的。

圖8 擴張型海灣時均濃度差和懸沙濃度的空間分布Fig. 8 Distribution of suspended sediment concentrationand difference in divergent embayment

圖9 擴張型海灣水面幅值、平衡剖面h數值結果和解析解結果對比Fig. 9 Comparison between numerical result and analytical result of surface elevation and equilibrium profile in divergent embayment

圖10 取水底地形為理論解(15)時速度幅值分布及其與數值計算的速度分布的對比Fig. 10 Distribution of velocity amplitude correspond to analytical solution(15)

4.3 矩形海灣

取矩形海灣長度為L=40 km，寬度為B=2 km，入射潮波和收縮型、擴張型海灣的相同。初始地形同樣采用平面斜坡。相比于收縮型海灣和擴張型海灣，矩形海灣明顯更快的達到平衡剖面，在模擬90年時剖面就已經不再變化，結果如圖11所示。圖中水底平衡時明顯的表現為平面斜坡形態(tài)，這和理論結果式（11）一致。圖中潮位幅值和速度幅值變化很小，與第2節(jié)的常數幅值假設相符合，圖12給出了有關對比結果以及水底剖面h理論與數值結果的對比。懸沙濃度差如圖13所示，平衡時濃度差也是非常小的。該圖中也給出了平衡地形所對應的時均懸沙濃度的理論結果（式（20）中）與數值模擬結果的對比，二者符合很好。

5 討論

本文的分析采用了短尺度海灣的假定，這使分析過程大為簡化，該假定的突出特點是潮流速度沿海灣長度的變化可以忽略，其僅隨時間做周期性變化，該流場所導致的潮周期平均海灣水底剖面自動成為平衡剖面，因為隨時間正弦變化的潮流速度可以滿足海灣內潮周期平均含沙量不隨時間變化的條件（式（8））。由上節(jié)解析解結果與數值解結果對比可見，基于以上簡化所得到的u和η以及平衡水深h和含沙量的理論結果仍然是可以和數值計算結果相符合的。下面對其中的原理，即對建立這些解析解所涉及的有關近似處理的合理性做進一步說明。

圖11 矩形海灣速度幅值分布(a)、水面幅值和地形剖面h演變(b)Fig. 11 Distribution of velocity amplitude (a), surface elevation and terrain profile evolution (b) in rectangular embayment

圖12 矩形海灣水面幅值、平衡剖面h數值結果和解析解結果Fig. 12 Numerical result and analytical result of surface elevation and equilibrium profile

圖13 矩形海灣時均懸沙濃度和時間濃度差的空間分布Fig. 13 Distribution of suspended sediment concentrationand difference in rectangular embayment

（1）關于海灣潮波水位和速度表達式的選擇。本研究對潮波速度采用了式（6）表達形式，且取其幅值為常數，對潮位采用了式（15）和式（16）的表達形式，為了說明這兩個表達式的合理性，除考慮連續(xù)方程（12）之外還需要考慮海灣內潮波的動量方程。在連續(xù)方程中忽略非線性項（乘積項），在動量方程也同樣地忽略非線性項，并忽略海灣橫向速度v, 因為其相對縱向速度u很?。ǎ瑒t可得線性化的水深平均連續(xù)方程和動量方程為（考慮海灣寬度變化并不改變動量方程的形式，參見Todeschini等[8]）

式中，φ是水底摩擦導致的水位與速度之間的相位差（摩擦相位角）；和是幅值。k在這里并不是潮波的波數，而是一個與相位差 φ一樣的待定量，事實上，這兩個量都是(是潮流波數)階小量，在最后的水位和速度的表達式中是可以忽略掉的，但在方程求解過程中是需要的，這可從下面推導中了解到。將以上兩式代入方程（27）和（28）可求得方程解，見附錄中式（A8）、（A9）和（A10）。圖 14 給出了收縮型海灣幅值原來的表達式（16）與新的表達式（A8）的變化曲線以及在取不同水底摩擦系數時的變化曲線。由圖中結果可見，原來的表達式與新的表達式的結果具有相同的變化趨勢，但原來的表達式更與數值模擬結果接近。的結果中時的結果更與數值模擬結果接近。

圖14 水面幅值和速度幅值的理論解（收縮型海灣）Fig. 14 Theoretical solution of surface elevation and velocity amplitude (convergent embayment)

（2）關于水面η解析表達式中參數。由以上理論推導可以看出，當速度u的表達式已知時，仍然需要水面η和水深h二者之一是已知的才可確定另外一個量（h已知可確定η，或者η已知可確定h）。這表明，在線性理論的范圍內，表達水動力和地形演變耦合的方程系統(tǒng)還不是封閉的（可證，即使考慮三維流動結果也是這樣），只有考慮非線性系統(tǒng)式（21）至式（26）才可以使系統(tǒng)自行封閉，這是因為線性解一般只能給出水底地形的增長趨勢，但不能給出該增長達到平衡時的狀態(tài)，后者屬于有限幅值效應，需要對地形演變進行非線性分析才能得到。在這種情況下，需要引入一個參數（即以上解中的參數）來使系統(tǒng)封閉。該參數是取決于海底平衡剖面的長度大小，所以事實上它也是海底平衡剖面所包含的參數，只不過本研究將其轉化為水面表達式中的參數。因為海灣潮位比水底地形更容易觀測和確定，所以這樣的轉化是把水底地形形態(tài)這樣一個不容易確定的問題轉化為水面變化這樣一個容易確定的問題。

（3）關于線性理論的限制。上面的理論推導和第2節(jié)的理論分析都忽略了η和u自身乘積和相互乘積的有關項，即忽略了非線性項，這導致所得結果不能考慮非線性效應。波浪非線性的大小在短波（深水波）情況是由波高與波長之比（波陡）來表征的，在長波（淺水波）情況是由波高與水深之比來表征的[14]。潮波屬于淺水波，所以其非線性的大小是由比值來衡量的，即在的幅值遠小于水深h的情況下，線性理論才適用。這使得本研究所得理論解不適合h較小的靠近岸線附近區(qū)域。數值模擬結果表明，這一區(qū)域存在水面爬坡形成的泥沙淤積，導致了存在突出的地形隆起，這會影響速度幅值空間分布和時均懸沙濃度的空間分布。

6 結論

本文采用了短尺度海灣（海灣長度遠小于潮汐波長）假定研究了海灣平衡剖面的確定問題，由于海岸潮汐波長尺度一般在500 km量級（半日潮），所以本文研究結果適合于長度為十幾或者幾十千米的海灣。基于該假定，理論分析和數值模擬研究結果給出了收縮型、擴張型和矩形3種典型海灣平面形態(tài)所對應的海灣平衡剖面，得到了這3種類型海灣平衡剖面的解析解及其對應的時均懸沙濃度的解析解，并與水深平均二維數值模型的計算結果進行對比。研究所得結論如下：

（1）對于長度遠小于潮流波長的短海灣，海灣內潮流速度幅值可以假定為常數。這一假定為本文數值計算結果和線性理論分析結果所證實：矩形海灣平衡地形的流場速度幅值沿空間是均勻分布的，收縮型和擴張型海灣的速度幅值在海灣的水深較大區(qū)域也基本是均勻分布的，但在靠近海岸的淺水區(qū)域會受到局部水底淤積隆起和較大的水底摩擦的影響，呈現偏離均勻分布的變化。

（2）基于海灣內潮流速度為常數的假定建立了短海灣平衡剖面的解析解和對應的時均懸沙濃度解析解。解析解對應的收縮型和擴張型海灣的水面升高沿海灣存在微小空間變化，這與等寬的矩形海灣的均勻水面有所不同。所建立的海灣水面解析表達式含有一個待定參數β，表征海灣寬度變化的影響（β的值與寬度參數的值符號相同，二者對收縮型海灣大于0，對擴張型海灣小于0，對矩形海灣等于0），所以，在實際應用中，可根據具體海灣的水面觀測結果來確定。本文利用數值模擬結果確定了這一參數的值，基于這一水面表達式而得到的海灣平衡剖面的解析解及其對應的時均懸沙濃度的解析解與數值結果符合很好。

（3）所得到的平衡剖面對收縮型海灣呈下凹形態(tài)，對擴張型海灣呈上凸形態(tài)，平面形態(tài)居于二者之間的矩形海灣的平衡剖面呈平面斜坡形態(tài)。

（4）平衡剖面對應的時均懸沙濃度沿海灣分布整體呈現緩慢升高趨勢，在靠近海岸區(qū)域升高變快。

以上研究結果不適合海灣長度接近或大于潮汐波長的長尺度海灣情況，如何對這類海灣的平衡剖面進行解析分析是今后研究的一個值得關注的課題，目前該方面的研究還僅有數值方面的結果（Todeschini等[8]）。

附錄：

1 海灣水位和速度的求解

把式（A8）代入方程（A6）可得速度幅值表達式

式（A8）和式（A9）中

進一步可由方程（A2）和（A4）求得波數k和相位角 φ，為了避免復雜的推導過程，這里采用近似處理，即直接采用水深為常數時的結果，后者可以容易在方程（A1）至(A4)中取水深為常數而得到，結果為

為了考慮實際水深是變化的，可將以上表達式中水深h取為當地水深（變水深）。