國立中央大學地震災害鏈風險評估及管理研究中心
大規模地震情境模擬對外資料

2021-東部地震情境模擬

        東部衝擊情境，預期成果係完成琉球隱沒帶情境震源參數擬定，並建立指標性特徵震源模型，依據擬定震源模型進行地震波模擬分析，最終提供災損分析使用之峰值地動加速度(peak ground acceleration, PGA)、峰值地動速度(peak ground velocity, PGV)、峰值地動位移(peak ground displacement, PGD)、0.3秒譜加速度(0.3-s of spectral acceleration, SA0.3)與 1.0 秒譜加速度(1.0-s of spectral acceleration, SA1.0)，其地動數值涵蓋空間範圍包含整個臺灣東部主要地區(宜蘭、花蓮及台東北段)。
        工作項目包含震源參數擬定、速度與構造模型建置、特徵震源模型建置、地震波模擬分析、不同數值模擬方法比較及地表地震動歷時與反應譜計算分析等工作。其中，參數與方法部分，模擬分析所需參數蒐集及擬定部分，陳述如何針對琉球隱沒帶地震蒐集相關文獻用以擬定震源模型參數；地形與速度構造模型建置部分，說明供後續情境地動模擬使用且覆蓋分析區域之地震波速度構造模型；特徵震源模型建置部分，考慮琉球隱沒帶地震之特徵規模及臺灣東部地區可能之地動影響模式，設定可能震源破裂模式，詳加擬定震源參數，並製表呈現；地震波模擬分析部份闡述使用方法；不同低頻震波數值模擬方法之比較。另，地震動歷時與反應譜計算分析部分，則針對地表地震動歷時與反應譜計算分析，提出涵蓋研究區域內之模擬成果。

PGA (gal) PGV (cm/s) PGD (cm) SA0.3 (g) SA1.0 (g)

參數與方法

此處闡述考量琉球隱沒帶情境地震而彙整之震源參數，並以此作為依據，以地震學、地體動力學與工程地震學角度，考慮對臺灣東部地區較具衝擊代表性之情境，擬定特徵震源模型。另對於速度、地形構造之數值模擬所需之模型建構進行考量。最終針對該單一特徵震源模型進行情境模擬地震動，模擬完整情境地震歷時，計算地震動參數，包含(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)、(4) 0.3 秒譜加速度(SA0.3)與(5) 1.0 秒譜加速度(SA1.0)等五項參數作為產出，以利後續災損推估小組推估臺灣東部地區震災情境。綜整上述說明，模擬分析程序詳見下圖 。

圖 2- 1 震源情境地震動模擬分析流程
圖 ● 震源情境地震動模擬分析流程

震源參數擬定

        臺灣東北部外海琉球海溝具引致大規模隱沒帶地震之潛勢，歷史上於 1920年，此區域卽發生規模達 8.0 之地震事件，Theunissen et al. (2010)提出該區域潛勢斷層 ISZ:板塊介面孕震區域(interpolate seismogenic zone)模式為最可能發生之1920 年地震與琉球隱沒帶地震之破裂型態。此外，Hsu et al. (2012)使用 2005 年至 2010 年間之連續全球定位系统(Global Positioning System, GPS)資料，推估琉球隱沒帶斷層幾何及滑移虧損(slip-deficit)，該區域具有發生震矩規模 7.5 至 8.7 的災害地震潛勢。本研究綜合考量上述論文研究成果，並根據國家地震工程中心SHACC level 3 計畫綜整近年臺灣海陸域斷層調查成果(PSHA SSHAC-3 SSC, 2019)，作為擬定琉球隱沒帶情境震源參數之依據。

圖 ● SSHAC Level 3 琉球隱沒帶震源參數模型

