馬尼拉地震帶海嘯對廣東沿海影響的數(shù)值研究
——潮汐與海嘯的非線性作用的影響

曹永港 劉長建 劉同木馮硯青 劉愉強(qiáng) 廖世智

1）國家海洋局南海調(diào)查技術(shù)中心，廣州市海珠區(qū)新港西路155號 510300

2）天津大學(xué)水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津市南開區(qū)衛(wèi)津路92號 300072

3）國家海洋局南海維權(quán)技術(shù)與應(yīng)用重點(diǎn)實(shí)驗(yàn)室，廣州市新港中路353號 510310

4）水沙科學(xué)與水災(zāi)害防治湖南省重點(diǎn)實(shí)驗(yàn)室，長沙市天心區(qū)萬家麗南路2段960號 410114

0 引言

海嘯是由水下地震、火山爆發(fā)、水下塌陷及滑坡等大地活動或氣象變化產(chǎn)生的破壞性海浪，其中，由水下地震產(chǎn)生的地震海嘯較為普遍，表現(xiàn)為海底地形急劇升降變動引起的海水強(qiáng)烈擾動。根據(jù)海嘯傳播至岸邊時間的長短，它可分為近海地震海嘯和遠(yuǎn)洋地震海嘯。海嘯的發(fā)生頻率雖然較低，可一旦發(fā)生，破壞力是巨大的。2004年12月26日印度尼西亞蘇門答臘島西北附近海域發(fā)生了MW9.1強(qiáng)烈地震，地震在印度洋沿岸引發(fā)罕見的地震海嘯并傳播到全球各大洋（Yeh et al，2007）。2011年3月11日日本東北部海域發(fā)生了MW9.0強(qiáng)烈地震，這是日本近海發(fā)生的最強(qiáng)烈的一次大地震，并引發(fā)了太平洋范圍的海嘯（王培濤等，2012）。據(jù)統(tǒng)計，每2年全球發(fā)生1次局地破壞性海嘯，每10年發(fā)生1次越洋大海嘯（于福江等，2011）。Jing等（2013）研究發(fā)現(xiàn)，地震引發(fā)海嘯的概率很高，而在過去約100年里，震級≥8.8的地震都觸發(fā)了海嘯。

我國處于太平洋西部，海區(qū)遼闊，海岸線長達(dá)18000km，近海海域大多位于環(huán)太平洋地震帶上，歷史上地震海嘯時有發(fā)生，是世界上最早記錄到地震海嘯的國家。據(jù)不完全統(tǒng)計，公元47～2004年，中國沿海共發(fā)生29次地震海嘯，其中，8～9次為破壞性海嘯。雖然自1949年新中國成立以來我國并未發(fā)生破壞性海嘯，但近期南海地震頻發(fā)，一旦誘發(fā)海嘯將給我國南部、東南部發(fā)達(dá)城市造成巨大損失（包澄瀾等，2005；楊馬陵等，2005；Liu et al，2007）。我國歷史上曾經(jīng)發(fā)生過多次海嘯，其中，臺灣地區(qū)周圍是海嘯的高發(fā)區(qū)域，其次是大陸架區(qū)域，低發(fā)區(qū)是渤海區(qū)域。在各省份中，按照有歷史記錄的海嘯分析，浙江的發(fā)生次數(shù)最多，高達(dá)45次，其次是江蘇、山東、上海、福建、臺灣、廣東。按照已確定的海嘯記錄分析，臺灣發(fā)生的次數(shù)最多，有8次，其次是山東、廣東和浙江。不論哪個數(shù)據(jù)都表明，我國是一個海嘯發(fā)生的危險區(qū)，而南海區(qū)域是中國受海嘯威脅最大的地方，這是因?yàn)閰嗡螎u西側(cè)和馬尼拉海溝一帶是地震活躍帶（潘文亮等，2009）。Li等（2016）認(rèn)為，常用的均勻滑動模型明顯低估了海嘯災(zāi)害。因此，研究建立廣東沿海海嘯預(yù)警機(jī)制，在海嘯發(fā)生后迅速預(yù)報海嘯到達(dá)時間、海嘯最大波高，可為海嘯防災(zāi)減災(zāi)提供有力支持，以期最大化地減少海嘯造成的損失。

