WAVEWATCH Ⅲ不同海冰源項的海浪模擬效果對比

苗琪，徐福敏，俞茂玲

( 1. 中國電建集團 華東勘測設計研究院有限公司，浙江 杭州 311122；2. 河海大學 港口海岸與近海工程學院，江蘇 南京210098)

1 引言

20 世紀50年代以來，作為北冰洋一部分的波弗特海（Beaufort Sea）近海大陸架被發現具有豐富的石油、天然氣儲量。隨著全球變暖，北極海冰儲量大量減少[1–3]，北極航道潛力巨大。北極地區開采所得的礦產利用北極航道經由波弗特海運往世界各地，較傳統方式可節約近一半時間，因此近年來隨著氣候變暖和亞洲天然氣需求激增，各國均對該海域的極地研究展現出極大關注。對北極風暴下波弗特海?馬更些河（Mackenzie River）河口水域的海浪特性進行分析研究，是油氣開采以及北極航道開發中的關鍵一環。

風暴下短期內海冰密集度、海冰厚度、冰緣線位置均劇烈變化的動態海冰，往往會給北極地區海浪的精確模擬帶來巨大困難。國外早期針對波弗特海海域的波浪模擬及海浪特性研究中[4–7]，在海冰的處理上多采用的是未經可靠驗證的模型或經驗公式，或著眼于夏季無海冰覆蓋的近岸區域，對海冰?海浪間的相互作用考慮較為粗糙，且關于北極風暴過程中海冰的短期運動變化、不同風暴下的海浪分布特征的研究依然較少，因此對海冰影響下波弗特海至馬更些河河口海域的海浪進行精確數值模擬，以及對如何在海浪數值模擬中正確計量海冰的影響開展研究具有一定的研究意義。

海浪模式是現今較為成熟的海浪計算及預報手段。第三代海浪模式WAVEWATCH Ⅲ在控制方程、程序結構、數值方法等很多方面都作出改進，因此近年來在海浪業務預報上得到了廣泛的應用[8–11]。在版本V5.16 中，WAVEWATCH Ⅲ模式已具備多種海冰對海浪的能量耗散源項(IC1?IC4，IS1?IS2)，能在極地地區海浪數值模擬中相對精確地考慮海冰的影響，因此近年來，一些學者開始利用WAVEWATCH Ⅲ模式開展波弗特海海域的海冰?海浪相互作用研究，如Erick 等[12]利用WAVEWATCH Ⅲ模式中的IC3源項對秋季波弗特海的波浪進行后報，并利用后報結果和浮標實測數據分析了不同冰種對海冰的能量耗散能力；Cheng 等[13]利用實測冰情、浪情率定了WAVEWATCH Ⅲ模式中IC3 源項的參數，并對波弗特海海域的海浪進行模擬。這些研究多針對特定海冰源項，而對于WAVEWATCH Ⅲ模式中的多種海冰源項缺乏詳細而完整的評述。據此，本文利用WAVEWATCH Ⅲ模式中的不同海冰源項分別對海冰影響下波弗特海夏季開冰期的暴風浪進行數值模擬，以期為極地地區海洋工程建設、海洋防災減災等提供一定的參考。

2 WAVEWATCH Ⅲ模式介紹

2.1 控制方程

WAVEWATCH Ⅲ模式中，球坐標系下關于動譜密度的控制方程通常表示為

式中，N為動譜密度；t為時間； ?x為水平拉普拉斯算子；k為波數矢量；θ為波向；S為源函數；σ為圓頻率；Cg為群速度；U為環境流速；d為水深；s為與波向相同的坐標；m為與s垂直的另一坐標。

2.2 海冰源項