地形與速度構造模型建置

        地震波傳模擬所需速度構造模型建置，係以臺灣地區速度構造研究(KuoChen et al., 2012)之三維速度分布數值，作為本研究模擬計算之模型建置基礎資料。該速度分布數據係藉由中央氣象局強震站及短週期地震站、中央氣象局及中央研究院寬頻地震站、大量布設臨時站所記錄到之大量地震波走時資訊，經由逆推 方法計算建構而成 。 參考三維速度分布數值所強調之構造不均質(heterogeneity)與側向變化(lateral variation)，提升了僅考量無方向性變化的一維速度模型可能影響地震模擬評估之不足。評估本先導研究情境模擬工作所需要之涵蓋空間範圍後，整合三維速度分布數值所建構之臺灣本島及外海區域速度構造模型。

        此外，因數值方法特性，本研究所採行之地震動模擬方法於低頻模擬時可考慮地形起伏造成的震波散射效應，例如震波傳遞至山峰及山脊處時所造成振幅放大，山谷處造成振幅衰減(Lee et al., 2009)。因地形效應額外產生之表面波與散射波，預期會影響測站處記錄之波形。本研究計畫以解析度為一弧秒(約 1.85 km)之ETOPO1 (Amante and Eakins, 2009)做為初始地形模型，考量前述臺灣速度構造模型後，經穩定性測試後，計算出震波模擬適合之模型解析，將地形模型重採樣為 300 m 取一格點，作為數值計算網格水平方向格距，建立數值網格模型，模型空間範圍考量納入涵蓋琉球隱沒帶情境震源與臺灣東部目標場域之範圍。

圖 2- 3 臺灣本島及外海區域地下速度構造模型
圖 ● 臺灣本島及外海區域地下速度構造模型
圖 2- 4 臺灣本島及外海區域數值地形模型
圖 ● 臺灣本島及外海區域數值地形模型
圖 2- 5 臺灣本島及外海區域數值地形與速度構造整合模型
圖 ● 臺灣本島及外海區域數值地形與速度構造整合模型

特徵震源模型建置

        考量琉球隱沒帶具轉折之線形與延伸破裂面，本完整詳實地震動模擬須簡化其琉球隱沒帶線形，合理考量符合線型之斷層長度與寬度，進而假定特徵震源模型。

        特徵震源模型之參數擬定原則，則採用日本地震調查研究推進本部之 Recipe(Irikura and Miyake, 2011)方法，取其特徵震源擬定之標準化流程作為依據，針對琉球隱沒帶情境地震破裂，探討長週期地震動對臺灣東部區域之影響性，假定斷層面上具單一地栓(asperity)為主要能量釋放位置，將地栓中心置於沿斷層走向之中線，地栓上緣與斷層上緣重合，並考慮震源(即破裂起始點)位於地栓上緣最外海處，考量斷層由海床往深處破裂琉球隱沒帶情境地震震源模型以平面與側向示意；模型之震源位置(即破裂起始點),斷層面之巨觀與微觀震源參數分列於下表。

圖 2- 6 琉球隱沒帶情境震源之特徵震源模型平面示意圖
圖 ● 琉球隱沒帶情境震源之特徵震源模型平面示意圖
圖 2- 7 琉球隱沒帶情境震源之特徵震源模型之側向示意圖
圖 ● 琉球隱沒帶情境震源之特徵震源模型之側向示意圖
表 2- 1 琉球隱沒帶情震模型震源位置參數
圖 ● 琉球隱沒帶情震模型震源位置參數
表 2- 2 琉球隱沒帶情震地震巨觀震源參數表
圖 ● 琉球隱沒帶情震地震巨觀震源參數表
表 2- 3 琉球隱沒帶情震地震微觀震源參數表
圖 ● 琉球隱沒帶情震地震微觀震源參數表

虛擬測站點位

點位分布設定

        本研究配合災損推估小組之評估網格設定，依據國家災害防救科技中心(National Science and Technology Center for Disaster Reduction, NCDR)格點解析度 500 公尺，確認分布範圍為東經 118.125 度至 122.005 度，北緯 21.895 度至 26.385 度，範圍涵蓋臺灣本島與離島範圍，共 132,712 網格點。下圖為 NCDR 網格示意圖，黑色格點即為網格點位。

圖 2- 8 國家災害防救科技中心(NCDR)500 公尺網格點位
圖 ● 國家災害防救科技中心(NCDR)500 公尺網格點位
本研究設定點位

        對應災損推估小組後續評估工作，空間範圍與解析上，本研究基於 NCDR 500 公尺網格點擬定虛擬測站，進行地震動模擬計算。所採行點位分布範圍為東經120.500 度至 122.005 度，北緯 22.1 度至 25.5 度，共 24,755 點位，空間點位選取範圍系考慮琉球隱沒帶情境震源所在區域、臺灣東部生活圈，及關鍵基礎設施進而擬定之範圍。分布於臺灣東部陸地區域之灰色格點為 NCDR 500 公尺網格點位，視為虛擬測站；藍框區域為情境地震動模擬範圍。

