基于全和諧型淺水方程的有限體積海嘯數(shù)值模型開發(fā)與應(yīng)用

周文，王培濤，王崗*，于福江，鄭金海，梁秋華

( 1. 海岸災(zāi)害及防護(hù)教育部重點(diǎn)實(shí)驗(yàn)室（河海大學(xué)），江蘇 南京 210098；2. 國(guó)家海洋環(huán)境預(yù)報(bào)中心，北京 100081；3. 拉夫堡大學(xué) 建筑與土木工程學(xué)院，英國(guó) 萊斯特郡 LE11 3TU)

1 引言

海嘯作為最具危險(xiǎn)性和破壞力的海洋災(zāi)害之一，給沿海地區(qū)造成了巨大的經(jīng)濟(jì)損失和人員傷亡。如1960 年5 月22 日智利瓦爾迪維亞（Valdivia）發(fā)生9.5 級(jí)大地震，這也是人類歷史上用儀器監(jiān)測(cè)到的最強(qiáng)地震，隨之引發(fā)的巨大海嘯，在當(dāng)?shù)氐暮[波波高達(dá)25 m，對(duì)附近一個(gè)居住著250 萬(wàn)人口的城鎮(zhèn)造成約1 萬(wàn)人死亡、500 億～700 億美元的經(jīng)濟(jì)損失。此外，該海嘯還影響了整個(gè)太平洋海域，導(dǎo)致了太平洋另一端的日本數(shù)百人喪生[1]。

為了應(yīng)對(duì)海嘯災(zāi)害，人們?cè)诤[監(jiān)測(cè)、預(yù)報(bào)和致災(zāi)機(jī)理等方面做了大量研究。由于人類史上災(zāi)難性的海嘯主要是由地震引起的，因此人們通過鋪設(shè)海底地震儀、投放海嘯浮標(biāo)和設(shè)置沿岸驗(yàn)潮站等建立海嘯監(jiān)測(cè)網(wǎng)。地震儀和觀測(cè)海底變形的大地傳感器等通過海底線纜將數(shù)據(jù)傳輸給陸上監(jiān)測(cè)中心。海嘯浮標(biāo)系統(tǒng)主要是美國(guó)國(guó)家海洋與大氣管理局（NOAA）下屬太平洋海洋環(huán)境實(shí)驗(yàn)室研制的海嘯深海評(píng)估和報(bào)告浮標(biāo)系統(tǒng)（Deep-ocean Assessment and Reporting of Tsunamis, DART），它是由海底壓力記錄儀、海表浮標(biāo)和地球同步氣象衛(wèi)星構(gòu)成的立體監(jiān)測(cè)系統(tǒng)[2]，目前已更新至第4 代。目前全球海域設(shè)置了60 多個(gè)DART和900 余個(gè)驗(yàn)潮站[3]，覆蓋全面的大洋與近岸監(jiān)測(cè)網(wǎng)為海嘯預(yù)警以及海嘯研究提供了大量可靠的原始數(shù)據(jù)。

隨著計(jì)算機(jī)和數(shù)值計(jì)算技術(shù)的發(fā)展，數(shù)值模擬已成為海嘯的主要研究方法并在海嘯預(yù)警系統(tǒng)中起到了關(guān)鍵作用。目前常見的海嘯模型可主要分為基于淺水方程的靜壓模型和考慮波浪動(dòng)水壓力的數(shù)值模型兩類。第1 類是基于淺水方程的數(shù)值模型，由于其計(jì)算效率高且在大多數(shù)模擬過程均能得到較好的結(jié)果，應(yīng)用最為廣泛。華盛頓大學(xué)LeVeque 教授團(tuán)隊(duì)在1994 年發(fā)布了用于求解雙曲型方程的軟件包Clawpack（Conservation Laws Package），且其版本一直在不斷完善與更新[4]。2010 年GeoClaw[5]作為Clawpack求解淺水方程的模塊用于海嘯模擬。GeoClaw 包含笛卡爾坐標(biāo)和球坐標(biāo)兩種形式，基于高分辨率激波捕捉的有限體積法進(jìn)行數(shù)值求解，同時(shí)采用自適應(yīng)網(wǎng)格嵌套技術(shù)并考慮了干濕邊界的模擬。第2 類常見的海嘯模型主要是基于Boussinesq 方程，考慮了波浪的頻散性。美國(guó)特拉華大學(xué)建立了基于完全非線性Boussinesq 型方程的波浪模型—FUNWAVE[6]。該模型此后進(jìn)一步更新為TVD[7]及GPU[8]版本。Yamazaki等[9]建立了直接求解Navier-Stokes 方程的非靜壓波浪模型模擬海嘯的傳播演化過程。Zhao 等[10]基于Green-Naghdi 理論建立了考慮海底動(dòng)態(tài)可變的非靜壓波浪模型模擬二維和三維的海嘯傳播過程。

雖然第2 類海嘯模型有更高的數(shù)值精度，但考慮到海嘯預(yù)報(bào)的時(shí)效性，目前海嘯模型中仍以淺水方程數(shù)值模型居多。早期的部分海嘯模型如MOST[11]等采用有限差分法（FDM）求解控制方程，其優(yōu)點(diǎn)是簡(jiǎn)單成熟，方便構(gòu)造高精度的計(jì)算格式，然而處理復(fù)雜網(wǎng)格時(shí)不夠靈活[12]。由于有限體積法（FVM）數(shù)學(xué)表達(dá)直觀、物理涵義清晰、滿足守恒規(guī)律且能方便地捕捉淺水方程中的間斷解[12]，已逐漸成為一種重要的數(shù)值模型方法。淺水方程的有限體積數(shù)值求解主要涉及源項(xiàng)與通量的處理以及動(dòng)邊界（干濕邊界）的準(zhǔn)確模擬。關(guān)于源項(xiàng)與通量的處理目前主要是通過不同的數(shù)值方法用以平衡有限體積法中的通量梯度和源項(xiàng)[13-15]。而動(dòng)邊界的準(zhǔn)確模擬目前主要基于Riemann近似解。HLL（Harten, Lax and van Leer）Riemann 近似求解器[16]在1983 年提出，它根據(jù)有限體積單元界面處左右波速的預(yù)測(cè)值將界面通量值的計(jì)算劃分為3 個(gè)區(qū)域，能夠很好地解決一維淺水流動(dòng)問題。在此基礎(chǔ)上，HLLC Riemann 近似求解器[17]考慮了二維問題中另一個(gè)方向的橫波影響，將單元界面通量值的計(jì)算劃分為4 個(gè)區(qū)域，用以簡(jiǎn)便和精確地求解二維淺水波動(dòng)問題。為了合理簡(jiǎn)潔地模擬動(dòng)邊界，Liang 和Borthwick[18]提出了全和諧型淺水方程，以自動(dòng)滿足控制方程通量與源項(xiàng)的守恒，并通過對(duì)地形重構(gòu)抑制單元內(nèi)的虛擬流動(dòng)，開發(fā)了一種高性能強(qiáng)穩(wěn)定性的淺水?dāng)?shù)值模型。該模型數(shù)值方法操作簡(jiǎn)便且自動(dòng)滿足單元體積各參量守恒，已被包括FUNWAVE 等數(shù)值模型廣泛采用[7]。

