基于GeoClaw海嘯數值模擬程序對2011年日本東北大地震有限斷層模型的驗證

胡 培, 景惠敏*, 楊 煒

基于GeoClaw海嘯數值模擬程序對2011年日本東北大地震有限斷層模型的驗證

胡 培1, 景惠敏1*, 楊 煒2

（1.寧波大學 機械工程與力學學院, 浙江 寧波 315211; 2.寧波大學科學技術學院, 浙江 寧波 315300）

由于海嘯波具有長周期、慢衰減的特征, 對于能夠激發海嘯的海底地震, 海嘯波資料一定程度上可以作為利用地震波資料和大地測量資料反演得到的震源有限斷層模型的補充數據, 以驗證該模型的合理性. 本文以2011年日本東北大地震海嘯期間反演的震源模型為例, 利用Okada彈性半無限空間理論將震源有限斷層模型轉化為垂向海底變形; 然后利用GeoClaw海嘯數值模型進行海嘯數值模擬, 在計算區域模擬實際觀測點的海嘯波波形數據; 并將模擬海嘯波波形數據與實際觀測值進行對比, 以驗證震源有限斷層模型的可靠性.

日本東北大地震; 有限斷層模型; GeoClaw數值模擬

2011年3月11日(UTC), 在日本東北部發生了M9.0級地震[1], 引發了海嘯災難, 造成大量房屋沖毀, 近2萬人失蹤與死亡. 該地震是日本記錄以來震級最大的一次地震[2]. 對于此次發震機制和破裂特征, 國際各地震研究所利用觀測網記錄的各種觀測數據進行研究, 根據不同的數據得到了各自的有限斷層模型[3-7], 結果顯示各斷層模型差異明顯. 目前, 震源機制測定可根據地震波資料和大地測量資料反演來實現[8]. 但對于特大地震而言, 這兩種資料反演的工作都存在不同程度上的欠缺[9]. 對于激發海嘯的特大地震, 其海嘯波具有長周期、慢衰減的特點, 從而可在一定程度上佐證各有限斷層模型的合理性. 一般來說, 我們可以借助成熟的海嘯數值模擬程序計算出海底變形引起的海嘯波[10], 將模擬波形與各觀測站觀測數據進行對比, 進而來論述各有限斷層模型的合理性[11].

在海嘯數值模擬方面, 國際上主要采用MOST[12]、TUNAMI-N1[13]、COMCOT[14]等數值模型進行海嘯數值模擬, 此類數值模擬均采用有限差分法求解二維淺水方程. 有限差分法通過泰勒展開式來構造微分方程中各階導數, 利用差分替代微商, 其離散得到的代數方程沒有實際的物理意義, 僅在數值上接近. 另外, 在進行網格計算時, 會產生數值頻散現象影響海嘯模擬結果.

近年來, GeoClaw模型[15]開始廣泛應用于海嘯的數值模擬, 該模型使用有限體積法求解非線性淺水方程. 有限體積法對控制方程進行積分, 然后采用高斯定理將擴散項和對流項的體積積分轉化為面積分, 得到的方程能夠表示控制容積內物理量的守恒. 這種守恒的特性在處理有間斷問題時(如海嘯)依舊具有合理性. 同時, 該模型使用自適應網格細化技術, 通過預先設定的閾值, 確定需要進行細化的計算區域, 將細化的網格嵌套在粗網格上, 網格層層嵌套, 對空間與時間的步長進行細化. 相比于其他數值模型的網格處理, 具有較高的計算效率, 同時還具有較為滿意的精度. 因此, 本文采用該模型進行海嘯數值模擬.

本文使用GeoClaw和30′精度海底地形數據對各有限斷層模型進行數值模擬, 期望通過使用合理的數值模型以及高分辨率的海底地形數據提供貼近實際的模擬條件, 使得各有限斷層模型的合理性論述結果更具有說服力.

1 海嘯數值模擬的構建

