ABSTRACT Title of Thesis: ASSESSING FAULT SLIP HAZARD IN TAIWAN USING SPACE GEODESY Kathryn Rose Robbins, Master of Science, 2022 Thesis Directed By: Assistant Professor Mong-Han Huang, Department of Geology Taiwan is a geologically complex region due to the continuous collision of the Eurasian Plate and the Philippine Sea Plate. This study aimed to quantify the interseismic crustal deformation of Taiwan and detail the island?s seismic hazard potential using space geodesy. Data were collected between 2016 and 2021 through C-band Copernicus Sentinel-1 synthetic aperture radar imagery and continuous GNSS data from Academia Sinica, Taiwan. I excluded major earthquake events within this time period and generated a dataset conssisting of interferometric synthetic aperture radar ground motion velocities with GNSS corrections and interpolated GNSS ground motion velocities. Then, utilizing this dataset, I performed a deformation rate analysis and error analysis. Next, I explored block modeling and used a total variation regularization approach to determine the reference block model that best reduced velocity residuals and minimized the number of independently rotating blocks. Results suggested that the Taipei Basin, Ilan Basin, Western Foothills, and Longitudinal Valley were experiencing increased total strain rate accumulation and, therefore, posed increased seismic hazard. ASSESSING FAULT SLIP HAZARD IN TAIWAN USING SPACE GEODESY by Kathryn Rose Robbins Thesis submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Master of Science 2022 Advisory Committee: Assistant Professor Mong-Han Huang, Chair Professor Vedran Leki? Professor Laurent Mont?si ? Copyright by Kathryn Rose Robbins 2022 Acknowledgements I would like to express my deepest gratitude to Dr. Mong-Han Huang for his patience, knowledge, and endless supply of coffee. Thank you for introducing me to the world and possibilities of remote sensing. Many thanks to Dr. Eileen Evans as well, for sharing her block model knowledge and time. I am also thankful to Paul and Figgy Franklin, for keeping my spirits high and providing me with boundless love and laughter. Sara Attayan, for knowing me better than I know myself and for giving selfless guidance. Alice and Allan Robbins, for encouraging me to follow my dreams. My lab and cohort friends, for making sure my studies were fueled by plentiful brunches. Finally, I would like to express gratitude to Pragnyadipta Sen, for helping me find my passion for geology. This research is supported by the National Science Foundation award EAR-2026099. ii Table of Contents Acknowledgements ....................................................................................................... ii Table of Contents ......................................................................................................... iii List of Tables ................................................................................................................ v List of Figures .............................................................................................................. vi List of Abbreviations ................................................................................................... ix Chapter 1: Introduction ................................................................................................. 1 1.1 Components of the Earthquake Cycle................................................................. 1 1.2 Geologic Overview of Taiwan ............................................................................ 2 1.3 Purpose of Study ................................................................................................. 6 1.4 Hypotheses & Significance of Study .................................................................. 7 Chapter 2: Data and Methods ....................................................................................... 8 2.1 Synthetic Aperture Radar and GNSS Data Collection ....................................... 8 2.2 Geodetic Technique: Interferometric Synthetic Aperture Radar ........................ 9 2.2.1 Limitations ................................................................................................. 12 2.3 Data Analysis Process for InSAR ..................................................................... 13 2.3.1 GNSS Time-Series Adjustment and Mean Velocity Estimation ............... 13 2.3.2 InSAR Scientific Computing Environment Stack Processor Module ....... 15 2.3.3 Miami InSAR Time-Series Software in Python (MintPy) ......................... 17 2.3.4 Time-Series Re-Model & Mean Velocity with GNSS Correction ............ 19 2.3.5 Merge Ascending and Descending from East and West Taiwan ............... 23 2.3.6 Convert GNSS-Corrected LOS InSAR & GNSS Velocities ..................... 23 2.3.7 Mean Velocity Gradient ............................................................................. 25 2.4 Deformation Rate Analysis ............................................................................... 26 2.4.1 Dilatation, Maximum Shear, and 2nd Invariant .......................................... 26 2.5 Error Analysis ................................................................................................... 30 2.5.1 Solving for RMS Misfit & ? ...................................................................... 30 2.5.2 Solving for Uncertainty.............................................................................. 31 2.5.3 Noise Structure Contributions in the Deformation Rate Analysis ............. 32 2.6 Geodetic Technique: Block Modeling .............................................................. 34 2.6.1 Limitations ................................................................................................. 37 2.7 Data Analysis Process for Block Model ........................................................... 37 2.7.1 Map Faults on QGIS .................................................................................. 37 2.7.2 Block Model Construction and Evaluation ................................................ 37 Chapter 3: Results ....................................................................................................... 41 3.1 InSAR Results ................................................................................................... 41 3.2 Deformation Rate Analysis Results .................................................................. 50 3.3 Error Analysis Results ...................................................................................... 54 3.4 Block Model Results ......................................................................................... 60 Chapter 4: Discussion ................................................................................................. 69 4.1 Combining Geodetic Techniques ...................................................................... 69 4.2 Evaluating Data and Hazard ............................................................................. 70 4.2.1 Data Comparison ....................................................................................... 70 4.2.2 Hazard Comparison ................................................................................... 73 iii 4.2.3 Block Model Comparisons ........................................................................ 80 4.3 Implications: Detected Potential Locked Faults ............................................... 83 Chapter 5: Conclusion................................................................................................ 88 Bibliography ............................................................................................................... 91 iv List of Tables Table 1: Relationship between pixel count and physical distance Table 2: Details on faults present in Taiwan v List of Figures Figure 1: Strain accumulation on fault over time Figure 2: Schematic block diagram of Taiwan (from Angelier et al., 2000, after Angelier et al., 1986) Figure 3: Earthquake map comparison (USGS, 2021) Figure 4: Tectonic and geologic layout of Taiwan Figure 5: Example of phase and amplitude from ENVISAT Figure 6: Ascending and descending flight paths and associated geometries Figure 7: Concept of phase change over time (based on Funning, 2021) Figure 8: GNSS time-series modeling and adjustment Figure 9: Processing of SAR phase values to InSAR Figure 10: Schematic of principal strains and corresponding strain tensor components Figure 11: Components of interseismic creep Figure 12: Schematic of block reduction technique (from Evans et al., 2015) Figure 13: Ascending and descending interferogram networks Figure 14: Ascending and descending MintPy mean velocities Figure 15: Uncorrected ascending and descending LOS InSAR velocities Figure 16: Corrected ascending LOS InSAR velocities, west and east Taiwan Figure 17: Corrected descending LOS InSAR velocities, west and east Taiwan Figure 18: Merged LOS InSAR velocities Figure 19: FIG dataset, final GNSS-corrected InSAR / interpolated GNSS velocities Figure 20: East-west, north-south, and vertical mean velocity gradients Figure 21: Deformation rate analysis with various smoothing windows vi Figure 22: Deformation rate analysis in 144-pixel resolution Figure 23: RMS misfit for time-series modeling Figure 24: RMS misfit for various ? values used in time-series modeling Figure 25: Uncertainty from 3-D velocity transformation Figure 26: Variance between each and all pixels in non-deforming region Figure 27: Downsampled semi-variogram and covariogram model Figure 28: Faults utilized in block model construction Figure 29: Block model construction with locking depths and dip angles Figure 30: Block model evaluation Figure 31: Slip rates associated with reference block model Figure 32: Uncertainty values associated with reference block model Figure 33: Change in descending LOS InSAR velocities between 2016-2021 and 1995-1999 (Huang et al., 2016) Figure 34: 2nd invariant of strain rate throughout time (Yu et al., 1997; Lin et al., 2010; Huang et al., 2016) Figure 35: Change in 2nd invariant of strain rate between 2016-2021 and 2003-2005 (Lin et al., 2010) Figure 36: TEMP PSHA2020 (Chan et al., 2020) main conclusions comparison Figure 37: Faults associated with block model construction and total slip rates from reference block model Figure 38: Evidence for active Taipei fault in Taipei Basin Figure 39: Evidence for potential fault in Western Foothills, Liuchia fault or extension of the Muchiliao ? Liuchia fault vii Figure 40: Evidence for potential fault in Western Foothills, related to Fengshan structure viii List of Abbreviations ASF ? Alaska Satellite Facility CGS ? Central Geologic Survey DEM ? Digital Elevation Model ECMWF ? European Centre for Medium-Range Weather Forecasts ENVISAT ? Environmental Satellite ERS - European Remote-Sensing ESA ? European Space Agency FIG ? Final InSAR / GNSS GAM ? Global Atmospheric Model GNSS ? Global Navigation Satellite System GPS ? Global Positioning System ISCE ? InSAR Scientific Computing Environment InSAR ? Interferometric Synthetic Aperture Radar JPL ? Jet Propulsion Lab LOS ? Line-of-Sight MLC ? Multi-Look Complex MintPy ? Miami InSAR Time-Series Software in Python MRV ? Mean Residual Velocity NASA ? National Aeronautics and Space Administration RMS ? Root Mean Square SAR ? Synthetic Aperture Radar SLC ? Single-Look Complex SRTM ? Shuttle Radar Topography Mission ix TEM PSHA - Taiwan Earthquake Model of Probabilistic Seismic Hazard Analysis TVR ? Total Variation Regularization USGS ? United States Geologic Survey WGS 84 ? World Geodetic System 1984 x Chapter 1: Introduction 1.1 Components of the Earthquake Cycle There are three components to the earthquake cycle: interseismic period, coseismic period, and postseismic period. During the interseismic period, strain varies over time and space. Spanning from years to centuries, strain accumulates on portions of the fault that are locked. This strain will continue to accumulate until a critical value is reached and the frictional strength along the fault plane fails. The seconds to minutes when there is movement along the fault plane is known as the coseismic period. The rupture zones containing the motion located within the Earth?s crust and/or mantle generate seismic waves in all directions. This rupture event enables the stress along the fault plane to readjust to a stable distribution. Following an earthquake event, there is postseismic activity where additional deformation resulting from relaxation can occur (Figure 1) (Kanamori, 1994). For example, the Mw 6.4 Hualien earthquake located in northeast Taiwan, which will be of interest to this study, ruptured on February 6th, 2018. A maximum slip of 1.6 m occurred during the coseismic period. Within approximately seven months of the mainshock, an additional 0.4 m ? 0.6 m of postseismic afterslip occurred as the fault relaxed (Zhao et al., 2020). 1 Figure 1. Plot displaying generalized strain accumulation on a fault as a function of time. 1.2 Geologic Overview of Taiwan Taiwan is geologically complicated as it is located at the intersection of two collision regimes (Figure 2). (i) There is the subduction of the Eurasian Plate beneath the Philippine Sea Plate due to part of the South China Sea lithosphere being attached to the Eurasian Plate, and (ii) there is the northward subduction of the Philippine Sea Plate beneath the Eurasian Plate. It is the latter subduction event (~5 Ma) that initiated the Taiwan Orogeny as the Luzon Arc, associated with the Philippine Sea Plate, introduced oblique collision with the Chinese Collision Margin, an accretionary prism, associated with the Eurasian Plate (Huang et al., 2017; Hsu et al., 2016). As the collision progresses, the subduction of the Eurasian Plate beneath the Philippine Sea Plate is captured between the Manila Trench and Luzon Arc, and the subduction of the Philippine Sea Plate beneath the Eurasian Plate is captured by the 2 Ryukyu Trench and Ryukyu Arc (Figure 2). As a result, in respect to the Penghu Islands off the west coast of Taiwan, the overall pattern of motion in Taiwan is fan-shaped with decreasing velocity from east to west. Figure 2. Schematic block-diagram illustrating the collisional kinematics of Taiwan. The asterisk represents the Longitudinal Valley Fault that acts as a boundary between the Eurasian Plate and the Philippine Sea Plate (from Angelier et al., 2000, after Angelier et al., 1986). As a result of the complex tectonic environment and fast plate convergence rate, Taiwan experiences increased seismic activity. The Gutenberg-Richter Relationship for all of Taiwan, which describes the relationship between time interval occurrence and earthquake magnitude, has a slope (i.e., b-value) lower than the global average of 3 1.0 (Wang et al., 2015). Subsequently, higher magnitude earthquakes are more common (Figure 3) in Taiwan compared to locations with an average b-value. Figure 3. U.S. Geologic Survey Earthquake Catalog detailing earthquakes Mw 4.5+ from July 1st, 2020 to 2021 in the east coast of the United States and Taiwan. Red lines indicate plate boundaries (USGS, 2021). Furthermore, the collisional tectonics of Taiwan have given rise to different geologic provinces: the Coastal Plain, the Western Foothills, the Hsuehshan Range, the Central Range, the Longitudinal Valley, and the Coastal Range (Ho, 1986) (Figure 4). The Coastal Plain is composed of alluvial deposits, and the Western Foothills and Hsuehshan Range are composed of marine sediments. The Central Range is composed of Tertiary metamorphism that increases in grade from west to east due to the collision orientation of the Luzon Arc, and the Coastal Range is composed of andesitic volcanic rock from the Luzon Arc. Stitching the Eurasian Plate and Philippine Sea Plate together, the Longitudinal Valley consists of Quaternary clastic fluvial sediments (Hsu 4 et al., 2009). An additional area of interest includes the Pingtung Plain in southwest Taiwan which is made of weak muds from the South China Sea rift (Hu et al., 2007). Figure 4. Tectonic and Geologic Layout of Taiwan. The red lines indicate locations of subduction, and the dashed red lines indicate potential locations of subduction. The white lines indicate geologic province boundaries: CP ? Coastal Plain, WF ? Western Foothills, HR ? Hsuehshan Range, CR ? Central Range, LV ? Longitudinal Valley, CoR ? Coastal Range. Additional regions of interest include: PP ? Pingtung Plain, TB ? Taipei Basin, IB ? Ilan Basin. 5 1.3 Purpose of Study Various geodetic techniques have been used to study the evolving crustal deformation of Taiwan. Using Global Navigation Satellite Systems (GNSS) measurements, Yu et al. (1997) and Lin et al. (2010) observe more than 80 mm/yr of collision between east Taiwan and the Penghu Islands ~30 km west of Taiwan. The GNSS data shows a clockwise rotation in northeast Taiwan that indicates post-rift opening of the Okinawa Trough as well as lateral extrusion (i.e., tectonic escape) in southwest Taiwan (Lacombe et al., 2001). However, the GNSS measurements from those studies only reveals first order surface deformation as they cannot provide high spatial sampling. In particular, they are unable to show creeping/locked fault traces that are on a smaller spatial scale, and it is difficult to determine interseismic fault locking depth without high spatial resolution crustal deformation data (Elliott et al., 2016). Alternatively, interferometric synthetic aperture radar (InSAR) provides high spatial resolution measurements of surface deformation (e.g., fault creeping/locking in west Taiwan) at a lower accuracy with cm-scale precision (B?rgmann et al., 2000; Elliott et al., 2016). Recent work (e.g., Huang and Evans, 2019; Tong et al., 2013; Weiss et al., 2020) attempts to combine both GNSS and InSAR to achieve a high spatial resolution with relatively high accuracy. This method is done by utilizing GNSS for the long- wavelength spatial deformation and InSAR for short-wavelength features. Huang and Evans (2019) estimate crustal deformation in southwest Taiwan using six years of InSAR and GNSS data. They characterize fault slip and locking depth 6 for the major fault system in southwest Taiwan. Weiss et al. (2020) uses this technique to estimate surface strain rate for the North Anatolian Fault, Turkey. In this study, I employ GNSS and InSAR to characterize interseismic crustal deformation from 2016 to 2021 and to utilize block modeling to quantify slip budget along faults to forecast seismic hazard potential in Taiwan. 1.4 Hypotheses & Significance of Study Along an active plate boundary, continuous plate convergence will cause energy stored in locked faults to eventually rupture and initiate seismic events. The quantification of interseismic fault slip budget can help evaluate seismic risks. For example, geodetic techniques can be used to estimate crustal strain accumulation in the deforming brittle crust to evaluate seismic potential along faults. Evaluating the trends and behaviors of interseismic slip budget will enable proper seismic hazard forecasts models to be developed. 7 Chapter 2: Data and Methods 2.1 Synthetic Aperture Radar and GNSS Data Collection Data were obtained from the European Space Agency?s (ESA) Sentinel-1 mission for the Copernicus initiative. This mission collects synthetic aperture radar (SAR) acquisitions with a wavelength of 56.7 mm (C-band) in the form of single-look complex (SLC) data (i.e., no image compression). The SLC products containing the SAR images were downloaded from the Alaska Satellite Facility (ASF), University of Alaska database through the National Aeronautics and Space Administration (NASA) Earth Data service. In this study, I used Sentinel-1 SAR acquisitions from ascending track 69 and descending track 105 from November 2016 to July 2021. The precise orbital data, which details the trajectories of the Sentinel-1 satellites' flight paths, were downloaded from ESA Science Hub. The digital elevation models (DEM) were downloaded from the NASA Jet Propulsion Lab?s (JPL) Shuttle Radar Topography Mission (SRTM) with 30 m resolution and 3-arc second (Farr et al., 2007). The SRTM data are currently stored in the United States Geologic Survey (USGS) Measures project. The DEM data were used to remove elevation contributions to phase in the InSAR images. Combining the orbital and DEM data enabled geocoding. That is, the differing Sentinel-1 satellites? flight paths were aligned to the DEM that utilized coordinates from spatial reference system World Geodetic System 1984 (WGS 84). The weather model used in the troposphere noise correction was 8 downloaded from European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 weather model products. The continuous GNSS time-series data were processed and maintained by the Central Geologic Survey (CGS), the Central Weather Bureau, the Ministry of Interior, Taiwan, the Global Positioning System (GPS) Laboratory at the Institute of Earth Science, and Academia Sinica, Taiwan. These data were downloaded from Academia Sinica, Taiwan. 2.2 Geodetic Technique: Interferometric Synthetic Aperture Radar Radar is a detection instrument used to map an objects distance or speed. SAR, which utilizes radar technology, collects the backscatter of electromagnetic waves to produce high resolution radar images. To perform a single acquisition, a radar positioned on the side of an object in orbit emits a sequence of radar pulses in the microwave-range at a specific look angle, the pulses are backscattered by the ground surface, and the amplitude and phase of the backscatter is recorded (B?rgmann et al., 2000) (Figure 5). Acquiring a stretch or scan of backscattered energy acts as a synthetic antenna of considerable size, which enables the image to be of greater resolution. Acquisitions can be gathered in either the ascending direction (south to north) or the descending direction (north to south) (Figure 6). This is known as the heading direction. Multiple acquisitions of the surface are acquired by repeated flybys over a period of time. If the ground experiences movement during this time, the backscattered signals comprised of phase and amplitude values can be used to quantify the surface displacement. The Wave Equation can be represented as (Equation 1): 9 ?(?, ?) = ?sin(?? ? ?? + ?) [1] where Y is vertical velocity of the wave, A is amplitude, k is wave number, x is horizontal position, ? is angular frequency, t is time, and ? is phase. Figure 5. Phase and amplitude of Environmental Satellite (ENVISAT) SAR acquisition dated 01-04-2007 (month-day-year) for southwest Taiwan. This acquisition is SLC where pixel resolution is 5 m east-west and 20 m north-south. As a result of the pixel resolution, there is distortion. 10 Figure 6. Ascending and descending flight paths detailing look angle and heading direction. InSAR uses SAR images to produce interferograms that reveal ground elevation information. The amplitude data from the SAR images are used to align the series of images taken over a period. The phases of these images are used to determine ground movement. By subtracting ?2 - ?1 for each pixel, the ground displacement moving towards or away from the satellite in the line-of-sight (LOS) direction can be calculated (Figure 7) (B?rgmann et al., 2000). If the whole area moves by the same amount, there will be no phase shift; therefore, InSAR only measures relative deformation. 11 Figure 7. Concept of phase difference produced from ground motion captured by SAR acquisitions (based on Funning, 2021). 2.2.1 Limitations The reliability of InSAR is limited by atmospheric structure (i.e., troposphere and ionosphere) and decorrelation between an image pair. Water vapor in the troposphere can refract radar waves which distorts the returned phase and subsequently skews the suggested amount of deformation. Longer wavelengths have a lower refraction index and are less susceptible to weather (Massonnet and Feigl, 1998). Due to the ionization properties of the ionosphere, this layer of the atmosphere is also able to distort radio waves (Massonnet et al., 1993). Decorrelation refers to the change in ground surface features (i.e., snow, vegetation, flooding) between image acquisitions. 12 This ground surface feature change is not related to movement but can cause a phase change or signal loss (B?rgmann et al., 2000). Additionally, InSAR only measures ground deformation in LOS, and it is not sensitive to north-south motions. 2.3 Data Analysis Process for InSAR 2.3.1 GNSS Time-Series Adjustment and Mean Velocity Estimation The GNSS time-series were manually adjusted using a constant offset to remove transient events (i.e., GNSS station malfunctions) at each GNSS station that were not associated with a known earthquake (Figure 8). 13 Figure 8. Example of GNSS time-series manual adjustment showing the original time-series of GNSS station C001 in Taiwan and the adjusted time-series of GNSS station C001. The adjustment is located at 2018.4 in the east-west motion time-series. I then used the following mathematical model to model the adjusted GNSS time-series for the east, west, and vertical component of each station (Equation 2). The model was solved as an inverse problem (see Section 2.3.4). 14 ????(?, ?, ?, ??) = ?1(?, ?, ?) + ?2(?, ?, ?)?? + ?3(?, ?, ?) sin(2???) + ?4(?, ?, ?) cos(2???) + ?5(?,?, ?) sin(4???) + ?6(?, ?, ?) cos(4???) [2] where m1 is a constant representing an adjustment in distance for the time-series, m2 is the linear trend of the station throughout time, m3 and m4 are the annual seasonality of the station throughout time, and m5 and m6 are the semi-annual seasonality of the station throughout time. 2.3.2 InSAR Scientific Computing Environment Stack Processor Module The InSAR products were processed and produced using InSAR Scientific Computing Environment (ISCE) software developed at NASA JPL Caltech (Rosen et al., 2012). Stack Processor is a module of the ISCE software package that enables SAR images to be combined to generate InSAR images (Figure 9) (Fattahi et al., 2017; Rosen et al., 2012) while automatically applying phase unwrapping using Snaphu (Chen and Zebker, 2002). First, I used the Stack Processor module to define various details of interest: SAR acquisitions, orbital data, DEM data, bounding box, auxiliary data (i.e., Sentinel- 1 instrument parameters), the number of adjacent SAR images to be processed (network of 3), and the start and end dates. Then, the SAR reference and target image pairs were coregistered pixel-to-pixel by amplitude to align the pixel locations in each image and determine the phase difference. Next, the images were downsampled as they were taken from SLC to a multi-look complex (MLC) (i.e., coarser pixel resolution) by compression, a phase filter was applied to enhance the phase signal, and phase 15 associated with the Earth?s curvature was taken into consideration as a simulated phase signal known as ?Earth flatten? was applied to remove additional fringes. After the module processed the SAR images, it unwrapped the resulting interferograms. Phase unwrapping unbound the phase values between -? and ? to determine the true surface deformation values. Figure 9. By differencing the phase values from different SAR acquisitions, a wrapped interferogram is produced that details ground motion using ?fringes.? Each fringe represents a phase change of 2? which is equivalent to 3 cm (C-band). Note: The SAR acquisitions are dated 01-04-2007 (month-date-year) and 03-13-2007 and were collected by ENVISAT. The phase values from t1 and t2 were in the form of SLC where pixel resolution was 5 m east-west and 20 m north-south as resolution is a function antenna size and 16 wavelength. When the phase values were differenced, they produced an image of ground motion in Taiwan. The interferogram must then be geocoded to properly organize the pixels that were previously aligned in capture order (i.e., descending captures easterly pixels first). At this point, the data were not yet fully geocoded and remained in radar coordinates. Furthermore, the Stack Processor module processed the ascending track data and the descending track data separately. 2.3.3 Miami InSAR Time-Series Software in Python (MintPy) I used MintPy (Yunjun et al., 2019) to generate InSAR time-series. MintPy applied tropospheric and ionospheric noise corrections, generated time-series, and determined LOS InSAR mean velocities for both the ascending and descending tracks. This software used a small baseline subsets approach to find the best-fitting time-series for the given interferograms while minimizing the implied velocities (Bernardino et al., 2002). To generate a displacement time-series and mean velocities from the phase- unwrapped interferograms for both the ascending and descending LOS InSAR data, the following steps were taken: (i) The network of interferograms was modified to exclude interferograms with average spatial coherences (i.e., pixel's phase similarity to surrounding pixels) less than predefined threshold value of 0.3. This was called a coherence-based modification and was optional. (ii) A reference pixel among the pixels with the highest average spatial coherences was randomly selected. This reference pixel was assumed to have minimal deformation and impact from phase delays. If specified, a reference region could also be defined. (iii) The stack of interferograms was examined and underwent a correction for unwrapping errors by enabling bridging (uses spatial 17 resolution of phase to detect and correct unwrapping error) and phase closure (uses pixel phase relationship of same pixel over time used to adjust phase cycle), and a water mask was applied to remove decorrelated bodies of water. This step was optional. (iv) The time-series was solved for by minimizing the interferometric phase residual between each pixel throughout time. During this step, a temporal coherence (i.e., stability of pixel throughout time) mask and a spatial coherence mask can be applied to remove pixels below the predefined threshold value. After masking, each pixel may have a different number of interferograms. To address this potential problem, a minimum number of interferograms required for inclusion in the final product was set. (v) A tropospheric delay correction was applied using Global Atmospheric Models (GAMs) data. I did not de-ramp the phase data as the deformation was of larger scale and the appearance of a ramp would likely be a true phase signal. (vi) A topographic phase residual correction caused by DEM error was applied. (vii) The remaining phase component that could neither be corrected nor determined as ground deformation was removed as noise. (viii) The linear slope that best fit the displacement time-series for each pixel was determined. (ix) For a post-processing step, the data were geocoded from radar coordinates. As a result, the displacement time-series of each pixel was determined, and the mean velocity of each pixel throughout time was estimated. MintPy was processed separately for the ascending track data and the descending track data. Furthermore, to avoid phase unwrapping errors due to steep mountain ranges and dense vegetation between west and east Taiwan that lacked stable persistent scatterers, MintPy processing was completed separately for both west and east Taiwan. This enabled west and east Taiwan to be examined relative to one another. The reference 18 point was set to both respective sides when processing MintPy to best capture pixel data. 2.3.4 Time-Series Re-Model & Mean Velocity with GNSS Correction For a more in-depth data analysis, a module developed by Huang and Evans (2019) was used to (i) re-model the time-series for the ascending and descending LOS InSAR velocities using a polynomial and (ii) re-estimate the ascending and descending LOS InSAR mean velocities to include a GNSS correction. To exclude low-quality pixels from the ascending track and descending track LOS InSAR data, the temporal and spatial coherence thresholds of the mean ascending and descending LOS InSAR velocities were determined through a trial-and-error method of pixel visibility. Temporal coherence refers to the stability of a pixel throughout time ? how similar the pixel phase is between acquisitions. The more stable a pixel, the higher the temporal coherence. Deformation reduces the temporal coherence within reason. Spatial coherence refers to the consistency of a pixel?s phase to surrounding pixels. Sharp phase changes between neighboring pixels may indicate an error. I used a temporal coherence threshold of 0.3 for west and east Taiwan, and a spatial coherence threshold of 0.4 for west Taiwan and east Taiwan. The average of the temporal and spatial coherence values above the defined thresholds was determined to be the final coherence value. This final value must be above the final predefined mask value of 0.35. Applying this mask excluded all pixels with values lower than the final coherence value. After the ascending and descending LOS InSAR velocities were masked, (i) the elevation of each pixel for each interferogram was defined using the DEM, (ii) the look 19 angles and heading directions for each interferogram were defined, (iii) the reference image was set to the first acquisition date for both the ascending and descending track interferograms, and (iv) the latitude and longitude data for the bounding box were linked to the module to geolocate each pixel in each interferogram. The west Taiwan ascending and descending LOS InSAR velocities were assigned a local reference region in the west, and the east Taiwan ascending and descending LOS InSAR velocities were assigned a local reference region in the east. The designated reference regions were at an area without known faults and minimal seasonal surface movement due to hydrologic cycles and human induced land subsidence. Subsequently, the reference regions were considered as stable regions with zero movement. Once the reference region was defined, the mean velocity of each pixel location throughout time was determined. To fit the time-series of each pixel for ascending and descending LOS InSAR velocities for both west and east Taiwan, I generated a mathematical model with a linear velocity term, annual periodic terms, and semi-annual periodic terms (Equation 3). The terms utilized match the general pattern of interseismic deformation anticipated in Taiwan. For example, the linear velocity term accounted for the overall mean velocity of the pixel, the annual periodic terms took into consideration the wet and dry seasons? influence on motion, and the semi-annual periodic terms considered sub-tropical precipitation events (e.g., monsoon vs. typhoon events). These terms for a pixel at location (x,y) were represented as: 20 ???(?, ?, ??) = ?1(?, ?) + ?2(?, ?)?? + ?3(?, ?)sin(2???) + ?4(?, ?) cos(2???) + ?5(?, ?) sin(4???) + ? ?? ?6(?, ?) cos(4???) + ?7(?, ?)?(?EQ) + ?8(?, ?) ?(?EQ) ln [1+ ( ? EQ)] [3] ? where m1 is a constant representing a constant adjustment for the time-series, m2 is the linear trend of the pixel throughout time, m3 and m4 are the annual seasonality of the pixel throughout time, m5 and m6 are the semiannual seasonality of the pixel throughout time, m7 is the Hualien earthquake (i.e., a notable earthquake within the timeframe of interest) displacement coefficient, H(tEQ) is the step function to remove the 2018 Hualien Earthquake coseismic event, m8 represents the postseismic period with a relaxation time of ? = 121 days (See Section 2.5.1 for more details). Now, I assume G matrix represents the mathematical model described in Equation 3, ? = ?????? [4] where ? is the data vector (LOS (x, y, t)), G is the mathematical model that relates the model parameters to the data (right hand side of Equation 3), and ????? is the model vector (m1, m2, m3, m4, m5, m6, m7, m8). This mathematical model enables the fitting of a time- series at each pixel location throughout time for both the ascending and descending track data for both west and east Taiwan. 21 I used a least squares inversion to solve for the mathematical model and estimate the coefficient of each term. This solving approach minimizes the sum of squares of the residuals (Equation 5). ????? = (???)????? [5] where ????? is the model vector, G is the mathematical model, GT is the transpose of the mathematical model, and ? is the data vector. Then, using the pre-processed time-series of each GNSS station, I applied a GNSS-correction to the mean velocities of each pixel derived from the calculated time- series. This correction applied the accuracy of GNSS to the high spatial resolution of InSAR. The GNSS-correction only considered the same time period as the InSAR data and was comparable to the LOS InSAR data as the displacements were projected onto the satellite look angle and heading direction. To apply the GNSS correction to the ascending and descending LOS InSAR velocities, a ramp model that best fit the InSAR and GNSS data velocity differences was constructed. The coefficients of the ramp model were solved for by inversion (Equation 5) with the velocity residuals as data. Removing the ramp from the uncorrected ascending and descending LOS InSAR velocities produced the GNSS- corrected ascending and descending LOS InSAR velocities. Additionally, the ascending and descending ramps for west and east Taiwan could be applied to the time- series for GNSS correction inclusion. 22 2.3.5 Merge Ascending and Descending from East and West Taiwan To begin merging the ascending and descending LOS InSAR velocities from west and east Taiwan, the low coherence (or low quality) pixels (e.g., pixels capturing water) from each dataset were masked out and set to 0. High coherence pixels were set to 1, and pixels that were high coherence in both datasets were set to 2. During the GNSS correction, west and east Taiwan were assigned the same reference location; therefore, here, they were merged without searching for a common reference region. For accuracy purposes, if there were overlapping real valued pixels from both datasets, the pixels from the east Taiwan dataset were kept while the pixels from the west Taiwan dataset were set to 0. This merged masking process was done for both the ascending and descending LOS InSAR velocities. Once the masks were created, the values of the real-valued pixels were utilized and datasets with ascending LOS InSAR velocities and descending LOS InSAR velocities for all of Taiwan were created. 2.3.6 Convert GNSS-Corrected LOS InSAR & GNSS Velocities GNSS-corrected LOS InSAR velocities and GNSS velocities were utilized to estimate 3-D deformation: east-west, north-south, and vertical motion. First, the GNSS velocities were interpolated to InSAR pixels using 2-D cubic interpolation in Matlab. The mesh size matched the pixel location and size of that from InSAR geocoded to the DEM. The inclusion of these velocities enabled a more accurate 3-D velocity field of Taiwan to be constructed as, for example, InSAR has poor sensitivity to north-south velocities and GNSS velocities are less sensitive to atmospheric phase delays. The GNSS-corrected LOS InSAR velocities and interpolated GNSS velocities were converted to 3-D deformation by relating the heading direction and look angle of 23 the satellites to the velocity data through an inverse problem in the form of Equation 4. The final velocity product was as follows (Equation 6): ???? cos??sin?? sin??sin? ? ?cos?? ??? ? cos??sin?? sin??sin?? ?cos? ? ? ? ????? = 1 0 0 [?? ] [6] ????? 0 1 0 ?? [????? ] [ 0 0 1 ] where data vector ? contains: LOSA,D the LOS velocity for the ascending and descending tracks and GNSSE,N,Z the interpolated GNSS velocities in east, north, and vertical, respectively. Matrix G contains: ?A and ?D the satellite heading direction for the ascending and descending tracks and ?A and ?D the satellite look angle of the ascending and descending tracks, respectively. This matrix relates the InSAR and GNSS velocities to their 3-D components. Model vector ????? contains the 3-D velocity outputs UE,N,Z. The linear inverse problem was solved for using a least squares inversion to minimize the sum of squares of residuals and to determine the best fit model (Equation 5). Additionally, in order to weigh each component of the output 3-D velocity dataset based on misfit, I incorporated a weighting matrix W into the least squares inversion (Equation 7). ????? = (????)??????? [7] 24 where matrix W (Equation 8) is used to weigh the data and is solved for during Section 2.5.1. ? ?2 ? 0 0 0 0 0 ? ?2 ? 0 0 0 ? = 0 0 ? ?2 ???? 0 0 ? [8] 0 0 0 ? ?2 ???? 0 ? [ 0 0 0 0 ? ?2????? ] where matrix W is the weighting matrix and ? for ascending, descending, and GNSS represents the misfit values produced from the inversion of Equation 3. Using the resulting 3-D velocity outputs, the Final InSAR and GNSS (FIG) dataset was created. This dataset includes: the weighted mean GNSS-corrected InSAR / interpolated GNSS velocity values with the associated uncertainties and the GNSS velocity values with the associated uncertainties. The GNSS velocity values were appended to the dataset for additional data point inclusion. Uncertainties are solved during Section 2.5.2. Additionally, a Reduced FIG dataset was created, which contained the values within the FIG dataset downsampled to every 10 pixels in both the x- and y-direction. The FIG and Reduced FIG datasets contain pixels that are 50 m x 50 m and 500 m x 500 m, respectively. 2.3.7 Mean Velocity Gradient To better visualize locations of increased strain accumulation with source unspecified resulting in a sharp change in velocity, I calculated the east-west, north- 25 south, and vertical mean velocity gradients of the FIG dataset in a raster image in the WGS 84 system. Additionally, prior to mean velocity gradient calculation, the velocity data were smoothed to distinguish noise signal from tectonic signal (see Section 2.5.3). ? 2 ? 2 ??,?,?(?, ?) = ?[ ??,?,?(?, ?)] + [ ??,?,?(?, ?)] [9] ?? ?? Where ??,?,?(?, ?) is east-west, north-south, or vertical mean velocity gradient, ? ??,?,?(?, ?) is the velocity difference (mm/yr) between neighboring pixels in the x-?? direction for east-west velocity, north-south velocity, and vertical velocity, and ? ??,?,?(?, ?) is the velocity difference (mm/yr) between neighboring pixels in the y-?? direction for east-west velocity, north-south velocity, and vertical velocity. Each pixel is 50 m x 50m, which is equivalent to 5000 cm x 5000 cm. 2.4 Deformation Rate Analysis 2.4.1 Dilatation, Maximum Shear, and 2nd Invariant The deformation tensor defines position change within a body due to external forces (Figure 10). Using the Reduced FIG dataset, I determined the deformation rate tensor every 1 km where the InSAR samples were located every 500 m. The calculated deformation rate tensors considered the nearest 30, 144, and 420 pixels (to be continued in Section 2.5.3). From this, I calculated dilatation (unit: yr-1), maximum shear (unit: yr-1), and 2nd invariant (unit: yr-1). Note: Dilatation, maximum shear, and 2nd invariant refer to their rate per year. The purpose of the deformation rate analysis was to quantify 26 2-D deformation fields across Taiwan. This analysis assumed that deformation fields were subject to variations in stress rather than strength (Fagereng and Biggs, 2019). The deformation rate tensor is defined as: ????? ????? ??? = [ ] ????? ????? 1 1 ????? ( ????? + ????? ) 0 ( ????? ? ??? ) = [ 2 ?? 1 ] + [ 2 1 ] ( ????? + ?????) ????? ? ( ??? ? ??? ) 02 2 ?? ?? ??? 1 ??? ?? 1 ?? ?? ( + ? ) 0 ( ? ? ? ) ?? 2 ?? ?? 2 ?? ?? = [ ] + [ ] [10] 1 ?? ( ? ?? ?? + ? ) ? 1 ?? ?? ? ( ? ? ? ) 0 2 ?? ?? ?? 2 ?? ?? where the deformation rate tensor, ???, is the sum of the strain rate (irrotational) matrix ??? ?? 1 ?? ?? 1and rotational rate matrix. = ??? ? ? ? ?? ?? , = ??? , and ( + ) = ( ??? + ?? ?? 2 ?? ?? 2 ?? ????? ). The off-diagonal terms in the rotational matrix are equal in quantity but change in sign. 27 Figure 10. Schematic of maximum (?1) and minimum (?2) principal strains and the corresponding strain tensor components (Dxx, Dyy, Dxy) influencing a square. This schematic does not consider rotation. Using components of the deformation rate tensor, I solved for dilatation, the overall change in volume due to deformation. Dilatation is the sum of principal strains, which are the eigenvalues of a strain rate tensor (Equation 11, 12). | A - ? ? I | = 0 [11] where A is the strain tensor, ? are the eigenvalues, and I is the identity matrix. The | | sign represents the determinant operation. ? = ?1 + ?2 [12] 28 where ? is dilatation and ?1 and ?2 are the maximum and minimum principal strains (i.e., eigenvalues). Then, maximum shear was solved to determine the factor in which deformation occurred in a specific direction (Equation 13). In this case, maximum shear (i.e., change in shape/angle) corresponds to the greatest shear at 45? to the principal strains. ?1? ?? 2??? = [13] 2 where ???? is maximum shear and ?1 and ?2 represent the maximum and minimum principal strains (i.e., eigenvalues). Invariants of the deformation rate tensor are properties that do not change under coordinate rotation. The 2nd invariant of strain rate determines the total strain rate accumulation of the area of interest, which highlights localities with increased seismic risk (Equation 14) (Pagani et al., 2021). It acts as a combination of both the dilatation (contraction and extension) and maximum shear stress. ? 2 2 1 2 ?2 = ??? + ??? + 2 [ (??? + ???)] [14] 2 29 where ? is 2nd2 invariant of strain rate and Dxx, Dyy, Dxy, Dyx are components of the symmetric strain rate tensor. ??? and ??? cannot be assumed to be of the same value as rotation, which does not address shape change and is not taken into consideration. 2.5 Error Analysis I incorporated an error analysis into the InSAR and GNSS velocity solutions by determining root mean square (RMS) misfit. As previously mentioned, the calculated misfit was used to define the weighting matrix, W (Equation 8), to properly weigh between the GNSS-corrected InSAR and the interpolated GNSS velocities. The uncertainties, inferred from misfit, produced by taking the LOS and GNSS velocities to east-west, north-south, and vertical were utilized to confirm consistent transformation and were appended to the FIG dataset for later usage. Furthermore, to distinguish tectonic signal from noise in the deformation rate analysis, I quantified distance-correlated noise structure using a semi-variogram and covariogram model for a region without known surface deformation. 2.5.1 Solving for RMS Misfit & ? When constructing the mathematical model (Equation 3) that best fits the velocity data, I calculated the RMS misfit to detail the misfit between the model and the observed velocity values (Equation 15). Specifically, the RMS misfit was calculated for the ascending and descending LOS InSAR data (every pixel at all times) and the east-west, north-south, and vertical GNSS data (every station at all times). Since west and east Taiwan were processed separately, the west and east Taiwan InSAR misfits were calculated separately and then merged. 30 1 2 ?(?, ?) = ? ???=1(??(?, ?) ? ?????? ?(?, ?)) ??? ? = 1,??; [15] where ?(x,y) is the RMS misfit of pixel (x,y), ?i is observed velocity data, and ?????? is best fit model velocity data, i is the index of an acquisition, and N is the total number of acquisitions. Additionally, to determine the relaxation time (?) for removing the postseismic contribution from the 2018 Hualien Earthquake in the mathematical model (Equation 3), I calculated the RMS misfit of all the pixels with a given ? value between 1 and 600 days in a 20-day step size: 1 ?(??) = ? ? ? ?=1(??(? 2 ?) ? ??????(??)) ??? ? = 1,??; j = 1,?600 (days); [16] ? where ?(??)is the RMS misfit of relaxation time (?) in j days, ?i is observed velocity data, and ?????? is best fit model velocity data, i is the index of an acquisition, and N is the total number of acquisitions. The relaxation time that produced the least amount of misfit was used in Equation 3. 2.5.2 Solving for Uncertainty The uncertainty of every pixel was solved by evaluating the misfit of the inversion that transformed the LOS InSAR and GNSS velocities to 3-D velocities. The 31 uncertainty values at each time were inferred from the velocity misfit values using a linear inverse problem in the ? = ?????? form (Equation 17): ??(?, ?) cos??(?, ?)sin??(?, ?) sin??(?, ?)sin??(?, ?) ?cos??(?, ?) ??(?, ?) ??? (?, ?) cos??(?, ?)sin??(?, ?) sin??(?, ?)sin??(?, ?) ?cos??(?, ?) ? ? ?????(?, ?) = 1 0 0 [????(?, ?)][17] ? ?????(?, ?) 0 1 0 ????(?, ?) [ ? ?????(?, ?) ] [ 0 0 1 ] where data vector ? contains: ?A,D (x,y) the misfit for the ascending and descending track pixels and ? GNSSE,N,Z (x,y) the interpolated GNSS misfit for east, north, and vertical, respectively. Matrix G contains: ?A and ?D the satellite heading direction for the ascending and descending tracks and ?A and ?D the satellite look-angle of the ascending and descending tracks, respectively. This matrix relates the GNSS-corrected InSAR / interpolated GNSS velocities to their 3-D components. Model vector ????? contains the 3-D velocity uncertainty estimates, UncE,N,Z, for east, north, and vertical components, respectively. The identity matrix is used to bring the GNSS misfit values through the inverse problem with no transformation as they are already in 3-D form. 2.5.3 Noise Structure Contributions in the Deformation Rate Analysis When calculating the deformation rate tensor in Section 2.4, various sizes of smoothing windows were utilized. To determine which level of smoothing best eliminated the noise structure contribution, I calculated a semi-variogram model and covariogram model of a non-deforming region for error estimation (Sudhaus and Jonsson, 2011). The semi-variogram was modeled from pixel variance with distance in the x- and y-direction (Equation 18) and suggested the use of an exponential equation 32 to model the covariogram (Equation 19). The covariogram estimated pixel spatial correlation with distance (i.e., covariance) (Sudhaus and Jonsson, 2009). The semi-variogram was defined as (Equation 18): ? ?(?) = ?2(1 ? ???) [18] where S(r) is the modeled semi-variogram between two pixels, ?2 is the variance, r is the distance between the two pixels, and ? is the characteristic wavelength of the transect. The modeled covariogram, produced from an exponential mathematical model, was defined as (Equation 19): ? ?(?) = ?2??? [19] where C(r) is the covariance between two pixels, ?2 is the variance, r is the distance between the two pixels, and ? is the characteristic wavelength of the transect. The modeled semi-variogram was solved for with an exponential, spherical, and gaussian mathematical model. The exponential model fit best and was subsequently utilized to model the covariogram. The unknown ?2 and ? in the covariogram model (Equation 19) were solved using an inverse problem with the FIG dataset velocities for a non-deforming region as data constraints. The covariogram model acted to estimate 33 the characteristic wavelength of correlation to quantify the assumption that variables closer in distance tend to be more similar (Watson et al., 2022; Hussain et al., 2016). Therefore, a smoothing window size that is smaller than the characteristic distance ? at which pixels were spatially correlated may display noise signals and not accurately capture the tectonic deformation influencing the region. Given that the deformation tensor was calculated every 1 km and each pixel is 500 m x 500 m in the Reduced FIG dataset, utilizing the nearest 30, 144, and 420 pixels resulted in a 1.5, 3.4, and 5.8 km radius of values being incorporated into the tensor, respectively (Table 1). Table 1. Conversion of pixel count to physical distance. Number of Nearest Pixels (area) Radius of Pixels Radius in km (A = ?r2) 30 3.09 1.5 144 6.77 3.4 420 11.56 5.8 2.6 Geodetic Technique: Block Modeling Block modeling enables a region of interest to be divided into blocks that are completely bound by segments (i.e., faults). This allows a plate-boundary scale area to be examined at the regional-scale, microplate level. In this study, block construction was based upon many lines of evidence: mean velocity gradients, active fault traces, geologic maps, and block geometries from published work (e.g., Ching et al., 2011; Chang et al., 2016). Block modeling uses inverse theory to apply observed velocities to a given construction to determine block rotation and translation and the difference of motion 34 between blocks to be interpreted as fault slip (Meade and Loveless, 2009). These blocks are subject to a variety of behaviors. Rigid translation and rotation occur when the block moves or rotates without causing internal deformation, whereas elastic behavior occurs when part of a fault segment slips and then causes internal elastic deformation within the block. The rigid translation and rotation between blocks can be considered as the long-term crustal motions between plates, where slip above or below the locking depth can be considered as coseismic or interseismic deformation, respectively (Figure 11). To estimate contributions from both rigid and elastic behaviors, block modeling takes into consideration a dislocation model (back-slip) that details the coseismic period defined by various degrees of elastic behavior (Figure 11). The dislocation model can constrain fault parameters given surface deformation (observation) and fault geometry. The outcome of subtracting a dislocation model (i.e., back-slip) from a completely rigid block motion produces a more realistic representation of interseismic slip and surface crustal deformation during this period. 35 Figure 11. Removing coseismic rupture represented by a dislocation model from a block model (i.e., long-term slip) that assumes complete rigidity at depth enables interseismic creep to be estimated. Blue represents a transect cutting across the fault. Red indicates creep. I used ?Blocks? (Meade and Loveless, 2009) with a Total Variation Regularization (TVR) algorithm developed by Evans et al. (2015), to determine the reference block solution to the observed velocities. The latter algorithm is unique in that it allows microplates to be grouped with neighboring microplates to identify the most influential faults and produce the simplest block model configuration that can sufficiently explain observations with given data uncertainty (Figure 12). Figure 12. Conceptualization of block reduction where the black lines represent faults, and the colored stars represent Euler pole locations for each block. The left-side panel is separated into Blocks 1, 2, and 3 with slip along the faults bounding them. The right- side panel combines Blocks 1 and 2 and assigns them the same Euler pole and rotation rate. There is no longer slip between them and they are considered to be one larger block (from Evans et al., 2015). 36 2.6.1 Limitations There are four main limitations to block modeling. First, the model assumes interseismic deformation is a combination of rigid block movements and the subtraction of (mostly shallower) elastic dislocations. However, the crust is neither rigid nor elastic. Second, it only takes into consideration horizontal velocities due to program constraints and conceptual limitations. This may not be ideal for a convergent plate boundary because in Taiwan the vertical component of motion is also significant. Third, block modeling does not consider faults that are not incorporated into the model construction. That is, block segments can be removed under regularization, but they cannot be added. Fourth, it does not consider the contribution from postseismic deformation, which is often ~10% of coseismic deformation (e.g., Zhao et al., 2017). 2.7 Data Analysis Process for Block Model 2.7.1 Map Faults on QGIS To define known and potential faults in Taiwan, the following sources were utilized: active fault maps from the CGS of Taiwan (CGS, 2012), geologic maps from the Chinese Petroleum Corporation (CPC, 1989a-d), Chen (2016), Chang et al. (2016), Ching et al. (2011), SRTM DEM of Taiwan, the FIG dataset, and the mean velocity gradients. 2.7.2 Block Model Construction and Evaluation The active and potential faults were traced on QGIS and could be uploaded into either the Blocks (Meade and Loveless, 2009) or Blocks TVR (Evans et al., 2015) 37 program in Matlab. Using the faults determined on QGIS as a guide, the block segments and blocks were delineated in the Segment Manager module. Each block segment must be a part of a bounded block and each block must be fully bound by block segments. Utilizing fault trace and information provided by Chan et al. (2020), Shyu et al. (2020), Shyu et al. (2016), and the CGS of Taiwan (CGS, 2012), each block segment was assigned a locking depth and dip angle. Following the delineation of (i) block segments and blocks and (ii) block segment locking depth and dip, either program could be used to apply linear block theory, which uses linear algebraic functions to decompose the given surface velocities into plate rotation and intrablock strain accumulation. The observed, input surface velocities were those from the FIG dataset. During linear block theory application, Blocks minimizes the L2 norm (least squares: minimize the sum of squares of differences) of the data residuals in order to find the calculated velocities that best fit the observed velocities (Meade and Loveless, 2009). Alternatively, Blocks TVR optimization algorithm minimizes the L2 norm (least squares: minimize the sum of squares of differences) of data residuals and the L1 norm (least absolute error: minimize the sum of absolute values of residuals) of block motion variations (Evans et al., 2015). The latter approach enables the grouping of blocks to reduce the number of microplates fitting the data. To produce a collection of potential models, Blocks TVR was used as it minimized both the residual velocities (L2 norm) and the variation in block rotation vectors (L1 norm) to various extents. To process the L1 norm, a Matlab toolbox called CVX: Disciplined Convex Programming was required. This is a modeling program that 38 solves convex optimization problems. The following equation was used to minimize and balance both the L2 and L1 norms (Equation 20): 1 min??2(????? ? ?)? + ???????? [20] 2 1 where ? controls the strength of regularization (similar to a smoothing parameter), W is the diagonal weighting matrix that takes into consideration velocity uncertainties, G is the mathematical model, ? is the data vector describing the interseismic east-west and north-south velocities, D is the linear differential operator that applies smoothing to the block model, and ???? is the vector of block rotation. The L2 component (Equation 20) produced calculated velocities that best fit the block construction with minimized data residuals. The calculated velocities were influenced by the diagonal weighting matrix W that was built from velocity misfit values (see Section 2.3.6). The L1 component (Equation 20) utilized the parameter ? to control the amount of regularization (i.e., control block motion variation). For example, increasing ? would add weight to the smoothing operator D and decrease the number of distinct microplates utilized to fit the given surface velocity data. However, both the L2 norm and L1 norm cannot be fully minimized simultaneously. For example, if the L2 norm was minimized to contain no data residuals, then the L1 norm would be maximized with maximum blocks. Alternatively, if the L1 norm were to be minimized to a singular block, then the L2 norm would be maximized with high data residuals due to the low number of degrees of freedom. 39 Therefore, Blocks TVR serves to display a range of balances between the L2 and L1 norm by processing multiple iterations of varying ? values (Evans et al., 2015). The collection of outputs was evaluated by examining the mean residual velocities (MRV) to the number of blocks and by applying a reduced chi-squared test. Both evaluations followed the same pattern but supplied different information. MRV quantified how the output block model velocity residuals compared to the input geodetic velocity uncertainties given block regularization. MRV did not scale for data uncertainty in the models. The reduced chi-squared test determined how the residual velocities of the observed and modeled velocities statistically compared to the FIG dataset uncertainty (Equation 21). This was done for both east-west and north-south velocities. The closeness was highest at a value of one. ? 2 ?2 = ?2 [21] ?? where ?2 is the reduced chi-squared value representing if there is a statistically significant difference between the observed and calculated values, ?? is the residual between observed (i.e., FIG dataset) and block modeled calculated velocities, and ?? is the FIG dataset uncertainty (inferred from the misfit in the polynomial equation fitting, Equation 3). ?2 > 1 indicates under-fitting, ?2 < 1 indicates over-fitting, and ?2 = 1 indicates best fit. 40 Chapter 3: Results This chapter shows results from the InSAR and GNSS data processing, deformation rate analysis, error analysis, and application of block modeling. 3.1 InSAR Results The SAR image acquisitions utilized to create the interferograms and the subsequent interferogram network collected by Sentinel-1 suggested a relatively moderate and uniform spatial coherence of 0.3. Each image acquisition was used in three interferogram pairs (Figure 13). As expected with improved satellite technology, the variation in perpendicular baseline was minimal with a maximum of 150 m for both the ascending and descending LOS InSAR tracks (Figure 13). Perpendicular baseline refers to the orbital tracks of repeated flybys with respect to the original flight track. 41 Figure 13. Perpendicular baseline with the average spatial coherence between an image acquisition and the selected three image acquisitions to produce pairs for the interferogram network. (A) Ascending interferogram network. (B) Descending interferogram network. 42 Both ascending and descending LOS InSAR mean velocities and time-series produced from MintPy showed negative motions in the central west coast of Taiwan, which suggested surface subsidence (Figure 14). However, from southwest to northeast, there was a gradual transition from positive to negative motion in the ascending track but negative to positive motion in the descending track. This opposite sign of motion between the ascending and descending tracks suggested a dominating horizontal motion. These results were encouraging of more critical analysis to determine specific horizontal and vertical directions of motion. The time-series re- model and mean velocity re-estimation with GNSS correction results suggested consistency between the InSAR and GNSS velocities that were further improved upon with the GNSS correction (Figure 15, 16, 17). With respect to the GNSS velocities, the ascending LOS InSAR velocities were underestimated by ~20 mm/yr, and the descending LOS InSAR velocities were overestimated by ~20 mm/yr. The error within the ascending and descending LOS InSAR velocities, likely due to the difference in reference region from GNSS, were corrected to produce the GNSS-corrected LOS InSAR velocities, and east and west Taiwan were merged (Figure 18). The ascending and descending GNSS-corrected LOS InSAR velocities were converted to east-west, north-south, and vertical velocities and combined with the interpolated GNSS velocities while taking into consideration proper weighting based on misfit (Figure 19). The final product suggested (i) westerly motion of up to 60 mm/yr in east Taiwan due to the northwestward collision of the Philippine Sea Plate into the Eurasian Plate, (ii) northly motion of up to 40 mm/yr along the Longitudinal Valley fault due to collision and southerly motion of up to 25 mm/yr in southwest 43 Taiwan due to tectonic escape, and (iii) subsidence of up to 40 mm/yr in west Taiwan due to water pumping and uplift of up to 20 mm/yr in central Taiwan due to collisional mountain-building events. Determining the east-west, north-south, and vertical mean velocity gradients from the FIG dataset revealed areas of increased strain rate accumulation focused along the Longitudinal Valley fault, the Western Foothills, and the central west coast of Taiwan (Figure 20). However, the increased strain rate accumulation located on the central west coast of Taiwan is likely not tectonic as this area is subject to groundwater pumping. MintPy Mean Velocity MintPy Mean Velocity Ascending, West and East Descending, West and East A B Flight direction with right-look Flight direction with right-look Figure 14. Mean velocities from (A) ascending and (B) descending LOS InSAR tracks produced by MintPy with coordinates in pixels (Yunjun et al., 2019). (A) Negative ascending LOS InSAR velocities indicate easterly motion and/or subsidence. (B) Negative descending LOS InSAR velocities indicate westerly motion and/or subsidence. The black squares indicate reference locations. One pixel is ~50 m x 50 m. West and east Taiwan were calculated separately but are merged above. The region impacted by the 2018 Hualien earthquake was removed (see white rectangle). 44 Figure 15. Uncorrected ascending LOS InSAR velocities overlain with GNSS velocities for (A) west and (B) east Taiwan. Negative ascending LOS InSAR velocities indicate easterly motion and/or subsidence. Uncorrected descending LOS InSAR velocities overlain with GNSS velocities for (C) west and (D) east Taiwan. Negative descending LOS InSAR velocities indicate westerly motion and/or subsidence. 45 Figure 16. GNSS-corrected ascending LOS InSAR velocities overlain with GNSS velocities for (A) west and (C) east Taiwan. Negative ascending LOS InSAR velocities indicate easterly motion and/or subsidence. Comparison of ascending LOS InSAR velocities to GNSS velocities for (B) west and (D) east Taiwan before and after the GNSS correction. 46 Figure 17. GNSS-corrected descending LOS InSAR velocities overlain with GNSS velocities for (A) west and (C) east Taiwan. Negative descending LOS InSAR velocities indicate westerly motion and/or subsidence. Comparison of descending LOS InSAR velocities to GNSS velocities for (B) west and (D) east Taiwan before and after the GNSS correction. 47 Figure 18. (A) Ascending and (B) descending GNSS-corrected LOS InSAR velocities of merged west and east Taiwan. (A) Negative ascending LOS InSAR velocities indicate easterly motion and/or subsidence. (B) Negative descending LOS InSAR velocities indicate westerly motion and/or subsidence. 48 Figure 19. The FIG dataset with (A) eastward, (B) northward, and (C) upward components of motion in mm/yr. Positive values correspond to (A) easterly motion, (B) northerly motion, and (C) uplift. The black lines indicate CGS of Taiwan?s active faults (CGS, 2012). Gray indicates regions of no data. 49 Figure 20. (A) East-west, (B) north-south, and (C) vertical mean velocity gradients derived from the FIG dataset. Increased mean velocity gradient indicates increase in strain rate accumulation. The black lines indicate CGS of Taiwan?s active faults (CGS, 2012). Gray indicates regions of no data. 3.2 Deformation Rate Analysis Results Suggested by the error analysis results (see Section 3.3), incorporating the nearest 144 pixels (smoothing window size) into the deformation rate tensor produced the highest resolution results that distinguished tectonic signal from noise (Figure 21, 22). Tectonic signals were of a minimum resolution of ~3.4 km. Deformation rate analysis results suggested the highest deformation rate characterized by both dilation and maximum shear along the Longitudinal Valley fault in east Taiwan. To a lesser extent, increased deformation was also observed in west Taiwan where the main fold- and-thrust belts (i.e., Western Foothills) are located (Figure 22). The 2nd invariant of strain rate (Figure 22) highlighted regions with greater total strain rate accumulation 50 and implied locked faults. In addition to the Longitudinal Valley fault, several faults in southwest Taiwan appeared to be locked as well. When calculating the deformation rate tensor used in dilation, maximum shear, and 2nd invariant, increasing the number of nearest pixels reduced the noise but subsequently reduced the ability to apply estimated values to specific faults and structures. 51 Figure 21. Deformation rate analysis with (A, B, C) dilatation rate, (D, E, F) maximum shear rate, and (G, H, I) 2nd invariant rate in (A, D, G) 30-pixel resolution, (B, E, H)144-pixel resolution, and (C, F, I) 420-pixel resolution. Positive dilatation values indicate contraction and negative values indicate expansion. High maximum shear values indicate increased shearing. 2nd invariant values describe both dilatation and maximum shear as total strain rate accumulation. Gray indicates regions of no data. 52 Figure 22. (A) Dilatation rate, (B) maximum shear rate, and (C) 2nd invariant rate of Taiwan with a smoothing window of 144 pixels. The black lines indicate active faults categorized by the CGS of Taiwan (CGS, 2012). Gray indicates regions of no data. 53 3.3 Error Analysis Results The RMS misfit estimates from the inverse problem (Section 2.3.4) for mean velocity models of each InSAR pixel and interpolated GNSS throughout time are shown in Figure 23. Misfit was generally < ~3 mm/yr for ascending LOS InSAR, descending LOS InSAR, east-west GNSS, and north-south GNSS. Misfit increased to ~5 mm/yr ?10 mm/yr for vertical GNSS. This increased misfit is likely due to that GNSS receives higher noise levels in vertical components and is more susceptible to atmospheric delay-related errors (Serpelloni et al., 2013). RMS misfits produced from varying ? values related to days of postseismic relaxation from the 2018 Hualien earthquake displayed minimal change in error estimation between the data and model. ? = 121 resulted in a minimized RMS misfit of 2.9245 mm/yr (Figure 24). I further estimated the uncertainty values (i.e., error) produced during the transformation of GNSS-corrected ascending and descending LOS InSAR velocities and interpolated east-west, north-south, and vertical GNSS velocities to the 3-D FIG dataset (see Section 2.3.6) (Figure 25). The east-west and north-south FIG dataset velocities produced the least amount of misfit of < ~2 mm/yr and < ~3 mm/yr, respectively. The vertical FIG dataset velocities displayed a larger range of uncertainties between ~3 mm/yr ? 10 mm/yr. This was likely due to the reliance on the inherently noisier InSAR data in the calculation as GNSS poorly resolves vertical motions and was to be weighted less. 54 Figure 23. RMS misfits produced from the transformation of observed (A) ascending and (B) descending LOS InSAR velocities and interpolated (C) east-west, (D) north-south, and (E) vertical GNSS velocities to modeled velocity values. 55 Figure 24. RMS misfits produced from utilizing various ? values (0 to 600 days in 20-day step sizes) in the transformation of observed velocities to modeled velocities. The yellow star indicates ? = 121, which is associated with the lowest RMS misfit estimation. 56 Figure 25. Uncertainty values produced from transforming GNSS-corrected ascending and descending LOS InSAR velocities and interpolated GNSS velocities to the FIG dataset with (A) east-west, (B) north-south, and (C) vertical velocities. 57 To best interpret the deformation rate analysis (see Section 2.4), the noise structure of a non-deforming region was examined using a semi-variogram model and covariogram model (Figure 26, 27). On the semi-variogram each point signified the variance between each and all pixels at distance. The points were particularly dense between a low variance of 0 to 10 as similarity was displayed. Binning (i.e., where the mean value at intervals is determined) of the points produced a downsampled semi- variogram that displayed the best-fitting, exponential model of variance among pixels. The downsampled semi-variogram was inverted to an exponential covariogram model that reached a critical wavelength at ~2.8 km. This was the maximum wavelength that pixels were spatially correlated due to noise. Therefore, utilizing the nearest 144 pixels (i.e., 3.4 km) in each deformation rate tensor calculation enabled tectonic signal to be distinguished from noise signal (Table 1). 58 Figure 26. Semi-variogram displaying variance as a function of separation distance for each and all pixels. The open red circles represent the relative locations along the x-axis used for binning the data. 59 Figure 27. Downsampled semi-variogram displaying variance between each and all pixels as a function of distance. Overlain is the best fit exponential model. Inverted from the semi-variogram is the exponential covariogram model. 3.4 Block Model Results The original block model construction consisted of 50 blocks with each block segment having an assigned locking depth and dip angle (Table 2; Figure 28, 29). Block segments were better constrained in west Taiwan due to the focus of previous research and study efforts. Utilizing Blocks TVR, I examined an output suite of block model constructions with ? from 1 to 2000 in step sizes of 5. The reference model was located at the first notable drop in independent blocks where the MRV was similar to the FIG dataset uncertainty of ~ 3 mm/yr, and the reduced chi-squared result suggested slight 60 under-fitting. The reference block model was produced with a ? value of 200 and 42 blocks (Figure 30). The subsequent MRV was 2.95 mm/yr. The removed block segments, from regularization, were in northern Taiwan (Figure 31). The modeled total slip rates produced by the reference block model highlighted faults in both west and east Taiwan (Figure 31). In Western Foothills, the Chiuchiungkeng fault (#19 in Table 2), Chushiang structure (39), and Gukeng structure (40) had an increased total slip rate between 40 mm/yr ? 55 mm/yr. South of these faults and structures, there was an increased slip of up to 60 mm/yr on various, sporadic faults. Within the Longitudinal Valley, the Longitudinal Valley fault (33) experienced a total slip rate between 40 mm/yr ? 50 mm/yr. Left-lateral strike-slip rates of up to 40 mm/yr were found within the Western Foothills on the Kaoping River structure (28) and Longchuan structure (42) and within the Longitudinal Valley on the Longitudinal Valley fault (33). Opposing strike-slip motions were located on the Northern Ilan structure (37) (i.e., left-lateral) and the Southern Ilan structure (38) (i.e., right-lateral), which indicated rotation of the Ilan basin due to the opening of the Okinawa Trough. This type of behavior is different than that observed on faults and structures bounding the Pingtung Plain, which is undergoing southwest tectonic escape. Additionally, right-lateral motion of up to 30 mm/yr was observed along an east-west trending transect through the Hsuehshan Range and Coastal Range. This may be linked to the placement of a broad shear zone. Convergent motion was most prominent along the Longitudinal Valley fault (33), up to 30 mm/yr, where the Philippine Sea Plate collides with the Eurasian Plate. 61 Divergent motion of up to 20 mm/yr was focused along the Northern Ilan structure (37) and Southern Ilan structure (38). The highest east-west velocity residuals (up to 10 mm/yr) produced from the reference block model construction were focused in the mountainous region of Taiwan and the Longitudinal Valley (Figure 32). These locations contain faults and structures that are less constrained and due for improved mapping. Chang et al. (2016) suggest residual velocities of up to 20 mm/yr in these regions. The highest north-south velocity residuals were in the Ilan Basin and the southern portion of the Coastal Range. The Ilan Basin displays an opening and rotating motion that is likely under the influence of more complicated fault and structure dynamics than what is currently understood. The southern portion of the Coastal Range is likely experiencing loosely defined fault branching. Table 2. Faults of Taiwan (CGS, 2012; Shyu et al., 2020) with associated locking depths and dip angles (Chan et al., 2020). Fault Name Number Locking Depth Dip Angle Shanchiao fault 1 14.0 48.2 Shuanglienpo structure 2 5.0 33.0 Yangmei structure 3 3.0 60.0 Hukou fault 4 10.0 30.0 Fengshan River strike-slip structure 5 13.8 85.0 Hsinchu fault 6 10.0 45.0 Hsincheng fault 7 12.9 30.0 Hsinchu frontal structure 8 10.0 30.0 Touhuanping structure 9 12.0 85.0 Miaoli frontal structure 10 10.0 30.0 Tunglo structure 11 3.5 30.0 East Miaoli structure 12 4.0 30.0 Shihtan fault 13 10.8 75.0 Sanyi fault 14 9.0 15.0 62 Tuntzuchiao fault 15 14.8 85.0 Changhua fault 16 12.0 22.1 Chelungpu fault 17 12.0 15.0 Tamaopu - Shuangtung fault 18 6.0 30.0 Chiuchiungkeng fault 19 12.0 30.0 Meishan fault 20 14.7 85.0 Chiayi frontal structure 21 12.0 15.0 Muchiliao - Liuchia fault 22 12.0 30.0 Chungchou fault 23 12.0 30.0 Hsinhua fault 24 15.0 85.0 Houchiali fault 25 5.0 45.0 Chishan fault 26 10.8 75.0 Hsiaokangshan fault 27 7.0 30.0 Kaoping River structure 28 12.3 75.0 Chaochou fault 29 11.1 75.0 Hengchun fault 30 15.0 75.0 Hengchun offshore structure 31 4.0 30.0 Milun fault 32 10.0 75.0 Longitudinal Valley fault 33 20.0 55.0 Central Range structure 34 20.0 45.0 Luyeh fault 35 4.0 37.5 Taimali coastline structure 36 10.6 75.0 Northern Ilan structure 37 9.4 60.0 Southern Ilan structure 38 11.3 60.0 Chushiang structure 39 3.0 55.0 Gukeng structure 40 12.0 85.0 Tainan frontal structure 41 12.0 18.8 Longchuan structure 42 12.0 60.0 Youchang structure 43 12.0 75.0 Fengshan Hills frontal structure 44 15.0 30.0 Fengshan structure 45 N/A N/A 63 Figure 28. Faults of Taiwan used in the block model construction (CGS, 2012; Shyu et al., 2020). 64 Figure 29. Original block model construction and associated (A) locking depths and (B) dip angles (Chan et al., 2020). 65 Figure 30. Evaluation of block model suite produced from TVR. Displayed are the MRVs, chi-squared values, and the number of blocks in the output models. The star indicates the reference block model. 66 Figure 31. (A) Total slip rates produced from reference block model construction with (B) positive strike-slip rates indicating right- lateral motion and (C) positive convergence rates indicating convergent motion. The removed block segments appear in black. 67 Figure 32. (A) East-west, (B) north-south and (C) total residual velocity values produced from the reference block model construction. The preserved block segments are in black, and the removed block segments are in gray. 68 Chapter 4: Discussion The main goals of this chapter are to characterize interseismic crustal deformation and quantify the slip budget to highlight seismic hazard potential in Taiwan. The deformation rate analysis produced from the Reduced FIG dataset suggested complementary, more detailed findings compared to previous research based on GNSS-only datasets (e.g., Chang et al., 2016; Ching et al., 2015). The results indicated that the highest velocity values and total strain rate accumulation suggested by the 2nd invariant of strain rate were focused in the Longitudinal Valley, the Western Foothills (including the Pingtung Plain), the Ilan Basin, and the Taipei Basin. The block model further examined these areas of interest and determined that faults located in the Longitudinal Valley and Western Foothills had the greatest total slip rate. 4.1 Combining Geodetic Techniques GNSS has provided high accuracy observations of crustal deformation in Taiwan since the first deployment in the early 1990s (Yu et al., 1997; Ching et al., 2007; Hsu et al., 2009; Lin et al., 2010). Stations are located roughly every tens of kilometers and can readily detect uplift, subsidence, and horizontal motions (Ching et al., 2011). Although Taiwan has one of the highest GNSS deployment densities in the world, GNSS measurements alone do not allow us to relate detected ground motions to specific faults and structures. Therefore, combining the high resolution of InSAR and the high accuracy of GNSS is worthwhile as it enables accurate ground motion velocities at high spatial resolutions to be estimated and potential locked faults to be detected (B?rgmann et al., 2000; Elliott et al., 2016). Potential locked faults can be 69 verified with additional lines of evidence (i.e., mean velocity gradient, 2nd invariant of strain rate, geologic observations, etc.). Identifying and confirming potentially active faults is important for quantifying fault slip budget and partitioning throughout the island and forecasting seismic hazard with increased accuracy. 4.2 Evaluating Data and Hazard 4.2.1 Data Comparison Huang et al. (2016) estimate crustal velocities using InSAR and CGNSS for west Taiwan. Broadly, I shared similar overall findings for tectonic escape in southwest Taiwan and land subsidence in west Taiwan (Figure 33). Southwest Taiwan displayed a uniform descending LOS InSAR velocity of ~-25 mm/yr which indicated either subsidence or westerly motion from 1995 to 1999 using European Remote-Sensing (ERS) satellite data and 2005-2008 using Envisat satellite data (Huang et al., 2016). My 2016-2021 findings suggested a sporadic increase and decrease in descending LOS InSAR velocity throughout southwest Taiwan by ~10 to 20 mm/yr (Figure 33). These changes in velocity may be due to stress field alterations from the Mw 6.3 and 6.4 earthquakes in 2010 and 2016, respectively. Similar findings were documented by Hsu et al. (2018). Alternatively, they could also be due to the presence of phase unwrapping errors of InSAR data during 1995-1999 (Huang et al., 2016), which may be due to the fewer number of SAR acquisitions (< 10 images per year) during that time period. That is, fewer SAR acquisitions could cause increased phase unwrapping errors in vegetated terrains such as in Taiwan. Additionally, the spatial and temporal baseline between pairs using ERS and Envisat (35 days minimum) was generally greater than those using 70 Sentinel-1 data (6 days minimum). A greater spatial baseline would cause lower coherence and could hinder the production of high phase-quality interferograms. Ultimately, InSAR during this time period may underestimate LOS InSAR velocities. Comparing descending LOS InSAR velocities along the central west coast of Taiwan between 1995-1999 (Huang et al., 2016) and 2016-2021, maximum subsidence was detected in Yunlin County and increased over time (Figure 32). This may be indicative of either increased groundwater pumping or exponentially growing groundwater pumping consequences in relation to soil compaction (Holzer, 1984). 71 Figure 33. The (C) difference between the (A) descending LOS InSAR velocities from 2016-2021 using InSAR- and GNSS-derived data and (B) descending LOS InSAR velocities from 1995-1999 using InSAR-derived data (Huang et al., 2016). Negative LOS InSAR velocities (brown) indicate either subsidence and/or westerly motion. (C) Negative LOS InSAR velocity change (red) indicates increase in either subsidence and/or westerly motion. The descending LOS InSAR velocities from Huang et al. (2016) have been adjusted to accommodate for the difference in reference point 72 4.2.2 Hazard Comparison Utilizing the interpolated GNSS velocity data from Lin et al. (2010), I calculated 2nd invariant rates over Taiwan both before (1990-1995) and after (2003- 2005) the 1999 MW 7.6 Chi-Chi earthquake (Figure 34). In comparison to the 2016- 2021 2nd invariant rates estimated from the Reduced FIG dataset, I observed similar patterns to that of the post-Chi-Chi earthquake with observations being limited due to noise associated with InSAR data collection (i.e., topography correlation, atmospheric noise, etc.). Similarly, I discerned contraction along the Longitudinal Valley Fault and within a band encompassing a section of the Western Foothills. This band may be influenced by the Pingtung Plain, which is experiencing tectonic escape. Compared to the 2005-2009 InSAR and CGNSS derived 2nd invariant rates from Huang and Evans (2019), the 2nd invariant rates during 2016-2021 show minimal variation (Figure 34C and 34D). The north-south trending region of increased total strain rate accumulation in the Western Foothills persists throughout time. Lin et al. (2010) suggest that large earthquakes and their postseismic deformation may influence island-wide patterns of strain accumulation. For example, there was a minimal similarity to the dilation pattern observed before and after the 1999 Chi-Chi earthquake. Moreover, differencing the 2nd invariant of strain rate from 2016- 2021 (InSAR- and GNSS-derived) and 2003-2005 (GNSS-derived) (Lin et al., 2010) suggested an increase in total strain rate accumulation of up to 20 mm/yr propagating northward in the Longitudinal Valley (Figure 35). Regions in the Western Foothills and along the southern boundary of the Taipei Basin displayed an increase in total strain 73 rate accumulation as well. However, the GNSS-derived 2nd invariant of strain rate values has decreased spatial sampling which may introduce inconsistencies when calculating the change in total strain rate accumulation. 74 Figure 34. 2nd invariant of strain rate throughout time highlighting before (A) 1990-1995 (Yu et al., 1997) and after (B) 2003-2005 (Lin et al., 2010), (C) 2005-2009 (Huang et al., 2016), (D) 2016-2021 the 1999 Chi-Chi earthquake. (A, B) 2nd invariant of strain rate from GNSS-derived data. (C) 2nd invariant of strain rate from InSAR- and CGNSS-derived data. (D) 2nd invariant of strain rate from InSAR- and GNSS-derived data. The black dots indicate GNSS station locations. 75 Figure 35. The difference between the (A) 2nd invariant of strain rate from 2016-2021 using InSAR- and GNSS-derived data and the (B) 2nd invariant of strain rate from 2003-2005 using GNSS-derived data (Lin et al., 2010). The black dots indicate GNSS station locations. (C) Blue indicates an increase in total strain rate accumulation given an increase in 2nd invariant of strain rate. 76 Chan et al. (2020) provide the 2020 Taiwan Earthquake Model (TEM) of Probabilistic Seismic Hazard Analysis (PSHA). The PSHA model was calculated based on three parts: seismic-hazard source model, ground motion model, and probabilistic calculation (Field, 2001). The seismic-hazard source model is a list defining earthquake scenarios using previous event details (i.e., magnitude, location, depth, seismicity rate). The ground motion model, or the basic attenuation relationship, specifies the ground motion levels given distance from a certain magnitude event. More complex models allow for the inclusion of rock type and source faulting style. They can also compute scenarios with combined fault-segments ruptures. The probabilistic calculation determines the distribution of possible ground-motion levels given an earthquake scenario to estimating the probability of exceedance of ground-motion levels over a certain span of time. The previous and original TEM PSHA was produced in 2015 (Wang et al., 2016) based on a seismogenic structure source database proposed by Shyu et al. (2016). Different from its predecessor, the TEM PSHA2020 improved upon the seismogenic structure database, incorporated the probability of rupture along multiple structures, included the effect of fault memory given the time-lapse from the previous rupture, updated the earthquake catalog used in defining area sources, smoothed the utilized seismicity rates, utilized more accurate ground motion prediction equations, and incorporated site effect (i.e., shear wave amplification given subsurface characteristics). The final analysis with site effect inclusion specifically suggested increased seismic hazard potential in the Taipei Basin, Ilan Basin, and Chianan Basin (Figure 36). Additionally, observations suggested increased hazard potential along the 77 Longitudinal Valley Fault and within the Western Foothills as well. Chan et al. (2020) suggest a 10% probability of 0.4 g, 1.4 g, and 0.75 g peak ground motion in 50 years for Taipei Basin, Ilan Basin, and Chianan Basin, respectively. Given the high-resolution results of this study, I was better able to define the specific faults, whether they were identified by CGS of Taiwan (CGS, 2012) or inferred by Shyu et al. (2020), for which the increased seismic hazard potential was associated. In the Taipei Basin, the 2nd invariant of strain rate did not observe clear total strain rate accumulation along the Shanchiao fault (1) (CGS, 2012). However, the findings revealed a potential locked fault extending along the basin?s southern boundary which has been inferred as the Hsindian fault (e.g., Chen, 2016). This inferred fault, when compared to DEM and relief data, does not appear to be an artifact of topography. In the Ilan Basin, both probabilistic seismic hazard for 10% in 50 years and 2nd invariant of strain rate highlighted increased seismic hazard potential along the Northern Ilan structure and the Southern Ilan structure. Within the Chianan Basin, the 2nd invariant of strain rate did not suggest increased total strain rate accumulation along the inferred Chiayi frontal structure (21) (Shyu et al., 2020). This structure showed a more prominent increase in seismic hazard potential when examining probabilistic seismic hazard for 2% in 50 years, rather than 10%. East of the Chiayi frontal structure, I observed increased total strain accumulation along mapped faults trending north-south in the Western Foothills (i.e, Chushiang structure, 39; Gukeng structure, 40; Chiuchiungkeng fault, 19). 78 Figure 36. Comparison of regions with increased seismic hazard suggested by Chan et al. (2020) and 2016-2021 2nd invariant of strain rate. From North to South: Taipei Basin, Ilan Basin, Chianan Basin. The black lines indicate active faults identified by the Central Geologic Survey of Taiwan (CGS, 2012). The blue lines indicate potential faults. 79 4.2.3 Block Model Comparisons When comparing the modeled total slip rates from the reference block model (Figures 31 and 37) to the total slip rates estimated by Chan et al. (2020), I observed encouraging similarities in relative velocities. The increased total slip rates corresponded with greater (i.e., above average) recorded total slip rates estimated by Chan et al. (2020). The Ilan Basin, the Northern Ilan structure (37) and the Southern Ilan structure (38) had recorded total slip rates of 3.29 mm/yr and 5.48 mm/yr (Chan et al., 2020), respectively. The modeled total slip rate for these structures was ~30 mm/yr. In the Longitudinal Valley and east Central Range, the Longitudinal Valley fault (33) and Central Range structure (34) had recorded total slip rates of 11.35 mm/yr and 7.28 mm/yr, respectively. The modeled total slip rates for this fault and structure were ~45 mm/yr. In contrast to the modeled total slip rates (Figure 31 and 37), the Chushiang structure (39) Gukeng structure (40), and Chiuchiungkeng fault (19) in the Western Foothills did not collectively have an above-average recorded total slip rate. Respectively, they had recorded total slip rates of 2.90 mm/yr, 0.94 mm/yr, and 4.66 mm/yr. Only the Chiuchiungkeng fault (19) had above-average slip rate. The modeled values suggested between 40 mm/yr ? 50 mm/yr of total slip rate. This indicates the Gukeng structure (40) experiences notably more motion than previously observed. Unexpectedly, the block segment trending east-west across the Hsuehshan Range and Central Range displayed increased total slip rate of up to 30 mm/yr 80 (Figures 31 and 37). This block segment was included due to necessary block division as recommended from previous block models (Ching et al., 2011, Chang et al., 2016), and did not correspond to any known faults. That is, a change in velocity vector direction required subdivision of north Taiwan (Ching et al., 2011). Potentially, this increased total slip rate is linked to a broad shear zone. While the total slip rate trends (i.e., both modeled and recorded rates) on the major faults and structures of interests did correspond in trend, there is overall discrepancy. This discrepancy may be a result of examining convergence on different timescales (Shyu et al., 2006; Shyu et al., 2016). The modeled total slip rates only account for convergence between 2016 and 2021. The recorded total slip rates (Chan et al., 2020) are the result of a combination of field investigations and geodetic studies. That is, the recorded total slip rates detail millennial rates of convergence rather than just a 5-year period. 81 Figure 37. Faults of Taiwan used in the block model construction (CGS, 2012; Shyu et al., 2020), and the total slip rates produced from the reference block model. The reference block model construction of Taiwan with 42 blocks shared similarities to that of Huang and Evans (2019). Concerning block model extent and evaluation, our reference block model included the entire island, rather than only the southwest section, and suggested a decrease in MRV from 3.6 mm/yr (Huang and Evans, 2019) to 2.95 mm/yr. The reduction in MRV is likely due to increased data points resulting in an updated block model construction and better quality datapoints due to a longer observational time period. Both block models expressed similar 82 maximum velocity residuals of up to 10 mm/yr (Huang and Evans, 2019) and 12 mm/yr. Both block models suggested up to 40 mm/yr of general strike-slip motion in Western Foothills and increased left-lateral strike-slip motion in the Pingtung Plain. The left-lateral motion related to the Pingtung Plain appears to be focused on different fault segments for Huang and Evans (2019) and the reference model. However, they both suggest a maximum left-lateral motion of ~40 mm/yr linked to the tectonic escape. When comparing the reference block model to the TEM PSHA2020, I found additional evidence for increased seismic hazard in the Ilan Basin. However, I did not find a total slip rate increase in the northern Taipei Basin and the Chianan Basin (see Figure 36 for locations). In the Ilan Basin, the reference block model suggested an increased total slip rate along the Northern Ilan structure and Southern Ilan structure with 20 mm/yr of left-lateral and right-lateral motion, respectively. This duality of motion suggests rotation and opening of the Ilan Basin. Hu et al. (2001) suggest that the motion is likely aligned with divergence given the basin?s proximity to the opening of the Okinawa Trough. The reference block model suggested that faults in the Western Foothills (i.e, Chushiang structure, 39; Gukeng structure, 40; Chiuchiungkeng fault, 19) are experiencing an increased total slip rate that is dominantly convergent. 4.3 Implications: Detected Potential Locked Faults Based on the FIG dataset, mean velocity gradients, and 2nd invariant of strain rate, I identified three potential active, locked faults in the north and southwest Taiwan 83 (Figure 38, 39, 40. Some of these potential faults have been previously identified, but they are considered inactive or unconfirmed active (Teng et al., 2001; Huang et al., 2007). These potential faults each maintain at least three lines of evidence. However, to fully confirm the activity of these faults, the following avenues should be explored as well: paleoseismological data, further field geology investigations (e.g., outcrops displaying faulting behaviors/features), and building structure damage due to fault creep. Figure 38. Evidence for the potential Taipei fault (i.e., Hsindian fault) along the southern boundary of the Taipei Basin. The black lines indicate active faults identified by the Central Geologic Survey of Taiwan (CGS, 2012). The red lines (solid and dashed) represent the potential fault. 84 Figure 39. Evidence for potential Liuchia Fault or extension of the Muchiliao ? Liuchia Fault (22) between Tainan City and Kaohsiung City in the Western Foothills. The black lines indicate active faults identified by the Central Geologic Survey of Taiwan (CGS, 2012). The red lines (solid and dashed) represent the potential fault. 85 Figure 40. Evidence for a potential fault between Tainan City and Kaohsiung City in the Western Foothills. The southern portion of the potential fault may be a part of the Fengshan structure (45). The black lines indicate active faults identified by the Central Geologic Survey of Taiwan (CGS, 2012). The red lines (solid and dashed) represent the potential fault. Along the southern boundary of the Taipei Basin, I identified the Taipei fault as a locked fault because of observed vertical velocity change and increased total strain rate accumulation. The Taipei fault has been identified as an inactive reverse fault in the southern portion of the basin (Teng et al., 2001; Huang et al., 2007), and poses risk to Taipei City: the capital of Taiwan and a populated urban center. Previous studies suggest the observed ground motion behavior is a result of groundwater pumping inducing soil compaction, aquifer deformation, and general subsidence (Chen et al., 2007). Although surface creep along the Taipei fault is observed, I cannot discern 86 whether it is surface subsidence due to groundwater pumping or fault creep. Future studies on seasonal variation of surface movement and how it relates to precipitation and groundwater discharge data may provide further insight into identifying the cause of surface creep. The potential locked Liuchia Fault or extension of the Muchiliao ? Liuchia Fault (22) (Shyu et al. 2020; Chan et al., 2020) and the potential locked fault between Tainan City and Kaohsiung City located in the Western Foothills were identified due to a sudden east-west velocity change (i.e., change in east-west velocity rate and mean velocity gradient) and an increase in total strain rate accumulation. Collectively, they pose risk to Tainan City and Kaohsiung City, the 4th and 2nd largest cities in Taiwan. The potential Liuchia Fault may be an extension of the Muchiliao ? Liuchia Fault (22), and, if so, then the Muchiliao ? Liuchia Fault with increased horizontal length and displacement area would be able to achieve events of different magnitudes and frequencies. The potential fault between Tainan City and Kaohsiung City, which is of unusual geometry, may be multiple faults, extensions or branches of previously mapped faults (e.g., Fengshan structure (45)), and/or reactivations of historic faults. 87 Chapter 5: Conclusion Taiwan is at the dynamic center of two rapidly converging plates. As a result, the island hosts a complex system of faults with increased earthquake activity that provides scientists with a unique natural laboratory to study active tectonics. In the past few decades, alongside routine lower magnitude earthquake events, Taiwan has experienced notable inland earthquakes including the 1999 Mw 7.6 Chi-Chi earthquake and 2016 Mw 6.4 MeiNong earthquake. Both events left an immense trail of destruction. Therefore, it is important to understand the fault and structure systems of the island, past and present crustal deformation velocities, and the consequences of slip budget deficits on individual faults. Combining the capabilities of GNSS and InSAR, I can more clearly define interseismic crustal deformation of Taiwan. Together, these geodetic techniques produce high accuracy velocities with high spatial resolution (i.e., mm-scale). The final InSAR and GNSS dataset (i.e., FIG dataset) reveals up to 60 mm/yr of westerly motion in east Taiwan due to the northwestward collision of the Philippine Sea Plate into the Eurasian Plate, (ii) up to 40 mm/yr of northly motion along the Longitudinal Valley fault due to collision and up to 25 mm/yr of southerly motion in southwest Taiwan due to tectonic escape, and (iii) up to 40 mm/yr of subsidence in west Taiwan due to water pumping and up to 20 mm/yr of uplift in central Taiwan due to collisional mountain-building events. Using these interseismic velocity estimates, I perform a deformation rate analysis that suggests increased deformation from both dilatation and maximum shear along the Longitudinal Valley fault and within the Western Foothills. The 2nd 88 invariant of strain rate highlights regions of greatest total stain rate accumulation and potential locked faults. The error analysis improves confidence in the study and its progressions. Misfit values from the inverse problem used to determine the best-fitting model of each InSAR pixel and interpolated GNSS throughout time produce average values of 3 mm/yr and maximum values of 10 mm/yr. Uncertainty values estimated during the transformation from LOS to 3-D deformation produce average values of 3 mm/yr and maximum values of 10 mm/yr. To improve data interpretation, I used an exponential model to simulate atmospheric noise structure in Taiwan from a region without known active faulting. I found that surface deformation features with a wavelength of ~2.8 km could be interpreted as a tectonic signal. By applying the calculated interseismic crustal deformation velocities and deformation rate analysis observations to block model construction and analysis, I determine a reference block model identifying 42 unique blocks with an MRV of 2.95 mm/yr. The reference model suggests total slip rates along faults and structures throughout all of Taiwan. I can interpret these total slip rates to forecast localities, faults, and structures with increased seismic hazard. However, this type of model does not consider fault slip deeper than the main detachment of Taiwan. As a result, this earthquake hazard model cannot account for deep events such as the 2010 Jia- Shian earthquake and the 2016 MeiNong earthquake. Future work involves combining the geodetic data (i.e., slip deficit accumulation rate) of this study with the Gutenberg-Richter Relationship of Taiwan (i.e., interseismic seismicity rate) to determine the maximum anticipated earthquake 89 for Taiwan. Additionally, extending block modeling capabilities to include vertical velocity observations will be of interest as Taiwan experiences both substantial horizontal and vertical motions. 90 Bibliography Angelier, J., Barrier, E., & Chu, H. T. (1986). Plate collision and paleostress trajectories in a fold-thrust belt: The foothills of Taiwan. Tectonophysics, 125(1-3), 161-178. doi:10.1016/0040-1951(86)90012-0. Angelier, J., Chu, H., Lee, J., & Hu, J. (2000). Active faulting and earthquake hazard: The case study of the Chihshang Fault, Taiwan. Journal of Geodynamics, 29(3-5), 151-185. doi:10.1016/s0264-3707(99)00045-9 Berardino, P., Fornaro, G., Lanari, R., and Sansosti, E. (2002), Monitoring based on small baseline differential SAR interferograms, IEEE Transactions on Geoscience and Remote Sensing, 40, 2375-2383. Bu?rgmann, R., Rosen, P. A., & Fielding, E. J. (2000). Synthetic Aperture Radar Interferometry to Measure Earth?s Surface Topography and Its Deformation. Annual Review of Earth and Planetary Sciences, 28(1), 169-209. doi:10.1146/annurev.earth.28.1.169. Central Geologic Survey (CGS), (2012). Active fault of Taiwan, MOEA. Retrieved from https://fault.moeacgs.gov.tw/TaiwanFaults_2009/ PageContent.aspx?type=C&id=5 Chan, C.-H., Ma, K.-F., Shyu, J. B. H., Lee, Y.-T., Wang, Y.-J., Gao, J.-C., Yen, Y.- T., & Rau, R.-J. (2020). Probabilistic seismic hazard assessment for Taiwan: TEM PSHA2020. Earthquake Spectra, 36, 137?159. https://doi.org/10.1177/8755293020951587 Chang, W. L., Ching, K. E., Lee, C. H., Lee, Y. R., & Lee, C. F. (2016). Earthquake potential of active faults in Taiwan from GPS observations and block modeling. Seismological Research Letters, 87(6), 1274?1286. https://doi.org/10.1785/0220160094 Chen, C. W., & Zebker, H. A. (2002). Phase unwrapping for large SAR interferograms: Statistical segmentation and Generalized Network Models. IEEE Transactions on Geoscience and Remote Sensing, 40(8), 1709? 1719. https://doi.org/10.1109/tgrs.2002.802453 91 Chen, C.-T., Hu, J.-C., Lu, C.-Y., Lee, J.-C., & Chan, Y.-C. (2007). Thirty-year land elevation change from subsidence to uplift following the termination of groundwater pumping and its geological implications in the metropolitan taipei basin, Northern Taiwan. Engineering Geology, 95(1-2), 30?47. https://doi.org/10.1016/j.enggeo.2007.09.001 Chen, W. S. (2016), Introduction to Taiwan geology (in Chinese), Geologic Society of Taiwan Ching, K.-E., Rau, R.-J., Lee, J.-C., & Hu, J.-C. (2007). Contemporary deformation of tectonic escape in SW Taiwan from GPS observations, 1995-2005. Earth and Planetary Science Letters, 262(3-4), 601?619. https://doi.org/10.1016/j.epsl.2007.08.017 Ching, K.-E., K. M. Johnson, R.-J. Rau, R. Y. Chuang, L.-C. Kuo, and P.-L. Leu (2011), Inferred fault geometry and slip distribution of the 2010 Jiashian, Taiwan, earthquake is consistent with a thick-skinned deformation model, Earth Planet. Sci. Lett., doi:10.1016/j.epsl.2010.10.021. CPC Corporation (1989a). Geologic map of western Taiwan, Chiayi sheet, Taiwan Petroleum Exploration Division, CPC Corporation scale 1:100,000. CPC Corporation (1989b). Geologic map of western Taiwan, Miaoli sheet, Taiwan Petroleum Exploration Division, CPC Corporation scale 1:100,000. CPC Corporation (1989c). Geologic map of western Taiwan, Taichung sheet, Taiwan Petroleum Exploration Division, CPC Corporation scale 1:100,000. CPC Corporation (1989d). Geologic map of western Taiwan. Tainan sheet, Taiwan Petroleum Exploration Division, CPC Corporation scale 1:100,000. Elliott, J. R., Walters, R. J., and Wright, T. J., (2016), The role of space-based observation in understanding and responding to active tectonics and earthquakes, Nat. Commun., 7, 13844, doi:10.1038/ncomms13844. Evans, E. L., Loveless, J. P., & Meade, B. J. (2015). Total variation regularization of geodetically and geologically constrained block models for the Western United States. Geophysical Journal International, 202(2), 713? 727. https://doi.org/10.1093/gji/ggv164 92 Fagereng, ?., & Biggs, J. (2019). New perspectives on ?geological strain rates? calculated from both naturally deformed and actively deforming rocks. Journal of Structural Geology, 125, 100?110. https://doi.org/10.1016/j.jsg.2018.10.004 Farr, T. G., et al. (2007), The shuttle radar topography mission, Rev. Geophys., 45, RG2004, doi:10.1029/2005RG000183. Fattahi, H., Agram, P., & Simons, M. (2017). A Network-Based Enhanced Spectral Diversity Approach for TOPS Time-Series Analysis. IEEE Transactions on Geoscience and Remote Sensing, 55(2), 777-786. doi:10.1109/tgrs.2016.2614925 Field, E. (2001). Probabilistic Seismic Hazard Analysis ( PSHA ) A Primer. Funning, G. (2021, August 23). Introduction to geophysical modeling of geodetic data. Lecture presented at 2021 InSAR Processing and Time-Series Analysis for Geophysical Applications: ISCE, ARIA-Tools, and MintPy short course. Ho, C.S. (1986). A Synthesis of the Geologic Evolution of Taiwan. Tectonophysics, vol. 125, no. 1-3, 1986, pp. 1?16., https://doi.org/10.1016/0040- 1951(86)90004-1. Holzer, T. L. (1984). Ground failure induced by ground-water withdrawal from unconsolidated sediment. Reviews in engineering geology, 6, 67-105. Hsu, Y.-J., Yu, S.-B., Simons, M., Kuo, L.-C., & Cheng, H.-Y. (2009). Interseismic crustal deformation in the Taiwan plate boundary zone revealed by GPS observations, seismicity, and earthquake focal mechanisms. Tectonophysics, 479(1-2), 4?18. https://doi.org/10.1016/j. tecto.2008.11.016 Hsu, W., Byrne, T. B., Ouimet, W., Lee, Y., Chen, Y., Soest, M. V., & Hodges, K. (2016). Pleistocene onset of rapid, punctuated exhumation in the eastern Central Range of the Taiwan orogenic belt. Geology, 44(9), 719-722. doi:10.1130/g37914.1 Hsu, Y.-J., Lai, Y.-R., You, R.-J., Chen, H.-Y., Teng, L. S., Tsai, Y.-C., et al. (2018). Detecting rock uplift across southern Taiwan mountain belt by integrated GPS and leveling data. Tectonophysics, 744, 275?284. https://doi.org/10.1016/j.tecto.2018.07.012 93 Hu, J.-C., Yu, S.-B., Angelier, J., & Chu, H.-T. (2001). Active deformation of Taiwan from GPS measurements and numerical simulations. Journal of Geophysical Research: Solid Earth, 106(B2), 2265?2280. https://doi.org/10.1029/2000jb900196 Hu, J., Hou, C., Shen, L., Chan, Y., Chen, R., Huang, C., . . . Nien, P. (2007). Fault activity and lateral extrusion inferred from velocity field revealed by GPS measurements in the Pingtung area of southwestern Taiwan. Journal of Asian Earth Sciences, 31(3), 287-302. doi:10.1016/j.jseaes.2006.07.020 Huang, S.-Y., Rubin, C. M., Chen, Y.-G., & Liu, H.-C. (2007). Prehistoric earthquakes along the Shanchiao Fault, Taipei Basin, Northern Taiwan. Journal of Asian Earth Sciences, 31(3), 265?276. https://doi.org/10.1016/j.jseaes.2006.07.025 Huang, M.-H., & Evans, E. L. (2019). Total variation regularization of geodetically constrained block models in Southwest Taiwan. Journal of Geophysical Research: Solid Earth, 124. https://doi.org/10.1029/2019JB018076 Hussain, E., Hooper, A., Wright, T. J., Walters, R. J., & Bekaert, D. P. (2016). Interseismic strain accumulation across the central North Anatolian Fault from iteratively unwrapped InSAR measurements. Journal of Geophysical Research: Solid Earth, 121(12). https://doi.org/10.1002/2016JB013108 Kanamori, H. (1994). Mechanics Of Earthquakes. Annual Review of Earth and Planetary Sciences, 22(1), 207-237. doi:10.1146/annurev.ea.22.050194.001231 Lacombe, O., Mouthereau, F., Angelier, J., & Deffontaines, B. (2001). Structural, geodetic and seismological evidence for tectonic escape in SW Taiwan. Tectonophysics, 333(1-2), 323?345. https://doi.org/10.1016/s0040- 1951(00)00281-x Lin, K.-C., Hu, J.-C., Ching, K.-E., Angelier, J., Rau, R.-J., Yu, S.-B., et al. (2010). GPS crustal deformation, strain rate, and seismic activity after the 1999 chi- chi earthquake in Taiwan. Journal of Geophysical Research, 115(B7), B07404. https://doi.org/10.1029/2009JB006417 94 Massonnet, D., Rossi, M., Carmona, C., Adragna, F., Peltzer, G., Feigl, K., & Rabaute, T. (1993). The displacement field of the Landers earthquake mapped by radar interferometry. Nature, 364(6433), 138-142. doi:10.1038/364138a0 Massonnet, D., & Feigl, K. L. (1998). Radar interferometry and its application to changes in the Earths surface. Reviews of Geophysics, 36(4), 441-500. doi:10.1029/97rg03139 Meade, B. J., & Loveless, J. P. (2009). Block modeling with connected fault-network geometries and a linear elastic coupling estimator in spherical coordinates. Bulletin of the Seismological Society of America, 99(6), 3124?3139. https://doi.org/10.1785/0120090088 Naik, S. P. (2017). Earthquake fault propagation, displacement and damages. Lecture presented at 8th International INQUA Meeting on Paleoseismology in New Zealand. Pagani, C., Bodin, T., M?tois, M., & Lasserre, C. (2021). Bayesian Estimation of Surface Strain Rates From Global Navigation Satellite System Measurements: Application to the Southwestern United States. Journal of Geophysical Research: Solid Earth, 126(6). doi:10.1029/2021jb021905 Rosen, P. A., Gurrola, E., Sacco, G. F., Zebker, H. (2012), The InSAR scientific computing environment, Proc. EUSAR, Nuremberg, Germany, 730?733. Serpelloni, E., Faccenna, C., Spada, G., Dong, D., & Williams, S. D. (2013). Vertical GPS ground motion rates in the euro-mediterranean region: New evidence of velocity gradients at different spatial scales along the Nubia- eurasia plate boundary. Journal of Geophysical Research: Solid Earth, 118(11), 6003?6024. https://doi.org/10.1002/2013jb010102 Shyu, J. B. H., Chuang, Y. R., Chen, Y. L., Lee, Y. R., & 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. Terrestrial, Atmospheric & Oceanic Sciences, 27(3). 95 Shyu, J. B., Yin, Y.-H., Chen, C.-H., Chuang, Y.-I., & Liu, S.-Ch. (2020). Updates to the on-land seismogenic structure source database by the Taiwan Earthquake Model (TEM) project for seismic hazard analysis of Taiwan. Terrestrial Atmospheric and Oceanic Sciences. 31, 469-478. 10.3319/TAO.2020.06.08.01. Sudhaus, H., & J?nsson, S. (2009). Improved source modelling through combined use of Insar and GPS under consideration of correlated data errors: Application to the June 2000 Kleifarvatn earthquake, Iceland. Geophysical Journal International, 176(2), 389?404. https://doi.org/10.1111/j.1365- 246x.2008.03989.x Sudhaus, H., & J?nsson, S. (2011). Source model for the 1997 Zirkuh earthquake (MW= 7.2) in Iran derived from Jers and ERS Insar Observations. Geophysical Journal International, 185(2), 676?692. https://doi.org/10.1111/j.1365-246x.2011.04973.x Teng, L. S., C. T. Lee, C. H. Peng, J. J. Chu, and W. F. Chen, (2001). Origin and geological evolution of the Taipei Basin, northern Taiwan. West. Pac. Earth Sci., 1, 115-142 Tong, X., Sandwell, D. T., & Smith-Konter, B. (2013). High-resolution interseismic velocity data along the San Andreas Fault from GPS and InSAR. Journal of Geophysical Research: Solid Earth, 118(1), 369-389. doi:10.1029/2012jb009442 USGS. (2021). Earthquake Catalog. Retrieved July 30, 2021, from https://earthquake.usgs.gov/earthquakes Wang, J.-H., Chen, K.-C., Leu, P.-L., & Chang, J.-H. (2015). B-values observations in Taiwan: A Review. Terrestrial, Atmospheric and Oceanic Sciences, 26(5), 475?492. https://doi.org/10.3319/tao.2015.04.28.01(t) Wang, Y.-J., Chan, C.-H., Lee, Y.-T., Ma, K.-F., Shyu, J. B. H., Rau, R.-J., & Cheng, C.-T. (2016). Probabilistic seismic hazard assessment for Taiwan. Terrestrial, Atmospheric and Oceanic Sciences, 27(3), 325-. doi:10.3319/TAO.2016.05.03.01(TEM) 96 Watson, A. R., Elliott, J. R., & Walters, R. J. (2022). Interseismic strain accumulation across the main recent fault, SW Iran, from Sentinel-1 Insar Observations. Journal of Geophysical Research: Solid Earth, 127(2). https://doi.org/10.1029/2021jb022674 Weiss, J. R., Walters, R. J., Morishita, Y., Wright, T. J., Lazecky, M., Wang, H., . . . Parsons, B. (2020). High?Resolution Surface Velocities and Strain for Anatolia From Sentinel?1 InSAR and GNSS Data. Geophysical Research Letters, 47(17). doi:10.1029/2020gl087376 Yu, S.-B., Chen, H.-Y., & Kuo, L.-C. (1997). Velocity field of GPS stations in the Taiwan area. Tectonophysics, 274(1- 3), 41?59. https://doi.org/10.1016/S0040-1951(96)00297-1 Yunjun, Z., Fattahi, H., & Amelung, F. (2019). Small baseline InSAR time series analysis: Unwrapping error correction and noise reduction. Computers and Geosciences, 133, [104331]. https://doi.org/10.1016/j.cageo.2019.104331 Zhao, Dezheng, et al. (2020). Multifault Complex Rupture and Afterslip Associated with the 2018 MW?6.4 Hualien Earthquake in Northeastern Taiwan. Geophysical Journal International, vol. 224, no. 1, pp. 416?434., https://doi.org/10.1093/gji/ggaa474. Zhao, B., B?rgmann, R., Wang, D., Tan, K., Du, R., & Zhang, R. (2017). Dominant controls of downdip afterslip and viscous relaxation on the postseismic displacements following the Mw7.9 Gorkha, Nepal, earthquake. Journal of Geophysical Research: Solid Earth, 122. https://doi.org/10.1002/2017JB014366 97