由于災(zāi)害性海嘯的影響通常波及整個(gè)大洋甚至全球范圍，本文基于球坐標(biāo)系下的全和諧型淺水方程建立了有限體積海嘯數(shù)值模型?？紤]到海嘯不僅會(huì)對(duì)橋墩、防波堤等近海建筑物產(chǎn)生影響，也會(huì)造成近岸大面積的陸地淹水，因此海嘯近岸爬高的精確模擬也是海嘯數(shù)值模擬的重要技術(shù)要求。沿海地區(qū)往往是經(jīng)濟(jì)較為發(fā)達(dá)的區(qū)域，不僅存在著復(fù)雜的自然地形還涉及大量的人工建筑物。本文采用地形重構(gòu)的方法，以模擬海嘯在此復(fù)雜地形上的近岸爬高過程。本文所建立的海嘯數(shù)值模型應(yīng)用于2015 年9 月16 日智利Mw8.3 級(jí)地震海嘯的模擬，以驗(yàn)證模型對(duì)越洋海嘯模擬預(yù)報(bào)的能力。

2 控制方程

越洋海嘯的傳播需要考慮地球的曲率和地轉(zhuǎn)偏向力的影響，因此控制方程選取考慮科氏力的球坐標(biāo)系下的非線性淺水方程，其矩陣形式為[19]

