高 峰, 朱川林, 李 昊, 唐友剛, 黃杰斯, 張慈珩, 劉 針, 姜云鵬
(1. 交通運輸部天津水運工程科學(xué)研究所,港口水工建筑技術(shù)國家工程實驗室,工程泥沙交通行業(yè)重點實驗室,天津 300456;2. 天津大學(xué)建筑工程學(xué)院,天津 300072; 3. 海南金海灣投資開發(fā)有限公司,海南 海口 570000)
中國是世界上海洋災(zāi)害最嚴重、最頻發(fā)的少數(shù)國家之一,而四面環(huán)海的海南省又是最易受到海洋災(zāi)害影響的地區(qū),對于地處海南省經(jīng)濟發(fā)展前沿的最大城市——??谑卸?,其瀕臨??跒乘蛩媾R的海洋災(zāi)害形勢十分嚴峻,超強臺風、風暴潮、極端大浪、海嘯、海平面持續(xù)上升、岸灘侵蝕加劇以及海岸生態(tài)環(huán)境變化加速等,均對沿海經(jīng)濟造成威脅相當嚴重,已成為沿海開發(fā)規(guī)劃實施和經(jīng)濟建設(shè)的制約因素之一。項目用海災(zāi)害風險分析主要是對項目災(zāi)害性風險的種類、項目用海對海洋功能的風險、項目用海對周邊利益相關(guān)者的風險進行論證和分析[1],進而提出避免災(zāi)害性風險發(fā)生的方案,或是災(zāi)害性風險發(fā)生后的應(yīng)急措施。因此,開展沿海工程海洋災(zāi)害風險評價分析,是有效提高災(zāi)害防御能力、控制風險和減少災(zāi)害損失的關(guān)鍵技術(shù)方法之一,對該地區(qū)近海涉水工程的建設(shè)具有重要意義。
目前,一般風險分析的程序是:風險識別、風險估計、風險評價、風險處理和風險決策[2]。來自海洋水動力的風險主要包括兩個方面:(1)一方面是用海項目自身引起的突發(fā)或緩發(fā)事件導(dǎo)致對海域動力環(huán)境造成的危害,發(fā)生于施工期居多,如波浪增加、流速加快、懸浮物擴散、船舶溢油等;(2)另一方面是由于自然災(zāi)害對海域使用項目造成的危害,發(fā)生于營運期居多,如遭受熱帶氣旋、風暴潮、海嘯等海洋災(zāi)害對工程造成的損壞。海南省四面環(huán)海、瀕臨南海,地理位置特殊,??跒澄挥诤D蠉u北部的海岸線上,分布于瓊北中部岬間海灣內(nèi),隔瓊州海峽與大陸的雷州半島相望,背依??谑屑俺芜~縣的馬村(老城)工業(yè)開發(fā)區(qū)。其所在海域的海洋水動力災(zāi)害主要包括:風暴潮、海浪、海嘯等。按照一般風險分析的程序,就近岸海洋動力風險分析研究來說,重點在于上述災(zāi)害的風險識別與估計,本文即通過實測資料分析結(jié)合數(shù)值模擬計算,確定上述各主要風險因子發(fā)生的強度大小和破壞情況。
登陸我國的臺風和強溫帶天氣過程往往造成風暴潮災(zāi)害,其成災(zāi)頻率高、致災(zāi)強度大,造成的人員和經(jīng)濟損失慘重。風暴潮災(zāi)害主要是由異常的風暴增水使潮位大幅升高而形成,致災(zāi)因子不僅包括風暴潮,還包括天文大潮、近岸浪及其三者之間的耦合作用。影響海南島北部的熱帶氣旋較多[3],是受風暴潮影響較為嚴重的區(qū)域之一,海水堆積引發(fā)的潮災(zāi)是??跒乘蛑饕Q鬄?zāi)害之一。根據(jù)相關(guān)資料統(tǒng)計,海口年均發(fā)生風暴潮約3次,而嚴重或特大風暴潮約每2年1次[4]。海口灣每年6—12月均有可能發(fā)生風暴潮,且集中發(fā)生在7—10月,尤以10月份最多。在對歷史風暴潮資料進行搜集整理的基礎(chǔ)上,本文研究基于現(xiàn)有的實測臺風資料通過將Monte Carlo法[5](又稱為隨機統(tǒng)計模型)與水動力模型相結(jié)合的方法實現(xiàn)對可能出現(xiàn)的臺風風暴潮進行數(shù)值模擬研究,從而推求出重現(xiàn)期為50、100以及500 a的極端高水位預(yù)測值。
1.1.1 計算方法 采用假定水體上下層均勻一致的垂向平均二維風暴潮數(shù)值模型作為近海臺風風暴潮的水動力模型,且二維模型突出了風暴潮預(yù)報與計算的最重要變量—風暴潮位,這也是預(yù)測臺風風暴潮引起的極端高水位的核心變量。模型采用大小模型嵌套方式,其中大模型的北側(cè)和西側(cè)以大陸岸線為陸地邊界,南面以北緯15.0°為外海邊界,東面以東經(jīng)117.9°為外海邊界,模型范圍覆蓋了北部灣、瓊州海峽、瓊州海峽以東及南海部分海域,這種規(guī)模的大尺度計算網(wǎng)格可以很好地對一次臺風影響過程進行數(shù)值模擬,其網(wǎng)格尺寸約為1 km,共有420×808個網(wǎng)格點,339 360個單元。同時,小模型南北兩側(cè)以大陸岸線為陸地邊界,南北向跨度最大處約180.0 km,東西向從108.7°E~110.9°E,跨度約為229.2 km。計算范圍及網(wǎng)格劃分如圖1(a)~(b)所示。
((a)大模型; (b)小模型; (c)1409號臺風模擬; (d)驗證結(jié)果。(a)Large scale model; (b) Local nested model; (c)Typhoon 1409 simulation; (d)Verification results. )
圖1 模型網(wǎng)格嵌套分布與1409號臺風模擬驗證情況
Fig.1 Distribution of model nested grid and verification result of No.1409 typhoon’s simulation
氣壓分布模型采用Holland模型,環(huán)流風速模型采用梯度風速公式,移行風速模型采用宮崎正衛(wèi)模型。由氣壓分布模型和風速模型組合成臺風風場的模擬方案,簡稱Holl模型[6]。典型的實測臺風依據(jù)美國FEMA使用的臺風歷史記錄數(shù)據(jù)處理方法—模擬圓法,提取在一定時期內(nèi)途徑以要模擬的城市為模擬點,以設(shè)定的距離為半徑的圓的所有臺風,以這些臺風的記錄作為研究對象[7-8]。在1949—2015年期間曾有151個臺風經(jīng)過以??跒碁橹行?,200 km為半徑的圓形研究區(qū)域,篩選登陸點中心最低氣壓在0~980 hPa范圍內(nèi)共計60個。在所有的60個實測臺風中,2014年7月18日經(jīng)過的“威馬遜”臺風強度最大,故利用鐵山港海區(qū)英羅灣附近氣象觀測站實測資料對該場臺風進行了風場驗證,“威馬遜”登陸后整個研究區(qū)域于2014年7月18日18:00的氣壓場與風場模擬情況可分別見下圖1(c)所示,模擬與實測的最大風速擬合良好圖1(d)。進而采用Monte Carlo模型生成100年內(nèi)可能出現(xiàn)的模擬臺風數(shù)據(jù)。由于臺風歷史資料采用的是1949—2015年間的60個典型臺風數(shù)據(jù),則500年內(nèi)應(yīng)當產(chǎn)生的臺風數(shù)為448個。對這些個可能出現(xiàn)臺風的參數(shù)進行模擬,可得所需的模擬臺風數(shù)據(jù)。采用隨機模型(Monte Carlo)隨機模型和水動力模型(Delft3D)相結(jié)合的方法得到水位過程線,并得出極值水位。首先采用Monte Carlo隨機模型,將篩選出的67年里60條臺風的氣壓和路徑分別作為獨立的元素隨機組合得到448個重組臺風。將重組臺風要素帶入已驗證的風暴潮水動力模型進行計算,其中臺風過程與天文潮位隨機組合,最后得到448個水位過程線,統(tǒng)計每條水位過程線中的最大水位值,用于后面的極端水位分析。
1.1.2 計算結(jié)果 通過以上過程,統(tǒng)計每條水位過程線中的最大水位值,部分較大水位過程線如圖2(a)所示。將統(tǒng)計出的448個最大水位值進行降序排列,計算得秀英站處最大水位序列的最高值達3.75 m。極端水位的回歸周期TR可按水位的排序參照式:TR=500/M估計,其中,M=8是448個最高水位的從高到低的排序,TR為重現(xiàn)期。例如,當M=1(對應(yīng)的是448個極端水位的最大值)時,這個水位的最大值對的重現(xiàn)期是TR1=500/1=500 a。相應(yīng)的極端水位所對應(yīng)的重現(xiàn)期可點繪成圖,且可由該圖確定出任一重現(xiàn)期所對應(yīng)的極端水位值。第五高水位3.02 m所對應(yīng)的是100年一遇的極端高水位值,重現(xiàn)期為500、100和50年所對應(yīng)的秀英站極端高水位預(yù)測值分別為3.75、3.02和2.68 m。
另外,通過采用重現(xiàn)期100年一遇氣壓P0=883 hPa,最大風速半徑R為33 km,臺風移速為35 km/h,采用假想路徑,對可能最大風暴潮增水(PMSS-ProbableMaximum Storm Surge)[9]進行了計算,計算結(jié)果表明最大增水達到4.42 m,其對應(yīng)增水過程見圖2(b)。
((a)水位過程線; (b)最大可能增水過程。(a)Stage hydrograph; (b)Maximum possible water increase process.)
圖2 各重現(xiàn)期秀英站臺風暴潮水位過程線和最大可能風暴增水過程圖
Fig.2 Storm surge water level process of different return period for Xiuying and the greatest possible storm surge elevation
海口灣屬臺風作用活躍區(qū)域,受臺風和大陸季風作用明顯,但該水域及周邊缺乏長期觀測站,現(xiàn)僅有玉苞角、白沙門兩個臨時站。影響該海域的臺風主要來自南海自東或東南經(jīng)廣東、海南登陸或從海南島南部海域向偏西方向移動的臺風。??跒乘蛉暧坷溯^少、以風浪為主,風浪頻率為76%~85%,涌浪頻率占14%~23%,常浪向為ENE,頻率為30.1%,次常浪向為NE,頻率為22.9%。波浪出現(xiàn)最少的方向為S~WSW。受季風風向影響,工程海域常浪向隨季節(jié)風而變化,東北季風期(11 月—次年3 月),常浪向為ENE,次常浪向為NE,季風轉(zhuǎn)換期(4—5 月,9—10 月)常浪向、次常浪向與東北季風期一樣,但常浪向出現(xiàn)頻率略低于冬季,次常浪向出現(xiàn)頻率略高于冬季;西南季風期(6—8 月)常浪向為NE。目前已知最大波浪由1409 號威馬遜臺風超強臺風產(chǎn)生,該超強臺風在??跒车臇|北側(cè)通過,極有利于西北向波浪的成長和傳入??跒?,在??跒吵霈F(xiàn)最大波浪,形成有效波高6.8 m,周期約7.8 s,該波高為偏西波向的異常特大值,在N-E-S方位的最大有效波高3.9 m,周期約6.3 s。
1.2.1 計算方法 該海域主要成風原因有兩種,一為季風、二為氣旋。因此,研究主要通過風場模擬計算波浪場,風場資料需綜合考慮季風和臺風兩種影響。具體研究思路如下:
(1)采用中國氣象數(shù)據(jù)網(wǎng)(www.data.cma.cn)公布的臺風路徑,結(jié)合歐洲中長期在南海海域的后報風場資料,選取每年影響工程海區(qū)的臺風及季風過程,對風場進行分析。
(2)根據(jù)多年風場數(shù)據(jù),采用SWAN模式建立大區(qū)域波浪計算模型,得到每年工程區(qū)附近外海深水區(qū)的年波浪極值,進而采用Pearson-Ⅲ型曲線適線法推算不同重現(xiàn)期波要素。
(3)根據(jù)外海不同重現(xiàn)期波浪條件,在自然地形條件下即工程建設(shè)前,利用SWAN模式建立小區(qū)域波浪計算模型,推算工程區(qū)設(shè)計波浪條件。
1.2.2 主要結(jié)果與分析 根據(jù)工程海區(qū)最近30年(1986—2015年)臺風和季風資料,通過臺風場模型和季風場數(shù)據(jù)分析計算工程海區(qū)大風過程,進而采用大范圍風浪模型計算工程海區(qū)風浪場,從中選取主要計算浪向ENE、NE、NNE、N、NNW與NW向的波浪年極值,分析工程外海深水區(qū)不同重現(xiàn)期波浪。提取??跒潮眰?cè)外海計算點(水深約-70 m)處的主要計算方向波浪重現(xiàn)期結(jié)果見下表1,其中前文所提到的60個典型臺風中,分別選取1522號“啟德”臺風與1117號“飛燕”臺風過境瓊州海峽時波高分布圖(見圖3),其中重現(xiàn)期100年有效波高為7.14 m,重現(xiàn)期50年有效波高為6.80 m。
表1 海口灣不同重現(xiàn)期有效波高結(jié)果Table 1 Results of different return period of significant wave height in Haikou-Bay /m
另外,依據(jù)上述波高結(jié)果,通過HUDSON公式對海口灣近岸工程可能采用斜坡式護岸或者防波堤的護面塊體重量進行了核算,如目前海口灣內(nèi)各人工島如普遍采用30 t護面塊體時,其對應(yīng)破壞波高可達8.78 m,當于重現(xiàn)期為1 400年波浪(年概率為0.071%),若以設(shè)計年限100年考慮,這一等級的塊體穩(wěn)定性保障率可達93%。同時,根據(jù)爬高和越浪的計算分析,極端高水位條件下有效波高超過7.9 m就會出現(xiàn)危害堤頂鋪面結(jié)構(gòu)安全的越浪量,該波浪條件對應(yīng)的重現(xiàn)期約為400年,既年出現(xiàn)頻率為0.25%。
圖3 典型臺風工程海域作用過程Fig.3 Functional process of typical typhoon in engineering sea area
隨著波浪由外海向近岸傳播,充分考慮了灣內(nèi)自然水深地形條件下波浪分布。采用SWAN模型對工程區(qū)的波浪分布情況進行計算分析。計算方向為ENE、NE、NNE、N、NNW與NW向。包含100年一遇高水位、設(shè)計高水位與設(shè)計低水位。其中,ENE最大有效波高為6.7 m(見圖4(a)),NE向最大有效波高5.8 m(見圖4(b)),NNE向為5.5 m,N向為3.8 m,NNW向為4.4 m,NW向為4.9 m。對??跒辰队绊戄^大的2個方向為ENE與NE向,重現(xiàn)期100年一遇高水位時H13%波高分布見圖4所示。
圖4 100年一遇高水位重現(xiàn)期100年有效波高分布結(jié)果Fig.4 Distribution of significant wave height for 100 years return period by a high water of 100 years return period
依據(jù)美國地質(zhì)調(diào)查局(USGS)提供的一份潛在海嘯地震斷裂帶風險圖[10],我國南海海域附近有三大俯沖帶被認為具有地震激發(fā)海嘯的潛在危險,分別為馬尼拉俯沖帶、琉球群島俯沖帶和北蘇拉威西島俯沖帶。其中,馬尼拉俯沖帶(又稱為馬尼拉海溝)對南海周圍的國家和地區(qū)會帶來很大的海嘯威脅,其位于亞歐大陸板塊和菲律賓海板塊交界處,南起菲律賓巴拉望島北端,沿菲律賓呂宋西部邊緣向北發(fā)展,北至臺灣島,總長度約1 000 km[11]。USGS根據(jù)馬尼拉海溝的方位特征和斷層幾何學(xué)將其劃分為6個斷裂帶[12],沿馬尼拉海溝分布如圖5(b)所示。
上述區(qū)域地震激發(fā)的海嘯波會使海南等我國南部海域暴露在海嘯波威脅之下,雖然歷史統(tǒng)計資料顯示幾次大的跨洋海嘯對我國東南沿海地區(qū)影響有限,但南海馬尼拉海溝地震帶形成的海嘯仍會對臺灣、海南、廣東和福建等沿海地區(qū)曾造成不同程度的影響,南海曾發(fā)生過約15次局地海嘯,造成了一定影響[13]。由于近幾年全球地質(zhì)運動較為活躍,特大地震海嘯發(fā)生頻率較高,馬尼拉海溝被認為是最有可能引發(fā)海嘯的潛在地震源?;谏鲜龇治?,仍不可忽視海嘯波的潛在威脅。
2.3.1 計算方法 本文采用淺水波數(shù)值模型GeoClaw,基于有限體積法對非線性淺水方程進行數(shù)值離散,考慮了海嘯波在近岸傳播的非線性作用,同時嵌入了自適應(yīng)網(wǎng)格技術(shù)可追蹤海嘯波高判斷網(wǎng)格加密與否。通過對馬尼拉海溝斷層特征的分析,Lau等[14]根據(jù)歷史上曾經(jīng)在附近區(qū)域發(fā)生過的地震和經(jīng)驗關(guān)系公式,確定了7.0、8.0級地震時E1~E3板塊的數(shù)值模擬參數(shù)(見表2),并基于Nguyen等[15]在USGS震源參數(shù)基礎(chǔ)上給出的震級為9.3級的震源參數(shù),此時對應(yīng)各板塊震源參數(shù)見表3,從而依此計算相應(yīng)的傳播對海南島周邊海域的影響。
( (a)南海附近俯沖帶;(b)馬尼拉海溝不同板塊分布。(a)Subduction zone near the south China sea。(b)Distribution of different plates in the Manila trench.)
圖5 中國南海附近的俯沖帶與馬尼拉海溝不同板塊產(chǎn)生的海底變形及位置
Fig.5 Sea bottom deformation and location of subduction zone and different tectonic plates from the Manila trench nearby China South Sea
利用Okada彈性半空間錯位理論模型,根據(jù)前文所述六塊斷裂帶參數(shù)生成相應(yīng)的海底變形可分析出海南島受E1、E2和E3板塊的影響較大。據(jù)此將對上述各板塊(震源深度15 km)分別發(fā)生7.0級和8.0級地震形成的海嘯進行數(shù)值模擬,進而估算對??跒乘虻挠绊懬闆r,還考慮了整個斷裂帶(E1~E6)發(fā)生破裂形成的特大地震海嘯時的極端情況。計算范圍從99°E~130°E,1°S~33°N,包含了我國東南沿海至南海海域以及菲律賓、越南等國家。
表2 震級分別在7.8和8.0時馬尼拉海溝E1~E3斷裂帶震源模擬參數(shù)Table 2 Source modeling parameters of E1 ~ E3 fault zone in Mw=7.8 and 8.0 respectively in the Manila trench
表3 Mw=9.3時,馬尼拉海溝E1~E6斷裂帶基本參數(shù)Table 3 Source modeling parameters of E1 ~ E6 fault zone in Mw=9.3 in the Manila trench
表4 各計算工況下工程區(qū)海嘯波特征值 Table 4 Characteristic value of tsunami under the different calculation conditions
2.3.2 結(jié)果與分析 數(shù)值模擬結(jié)果如圖6所示,分別給出了馬尼拉海溝北部3個板塊(E1、E2、E3)分別發(fā)生里氏7.0、8.0級海底地震形成的海嘯波。各工況下工程區(qū)的海嘯波特征值如表4所示。從走時圖中可以看到,海嘯波大概經(jīng)歷4~5 h即可到達??跒掣浇S?。對于單一板塊的一般震級地震,海嘯波到達海口附近波幅都小于4 cm,波長在50 km左右,周期1.6~2 h,海嘯持續(xù)時間(第一個大波的作用時間)約為0.7~1.1 h,海嘯波行進速度為8.1 m/s。對于整個斷裂帶發(fā)生9.3級大地震引發(fā)的極端海嘯,海嘯波到達工程區(qū)附近海域時最大波幅為2.25 m,波長115 km,周期3.4 h,海嘯持續(xù)時間1.3 h,海嘯波行進速度為9.4 m/s,該情況將對海口灣會產(chǎn)生影響,特別是岬灣形態(tài)的海口灣近岸,更易形成能量集中,應(yīng)予以關(guān)注。
圖6 各板塊分別在7.0和8.0級地震生成的海嘯波走勢圖Fig.6 The tsunami wave distribution of each plate in Mw=7.0 and 8.0 generated by earth quake
海洋水動力風險是用海風險的重要組成部分,關(guān)鍵在于風險識別和風險估測。本文根據(jù)對海口灣水域可能出現(xiàn)的風暴潮、極端波浪、海嘯等問題的研究,對該海域水動力風險因素進行了研究分析,得到風暴潮不同重現(xiàn)期增水水位,并在確定外海重現(xiàn)期100 a有效波高為7.14 m,重現(xiàn)期50 a有效波高為6.8 m的基礎(chǔ)上,計算了??跒辰恫煌ㄏ虻淖畲蟛ǜ撸瑢硟?nèi)近岸斜坡式防護工程采用護面塊體的重量校核。針對潛在海嘯威脅最大的馬尼拉海溝,進行了多種工況數(shù)值模擬,除9.3級大地震引發(fā)極端海嘯導(dǎo)致海嘯波達2.25 m之外,7、8級海嘯對于海口灣所受影響較小。這些計算結(jié)果反映海口灣沿岸在不同水位下的漫灘情況,對??跒逞匕吨匾墓こ添椖恳约胺罏?zāi)規(guī)劃的制定,具有參考意義。