尹 超 , 黃海軍 王道儒, 謝 琳
(1. 中國科學(xué)院海洋研究所, 山東 青島 266071; 2. 中國科學(xué)院大學(xué), 北京 100049; 3. 海南省海洋與漁業(yè)科學(xué)院, 海南 ???571126)
臺風(fēng)是對海南影響較大的氣象災(zāi)害, 特別是水產(chǎn)養(yǎng)殖行業(yè)。海南臨高縣后水灣網(wǎng)箱養(yǎng)殖密度和產(chǎn)值均較高, 是臨高縣支柱產(chǎn)業(yè), 占全縣GDP一半以上。2011年, 遭受臺風(fēng)“納沙”(NESAT)和“尼格” (NALGAE)襲擊, 受損2050口網(wǎng)箱, 占全省的78.8%, 直接經(jīng)濟損失7.63億元; 2014年, 遭受臺風(fēng)“威馬遜”(RAMMASUN)及“海鷗”(KALMAEGI)襲擊, 受損1 383口網(wǎng)箱, 占全省的45%。到了2018年, 海南省開始推進(jìn)漁業(yè)產(chǎn)業(yè)結(jié)構(gòu)調(diào)整, 在陵水縣新建3個大型智能化養(yǎng)殖網(wǎng)箱, 從水深20米以淺近岸區(qū)域拓展到50米以深的深遠(yuǎn)海區(qū)域, 從養(yǎng)殖裝備直徑20~ 30米量級拓展到近100米量級。隨著養(yǎng)殖產(chǎn)業(yè)升級, 為了維護(hù)生命和財產(chǎn)安全, 需要加強近岸海域的海浪預(yù)報系統(tǒng)的建設(shè)以提供科學(xué)依據(jù), 提高海洋防災(zāi)減災(zāi)的能力。
影響海浪預(yù)報的最重要因素是風(fēng)場預(yù)報的準(zhǔn)確性[1]。目前主要的兩大氣象業(yè)務(wù)預(yù)報系統(tǒng)分別是歐洲中期天氣預(yù)報中心的全球模式ECMWF(European Centre for Medium-range Weather Forecasts), 和美國全球預(yù)報系統(tǒng)GFS(Global Forecast System)。GFS是由美國國家海洋和大氣管理局NOAA(National Oceanic and Atmospheric Administration)推出的全球數(shù)值天氣預(yù)報計算模式, 該模式每天分別于世界時00時、06時、12時和18時發(fā)布四次預(yù)報結(jié)果, 預(yù)報時效長達(dá)16天。從2015年開始空間分辨率已經(jīng)提高到了0.25°, 時間分辨率隨著預(yù)報時效增加而降低, 5天內(nèi)分辨率為1小時, 5天至10天分辨率為3小時, 10天至16天分辨率降為12小時。
目前國內(nèi)外的業(yè)務(wù)化海浪數(shù)值預(yù)報主要基于第3代海浪模式, 廣泛應(yīng)用的主要有WAM(Wave Modeling), WAVEWATCH Ⅲ和SWAN(Simulating Waves Nearshore)等。美國氣象局采用WAVEWATCH Ⅲ, 歐洲氣象中心采用WAM, 英國氣象局采用WAM和SWAN等。我國的國家海洋環(huán)境預(yù)報中心用WAM和WAVEWATCH Ⅲ模型進(jìn)行大洋的海浪預(yù)報, 用WAVEWATCH Ⅲ和SWAN嵌套模型進(jìn)行中國海的海浪預(yù)報[2]。通常, WAVEWATCH Ⅲ模式適用于深水大洋波浪計算, SWAN模式更適和近岸波浪計算, WAVEWATCH Ⅲ和SWAN模式的嵌套常用于跨尺度的波浪計算[3]。
臺風(fēng)浪往往伴有風(fēng)暴潮, 風(fēng)暴潮會引起海岸增水[4,5], 水位的變化直接影響近岸波高的大小[6], 而波浪的輻射應(yīng)力的又會引起水位的改變[7]。潮汐引起的水位變化也會影響風(fēng)暴潮和海浪的增長。網(wǎng)箱養(yǎng)殖區(qū)往往設(shè)置在半封閉海灣中, 以躲避強風(fēng)浪的直接作用。但其幾何特征有利于風(fēng)暴潮增水, 以至于增大了臺風(fēng)浪風(fēng)險。所以, 近岸的風(fēng)浪潮流動力過程存在著強相互作用, 而浪潮流耦合模型適用于近岸動力環(huán)境問題的解決[8]。Dietrich耦合了ADCIRC(Advanced circulation model)和SWAN模式研究了颶風(fēng)Gustav (2008)過程的臺風(fēng)浪和風(fēng)暴潮特征[9]; 馮興如等基于ADCIRC+SWAN模型研究了浙江和福建海域臺風(fēng)浪變化特征[10]。
本研究通過海南島近岸養(yǎng)殖區(qū)臺風(fēng)浪精細(xì)化預(yù)報系統(tǒng)的建設(shè), 提高網(wǎng)箱養(yǎng)殖環(huán)境監(jiān)測預(yù)警水平。為漁業(yè)工作者投喂作業(yè)和維護(hù)設(shè)備提供合理的時間窗口, 并降低其勞動風(fēng)險; 為養(yǎng)殖企業(yè)的生產(chǎn)決策和地方政府的應(yīng)急管理服務(wù)提供科學(xué)合理的依據(jù)。
ADCIRC模式由Luettich和Westerink于1992年開發(fā)并不斷完善發(fā)展[11]。該模式采用了有限元方法, 在二維條件下應(yīng)用了深度平均正壓和淺水方程, 通過求解廣義波連續(xù)性方程(Generalized Wave Continuity Equation, GWCE), 水深和動量方程可以表述為:
式中: U和V是垂向積分平均流速; f為科氏力參數(shù); ζ為海表面自由高度; Ps為海水自由表面大氣壓; ρ0是海水密度; η為牛頓引潮勢; τs,x, τs,y為海表面應(yīng)力; H=ζ+h為水深加自由表面高度; τb,x, τb,y為底切應(yīng)力; g為重力加速度; M, D和B分別代表了側(cè)向應(yīng)力梯度、動量耗散和斜壓梯度。為了滿足CFL條件, 時間步長取為1 s。ADCIRC模式帶有干濕網(wǎng)格判斷, 最小水深設(shè)為0.1 m。
SWAN模式包含了較多的近岸波浪物理過程, 例如折繞射、淺化、白冠破碎、底摩擦、波-波相互作用等。笛卡兒坐標(biāo)下控制方程表述為:
在耦合過程中, ADCIRC和SWAN使用同一套非結(jié)構(gòu)網(wǎng)格。所以在數(shù)據(jù)傳遞過程中不需要再調(diào)用計算資源進(jìn)行插值和數(shù)據(jù)分配, 保證了計算效率和穩(wěn)定性。因為SWAN模式的有更大的時間步長, SWAN和ADCIRC模式的數(shù)據(jù)交換時間間隔設(shè)定為SWAN時間步長[9]。由ADCIRC模式計算的風(fēng)速、海表面水位和流速每10分鐘傳遞給SWAN模式, 而SWAN模式將輻射應(yīng)力梯度傳回給ADCIRC模式。
水深地形采用美國國家海洋大氣管理局NOAA提供的ETOPO1數(shù)據(jù)集, 該數(shù)據(jù)為柵格化數(shù)據(jù), 分辨率為1′×1′。由于ETOPO1對近岸水深的觀測誤差較大, 近岸水深數(shù)據(jù)采用中華人民共和國海軍航海保證司令部出版的海圖, 分辨率從1︰150 000加密到1︰50 000再到重點海灣的1︰25 000(見圖1)。
近岸養(yǎng)殖區(qū)臺風(fēng)浪耦合數(shù)值預(yù)報模型計算區(qū)域范圍是105.5°—119.7°E和 10.8°—24.7°N, 包括海南島和南海大部分的海域。由于海南島海岸島嶼眾多, 為了精確描述岸線特征, 達(dá)到精細(xì)化預(yù)報要求, 我們從Google Earth提取了整個海南島的岸線, 使沿岸區(qū)域分辨率達(dá)到了100 m。網(wǎng)格分辨率從南海深海到陸架到近岸再到重點海灣逐漸提高(見圖2)。模型在開邊界處加入了俄勒岡州立大學(xué)開發(fā)的TPXO模型(http: //volkov.oce.orst.edu/tides/)的8個主要天文分潮(P1、Q1、K1、O1、K2、N2、M2和S2)。
圖1 海南島近岸地形分布圖 Fig. 1 Topography around Hainan Island
圖2 模式計算域網(wǎng)格 Fig. 2 Unstructured grid of the model domain
2014年09號臺風(fēng)“威馬遜”生成于西北太平洋, 2014年7月11日由熱帶低壓升級為熱帶風(fēng)暴, 7月14日升級為強熱帶風(fēng)暴, 并于下午5時升級為臺風(fēng)。7月15日“威馬遜”在菲律賓阿爾拜省拉普拉普登陸, 7月16日臺風(fēng)向西北偏西移動, 進(jìn)入南海海域。17日, 進(jìn)入南海北部, 形成為強臺風(fēng), 其路徑見圖3。臺風(fēng)“威馬遜”于7月18日15時30分在海南省文昌市翁田鎮(zhèn)登陸, 中心最大風(fēng)力達(dá)到了60 m/s, 中心氣壓為910 hPa, 10級風(fēng)圈半徑180 km, 7級風(fēng)圈半徑300 km。于19日7時10分在廣西防城港市再次登陸, 之后強度逐漸減弱, 于20日5時進(jìn)入云南境內(nèi), 并減弱為熱帶低壓[13]。超強臺風(fēng)“威馬遜”對海南、廣東和廣西沿海地區(qū)造成的直接經(jīng)濟損失達(dá)265.5億元。
圖3 海洋觀測站(藍(lán)點)和臺風(fēng)“威馬遜”路徑(紅色實線) Fig. 3 Oceanographic stations (blue dots) and track of the RAMMASUN cyclone (red line)
潮位和波浪觀測數(shù)據(jù)由國家海洋局南海分局??谥行恼咎峁? 站位如圖3所示。潮位數(shù)據(jù)每小時一次, 波浪數(shù)據(jù)每3小時一次, 當(dāng)臺風(fēng)強度增強則加密到1小時一次。??谡静ɡ藬?shù)據(jù)因為儀器原因, 在接近波高最大值時停止工作, 缺失數(shù)據(jù)。
強迫場方面使用CFSv2(Climate Forecast System Version 2)再分析數(shù)據(jù), 包括海面10 m風(fēng)場和海表面氣壓場。CFSv2系統(tǒng)由美國環(huán)境預(yù)報中心NCEP(National Centers for Environmental Prediction)開發(fā), 它是全耦合模型包括了大氣、海洋、陸地和海冰模塊, 于2011年投入業(yè)務(wù)化運行。該強迫場時間分辨率為1小時, 空間分辨率為0.2°×0.2°。
模型計算時間從2014年7月9日18時至20日6時。模式的輸出變量設(shè)定含有效波高、譜峰周期、平均周期、波向, 輸出頻率為每小時一次。為檢驗?zāi)J降慕Y(jié)果, 將計算的潮位和有效波高結(jié)果與臺風(fēng)“威馬遜”期間觀測波高進(jìn)行對比, 對比結(jié)果見圖4。由近岸觀測站潮位對比(圖4a和b)可知, 模擬的潮位在東方站和實測值基本一致, 海口站略低且峰值提前, 原因可能是再分析風(fēng)圈結(jié)構(gòu)和實際風(fēng)圈有差別。由近岸觀測站有效波高對比圖(圖4c和d) 可知, 有效波高增長趨勢一致性較好, 最高值基本在同一時刻。海口站和東方站模擬有效波高值均略大于實測值。
2018年7月19日在海南島近岸生成了熱帶風(fēng)暴, 其過程持續(xù)到7月25日。圖5給出的是GFS預(yù)報的熱帶風(fēng)暴過程的風(fēng)速場和氣壓場分布圖, 預(yù)報時間均為2018年7月22日18時(UTC)。左圖為2018年7月21日0時(UTC)發(fā)布的預(yù)報數(shù)據(jù), 該組數(shù)據(jù)的預(yù)報時效約48小時。右圖為2018年7月22日0時(UTC)的發(fā)布數(shù)據(jù), 該組數(shù)據(jù)的預(yù)報時效約24小時。圖5a 和c顯示, 熱帶風(fēng)暴中心大體位于瓊州海峽中部, 海南島東部風(fēng)速較大, 能達(dá)到18 m/s。圖5b和d顯示, 預(yù)報的熱帶風(fēng)暴中心向西偏移, 大體位于海南島西北部海域, 海南島南部風(fēng)速較大, 達(dá)到了16 m/s, 比21日0時發(fā)布的預(yù)報數(shù)據(jù)極值偏小。
圖4 1409號臺風(fēng)“威馬遜”期間沿岸各站潮位及有效波高觀測值和模擬值 Fig. 4 Computed results of sea surface elevations and significant wave heights (red lines) vs. observations (black dots) during Typhoon RAMMASUN
圖5 GFS預(yù)報的2018年7月22日18時(UTC)海南島附近風(fēng)速場和氣壓場 Fig. 5 Wind and surface air pressure fields predicted by GFS at 18: 00 UTC on July 22, 2018 near Hainan Island
圖6顯示的是模擬的2018年7月22日18時(UTC)海南島附近波浪場, 圖6a采用2018年7月21日0時(UTC)發(fā)布的預(yù)報風(fēng)速場和氣壓場, 而圖6b采用2018年7月22日0時(UTC)發(fā)布的預(yù)報風(fēng)速場和氣壓場。有效波高的分布與圖5風(fēng)速強度分布基本吻合。圖6a波浪高值區(qū)域在海南島東部海域, 最大值達(dá)到了5 m。圖6b高值區(qū)域在海南島南部海域, 相比于圖6a最大值減少了大概1 m。
將48小時、24小時預(yù)報值和海南陵水灣波浪觀測數(shù)據(jù)進(jìn)行對比(見圖7)。時間范圍選取7月20日16時-24日2時, 其中紅線代表48小時預(yù)報值、藍(lán)線代表24小時預(yù)報值、黑色點代表實測值。48小時預(yù)報最大有效波高為3.2 m左右, 24小時預(yù)報最大有效波高為2.15 m左右, 二者在近岸的最大有效波高差值接近1.0 m。在7月23日4時之前, 48小時預(yù)報的有效波高較24小時預(yù)報值偏大, 之后較24小時預(yù)報值偏小且和實測值更為吻合。48小時預(yù)報結(jié)果存在極值偏高的現(xiàn)象, 且波高增長時間及區(qū)間提前。24小時預(yù)報有效波高在波浪增長、整體趨勢上和實測值更為吻合, 7月23日9時之后的數(shù)值較實測值偏大。對比圖5可以看出, 48小時預(yù)報最大風(fēng)速中心區(qū)域在海南島東海岸, 陵水波浪觀測站受其影響較大, 風(fēng)速比24小時預(yù)報值偏大, 導(dǎo)致波浪計算結(jié)果也較大。
圖6 模擬的2018年7月22日18時(UTC)海南島附近波浪場 Fig. 6 Simulated significant wave heights at 18: 00 UTC on July 22, 2018 near Hainan Island
圖7 熱帶風(fēng)暴期間48小時和24小時預(yù)報的有效波高 Fig. 7 The 48 h and 24 h forecasts of significant wave heights during the tropical storm
最后對有效波高進(jìn)行定量檢驗評估, 評估參數(shù)包括平均絕對誤差(MAE)、平均相對誤差(MRE)、均方根誤差(RMSE)和皮爾遜相關(guān)系數(shù)(R)[14-16]。其中, 平均絕對誤差、平均相對誤差和均方根誤差體現(xiàn)了觀測值同實測值的平均偏離情況, 皮爾遜相關(guān)系數(shù)則能反映觀測值同實測值變化趨勢的密切程度。檢驗評估的計算公式如下所示:
有效波高的預(yù)報值與觀測值的檢驗統(tǒng)計結(jié)果如表1所示。24小時和48小時預(yù)報的MAE均小于0.3 m, 24小時的MAE為0.22 m, 相對于48小時預(yù)報值提高了23%。48小時預(yù)報值MRE為20.75%, 而24小時預(yù)報值MRE為17.0%。48 小時預(yù)報值RMSE為0.439 m, 24小時預(yù)報值的RMSE提高到了0.272 m, 提高幅度較大。24小時和48小時的相關(guān)系數(shù)R均大于0.8, 24小時預(yù)報值略有提高達(dá)到了0.87。四項統(tǒng)計指標(biāo)24小時預(yù)報值均好于48小時的結(jié)果, 其中RMSE和MAE指標(biāo)提高較大。
表1 預(yù)報值與觀測值的統(tǒng)計檢驗結(jié)果 Tab. 1 Statistical analysis of forecast results and observations
本文根據(jù)近岸網(wǎng)箱養(yǎng)殖區(qū)臺風(fēng)浪精細(xì)化預(yù)報的工作需求, 基于非結(jié)構(gòu)網(wǎng)格的ADCIRC+SWAN耦合模型, 采用GFS風(fēng)場和氣壓場數(shù)據(jù), 建立了適用于海南近岸的海浪預(yù)報系統(tǒng)。利用近岸觀測資料對2014年7月“威馬遜”臺風(fēng)過程進(jìn)行后報驗證, 在此基礎(chǔ)上對2018年7月海南島近岸海域熱帶風(fēng)暴過程進(jìn)行海浪預(yù)報及結(jié)果統(tǒng)計分析檢驗。主要結(jié)論如下:
(1) 采用高分辨率CFSv2風(fēng)場和氣壓場作為模式輸入, 在最大波高與波浪增長衰減過程兩個方面, 計算的結(jié)果和近岸觀測一致性較好。模擬的有效波高值略大于實測值。
(2) GFS預(yù)報的風(fēng)場與氣壓場每6小時不斷更新, 熱帶氣旋的路徑、中心氣壓和中心強度均不斷變化。預(yù)報時效對海浪預(yù)報有直接影響。
(3) 熱帶風(fēng)暴過程的預(yù)報結(jié)果對比發(fā)現(xiàn), 48小時預(yù)報的最大波高比實測偏大較多, 24小時預(yù)報的波浪增長衰減過程和實測值更為吻合。24小時預(yù)報的平均相對誤差達(dá)到了17%, 均方根誤差為0.272 m, 24小時和48小時預(yù)報的相關(guān)系數(shù)均大于0.8。
通過對預(yù)報結(jié)果的統(tǒng)計檢驗, 該模型的預(yù)報精度可以滿足近岸養(yǎng)殖區(qū)臺風(fēng)浪預(yù)報業(yè)務(wù)的需求??梢詾轲B(yǎng)殖企業(yè)的生產(chǎn)決策和地方政府的應(yīng)急管理服務(wù)提供科學(xué)合理的依據(jù)。