WAVEWATCH Ⅲ模式中海冰損耗作用最傳統的處理方式是采用經驗方法[14]，近似地將海冰密集度大于某臨界值（由用戶指定）的網格點視作陸地，小于某臨界值的視為不受海冰影響的開闊水域，兩者之間進行線性插值，并通過修改控制方程使波能在每個海冰網格點處進行百分比傳遞，該方法記為IC0。WAVEWATCHⅢ V5.16 版本中引入了4 種海冰損耗源項，分別以IC1?IC4 來表示：第一種海冰損耗源項（IC1）對所有的波浪成分都采用相同的耗散率，是一種不依賴于頻率的簡單損耗過程，在計算時需要輸入海冰的密集度以及損耗率；第二種海冰損耗源項（IC2）將海冰視作一個連續的彈性薄盤，能量損耗由冰下邊界層的摩擦產生[15–17]，在計算時需要輸入海冰密集度、海冰有效剪切模量、黏度；第三種海冰損耗源項（IC3）認為海冰是一個黏滯彈性層[18]，對損耗的計算需要輸入海冰的密集度、厚度、有效黏度系數、密度和有效剪切模量；第四種海冰損耗源項（IC4）對能量損耗給出了幾組簡單的依賴于頻率的經驗或參數化的公式[19–24]，從而使得源項能夠簡單靈活卻有效地再現冰?浪相互作用過程中海浪能量的衰減。

散射源項有兩種，分別以IS1 和IS2 表示：第一種散射源項（IS1）是比較簡單的海冰對海浪的守恒作用，僅考慮海冰密集度，海冰只散射海浪能量，但是不耗散海浪的能量；第二種散射源項（IS2）依賴于浮冰尺寸，并且考慮了海浪作用下海冰的破碎從而能夠不斷更新最大浮冰尺寸[24]。

3 驗證資料來源及模型設置

3.1 驗證浮標資料

2014年7月 末，陸 緣 海 冰 帶（Marginal Ice Zone，MIZ）項目由美國海軍在波弗特海/楚科奇海開展（http://www.apl.washington.edu/project/project.php?id=miz），范圍為72°～79°N，137°～170°W。項目過程中共投放3 個配備定位系統的SWIFT（Surface Wave Instrument Floats with Tracking）漂移浮標，并于9月底回收。SWIFT10 和SWIFT11 浮標于7月27 日在阿拉斯加北部約100km 的近岸處從RV“Ukpik”號極地考察船上投放，向西北方向漂移。在9月1 日之前兩者漂移路徑大致相同，然后SWIFT10 浮標開始進入高密集度海冰區域，直至9月15 日。SWIFT15 浮標在8月5?17 日期間由R/V “Araon”號極地考察船投放并回收，且全程被海冰環繞。各浮標漂移路徑見圖1。

由于SWIFT15 浮標全程被高密集度海冰包圍，幾乎處于擱淺狀態，大部分時間測得有效波高均為0m，最大有效波高不足0.7m，故本次研究僅選取SWIFT10和SWIFT11 兩個浮標用作后續驗證。

圖1 SWIFT 浮標漂移路徑Fig. 1 The track of SWIFT drifting buoys

3.2 模型設置

2014年8?9月有多場風暴襲擊波弗特海至馬更些河河口水域，且影響范圍較廣。研究區域過大，通常會導致計算效率過于低下，故研究采用自嵌套的方案，即先用分辨率較低的網格對較大的海域進行計算，然后將計算所得波浪譜作為邊界輸入到內層分辨率更高、范圍更小的模型中進行計算，以提高計算效率。嵌套外層包括整個東西伯利亞海、楚科奇海、波弗特海和部分位于西北太平洋的白令海，內層為波弗特海至馬更些河口外水域，據此建立自波弗特海至馬更些河口的WAVEWATCH Ⅲ兩級嵌套模型，由外層的計算獲取傳播至波弗特海的所有低頻涌浪成份，輸入到第二層中，然后在存在大量海冰的波弗特海進行風?冰?浪相互作用模擬研究，并探究WAVEWATCH Ⅲ對海冰影響下的海浪場模擬的準確性。

嵌套計算時，中低頻截斷頻率取為0.04118 Hz，高頻截斷頻率取為0.7904 Hz，以對數分布fn+1=1.1fn劃分為32 個頻率網格，外層譜方向劃分為30 個波向，分辨率為12°，內層譜方向劃分為36 個波向，分辨率為10°。

3.2.1 計算區域和水深