圖 2- 9 本研究設定點位與震波模擬範圍對照圖
圖 ● 本研究設定點位與震波模擬範圍對照圖

地震動歷時與反應譜計算分析

依前述地震波模擬分析方法，考慮琉球隱沒帶情境地震之特徵震源模型，分別進行低頻與高頻模擬，進而以混合法合成於工程基盤之寬頻加速度地震歷時，爾後分別積分一次得速度歷時、積分二次得位移歷時，各地震歷時之最大振幅處即分別為(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)等地動數值。另一方面，選取加速度地震歷時作為輸入地震歷時，考慮 5%阻尼計算該地震歷時之反應譜(response spectrum)，爾後分別取週期 0.3 秒與週期 1.0 秒之譜加速度作為反應譜地動數值。本研究所模擬之地震動數值敘述如下：

考慮震源位於地栓上緣最外海處，斷層由海床往深處破裂，因斷層破裂的方向指向臺灣本島花蓮方向，破裂的方向性效應強烈，此起情境引致臺灣花蓮區域地震動甚鉅，PGA 分布顯示最大值為 594.4 cm/s2，出現在花蓮縣豐濱鄉一帶；PGV 分布，最大值為 141.5 cm/s，僅在非常少數區域出現此極端值，而豐濱鄉、瑞穗鄉、光復鄉及鳳林鎮區域皆有超過 80 cm/s 的 PGV 分布；該區域附近具有最大對應之震度為 6 強；PGD 分布在秀林鄉、萬榮鄉南部及豐濱鄉沿海一帶有超過 75cm 的分布，而最大值達到 98cm；SA0.3 分布與 SA1.0 分布之最大值則分別為 1.01 g 與 0.73 g。

圖 3-1 琉球隱沒帶情境地震模擬臺灣東部地區 PGA 分布(圖 a)與PGV 分布(圖 b)
圖 ● 琉球隱沒帶情境地震模擬臺灣東部地區 PGA 分布(圖 a)與 PGV 分布(圖 b)
圖 3-2 琉球隱沒帶情境地震模擬臺灣東部地區震度分布(圖 a)與 PGD 分布(圖 b)
圖 ● 琉球隱沒帶情境地震模擬臺灣東部地區震度分布(圖 a)與 PGD 分布(圖 b)
圖 3-3 琉球隱沒帶情境地震模擬臺灣東部地區之 SA0.3 分布(圖 a)與 SA1.0 分布(圖 b)
圖 ● 琉球隱沒帶情境地震模擬臺灣東部地區之 SA0.3 分布(圖 a)與 SA1.0 分布(圖 b)

成果總結

        過去，琉球隱沒帶地震事件曾引致地震災害，如 1920 年琉球隱沒帶地震事件、1986 花蓮外海地震事件等，皆造成傷亡與建物毀損。本研究計畫擇定琉球隱沒帶作為假想目標，考慮單一地栓位置與震源向西破裂起始位置、震波傳遞等因素，擬定單一震源破裂情境，針對臺灣東部地區以混合法進行地震動模擬，完成情境模擬地動參數，包含峰值地動加速度、峰值地動速度、峰值地動位移、0.3 與1.0 秒譜加速度等，供後續災損推估小組進行評估。研究成果顯示，該震源模式對臺灣東部區域引致高地動量值，尤其是花蓮縣豐濱鄉、瑞穗鄉、光復鄉及鳳林鎮附近區域，具有 6 強以上震度；且引起之地動值兼具短周期及長周期之地震動，可能引致較大規模之地震災害。

        然而，本研究僅考慮單一地栓震源破裂情境，倘若琉球隱沒帶破裂，不全然為本研究所假想之震源破裂型態以及推估規模。故本研究建議如下：

  1. 針對單一地栓震源破裂面，未來可考慮不同地栓位置之滑移量分布模式、震源破裂位置等多種破裂情境，對各情境產出之地震動參數進行統計分析，有效量化地震動中值與極值，估算地震動範圍。
  2. 對於分析地點進行場址放大因子研究，提供更高空間解析之準確場址放大，藉以提升地動空間解析。
  3. 有限差分法或譜元素法皆為低頻三維地動模擬有效工具，透過量化比較指出，兩者計算結果相似；但須注意操作過程使用之參數及流程之正確性。

