基于二維淺水方程和SWAN模型的風(fēng)暴潮模擬研究

于英霞，張秀姬，宋懷輝

(1.河南科技大學(xué) 土木工程學(xué)院，河南 洛陽 471023；2.洛陽理工學(xué)院 土木工程學(xué)院，河南 洛陽 471023)

0 引言

隨著全球氣候變化以及海洋開發(fā)利用的進(jìn)一步深入，風(fēng)暴潮災(zāi)害的頻發(fā)程度和危害程度也將呈繼續(xù)上升趨勢[1]。目前，風(fēng)暴潮災(zāi)害已位居中國海洋災(zāi)害首位[2]。上海由于其特殊的地理位置，常遭受臺風(fēng)引起的強(qiáng)風(fēng)暴潮侵襲。因此，建立一個高精度的風(fēng)暴潮數(shù)值預(yù)報模型，對于其海岸保護(hù)系統(tǒng)的設(shè)計、風(fēng)險評估和應(yīng)急疏散規(guī)劃都具有非常重要的意義。

目前，運用最廣泛的風(fēng)暴潮模式主要有英國的海模型(sea model)[3]、美國的颶風(fēng)暴潮模式(sea,lake and overland surges from hurricanes, SLOSH)[4-5]、有限體積海岸海洋模型(finite-volume coastal ocean model, FVCOM)[6-7]和有限元海洋模式(advanced circulation model,ADCIRC)[8]。中國從20世紀(jì)80年代開始，利用引進(jìn)的風(fēng)暴潮模型進(jìn)行大量的數(shù)值模擬研究和海洋預(yù)報業(yè)務(wù)。文獻(xiàn)[9]利用FVCOM對發(fā)生在廈門灣附近海域的“莫蘭蒂”臺風(fēng)引起的風(fēng)暴潮進(jìn)行模擬，文獻(xiàn)[10]利用ADCIRC對上海地區(qū)歷史風(fēng)暴潮進(jìn)行大量數(shù)值模擬研究，雖然兩者都對增水分布、風(fēng)生流場等風(fēng)暴潮特征進(jìn)行了較深入研究，但未考慮波浪的作用。波浪與風(fēng)暴潮共存于同一水體，相互作用又相互影響，在具有寬闊大陸架的區(qū)域，波浪可使水位增加5%～15%，在具有狹窄大陸架或陡坡的區(qū)域，水位可增加35%～40%[11]。文獻(xiàn)[12]利用ADCIRC和淺水波浪數(shù)值模型(simulating waves nearshore,SWAN)建立耦合模型,對珠江口風(fēng)暴潮、近岸浪開展數(shù)值模擬與預(yù)報，文獻(xiàn)[13]利用FVCOM和SWAN建立耦合模型,對汕尾沿岸海域影響較大的兩次臺風(fēng)風(fēng)暴潮過程進(jìn)行模擬，分析浪流耦合作用對汕尾風(fēng)暴潮增水的影響，雖然模擬精度都有所提高，但受制于引進(jìn)的風(fēng)暴潮模型，未能充分考慮波浪與風(fēng)暴潮的相互作用。近年來，國內(nèi)學(xué)者自主建立的一些風(fēng)暴潮模型也很少考慮波浪的影響作用[14-16]。另一方面，研究表明不同地區(qū)波浪對風(fēng)暴潮的影響有很大差異[17-19]，而上海沿海地區(qū)地處長江口和杭州灣交匯處，風(fēng)、浪、潮的相互作用復(fù)雜，但相關(guān)研究較少[20]。

本文基于二維淺水方程建立風(fēng)暴潮模型，利用高分辨率的三角形非結(jié)構(gòu)網(wǎng)格與SWAN建立耦合模式。該模式充分考慮波浪的輻射應(yīng)力與表面風(fēng)應(yīng)力作用，同時考慮風(fēng)暴潮過程對波浪的影響。通過對臺風(fēng)溫妮(Winnie)期間上海沿海地區(qū)復(fù)雜的風(fēng)暴潮增水過程進(jìn)行模擬分析，討論了波浪對該地區(qū)風(fēng)暴潮數(shù)值模擬的影響，對提高區(qū)域風(fēng)暴潮預(yù)報精度以及潮流-波浪耦合研究均具有一定的實際應(yīng)用價值。

1 臺風(fēng)Winnie簡況

