林澤輝

(惠州市惠陽區(qū)永良堤圍管理所，廣東 惠州 516007)

我國東南沿海地區(qū)氣象水文條件復(fù)雜，水文方面，海水水位受天體引潮力作用而周期性漲落。氣象方面，西太平洋熱帶地區(qū)氣壓受到擾動作用成渦，并不斷發(fā)展移動，這些氣旋登陸后會對人類社會造成巨大的破壞。另一方面，臺風(fēng)所造成的大氣壓強(qiáng)場變化也會導(dǎo)致海面異常起伏，產(chǎn)生風(fēng)暴潮，當(dāng)風(fēng)暴潮恰好疊加上天文大潮，海平面會巨幅上升，這對海岸防護(hù)、河口擋潮等水利工程帶來了巨大挑戰(zhàn)。因此研究臺風(fēng)過程前后的波浪以及增水具有重要意義。文章利用…模型對深圳地區(qū)水域在臺風(fēng)作用下的水動力特征進(jìn)行了數(shù)值模擬研究。

1 研究方法

1.1 區(qū)域概況

深圳市海岸線長約250km，東起大亞灣，西至珠江口，增水>50cm的風(fēng)暴潮過程年平均出現(xiàn)2.1次。沿岸分布有機(jī)場、重點(diǎn)港口、大型住宅區(qū)、核電站等大型企業(yè)，人口密度極高。

深圳市年平均風(fēng)速2.6m/s，冬季各月風(fēng)速較大為3.0m/s左右，夏季各月較小約2.0m/s。大風(fēng)日數(shù)(風(fēng)速≥17m/s，即風(fēng)力相當(dāng)于8級或8級以上)年平均為7.3d。全年各月均出現(xiàn)大風(fēng)，其中以7-9月居多，共占全年大風(fēng)日數(shù)的44%，其次為6月和10月，由于受臺風(fēng)或熱帶低壓影響，常出現(xiàn)大風(fēng)，因此，大風(fēng)出現(xiàn)多的月份，也是臺風(fēng)季節(jié)。

1.2 數(shù)據(jù)資料

影響風(fēng)暴增水量級的因素有很多，包括臺風(fēng)氣壓場分布、臺風(fēng)風(fēng)場的方向與大小以及近岸地形變化因素。在從我國東南沿海登陸的臺風(fēng)里，影響到陽江至汕尾海域的臺風(fēng)均可引起深圳海域的風(fēng)暴潮增水。目前在深圳海域，國家海洋局共設(shè)有兩個長期驗潮站，分別是赤灣和鹽田站，緊鄰深圳大亞灣岸線的有位于惠州市的惠州海洋站，各站位置見圖1。收集上述站點(diǎn)自建站以來臺風(fēng)期間的潮位資料統(tǒng)計歷史風(fēng)暴潮增水情況。廣東省水文局1964年在深圳市赤灣港內(nèi)建設(shè)有長期水文觀測站，1986年國家海洋局南海分局在相距不到1km處建設(shè)了赤灣海洋站，并于1986年開始投入業(yè)務(wù)化觀測。鹽田海洋站位于大鵬灣鹽田港海域，由國家海洋局南海分局建設(shè)和管理，2003年投入使用?；葜莺Ｑ笳疚挥诖髞啚碁稠斘鱾?cè)荃灣港區(qū)，由國家海洋局南海分局建設(shè)和管理，2006年建成投入使用，他們的分布位置如圖1所示。

圖1 資料站位位置示意圖

1.3 研究方法