式中，θ和φ分別為地球的經(jīng)度和緯度；t為時(shí)間；向量U為變量；F和G分別是θ和φ方向的通量；S為源項(xiàng)；具體為

式中，r為地球半徑；ξ為相對(duì)于靜止海平面的波面位移；hs為靜水深；h=ξ+hs為總水深；uθ和uφ分別為θ和φ方向沿水深平均的速度。源項(xiàng)中第2、第3 項(xiàng)分別對(duì)應(yīng)科氏力和摩阻項(xiàng)，其中f=2ωsinφ為科氏力系數(shù)；ω為地球繞地軸自轉(zhuǎn)的角速度；Cf=gn2/h1/3為底摩阻系數(shù)，n是Manning 系數(shù)。

為了平衡通量和源項(xiàng)，本文利用Liang 和Borthwick[18]提出的全和諧型淺水方程表達(dá)形式改寫源項(xiàng)式（5）中波面的空間導(dǎo)數(shù)。以θ方向?yàn)槔?，?/p>

式中，zb為基準(zhǔn)線上的海底地形高程；η=zb+h為波面高程。

在此不考慮zb和hs隨時(shí)間t改變，從而有

因此式（2）至式（5）的最終形式為

式（9）至式（10）將地形梯度分離至方程右側(cè)的源項(xiàng)中用以平衡方程左側(cè)通量的變化，從數(shù)學(xué)表達(dá)上保證了方程的守恒形式。

此外，為方便模型用于小范圍特別是物理實(shí)驗(yàn)的模擬，本文在建立球坐標(biāo)模型的同時(shí)也開發(fā)了笛卡爾坐標(biāo)海嘯模型。由于兩套方程的數(shù)值求解方法完全一致，后文將統(tǒng)一介紹。

3 數(shù)值方法

本文基于均勻矩形空間單元，利用Godunov 格式的有限體積法[18]求解方程矩陣?；贛USCL-Hancock 格式，采用HLLC Riemann 近似求解器[17]計(jì)算單元界面上的流體通量，保證數(shù)值格式在時(shí)間和空間上均為二階精度。

3.1 MUSCL-Hancock 格式

3.1.1 預(yù)測(cè)步

首先在時(shí)間上求解半個(gè)時(shí)間步（Δt/2）對(duì)應(yīng)的流體變量預(yù)測(cè)結(jié)果

式中，上標(biāo)n代表第n時(shí)刻；i、j為有限體積單元標(biāo)號(hào)；Δθ、Δφ表示球坐標(biāo)系水平方向的單元大??；fe、fw，gn、gs分別代表單元東、西、北、南4 個(gè)方向的流體界面通量。單元界面的流體變量利用MUSCL 線性重構(gòu)計(jì)算[20]，以單元東西界面為例

此外，源項(xiàng)中海底地形底坡項(xiàng)直接利用中心差分法離散，科氏力項(xiàng)直接代入計(jì)算即可，而摩阻項(xiàng)后文詳細(xì)介紹。

3.1.2 校正步

在校正步將變量遞進(jìn)至完整的時(shí)間Δt，具體為

式中，F(xiàn)E、FW、GN、GS分別為求解Riemann 問題得到的單元四周的流體界面通量，單元界面處Riemann 變量的求解依然利用MUSCL 線性重構(gòu)[20]。

3.2 干濕邊界

為了準(zhǔn)確模擬波浪在陸地和海洋邊界處的運(yùn)動(dòng)，本文采用Liang[21]提出的地形重構(gòu)方法。以單元東界面為例，為了確保計(jì)算過程中單元的水深恒為正，則地形重構(gòu)為

式中，上標(biāo)L 和R 分別表示界面的左右側(cè)。對(duì)應(yīng)的水深為

綜合式(17)至式(19)可求得Riemann 變量η和q。

為了使干濕邊界的處理不失一般性，考慮如圖1干濕單元相鄰的情形。單元(i,j)的東界面通量的計(jì)算仍以地形為基線，而西界面的計(jì)算基線為水面線，假設(shè)水面靜止，東側(cè)地形高于西側(cè)則水面不會(huì)發(fā)生流動(dòng)，而數(shù)值計(jì)算時(shí)單元(i,j)內(nèi)會(huì)引入從東至西的通量。為了避免虛假通量的引入，進(jìn)一步修正zb與η為