1997年第11號臺風(fēng)Winnie于8月10日8時在關(guān)島東北偏東洋面上生成，即北緯15.4°、東經(jīng)153.8°，當(dāng)時中心氣壓995 hPa，近中心最大風(fēng)力為8級(20 m/s)，朝西北偏西方向移動，14時發(fā)展成強(qiáng)熱帶風(fēng)暴。11日8時到達(dá)北緯16.3°、東經(jīng)150.1°時發(fā)展為臺風(fēng)，當(dāng)時中心氣壓975 hPa，近中心最大風(fēng)力12級(風(fēng)速33 m/s)，并繼續(xù)朝西北偏西方向移動。12日20時到達(dá)北緯18.1°、東經(jīng)144.8°時達(dá)到鼎盛時期，當(dāng)時中心氣壓920 hPa，近中心最大風(fēng)速達(dá)60 m/s，并繼續(xù)向西北偏西方向移動，向浙江中南部靠近，18日晚21:30在浙江溫嶺登陸，登陸時中心氣壓955 hPa，近中心最大風(fēng)速為40 m/s，臺風(fēng)登陸后中心氣壓上升，近中心風(fēng)力減弱，于19日8時降為強(qiáng)熱帶風(fēng)暴，14時降為熱帶風(fēng)暴，后經(jīng)浙江北部、安徽、山東、遼寧省于21日2時減弱為低壓。

8月18日至8月20日上海地區(qū)遭受此臺風(fēng)影響嚴(yán)重。臺風(fēng)Winnie(199711號)影響上海期間，適逢農(nóng)歷七月半天文大潮(一年中天文潮最高時段)、風(fēng)潮遭遇，沿杭州灣、長江口和黃浦江干流各水文站均出現(xiàn)有記錄以來的最高潮位。杭州灣沿岸的各站比原歷史記錄最高潮位抬高42～64 cm，其中金山嘴站達(dá)6.57 m；沿長江口各站抬高13～36 cm，其中外高橋站達(dá)5.99 m；沿黃浦江各站抬高24～50 cm，其中吳淞站達(dá)5.99 m，黃浦公園站為5.72 m，吳涇站為4.82 m，米市渡站為4.27 m。

2 風(fēng)暴潮耦合模式的建立

2.1 風(fēng)暴潮耦合模式

風(fēng)暴潮過程中潮流、風(fēng)生流組成的復(fù)合流場的運動規(guī)律可以用守恒性的淺水方程來描述。本文應(yīng)用守恒性高、區(qū)域適應(yīng)性強(qiáng)、分辨率高的非結(jié)構(gòu)網(wǎng)格有限體積法離散淺水方程，運用HLLC格式計算界面通量，利用動邊界跟蹤的模式來處理干、濕交替的邊界問題，建立了二維風(fēng)暴潮水動力數(shù)值預(yù)報模型[20]。

采用SWAN模擬風(fēng)浪在向岸邊傳播過程中各項因素對風(fēng)浪形態(tài)的影響。SWAN 能夠系統(tǒng)地模擬波浪產(chǎn)生過程、耗散過程及傳播過程，并且源代碼公開，被廣泛應(yīng)用于近海岸、湖泊及河口等淺水地區(qū)的波浪研究。

圖1 風(fēng)暴潮和波浪耦合數(shù)值模式流程圖

本文通過運用高分辨率非結(jié)構(gòu)三角形網(wǎng)格將風(fēng)暴潮模型與SWAN進(jìn)行雙向耦合，能夠更好地反映風(fēng)暴潮與波浪之間的非線性相互作用。耦合計算過程中，每個循環(huán)時間間隔內(nèi)，風(fēng)暴潮模型首先運行，它將上一個時間間隔中SWAN模型內(nèi)插波浪參數(shù)所得到的表面風(fēng)應(yīng)力和輻射應(yīng)力加入到包括天文潮在內(nèi)的風(fēng)暴潮過程的水動力學(xué)計算中，得到水位和流速。在這個時間間隔的最后一個時間步長中，風(fēng)暴潮模型將該間隔內(nèi)計算得到的風(fēng)場、水位和流場數(shù)據(jù)作為計算條件代入到SWAN，然后SWAN利用這些數(shù)據(jù)計算有潮流影響的臺風(fēng)浪過程，并計算出新的表面風(fēng)應(yīng)力和輻射應(yīng)力，從而傳遞給風(fēng)暴潮模型用于下一個時間間隔潮流場的計算并繼續(xù)反饋給SWAN，實現(xiàn)雙向耦合的循環(huán)迭代。耦合模式的計算流程如圖1所示。

2.2 網(wǎng)格及模式參數(shù)

