李文博，李煊，李慶杰，陳默，王立鵬

（1.山東省海洋生態(tài)環(huán)境與防災減災重點實驗室，山東 青島 266100；2.國家海洋局北海預報中心，山東 青島 266100；3.海軍航空大學，山東 煙臺 264000；4.自然資源部煙臺海洋中心，山東 煙臺 264000；5.國家海洋局北海海洋技術保障中心，山東 青島 266100）

1 引言

臺風過程中伴隨的大風、巨浪和風暴潮等海洋災害會給人類海上活動及海岸設施帶來巨大威脅，極端天氣條件下海浪的研究及預報具有非常重要的意義。海表面拖曳系數(shù)是在高風速下進行海浪數(shù)值模擬的關鍵因子之一，選取合理的拖曳系數(shù)將直接決定海浪數(shù)值預報或后報的精準度[1-2]。

在海氣相互作用中，海面風應力是海浪產生的主要驅動力之一，海氣之間的動量通量通常由海表面拖曳系數(shù)或者海面空氣動力學粗糙長度來表征。最初研究認為海表面拖曳系數(shù)是常量，隨著研究的深入，人們逐漸發(fā)現(xiàn)海表面拖曳系數(shù)與風速之間存在一定的數(shù)學關系。GARRATT[3]通過研究給出了Charnock無量綱粗糙度表達式中常數(shù)α的取值以及海表面拖曳系數(shù)與海表面10 m 高度處風速（U10）的參數(shù)關系；SMITH[4]用渦相關法分析得出海表面拖曳系數(shù)隨風速（6～22 m/s）的變化規(guī)律；WU[5]提出了海表面拖曳系數(shù)的線性參數(shù)化方案，并認為在颶風條件下海表面拖曳系數(shù)隨風速的變化仍遵從線性關系。LARGE 等[6]利用雷諾動量耗散法對海表面動量通量進行了估算，認為海表面拖曳系數(shù)在風速為 4～11 m/s 時為常數(shù)，大于 11 m/s 后與U10呈線性增長關系。GEERNAERT[7]獲取了北海的觀測數(shù)據，利用渦相關法并考慮波齡的影響進行分析，得出了4～24 m/s 風速條件下海表面拖曳系數(shù)的線性表達式。以上研究均認為海表面拖曳系數(shù)或動力學粗糙長度隨著風速增加呈線性增長。然而近年來觀測資料、實驗室數(shù)據以及數(shù)值模擬結果均表明高風速下海表面拖曳系數(shù)的增長有所減緩，甚至出現(xiàn)減小的趨勢。ALAMARO 等[8]通過實驗得出在風速達到25 m/s 后海表面拖曳系數(shù)隨風速的增大而減??；MAKIN[9]利用Barenblatt提出的有限飽和沫滴懸浮層理論，建立了高風速下海表面動力學粗糙度模型，并指出飛沫水滴組成的飛沫邊界層是海表面拖曳系數(shù)在高風速下發(fā)生衰減的原因；HOLTHUIJSEN等[10]根據觀測指出海表面拖曳系數(shù)在風速達到40 m/s 后開始隨風速的增加而減小，當風速達到80 m/s時海表面拖曳系數(shù)將降為零。

在臺風浪的模擬過程中，風應力的準確估計及其對海洋上層影響的表征十分重要。海表面拖曳系數(shù)是海氣動量通量參數(shù)化研究的關鍵，因此本文著眼于高風速下海表面拖曳系數(shù)的參數(shù)化方案研究，期望得到高風速下更為準確且便于數(shù)值模式應用的海表面拖曳系數(shù)表達式，以提高臺風過程中數(shù)值模擬預報的準確性。

2 研究方法

2.1 海浪數(shù)值模式

