基于全和諧型淺水方程的有限體積海嘯數值模型開發與應用

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

( 1. 海岸災害及防護教育部重點實驗室（河海大學），江蘇 南京 210098；2. 國家海洋環境預報中心，北京 100081；3. 拉夫堡大學 建筑與土木工程學院，英國 萊斯特郡 LE11 3TU)

1 引言

海嘯作為最具危險性和破壞力的海洋災害之一，給沿海地區造成了巨大的經濟損失和人員傷亡。如1960 年5 月22 日智利瓦爾迪維亞（Valdivia）發生9.5 級大地震，這也是人類歷史上用儀器監測到的最強地震，隨之引發的巨大海嘯，在當地的海嘯波波高達25 m，對附近一個居住著250 萬人口的城鎮造成約1 萬人死亡、500 億～700 億美元的經濟損失。此外，該海嘯還影響了整個太平洋海域，導致了太平洋另一端的日本數百人喪生[1]。

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

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

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

由于災害性海嘯的影響通常波及整個大洋甚至全球范圍，本文基于球坐標系下的全和諧型淺水方程建立了有限體積海嘯數值模型。考慮到海嘯不僅會對橋墩、防波堤等近海建筑物產生影響，也會造成近岸大面積的陸地淹水，因此海嘯近岸爬高的精確模擬也是海嘯數值模擬的重要技術要求。沿海地區往往是經濟較為發達的區域，不僅存在著復雜的自然地形還涉及大量的人工建筑物。本文采用地形重構的方法，以模擬海嘯在此復雜地形上的近岸爬高過程。本文所建立的海嘯數值模型應用于2015 年9 月16 日智利Mw8.3 級地震海嘯的模擬，以驗證模型對越洋海嘯模擬預報的能力。

2 控制方程

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

式中，θ和φ分別為地球的經度和緯度；t為時間；向量U為變量；F和G分別是θ和φ方向的通量；S為源項；具體為

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

為了平衡通量和源項，本文利用Liang 和Borthwick[18]提出的全和諧型淺水方程表達形式改寫源項式（5）中波面的空間導數。以θ方向為例，有

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

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

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

式（9）至式（10）將地形梯度分離至方程右側的源項中用以平衡方程左側通量的變化，從數學表達上保證了方程的守恒形式。

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

3 數值方法

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

3.1 MUSCL-Hancock 格式

3.1.1 預測步

首先在時間上求解半個時間步（Δt/2）對應的流體變量預測結果

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

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

3.1.2 校正步

在校正步將變量遞進至完整的時間Δt，具體為

式中，FE、FW、GN、GS分別為求解Riemann 問題得到的單元四周的流體界面通量，單元界面處Riemann 變量的求解依然利用MUSCL 線性重構[20]。

3.2 干濕邊界

為了準確模擬波浪在陸地和海洋邊界處的運動，本文采用Liang[21]提出的地形重構方法。以單元東界面為例，為了確保計算過程中單元的水深恒為正，則地形重構為

式中，上標L 和R 分別表示界面的左右側。對應的水深為

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

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

圖1 地形重構示意Fig. 1 Local topography reconstruction

其中

經過上述修正，東西界面的計算基線得到統一，保證了通量計算的平衡性。

3.3 底摩阻

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

式中，Sf為摩阻項，單寬流量q=[huθ huφ]T，其隨時間演進的計算為

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

其物理意義是摩阻作用至多會使流體停止運動而不產生逆向流動。

3.4 邊界條件和穩定性準則

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

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

對于閉合邊界、即全反射邊界，有

為了保證數值計算的穩定，時間步Δt的確定遵循Courant-Friedrichs-Lewy (CFL)準則

虛擬網格

式中，C為Courant 系數；ΔSθ和ΔSφ分別為經度與緯度方向的單元空間長度。

3.5 海嘯源模型

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

4 2015 年智利地震海嘯數值模擬

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

4.1 海嘯產生階段模擬

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

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

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

4.2 智利海域區域海嘯模擬

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

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

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

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

表2 智利近岸測站海嘯先導波到達時間和波幅模擬與實測比較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 泛太平洋海嘯模擬

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

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

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

圖7 海嘯先導波到達時間等值線圖Fig. 7 Contours of the arrival time for the leading tsunami wave

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

表3 DART 浮標海嘯先導波到達時間和波幅模擬與實測比較Table 3 Comparisons of simulation and observation of the arrival time and amplitude for the leading wave at DART buoys

圖8 DART 海嘯浮標實測波面過程與模擬結果對比Fig. 8 Comparisons of the free surface elevation between DART buoys and simulations

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

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

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

5 結論

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

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

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