本文研究范圍為上海沿海地區(qū)，涵蓋了小尺度的河道(黃浦江)，中尺度的長江口和杭州灣，以及大尺度的東海部分海域。由于研究區(qū)域尺度跨越很大，為提高預(yù)報精度，設(shè)定了一個相對較小但合理的研究區(qū)域。其中，黃浦江和長江上游邊界分別取至米市渡和徐六涇，杭州灣邊界取至乍浦，外海的北邊界取至呂泗以北40 km處，南邊界至象山，東邊界取至吳淞口以北200 km處。模型的計算區(qū)域約為120.9°E～123.6°E，29.5°N ～33.0°N，該區(qū)域水深分布圖及長江口的局部放大圖分別如圖2a和圖2b所示。

(a) 計算區(qū)域

圖3 計算區(qū)域三角形非結(jié)構(gòu)網(wǎng)格及站點分布

為了更好地適應(yīng)上海沿海地區(qū)河網(wǎng)密布及地形復(fù)雜的特點，本文采用三角形非結(jié)構(gòu)網(wǎng)格進(jìn)行研究區(qū)域網(wǎng)格劃分。耦合模式的網(wǎng)格圖和站點分布見圖3，共包含3 162個網(wǎng)格結(jié)點，5 634個網(wǎng)格單元。網(wǎng)格分辨率在外?？拷_邊界處約為8 km，隨著向近岸地區(qū)的靠近，網(wǎng)格逐漸加密，到長江口的淺水區(qū)域只有幾百米，而黃浦江僅為幾十米。

模型的初始條件為水位和流速。它們對外界動力的響應(yīng)十分迅速，因此假設(shè)模型計算開始時海面是靜止的，即初始時刻所有單元的流速均為0，水位取平均海平面高程。外海開邊界利用全球天文潮位預(yù)報系統(tǒng)TPXO7.0中的8個主要分潮(M2、S2、K2、N2、K1、O1、P1和Q1)的調(diào)和常數(shù)來給定。對于黃浦江和長江口的上游邊界，分別選用米市渡站和徐六涇站的實測數(shù)據(jù)；對于杭州灣的開邊界，選用乍浦站的實測數(shù)據(jù)。

3 風(fēng)暴潮的模擬與分析

利用建立的風(fēng)暴潮和波浪的耦合模式，對臺風(fēng)Winnie(199711號)在上海近海海域引起的風(fēng)暴潮增水情況進(jìn)行了模擬。為了分析風(fēng)暴潮數(shù)值模擬中波浪的作用，分別采用不考慮波浪對風(fēng)暴潮影響的純風(fēng)暴潮模式(記為“非耦合模式”)、考慮波浪和風(fēng)暴潮相互作用的耦合模式(記為“耦合模式”)進(jìn)行計算，并將兩種模式的計算結(jié)果與實測值進(jìn)行對比分析。

3.1 潮位的模擬與驗證

為了分析臺風(fēng)Winnie影響下的潮水位和流速變化過程，計算了1997年8月17日2時至8月20日23時的風(fēng)暴潮過程。為了獲取準(zhǔn)確的初始條件，模型啟動24 h后才開始預(yù)報。

在風(fēng)暴潮對上海影響嚴(yán)重的69 h模擬中，計算出風(fēng)暴潮和波浪耦合作用下的潮位過程線，與觀測潮位過程的相位一致，水位增幅基本接近。橫沙和吳淞兩個測站的潮位對比結(jié)果分別見圖4a和圖4b。最大風(fēng)暴潮增水出現(xiàn)在8月19日0點，計算值與實測值吻合較好，說明耦合模式在風(fēng)暴潮增水計算方面具有較高精度。

(a) 橫沙

3.2 潮流場的變化特征分析

非耦合模式和耦合模式計算所得臺風(fēng)Winnie最高潮位時刻長江口潮流場分別見圖5a和圖5b，杭州灣潮流場分別見圖6a和圖6b。非耦合模式下長江口與杭州灣交匯處產(chǎn)生了3.16 m/s的最大流速，耦合模式下該區(qū)域最大流速為3.28 m/s，比未考慮波浪作用的情況增速0.12 m/s，且流速超過3 m/s的區(qū)域略微增大。同時，長江口入海處流速分布也發(fā)生一些改變，北槽中西側(cè)和佘山附近的流速都增大了0.2 m/s。非耦合模式下杭州灣開邊界處有高達(dá)3 m/s往東北方向的流速，且區(qū)域較大。而耦合模式下該處最大流速僅有2.2 m/s。附近區(qū)域流速在迎面而來的潮流影響下迅速轉(zhuǎn)向，且速度不到1 m/s。但上部區(qū)域在波浪的影響下流速由不到0.8 m/s增速至1.2 m/s左右。流速等值線的分布也發(fā)生了較大改變，流速大于2.0 m/s的區(qū)域雖然減少，但1.0～1.4 m/s流速區(qū)域增大很多。