圖1 地形重構(gòu)示意Fig. 1 Local topography reconstruction

其中

經(jīng)過上述修正，東西界面的計(jì)算基線得到統(tǒng)一，保證了通量計(jì)算的平衡性。

3.3 底摩阻

在利用MUSCL-Hancock 格式數(shù)值求解后進(jìn)一步利用分裂格式(splitting point-implicit scheme)對(duì)摩阻項(xiàng)離散[22]，其等價(jià)于求解以下常微分方程

式中，Sf為摩阻項(xiàng)，單寬流量q=[huθ huφ]T，其隨時(shí)間演進(jìn)的計(jì)算為

式中，D為隱式系數(shù)，設(shè)F=(Sf/D)n。為了確保干濕邊界附近數(shù)值計(jì)算的穩(wěn)定性，對(duì)F設(shè)置閾值

其物理意義是摩阻作用至多會(huì)使流體停止運(yùn)動(dòng)而不產(chǎn)生逆向流動(dòng)。

3.4 邊界條件和穩(wěn)定性準(zhǔn)則

如圖2 以某個(gè)方向第一個(gè)單元為例，邊界處虛擬單元的變量利用相鄰單元的值近似處理，單寬流量q的取值根據(jù)邊界類型而變化。對(duì)于開放邊界，有

圖2 邊界條件示意圖Fig. 2 Sketch of boundary condition

對(duì)于閉合邊界、即全反射邊界，有

為了保證數(shù)值計(jì)算的穩(wěn)定，時(shí)間步Δt的確定遵循Courant-Friedrichs-Lewy (CFL)準(zhǔn)則

虛擬網(wǎng)格

式中，C為Courant 系數(shù)；ΔSθ和ΔSφ分別為經(jīng)度與緯度方向的單元空間長(zhǎng)度。

3.5 海嘯源模型

由于歷史上災(zāi)害性海嘯幾乎都是由地震所引起的，所以海嘯源模型主要考慮地震引起的海水波動(dòng)過程，它是海嘯數(shù)值計(jì)算的基礎(chǔ)，其適用性直接影響海嘯波的傳播及海嘯與海岸間相互作用。自Steketee[23-24]把位錯(cuò)理論引入地震學(xué)后，許多地震學(xué)者針對(duì)不同的斷層類型發(fā)展了不同的位錯(cuò)理論，使得計(jì)算地球同震形變的方法得到不斷發(fā)展。Kajiura[25]首次提出可以將靜態(tài)的海底位移轉(zhuǎn)化成海嘯產(chǎn)生階段自由表面的初始邊界條件，使得地震海嘯數(shù)值模擬過程進(jìn)一步簡(jiǎn)化。為此，利用斷層模型計(jì)算地震引發(fā)的海床位移量，從而估算出海底形變引起初始水面高度來初始化海嘯數(shù)值模型已成為海嘯模擬的慣用方法。目前，Mansinha 和Smylie[26]以及Okada[27]基于彈性錯(cuò)移理論發(fā)展的兩套斷層模型被廣泛應(yīng)用，大量的研究和應(yīng)用實(shí)例表明此類模型對(duì)逆沖斷層地震海嘯源計(jì)算具有較好的適用性[28-30]。Okada 模型[27]計(jì)算海底同震位移需要輸入7 個(gè)震源參數(shù)，分別為震中位置、震源深度、斷層走向角、傾角、滑動(dòng)角、滑移量以及斷層破裂的尺度（斷層長(zhǎng)度及寬度）。震源基本信息（震中位置、震源深度）可由地震速報(bào)系統(tǒng)測(cè)定，而震源機(jī)制角可由W-phase 反演系統(tǒng)得到。震源的尺度和滑移量可參考震級(jí)與破裂尺度經(jīng)驗(yàn)關(guān)系率估計(jì)[31]。需要說明的是，上述方法在海底較為平緩且為典型的逆沖斷層條件下較為適用。當(dāng)?shù)卣鹨鸲钙麓蠓秶乃竭\(yùn)動(dòng)或者地震引起斷層水平運(yùn)動(dòng)比垂向運(yùn)動(dòng)大得多的走滑斷層時(shí)應(yīng)該考慮水平位移在海嘯產(chǎn)生過程中的貢獻(xiàn)。