總體上講, 海嘯的數值模擬需要經過三個階段: 海嘯激發階段, 海嘯傳播階段和上岸階段. 海嘯激發階段即海底地震發生時, 引發海底斷層發生大規模變形, 使得海底表面產生一個突發的變形, 使海平面發生抬升. 與海嘯的越洋傳播相比, 震源斷層的變形在極短時間內完成, 因此在海嘯數值模擬中, 通常將海底變形與海平面水位的抬升作同步處理. 海面抬升水位在重力作用下向四周擴散, 并發展成重力長波. 由于海嘯傳播的大尺度特征, 模擬海嘯傳播過程中應考慮地球曲率等因素. 上岸階段即海嘯傳播至海岸淺水地段, 海嘯波波長急劇短縮, 海嘯波波高急劇增加, 形成巨大的“水墻”. 此階段由于水動力過程的非線性增強, 需要配合高分辨率地形數據進行計算以提升數值模擬的準確性, 而由于高分辨率地形數據的缺乏, 海嘯數值模擬的準確度受到了較大的限制. 本文使用30′精度的SRTM30地形數據(Shuttle Radar Topography Mission, http://www2.jpl.nasa.gov/srtm/)進行海嘯數值模擬, 期望得到一個較為合理的模擬結果.

1.1 海嘯源

2011年3月11日日本東北地震發生之后, 國際地震研究機構根據初步震源機制、斷層模型等參數確定了各自的震源有限模型[10]. 其中, 包括美國國家地質調查局提供的平均位錯模型、USGS瞬時破裂模型和USGS分時破裂模型, 中國地震局地球物理研究所提供的IGPCEA模型, 中國科學院地質與地球物理研究所提供的IGCAS模型以及日本國土地理院提供的GSI模型. 其具體參數見表1.

表1 各地區地震研究機構有限斷層模型的參數

圖1 IGCAS模型斷層分布圖

圖2 IGCAS模型海底變形圖

1.2 傳播模型

GeoClaw模型是Conservation Laws Package (ClawPack)軟件用于模擬地震海嘯的模塊. 其實現的技術路徑可分為兩步: 第一步, 利用Okada彈性半無限空間理論將有限斷層模型轉化為海底垂向變形, 進一步轉化為海面的初始水位; 第二步, 使用有限體積法求解二維淺水方程, 并考慮海底底部摩擦對于海嘯波傳播的影響, 同時引入了非線性限制器來抑制數值計算過程中產生的非物理振蕩, 使其在空間和時間都達到了二階精度, 并且采用自適應網格技術, 對海嘯波進行追蹤, 并在波高變化劇烈的區域進行網格加密, 來達到減少計算量及提高計算效率的目的.

模型的控制方程的守恒形式如下:

式中:為曼寧系數, 取0.025.

1.3 模型設置

本文海嘯波傳播模型計算區域為60°E～60°W, 80°S～80°N. 將計算區域劃分為以1°為間隔的粗網格, 共生成38400個網格. 并設置二級、三級精細化網格, 精度分別為30′和5′. 采用自適應網格技術, 一般模擬時使用1°精度網格; 當模擬的過程值達到閾值時, 切換至高分辨率網格, 繼續網格推進, 以達到優化計算效率的目的. 同時, 采用精度為30′的STRM30地形數據, 以提高模擬海嘯上岸階段的準確率.

模擬所用的海嘯觀測數據由NOAA(National Oceanic and Atmospheric Administration, https:// ngdc.noaa.gov/hazard/dart/2011honshu_dart.html)提供, 圖3為NOAA提供的2011年日本東北大地震海嘯傳播時間圖. 圖中表示了海嘯在不同時刻所在的位置, 圖中標有數字的方塊為該區域中記錄到海嘯波的浮標, 用于記錄浮標點處不同時刻的水位高度, 同時NOAA給出了各觀測點處的擬合潮汐分量與殘差值. 其中, 擬合潮汐分量為僅在潮汐影響下不同時刻的水位高度; 殘差值為觀測數據與擬合潮汐分量之差, 即除潮汐影響之外不同時刻的水位高度變化.

圖3 2011年日本東北大地震海嘯傳播時間圖

本文將使用殘差值與模擬海嘯波對比來論述模擬海嘯的合理性. 為保證模擬海嘯結果的合理性, 應使用多點浮標數據與模擬海嘯進行比對, 同時對各浮標殘差值進行分析篩選, 篩選原理為觀察海嘯波到達該浮標點之前是否存在異常值. 以46402點、46407點為例, 圖4為46402點與46407點的殘差時序圖, 46402點在海嘯波到達之前, 殘差在-0.05m附近波動, 當海嘯波到達該潮汐站時, 殘差為0.08538m. 此時可以認為當海嘯波傳到該浮標點時存在其他因素對海面高程產生影響, 此點與模擬海嘯數據差異較大的可能性應是其他因素導致實際觀測數據異常, 而不是模擬結果有問題. 46407點在海嘯波到達之前, 殘差數值在零附近波動, 當海嘯波到達該浮標, 殘差值產生明顯變化. 則可以認為在海嘯波傳至該點前并無其他因素對其附近海面產生影響, 所記錄殘差即由海嘯波產生. 此類浮標數據可與模擬數據進行對比.

使用上述篩選準則對上述浮標進行篩選, 并最終確定18個浮標, 同時在計算區域內設置對應的監測點. 浮標位置信息詳見表2.

表2 所選中的浮標列表

圖4 46402點(a)與46407點(b)的殘差時序圖

圖6 各觀測點記錄的海嘯波傳播時程曲線

2 海嘯數值模擬

2.1 IGCAS模型海嘯數值模擬

本文共使用6種地震有限斷層模型來進行海嘯數值模擬, 本節以IGCAS模型進行闡述.

IGCAS斷層模型在海底破裂引發了海平面的抬升, 抬升水位受到重力影響向四周擴散形成海嘯波. 同時當模擬海嘯傳播時, 設置在浮標坐標處的監測點記錄海面水位的變化. 圖5、圖6分別為IGCAS模型激發的海嘯在2h內的傳播過程以及各監測點所記錄的海面水位時序變化情況.

圖6中, 顏色較淺曲線表示各浮標點記錄的殘差值隨時間的變化, 顏色較深曲線則為模擬海嘯傳播時各監測點所記錄的海面隨時間的變化, 橫坐標以海嘯激發的時刻為0, 縱坐標以海平面為基準海面水位. 通過觀察各個浮標點觀測值與模擬值的對比并對其進行分析: (1)21413點與21419點的觀測值在0時附近出現異常值. 兩點距離震源有一定距離, 而異常值出現在海嘯激發時, 因此可基本排除是由海嘯波傳播引起的海平面變化; (2)各監測點所記錄的模擬海嘯波的波形基本上都與浮標點所記錄的殘差數據吻合. 同時對于接近震源的監測點(21413點、21414點、21415點等)模擬海嘯波的到時晚于實際海嘯波, 而對于距離震源較遠監測點, 模擬海嘯波的到時早于實際海嘯波. 此現象在32413點與51406點表現最為明顯. 因此可以認為使用IGCAS模型在進行模擬海嘯時, 海嘯波波速略大于實際海嘯波的波速. (3)使用IGCAS模型進行海嘯數值模擬時, 模擬海嘯可以較好地匹配實際海嘯的波高與到時. 但是當海嘯波穿過該浮標點之后, 實際觀測數據中依舊存在較大幅度的海面波動, 而模擬海嘯的海嘯波在穿過監測點之后, 海面迅速趨于水平狀態. 造成這種情況的原因主要是由于在設置邊界條件時, 采用了開放式邊界, 因此無法考慮海嘯波傳播至固體界面時產生的反射波對海面的影響.

表3 各監測點模擬波高與實測波高差值

注: (1)絕對差(m)=觀測波高-模擬海嘯波高; (2)相對差(%)=(觀測波高-模擬海嘯波高)/實測殘差.

2.2 各有限斷層模型誤差分析

分別將各個震源模型模擬得到的監測點海面波動時序圖與對應浮標數據進行對比, 均可在浮標數據時序圖中匹配到相應形狀. 模擬結果和實測數據的接近程度可以在一定程度上反映模擬程序的可靠性, 形狀的匹配程度一般由海嘯傳播至觀測點時的波高與到時所控制, 而通過觀察各觀測點的模擬結果與浮標數據的對比, 并沒有發現模擬海嘯波到時嚴重偏離實際海嘯波到時的現象. 因此, 研究各斷層模型海嘯模擬結果可靠性的問題便簡化成了通過對比各觀測點模擬海嘯波高的誤差來分析各斷層海嘯模擬結果的可靠性問題. 表3即為各監測點上觀測波高與模擬波高之間的絕對差以及它們的相對誤差.

圖7 各觀測點模擬波高與實測波高絕對差

圖8 各觀測點模擬波高與實測波高相對差

為更直觀地了解數據變化趨勢, 將上述觀測點根據海嘯波到時進行排序并在圖中表示, 圖7為各觀測點模擬波高與實測波高絕對差圖, 圖8為各觀測點模擬波高與實測波相對差圖

通過對上述的數據與圖形進行分析評價各有限斷層的合理性, 最終得出如下結果:

(1)從圖7可以發現, 各觀測點的絕對誤差總體還是隨著傳播時間的增加變小. 當海嘯波向四周海域傳播時, 海嘯的能量也開始進行消耗, 最直觀的表現即海嘯波波高的降低. 若觀測點處無其他干擾因素, 一般情況下傳播時間長的觀測點模擬波高與觀測波高的差值普遍較小. 由于較遠觀測點殘差較低, 在進行相對誤差分析時, 若存在其他因素干擾, 則極易在相對誤差中放大. 因此, 就出現了圖8的情況, 即在絕對誤差較小的點處, 卻有著較大的相對誤差.

(2)平均位錯模型與GSI模型相較于其他模型, 在準確性上存在較大差距. 通過簡單分析可知, 平均位錯模型得到的震級雖然與現實接近, 但是該模型只存在一個斷層, 由斷層破裂引起的海底變形相對平均, 即地震能量并不集中, 因此模擬得到的海嘯波波高普遍低于實際觀測值. 對于GSI模型而言, 本身計算得到的地震矩小于實際情況, 因此不難發現該模型的模擬波高均小于觀測波高. 雖然, 兩個模型的準確性有待提升, 但是不可否認這兩個模型在海嘯波的到時以及各監測點的海平面變化趨勢與實際觀測的海平面變化具有一定的吻合. 而使用這些模型(尤其是平均位錯模型)進行海嘯預警時, 在預測海嘯波傳播位置的確定性方面具有一定的參考價值.

(3)IGCAS模型的模擬結果與實際觀測數據較為接近, 其18個監測點中僅存在2個點數值異常. 在32413點, 所有斷層模型均具有較大的誤差, 并且該點距離震源距離是所有18個點中最遠的, 本身距離越遠, 模擬結果的準確率就容易下降, 但同時也不能排除其他干擾因素導致該點海面變化異常的可能性. 在51407點, 實際波高與IGCAS模型的模擬波高差為12.2cm, 相對誤差為46.47%. 而IGPCEA模型、USGS瞬時/分時模型在該點的相對誤差分別為26.07%、－4.38%、－9.12%. 這3個模型在該點模擬結果表現良好, 可排除其他因素對該浮標點檢測數據造成影響. 因此, 可認為IGCAS模型在51407點處的模擬結果確實存在較大誤差. 所以, 當使用該模型進行海嘯數值模擬時, 應警惕可能出現的異常觀測點.

(4)通過USGS瞬時破裂模型與分時破裂模型的模擬結果的對比, 我們可以初步得到以下結論: 對于范圍大、持續時間久的海嘯傳播過程, 分別使用相同斷層參數的瞬時破裂模型與分時破裂模型進行海嘯數值模擬, 模擬的海嘯波在遠場海域基本無影響, 在距離震源較近的觀測點(21413點)處雖存在一定差異, 而相較于本身的波高來說, 差異也在可接受范圍內. 此結論在本次日本東北地震海嘯模擬中成立, 但在更大范圍內的適用性還需要后續更多驗證.

(5)IGPCEA模型的模擬結果與實際數據接近. 雖然在18個監測點中, 51406點與32413點模擬數據與實際觀測數據對比有一定差距, 考慮到這兩個觀測點與震源的距離較遠, 結合其他模型均在該點處有較大相對誤差, 因此我們可以初步認為在這兩點出現了其他的因素影響海面水位的變化. 因此, 初步可以認為使用IGPCEA模型進行的海嘯數值模擬時, 其模擬結果具備較好的準確性.

3 總結

本文對現有日本東北大地震斷層模型進行海嘯數值模擬, 以進一步論述各斷層模型的合理性. 此方法的討論雖早已存在, 但因為受到當時海嘯數值模型以及地形數據精度等限制, 模擬得到的結果不盡人意. 現今, 筆者使用GeoClaw海嘯數值模型與更高精度的地形數據進行海嘯數值模擬, 期望獲得更為精確的結果. 通過各斷層模擬結果之間的橫向對比, 最終, IGPCEA模型的模擬結果與實際觀測結果更為接近, 表明此斷層模型相較于其他斷層模型較為合理.

本次的研究依舊存在不少需要加強之處. 在海嘯的數值模型上, 對數值模型的邊界設置僅采用了開放式邊界, 而忽略了海嘯波傳播至固體邊界產生的反射波對海面的影響. 另外, 在進行此類越洋海嘯的數值模擬時, 忽略了科里奧利力對海嘯傳播的影響. 在觀測數據選擇上, 直接使用了殘差值作為海嘯觀測數據, 而忽略了其他因素(如風暴潮、颶風等)對于殘差值的影響. 期望在未來的研究中, 通過對海嘯數值模型的優化與升級得到更加具有說服力的模擬結果.

[1] JMA. The 2011 off the Pacific coast of Tohoku Earthquake 28th report[R]. 2011.

[2] USGS. Largest earthquakes in the world since 1900[R]. 2011.

[3] Hayes G P. Rapid source characterization of the 2011M9.0 off the Pacific coast of Tohoku Earthquake[J]. Earth, Planets and Space, 2011, 63(7):529-534.

[4] Hayes G P, Earle P S, Benz H M, et al. 88 hours: The U.S. geological survey national earthquake information center response to the 11 March 2011M9.0 Tohoku Earthquake[J]. Seismoligical Research Letters, 2011, 82(4):481-493.

[5] 張勇, 許力生, 陳運泰. 2011年3月11日日本本州東海岸近海地震破裂過程快速反演結果v1[EB/OL]. [2021- 05-16]. 中國地震局地球物理研究所. http://ddsep.cea- igpac.cn/rp/2011/20115e74367081165e565e5672c4e1c6d775cb88fd16d7757309707783488fc77a0b5fe901f53cd6f147ed3679c.

[6] Wang W, Hao J, Yao Z. Preliminary result for rupture process of Mar.11, 2011,M8.9 earthquake, near the East Coast of Honshu, Japan [EB/OL]. [2021-05-22]. http:// www.igg.cas.cn/xwzx/zhxw/201103/t20110312_3082869.html.

[7] Imakiire T, Kobayashi T. Crustal deformation and fault model of the 2011 off the Pacific coast of Tohoku Earthquake[J]. Bulletin of the Geospatial Information Authority of Japan, 2011, 2011:70-79.

[8] 王振, 孟國杰, 橫田佑助, 等. 利用1-Hz GPS波形數據反演2011年日本東北大地震震源破裂過程[J]. 地震, 2013, 33(3):13-23.

[9] 周正陽, 張勇. 利用海嘯波數據反演特大海嘯地震的地震矩張量[C]//2015中國地球科學聯合學術年會論文集(十三)——專題37地球氣候系統歷史、專題38強震機理、孕育環境與地震活動性分析、專題39大數據時代地球物理信息學及其應用. 北京, 2015:2.

[10] 景惠敏, 張懷, 吳忠良, 等. 利用海嘯數值模擬結果進行海底地震有限斷層模型驗證[J]. 地震, 2013, 33(4): 207-213.

[11] 藝帆, 安超. 海嘯反演對于斷層位置不確定性的敏感性分析[C]//2020年中國地球科學聯合學術年會論文集(十五)——專題四十三: 海洋地球物理、專題四十四: 海嘯及海嘯預警研究、專題四十五: 電磁地球物理學研究應用及其新進展. 重慶, 2020:27.

[12] Titov V V, Gonzalez F I. Implementation and testing of the method of splitting tsunami (MOST) model. [EB/OL]. [2021-05-22].https://www.pmel.noaa.gov/pubs/PDF/tito 1927/tito1927.pdf.

[13] Imamura F, Shuto N, Goto C. Numerical simulations of the transoceanic propagation of tsunamis[C]//Proceedings of the Sixth Congress Asian and Pacific Regional Division, IAHR. Kyoto, Japan. 1988:265-272.

[14] 潘文亮, 王盛安. COMCOT數值模式的介紹和應用[J]. 海洋預報, 2009, 26(3):45-52.

[15] George D L, LeVeque R J. Finite volume methods and adaptive refinement for global tsunami propagation and local inundation[J]. Science of Tsunami Hazards, 2006, 24(5):319.

[16] Okada Y. Surface deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 1985, 75(4):1135-1154.

[17] Hanks T C, Kanamori H. A moment magnitude scale[J]. Journal of Geophysical Research: Solid Earth, 1979, 84 (B5):2348-2350.

（責任編輯 章踐立）

Verification of finite fault model of the 2011 Tohoku Tsunami based on GeoClaw tsunami simulation model

HU Pei1, JING Huimin1*, YANG Wei2

( 1.Faculty of Mechanical Engineering & Mechanics, Ningbo University, Ningbo 315211, China;2.College of Science & Technology Ningbo University, Ningbo, 315300, China )

Tsunami waves are characterized by long periods and slow subsidence. For the submarine earthquake triggering a tsunami, the tsunami waveform data can be used as supplementary data for verifying the rationality of the finite fault model obtained by inversion of seismic wave data and geodetic data. Taking the earthquake source of the 2011 Tohoku tsunami as an example, we transform the finite fault model into vertical variation in the seabed topography according to Okada elastic half-space theory. Then we use the GeoClaw tsunami model to perform the tsunami simulation. Some monitoring points in the calculation domain are set on those points which can be utilized to obtain the tsunami waveform data. By comparing the simulated tsunami waveform data with the actual observation data, we analyze the reliability of the models and verify the effectiveness of most reasonable finite fault model of the earthquake.

Tohoku earthquake; finite fault model; GeoClaw numerical simulation

P315

A

1001-5132（2022）03-0106-09

2021?10?08.

寧波大學學報（理工版）網址: http://journallg.nbu.edu.cn/

國家自然科學基金青年科學基金（41304072）.

胡培（1993－）, 男, 浙江寧波人, 在讀碩士研究生, 主要研究方向: 海嘯數值模擬. E-mail: hup93@outlook.com

通信作者:景惠敏（1984－）, 女, 山東德州人, 博士, 主要研究方向: 地震海嘯的數值模擬. E-mail: jinghuimin@nbu.edu.cn