海浪數(shù)值模擬和預報工作最早開始于20 世紀40 年代，經過前人的不斷探索與技術的發(fā)展，海浪數(shù)值模擬預報模式已經從第一代發(fā)展到第三代。荷蘭Delft理工大學（Delft University of Technology）對第二代海浪模式進行了改進和優(yōu)化，并研制開發(fā)出第三代近岸海浪數(shù)值模式（Simulation WAves Nearshore，SWAN）。經過多年的改進，SWAN 模式趨于成熟，目前已被廣泛使用。SWAN 模式采用全隱式有限差分格式和基于能量守恒原理的平衡方程，它將波浪進入淺水后的傳播特點考慮在內的同時保證計算過程中時空網格能夠相對穩(wěn)定地進行，而且計算空間網格和時間步長上不會受到牽制。在平衡方程的各源項中，SWAN 模式不僅考慮了風輸入項、四波相互作用、破碎和摩擦項等，還考慮了波浪的深度破碎作用和三波相互作用。該模式包括了當前海浪預報的較新成果[11-13]。

2.2 海表面拖曳系數(shù)

海氣相互作用的一個重要形式是海氣界面的動量交換過程，風速脈動會產生動量通量，海表面風應力是驅動上層海洋環(huán)流和風浪的主要因素，因此，如何表征風應力對海氣界面的影響至關重要。在海氣相互作用的研究中，通常開展海氣動量通量的參數(shù)化研究，并用塊體公式來計算海面風應力。在這一研究過程中，海氣動量通量一般由海表面拖曳系數(shù)或者海面空氣動力學粗糙長度進行參數(shù)化。公式如下：

式中：Cd,z為海表面Z高度處拖曳系數(shù)；Uz為海表面Z高度處平均風速；U0為海表面處平均風速，通常取值為0，故上式可寫為：

拖曳系數(shù)描述的是動量通量在海-氣界面的傳輸程度。在流體動力學中，阻力系數(shù)即表面拖曳系數(shù)Cd是用于量化流體環(huán)境（例如空氣或水）中物體阻力的無量綱量，較低的阻力系數(shù)表示物體具有較少的空氣動力學或流體動力學阻力，阻力系數(shù)總是與特定表面積相關聯(lián)。Cd公式如下：

式中：Fd是阻力；ρ是流體的質量密度；u是物體相對于流體的運動速度；A是作用區(qū)域。任何對象的阻力系數(shù)均受表面摩擦阻力和寄生阻力兩個基本因素影響，其中寄生阻力是在拖曳物體穿過大氣中流體介質過程中由氣態(tài)介質的空氣動力阻力引起。Cd并不是常數(shù)，而是作為流速、流動方向、物體位置、物體尺寸、流體密度和流體粘度的函數(shù)。流速、流體粘度和物體尺寸也被結合到被稱為雷諾數(shù)（Re，慣性力和粘性力的比值）的無量綱量中，以便量化給定條件的流動，Cd是Re的相關量。在可壓縮流體中，Cd同時也是馬赫數(shù)（Ma，表示超過邊界的流速與聲音的局部速度的比率的無量綱量）的函數(shù)。對于某些形狀，Cd只取決于Re、Ma和流動方向，當Ma很低時，阻力系數(shù)與Ma無關。

基于Cd與U10具有相關性這一前提，依據有限的實測資料和不同的參數(shù)化方法，前人給出了低風速條件下Cd隨風速線性變化的關系。圖1為不同研究者獲得的Cd與U10的線性關系圖。

圖1 Cd與U10線性關系圖Fig.1 The linear relationship between Cd and U10

低風速下Cd線性化參數(shù)方案基本成型，但高風速下Cd與U10的關系仍是需要繼續(xù)探討的問題。近年來基于外海觀測數(shù)據、實驗數(shù)據和數(shù)值模型數(shù)據的研究都顯示出高風速下Cd隨風速的增長趨于平緩，并且在風速的持續(xù)增長下Cd將出現(xiàn)減小的趨勢。究其原因，可能是由于白帽破碎和飛沫造成的。

在SWAN 模式中，默認使用WU[5]得出的關于風速變化的Cd表達式：

2.3 海表面拖曳系數(shù)的參數(shù)化方案

COARE 3.0[14]模型是國際上廣泛采用的海氣動量通量計算模型，其中海表面動力粗糙度長度有3種參數(shù)化方案，分別由 YELLAND 等[15]、TAYLOR等[16]和OOST 等[17]提出，本文分別稱為 YT96 方案、TY01方案和O02方案。