現(xiàn)有的海嘯計算模式中，一般不含有潮汐模塊。這是由于通常不考慮潮汐與海嘯的非線性作用，將計算出的海嘯波高與潮位值進(jìn)行線性疊加來得到關(guān)心站點(diǎn)的水位。但是，對于近岸淺水域，潮汐與海嘯的非線性作用會明顯增大，故影響海嘯波的波幅及到達(dá)時間。本文利用COMCOT模型和全球潮汐模式TPXO7.2，建立廣東沿海天文潮與海嘯的耦合數(shù)學(xué)模型，針對馬尼拉地震帶海嘯，通過設(shè)計海嘯震源計算，分析廣東省沿海海嘯風(fēng)險。

1 數(shù)值模型

COMCOT模型（Cornell Multi-grid Coupled Tsunamimodel）是由Cornell大學(xué)土木與環(huán)境工程系Philip Liu研究組開發(fā)基于淺水長波方程的海嘯數(shù)值計算模式，該模式已被多次用來模擬歷史海嘯事件，如1992年印度尼西亞Flores Islands海嘯（Liu et al，1995）、2003年Algeria海嘯（Wand et al，2005）、2004年印度洋大海嘯（Wand et al，2006）以及2011年日本海嘯（應(yīng)超等，2015；Jing et al，2012），均有效模擬了海嘯生成、傳播、爬高和淹沒的整個過程。模式可根據(jù)研究范圍，靈活選擇球面與直角坐標(biāo)系下的線性與非線性淺水方程，所有方程均采用顯式蛙跳有限差分法進(jìn)行離散。

1.1 基本方程

球坐標(biāo)下，考慮科氏力忽略海底摩擦，線性淺水方程為

式中，η為相對于平均海平面的自由表面位移；R為地球半徑；H＝η＋h為總的水深，h為靜水深；P為沿緯度單位寬度的通量；Q為沿經(jīng)度單位寬度的通量；φ、ψ分別為經(jīng)度和緯度；f為科氏力系數(shù)；g為重力加速度。海嘯傳播至近岸，采用笛卡爾坐標(biāo)下非線性淺水波方程，并考慮底摩擦效應(yīng)，其連續(xù)方程和動量方程為

式中，P、Q分別為x和y方向上的體積通量；n為曼寧粗糙系數(shù)。

1.2 數(shù)值格式

對式（1）～（5）采用了蛙跳有限差分格式進(jìn)行離散，非線性對流項(xiàng)采用迎風(fēng)向上格式離散。圖1為交錯格式示意圖。由圖1可見，波高和體積通量在空間和時間上都是交錯的，波高及水深位于網(wǎng)格中心，體積通量位于網(wǎng)格線上。數(shù)值格式的截斷誤差為o（Δx2，Δy2，Δt2），在時空上均有二階精度。COMCOT是比較成熟的海嘯數(shù)值模型，具體算法詳見文獻(xiàn)（Wang，2009）。

圖1 交錯格式示意圖

1.3 邊界條件

海嘯初始條件的確定是通過輸入地震斷層參數(shù)，由彈性斷層模型計算得到。模式中可選擇Mansinha和Smylie的彈性半空間錯移模型和Okada的理論模型，這也是現(xiàn)階段常用的2種彈性斷層模型。彈性斷層模型基于彈性錯移理論建立，主要利用斷層錯動資料計算特定點(diǎn)處的位移量。例如，斷層在j方向上錯動Δui，對整個矩形斷層范圍積分，即得到

式中，δjk為克羅內(nèi)克函數(shù)變量；vk為∑面向外的垂直向量；μ、λ為拉梅常數(shù)；uji為受到破裂面在j方向的單位應(yīng)力作用而在地表i方向產(chǎn)生的位移。

COMCOT模型中的斷層模型共輸入9個參數(shù)，分別為震中經(jīng)度、緯度、震源深度、斷裂長度、斷裂寬度、滑動量、走向角、傾角、滑移角。斷層模型示意圖見圖2。

