國立中央大學地震災害鏈風險評估及管理研究中心
大規模地震情境模擬對外資料
2020-中洲構造特徵震源地動模擬
本研究「臺灣南部都會區特定斷層潛勢之大規模地震情境模擬地震動評估」作為輔助執行「大規模地震衝擊評估模型強化與應用面建構」研究計畫，採用震源情境地震動模擬分析流程進行震源參數擬定與地震動模擬計算，本期(109 年度)考量位於台灣南部、具高發震潛勢之中洲構造錯動破裂，影響地表面地震動之情境，並納入(1) 極端地震動情境與(2) 中值地震動情境作為考量，最終提供可作為災損分析輸入參數之峰值地動加速度(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)，其地震動數值涵蓋空間範圍包含台南都會全區。
參數與方法
此處闡述考量中洲構造破裂之情境地震而彙整之震源參數，並以此作為依據，以地震學、地體動力學與工程地震學角度，考慮對台南都會區較具衝擊之情境，擬定特徵震源模型。另對於速度、地形構造之建構進行考量。最終針對該單一特徵震源模型進行情境模擬地震動，模擬完整情境地震歷時，計算地震動參數，包含(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)、(4) 0.3 秒譜加速度(SA0.3)與(5) 1.0 秒譜加速度(SA1.0)等五項參數作為產出，以利後續災損推估小組推估台南都會區震災情境。綜整上述說明，模擬分析程序詳見下圖。
中洲構造概述
1949 年中國石油公司透過重力探測，發現台南平原地區有多處重力異常區域，後續透過震測資料、航空攝影測量與重力測勘等技術，證實了台南台地與中洲階地有兩處重力異常區，同時，鑽井資料也證實了地下構造的存在(Hsieh,1972; Sun, 1964)。在地貌特徵的部分，透過光學雷達(Light Detection and Ranging, LiDAR)資料繪製而成的三維地形模型來看，可看出台南台地兩側有明顯高低落差，且台南台地東側明顯比西側陡峭(經濟部中央地質調查所, 2013)，Hsieh (1972)推測可能為泥貫入體形成了此處的背斜構造。此後，Chen and Liu(2000)透過分析台南地區的抬升與沉降速率，發現在全新世時期，抬升速率最快的區域在台南台地與中洲階地兩處，分別 5 mm/yr 與 7 mm/yr，雖然同為抬升速率最快的區域，但 Chen and Liu (2000)認為這兩區域的成因並不相同，從地形判釋上來看，台南台地的成因可能受到泥貫入體所影響，而中洲構造的部分，因其抬升速率有向東增加的趨勢，因此，Chen and Liu (2000)推論，中洲階地下方的中洲構造才是造成抬升的主要機制。雖然中洲構造不屬於中央地質調查所於 2012 年公佈的活動斷層之一，但因其仍具有孕震潛勢，故被台灣地震模型(Taiwan Earthquake Model, TEM)視為孕震構造。
台灣位處於歐亞板塊與菲律賓海板塊交界帶，頻繁的碰撞使得台灣地震活動頻仍，為了有效降低地震災害所帶來的風險，TEM 致力於建立台灣的地震參數資料庫，不僅關注已有明顯地表破裂痕跡或明顯位移量的活動斷層，同時也針對具有孕震潛能，但僅有地貌特徵且缺乏地表破裂的孕震構造進行孕震構造參數的收集，以評估台灣的地震災害風險。
中洲構造位於台南台地與大灣低地東側，北起自台南市官田區，南至台南市關廟區，為一長度為 29.2 公里，傾角為 30°，深度 12 公里的逆斷層(Shyu et al., 2016; Shyu et al., 2020)(藍線為中州構造位置所在處，黃色圓圈為台南市政府所在位置。
Sun (1964)透過航測照片，判釋出台南至高雄間許多的褶皺構造與斷層構造，其中便包括了中洲背斜，而從 LiDAR 測量資料繪製而成的三維地形模型來看，中洲構造兩側地形高低落差極大，也可推測出中洲構造存在的位置(經濟部中央地質調查所, 2013)，並進而繪製出中洲構造線型(Shyu et al., 2016)。Chen and Liu (2000)經鑽井資料來計算出中洲上升側的抬升速率構造約為 5-8mm/yr。Shyu et al. (2016)透過前人研究與地貌特徵來估算長期滑移速率，得出中洲構造平均滑移速率約 12.2 mm/yr (9.02mm/yr-18.71mm/yr)，再現週期約為 50-190 年，破裂面積預估可達 852 平方公里 (Shyu et al., 2020)。在最大地震矩規模方面，Shyu et al. (2020)採用了兩種不同的公式來做計算，基於全世界歷史地震的震源參數，推導出的回歸式所得出的 MW 為 6.97 (Wells and Coppersmith, 1994)，而台灣位處於板塊交界帶，多以逆衝斷層為主，考慮到台灣特有的條件後所推導出的回歸式則是 MW 為 6.93 的數值 (Shaw, 2009; Yen and Ma, 2011)。
TEM 也依據泊松模型(Poisson process)計算各孕震構造的發震機率，在2016 年公布的未來 30 年台灣孕震構造之發震機率圖中(Wang et al., 2016)，中洲構造的發震機率最大規模為 MW6.9，最大規模破裂機率為 24%，且其位置接近人口稠密的台南都會區，對於鄰近區域及南部科學園區影響甚大，故本研究將以中洲構造作為地震動預估與災損評估的主要對象。
地形與速度構造模型建置
地震波傳模擬所需速度構造模型建置，係以近期台灣地區速度構造研究(Kuo-Chen et al., 2012)之大尺度三維速度分佈數值，作為本研究模擬計算之模型建置基礎資料。該速度分布數據係藉由中央氣象局強震站及短週期地震站、中央氣象局及中央研究院寬頻地震站、大量布設臨時站、海底地震儀所記錄到之大量地震波走時資訊，經由逆推方法計算建構而成。本研究另行檢核三維速度構造各位置與深度點位之 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)，提升了僅考量無方向性變化的一維速度模型可能影響地震模擬評估之不足，並考量細緻台北盆地三維速度分佈數值，精進盆地區域波傳速度分佈之解析力，合理考量情境模擬工作所需涵蓋空間範圍後，整合三維速度分佈數值與地表地形之震波傳遞模擬數值網格如圖所示。
特徵震源模型建置
考量中洲構造線型於新市區與山上區交界有較大轉折，根據斷層動力學模式，轉折處可能為地栓(asperity)區域，於地震發生時可能引致較大地表位移、輻射較大地震波能量。本研究進行地震動模擬評估，合理化考量計算資源使用，簡化中洲構造線型並考量符合線型之斷層長度與寬度，設定破裂模式，進而假定特徵震源模型。
特徵震源模型之參數擬定原則，採用日本地震調查研究推進本部之 Recipe(Irikura and Miyake, 2011)方法，取其特徵震源擬定之標準化流程作為依據，針對中洲構造情境地震破裂，探討地震動對台南都會區影響性，並假定斷層面上具不同面積分布之雙地栓(asperity)特徵，該雙地栓為主要地震波能量釋放位置，較大面積之地栓鄰近中洲構造最大轉折處，較小面積之地栓則位於中洲構造南段。考量地栓分佈後，設定兩種震源情境探討對台南都會區之地震動影響：
- 情境破裂模型 01：大小面積二地栓上緣與地表重合，二地栓上緣中點分別位於地表簡化線型三等分點處，破裂起始位置(即震源位置)位於較大地栓角落靠近斷層面中央處，係為 5.30 公里深。此情境考量地震主要能量釋放位置位於地表，為極端地震動情境；
- 情境破裂模型 02：大小面積二地栓中心分別位於全斷層面深度中線之三等分點處，破裂起始位置(即震源位置)位於較大地栓單邊中點靠近斷層面中央處，係為 5.90 公里深，此設計情境為考量中值地震動情境。
此外，本研究假設二種情境破裂模型之震源位置皆靠近斷層面中央，考量破裂方向為均向(isotropic)，無顯著方向性效應(directivity effect)。
二組中洲構造情境破裂震源模型以平面示圖見下圖；斷層面之巨觀與微觀震源參數則分列於下表。
虛擬測站點位
NCDR 點位
本研究配合災損推估小組之評估網格，配合國家災害防救科技中心(National Science and Technology Center for Disaster Reduction, NCDR)格點解析度為 500 公尺，分佈範圍為東經 118.125 度至 122.005 度，北緯 21.895 度至 26.385 度，範圍涵蓋臺灣本島與離島範圍，共 132,712 網格點。圖為 NCDR 網格示意圖，黑色格點即為網格點位。
本研究設定點位
對應災損推估小組後續評估工作，空間範圍與解析上，本研究基於 NCDR 500公尺網格點擬定虛擬測站，進行地震動模擬計算。所採行點位分佈範圍為東經120.005 度至 121.700 度，北緯 22.850 度至 23.45 度，共 14,431 點位，空間點位選取範圍系考慮中洲構造情境震源所在區域以及台南都會區範圍進而擬定。如下圖所示，分佈於台南都會區之黑色格點為 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 and Turkel, 1976)與 optimizedDispersion Relation Preserving (DRP/opt, e.g. Tam and Webb, 1993; Hixon, 1997)能有效降低數值頻散(numerical dispersion)效應。同時，該演算法採取完美匹配吸收層(perfectly matched layer, PML)做為吸收邊界條件(absorbing boundary condition, e.g. Berenger, 1994; Marcinkovich and Olsen, 2003)，能有效降低人為設定模型邊界產生之反射波。
為考慮地表地形或地下構造反射面效應，Zhang et al. (2012)透過坐標轉換，將地表地形與地下構造界面以曲線性坐標系統(curvilinear coordinate system)描述，進而使物理空間中各曲線性坐標位置
透過 Jacobian 變換至電腦計算所使用的直角坐標系統 
。
高頻模擬：隨機式方法
針對需求所擬定強地動模擬，須計算 PGA、SA0.3 等高頻地震動參數，模擬頻段須達 10 Hz，上述之低頻模擬方法無法反應 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)，詳細模擬參數列於下表。
混合法合成寬頻地震歷時
考量前述地栓滑移量分佈之隱沒帶震源情境特徵震源模型，並採行之 14,431 網格點位，以三維有限差分法與隨機式分析方法分別模擬各點位之低頻與高頻地震波，並以剪力波波速輔以互相關(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)。
本研究以叢集式電腦進行平行運算，若以中洲構造單一情境破裂情境為例，中洲構造之特徵震源模型共 4,380 子斷層與地表面 14,431 虛擬測站，其運算量共 63,207,780 次點震源模擬，以 64 計算核心進行計算，低頻模擬約莫耗時 5 小時，高頻模耗時約 3 小時，混合法與場址修正耗時約 3 小時。後續若對單一斷層進行多次情境模擬，仍須審慎評估計算量與計算時間，合理使用並擴增計算資源使計算效益最大化。
地震動參數分布與反應譜計算分析
地震動參數
依前述地震波模擬分析方法，考慮中洲構造破裂情境地震之特徵震源模型，分別進行低頻與高頻模擬，進而以混合法合成於工程基盤之寬頻加速度地震歷時，爾後分別積分一次得速度歷時、積分二次得位移歷時，各地震歷時之最大震幅處即分別為(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)等地動數值。另一方面，選取加速度地震歷時作為輸入地震歷時，考慮 5%阻尼計算該地震歷時之反應譜(response spectrum)，爾後分別取週期 0.3 秒與週期 1.0 秒之譜加速度作為反應譜地動數值；此外，另以氣象局新制震度量級同時考量 PGA 與 PGV，呈現模擬結果。各組模型引致台南都會區地震動分述如下：
情境破裂模型 01：此情境考量地震主要能量釋放位置位於地表，為極端地震動情境，下圖分別為 PGA、震度、PGV、PGD、SA0.3 與 SA1.0 地震動參數於空間分佈。模擬震度逾 7 級的區域皆位於地表斷層線型附近偏上盤側，包含官田、大內、善化、山上、新市、新化、關廟等行政區。此外，PGA 與 PGV地震動參數之空間分布呈現震源區往台南西北與西南方向有較高地震動數值，該數值分布與地栓位置以及情境破裂帶朝向斷層面 up-dip 方向有關；PGD 地震動參數分布多反應長週期地震動，模擬 PGD 呈現 S 波震源輻射型態(radiation pattern)分布；SA0.3 地震動參數分布則反應相對短週期地震波響應，量值高區多位於斷層線型偏上盤側，與 PGA 分布相似；SA1.0 地震動參數分布則反應相對長週期地震波響應，除斷層線型偏上盤側為量值高區外，震源往西南方向亦有明顯較大譜加速度值。
情境破裂模型 02：考量雙地栓位於斷層面中心附近，此設計情境為考量中值地震動情境，整體地震動參數量值較情境破裂模型 01 為小，分別為 PGA、震度、PGV、PGD、SA0.3 與 SA1.0 地震動參數於空間分佈。模擬震度逾 7 級的區域皆位於地表斷層線型附近偏上盤側，包含官田、大內、關廟等三行政區，分布範圍較情境破裂模型 01 之地震動有限且近地栓位置；PGV 地震動參數分布亦與地栓位置有所關連，呈現地栓區域有較大地震動量值；PGD 地震動參數分布亦如同情境破裂模型 01，與 S 波震源輻射型態較有關聯；SA0.3 與 SA1.0地震動參數量值高區多分布於斷層線型沿線，形成二狹長帶狀分布，亦與設定地栓位置有所關連，而 SA1.0 亦顯示震源往西南方向亦有明顯較大譜加速度值。
上述模擬成果顯示，在未考慮顯著震源破裂方向性效應情況下，地栓分布位置與 S 波輻射型態扮演地震動參數分布極為重要的角色；此外，地栓於地下構造中的深度位置，亦影響於地表面的地震動參數分布。
震度分布範圍
作為防災策略擬定、更甚者救災考量，彙整單一震源多重破裂情境模擬之地震動高區聯集與交集，可更加明確判斷地震動分佈，後續亦較容易應用於災防策略推演。下圖為震源模型於 NCDR 網格點上所引致之特定模擬震度之聯集與交集，因情境破裂模型 01 與 02 主要差異來自於地栓深度分布，於(1)震源破裂起始位置相同且(2) 破裂速度一致之情況下，二者模型產出之地震動參數分布特徵相近，而量值分布主要由地栓位置與 S 波震源輻射型態所控制，故細究二者震度分布範圍，「震度交集」分布呈現由模擬模型 02 所主導，而「震度聯集與交集」所形成的分布則呈現由模擬模型 01 所主導。
根據各圖呈現結果，兩種情境將引致台南全區達震度 5 弱；由於民國 109 年起，我國中央氣象局更新震度計算方式並同時採用 PGA 與 PGV，因此在震度 5 強分布區域明顯為反應 PGV 分布，與震源 S 波輻射型態相近；震度 6 弱分布區域明顯靠近台南都會區西部，雖模型 01 與 02 無顯著震源破裂方向性，但由於模型 01 地栓位於近地表，發震起始點(即震源)往地栓破裂之方向帶有 up-dip 分量方向性，故有部分方向性效應凸顯於震度 6 弱分布，加成 S 波輻射型態影響；震度 6 強交集區域顯示模型 01 與模型 02 於地栓附近為引致震度高區，震度 6 強聯集區域顯示模型 01 之 up-dip 分量方向性效應造成高震度往台南西部延伸；震度 7則顯示高震度區即地栓往地表方向投影、位於中洲構造線型處。此外，因台南都會區西部下方為深厚沉積層，雖大尺度三維地下速度構造無法完全反應該沉積層對地震波放大效應之貢獻，但對於長週期地震動參數，如 PGV、PGD、SA1.0 等，仍具相當程度貢獻，有待後續釐清沉積層對地震動之影響性。
成果總結
本研究聚焦於台南都會區情境地震動模擬評估，並合理化考量計算資源使用，設定二種不同破裂模式，以期有效量化不同震源模式造成台南都會區之潛勢地震動，並可近一步作為災損推估之關鍵輸入資訊。盤點本計畫關鍵執行工作項目，包含：
- 中洲構造震源參數彙整：彙整台灣地震模型組織近期對於中洲構造之工作成果；
- 地形與速度構造模型建置：彙整大尺度三維速度分佈數值，包含 P 波與S 波三維速度分佈數值，以及 GDEM 高解析數值地形模型；
- 數值網格穩定測試：依彈性波理論與數值穩定條件，嚴格計算數值網格對於模擬結果之穩定性；
- 特徵震源模型建置：設定中洲構造破裂情境，透過日本 Recipe 方法，擬定二種破裂模式並建立其巨觀與微觀震源參數；
- 透過數值模擬方法進行各地震動參數計算與比較。
基於上述參數率定與模型建立流程，考量系統性震源參數設置與建構台南都會區之地下構造模式與數值地形網格，實現中洲構造破裂之情境地震模擬，綜觀上述階段性模擬成果，顯示精緻化震源與地下構造參數，對震波傳遞模擬有直接助益，逐步實現以地震波數值模擬方法量化地震動之可行性。藉由本研究產出之地震動參數作為災損推估模式之關鍵輸入參數，近程目標係為檢核與精進目前災損推估模式之可行手段，協助檢核現今國內災損推估與規範採用之建物易損性曲線(fragility curve)是否如預期，長遠目標則可應用於緊急應變與災損策略之擘劃與修訂，以及地震災害應變演練，防範地震災害於未然。
然而，本研究目前僅考慮二種震源破裂情境，倘若中洲構造發生實際破裂並引致地震，不全然為本研究所假想之震源破裂型態以及推估規模；此外，誠如「3.2震度分布範圍」一節所述，大尺度三維速度構造並非完全能反應沉積層以及地表極淺層場址效應對於台南都會區之地震動影響，以及土壤非線性效應對於地震動之折減。基於上述研究議題，本研究報告建議如下：
- 針對地震動情境模擬之震源設定，若要詳實考量地震動極端值與中間值分布，仍待考量不同地栓分布型態，包含地栓大小尺度、地栓分配位置等；此外，設定不同震源破裂起始位置，並考量不同破裂速以及各子斷層之震源時間函數，意函納入震源破裂方向性效應因素。藉由探討上述震源特性，由運動學方法模擬分析震源動力學模式，可進一步掌握震源效應對於目標區域之分布型態影響性。
- 建構目標區域沉積層至地表極淺層覆土之地下構造物性參數，並探討目標區域之頻率相依(frequency-dependent)之地震基盤地震動轉換至工程基盤地震動、淺層場址效應造成震波振幅放大，以及因土壤非線性效應所造成震波振幅折減，對於提升地震動空間解析力，或更進階之土壤液化潛勢分析有相當程度的助益。
參考文獻
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.
Chen, Y.-G., and Liu, T.-K. (2000). Holocene uplift and subsidence along an active
tectonic margin southwestern Taiwan, Quat. Sci. Rev., 19(9), 923-930.
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.
Hsieh, S. H. (1972). Subsurface geology and gravity anomalies of the Tainan and
Chungchou structures of the coastal plain of southwestern Taiwan (臺灣西南部濱海平原台南及中洲兩構造之地下地質與重力異常之研究), Petroleum
Geology of Taiwan, 10, 323–338.
Hixon, R. (1997). On increasing the accuracy of MacCormack schemes for aero acoustic
applications, AIAA Paper, doi: 10.2514/6.1997–1586.
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., Jenkins, D. M., Mechie, J., Roecker, S. W., Wang, C. Y., and
Huang, B. S. (2012),Seismic evidence for the α-β quartz transition beneath Taiwan from Vp/Vs tomography, Geophys. Res. Lett., 39, L22302, doi:10.1029/2012GL053649
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
Shaw, B. E. (2009). Constant stress drop from small to great earthquakes in magnitude area scaling, Bull. Seism. Soc. Am., 99(2A), 871–875.
Shyu, J. B. H., Chuang, Y. R., Chen, Y. L., Lee, Y. R., and Cheng, C. T. (2016). A New On Land Seismogenic Structure Source Database from the Taiwan Earthquake Model (TEM) Project for Seismic Hazard Analysis of Taiwan, Terr. Atmos. Oceanic. Sci., 27(3), 311–323.
Shyu, J. B. H., Yin, Y. H., Chen, C. H., Chuang, Y. R., and Liu, S. C. (2020). Updates to the
on-land seismogenic structure source database by the Taiwan Earthquake Model (TEM) project for seismic hazard analysis of Taiwan, Terr. Atmos. Oceanic. Sci.,31(4), 469–478.
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).
Sun, S. C. (1964). Photogeologic study of the Tainan-Kaohsiung coastal plain area. (臺
灣省臺南至高雄間平原區域航照地質之研究), Petroleum Geology of Taiwan,3, 39-51.
Tam, C. and Webb, J. C. (1993). Dispersion-Relation-Preserving finite difference
schemes for computational acoustics, J. Comput. Phys., 107, 262–281.
Wang, Y. J., Chan, C.-H., Lee, Y.-T., Ma, K.-F., Shyu, J. B. H., Rau, R.-J., and Cheng, C.-T.
(2016). Probabilistic seismic hazard assessments for Taiwan, Terr. Atmos. Ocean. Sci., 27(3), 325–340.
Wells, D. L., and Coppersmith, K. J. (1994). New empirical relationships among
magnitude, rupture length, rupture width, rupture area, and surface displacement, Bull. Seism. Soc. Am., 84(4), 974–1002.
Yen, Y. T., and Ma, K. F. (2011). Source-scaling relationship for M4.6–8.9 earthquakes,
specifically for earthquakes in the collision zone of Taiwan, Bull. Seism. Soc. Am.,101(2), 464–481.
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 頁，臺灣臺北。
經濟部中央地質調查所. (民 102)。重要活動斷層調查特性研究–近斷層高精度地
形資料之判釋與分析(3/4)，經濟部中央地質調查所報告第 102-08 號。