早期人們主要通過經(jīng)驗總結(jié)歸納的方法研究風(fēng)暴潮，觀察臺風(fēng)天時海塘外潮位變化研究風(fēng)暴潮的成因。20世紀(jì)中期，隨著計算機(jī)技術(shù)的發(fā)展，人們開始利用數(shù)值模擬研究風(fēng)暴潮。1956年W·Hansen首次利用計算機(jī)對歐洲北海的最大風(fēng)暴潮增水進(jìn)行了研究。隨著數(shù)值模擬技術(shù)的進(jìn)步，Jelesnianski建立了研究風(fēng)暴潮的SPLASH模式以及基于二維水動力的SLOSH模式，在不斷地摸索研究后，該模式的預(yù)報精度不斷提高。近年來，不斷發(fā)展的模型計算實現(xiàn)了對復(fù)雜地形影響的考慮，提高了計算的準(zhǔn)確性[1]。其次，有關(guān)風(fēng)暴增水與天文潮的非線性擬合機(jī)制也變得更加清晰，這也提高了風(fēng)暴潮數(shù)學(xué)模型的準(zhǔn)確性。我國自20世紀(jì)60年代也開展了關(guān)于風(fēng)暴潮的相關(guān)研究。馮士筰所著的《風(fēng)暴潮導(dǎo)論》系統(tǒng)闡述了風(fēng)暴潮的概念，形成機(jī)制，并介紹了早期的風(fēng)暴潮數(shù)值預(yù)報方法。趙長進(jìn)等利用ADCIRC潮流與SWAN波浪耦合模式對長江口附近海域的風(fēng)暴潮過程進(jìn)行了模擬。黃本勝等用SELFE模型構(gòu)建了南海—珠江河口雙重嵌套模型，研究了“山竹”臺風(fēng)過程中珠江河口的風(fēng)暴增水與波浪增水分布。

在波浪模擬方面，目前工程界使用較為廣泛的數(shù)學(xué)模型有兩類，一類是波浪流體動力學(xué)模型，另一類是動譜密度方程。緩坡方程模型與Boussinesq方程模型是主要的波浪流體力學(xué)模型。經(jīng)典的Boussinesq僅僅適用于淺水，并不適用于水深變化大的水域和不規(guī)則波的情況。因此在工程中適用范圍不廣。動譜密度方程的基礎(chǔ)是能量平衡理論，在模擬過程中可提高海岸波浪場的動力作用過程的計算效率。該方程以動譜密度來描述波浪變形，考慮了對波-波非線性相互作用、波流相互作用、波浪破碎和底部摩擦引起的能量耗散等過程，通過將源項添加到方程中，可以全面考慮動力作用的復(fù)雜過程，使其能完成大尺度區(qū)域的波浪產(chǎn)生、成長、傳播和變形的計算過程。其中對物理源項的處理和研究是波作用量守恒的關(guān)鍵，波作用量守恒方程發(fā)展到現(xiàn)在已經(jīng)到了第三代，逐漸開發(fā)出WAVEWATCH、SWAN、MIKE21SW、TOMAWAC等模型。

1.4 模型設(shè)置

1.4.1 模型介紹

文章采用的是DELFT3D-FLOW與SWAN的耦合模型，DELFT3D是一款由荷蘭DELTRA開發(fā)的開源水動力模擬軟件，在風(fēng)暴潮模擬方面已經(jīng)廣泛應(yīng)用。

Delft3D-Flow是基于淺水方程和Boussinesq假設(shè)，采用有限差分法求解斜壓Navier-Stokes方程的模型，模型控制方程組分別為：

質(zhì)量守恒方程：

(1)

式中：ζ為水位(總水深)；d為到參考平面的當(dāng)前水深(凈水深)；U、V-x、y為方向上的速度分量；Q為單位面積上質(zhì)量強(qiáng)度。

動量守恒方程：

(2)

(3)

式中：f為科氏力參數(shù)；g為重力加速度；νh為動態(tài)水平渦流黏度；τsx,τsy為x、y方向上的海平面風(fēng)壓分量；τbx,τby為x、y方向上的底部的剪應(yīng)力分量；ρ0為參考密度；ρ'為不規(guī)則密度。

SWAN模型的控制方程為：

(4)

式中：Cx，Cy，Cσ，Cθ為二維動譜密度在x、y、σ、θ為空間上的傳播速度；σ為相對頻率；θ為波向；S為能量源項。

1.4.2 模型設(shè)置

1)模型范圍：

文章所采用的模型囊括整個珠江三角洲區(qū)域，見圖2。

圖2 模型區(qū)域示意圖

2)地形水深：