（1）YT96方案僅考慮了風速的影響：

式中：α是一個通用常數(shù)或者無量綱化粗糙長度；v是運動粘度；Z0為海表面動力粗糙度長度;u*為摩擦速度。

（2）TY01 方案在風速的基礎上同時考慮了波浪的影響：

式 中 ：HS為 有 效 波 高 ，0.015u10N)；LP為主波波長，tw= 0.729u10N;HS/LP近似為主波波陡；v為運動粘度；u10N為海面10 m風速。

（3）O02方案也考慮了波浪的影響，但是選取了與TY01方案不同的波浪參數(shù)：

式中：LP為主波波長；CP為主波向速度為波齡；v為運動粘度。

Cd的參數(shù)化方法通常有3種：渦相關法、廓線法和慣性耗散法。廓線法內容為：在高風速條件下達到穩(wěn)定的中性層結大氣，大氣邊界層的垂直風廓線可以由以下對數(shù)關系式表達：

式中：u(z)表示海表面Z高度處的風速；u*為摩擦速度；k為卡曼系數(shù)，一般取值為0.4；z0為海面粗糙度。在流體力學中，摩擦速度通常用來與真實速度進行比較，摩擦速度可以用來重寫剪切應力。根據摩擦速度的定義：

式中：τ 表示流體的任意層中的剪切應力；ρ表示流體的密度。結合式（2），拖曳系數(shù)可以寫為：

再結合式（8）可以得到海表面拖曳系數(shù)與海面粗糙度的關系為：

為研究高風速下Cd與風速的關系，本文使用美國國家浮標數(shù)據中心（National Data Buoy Center，NDBC）的高風速實測數(shù)據作為研究資料。本文選用了9 個颶風過程期間8 個浮標的12 份浮標數(shù)據。所選浮標的有效觀測數(shù)據中，最大風速達到40.9 m/s，風速在25～41 m/s 的有效數(shù)據超過25 個，這表明高風速區(qū)域海表面拖曳系數(shù)的研究范圍有了很大拓展，數(shù)據的準確度也有一定保證。根據TY01和O02方案并結合式（8）和式（10），可分別計算出兩種方案下的Cd，在此基礎上對Cd與U10的散點關系進行二項式擬合，即可得出兩種方案下Cd與U10的參數(shù)化公式。由于YT96 方案本身即為與風速相關的關系式，故不再對其進行擬合，將YT96 方案作為新的參數(shù)化關系的對比方案之一。圖2 和圖3 為TY01 和O02 方案計算所得的Cd及其擬合函數(shù)的曲線圖。

圖2 TY01方案所得Cd及其擬合曲線Fig.2 Cd and its fitting curve obtained from the TY01 scheme

圖3 O02方案所得Cd及其擬合曲線Fig.3 Cd and its fitting curve obtained from the O02 scheme

所得擬合函數(shù)曲線中，O02 方案中高風速區(qū)的Cd隨風速增大而增大的趨勢有所減緩，而TY01 方案在風速增大到一定程度后Cd的增長更為顯著。根據所得散點圖進行分段二項式擬合，得到兩個Cd與U10的參數(shù)化公式：

TY01方案為：

對比兩種方案所得的Cd與U10的散點圖及參數(shù)化公式，統(tǒng)計參數(shù)包括誤差平方和（Sum of Squares for Error，SSE）、相關系數(shù)（R2）、調整后的R2以及均方根誤差（Root Mean Square Error，RMSE），結果見表1。通過綜合比較分析可知，TY01 方案下由U10計算Cd存在兩個問題：Cd在高風速下的計算結果數(shù)值過大，有失準確；Cd與U10的變化趨勢在高風速處陡增，與實際相反。故本文淘汰TY01 方案所得的參數(shù)化公式，選擇O02 方案的擬合結果（以下稱新O02 方案）進行下一步的實測對比驗證。實測數(shù)據選自中國南海海洋觀測平臺，站點位置距離海岸約6.5 km，測量數(shù)據范圍全長67 m，水下14 m。