本文建立的WAVEWATCH Ⅲ兩級自嵌套海浪模型網格采用球面坐標，各層范圍為：（1）嵌套Ⅰ區：整個東西伯利亞海、波弗特海和部分位于西北太平洋的白令海，范圍為58°～78°N，140°E～110°W，空間分辨率為0.25°×0.25°，網格節點數為441×81；（2）嵌套Ⅱ區：波弗特海至馬更些河河口外水域，范圍為65°～75°N，125°～175°W，空間分辨率為5′×5′，網格節點數為601×121（圖2）。

圖2 兩級嵌套模型研究區域示意圖Fig. 2 Two-level nested computational domain

模型水深數據采用美國國家海洋與大氣管理局（NOAA）的ETOPO1 全球地形數據集，其空間分辨率為1′×1′，在各級嵌套區域分別插值至各網格點處。

3.2.2 驅動風場

模型風場數據采用CCMP（Cross-Calibrated Multi-Platform）V2.0 風場，它是CCMP 風場的升級版本，在繼承CCMP 風場高時間、空間分辨率的基礎上（時間分辨率為6 h，空間分辨率為0.25°×0.25°），又將數據覆蓋時間從2011年12月延長至今。對2014年8月1 日00 時 至9月30 日18 時 的CCMP V2.0 風 場 進 行空間線性插值，作為模型驅動風場。

鑒于CCMP V2.0 風場在相關研究中應用暫時較少，將CCMP V2.0 風場的風速與SWIFT 漂移浮標的實測風速進行對比，以驗證其可靠性。將CCMP 風速分別插值至每個時刻兩個浮標所在位置處，對比結果如圖3 所示。可以看出，兩條風速曲線趨勢和變化幾乎完全一致，數據吻合較好，證明研究所用的CCMP V2.0 風場風速與波弗特海的風速吻合較好，能夠較好地反映北極風暴過程。

3.2.3 海冰資料

海冰密集度數據資料來源有兩種。第一種為NOAA 最優插值海面溫度分析數據集(NOAA Optimum Interpolation 1/4 Degree Daily Sea Surface Temperature Analysis, Version 2)，被用于嵌套模型的外層，即從白令海峽至波弗特海的整個海域。該數據集由NOAA/NCDC 提供，空間上覆蓋全球，分辨率為0.25°×0.25°，時間上從1981年9月1 日至今，分辨率為1 d。第二種為MASAM2 逐日4km 北極海冰密集度數據集（MASAM2: Daily 4km Arctic Sea Ice Concentration,Version 1），被用于模型的第二層，即波弗特海至馬更些河河口。該數據集時間覆蓋范圍為2012年7月至今，分辨率為1 d，空間覆蓋范圍為29.08°～90°N，環全球經度，空間分辨率為4km×4km。

海冰厚度數據采用美國國家環境預報中心（NCEP）氣候預報系統2.0 版本（CFSV2）的逐時時間序列產品（NCEP Climate Forecast System Version 2（CFSV2）Selected Hourly Time-Series Products），范圍覆蓋全球，有0.2°×0.2°、0.5°×0.5°、1.0°×1.0°和2.5°×2.5°的空間分辨率可供選擇，時間覆蓋為2011年1月1 日至今，分辨率為6 h。CFSV2 海冰厚度數據相比實測數據往往偏大[25–26]，故在本次模擬中對其乘以0.6 的折減。

3.2.4 參數設置

計算時，對于海冰的散射作用，由于計算海域中浮冰直徑相對波長來說較小（可達2m，但通常小于1m），能帶來的散射作用有限，且WAVEWATCH Ⅲ模式中的海冰散射模塊尚不能很好地體現海冰對海浪能量的衰減作用，故本次計算不考慮海冰的散射作用。對于海冰對海浪能量的衰減作用，分別采用IC0、IC1、IC2、IC34 種機理進行考慮，其中IC0 額外設置海冰損耗上下限為50%～50%及25%～75%兩種方案。其中50%～50%方案為模型的缺省設置方案，而25%～75%方案為一些學者采用過的改進方案。

圖3 CCMP V2.0 風場風速驗證Fig. 3 Validation of CCMP V2.0 wind speed