4 2015 年智利地震海嘯數(shù)值模擬

2015 年9 月16 日22:54:32（UTC）位于納斯卡板塊和南美洲板塊之間的南美洲俯沖帶發(fā)生Mw8.3 級(jí)地震。震中位于31.57°S, 71.67°W，距離智利伊亞佩爾市以西48 km，震源深度為22.4 km。該地震引發(fā)了太平洋范圍的中等強(qiáng)度海嘯，造成近場(chǎng)地區(qū)30 多人傷亡，經(jīng)濟(jì)損失高達(dá)6 億美元[32]。

4.1 海嘯產(chǎn)生階段模擬

考慮海嘯預(yù)警業(yè)務(wù)化需求，通常假定海底地震具有一致震源機(jī)制特征，采用點(diǎn)源、單斷層面和平均滑動(dòng)量的假定以簡(jiǎn)化斷層破裂的復(fù)雜性，忽略了斷層破裂的局地特征以便最短時(shí)間內(nèi)給出相對(duì)合理的震源破裂特征。自Kanamori 和Rivera[33]提出基于W 震相測(cè)定地震矩張量的方法后，全球多個(gè)預(yù)警中心基于W 震相研制出了中強(qiáng)震矩心矩張量反演系統(tǒng)。基本實(shí)現(xiàn)在震后8～15 min 準(zhǔn)確獲取全球Mw≥6.5 地震的震源機(jī)制解，從而有效提高了海嘯定量預(yù)警及時(shí)性。表1 為2015 年智利Mw8.3 級(jí)地震斷層基本參數(shù)，其中地震斷層破裂尺度和滑移量依據(jù)Wu 等[31]關(guān)系率估算所得，即地震斷層長(zhǎng)度約為212 km，寬度約為79 km，平均滑移量6.3 m。斷層尺度量級(jí)與記錄的深海海嘯波長(zhǎng)對(duì)比相當(dāng)，估算的結(jié)果可信。選擇USGS 快速W-phase 解作為該事件的震源機(jī)制解，從震源機(jī)制解可知該地震斷層為典型的低傾角逆沖斷層，因此可利用Okada 彈性位錯(cuò)理論模型[27]計(jì)算海嘯傳播模型所需的海表面初始形變（圖3）。從海表面初始形變分布可知主要的形變呈輕微NNE 向的扁狀橢圓形分布，緯度方向跨度約2°，經(jīng)度方向跨度約0.5°。

圖3 地震引起的海表面初始變形Fig. 3 Initial sea surface deformation induced by the earthquake

表1 智利地震斷層參數(shù)Table 1 Fault parameters of Chile earthquake

4.2 智利海域區(qū)域海嘯模擬