表1 TY01和O02兩個方案擬合函數(shù)對比統(tǒng)計結果Tab.1 Comparison of the statistical results of the fitting functions obtained by the TY01 and O02 schemes

將新O02 方案與YT96 方案以及根據實測u*和U10得出的Cd進行對比，結果見圖4。從圖中可以看出兩個函數(shù)曲線與實測結果的變化趨勢一致，數(shù)量級相同，故新O02 方案能夠合理地反映Cd與U10的變化關系。

圖4 YT96、新O02方案與實測數(shù)據對比圖Fig.4 Comparison chart between YT96 and new O02 schemes and measured data

新 O02 方案的R2為 0.997 0，SSE 和 RMSE 分別為2.273 3 × 10-5和8.250 1 × 10-5，數(shù)量級均為10-5，表明該方案下的參數(shù)化擬合函數(shù)與實測數(shù)據的擬合效果良好。

3 參數(shù)化拖曳系數(shù)在臺風浪數(shù)值模擬中的驗證

3.1 颶風“卡特里娜”和“瑞塔”介紹

風暴“卡特里娜”于 2005 年 8 月 23 日（北京時，下同）在巴哈馬群島附近生成，24 日增強為颶風后，于佛羅里達州南部以1 級颶風強度登陸，數(shù)小時內其強度迅速增強為5 級，近中心最高持續(xù)風速為150 n mile/h，29日在密西西比河口登陸后強度減弱為3級，之后系統(tǒng)加速向東北移動，31日在俄亥俄州轉化為溫帶氣旋。風暴“瑞塔”于2005 年9 月18 日在特克斯與凱科斯群島形成，20 日早上升級為1 級颶風，21 日進一步升級為5 級颶風，24 日在德克薩斯州薩賓帕斯與路易斯安那州卡梅倫之間的海岸地區(qū)登陸。雖然颶風“瑞塔”登陸時已減弱為3級颶風，但其最大風速超過190 km/h，給沿海一些中小城鎮(zhèn)造成了嚴重的水災和不同程度的財產損失。

3.2 海浪模型的設置和數(shù)據來源

本文選定的計算模擬區(qū)域為79°～98°W，18°～32°N。水深采用全球水深地形模型ETOPO1 數(shù)據，利用非結構網格，網格邊界分辨率是1/20°。采用美國國家海洋和大氣管理局（National Oceanic and Atmospheric Administration，NOAA）提供的浮標數(shù)據作為實測數(shù)據進行對比驗證。根據颶風路徑及NDBC 所提供的浮標位置找到6 個有效浮標數(shù)據，浮標經緯度位置及水深見表2。

表2 浮標信息表Tab.2 Buoys information table

本文SWAN 模式中的風場采用模型風場與CCMP（Cross-Calibrated，Multi-Platform）風場作為背景風場構造的合成風場，模型風場由圓形對稱風場Jelesnianski-II 模型表達式得出[18-19]，背景風場采用經過線性插值后的CCMP 風場，時間步長為1 h，空間分辨率為5′×5′。采用隨風速變化的權重系數(shù)e構造合成風場，構造公式為：

式中：Vm為模型風場；Ven為背景風場；c是一個與臺風影響范圍有關的參數(shù)；系數(shù)n一般取9 或10；r表示計算點到臺風中心的距離；Rmax為最大風速半徑。合成風場綜合了模型風場與背景風場的優(yōu)勢，在颶風中心附近模型風場中起決定性作用，而在臺風外圍背景風場中起主導作用。

誤差分析表明，颶風“卡特里娜”合成風場與實測風場的平均絕對誤差（Mean Absolute Error，MAE）為 1.93 m/s,平均相對誤差（Mean Relative Error，MRE）為20.4%；颶風“瑞塔”合成風場與實測風場的 MAE 為 1.78 m/s，MRE 為 18.8%。上述數(shù)據說明合成風場基本能夠反映真實風場的結構和變化。

3.3 數(shù)值實驗結果及分析