由于應用于該大范圍海域時，缺乏海冰的剪切模量、黏滯系數等實測屬性參數，故除海冰密集度、厚度由衛星資料提供外，其余海冰相關的模式參數均取為模型缺省值。其他源項設置還包括風能輸入和白浪耗散項選取Ardhuin 等[27]的WAM4 方案，四波非線性相互作用項選取DIA 方法[28]，三波非線性相互作用項選取LTA 方法[29]，底摩擦項選取JONSWAP 底摩擦公式[30]，水深變淺引起的波浪破碎項選取Battjes-Janssen 的方法[31]，參數均為默認值。

4 模擬結果及對比分析

4.1 模擬結果

由于SWIFT 浮標屬于漂移浮標，所處位置隨時間在不斷變化，需利用WAVEWATCHⅢ模式中的沿軌跡輸出功能輸出兩個SWIFT 浮標的軌跡在對應時刻的波浪要素以進行對比驗證。WAVEWATCH Ⅲ模式的沿軌跡輸出功能僅能輸出一維海浪譜，利用公式編程由離散的海浪譜換算有效波高（其中m0為海浪譜的零階譜矩），從而與實測數據進行驗證比較。圖4 為不同海冰方案下兩個SWIFT 浮標處的有效波高模擬結果及浮標實測結果。

由圖4 可見，在9月1 日之前，由于SWIFT10 浮標和SWIFT11 浮標均處于開闊水域中，海冰對兩浮標唯一的影響為當風向背離海冰區時風區長度的不同，但該影響較為微弱，故此時不同的海冰耗散機理對浮標所在位置的海浪模擬并無明顯影響，5 種機理結果并無明顯區別，部分方案結果幾乎完全一致，因此部分方案的結果被遮擋，但各方案結果均較為理想。

9月1?14 日，SWIFT10 浮標實測有效波高開始迅速衰減，從浮標掛載的攝像頭錄像可以發現，此時SWIFT10 浮標已被海冰所包圍（圖5），海面上大面積覆蓋的海冰對波浪能量產生劇烈的衰減作用。這段時間內，IC0（25%）、IC2、IC3 有效波高模擬結果的趨勢與實測數據相同，但數值明顯偏高，IC1 趨勢與實測相同，且數值相對較為接近。IC2 模擬結果與IC1 非常接近，幾乎被完全遮擋（從下文誤差分析中也可以看出兩方案差距較小，SWIFT10 浮標處IC2 誤差略微高于IC1，而在SWIFT11 浮標處略低于IC1）。IC0（50%）在該時段內波要素始終為0。

9月14 日以后，SWIFT10 浮標重新回到開闊水域，IC0（50%）方案仍圍繞觀測值出現劇烈上下波動，而其余方案均較為理想。與SWIFT10 浮標不同，SWIFT11 浮標全程位于冰密集度較低的開闊海域內，從圖4b 可以看出該浮標的有效波高模擬結果趨勢與9月14 日后的SWIFT10 浮標結果類似，即IC0（50%）方案出現異常波動，而其余方案均擬合較好。

IC0（50%）的參數設置，根據其原理，是一種不連續的處理方式，會對海冰密集度大于50%區域的海浪能量產生劇烈衰減，而對海冰密集度小于50%的區域衰減較微弱。根據衛星及遙感數據，浮標所在的波弗特海域在計算時間段內海冰密集度多為40%～70%，這就導致在該密集度范圍下，IC0（50%）方案對密集度變化過于敏感，結果極易產生突變及波動。在9月1?14 日，海冰密集度大于50%，故該時段內IC0（50%）方案模擬結果波高與周期均為0。

圖4 不同海冰方案下兩個SWIFT 浮標有效波高對比Fig. 4 Comparison of simulated significant wave height by two SWIFT bouys of different ice source terms

IC0（25%）方案雖然也是一種不連續的處理方式，會對海冰密集度大于75%區域的海浪能量產生劇烈衰減，而對海冰密集度小于25%的區域衰減較微弱，但其能在兩密集度范圍之間均勻過渡。因而雖然機理相同，但并未出現類似IC0（50%）方案的波高突變。

4.2 不同海冰模式下模擬結果誤差分析