智利海域區(qū)域海嘯模擬范圍為60°S～5°N、65°～110°W，地形數(shù)據(jù)來自ETOPO1 (https://www.ngdc.noaa.gov/mgg/global/global.html)，其中以平均海平面為基準(zhǔn)，海域水深為負(fù)、陸地為正。模型采用分辨率為1′的均勻網(wǎng)格，時(shí)間步長(zhǎng)Δt=2 s，整個(gè)模擬過程為8 h（圖4）。Manning 系數(shù)在水深小于100 m 時(shí)取0.025，其余則取0。計(jì)算域四周設(shè)置為開放邊界。本文模擬采用CPU 為Intel(R) Core(TM) i7-7500U 的筆記本電腦，整個(gè)模擬耗時(shí)50.5 h。

為了驗(yàn)證結(jié)果，本文選取了14 個(gè)智利近岸測(cè)站的實(shí)測(cè)數(shù)據(jù)（各測(cè)站位置見圖4）。模擬結(jié)果與測(cè)站實(shí)測(cè)比較見圖5，各測(cè)站第一個(gè)測(cè)到的海嘯先導(dǎo)波時(shí)間與波幅比較見表2。模擬與實(shí)測(cè)的海嘯先導(dǎo)波到達(dá)時(shí)間及波幅均吻合良好。智利海域的海嘯持續(xù)時(shí)間較長(zhǎng)，且實(shí)測(cè)相較于模擬結(jié)果呈現(xiàn)出更多的微小波動(dòng)現(xiàn)象，這可能是由于近岸劇烈變化的地形及岸線影響。此外，COQUIMBO 測(cè)站所測(cè)最大波高達(dá)4.75 m，最大波不是先導(dǎo)波，且先導(dǎo)波后的幾個(gè)波幅基本超過2.5 m，大于模擬結(jié)果。其可能的原因一方面是該測(cè)站位于震源附近，受復(fù)雜的地震影響，概化的海嘯源模型沒有充分反演海嘯產(chǎn)生階段一系列復(fù)雜的過程。另一方面是該測(cè)站所處位置地形變化劇烈，模型采用的網(wǎng)格未能精細(xì)刻畫地形。

圖4 智利區(qū)域海嘯和泛太平洋海嘯模擬范圍及DART 浮標(biāo)與智利近岸測(cè)站位置Fig. 4 Domains of Chile regional tsunami and the Pacific tsunami and the location of DART buoys and coastal tide-gauge stations

圖5 智利近岸測(cè)站波面過程與模擬對(duì)比Fig. 5 Time histories of observation and simulation at Chile coastal tide-gauge stations

表2 智利近岸測(cè)站海嘯先導(dǎo)波到達(dá)時(shí)間和波幅模擬與實(shí)測(cè)比較Table 2 Comparisons of simulation and observation of the arrival time and amplitude for the leading wave at Chile coastal tidal-gauge stations

4.3 泛太平洋海嘯模擬

泛太平洋海嘯模擬區(qū)域?yàn)?5°S～70°N，100°E～60°W，地形數(shù)據(jù)同樣來自ETOPO1，采用空間分辨率為4′的均勻網(wǎng)格。Manning 系數(shù)與4.2 節(jié)區(qū)域海嘯相同，計(jì)算域四周均為開放邊界。模擬時(shí)間步長(zhǎng)Δt=5 s，整個(gè)模擬過程為25 h，同樣采用CPU 為Intel（R）Core（TM）i7-7500U 的筆記本電腦，用時(shí)28.0 h。

地震發(fā)生后不同時(shí)刻的瞬時(shí)波面過程見圖6，第一個(gè)出現(xiàn)的海嘯先導(dǎo)波到達(dá)時(shí)間見圖7。海嘯傳播初期基本以震中為圓心向四周擴(kuò)散（圖6a）。當(dāng)海嘯經(jīng)過海脊和海底山等地形時(shí)，由于其在水深較淺的地區(qū)傳播速度變慢，導(dǎo)致海嘯波峰線呈現(xiàn)V 形，即在海脊和海底山附近的海嘯波傳播較慢而在深水處傳播較快。由于太平洋NE 側(cè)水深較SW 一側(cè)深，所以海嘯波峰線逐漸從NE-SW 方向變化至SN 向（圖6d），其對(duì)應(yīng)的先導(dǎo)波到達(dá)時(shí)間也呈現(xiàn)出NE 側(cè)快于SW 側(cè)的特點(diǎn)（圖7）。此外，太平洋SW 側(cè)存在眾多海脊、海底山及海島，導(dǎo)致此處的海嘯波較NE 側(cè)呈更加復(fù)雜的波面過程。

圖6 地震后越洋海嘯瞬時(shí)波面分布Fig. 6 Snapshots of simulated tsunami wavefields after the earthquake

圖7 海嘯先導(dǎo)波到達(dá)時(shí)間等值線圖Fig. 7 Contours of the arrival time for the leading tsunami wave

進(jìn)一步選取20 個(gè)分布于環(huán)太平洋沿岸及位于太平洋中部的DART 浮標(biāo)（見圖4）驗(yàn)證本文所建立的數(shù)值模型，其中波面過程對(duì)比見圖8、海嘯先導(dǎo)波的到達(dá)時(shí)間和波幅對(duì)比見表3。由于從智利至太平洋東北側(cè)的海域地形相對(duì)規(guī)則，海嘯最大波通常為其先導(dǎo)波，且波形規(guī)則、持續(xù)時(shí)間較短。該海域的海嘯模擬結(jié)果與實(shí)測(cè)吻合較好。相比于東北側(cè)，太平洋西南側(cè)的地形復(fù)雜，存在著海脊、海底山和海島等地形。海嘯在此受海底地形折射、繞射等影響，自由水面呈較為復(fù)雜的波動(dòng)過程，海嘯最大波不再是先導(dǎo)波。該海域本文模擬的前幾個(gè)波動(dòng)過程與實(shí)測(cè)基本吻合，但后續(xù)的波動(dòng)過程無論在振幅還是相位上均存在一定差別。這可能是由于模型采用的網(wǎng)格未能精細(xì)刻畫出該區(qū)域復(fù)雜的地形所致。