模擬范圍內(nèi)的地形水深采用的是中華人民共和國海事局2016年繪制的海圖插值數(shù)據(jù)。將地形水深數(shù)據(jù)轉(zhuǎn)化到統(tǒng)一基面(平均海平面)后，用IDW插值方法插值到網(wǎng)格點(diǎn)上，得到模型的水深地形。

3)開邊界設(shè)置：

開邊界采用天文潮水位，主要考慮M2、N2、S2、K2、K1、O1、P1、Q1等八個天文分潮以及MS4、M4、M6等淺水分潮共11個分潮的作用。分潮的調(diào)和分析常數(shù)由2016-2019年的潮位實測數(shù)據(jù)經(jīng)MATLAB中的T_TIDE工具包調(diào)和分析所得。模型上部采用流量邊界，輸入數(shù)據(jù)為梧州、石角、博羅三個水文站的2016-2019洪季平均流量。

4)底摩擦設(shè)置：

模型根據(jù)水深范圍設(shè)置不同大小的糙率，水深越大糙率越小。

1.4.3 風(fēng)場輸入

文章的風(fēng)場采用美國國家大氣研究中心(NCAR)等機(jī)構(gòu)共同開發(fā)的中尺度大氣預(yù)報系統(tǒng)。WRF模式適用于臺風(fēng)模擬，可以通過嵌套網(wǎng)格運(yùn)算來達(dá)到高精度模擬的目的。本次采用了WSM6微物理方案、KainFritsch卷積云方案、Monin-Obukhov表面層方案、Noah陸面層方案、MYJ大氣邊界層方案、RRTMG輻射方案模擬臺風(fēng)。

2 模型驗證

2.1 風(fēng)場驗證

首先利用實測資料對所建立的風(fēng)暴潮-波浪耦合模型進(jìn)行率定驗證。選取1713號臺風(fēng)(2017年第13號臺風(fēng)—天鴿，其移行路徑如圖3所示)作為模型邊界條件，為了證明WRF模式的準(zhǔn)確性其在風(fēng)暴潮模擬方面的優(yōu)勢，文章也用其他論文常用的經(jīng)驗風(fēng)場模型做了驗證對比。經(jīng)驗風(fēng)場采用的是宮崎正衛(wèi)模型作為移動風(fēng)場，HOLLAND模型作為梯度風(fēng)場，他們的公式分別為：

圖3 “天鴿”臺風(fēng)移行路徑

(5)

(5)

式中：B為HOLLAND模型氣壓參數(shù)；f為科氏力參數(shù)，與模擬區(qū)域緯度有關(guān)；R為最大風(fēng)速半徑，與中心氣壓有關(guān)。

圖4為WRF模型模擬風(fēng)場值以及經(jīng)驗風(fēng)場模擬值與實測風(fēng)場的對比，可以看出，WRF模型模擬的風(fēng)場在峰值以及過程上都比經(jīng)驗風(fēng)場模型更貼近實際，因此利用WRF模型生成的臺風(fēng)風(fēng)場作為風(fēng)暴潮-波浪耦合模型的邊界條件更為合理。

圖4 WRF風(fēng)場及模型風(fēng)場與實測風(fēng)場對比

2.2 風(fēng)暴潮及波浪率定驗證

以1713號臺風(fēng)(2017年第13號臺風(fēng))為邊界條件分別驅(qū)動珠江口風(fēng)暴潮模型以及珠江口波浪-風(fēng)暴潮耦合模型，得到如圖5所示的赤灣站與鹽田站水位驗證曲線，可以看出，考慮波浪后的水位過程曲線更符合實際，最大誤差絕對值為13cm，平均相對誤差為11.3%，標(biāo)準(zhǔn)差為2.94，符合模擬精度要求。

圖5(a) 赤灣站不同模型下水位曲線

圖5(b) 鹽田站不同模型下水位曲線

3 深圳市沿岸漫堤風(fēng)險評估

3.1 可能最大風(fēng)暴潮模型計算

為了尋求可能影響深圳市的最大風(fēng)暴潮對深圳市沿岸的影響，采用概率論法進(jìn)行計算可能最大臺風(fēng)中的近中心最大風(fēng)速Vmax，臺風(fēng)中心氣壓P0。

