國立中央大學地震災害鏈風險評估及管理研究中心
大規模地震情境模擬對外資料
2019-山腳斷層特徵震源地動模擬
鑒於先導型計畫逐步落實分析陸域山腳斷層與海域琉球隱沒帶對北台灣都會區之影響性，民國 108 年起，透過「科技部災害防救科技創新服務方案」進行「山腳斷層多重震源情境模擬」研究計畫乙案，為精進前期研究成果，納入更精細之地表地形模型與台北盆地幾何，並擇定可能影響台北盆地與鄰近區域之山腳斷層情境地震事件，整合出具代表性之震源巨觀參數，並分析評估特定震源特徵模型之於盆地區域之地震動分佈，作為提供後續災損分析評估之科學性基礎輸入資訊。本研究「山腳斷層特徵震源地動模擬」作為輔助執行「山腳斷層多重震源情境模擬」研究計畫，採用震源情境地震動模擬分析流程進行震源參數擬定與地震動模擬計算，預期成果係完成山腳斷層情境震源參數擬定，並建立指標性特徵震源模型，依據擬定震源模型考量台北盆地與鄰近區域地下構造進行地震波模擬分析，最終提供可作為災損分析輸入參數之峰值地動加速度(peak groundacceleration, 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)，其地震動數值涵蓋空間範圍包含台北盆地與鄰近區域。
參數與方法
此處闡述考量山腳斷層破裂之情境地震而彙整之震源參數，並以此作為依據，以地震學、地體動力學與工程地震學角度，考慮對台北盆地與鄰近區域較具衝擊之情境，擬定特徵震源模型。另對於速度、地形構造之建構進行考量。最終針對該單一特徵震源模型進行情境模擬地震動，模擬完整情境地震歷時，計算地震動參數，包含(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)、(4) 0.3 秒譜加速度(SA0.3)與(5) 1.0 秒譜加速度(SA1.0)等五項參數作為產出，以利後續災損推估小組推估台北盆地區域震災情境。綜整上述說明，模擬分析程序詳見下圖。
震源參數擬定
本研究基於台灣地震模型(Taiwan Earthquake Model, TEM)組織於 2015 發表之全台 38 條孕震構造，以及中央地質調查所 2012 年公布之台灣活動斷層分佈等資訊，以北台灣之重要控制斷層－山腳斷層作為大規模地震情境模擬之假想對象。山腳斷層陸域部分，TEM 版本與地調所版本資訊無較大歧異，線型分佈類似，山腳斷層為呈北北東走向之正斷層；此外，TEM 組織以高解析度地形資料審視構造線位置，對山腳斷層線型轉折與連續性有較精細之描繪。本研究旨於考量山腳斷層引發地震之情境對台北盆地區域之影響性研究，借鏡 106 年執行「大規模地震模擬情境」同以山腳斷層破裂情境為例，參考並彙整 TEM 版本與地調所版本山腳斷層資訊與線型作為情境模擬依據。
地形與速度構造模型建置
地震波傳模擬所需速度構造模型建置，係以近期台灣地區速度構造研究(Kuo-Chen et al., 2012)之大尺度三維速度分佈數值，作為本研究模擬計算之模型建置基礎資料。該速度分佈數據係藉由中央氣象局強震站及短週期地震站、中央氣象局及中央研究院寬頻地震站、大量布設臨時站所記錄到之大量地震波走時資訊，經由逆推方法計算建構而成。此外，針對台北盆地區，另行考量盆地幾何形狀與震波傳遞速度分佈，本研究援引 Lee et al. (2008)所建構之台北盆地三維速度分佈數值，P 波與 S 波速度的深度剖面分別如下圖所示。本研究為通盤考慮台北盆地沉積層、盆地效應、深部構造對地震動之影響性，結合前述台灣地區大尺度三維速度分佈數值與台北盆地三維速度分佈數值，精進數值模擬對台北盆地的解析度。由於大尺度三維速度分佈數值僅有 P 波速度資料，因此大尺度 S 波三維速度分佈數值根據彈性波傳遞理論，假設地下速度構造特性為泊松固體(Poisson Solid)，由 P 波與 S 波速比值為√3 換算而得；盆地區域則考慮後續地震波傳遞模擬計算至工程基盤，設定最小 S 波波速為 760 m/s。結合不同尺度三維速度分佈後，根據彈性波傳遞理論，檢核空間各位置 P 波與 S 波速比值須滿足大於或等於√2，使波動方程式成立方可進行數值計算求解。
因數值方法特性，本研究採行之地震動模擬方法於低頻模擬時可考慮地形起伏造成的震波散射效應，例如震波傳遞至山峰及山脊處時所造成震幅放大，山谷處造成震幅衰減(Lee et al., 2009)。因地形效應額外產生之表面波與散射波，預期會影響測站處記錄之波形。本研究參考美日合作產出之 GDEM v002 地形模型(NASA/METI/AIST/Japan Spacesystems, and U.S./Japan ASTER Science Team, 2009)，該模型之解析度為 30 m，考量前述北台灣速度構造模型後，經穩定性測試，計算出震波模擬適合之模型解析，將地形模型重採樣為 50 m 取一格點，作為數值計算網格水平方向格距，建立數值網格模型，該格距可確保若震波模擬至 1.5 Hz，單一波長以 10 個數值網格點描述，可有效避免數值頻散(numerical dispersion)效應。本研究參考大尺度三維速度分佈數值，強調並納入構造不均質(heterogeneity)與側向變化(lateral variation)，提升了僅考量無方向性變化的一維速度模型可能影響地震模擬評估之不足，並考量細緻台北盆地三維速度分佈數值，精進盆地區域波傳速度分佈之解析力，合理考量情境模擬工作所需涵蓋空間範圍後，P波與S波整合三維速度分佈數值與地表地形之震波傳遞模擬數值網格如下圖所示。
特徵震源模型建置
考量山腳斷層線型約略於現台北捷運北投機廠位置有所轉折，若以機廠為參考點，往東北方向稱之北段，往西南方向稱之南段，本研究聚焦於精緻化台北盆地區域地震動模擬評估，合理化考量計算資源使用，擇定山腳斷層南段破裂作為震源破裂假想情境，合理考量符合線型之斷層長度與寬度，設定不同破裂模式，進而假定特徵震源模型。
特徵震源模型之參數擬定原則，則採用日本地震調查研究推進本部之 Recipe(Irikura and Miyake, 2011)方法，取其特徵震源擬定之標準化流程作為依據，針對山腳斷層南段情境地震破裂，探討地震動對台北盆地區域影響性，並假定斷層面上具單一地栓(asperity)為主要能量釋放位置，將地栓中心置於沿斷層走向之中線，地栓上緣與斷層上緣重合，並考慮震源(即破裂起始點)位於地栓上緣或下緣之四個不同位置，產製四組震源模型，意即評估不同斷層破裂方向性效應(directivity effect)對台北盆地與鄰近區域之衝擊：
模型 01：震源位於地栓上緣最南端，考量斷層由南往北、朝台北盆地破裂；
模型 02：震源位於地栓上緣中間處，考量斷層朝台北盆地破裂；
模型 03：震源位於地栓上緣最北端，考量斷層由北往南、朝台北盆地破裂；
模型 04：震源位於地栓下緣中間處，考量斷層由深部往淺部、背向台北盆地破裂。
虛擬測站點位
NCDR 點位
本研究配合災損推估小組之評估網格，配合國家災害防救科技中心(National Science and Technology Center for Disaster Reduction, NCDR)格點解析度為 500 公尺，分佈範圍為東經 118.125 度至 122.005 度，北緯 21.895 度至 26.385 度，範圍涵蓋臺灣本島與離島範圍，共 132,712 網格點。下圖為 NCDR 網格示意圖，黑色格點即為網格點位。
本研究設定點位
對應災損推估小組後續評估工作，空間範圍與解析上，本研究基於 NCDR 500 公尺網格點擬定虛擬測站，進行地震動模擬計算。所採行點位分佈範圍為東經 121.325 度至 121.675 度，北緯 24.925 度至 25.275 度，共 4,387 點位，空間點位選取範圍系考慮山腳斷層情境震源所在區域、台北盆地與鄰近區域進而擬定之範圍。如下圖所示，分佈於台北盆地區域之黑色格點為 NCDR 500 公尺網格點位，視為虛擬測站；紅框區域為情境地震動三維有限差分法數值模擬範圍。
地震波模擬分析
低頻模擬：三維有限差分法
本研究之低頻地震波模擬方法，採用有限差分法(finite-difference method)，於三維地形及地下速度構造模型設定中，透過擬定假想地震之情境(特徵震源模型)進行模擬，評估低頻地震波對特定空間點位、工址或測站產生之強地動影響。
震波傳遞之模擬方法係基於物理機制下，考慮地震源與震波傳遞路徑中各種耦合效應，計算特定點位的震動歷時反應，其中震源項包含地震震源機制、斷層幾何與破裂型態，而路徑項則包含不同尺度的地下構造模型(包含波速、密度與衰減因子等)。以有限差分法進行震波模擬，主要係由求解質點運動方程與組構關係隨時間之變化，即為模擬質點受力後於空間中的運動行為與震波隨時空變化之傳遞過程。有限差分法為基於網格運算之數值方法，即為求網格系統中某格點位置之差分運算，運算過程僅需鄰近格點即可計算，因此可透過訊息傳輸介面(Message Passing Interface, Gropp et al., 1996)將數值方法改寫，並在叢集式電腦(PC cluster)進行平行化計算 (parallel computation)，進而提高效率或進行更為高頻的模擬工作。
本研究採用 Zhang et al. (2012)之牽引力鏡像有限差分法(traction-image finite difference method)可考慮(1) 三維地下構造與(2) 地表地形，進行震波模擬。該方法採用 MacCormak 格式(MacCormack, 1969; Gottlieb & Turkel, 1976)與 optimized Dispersion Relation Preserving (DRP/opt, e.g. Tam & Webb, 1993; Hixon, 1997)能有效降低數值頻散(numerical dispersion)效應。同時，該演算法採取完美匹配吸收層(perfectly matched layer, PML)做為吸收邊界條件(absorbing boundary condition, e.g. Berenger, 1994; Marcinkovich & Olsen, 2003)，能有效降低人為設定模型邊界產生之反射波。
為考慮地表地形或地下構造反射面效應，Zhang et al. (2012)透過坐標轉換，將地表地形與地下構造界面以曲線性坐標系統(curvilinear coordinate system)描述，進而使物理空間中各曲線性坐標位置
透過 Jacobian 變換至電腦計算所使用的直角坐標系統
。
高頻模擬：隨機式方法
針對需求所擬定強地動模擬，須計算 PGA、SA0.3 等高頻地震動參數，模擬頻段須達 10 Hz，2.6.1 節之低頻模擬方法無法反應 1-10 Hz 頻段之高頻地震動時間歷時，故本研究計畫選用 EXSIM (Atkinson and Assatourians, 2015)隨機式方法進行高頻地震動模擬，補足人工合成地震歷時之高頻響應。
EXSIM 係根據使用者給定之地震矩，以經驗公式給定斷層長寬形成一斷層面，而後離散化為一系列之子震源(subsource)，其中的每一個子震源皆遵循點震源隨機模型(point-source stochastic model)，依給定的地震歷時(duration)與剪力波輻射型態以隨機高斯噪訊(random Gaussian noise)模擬強地面加速度地動歷時。EXSIM 亦引入 Brune model，於傅氏譜(Fourier spectrum)考慮地震矩、應力參數、經驗衰減參數於頻率域之衰減效應，經快速傅立葉變換(fast Fourier transform, FFT)至時間域後，根據破裂速度疊加各子震源，並考慮場址效應修正函數，形成時間序列。本研究亦參考 Lee et al. (2015)根據場址 VS30 條件分類推估地震歷時之經驗公式，整合於高頻模擬計算
EXSIM 為相當具代表性之隨機式強地動模擬方法，除南加洲地震中心(Southern California Earthquake Center, SCEC)已對 EXSIM 實行相關檢核(e.g.Atkinson and Assatourians, 2015)外，EXSIM 模擬結果亦已與隱沒帶地震事件檢核(e.g. Atkinson and Macias, 2009)。
本研究以 EXSIM 模組進行隨機法強地動模擬時，給定山腳斷層南段破裂特徵震源模型之滑移量分佈作為輸入震源模型，並針對模擬範圍，考慮適用臺灣區域之隨機式模型(Sokolov et al., 2006; Sokolov et al., 2009; Huang et al., 2017)，詳細模擬參數列於下表。
混合法合成寬頻地震歷時
考量前述地栓滑移量分佈之隱沒帶震源情境特徵震源模型，並根據採行之 4,387 網格點位，以三維有限差分法與隨機式分析方法分別模擬各點位之低頻與高頻地震波，並以剪力波波速輔以互相關(cross-correlation)方法修正各點位之震波到時(arrival-time)後疊加低頻與高頻地震波，以工程常用之 VS=760 m/s 基盤面做為地震動計算基準，以混合法合成寬頻地震歷時至該基盤面，後續便以該寬頻震動歷時進行 5 項地動參數計算，包含(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)、(4) 0.3 秒譜加速度(SA0.3)與(5) 1.0 秒譜加速度(SA1.0)。
本研究以叢集式電腦進行平行運算，若以山腳斷層單一情境破裂情境為例，山腳斷層之特徵震源震源模型共 1,320 子斷層與地表面 4,387 需擬測站，其運算量共 5,790,840 次點震源模擬，以 192 計算核心進行計算，低頻模擬約莫耗時 15小時，高頻模耗時約 2 小時，混合法與場址修正耗時約 2 小時。後續若對單一斷層進行多次情境模擬，仍須審慎評估計算量與計算時間，合理使用並擴增計算資源使計算效益最大化。
地震動歷時與反應譜計算分析
地震動參數
依前述地震波模擬分析方法，考慮山腳斷層南段破裂情境地震之特徵震源模型，分別進行低頻與高頻模擬，進而以混合法合成於工程基盤之寬頻加速度地震歷時，爾後分別積分一次得速度歷時、積分二次得位移歷時，各地震歷時之最大震幅處即分別為(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)等地動數值。另一方面，選取加速度地震歷時作為輸入地震歷時，考慮 5%阻尼計算該地震歷時之反應譜(response spectrum)，爾後分別取週期 0.3秒與週期 1.0 秒之譜加速度作為反應譜地動數值；此外，近期氣象局更新震度量級，預訂於民國 109 年元月起實施新案，原由 PGA 訂定 0-7 級震度，現改為由PGA 與 PGV 依條件訂定 0-7 級震度，5 級與 6 級再各細分為 5 弱、5 強、6 弱、6強，便於與災損有更加對應關係，本研究計畫響應氣象局新制震度，模擬震度以新制震度量級呈現。各組模型引致北臺灣地震動分述如下：
模型 01：震源位於地栓上緣最南端，考量斷層由南往北、朝台北盆地破裂，下圖分別為 PGA、震度、PGV、PGD、SA0.3 與 SA1.0 地震動參數於空間分佈。模擬震度逾 7 級的區域包含新北市五股區、蘆洲區、泰山區，新莊區與三重區則有部分區域逾 7 級。此外，從 PGA、PGV、SA0.3 與 SA1.0 等地震動參數之空間分佈，可明顯看出震源破裂的方向性效應，於山腳斷層南段北邊與往北方向有較大地震動，山腳斷層南段南邊與往南方向地震動較小。另外，PGD 無明顯反應地震破裂方向性，但反應 S 波震源輻射型態(radiation pattern)以及盆地區域對地震動之影響。
模型 02：震源位於地栓上緣中間處，考量斷層朝台北盆地破裂，朝台北盆地破裂，分別為 PGA、震度、PGV、PGD、SA0.3 與 SA1.0 地震動參數於空間分佈。模擬震度逾 7 級的區域包含新北市五股區、蘆洲區、泰山區，新莊區與三重區則有部分區域逾 7 級。此特徵震源模型因震源位於斷層面上緣中間位置，故於各地震動參數之空間分佈無法明顯看出震源破裂的方向性效應，但因該區域地下構造效應特性與地栓效應耦合，PGA、PGV、SA0.與 SA1.0 等地震動參數於山腳斷層南段南北邊均形成地震動高區，且北邊分佈範圍較難邊為大。此外，PGD 雖無明顯反應地震破裂方向性，但反應 S 波震源輻射型態以及盆地區域對地震動之影響。
模型 03：震源位於地栓上緣最北端，考量斷層由北往南、朝台北盆地破裂，分別為 PGA、震度、PGV、PGD、SA0.3 與 SA1.0 地震動參數於空間分佈。模擬震度逾 7 級的區域包含新北市五股區、蘆洲區、泰山區，新莊區與三重區則有部分區域逾 7 級。此外，類似於模型 01，從 PGA、PGV、SA0.3 與SA1.0 等地震動參數之空間分佈，可明顯看出震源破裂的方向性效應，於山腳斷層南段南邊與往南方向有較大地震動，山腳斷層南段北邊與往北方向地震動較小，恰與模型 01 所模擬之地震動分佈相反。另外，PGD 無明顯反應地震破裂方向性，但反應 S 波震源輻射型態以及盆地區域對地震動之影響。
模型 04：震源位於地栓下緣中間處，考量斷層由深部往淺部、背向台北盆地破裂，分別為 PGA、震度、PGV、PGD、SA0.3 與 SA1.0 地震動參數於空間分佈。模擬震度逾 7 級的區域包含新北市五股區、蘆洲區、泰山區，新莊區與三重區則有部分區域逾 7 級。此特徵震源模型因震源位於地栓下緣中間位置，但由於震源破裂為由深部往淺部破裂，故由 PGA、PGV、PGD、SA0.3 與SA1.0 等地震參數皆可於斷層下盤方向可見較大地震動，此模型亦可為斷層破裂方向性效應之說明案例。此外，同前述，PGD 除了於斷層下盤方向顯示較大地震動(如新北市八里區)，亦反應 S 波震源輻射型態以及盆地區域對地震動之影響。
震度分布範圍
作為方災策略擬定、更甚者救災考量，彙整單一震源多重破裂情境模擬之地震動高區聯集與交集，可更加明確判斷地震動分佈，後續亦較容易應用於災防策略推演。下圖為震源模型於 NCDR 網格點上所引致之特定模擬震度之聯集與交集(震度≥5 弱、震度≥5 強、震度≥6 弱)，圖中清楚可見山腳斷層線型沿線所通過行政區皆有震度≥6 弱之地震動潛勢，且由於斷層位態(走向與傾角方向)朝向台北盆地，另受盆地效應影響，亦造成台北市西邊行政區皆有較大震度；此外，震度 6 強與震度 7之交集區域，可明顯見到本研究所設置地栓、近山腳斷層南段轉折之區域(新北市五股區、蘆洲區、泰山區、板橋區)為須重點關注區域，各模型皆造成這些行政區有 6 強以上之地震動潛勢；而不同破裂模式將造成南往桃園市龜山區、北往新北淡水區、三芝區、石門區皆有可能有較大模擬震度發生之情事。
地震波於盆地傳播效應
本研究考量細緻化台北盆地及鄰近區域地下構造模式，與前期之構造模型相比較，精進原有之地下構造模型，並採用更高解析度之地形模型。本研究再行採用前期計畫山腳斷層南段破裂情境之特徵震源模型，進行低頻模擬計算，爾後與前期產製之地震歷時波形比較，詳實量化有無盆地效應對地震波形之影響
性。
為求量化不同區域盆地效應之影響性，選取比較的測站有三，包含：(A)位於台北盆地西緣，且最靠近山腳斷層地表線型，該測站位於台北盆地基盤最深區域正上方；(B) 位於盆地中央，預期震波於盆地中複反射與共振效應最為強烈；(C)位於盆地東緣，較為遠離震源區，且位於台北盆地基盤較淺與平緩區域。下圖為 A、B、C 三站模擬地震歷時，顯示台北盆地構造對地震波有明顯放大效應，S 波震幅明顯放大(測站 A 與 B)，物理因素為松山層與基盤之波傳速度與密度差異形成較大的波阻抗對比(impedance contrast)，是影響地震波放大之主因；台北盆地之碟形幾何形貌呈現基盤面於盆地西緣較深、往東較淺之分佈形態，亦使地震波於盆地中造成來回震盪與建設性干涉，加劇地震波於盆地內放大效應，並且拉長地震延時(duration)(測站 B)；此外，台北盆地四周由山脈與台地圍繞，地形效應所產生之散射波亦為貢獻放大效之次要因素(測站 A、B 與 C)。綜觀上述階段性模擬成果，顯示精緻化震源與地下構造參數，對震波傳遞模擬有直接助益，逐步實現以模擬方法量化地震動之可行性。
成果總結
本研究聚焦於精緻化台北盆地區域地震動模擬評估，合理化考量計算資源使用，擇定山腳斷層南段破裂作為震源破裂假想情境，並設定不同破裂模式，以期有效量化臨近台北盆地區域之情境地震引致潛勢地震動，並可近一步作為災損推估之關鍵輸入資訊。盤點本計畫關鍵執行工作項目，包含：
- 山腳斷層震源參數彙整：彙整台灣地震模型組織近期對於山腳斷層之工作成果，並對山腳斷層震源參數進行更新；
- 地形與速度構造模型建置：彙整大尺度三維速度分佈數值、台北盆地 P波與 S 波三維速度分佈數值，以及 GDEM 數值地形模型；
- 數值網格穩定測試：依彈性波理論與數值穩定條件，嚴格計算數值網格對於模擬結果之穩定性；
- 特徵震源模型建置：設定山腳斷層南段破裂情境，透過日本 Recipe 方法，擬定四種破裂模式並建立其巨觀與微觀震源參數；
- 透過數值模擬方法進行各地震動參數計算與比較；
- 與前期成果進行有無盆地效應之地震波形比較。
基於上述參數率定與模型建立流程，考量系統性震源參數設置與細緻化台北盆地及鄰近區域地下構造模式，實現山腳斷層南段破裂之情境地震模擬，綜觀上述階段性模擬成果，顯示精緻化震源與地下構造參數，對震波傳遞模擬有直接助益，逐步實現以模擬方法量化地震動之可行性。藉由本研究產出之地震動參數作為災損推估模式之關鍵輸入參數，近程目標係為檢核與精進目前災損推估模式之可行手段，協助檢核現今國內災損推估與規範採用之建物易損性曲線(fragility curve)是否如預期，長遠目標則可應用於緊急應變與災損策略之擘劃與修訂，以及地震災害應變演練，防範地震災害於未然。
然而，本研究尚僅考慮單一地栓震源破裂情境，倘若山腳斷層南段破裂，不全然為本研究所假想之震源破裂型態以及推估規模；此外，場址效應對於臺北市區之地震動放大、土壤非線性效應對於地震動之折減，兩者之耦合關係仍待釐清。故本研究報告建議如下：
- 針對單一震源，仍須考量不同地栓分佈與震源破裂起始位置，進階量化震源破裂方向性效應，針對更多組之震源破裂情境進行模擬，並對產出之地震動參數進行統計分析，有效量化地震動中值與極值，估算地震動範圍。
- 詳細化目標地區之場址放大效應與土壤非線性效應，高解析之場址特性對於提升地震動空間解析有相當助益。
參考文獻
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.
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.
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., Chen, H. W., Liu, Q., Komatitsch, D., Huang, B. S., and Tromp, J. (2008). Three dimensional simulations of seismic-wave propagation in the Taipei Basin with realistic topography based upon the spectral-element method. Bull. Seism. Soc. Am., 98(1), 253-264. doi:10.1785/0120070033.
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.
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.
NASA/METI/AIST/Japan Spacesystems, and U.S./Japan ASTER Science Team (2009).
ASTER Global Digital Elevation Model [Data set]. NASA EOSDIS Land Processes DAAC. Accessed 2019-12-03 from https://doi.org/10.5067/ASTER/ASTGTM.002
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).
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.
行政院災害防救辦公室(民 105)。105 年災害防救白皮書，行政院，191 頁，臺灣臺北。
陳亮全、劉楨業、陳宏宇、柯孝勳、張國鎮、許健智、林正洪、馬國鳳與劉兆漢
(中央研究院大規模地震災害防治策略研議小組) (民 104)。大規模地震災害防治策略建議書，中央研究院，44 頁，臺灣臺北。