表3 DART 浮標(biāo)海嘯先導(dǎo)波到達(dá)時(shí)間和波幅模擬與實(shí)測(cè)比較Table 3 Comparisons of simulation and observation of the arrival time and amplitude for the leading wave at DART buoys

圖8 DART 海嘯浮標(biāo)實(shí)測(cè)波面過程與模擬結(jié)果對(duì)比Fig. 8 Comparisons of the free surface elevation between DART buoys and simulations

本文模擬的海嘯先導(dǎo)波到達(dá)時(shí)間與實(shí)測(cè)結(jié)果隨著傳播距離增加其誤差也逐漸增加，這可能是由于海嘯波的頻散效應(yīng)所致。Glimsdal 等[34]提出一個(gè)與傳播距離和時(shí)間呈正相關(guān)的“頻散時(shí)間”以量化海嘯的頻散效應(yīng)。這與本文模擬結(jié)果所顯示的隨著海嘯傳播距離和時(shí)間的增加先導(dǎo)波到達(dá)時(shí)間的誤差逐漸增加是一致的。此外，實(shí)測(cè)資料顯示在海嘯較大波動(dòng)到達(dá)之前會(huì)有一段幅度較小的減水過程，Watada 等[35]認(rèn)為其與海底的彈性形變有關(guān)。由于本文的模型沒有考慮該因素的影響，未能真實(shí)反演這一物理現(xiàn)象。

圖9 為智利海嘯在太平洋的最大波幅分布。在智利附近海域，海嘯最大波幅主要受斷層走向影響，呈現(xiàn)出由智利沿岸向西擴(kuò)散的趨勢(shì)。在距離震源較遠(yuǎn)的區(qū)域，海嘯最大波幅呈現(xiàn)由智利向日本延伸的趨勢(shì)，且基本與海脊走向一致，這可能是海嘯在海脊上形成所謂的海脊俘獲波所致[36]。

圖9 智利海嘯在太平洋的最大波幅分布Fig. 9 Distribution of maximum wave amplitude in the Pacific for Chile tsunami

5 結(jié)論

本文基于球坐標(biāo)系下的全和諧型淺水方程，采用MUSCL-Hancock 格式，并利用HLLC Riemann 近似求解器計(jì)算單元界面上的流體通量，最終建立了二階精度的Godunov 格式有限體積海嘯數(shù)值模型。由于近岸局部地區(qū)海嘯淹沒爬坡的精細(xì)化模擬和受災(zāi)評(píng)估中必將涉及到干濕邊界的處理；同時(shí)，良好的干濕邊界處理還有益于維持越洋海嘯模擬的穩(wěn)定性，因此，本文建立的模型包含了處理干濕邊界的功能。在此基礎(chǔ)上，本文利用彈性半空間位錯(cuò)理論得到的初始自由水面模擬了2015 年9 月16 日智利Mw8.3 級(jí)地震引發(fā)的區(qū)域海洋和泛太平洋海嘯傳播過程。分別通過與智利近岸14 個(gè)測(cè)站和環(huán)太平洋20 個(gè)DART 浮標(biāo)實(shí)測(cè)數(shù)據(jù)比較，進(jìn)一步驗(yàn)證了模型模擬實(shí)際海嘯過程的能力。

本文泛太平洋海嘯的模擬采用空間分辨率為

4′的均勻網(wǎng)格，在開闊海域和水深變化較緩的區(qū)域模擬結(jié)果與實(shí)測(cè)吻合良好，然而在一些島嶼及近岸受復(fù)雜地形及岸線影響的區(qū)域有一定誤差，主要是由于模擬所采用的網(wǎng)格空間分辨率未能精細(xì)刻畫這些因素的影響。為了兼顧實(shí)際模擬中的精度與效率，今后將在控制方程中考慮海水的壓縮性、地震導(dǎo)致的海底彈性形變等因素的影響并結(jié)合多層嵌套與并行技術(shù)進(jìn)一步提高模型的模擬能力。