利用中國氣象局(CMA)1949-2015年熱帶氣旋最佳路徑集(Best-track)數(shù)據(jù)，選取出研究區(qū)域400km直徑范圍內(nèi)對本研究區(qū)域有重要影響的臺風(fēng)，以所選取臺風(fēng)的最小臺風(fēng)中心氣壓值作樣本。采用極值I型分布計算1000年一遇的臺風(fēng)中心氣壓值作為可能最大臺風(fēng)的中心氣壓，計算結(jié)果如圖6(a)所示，即深圳區(qū)域1000a遇的臺風(fēng)中心氣壓P0=880hPa。

圖6(a) 影響深圳的臺風(fēng)中心最低氣壓極值I型分布

圖6(b) 影響深圳的臺風(fēng)中心最大風(fēng)速極值I型分布

利用中國氣象局(CMA)1949-2015年熱帶氣旋最佳路徑集(Best-track)數(shù)據(jù)，取深圳市海灣區(qū)域400km直徑范圍內(nèi)歷年路經(jīng)本海區(qū)的臺風(fēng)最大風(fēng)速值作樣本，對于沒有臺風(fēng)經(jīng)過研究區(qū)域的年份，該年Vmax取為臺風(fēng)Vmax系列中的平均值。采用極值I型分布計算1000a一遇的最大風(fēng)速值作為可能最大臺風(fēng)的最大風(fēng)速值，計算結(jié)果如圖6(b)所示，深圳區(qū)域1000a一遇的臺風(fēng)中心最大風(fēng)速Vmax=85m/s。另外，臺風(fēng)路徑選取為東南向西北由大鵬灣登陸的路徑。

輸出堤前代表點(diǎn)(具體見圖7)的最高潮位與堤頂高程進(jìn)行比對，了解沿岸淹沒的情況，由于模式計算結(jié)果是相對平均海平面的，為便于比對，需要根據(jù)當(dāng)?shù)仄骄Ｆ矫媾c85高程的關(guān)系，將計算結(jié)果轉(zhuǎn)成以85黃?；鏋槠鹚忝?，其中西段提防采用赤灣站高程關(guān)系，即計算結(jié)果加上0.54m；東段提防采用鹽田站高程關(guān)系進(jìn)行轉(zhuǎn)換，即計算結(jié)果加上0.66m。表1列出堤前統(tǒng)計的最高潮位及堤頂高程。從代表點(diǎn)的最高潮位和海堤高程比較來看，在可能最大熱帶氣旋影響下，深圳沿岸均可能發(fā)生漫堤。

圖7 堤前代表點(diǎn)位置示意圖

表1 可能最大風(fēng)暴潮漫堤分析 (單位：m，基面：85高程)

續(xù)表1 可能最大風(fēng)暴潮漫堤分析 (單位：m，基面：85高程)

3.2 海堤防御標(biāo)準(zhǔn)

對于大部分沿海地區(qū)來說，1000a一遇的標(biāo)準(zhǔn)是過高的。但是深圳作為特大型城市，其風(fēng)暴潮安全具有重要意義，因此考慮1000a一遇情形下的風(fēng)暴潮安全是有必要的，在此情形下得出的淹沒情況有利于應(yīng)急管理部門采取相應(yīng)的防范措施。

4 結(jié) 論

1)利用WRF中長尺度氣象預(yù)報模式作為波浪-風(fēng)暴潮耦合模型的風(fēng)場邊界條件輸入，提高了風(fēng)暴潮模型的精確度

2)在珠江口波浪-風(fēng)暴潮模型驗證實測資料準(zhǔn)確的基礎(chǔ)上，利用極值I型分布確定1000a一遇臺風(fēng)參數(shù)對影響深圳市海岸堤防未來發(fā)生的最大可能風(fēng)暴潮進(jìn)行計算分析。文章研究成果可為深圳市風(fēng)暴潮安全防范提供重要技術(shù)支撐，同時也為海洋水利防災(zāi)應(yīng)急部門提供決策參考。