(a) 非耦合模式

(a) 非耦合模式

從流場等值線的分布與流速極值發(fā)生的位置和大小來看，波浪作用在長江口區(qū)域?qū)α鲌龅挠绊戄^小，而在杭州灣其作用卻十分顯著。最高潮位時刻的風(fēng)場是由東往西的方向，表面風(fēng)應(yīng)力隨風(fēng)拖曳力系數(shù)的增大而增強(qiáng)，致使流向長江口和杭州灣的潮流流速增大。然而,長江口水深較淺，波浪在向長江口傳播時，水深變淺，再加上床面摩阻力的影響，在淺水中產(chǎn)生變形和破碎，波浪破碎以后，造成波高的減小，波浪作用也相應(yīng)減小，所以該區(qū)域潮流場大部分區(qū)域變化較小。而杭州灣的水深較長江口相對較大，岸線較平緩，波浪在向岸傳播過程中，引起波高在傳播過程中的空間變化較長江口小，波浪的空間變化也相對較小，因此傳播至杭州灣的波浪能量損失較小，在其影響下流入杭州灣的潮流流速增大。但最高潮位時刻該潮流正巧遇上杭州灣的退潮，這就造成杭州灣開邊界處退潮的最大流速大大降低。

3.3 最高潮位變化特征分析

非耦合模式和耦合模式計算所得最高潮位時刻長江口潮位分布分別見圖7a和圖7b，杭州灣潮位分布分別見圖8a和圖8b。最高潮位時刻長江口潮位升高明顯，北支支口涌入的潮量在其最狹窄段的下游形成6 m以上的高潮位，且5 m以上的高潮位范圍也從北上支口漫延至橫沙。杭州灣潮位變化不大，僅4.4～5.2 m的潮位范圍有所縮小。這是由于耦合模式潮流場的影響，使得潮水在北支支口大量匯集，而杭州灣的潮水則往長江口退去。在受波浪影響的表面風(fēng)應(yīng)力作用下，從外向杭州灣推進(jìn)的潮流增強(qiáng)，使得杭州灣高潮位等值線向開邊界處推移。

(a) 非耦合模式

(a) 非耦合模式

選取吳淞和橫沙兩個潮位站作為驗證點，將兩種模式的計算結(jié)果與實測值進(jìn)行對比，結(jié)果分別如圖9a和圖9b所示。與非耦合模式相比，耦合模式兩個站點的計算潮位都與實測潮位更吻合，特別是最高潮位，兩個站點的最高潮位計算值都更接近實測值。

(a) 橫沙

吳淞、橫沙、北槽中和長興4個潮位站兩種模式最高潮位誤差對比結(jié)果見表1。由表1可知：非耦合模式4個站點最高潮位的相對誤差為8.8%～10.1%，耦合模式預(yù)報精度明顯有提高，相對誤差為6.2%～7.9%，平均相對誤差為7.2%，較非耦合模式提高2.3%。這是因為波浪造成的水面小尺度的波動改變了水面的表觀粗糙度，同時波浪破碎會消耗大量的波能，輻射應(yīng)力也隨之而變，所以耦合模式能更好表現(xiàn)臺風(fēng)對大尺度潮流場的驅(qū)動力，從而改善了潮位增水。

表1 耦合模式和非耦合模式最高潮位誤差對比

4 結(jié)論

(1)構(gòu)建的風(fēng)暴潮與波浪耦合模式，充分考慮了風(fēng)暴潮與波浪的相互作用，并能有效模擬出臺風(fēng)Winnie引起的風(fēng)暴潮潮位變化過程，經(jīng)驗證4個站點最高潮位平均相對誤差僅7.2%，相對于未考慮波浪作用的純風(fēng)暴潮模型模擬精度提高了2.3%。

(2)波浪對上海沿海地區(qū)風(fēng)暴潮潮流場作用受地形影響較大，最高潮位時刻，長江口潮流場受波浪影響較小，而杭州灣則影響較大。

(3)臺風(fēng)Winnie影響期間，波浪對上海地區(qū)不同水域風(fēng)暴潮最高潮位影響也有區(qū)別，由于波浪和天文潮的綜合作用，潮水向長江口匯集，造成長江口潮位明顯增加。