參考文獻

Amante, C. and Eakins, B. W. (2009). ETOPO1 1 arc-minute global relief model:
         procedures, data sources and analysis, NOAA Technical Memorandum, NESDIS NGDC-24, 19 pp.
Atkinson, G. M and Assatourians, K. (2015). Implementation and validation of EXSIM
         (a stochastic finite-fault ground-motion simulation algorithm) on the SCEC broadband platform, Seismol. Res.. Lett., 86(1), 48–60, DOI:10.1785/0220140097.
Atkinson, G. M. and Macias, M. (2009). Predicted ground motions for great interface
         earthquakes in the Cascadia subduction zone, Bull. Seism. Soc. Am., 99(3), 1552–1578.
Berenger, J. A., (1994). Perfectly matched layer for the absorption of electromagnetic
         waves, J. Comput. Geophys., 114, 185–200.
Chaljub, E. (2009). Spectral element modeling of 3D wave propagation in the Alpine
          valley of Grenoble, France, in ESG 2006, Third Intl. Symposium on the Effects of Surface Geology on Seismic Motion, P.-Y. Bard, E. Chaljub, C. Cornou, F. Cotton, and P. Guéguen (Editors), LCPC Editions, ISSN 1628-4704, Vol. 2, 1467–1473.
Chaljub, E., Moczo, P., Tsuno, S., Bard, P. Y., Kristek, J., Käser, M., … & Kristekova, M.
          (2010). Quantitative comparison of four numerical predictions of 3D ground motion in the Grenoble Valley, France. Bulletin of the Seismological Society of America, 100(4), 1427-1455.
Gottlieb, D. and Turkel, E. (1976). Dissipative 2-4 methods for time-dependent
           problems, Math. Comput., 30, 703–723.
Gropp, W., Lusk, E., Doss, N., and Skjellum, A. (1996). A high-performance, portable
           implementation of the MPI message passing interface standard, Parallel Comput., 22(6), 789–828.
Hixon, R. (1997). On increasing the accuracy of MacCormack schemes for aero acoustic
           applications, AIAA Paper, doi: 10.2514/6.1997–1586.
Hsu, Y. J., Ando, M., Yu, S. B., and Simons, M. (2012). The potential for a great
            earthquake along the southernmost Ryukyu subduction zone, Geophys. Res. Lett.,39, L14302, doi:10.1029/2012GL052764.
Huang J. Y., Wen, K. L., Lin, C. M., Kuo, C. H., Chen, C. T., and Chang, S. C. (2017). Site
            correction of a high-frequency strong-ground-motion simulation based on an empirical transfer function, J. Asian Earth Sci., 138, 399–415.
Irikura, K. and Miyake, H. (2011). Recipe for predicting strong ground motion from
             crustal earthquake scenarios, Pure Appl. Geophys., 168, 85–104,doi:10.1007/s00024-010-0150-9.
Käser, M., M. Dumbser, and J. de la Puente Alvarez (2006). An efficient ADER-DG
             method for 3-dimensional seismic wave propagation in media with complex geometry, in ESG 2006, Third Intl. Symposium on the Effects of Surface Geology on Seismic Motion, P.-Y. Bard, E. Chaljub, C. Cornou, F. Cotton, and P. Guéguen(Editors), LCPC EditionsVol. 1, pp. 455–464, Grenoble.
Komatitsch, D., & Vilotte, J. P. (1998). The spectral element method: an efficient tool
             to simulate the seismic response of 2D and 3D geological structures. Bulletin of the seismological society of America, 88(2), 368-392.
Komatitsch, D., & Tromp, J. (1999). Introduction to the spectral element method for
             three-dimensional seismic wave propagation. Geophysical journal international,