邊界條件的設(shè)定中，水邊界設(shè)為開邊界。選擇線性淺水波方程時，海陸邊界設(shè)為垂直反射邊界；選擇非線性淺水波方程時，采用移動邊界方案，即考慮岸灘的干濕變化。本文通過源碼修改在模型中加入潮汐強(qiáng)迫邊界條件（李林燕等，2012；趙鑫等，2016），外海潮位邊界由全球潮汐模型（TPXO7.2）求得，該模型通過10個分潮推算天文潮位，包含8個主要分潮M2、S2、K1、O1、N2、P1、K2、Q1以及2個長周期分潮Mf和Mn，基本能夠構(gòu)造出外海深水處真實(shí)的天文潮過程

圖2 斷層示意圖

式中，ζ0為邊界處的潮位；ζp為邊界處靜壓水位；i為1～10，分別對應(yīng)上述分潮；Ai、αi分別為分潮在3條邊界處的振幅和遲角；ωi為分潮的角頻率。

2 海嘯模型驗(yàn)證

2.1 計算區(qū)域與網(wǎng)格布置

利用日本“311”海嘯對COMCOT模型進(jìn)行驗(yàn)證。模型采用2層嵌套網(wǎng)格，第1層網(wǎng)格的范圍是5°S～62°N、105.0°～179.5°E，計算范圍包括東海、臺灣海峽、南中國海、日本及部分西太平洋海域，網(wǎng)格大小為4′，控制方程選擇球坐標(biāo)下的線性淺水長波方程，忽略摩擦作用；第2層網(wǎng)格的范圍是14.5°～28.0°N、110°～126°E，網(wǎng)格大小為0.8′（圖3）。

圖3 網(wǎng)格范圍

2.2 震源參數(shù)設(shè)置

日本“311”地震海嘯發(fā)生后，Tang等（2012）利用D21418浮標(biāo)站在震后0.42～0.58h和D21401在0.80～1.23h之間的監(jiān)測數(shù)據(jù)，反演出初始海面位移有10m的抬升和3m的下沉，震源可由6塊沿著日本海溝的100km×50km的單位震源組合而成，在海嘯到岸前5h計算出了美國海岸線30多個城市的海嘯波爬高和淹沒情況，模擬值與后來的觀測值較為符合。Wei等（2013）也利用同樣的斷層設(shè)置，模擬計算了海嘯在日本近岸的爬高和淹沒，結(jié)果表明，模擬計算的茨城縣和青森縣的淹沒范圍準(zhǔn)確率達(dá)到85.5%。以上述National Oceanic and Atmospheric Administration的斷層設(shè)置作為海嘯模型的輸入條件，模擬此次海嘯波的產(chǎn)生和傳播，并輸出在日本近岸的DART浮標(biāo)點(diǎn)以及我國近海驗(yàn)潮站的波高序列。地震斷層參數(shù)設(shè)計見表1。

表1 地震斷層參數(shù)

圖4 日本東部海嘯浮標(biāo)分布

2.3 浮標(biāo)站驗(yàn)證

本研究收集海嘯發(fā)生后日本東部編號為21401＃（42.617°N，152.583°E）、21413＃（30.515°N，152.117°E）、21418＃（38.688°N，148.769°E）、21419＃（44.455°N，155.736°E）的海嘯浮標(biāo)站水位資料進(jìn)行驗(yàn)證，浮標(biāo)分布如圖4所示，模型驗(yàn)證如圖5所示。在4個浮標(biāo)點(diǎn)處，首波到岸時間的模擬值與觀測值之差均小于2min，21401＃和21419＃浮標(biāo)站，波高的模擬值與觀測值之差小于2cm，浮標(biāo)點(diǎn)處的模擬值與實(shí)測值較為符合；在21418＃h和21413＃浮標(biāo)站，波高的計算值與觀測值的誤差分別為36、21cm，誤差相對較大，這與海底地形誤差、驗(yàn)潮站局部海浪初值累積誤差等有關(guān)。

圖5 浮標(biāo)點(diǎn)處模擬值與觀測值對比

3 計算結(jié)果分析