為更直觀地定量比較不同機理對有效波高的模擬水平，確定該海域最佳的海冰機理，首先分別對不同機理的模擬結果和實測數據進行顯著性水平檢驗，假定模擬結果與實測數據之間相關性為0，計算得p值均小于0.05，證明存在相關性。其次，選取均方根誤差（Root Mean Square Error，RMSE）、相關系數（Correlation Coefficient，CC）和絕對平均誤差（Mean Absolute Error，MAE）3 個參數來統計分析實測值和計算值之間的誤差，分別定義為

圖5 9月1?14 日環繞SWIFT10 的海冰Fig. 5 Sea ice pictures taken around SWIFT10 from September 1st to 14th

表1 和表2 分別為SWIFT10 和SWIFT11 浮標有效波高模擬值與浮標值的相關系數、均方根誤差及絕對平均誤差。

由表1 和表2 可以看出，各方案下SWIFT11 浮標處的有效波高與實測數據之間的均方根誤差和絕對平均誤差普遍小于SWIFT10，而相關系數高于SWIFT10，證明WAVEWATCHⅢ模式對有效波高的模擬效果在遠離冰緣、水域相對開闊時比在被海冰控制下更好。此外，從SWIFT10 浮標處的有效波高結果可以看出，在離海冰較近，受其影響較大的水域，IC1 方案相較其他方案而言，模擬結果誤差更小，相關系數更高。

表1 SWIFT10 有效波高模擬值與浮標值的均方根誤差（RMSE）、相關系數（CC）和絕對平均誤差（MAE）Table 1 RMSE, CC and MAE of observed and simulated significant wave heights in SWIFT10

表2 SWIFT11 有效波高模擬值與浮標值的的均方根誤差（RMSE）、相關系數（CC）和絕對平均誤差（MAE）Table 2 RMSE, CC and MAE of observed and simulated significant wave heights in SWIFT11

盡管已有相關研究[25]顯示，在有較為完備的海冰屬性參數時，IC3 源項表現出的對海浪的能量衰減特征更加符合實際情況，但在本次研究中，由于研究區域范圍較大，并無能力獲得除海冰密集度、厚度以外的其他海冰數據，如剪切模量、黏性系數等，而這些時空均變化的海冰數據恰恰是IC3 源項必不可少的輸入參量，均取為缺省值顯然不合理。相反，IC1 源項對海冰數據要求較低，因此在應用于大范圍水域、資料匱乏的前提下，計算結果反而與實測數據更為接近。

本文所采用參數的模擬結果并非最優解，其他海冰?波浪源項在海冰數據充足、參數調試合理的情況下可能會得到更優的結果。但是在以下限定條件下：（1）應用于北極波佛特海這一特定海域的有效波高模擬，（2）缺乏海冰剪切模量、黏性系數等實測參數，（3）模型參數均采用默認值，使用IC1 海冰源項來計量海冰對海浪的衰減作用最為合理，能夠為后續淺水地區的進一步海浪嵌套模擬提供更加精準的邊界條件。

5 結論

建立典型北極風暴下的波弗特海多重嵌套WAVEWATCH Ⅲ海浪模型，利用WAVEWATCH Ⅲ模式中的不同海冰損耗源項設置不同海冰方案，并利用浮標實測數據對不同海冰源項在秋季波弗特海的海浪模擬效果進行對比研究。結論如下：

（1）在受海冰影響較小的開闊水域，各海冰機理的模擬結果區別較小，且均驗證良好；

（2）當海冰實際存在時，若在海浪模擬中不考慮海冰的能量衰減作用，將海冰密集度較高的地方直接視為陸地，則當海冰密集度大于某一閾值時，有效波高模擬結果會出現大量的向0 突變，結果不合理。因此在該種情況下，使用計入海冰對海浪的能量衰減作用的IC1、IC2、IC3 源項更為合理；

（3）在缺乏大范圍海冰實測數據、模型其他參數采用默認值的前提下對波弗特海至馬更些河河口這一特定水域進行海浪模擬時，IC1 海冰源項能夠更為準確地表達該海域特定的冰情、冰況對海浪的能量衰減特征，能夠為后續淺水地區的進一步海浪嵌套模擬提供更加精準的邊界條件。