139(3), 806-822.
Kristeková, M., Kristek, J., & Moczo, P. (2009). Time-frequency misfit and goodness-of fit criteria for    quantitative comparison of time signals. Geophysical Journal International, 178(2), 813-825.
Kristek, J., P. Moczo, and P. Pazak (2009). Numerical modeling of earthquake motion in
            Grenoble basin, France, using a 4th-order velocity- stress arbitrary discontinuous staggered-grid FD scheme, in ESG 2006, Third Intl. Symposium on the Effects of Surface Geology on Seismic Motion, P.-Y. Bard, E. Chaljub, C. Cornou, F. Cotton,
and P. Guéguen (Editors), LCPC Editions, ISSN 1628-4704, Vol. 2, p. 1517–1526.
Kuo-Chen, H., Wu, F. T., and Roecher, S. W. (2012). Three-dimensional P velocity
            structures of the lithosphere beneath Taiwan from the analysis of Taiger and related seismic data sets, J. Geophys. Res., B6.
Lee, S. J., Komatitsch, D., Huang, B. S., and Tromp, J. (2009). Effects of topography on
             seismic-wave propagation: an example from northern Taiwan, Bull. Seism. Soc. Am., 99, 314–325.
Lee, Y. T., Ma, K. F., Wang, Y. J., and Wen, K. L. (2015). An empirical equation of effective
             shaking duration for moderate to large earthquakes, Nat. Hazards, 75(2), 1779–1793.
MacCormack, R. W., (1969). The effect of viscosity in hypervelocity impact cratering,
              AIAA Paper, 69–354.
Maeda, T., Takemura, S., & Furumura, T. (2017). OpenSWPC: an open-source integrated
               parallel simulation code for modeling seismic wave propagation in 3D heterogeneous viscoelastic media. Earth, Planets and Space, 69(1), 1-20.
Marcinkovich, C. and Olsen, K. (2003). On the implementation of perfectly matched
                layers in a three-dimensional fourth-order velocity-stress finite difference scheme, J. Geophy. Res., 108, doi: 10.1029/2002JB002235. issn: 0148-0227.
PSHA SSHAC-3 SSC, 2019, Development of the Hazard Input Document for Taiwan
               using SSHAC Level 3 Methodology. Volume 2: SSC Technical Report, Taiwan Power Company.
Sokolov, V., Loh, C. H., and Jean, W. Y. (2006). Strong ground motion source scaling and
               attenuation models for earthquakes located in different source zones in Taiwan,4th International Conference on Earthquake Engineering, , Paper No. 003.
Sokolov, V., Wen, K. L., Miksat, J., Wenzel, F., and Chen, C. T. (2009). Analysis of Taipei
                basin response for earthquakes of various depths and locations using empirical data, Terr. Atmos. Ocean. Sci., 20, 687–702, doi: 10.3319/TAO.2008.10.15.01(T).
Stupazzini, M. (2009). 3D Ground Motion Simulation of the Grenoble Valley by
                GeoELSE, in ESG 2006, Third Intl. Symposium on the Effects of Surface Geology on Seismic Motion, P.-Y. Bard, E. Chaljub, C. Cornou, F. Cotton, and P. Guéguen(Editors), LCPC Editions, ISSN 1628-4704, Vol. 2, 1551–1560.
Stupazzini, M., R. Paolucci, and H. Igel (2009). Near-fault earthquake ground motion
                simulation in the Grenoble Valley by a highperformance spectral element code,Bull. Seismol. Soc. Am. 99, 286–301.
Tam, C. and Webb, J. C. (1993). Dispersion-Relation-Preserving finite difference
               schemes for computational acoustics, J. Comput. Phys., 107, 262–281.
Theunissen, T., Font, Y., Lallemand, S., and Liang, W. T. (2010). The largest
               instrumentally recorded earthquake in Taiwan: Revised location and magnitude, and tectonic significance of the 1920 event, Geophys. J. Int., 183, 1119–1133, doi:10.1111/j.1365-246X.2010.04813.x.
Zhang, W., Zhang, Z. G., and Chen, X. F. (2012). Three-dimensional elastic wave
               numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids, Geophys. J. Int., 190, 358–378.
行政院 (2016) 民國 105 年災害防救白皮書。
中央研究院 (2015) 大規模地震災害防治策略建議書。中央研究院報告 No.13