在我國南海海域附近，3大俯沖帶被認(rèn)為具有地震激發(fā)海嘯的潛在危險，分別為馬尼拉俯沖帶、琉球群島俯沖帶和北蘇拉威西島俯沖帶。馬尼拉海溝位于亞歐大陸板塊和菲律賓海板塊交界處，南起菲律賓巴拉望島北端，沿菲律賓呂宋西部邊緣向北發(fā)展，北至臺灣島，總長度約1000km。亞歐大陸板塊以70mm／a的速度向菲律賓海板塊之下俯沖，俯沖板塊與上覆板塊之間的匯聚擠壓作用持續(xù)了相當(dāng)長的一段時間，當(dāng)集聚的應(yīng)力釋放，便會發(fā)生海底地震，地震激發(fā)的海嘯波會使臺灣、福建、廣東、海南等沿海地區(qū)毫無遮攔地暴露在海嘯波的威脅之下。

本文通過使用COMCOT模型模擬了馬尼拉海溝地震后產(chǎn)生的海嘯對廣東省附近海域的影響，模型采用2層嵌套網(wǎng)格，網(wǎng)格范圍如圖5所示。美國地質(zhì)調(diào)查局根據(jù)海溝方位角和斷層幾何學(xué)將馬尼拉海溝分為6個斷裂帶，震源位置如圖6所示，地震震級為MW9.3（Phuong et al，2013），震源斷層參數(shù)見表2。

圖7為馬尼拉海溝發(fā)生地震后海嘯波傳播和先導(dǎo)波到達(dá)時間圖。馬尼拉海溝因地震引發(fā)的海嘯波3h之后傳播到廣東省沿岸海域附近。COMCOT可以得到外海島嶼附近海域更精細(xì)的模擬結(jié)果，由COMCOT第2層模型結(jié)果可以得到廣東省附近海嘯引發(fā)的最大增水值。

圖8為馬尼拉海溝地震引發(fā)的海嘯對廣東省附近海域影響的最大增水圖。由圖8可知，馬尼拉海溝處地震引發(fā)的海嘯由深海區(qū)傳入淺水區(qū)后，在淺水效應(yīng)的作用下，在東沙群島附近海嘯波有明顯的增高，最大增水達(dá)到11.95m，海嘯波傳到廣東省沿岸后最大增水約為4.87m。海嘯波的速度近似為（g H）1／2，當(dāng)海嘯波與高潮位重合時，耦合計算總水深H增大，海嘯傳播速度增大，傳播時間變短。若不考慮潮汐、海嘯的計算結(jié)果，在時間上會有所延遲，約15min左右，沿海各地受地形的影響，會略有不同。表3為在馬尼拉海溝地震影響下廣東各海洋站增、減水值及對應(yīng)時間。

圖6 馬尼拉海溝斷裂帶位置

表2 馬尼拉海溝斷層參數(shù)

為了更加詳細(xì)地了解海嘯波對廣東省沿岸的影響，選取廣東省沿岸的海門、甲子、汕尾、惠州港、香港、桂山島、廣州、大萬山、珠海港、海陵山島（閘坡）、湛江、下泊等處，得到其水位隨時間的變化序列，從而分析海嘯波對廣東省沿岸海域的影響。敏感點(diǎn)位置見圖9，受影響的時間序列見圖10。在各海洋站的統(tǒng)計中，東沙島和大萬山的最大增水相對較為顯著，分別為11.95、4.87m；其次是在珠海港和甲子附近，最大增水達(dá)到4m左右，惠州港附近3.76m，海門、汕尾、井岸、香港最大增水為2.0～2.5m，海陵山島（閘坡）、汕頭、茂石化（水東港）、東澳島、桂山島、西葛、潮州港（三百門）最大增水為1.5～2.0m，大鵬灣（鹽田港）、大亞灣、湛江、南澳島（云澳灣）、馬鞭洲（廣石化）、橫山、三灶島、燈籠山、流沙最大增水為1.0～1.5m；隨后較弱的為烏石、下泊、下港、博賀、電白、硵洲島（北港）、北津、海安，最大增水為0.5～1.0m；澳門、珠海（九洲港）、內(nèi)伶仃島、黃埔、橫門、珠海（香洲）、北街、南沙（水牛頭）、廣州、海沁、深圳機(jī)場（油碼頭）、上川島（三洲灣）、舢舨洲、蛇口（赤灣）最大增水為0.19～0.50m，相對稍弱。由圖10可知，馬尼拉海溝地震引發(fā)的海嘯約2.5h后到達(dá)廣東沿岸，海嘯波到達(dá)敏感點(diǎn)區(qū)附近后，影響在1h內(nèi)達(dá)到最大，隨后影響逐漸變小。內(nèi)伶仃洋因外海星羅棋布的島嶼阻攔和淺水地形影響，受海嘯波的影響微弱，海嘯波高相對較小。不到5h，最大增水覆蓋了廣東大部分沿海地區(qū)，而對于廣州、湛江以西的影響則較為微弱。