利用新O02 方案分別對颶風“卡特里娜”與“瑞塔”的臺風浪過程進行數(shù)值模擬，并與YT96 方案和SWAN 模式默認方案進行對比。圖5 為颶風“卡特里娜”過程有效波高實測數(shù)據與模擬值的對比，圖6為颶風“瑞塔”過程有效波高實測數(shù)據與模擬值的對比。根據圖像可知：

圖5 颶風“卡特里娜”過程有效波高實測數(shù)據與模擬值的對比Fig.5 Comparison of measured and simulated values of significant wave height during hurricane"Katrina"

圖6 颶風“瑞塔”過程有效波高實測數(shù)據與模擬值的對比Fig.6 Comparison of measured and simulated values of significant wave height during hurricane"Rita"

（1）對于颶風“卡特里娜”臺風浪過程，YT96 方案在浮標42001 和42002 處的模擬結果更好，而在浮標 42019、42035、42039 和 42040 處，新 O02 方案的參數(shù)化關系對臺風浪的模擬更接近實測值。

（2）對于颶風“瑞塔”臺風浪過程，YT96 方案在浮標42002 和42019 處的模擬結果更好，而在浮標42001、42035、42039和42040處，新O02方案的參數(shù)化關系對臺風浪的模擬更接近實測值。

（3）綜合來說，新O02 方案參數(shù)化關系的模擬效果優(yōu)于SWAN 模式默認方案，而與YT96 方案的模擬結果各有千秋。

通過颶風“卡特里娜”和“瑞塔”海浪數(shù)值模擬有效波高與實測數(shù)據的對比發(fā)現(xiàn)，新O02 方案中海浪數(shù)值模擬結果與實測浮標數(shù)據的趨勢較為一致，且與實測結果符合較好，說明新的拖曳系數(shù)參數(shù)化方案能夠滿足海浪數(shù)值模擬的應用要求。有效波高的對比檢驗利用MAE 和MRE 參數(shù)進行，檢驗結果見表3和表4。

表3 颶風“卡特里娜”和“瑞塔”有效波高模擬值與實測數(shù)據平均絕對誤差Tab.3 MAE of significant wave height between simulated values and measured data during hurricanes"Rita"and"Katrina"

表4 颶風“卡特里娜”和“瑞塔”有效波高模擬值與實測數(shù)據平均相對誤差Tab.4 MRE of significant wave height between simulated values and measured data during hurricanes"Rita"and"Katrina"

經過計算，在YT96、SWAN和新O02方案中，有效波高模擬值和實測數(shù)據對比統(tǒng)計的MAE 分別為0.40 m、0.50 m 和0.38 m，MRE 分別為9.4%、12.2%和9.0%。

4 結論

本文利用9 次颶風過程中的12 份浮標實測資料，采用 COARE 3.0 算法得到基于 TY01 和 O02 兩種考慮波浪方案下的海表面拖曳系數(shù)，通過二項式擬合得到兩種系數(shù)關于風速的參數(shù)化關系式，其中TY01 方案的系數(shù)結果在高風速下未表現(xiàn)出“飽和”特征，O02方案的結果符合已有的觀測規(guī)律，并且其與南海觀測平臺實測數(shù)據符合良好，因此選用基于O02的新參數(shù)化方案用于數(shù)值模擬研究。

為進一步驗證高風速下海表面拖曳系數(shù)參數(shù)化表達式的準確性，本文選取了颶風“卡特里娜”和“瑞塔”兩個過程進行海浪數(shù)值模擬。從模擬結果與實測數(shù)據的對比中可以發(fā)現(xiàn)，新O02 參數(shù)化海表面拖曳系數(shù)的模擬結果與NDBC 浮標實測數(shù)據符合更好。新的參數(shù)化關系方案的MRE 與SWAN 模式默認方案相比提高了3.2%，而與YT96 方案相比提高了0.4%，所以新的參數(shù)化關系方案能更準確地描述海表面拖曳系數(shù)在高風速下的變化，可以更好地用來進行高風速條件下的海浪數(shù)值模擬。