圖7 馬尼拉海溝地震海嘯波傳播及先導(dǎo)波到達(dá)時間

圖8 馬尼拉海溝地震引發(fā)的海嘯對廣東省沿岸海域最大增水圖

表3 在馬尼拉海溝地震影響下廣東各海洋站增、減水值及對應(yīng)時間

圖9 站位示意圖

圖10 馬尼拉海溝地震后典型站位水位時間序列

4 結(jié)語

本文基于COMCOT模型建立了廣東沿海天文潮與海嘯耦合數(shù)學(xué)模型，對日本“311”海嘯進(jìn)行了良好地驗(yàn)證，同時模擬了馬尼拉海溝地震后產(chǎn)生的海嘯對廣東附近海域的影響，對廣東沿岸由東向西的海洋站（海門至下泊等處）數(shù)據(jù)進(jìn)行了分析，主要研究結(jié)果如下。

（1）本文建立的海嘯耦合數(shù)學(xué)模型能夠較好地模擬日本“311”海嘯，首波到岸時間的模擬值與觀測值之差均小于2min，波高的計算結(jié)果與實(shí)測值對比良好，兩者的差異與海底地形誤差、驗(yàn)潮站局部海浪初值累積誤差等有關(guān)。本模型為準(zhǔn)確評估馬尼拉海溝地震后產(chǎn)生的海嘯對廣東沿海的影響提供了研究基礎(chǔ)。

（2）馬尼拉海溝地震引發(fā)的海嘯約2.5h后到達(dá)廣東沿岸，海嘯波到達(dá)敏感點(diǎn)區(qū)附近后，影響在1h內(nèi)達(dá)到最大，隨后影響逐漸變小。內(nèi)伶仃洋因外海星羅棋布的島嶼阻攔和淺水地形影響，受海嘯波的影響微弱，海嘯波高相對較小。

（3）基本上5h以內(nèi)，最大增水覆蓋了廣東大部分沿海地區(qū)，對近海工程將造成影響。東沙島和大萬山的最大增水相對較為顯著，分別為11.95、4.87m；其次是在珠海港和甲子附近最大增水達(dá)到4m左右，惠州港附近3.76m，海門、汕尾、井岸、香港最大增水為2.0～2.5m。而對于廣州，湛江以西的影響則較為微弱。

（4）由于南海海嘯發(fā)生的歷史記錄尚少，地震數(shù)據(jù)亦有所欠缺，震源機(jī)制解不能有效精確，故計算時間段無法準(zhǔn)確判定。線性疊加潮汐則比耦合潮汐的海嘯波計算結(jié)果稍安全。同時，潮汐對海嘯的影響敘述篇幅太長，這也是后續(xù)的研究工作，海嘯波抵達(dá)時間并無多大差異，主要差異體現(xiàn)在海嘯最大波高，將來可考慮多年最大、最小潮差時間段，一年當(dāng)中不同月份的最大、最小潮差時間段，以及受到東北季風(fēng)和西南季風(fēng)影響，甚至受到極端天氣影響情況下的模擬，以構(gòu)建相應(yīng)的海嘯源數(shù)據(jù)庫，建立更加完善的南海海嘯耦合數(shù)學(xué)模型。