# Monitoring interseismic activity on the Ilan Plain (NE Taiwan) using Small Baseline PS-InSAR, GPS and leveling measurements: partitioning from arc-continent collision and backarc extension

Monitoring interseismic activity on the Ilan Plain (NE Taiwan) using Small Baseline PS-InSAR, GPS... Abstract The Ilan Plain, located in Northeast Taiwan, represents a transition zone between oblique collision (between the Luzon Arc and the Eurasian Plate) and backarc extension (the Okinawa Trough). The mechanism for this abrupt transition from arc-continent collision to backarc extension remains uncertain. We used Global Positioning System (GPS), leveling and multi-interferogram Small Baseline Persistent Scatterer Interferometry (SBAS-PSI) data to monitor the interseismic activity in the basin. A common reference site was selected for the data sets. The horizontal component of GPS and the vertical measurements of the leveling data were converted to line-of-sight (LOS) data and compared with the SBAS-PSI data. The comparison shows that the entire Ilan Plain is undergoing rapid subsidence at a maximum rate of −11 ± 2 mm yr−1 in the LOS direction. We speculate that vertical deformation and anthropogenic activity may play important roles in this deformation. We also performed a joint inversion modeling that combined both the DInSAR and strong motion data to constrain the source model of the 2005 Ilan earthquake. The best-fitting model predicts that the Sansing fault caused the 2005 Ilan earthquake. The observed transtensional deformation is dominated by the normal faulting with a minor left-lateral strike-slip motion. We compared our SBAS-PSI results with the short-term (2005–2009) groundwater level changes. The results indicate that although pumping-induced surface subsidence cannot be excluded, tectonic deformation, including rapid southward movement of the Ryukyu arc and backarc extension of the Okinawa Trough, characterizes the opening of the Ilan Plain. Furthermore, a series of normal and left-lateral strike-slip transtensional faults, including the Choshui and Sansing faults, form a bookshelf-like structure that accommodates the extension of the plain. Although situated in a region of complex structural interactions, the Ilan Plain is primarily controlled by extension rather than by shortening. As the massive, pre-existing Philippines–Ryukyu island arc was pierced by the Philippine Sea Plate, the Ilan Plain formed as a remnant backarc basin on the northeastern corner of Taiwan. Radar interferometry, Asia, Seismicity and tectonics 1 INTRODUCTION Arc-continent collision and backarc extension are entirely opposite kinematic phenomena in classical orogenic theory. Arc-continent collisions cause mountain building via compressional structures, which typically occur within suture belts that form between converging plates. For example, Taiwan formed from the collision between the Philippine Sea Plate and the Chinese continental shelf (e.g. Huang et al.1997; Lee et al.2015). In contrast, backarc extension represents tectonic collapse that results from the relaxation of tension, the opposite of arc-continent collision. Backarc extension frequently occurs in interior continental areas behind collisional margins. For example, backarc basins are present in the central and western Mediterranean (e.g. Carminati et al.1998, 2012). Northeast Taiwan is one of the few areas on Earth that is currently experiencing two different stress regimes (Hu et al.1996, 2002). On the one hand, the Philippine Sea Plate is obliquely colliding with the southeastern margin of the Chinese continental shelf at a rate of 82 mm yr−1 with respect to the Eurasia Plate (Fig. 1, Yu et al.1997, 1999; Lin et al.2010; Tsai et al.2015). This convergence not only causes crustal thickening along several main mountain ranges in Taiwan (e.g. the Hsüehshan Range, Clark et al.1993; the Central Range, Hovius et al.2000) but also pushes the entire island northwestward (GPS monitoring, e.g. Yu et al.1997; Lin et al.2010) toward the interior of Eurasia (Suppe 1984, 1987; Teng 1990; Sibuet & Hsu 2004). Northeast Taiwan represents an exception to these trends in that this region is experiencing rapid southeastward extrusion (Angelier et al.2009; Hou et al.2009) along with the Ryukyu Islands (Fig. 1; Yu et al.1999; Nakamura 2004). Figure 1. View largeDownload slide Tectonic framework of Taiwan and the surrounding area. The Philippine Sea Plate is moving northwestward (directed N304.1°W) at a rate of 78 mm yr−1 in the MORVEL plate model (DeMets et al.2010). The Ryukyu arc-trench system is retreating to the south with a migration velocity of ∼30–40 mm yr−1. The GPS vectors are referenced to the fixed Chinese continental margin (permanent fixed GPS site SHAO in Shanghai, China; GPS data from Yu et al.1999; Nakamura 2004). LS: Lishan fault and LV: Longitudinal Valley fault. Figure 1. View largeDownload slide Tectonic framework of Taiwan and the surrounding area. The Philippine Sea Plate is moving northwestward (directed N304.1°W) at a rate of 78 mm yr−1 in the MORVEL plate model (DeMets et al.2010). The Ryukyu arc-trench system is retreating to the south with a migration velocity of ∼30–40 mm yr−1. The GPS vectors are referenced to the fixed Chinese continental margin (permanent fixed GPS site SHAO in Shanghai, China; GPS data from Yu et al.1999; Nakamura 2004). LS: Lishan fault and LV: Longitudinal Valley fault. One possible cause of this transition in Northeast Taiwan is a subduction polarity reversal of the Philippine Sea Plate (e.g. Kim et al.2005; Huang et al.2012; Ustaszewski et al.2012). The northern margin of the Philippine Sea Plate is subducting northward beneath the oceanic crust of the Okinawa Trough (Kuo-Chen et al.2015). Backarc basins typically form continent ward of plate collisions, for example, the Okinawa backarc basin (e.g. Sibuet et al.1987; Sibuet & Hsu 1997; Wang et al.1999; Shinjo & Kato 2000). To the south, however, the southwestern margin of the Philippine Sea Plate is obducting onto and overriding the South China Sea basin along the Manila trench (e.g. Seno 1977; Kato et al.2003; Ku & Hsu 2009; Wang & Bilek 2014; Hsiung et al.2015). Consequently, two entirely different tectonic domains have formed within a short distance (<100 km) of this margin in the vicinity of Northeast Taiwan. To the east and south, these domains consist of an extensional regime associated with the westernmost edge of the Okinawa backarc extension and coupled to the rapid southward rollback of the Ryukyu trench. A compressional domain to the west is represented by crustal shortening and thrusting along a major plate boundary and/or reverse faults, represented by the Lishan (LS) fault (Lee et al.1997), which exhibits overthrusting deformation behaviour (Kuo-Chen et al.2015). Compared to other deformation phenomena, interseismic strain accumulation is small in magnitude and is distributed over large spatial distances and temporal spans. The InSAR technique may have the potential to produce regional-scale deformation maps that can identify all active structures and measure the rate of strain accumulation across these active structures. However, short-term monitoring of these tiny but seismically active faults is not easy, especially for active intracontinental structures. Additionally, the expected low-amplitude deformations are comparable in scale to contamination by other errors, such as orbital ramps, temporal variations in the stratified troposphere (e.g. Li et al.2009; Jolivet et al.2011) and topographic effects. Quantitative analysis of interseismic deformation is imperative to understanding the tectonics of an area. The complementary strengths of InSAR, GPS and leveling measurements allow interseismic displacements to be inferred by appropriately combining measurements. Although the displacement fields of the Ilan Plain have been constructed from Permanent Scatterers SAR interferometry (PS-InSAR, e.g. Kang et al.2015) or GPS data (e.g. Angelier et al.2009; Hou et al.2009), a complete and spatially continuous interseismic deformation field associated with the neotectonic activity of the basin is not yet available. InSAR measurements are spatially continuous but are sensitive to only the surface displacement in the line-of-sight (LOS) direction of the radar (e.g. Wright et al.2004). GPS observations provide 3-D velocity vectors but are available only in limited quantities. Leveling measurements face the same problem. Additionally, these three data sets have different reference frames; for instance, leveling (Yang et al.2003; Chen et al.2011) and GPS measurements are relative to a tidal gauge in Keelung along the coast of North Taiwan and the S01R station on Paisha Island (e.g. Angelier et al.2009; Hou et al.2009) and the LOS mean rate is referenced to a ground-to-satellite radar system (Kang et al.2015). In this paper, we monitor interseismic deformation in the Ilan Plain by appropriately integrating InSAR, leveling and GPS measurements in the area. The different reference frames of these techniques are united into a common reference frame. A newly developed method called the Full Rank Matrix Small Baseline Subset InSAR processor (FRAM-SBAS, Biggs et al.2007; Li et al.2009) is applied to produce Small Baseline Persistent Scatterer Interferometry (SBAS-PSI). We compare the horizontal movement based on the GPS data and the vertical movement from the leveling with the movement tendencies of Persistent Scatterers (PS). We thereby characterize which patterns of deformation control the opening of the Ilan Plain, especially focusing on possible subsurface structures and the associated accommodation mechanism between arc-continent collision and backarc extension. 2 SEISMOTECTONIC SETTINGS The Ilan Plain of Northeast Taiwan is situated between two main mountain ranges on Taiwan Island: the Hsüehshan Range to the north and the Backbone Range to the south. Two seismic zones surround the plain. An area known as the Okinawa seismic belt occupies the area between the eastern coast of the Ilan Plain and the Okinawa Trough (seismic zone I; Fig. 2, Kao & Rau 1999). This seismic zone includes a series of normal faults and grabens that generally strike in an east–west direction (Kao & Rau 1999; Wang et al.2000; Huang et al.2012). The second seismic zone, known as the ‘Collision seismic zone’, lies just northeast offshore of Hualien County, Taiwan (Lo et al.2013), and forms an elbow-shaped concentration of earthquakes that occur at depths of approximately 20–30 and 30–50 km (yellow and red points in Fig. 2). The colliding seismic zone represents the location where the reversal of subduction polarity begins. The subduction of the Philippine Sea Plate shifts from northward to northwestward at this hinge point (Huang et al.2012). Figure 2. View largeDownload slide Seismic activity in Northeast Taiwan. The beach balls represent the focal mechanism solutions for earthquake events (red). All the focal mechanism solutions are calculated from the Broadband Array in Taiwan for Seismology (BATS) catalogue. The background seismicity is colour coded by focal depth (the events are from the BATS catalogue). The light grey lines represent the boundary of lithology. The black dotted lines ‘w1’ (Shyu et al.2005) and ‘w2’ (Lai et al.2009; Kang et al.2015) represent two possible northward projections of the Lishan fault. (I) Okinawa Seismic Zone and (II) Collision Seismic Zone. Figure 2. View largeDownload slide Seismic activity in Northeast Taiwan. The beach balls represent the focal mechanism solutions for earthquake events (red). All the focal mechanism solutions are calculated from the Broadband Array in Taiwan for Seismology (BATS) catalogue. The background seismicity is colour coded by focal depth (the events are from the BATS catalogue). The light grey lines represent the boundary of lithology. The black dotted lines ‘w1’ (Shyu et al.2005) and ‘w2’ (Lai et al.2009; Kang et al.2015) represent two possible northward projections of the Lishan fault. (I) Okinawa Seismic Zone and (II) Collision Seismic Zone. In addition to the Longitudinal Valley fault, which acts as a plate boundary between the Philippine Sea Plate and the Eurasian Plate (e.g. Yu & Kuo 2001; Chen et al.2007), another large-scale boundary fault, the LS fault, separates the Hsüehshan Range from the Backbone Range to the west of the Ilan Plain. The LS fault is a steeply dipping reverse fault with dip directions that change from eastward to westward along strike (Lee et al.1997). This structure separates two major geological structural units within the Taiwan orogeny: the Mio-Oligocene Daliao Group, which mainly occurs in the Hsüehshan Range, and the Oligocene Chiayang Group, which mainly occurs in the Backbone Range (Fig. 3, Teng 1990; Lee et al.1997). How far the LS fault extends to the north is a matter of debate. Some researchers have argued that the LS fault extends northwest along the northwestern margin of the Ilan basin (route ‘w1’ in Fig. 2; Shyu et al.2005). The fault becomes a normal fault upon entering the Ilan Plain (Shyu et al.2005). Other models have stated that the LS fault bends to the south and may connect to the east–west-striking Choshui fault at depth when it enters the Ilan Plain (route ‘w2’ in Fig. 2; Kang et al.2015). Figure 3. View largeDownload slide (a) Geological map (Ho 1986, 1988) and cross-sections of the Ilan Plain. The thin black contour lines represent the thickness of the fluvial sediments. (b) Lines 1 and 2 show transects based on exploratory seismic profiles, as reported in Chiang (1976). CS: Choshui fault; HZ: Huzing fault; IL: Ilan fault; JX: Jiaoxi fault; KS: Kengssu fault; NT: Niuto fault; TC: Toucheng fault; SS: Sansing fault and YC: Yaichein fault. Figure 3. View largeDownload slide (a) Geological map (Ho 1986, 1988) and cross-sections of the Ilan Plain. The thin black contour lines represent the thickness of the fluvial sediments. (b) Lines 1 and 2 show transects based on exploratory seismic profiles, as reported in Chiang (1976). CS: Choshui fault; HZ: Huzing fault; IL: Ilan fault; JX: Jiaoxi fault; KS: Kengssu fault; NT: Niuto fault; TC: Toucheng fault; SS: Sansing fault and YC: Yaichein fault. The Ilan Plain is a typical equilateral triangular basin, with its eastern coast facing the West Pacific Ocean. The basin is filled with 500–2000 m of fluvial sediments (Liu 1995; Dadson et al.2003). Consequently, no active structures can be observed on the surface. In terms of sedimentation rate, a previous study (Lin et al.2004) based on one well (the Wuyuan site) south of the Nanyang river suggested that the region featured a very rapid sedimentation rate of 21.7 mm yr−1 between 3400 and 3100 yr B.P., a slow sedimentation rate of ∼1.7 mm yr−1 during the next 2000 yr and an intermediate sedimentation rate of 6.9 mm yr−1 during the last 1000 yr. Two previous studies (Rau et al.2008; Ching et al.2011) proposed high slip rates (>20 mm yr−1) along the two boundary faults (Toucheng and Choshui) based on block models using campaign GPS data. However, possible seismogenic and buried structures were not taken into account in their models. Thus, these results may overestimate the seismic potential of these two boundary faults. The basement of the plain consists of consolidated basal layers that are overlain by shallow, uncemented fluvial sediments (Liu 1995). Previous seismic exploration transects revealed that the basement has an asymmetric washbowl-like shape (Fig. 3a, Chiang 1976). The thickest sediment occurs in an inlet where a river flows into the western Pacific (the 1600 m basement contour in Fig. 3a). However, the current subsidence centre seems to have shifted to the southern bank of the Lanyang River (Kang et al.2015). Chiang (1976) identified a series of blind faults in the basement. From north to south, these faults are the Toucheng (TC), Jiaoxi (JX), Niuto (NT), Huzing (HZ), Yaichein (YC), Ilan (IL), Kengssu (KS), Choshui (CS) and Sansing (SS) faults (fault names are based on Kang et al.2015). Chiang (1976) interpreted these faults as high-angle thrust faults (Fig. 3b) according to the stratal contact relationships but could not evaluate their recent activity. Kang et al. (2015) applied PS-InSAR technology to analyse the relationships with the basement faults and regional geodynamic setting. These authors revealed that these reverse faults could change into normal faults because the plain is experiencing subsidence. Lai et al. (2009) inverted the coseismic rupture process of the 2005 Ilan earthquake doublet with GPS and seismicity data sets. One 15-km-long fault was recognized between the Choshui and Sansing faults and was considered to be the seismogenic structure responsible for the 2005 Ilan earthquake doublet. This active structure exhibits normal faulting with a left-lateral sense of shear. 3 GPS AND LEVELING MEASUREMENTS In terms of GPS and leveling measurements, Taiwan is one of the most monitored places on Earth. The Ilan Plain contains eight permanent GPS stations (provided by Lin et al.2010), 29 temporary GPS monitoring points (provided by the Central Geological Survey) and 99 leveling points from four leveling lines in the leveling measurement project known as the 2001 Taiwan Vertical Datum (Yang et al.2003; Chen et al.2011). The time spans of the GPS and leveling measurements that we used in this study differ: the GPS data were collected from 2004 to 2013, and the leveling data were collected from 2000 to 2008 (Chen et al.2011). However, both types of data were converted to average velocities to determine the relatively stable long-term movement trends, for example, on a 10-yr scale. 3.1 Horizontal and vertical displacement fields We used GPS data collected from 2004 to 2013. Some workers have suggested that the Ilan Plain experiences an escape movement to the southeast (e.g. Angelier et al.2009; Hou et al.2009). The plain generally extrudes towards S139°E (Fig. 4a). A maximum velocity of ∼44 mm yr−1 occurs around the Su’ao coast in the southeastern corner of the plain. The velocity vectors in the northern half of the plain are generally slower than those in the southern half. The northern velocities range from 0.6 to 5.9 mm yr−1 towards 15°–160°, and the southern velocities range from 2.7 to 44 mm yr−1 towards 60°–139° (Hou et al.2009). Figure 4. View largeDownload slide (a) Horizontal velocity field relative to the Paisha station, S01R, which is located on Penghu Island along the Chinese continental margin (2004–2013). The white squares mark permanent GPS stations, and the black triangles markcampaign-mode GPS arrays. The velocity vectors are scaled according to a 95 per cent confidence interval. (b) Vertical velocity field (2000–2008) for the Ilan Plain according to leveling data. The elevation was calculated relative to the Keelung tidal gauge on Taiwan's northern coast (star in inset). The elevation estimates have a mean standard deviation of ± 1.64 mm yr−1 (Chen et al.2011). Figure 4. View largeDownload slide (a) Horizontal velocity field relative to the Paisha station, S01R, which is located on Penghu Island along the Chinese continental margin (2004–2013). The white squares mark permanent GPS stations, and the black triangles markcampaign-mode GPS arrays. The velocity vectors are scaled according to a 95 per cent confidence interval. (b) Vertical velocity field (2000–2008) for the Ilan Plain according to leveling data. The elevation was calculated relative to the Keelung tidal gauge on Taiwan's northern coast (star in inset). The elevation estimates have a mean standard deviation of ± 1.64 mm yr−1 (Chen et al.2011). The leveling results show that the Ilan Plain is presently experiencing rapid subsidence. According to the two E-W-trending profiles that transect the plain and are roughly divided by the 121°36΄E meridian, an obvious velocity shift appears between mountainous areas and the relatively flat plain. The vertical vectors increase from north to south along the N–S-striking survey transects. Several urban areas (e.g. Ilan County and the town of Luotung) show rapid subsidence. Near the Backbone Range, both the southern boundary of the basin and the Su’ao coast exhibit uplift. 3.2 Principal dilation and shearing strain rate In contrast to velocity data (e.g. Fig. 4), the strain rate tensor is independent of the reference frame. However, the maximum shearing and dilation strain rate from GPS data sets could reveal crustal deformation for an active fault during interseismic periods and its possible connection to seismic hazard potential (Ward 1994). We derived strain rates from the velocity estimates according to the procedure described by Shen et al. (1996). This method, which follows previous methods (e.g. Frank 1966; Prescott 1976), interpolates the strain rates from local geodetic measurements through continuous functions within the entire network (for details, see Lin et al.2010). The density of the GPS sites in this study was insufficient; therefore, we created triangles and calculated their centre strain rate. Although this method also applies interpolation method, however, the calculated strain rate is relatively more accurate, because the strain rate in the centre is determined based only on neighbouring real GPS sites. Our results indicate that NW–SE extension is occurring widely in the Ilan Plain (black arrows in Fig. 5a). The largest extensional rate is located on the southern bank of the Langyang River, with a strain rate of 2.66 μstrain yr−1 (μ stands for 10e−6; numbered (1) in Fig. 5a) in the NW–SE direction and a shortening rate of 1.28 μstrain yr−1 (numbered (2) in Fig. 5a) in the NE–SW direction (azimuths 132°–150°; Lin et al.2010). This extensional dilation strain rate pattern is consistent with the description of the coseismic displacement field of the 2005 Ilan earthquake doublet (Lai et al.2009) and the opening direction range of the plain (e.g. Yu et al.1997). Furthermore, this pattern indicates clockwise rotation, as shown by the variation in the direction of the extensional axis from north to south (Fig. 5a), which is also consistent with the opening direction of the Ilan Plain. Figure 5. View largeDownload slide (a) Horizontal principal axes of the strain rates. The amount and direction of the principal strain rates are shown by red (shortening) and black (extension) arrows. The green triangles are the locations of the GPS sites. The Ilan, Choshui and Sansing faults are three major active normal faults inferred from the maximum principal strain rate. Their locations are shifted to the right (dotted lines). From north to south, the maximum extensional axis seems to rotate in a clockwise manner. (b) Left-lateral shearing axes of the strain rates. The amount and direction of the shear strain rates are represented by the gold-coloured bold lines. The above three normal faults exhibit left-lateral senses of shear Figure 5. View largeDownload slide (a) Horizontal principal axes of the strain rates. The amount and direction of the principal strain rates are shown by red (shortening) and black (extension) arrows. The green triangles are the locations of the GPS sites. The Ilan, Choshui and Sansing faults are three major active normal faults inferred from the maximum principal strain rate. Their locations are shifted to the right (dotted lines). From north to south, the maximum extensional axis seems to rotate in a clockwise manner. (b) Left-lateral shearing axes of the strain rates. The amount and direction of the shear strain rates are represented by the gold-coloured bold lines. The above three normal faults exhibit left-lateral senses of shear Fig. 5(b) shows the distribution of the left-lateral shear strain rate in the Ilan Plain. We identified a series of E–W-striking discontinuous shear strain rates in the plain, indicating that some faults that underlie the plain may exhibit a left-lateral sense of shear. Considering the fault identified in seismic reflection profiles (Chiang 1976), we propose that lines (1)–(3) may correspond with the Ilan, Choshui and Sansing faults, respectively. A maximum sinistral strain rate of ∼0.7 μstrain yr−1 appears along the eastern segment of the Sansing fault (number (3) in Fig. 5b). Likewise, the epicentre locations of the 2005 Ilan earthquake doublet (Lai et al.2009) are consistent with the Sansing fault trace. 4 SBAS-PSI MONITORING 4.1 Methodology GPS and leveling data are highly accurate in terms of their respective horizontal and vertical components. However, these data are inherently point based. Current SAR interferometry analysis techniques can provide high-density measurements for monitoring land surface deformation. Among these techniques, such as stacking (e.g. Wright et al.2001; π-RATE, Biggs et al.2007; Elliott et al.2008; Wang et al.2009; Wang & Wright 2012), PS-InSAR technology (Ferretti et al.2001; Hooper et al.2012) and SBASs (Berardino et al.2002), PS-InSAR uses the movement of pixels or points that are coherent throughout the data set to determine the deformation in a particular area. StaMPS is a relatively mature platform that integrates commonly used approaches, such as PS-InSAR and SBAS algorithms (e.g. Hooper et al.2012). Indeed, StaMPS makes no a priori assumptions. However, this method does not have the ability to reduce noise from atmospheric delay, orbital ramps and topography. Furthermore, this method must resolve phase ambiguities in interferogram stacking either by searching a pre-defined solution space or by sparse phase-unwrapping methods. The efficiency and success of phase unwrapping cannot be guaranteed. Finally, the StaMPS software suite works well only in dry, rocky and urban areas. The rainy, wet and cloudy Taiwan Island would definitely introduce more atmosphere noise. In this study, we used a newly developed multitemporal SBAS-PSI method to detect the interseismic deformation in the Ilan Plain. This method is an effective tool for reducing atmospheric artefacts (e.g. Remy et al.2003; Cavalié et al.2008; Elliott et al.2008; Li et al.2009; Jolivet et al.2011, 2012, 2013, 2014) and provides more precise deformation signals (e.g. Biggs et al.2007). We applied the FRAM-SBAS method. FRAM-SBAS is a software suite that combines the open-source code of ROI_PAC (v3.1; Rosen et al.2004), The Generic Mapping Tools (Wessel & Smith 1998; Wessel et al.2013) and an inversion algorithm from Li et al. (2009). We constructed a master–slave network or small baseline image pair under the following conditions to minimize the spatial and temporal decorrelation. First, the perpendicular baselines of master–slave pairs were set to be smaller than 500 m. Second, each acquisition node in the network was given at least four link pairs, meaning that each node has a minimum of four connections with other nodes. Finally, we achieved 170 master–slave pairs (Fig. 6c). Coherent points were extracted based on coherent images; the points that were selected as PS shad coherent values that were greater than 0.38 with a standard deviation of <2.0 mm yr−1. We used the 3-D spatio-temporal approach that was described by Hooper (2008) to unwrap the phase measurements. Figure 6. View largeDownload slide (a) Coverage area of Track 461 for the Ilan Plain (yellow colour in map). (b) A total of 170 pairs were produced under the condition of pre-defined small perpendicular baselines. Figure 6. View largeDownload slide (a) Coverage area of Track 461 for the Ilan Plain (yellow colour in map). (b) A total of 170 pairs were produced under the condition of pre-defined small perpendicular baselines. 4.2 Envisat/ASAR images We used 32 C-band Advanced Synthetic Aperture Radar (ASAR) radar images that were acquired between 2003 November 22, and 2009 June 13, from Environmental Satellite (Envisat) on descending Track 461 (Fig. 6a). All the radar images covered most of the study area (Fig. 6a). The time intervals between the images from 2003 and 2004 and the images from 2008 and 2009 exceeded one year (Fig. 6b). The topographic phase was removed using the Shuttle Radar Topography Mission 3-arcsec Digital Elevation Model (DEM; Farr et al.2007). 5 INTERSEISMIC SURFACE DISPLACEMENT FIELDS IN THE ILAN PLAIN 5.1 Small baseline PSI results We obtained a total of 39 881 high-quality SBAS-PSI pixels from 170 interferogram pairs. All the standard deviations of the PSs were <2.0 mm yr−1. Because the study area includes many densely populated areas, the average data density reached approximately 100 PSs km−2. The LOS mean displacement rate in the Ilan Plain was limited to a range of −11 to 7.5 mm yr−1 (Fig. 7a). A more detailed analysis is provided in Section 5.3. The uncertainties were larger in mountainous regions and in the mountain–plain transition zones (Fig. 7b). We speculate that the large transition from uplift to subsidence and the associated unwrapping error are responsible. Figure 7. View largeDownload slide (a) SBAS-PSI velocity field along the LOS for the Ilan Plain relative to the common reference YILN (indicated by the asterisk). The squares (permanent GPS stations) and diamonds (campaign GPS sites) display the comparison results between the GPS and SBAS-PSI data. The search radius for the PSs was set to 2 km surrounding each GPS station. (b) Estimated standard deviation of the LOS velocity. Figure 7. View largeDownload slide (a) SBAS-PSI velocity field along the LOS for the Ilan Plain relative to the common reference YILN (indicated by the asterisk). The squares (permanent GPS stations) and diamonds (campaign GPS sites) display the comparison results between the GPS and SBAS-PSI data. The search radius for the PSs was set to 2 km surrounding each GPS station. (b) Estimated standard deviation of the LOS velocity. 5.2 Comparison with GPS data GPS, leveling and SBAS-PSI monitoring occurs in different reference frames. The GPS data were relative to the Paisha station (S01R), a permanent GPS site on Penghu Island. The leveling data were referenced to the Keelung tidal gauge on the northern coast of Taiwan. The SBAS-PSI velocity field was referenced to a ground-satellite system. We selected a common reference site to compare these three data sets. Generally, the vertical motion of the GPS measurement is not precise than that of the leveling. Thereafter, when we do the conversion from the eastward, northward, vertical velocity onto the LOS direction, we prefer to select the vertical motion provided by a leveling site that is much closed to a permanent/campaign GPS or under our searched radius. If there will not exist such a leveling site or exceed our searching radius, we have to use the vertical measurement of the GPS site. A permanent GPS station situated to the south of the Ilan County, the YILN was selected (snow sign in Fig. 7a). We achieved a total of 3069 PSs that could fall into a searching circle (2 km as the searching radius). Furthermore, the comparison result is improved by using the new common reference. This site also contains a leveling surveying site, numbered 9039. We unified the references for all three data sets. The eastward, northward and vertical motions were converted into the LOS direction by multiplying them by a full unit vector (Fialko et al.2001). Generally, the vertical motion of the GPS measurement is as not precise as that of the leveling measurements. Hence, during the conversion of the eastward, northward, and vertical velocities into the LOS direction, we used the vertical motion from a leveling site if one was close to the permanent/campaign GPS. If no such leveling site existed or exceeded our search radius, we used the vertical measurement from the GPS site. The descending Envisat/ASAR satellite flies with a mean azimuth angle of 12.41° and with a mean incidence angle of 23°, resulting in an estimated unit vector of (0.3724, −0.0880, 0.9239). The SBAS-PSI density is much higher than that of the GPS stations; hence, the average SBAS-PSI LOS velocity was estimated around each GPS station. We searched all the PSs within a certain distance from the common reference site. Then, we calculated their mean LOS displacement velocities. Here, the search distance is 2 km and the average LOS velocity of the YILN is −0.73 mm yr−1. At the same time, the converted eastward, northward and vertical component of the GPS data into the LOS direction is 0.01 mm yr−1. Thus, all the LOS velocities could be shifted into the common reference frame based on the YILN site by adding 0.74 mm yr−1. We compared all the campaign (displayed by diamonds) and permanent (represented by squares) GPS rate conversion results (triangles) with the PS pixels (coloured circles). This process rested on the assumption that the precision of the vertical GPS measurements was acceptable. The comparison results were satisfactory in the centre of the basin but became worse when approaching the boundary of the Ilan Plain, which may be related to the fact that the campaign and permanent GPS measurements were calculated from two surveying approaches with different measuring cycle periods, precisions, consistencies and stabilities. 5.3 Interseismic activity of the Ilan Plain The 1-D LOS measurement is a composite vector because both horizontal and vertical motions contribute to it. Therefore, the method that combines InSAR data with spatially sparse GPS observations has been widely applied in distinguishing the individual contributions from the horizontal and vertical components, for example, isolating the vertical rates associated with non-tectonic processes (e.g. groundwater and oil pumping, Bawden et al.2001) and/or extracting possible tectonic signals (e.g. significant uplift in the San Francisco Bay Area, Bürgmann et al.2006), by removing the contributions of the GPS-derived horizontal velocity field from the InSAR range-change rates. The Ilan Plain is experiencing not only horizontal extrusion but also rapid subsidence. Vertical deformation is always the result of both non-tectonic (including landslides, and seasonal and long-term groundwater level changes) and tectonic deformations (sediment settling, normal faulting and time-dependent earthquake cycle deformation). Thus, when the horizontal components of the GPS data (Ve and Vn) and the vertical measurements of the leveling data (Vu) are converted into the LOS direction (LOSconv) for comparison with the SBAS-PSI measurements, we are able to constrain the relative contributions of the southeastward extrusion and vertical motions in the opening process of the Ilan Plain. The horizontal components should represent purely tectonic deformation (e.g. the lateral extrusion of the block). Unfortunately, the vertical signals could not be easily isolated because they represent a combination of tectonic movement and the effects of anthropogenic activities, such as water pumping. We use the groundwater data from three major wells (ZX, WJ and LZ) to characterize how the long-term deformation of the Ilan Plain controlled by tectonic or non-tectonic signals (Fig. S1, Supporting Information). Our analysis suggests that surface elevation change is the combined effect of multiple aquifers at the same location; thus, pumping-induced surface subsidence cannot be excluded (Fig. S2, Supporting Information). This result is different from the suggestion proposed by Kang et al. (2015), the detailed discussion of pumping-induced deformation in the Ilan Plain is listed in Supporting Information Appendix S1. We multiplied the horizontal (Ve and Vn) and vertical (Vu) measurements by unit vectors of (0.3724, −0.0880, 0) and (0, 0, 0.9239), respectively. The sum of the vertical and horizontal conversion results should theoretically be equal to the average displacement velocity of the PSs (eq. 1).   $${\rm{LO}}{{\rm{S}}_{{\rm{conv}}}} = \Delta ({\rm{V}}{{\rm{u}}_{{\rm{conv}}}}) + \Delta ({\rm{V}}{{\rm{h}}_{{\rm{conv}}}}) = {\rm{LO}}{{\rm{S}}_{{\rm{SBAS}}}}$$ (1)  $${\rm{RMS}} = {\rm{LO}}{{\rm{S}}_{{\rm{SABS}}}} - {\rm{LO}}{{\rm{S}}_{{\rm{conv}}}}$$ (2) where Δ(Vuconv) and Δ(Vhconv) are the LOS conversion results from the horizontal and vertical components, respectively; LOSconv and LOSSBAS are the sum of the LOS conversion results from the horizontal and vertical components and the SBAS-PSI mean rate change, respectively and RMS is the residual error of LOSSBAS and LOSconv. 5.3.1 Boundary normal faults—Zones (1) and (6) The Ilan triangular basin is bounded by two boundary faults located along the southern slope of the Hsüehshan Range and the southern margin of the plain. These faults are typical normal faults, as inferred from a series of fault scarps observed in the field (e.g. Shyu et al.2005; Kang et al.2015). In our LOS mean displacement rate, these two flanking faults display decreasing LOS rates with respect to the YILN site (Fig. 8a). Warm-toned PS points appear in both the hanging and footwalls of the Hsüehshan Range Front (HSRF) and Backbone Range Front faults (BBRF), which are coloured orange and yellow, respectively (see Fig. 7a). Figure 8. View largeDownload slide (a) SBAS-PSI LOS velocity profile that crosses the Ilan Plain from north to south. The projected SBAS-PSI LOS velocities are represented by blue circles, yellow (campaign-mode GPS array) and blue diamonds (permanent GPS array), and magenta triangles (leveling stations). The blue line represents the topographic relief along a profile based on 40-m DEM data in Taiwan. The GPS and leveling data were collected from within a 5-km-wide swath along the cross-section. The SBAS PS pixels were collected from within a 2-km-wide swath. HSRF: Hsüehshan Range Front fault and BBRF: Backbone Range Front fault. (b) Velocity variation trend. The vertical thin black lines represent a set of possible faults. IL: Ilan fault; CS: Choshui fault and SS: Sansing fault. Figure 8. View largeDownload slide (a) SBAS-PSI LOS velocity profile that crosses the Ilan Plain from north to south. The projected SBAS-PSI LOS velocities are represented by blue circles, yellow (campaign-mode GPS array) and blue diamonds (permanent GPS array), and magenta triangles (leveling stations). The blue line represents the topographic relief along a profile based on 40-m DEM data in Taiwan. The GPS and leveling data were collected from within a 5-km-wide swath along the cross-section. The SBAS PS pixels were collected from within a 2-km-wide swath. HSRF: Hsüehshan Range Front fault and BBRF: Backbone Range Front fault. (b) Velocity variation trend. The vertical thin black lines represent a set of possible faults. IL: Ilan fault; CS: Choshui fault and SS: Sansing fault. Zone (1) exhibits similar movement trends in the leveling (e.g. 9033–9035) and PSI data sets (Fig. 8a). However, the density of GPS data is quite low in this region. We could not determine the deformation from the horizontal components. The vertical velocity decreased by ∼2 mm yr−1 from the Hsüehshan Range to the Ilan Plain (Fig. 8b). We speculate that the HSRF is a normal fault that dips to the south. Relative tothe YILN site, the region to the south of the BBRF appears to have an LOS rate that is ∼8.5 mm yr−1 faster than the area to the north (Fig. 8b). However, this tendency could not be verified by any other geodetic data set because no GPS or leveling survey points exist in the Backbone Range. Although the BBRF is generally considered to be a normal fault, we attribute this shortening tendency to the horizontal components because the GPS velocities indicate that the southern corner of the Ilan Plain is moving laterally to the southeast. Meanwhile, the radar satellites monitor the land deformation with a right-side-looking pattern, facing to the west, which shortens the baseline. 5.3.2 Interseismic activity in the Ilan Plain—Zones (2)–(5) Because the Ilan basin contains many large population centres, the tectonic deformation signals are likely juxtaposed with anthropogenic ground subsidence, for example, the massive pumping of groundwater. Therefore, we could not directly correlate every appearance of rate gradients to an active structure. Below, we describe several zones of movement. In Zone (2), velocity gradients of ∼5 and ∼6 mm yr−1 were measured in the SBAS-PSI and leveling data sets relative to the YILN site (Fig. 8b). The horizontal velocity is separated based on the permanent GPS station ILAN and decreases to the north and south. However, the decreasing rates of −3 and 3 mm yr−1 in these two directions are canceled out (Fig. 8b). Thus, the horizontal motion does not contribute to the LOS mean rate. Across the area north of the ILAN station, only the vertical movement contributes to the LOS displacement mean rate. In Zone (3), the GPS, leveling and SBAS data sets show similarities: the values of each decrease from south to north (Fig. 8a). The horizontal conversion rate change decreases by 3 mm yr−1 from the YILN station to the G250 GPS sites (Fig. 8b). However, no dramatic rate changes were observed in the leveling survey data from 7062 to 7063 (Fig. 8a). Thus, the horizontal subcomponents contributed more to the movement to the north of the YILN GPS station. Zone (4) contains the Lanyang River; thus, no PS pixels exist. Additionally, the GPS and leveling data are sparse. In Zone (5), an area of subsidence extends from the town of Luotung to the eastern coastline (blue zone in Fig. 7a). The maximum subsidence rate reaches −11 mm yr−1 relative to the YILN site along the LOS direction (Fig. 7a). Additionally, the location of this subsidence is consistent with the subsidence observed in the StaMPS data (Kang et al.2015). First, the northern and southern boundaries of the subsidence centre correspond to the Choshui and Sansing faults, respectively (Fig. 8b). The NW–SE extension along these faults is confirmed by the coseismic displacement field of the 2005 Ilan earthquake doublet (Lai et al.2009). On the other hand, the traces of these faults are not obvious in map view (Fig. 7a). The interseismic displacement fields throughout the Ilan Plain have been summarized in Table 1. Most of the regions are controlled by vertical activity (V). Based on the above analysis, we speculate that the subsidence might be controlled by a series of normal faults in the subsurface basement, although their deformation signals are subtle and mixed with human activity. On the other hand, the horizontal component begins to contribute more to the LOS composite velocity in Zone (3) farther to the south. According to the YILN reference site, Zone (3) is also the location where the horizontal velocity vectors begin to change their movement, changing from directed toward the YILN site to a clockwise rotation. The rotational pattern is consistent with that of the backarc extension of the Okinawa Trough. Table 1. Movement tendency comparison between the converted results (LOS) of the horizontal (GPS) and vertical deformation (leveling) and the PSI data. Zone number  LOSSBAS  LOSconv_V  LOSconv_H  LOSconv_SUM  RMS  MD  DS  (1)  3  3  NaN  3  0  V.  NF  (2)  5  6  0  6  −1  V.  NF  (3)  −2.5  0  −3  −3  −0.5  H.  Ext.  (4)  NaN  NaN  NaN  NaN  NaN  NaN  NaN  (5)  4.5  4  −2  2  2.5  V. + H.  NF  (6)  8.5  NaN  NaN  NaN  NaN  V. + H.  NF + Ext.  Zone number  LOSSBAS  LOSconv_V  LOSconv_H  LOSconv_SUM  RMS  MD  DS  (1)  3  3  NaN  3  0  V.  NF  (2)  5  6  0  6  −1  V.  NF  (3)  −2.5  0  −3  −3  −0.5  H.  Ext.  (4)  NaN  NaN  NaN  NaN  NaN  NaN  NaN  (5)  4.5  4  −2  2  2.5  V. + H.  NF  (6)  8.5  NaN  NaN  NaN  NaN  V. + H.  NF + Ext.  Notes: LOSconv_V: rate change in the vertical velocity conversion into the LOS direction; LOSconv_H: rate change in the horizontal velocity conversion into the LOS direction; LOSconv_SUM: sum of LOSconv_V and LOSconv_H; RMS: difference between LOSSBAS and LOSconv_SUM; MD: main deformation; V.: vertical deformation; H.: horizontal deformation; DS: deformation sources; NF: normal faulting; Ext.: lateral extrusion and NaN: no data. View Large 6 INFLUENCE OF THE 2005 ILAN DOUBLET COSEISMIC DEFORMATION The SBAS-PSI mean displacement rate results include one earthquake, known as the 2005 Ilan earthquake doublet that occurred during the studied interseismic period. Lai et al. (2009) simulated the coseismic displacement and suggested that the maximum coseismic slip was less than 50 mm. Although this amount of coseismic displacement could not seemingly contaminate our 6-yr mean interseismic displacement velocity, some scientific questions remain unresolved. For instance, what is the difference between the coseismic and interseismic displacement patterns? Did any strike-slip shear movement associated with the westward propagation of seismic events occur during the coseismic period? Furthermore, no study has applied the DInSAR method to study on the above-mentioned questions. Here, we use a pair of two radar images (captured before and after the 2005 Ilan earthquake doublet) to observe the coseismic displacement associated with the 2005 Ilan earthquake doublet. We also use the strong motion data provided by Lai et al. (2009). 6.1 Coseismic measurements from DInSAR The coseismic deformation due to the 2005 Ilan earthquake was measured using two descending Envisat/ASAR radar images. These two images were obtained three months before and after the 2005 Ilan earthquake doublet (Table 2). One was recorded on 2004 December 11, and the other was recorded on 2005 June 4. Table 2. Envisat/ASAR radar image pair used in this study. Number  Acquisition time  Path  Mode  Orbit  B⊥  Incidence  Ifg1  2004 December 11–2005 June 4  461  Envisat  Desc.  204 m  23.6°  Number  Acquisition time  Path  Mode  Orbit  B⊥  Incidence  Ifg1  2004 December 11–2005 June 4  461  Envisat  Desc.  204 m  23.6°  View Large We used the JPL/Caltech ROI_PAC code (Rosen et al.2004) to process these radar images to produce an interferogram with the shorter perpendicular baselines (∼150 m, Table 2). The coherence was further improved using multilooking methods, typically with 16 ranges. The topographic phase was removed using a simulated phase from the 3 arcsec (∼90 m) SRTM DEM (Farr et al.2007). The phase filtering (Goldstein & Werner 1998), phase unwrapping (e.g. using a minimum cost-flow approach, Chen 2001) and phase-to-displacement conversion and coherence estimation are the same as those of the conventional strip map interferometry (e.g. Fialko 2004; Wright et al.2004). Finally, we obtained a coseismic surface deformation map (Fig. 9). Figure 9. View largeDownload slide Interference fringe map for the coseismic deformation maps of the 2005 Ilan earthquake doublet along the line of sight (LOS). The usage of the ROI_PAC combines the two earthquake events in the Ilan earthquake doublet into one event. The two black beach balls stand for the focal mechanism solutions of the GCMT. Figure 9. View largeDownload slide Interference fringe map for the coseismic deformation maps of the 2005 Ilan earthquake doublet along the line of sight (LOS). The usage of the ROI_PAC combines the two earthquake events in the Ilan earthquake doublet into one event. The two black beach balls stand for the focal mechanism solutions of the GCMT. The coseismic deformation is mainly concentrated on the northern wall of the Sansing fault. Three interference fringes appear to surround its epicentre (Fig. 9). The coherence fringes only appear on the landward side. We therefore speculate that the Sansing fault may extend eastward to the oceanic bottom. Previous research indicated that no surface rupture occurred during the Ilan earthquake doublet (e.g. Liang et al.2005). 6.2 Joint inversion for the 2005 Ilan earthquake To obtain a better first-order description of the 2005 Ilan doublet earthquake source model, we use the SAR interferometry data together with the 13 static strong motion data published by Lai et al. (2009) to jointly inverse the coseismic slip distribution. 6.2.1 Finite-source inversion model We apply a finite fault model developed by Furuya et al. (2010) that can reasonably account for the observed deformation by DInSAR. We assume that the ground displacements are caused by distributed slip on multiple rectangular fault planes in a homogeneous elastic half-space. Using Okada's analytical solutions for the ground displacements due to dislocation sources (Okada 1992), we infer the location, size (length and width), dip angle, strike and slip distribution of the modeled faults. As already noted, the detected ground displacements from the DInSAR data show that the seismogenic structure for the 2005 Ilan earthquake can be simply assigning a single fault plane. Therefore, the fault plane is discretized into many rectangular subfault patches, each with a size of 1 × 1 km along the strike and dip directions, respectively. In order to improve the calculation efficiency, we downsampled the DInSAR data using a 2-D quantization algorithm known as quadtree (Welstead 1999). The number of pixels is reduced from 400 759 to 421. The initial inversions are carried out assuming that there are only left-lateral strike slip with normal slip components of the main shock according to the focal mechanism solutions (e.g. Broadband Array in Taiwan for Seismology, BATS; Global Centroid Moment Tensor Catalogue, GCMT). Therefore, the rake angle is set to vary in the range of [−65°, 65°] and inverse for the fault slip based on assumed fault geometry model. Laplace smoothing constraints of the adjacent fault patches are applied to avoid the abrupt variations of fault slip. The constraint of no significant slip at the fault edges except at the ground surface is applied to further regularize the inversion problem of fault parameters. The best-fitting model predicts that the coseismic rupture occurs along a fault plane with a strike direction of N265.9°W, a dip angle of 36.2°E dipping to the NW and the mean rake angle is −54.1°. The striking direction we searched is consistent with that from the focal mechanism solutions of the GCMT and BATS (see table 1 in Lai et al.2009). And the dip angle is a little smaller than the above-mentioned earthquake catalogues, but similar with the model proposed by Lai et al. (2009). Besides, the rake angle of −54.1° indicates that the seismogenic fault has much more dip-slip component dominated by normal faulting than we previously expected. Because better spatial coverage of the DInSAR data than those form GPS or strong motion, we expect could improve the previous model proposed by Lai et al. (2009). Significant strike-slip and dip-slip distributions are concentrated on the fault patches at depth of 2–5.5 and 3–8 km, respectively. And the main slip zone was located 15 km at length (Fig. 10). The maximum normal faulting and left-lateral shearing are 24 and 31 cm, respectively. Besides, the fault tip does not propagate to the ground surface. It is consistent with the truth that no surface rupture was reported after the earthquake. The seismic moment release of our preferred fault model is 1.80 × 1018 N·m, equivalent to an Mw 6.1 seismic event. This magnitude is larger than those from the of earthquake catalogues such as the GCMT (Mw 5.7) and the BATS (Mw 5.5). The difference could be attributed to the truth that we processed the Ilan doublet as only one seismic event. Figure 10. View largeDownload slide Coseismic slip distribution along the fault plane of the 2005 Ilan earthquake. (a) Strike-slip; (b) dip-slip and (c) slip magnitude. The black arrows indicate the slip vectors. Figure 10. View largeDownload slide Coseismic slip distribution along the fault plane of the 2005 Ilan earthquake. (a) Strike-slip; (b) dip-slip and (c) slip magnitude. The black arrows indicate the slip vectors. 6.2.2 Forward modeling We also used the above optimal fault geometry parameters to inverse the surface displacement. Then, the observed DInSAR results are compared with the simulated ones and the residuals between them are calculated. The horizontal and vertical slip vectors of the strong motion stations are also compared separately. Our simulated DInSAR pattern is roughly consistent with the observed result (see Figs 11a and b). The maximum coseismic displacement appears at the same zone where our interseismic SBAS-PSI measurements have shown (Fig. 8a). Both the left-lateral strike slip and the normal faulting along the Sansing fault determine a faring-away moving tendency along the LOS direction. And the residuals (Fig. 11c) between the observed and the simulated are less than 1.0 cm, especially in the maximum deformation zone (Fig. 11c). Figure 11. View largeDownload slide Forward modeling results of the coseismic displacement field of the 2005 Ilan doublet. (a) Observed; (b) Simulated and (c) Residual of the DInSAR. Comparison with strong motion data. (d) Horizontal and (e) Vertical displacements. The blue and red arrows are used to represent for the observed and simulated slip vectors. The stars stand for the epicentre locations of the 2005 Ilan doublet (Lai et al.2009). The black line represents the Sansing fault proposed seismogenic structure for the 2005 Ilan earthquake. Figure 11. View largeDownload slide Forward modeling results of the coseismic displacement field of the 2005 Ilan doublet. (a) Observed; (b) Simulated and (c) Residual of the DInSAR. Comparison with strong motion data. (d) Horizontal and (e) Vertical displacements. The blue and red arrows are used to represent for the observed and simulated slip vectors. The stars stand for the epicentre locations of the 2005 Ilan doublet (Lai et al.2009). The black line represents the Sansing fault proposed seismogenic structure for the 2005 Ilan earthquake. The horizontal coseismic displacement also shows an NW–SE extension pattern (Fig. 11d). The northern wall of the Sansing fault was moving northwestward; while the southern wall was moving southeastward (Fig. 11d). This pattern is consistent with the model proposed by Lai et al. (2009) model. Except the ILA008 station located at the Ilan Plain was displaying an uplift motion, almost all the other strong motion stations were experiencing a subsidence movement during the 2005 Ilan earthquake (Fig. 11e). Either the horizontal or the vertical coseismic deformation pattern is also consistent with our interseismic SBAS-PSI results (e.g. Fig. 8a) for the Ilan Plain. The modeling results have verified our knowledge achieved through the interseismic monitoring like the GPS, leveling and SBAS-PSI. There exists a set of parallel transtensional faults distributed underlying the Ilan Plain, for example, the Ilan, the Choshui and Sansing faults. They are arranged like a bookshelf-like structure to accommodate the clockwise rotation and the NW–SE opening in the basin. Although the subsidence tendency indicated by the interseismic velocity field by SBAS-PSI could not be only attributed to view tectonic deformation (e.g. anthropogenic activity can also cause land subsidence, see Supporting Information Appendix S1), however, both the co- and interseismic measurements can give us the reassessment of seismic hazards according to the seismogenic structures and strain accumulation across the major faults in the Ilan Plain. 7 DISCUSSION 7.1 Formation of the Ilan Plain and its associated transition mechanism In this paper, we analysed the possible contributionsof horizontal or vertical deformation to the interseismic activity of the Ilan Plain. The entire basin is experiencing fast-lateral extrusion to the southeast and rapid subsidence. Our results from integrated GPS, leveling and SBAS-PSI monitoring indicate that the extension of the basin is mainly controlled by two normal boundary faults on the northwestern and southern margins, the HSRF and the BBRF, respectively. We speculate that these two normal faults differ in their deformation behaviour. The BBRF is likely a passive boundary fault responsible for tearing the Backbone Range away from the Hsüehshan Range and reorienting the front of this range into an E–W direction (Liang et al.2005; Huang et al.2012). Therefore, the Ilan Plain is a typical intermountain basin bounded by the Backbone and Hsüehshan Ranges. Furthermore, the Ilan Plain is situated in the westernmost extension of the Okinawa Trough. The backarc extension of the Okinawa Trough might extend to Taiwan, thereby controlling the opening of the Ilan Plain (Sibuet et al.1987). The Ilan Plain appears to be undergoing decoupled extrusion from the northwest to the southeast at the surface and from shallow depths to the basement because the basin is underlain by basement rocks and an overlying sedimentary cover with two different compositions. The GPS monitoring (e.g. Angelier et al.2009; Hou et al.2009) has shown that the unconsolidated shallow sediments are being extruded laterally and exhibit ductile flow behaviour. At greater depths, however, the lower rigid basal units have been deformed in a typical continuous style and have behaved in an elastic manner. This configuration may explain why there is clockwise rotation around various Euler poles normal to the ground surface (Angelier et al.2009). GPS measurements have shown that an irregular clockwise rotation is observed throughout the Ilan Plain relative to a fixed Eurasia plate (‘irregular’ meaning that rotation occurs around different Euler Poles, Angelier et al.2009; Hou et al.2009). Additionally, the principal shear rate results (Fig. 5b) and the coseismic stress transfer phenomena along the Sansing fault (e.g. Liang et al.2005; Lai et al.2009) indicate that a set of left-lateral structures exist at depth in the basin. Our SBAS-PSI results also indicate that a series of normal faults are present within the basement (see Figs 7a and 8). The basement below the plain is therefore cut by a series of left-lateral normal faults, for example, the Ilan, the Choshui and the Sansing faults (Fig. 12b). These transtensional faults together produce a bookshelf-like structure (Fig. 12a) associated with differential rotation among several second-order blocks and surface subsidence. This movement pattern plays an important role in accommodating the asymmetric transition from arc-continent collision represented by the LS boundary fault and the backarc extension occurring in the Okinawa Trough (Fig. 12b). Figure 12. View largeDownload slide (a) Cartoon map explaining the left-lateral shear senseof the subsurface faults (modified from Liang et al.2005). We propose a bookshelf-like structure. (b) 3-D topographic map of the possible underground structures in the Ilan Plain (depths are not to scale). The opening of the basin is controlled by the bending of the main mountain ranges, the southward retreat of the Ryukyu arc and backarc extension in the Okinawa Trough. A series of transtensional faults is present in the basement. The black arrows indicate clockwise rotation at shallower depths and the extension of the basement. HSRF: Hsüehshan Range Front fault and BBRF: Backbone Range Front fault. Figure 12. View largeDownload slide (a) Cartoon map explaining the left-lateral shear senseof the subsurface faults (modified from Liang et al.2005). We propose a bookshelf-like structure. (b) 3-D topographic map of the possible underground structures in the Ilan Plain (depths are not to scale). The opening of the basin is controlled by the bending of the main mountain ranges, the southward retreat of the Ryukyu arc and backarc extension in the Okinawa Trough. A series of transtensional faults is present in the basement. The black arrows indicate clockwise rotation at shallower depths and the extension of the basement. HSRF: Hsüehshan Range Front fault and BBRF: Backbone Range Front fault. 7.2 Is the Ilan Plain a remnant backarc basin? Several models have been proposed to explain why and how the sudden shift from shortening to extension occurred in Northeast Taiwan. Hsu & Sibuet (1995) and Sibuet & Hsu (1997, 2004) suggested that oblique collision between the Philippine Sea Plate and the Eurasian Plate pushed Taiwan to the northwest (Fig. 13a). The island then acted as wedge, piercing and splitting the pre-existing large-scale Philippine–Ryukyu arc chain (Fig. 13a), which is now expressed in the northeastern and southwestern corners of Taiwan. The Okinawa Trough survives as a backarc basin that is connected with the northeastern corner of the island and includes the Ilan Plain (Figs 13b and c). Teng et al. (2000, 2001) and Shyu et al. (2005) regarded Northeast Taiwan as a post-collisional area. This region appears to no longer experience compression associated with the plate collision and has transitioned into a stage of crustal extension and rifting. To the east of Taiwan's eastern coast, the south-facing Ryukyu arc system is experiencing rapid southerly retreat. This motion has resulted in the opening and translation of the Okinawa Trough, including the westernmost region of extension, that is, the Ilan Plain, to the south. Huang et al. (2012) studied the seismicity and stress patterns in the transitional domain between Ilan and Hualien and found that the fault types change from thrust faults to strike-slip faults to normal faults under the Ilan Plain from south to north. This change is also considered to be caused by the bending of the Central Range belt (Huang et al.2012). Figure 13. View largeDownload slide Cartoon figure illustrating the possible formation history of the Ilan Plain. (a) A pre-existing Ryukyu-Philippine arc chain (RPAC) and the location of Taiwan prior to collision are shown with dotted and solid lines, respectively. TW: Taiwan Island and PSP: Philippine Sea Plate. (b) After arc-continent collision. The Ilan Plain and the Okinawa Trough are located behind the collision zone. LV: Longitudinal Valley fault and LS: Lishan fault. (c) 3-D topographic diagram of Taiwan and its surroundings. The southwestern Okinawa Trough is connected to the Ilan Plain offshore. The LS fault is interpreted to be a suture zone within the former Okinawa-Ilan Plain backarc basin system. The LV fault is the suture zone between the former Ryukyu and Luzon arcs.The blue dotted lines represent the possible evolution history (stages I–III) of the large-scale Ryukyu-Philippine arc chain (modified from Hsu & Sibuet 1995) as it was pierced and ruptured by Taiwan. Figure 13. View largeDownload slide Cartoon figure illustrating the possible formation history of the Ilan Plain. (a) A pre-existing Ryukyu-Philippine arc chain (RPAC) and the location of Taiwan prior to collision are shown with dotted and solid lines, respectively. TW: Taiwan Island and PSP: Philippine Sea Plate. (b) After arc-continent collision. The Ilan Plain and the Okinawa Trough are located behind the collision zone. LV: Longitudinal Valley fault and LS: Lishan fault. (c) 3-D topographic diagram of Taiwan and its surroundings. The southwestern Okinawa Trough is connected to the Ilan Plain offshore. The LS fault is interpreted to be a suture zone within the former Okinawa-Ilan Plain backarc basin system. The LV fault is the suture zone between the former Ryukyu and Luzon arcs.The blue dotted lines represent the possible evolution history (stages I–III) of the large-scale Ryukyu-Philippine arc chain (modified from Hsu & Sibuet 1995) as it was pierced and ruptured by Taiwan. Although Northeast Taiwan experienced uplift during the early stages of the orogenic episode of the island (e.g. Teng & Lin 2004), studies using GPS measurements (e.g. Rau et al.2008; Lin et al.2010), seismotectonic research (e.g. Wu et al.2010; Huang et al.2012) and lithospheric tomography (Kim et al.2005; Wu et al.2010) have all indicated that post-collisional extension rather than compression has dominantly affected the upper crust in Northeast Taiwan. Northeast Taiwan also experiences a considerable number of earthquakes that are associated with normal faults at intermediate depths (Wu et al.2008, 2010). A combination of LOS mean velocity data, leveling measurements and GPS data from the Ilan Plain allows us to integrate interpretations of structures at the surface and at depth and to further clarify the stress regime in Northeast Taiwan. 8 CONCLUSIONS The Ilan Plain experiences a greater degree of extensional stress than compressional stress. This extensional stress arises from both the backarc extension of the Okinawa Trough and the rapid rollback retreat of the Ryukyu Islands. The opening of the basin is controlled by two large-scale normal faults: the HSRF and BBRF normal faults on the northwestern and southern flanks, respectively. Additionally, several E–W-striking transtensional structures, the Choshui and Sansing faults, are distributed in the plain and form a bookshelf-like structure. The LS fault has transformed from a thrust fault into a normal fault. The Ilan Plain may also be a remnant backarc basin, similar to the Okinawa backarc basin. The continuous collision of the Philippine Sea Plate and the Chinese continental shelf plowed Taiwan into the margin, splitting the large-scale Ryukyu–Philippine arc chain and leaving the Ilan Plain and Okinawa Trough as two remnant backarc basins that are connected at the northeastern corner of Taiwan. Acknowledgements The authors are grateful to the editor and reviewers for constructive comments and suggestions. This work was supported by the National Natural Science Foundation (grant nos 41461164002, 41502219); Research Grant from Institute of Crustal Dynamics, China Earthquake Administration (grant nos ZDJ2017-15, ZDJ2017-29). This study also benefits the supports from the Silk Road Project (II): Tibet/Himalaya vs. Caucasus/Iran Orogenic Belts form Excellent Research Projects of National Taiwan University (grant no. 101-105R891204) and Integration Tectonic Research in East Asia from Ministry of Science and Technology in Taiwan under the grant 102-2116-M-002-021. We thank the European Space Agency (ESA) for providing the 32 Envisat/ASAR radar images. We are grateful to Professor Jiankun He and Dr Wenhao Shen, Dr. Xin Zhou for the valuable discussions. REFERENCES Angelier J. et al., 2009. Does extrusion occur at both tips of the Taiwan collision belt? Insights from active deformation studies in the Ilan Plain and Pingtung Plain regions, Tectonophysics , 466, 356– 376. Google Scholar CrossRef Search ADS   Bawden G.W., Thatcher W., Stein R.S., Hudnut K.W., Peltzer G., 2001. Tectonic contraction across Los Angeles after removal of groundwater pumping effects, Nature , 412, 812– 815. Google Scholar CrossRef Search ADS PubMed  Berardino P., Fornaro G., Lanari R., Sansosti E., 2002. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms, IEEE Trans. Geosci. Remote Sens. , 40, 2375– 2383. Google Scholar CrossRef Search ADS   Biggs J., Wright T., Lu Z., Parsons B., 2007. Multi-interferogram method for measuring interseismic deformation: Denali Fault, Alaska, Geophys. J. Int. , 170, 1165– 1179. Google Scholar CrossRef Search ADS   Bürgmann R., Hilley G., Ferretti A., Novali F., 2006. Resolving vertical tectonics in the San Francisco Bay Area from permanent scattererInSAR and GPS analysis, Geology , 34, 221– 224. Google Scholar CrossRef Search ADS   Carminati E., Wortel M.J.R., Meijer P.T., Sabadini R., 1998. The two-stage opening of the western–central Mediterranean basins: a forward modeling test to a new evolutionary model, Earth planet. Sci. Lett. , 160, 667– 679. Google Scholar CrossRef Search ADS   Carminati E., Lustrino M., Doglioni C., 2012. Geodynamic evolution of the central and western Mediterranean: tectonics vs igneous petrology constraints, Tectonophysics , 579, 173– 192. Google Scholar CrossRef Search ADS   Cavalié O., Lasserre C., Doin M.-P., Peltzer G., Sun J., Xu X., Shen Z.-K., 2008. Measurement of interseismic strain across the Haiyuan fault (Gansu, China), by InSAR, Earth planet. Sci. Lett ., 275, 246– 257. Google Scholar CrossRef Search ADS   Chen C.W., 2001. Statistical-Cost Network-Flow approaches to two-dimensional phase unwrapping for radar interferometry, PhD thesis , Stanford University, Stanford, CA, USA. Chen W.-S. et al., 2007. Late Holocene paleoearthquake activity in the middle part of the Longitudinal Valley fault, eastern Taiwan, Earth planet. Sci. Lett ., 264, 420– 437. Google Scholar CrossRef Search ADS   Chen K., Yang M., Huang Y.-T., Ching K.-E., Rau R., 2011. Vertical displacement rate field of Taiwan from geodetic levelling data 2000–2008, Surv. Rev. , 43, 296– 302. Google Scholar CrossRef Search ADS   Chiang S.-C., 1976. Seismic prospecting in the Ilan Plain, Min. Tech. , 14, 215– 221. Ching K.-E., Hsieh M.-L., Johnson K.-M., Chen K.-H., Rau R.-J., Yang M., 2011. Modern vertical deformation rates and mountain building in Taiwan from precise leveling and continuous GPS observations, 2000–2008, J. geophys. Res ., 116, B08406, doi:10.1029/2011JB008242. Clark M.B., Fisher D.M., Lu C.-Y., Chen C.-H., 1993. Kinematic analyses of the Hsüehshan Range, Taiwan: a large-scale pop-up structure, Tectonics , 12, 205– 217. Google Scholar CrossRef Search ADS   Dadson S.J. et al., 2003. Links between erosion, runoff variability and seismicity in the Taiwan orogen, Nature , 426, 648– 651. Google Scholar CrossRef Search ADS PubMed  DeMets C., Gordon R.G., Argus D.F., 2010. Geologically current plate motions, Geophys. J. Int. , 181, 1– 80. Google Scholar CrossRef Search ADS   Elliott J.R., Biggs J., Parsons B., Wright T.J., 2008. InSAR slip rate determination on the Altyn Tagh Fault, northern Tibet, in the presence of topographically correlated atmospheric delays, Geophys. Res. Lett. , 35, L12309, doi:10.1029/2008GL033659. Farr T.G. et al., 2007. The shuttle radar topography mission, Rev. Geophys. , 45, RG2004, doi:10.1029/2005RG000183. Ferretti A., Prati C., Rocca F., 2001. Permanent scatterers in SAR interferometry, IEEE Trans. Geosci. Remote Sens. , 39, 8– 20. Google Scholar CrossRef Search ADS   Fialko Y., Simons M., Agnew D., 2001. The complete (3-D) surface displacement field in the epicentral area of the 1999 MW7.1 Hector Mine Earthquake, California, from space geodetic observations, Geophys. Res. Lett . 28, 3063– 3066. Google Scholar CrossRef Search ADS   Fialko Y., 2004. Probing the mechanical properties of seismically active crust with space geodesy: study of the coseismic deformation due to the 1992 Mw7. 3 Landers (southern California) earthquake, J. geophys. Res ., 109, B03307, doi:10.1029/2003JB002756. Frank F., 1966. Deduction of earth strains from survey data, Bull. seism. Soc. Am. , 56, 35– 42. Furuya M., Kobayashi T., Takada Y., Murakami M., 2010. Fault source modeling of the 2008 wenchuan earthquake based on ALOS/PALSAR data, Bull. seism. Soc. Am. , 100, 2750– 2766. Google Scholar CrossRef Search ADS   Goldstein R.M., Werner C.L., 1998. Radar interferogram filtering for geophysical applications, Geophys. Res. Lett. , 25, 4035– 4038. Google Scholar CrossRef Search ADS   Ho C.-S., 1986. Geological Map of Taiwan. 1: 500,000. Cent. Geol. Surv. , MOEA, Taipei. Ho C.-S., 1988. An Introduction to the Geology of Taiwan: Explanatory Text of the Geology Map of Taiwan , 2nd edn, Ministry of Economic Affairs, Taipei, 192 pp. Hooper A., 2008. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches, Geophys. Res. Lett. , 35, L16302, doi:10.1029/2008GL034654. Hooper A., Bekaert D., Spaans K., Arıkan M., 2012. Recent advances in SAR interferometry time series analysis for measuring crustal deformation, Tectonophysics , 514–517, 1– 13, doi:10.1016/j.tecto.2011.10.013. Google Scholar CrossRef Search ADS   Hou C.-S. et al., 2009. The crustal deformation of the Ilan Plain acted as a westernmost extension of the Okinawa Trough, Tectonophysics , 466, 344– 355. Google Scholar CrossRef Search ADS   Hovius N., Stark C.P., Chu H.-T., Lin J.-C., 2000. Supply and removal of sediment in a landslide-dominated mountain belt: Central Range, Taiwan, J. Geol. , 108, 73– 89. Google Scholar CrossRef Search ADS PubMed  Hsiung K.H., Su C.C., Yu H.S., Chang J.H., 2015. Morphology, seismic characteristics and development of the sediment dispersal system along the Taiwan–Luzon convergent margin, Mar. Geophys. Res ., 36, 293– 308. Google Scholar CrossRef Search ADS   Hsu S.K., Sibuet J., 1995. Is Taiwan the result of arc-continent or arc-arc collision?, Earth planet. Sci. Lett ., 136, 315– 324. Google Scholar CrossRef Search ADS   Hu J.-C., Angelier J., Lee J.-C., Chu H.-T., Byrne D., 1996. Kinematics of convergence, deformation and stress distribution in the Taiwan collision area: 2-D finite-element numerical modelling, Tectonophysics , 255, 243– 268. Google Scholar CrossRef Search ADS   Hu J.-C., Yu S.-B., Chu H.-T., Angelier J., 2002. Transition tectonics of northern Taiwan induced by convergence and trench retreat, Geol. Soc. Am. Spec. Pap. , 358, 149– 162. Huang C.Y., Wu W.Y., Chang C.P., Tsao S., Yuan P.B., Lin C.W., Xia K.Y., 1997. Tectonic evolution of accretionary prism in the arc-continent collision terrane of Taiwan, Tectonophysics , 281, 31– 51. Google Scholar CrossRef Search ADS   Huang H.H., Shyu J.B.H., Wu Y.M., Chang C.H., Chen Y.G., 2012. Seismotectonics of northeastern Taiwan: kinematics of the transition from waning collision to subduction and postcollisional extension, J. geophys. Res ., 117, B01313, doi:10.1029/2011JB008852. Jolivet R., Grandin R., Lasserre C., Doin M.-P., Peltzer G., 2011. Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data, Geophys. Res. Lett. , 38, L17311, doi:10.1029/2011GL048757. Jolivet R. et al., 2012. Shallow creep on the Haiyuan Fault (Gansu, China) revealed by SAR Interferometry, J. geophys. Res ., 117, B06401, doi:10.1029/2011JB008732. Jolivet R., Agram P.S., Lin N.Y., Simons M., Doin M.-P., Peltzer G., Li Z., 2014. Improving InSAR geodesy using global atmospheric models, J. geophys. Res. , 119, 2324– 2341. Google Scholar CrossRef Search ADS   Jolivet R., Lasserre C., Doin M.-P., Peltzer G., Sun J., Rong D., Shen Z., Xu X., 2013. Creep spatial and temporal evolution revealed by InSAR time series analysis, ESA Spec. Publ. , 704, 90. Kang C.-C., Chang C.-P., Siame L., Lee J.-C., 2015. Present-day surface deformation and tectonic insights of the extensional Ilan Plain, NE Taiwan, J. Asian. Earth. Sci. , 105, 408– 417. Google Scholar CrossRef Search ADS   Kao H., Rau R.-J., 1999. Detailed structures of the subducted Philippine Sea plate beneath northeast Taiwan: a new type of double seismic zone, J. geophys. Res. , 104, 1015– 1033. Google Scholar CrossRef Search ADS   Kato T., Beavan J., Matsushima T., Kotake Y., Camacho J.T., Nakao S., 2003. Geodetic evidence of back-arc spreading in the Mariana Trough, Geophys. Res. Lett. , 30(12), 1625, doi:10.1029/2002GL016757. Kim K.H., Chiu J.M., Pujol J., Chen K.C., Huang B.S., Yeh Y.H., Shen P., 2005. Three-dimensional VP and VS structural models associated with the active subduction and collision tectonics in the Taiwan region, Geophys. J. Int. , 162, 204– 220. Google Scholar CrossRef Search ADS   Ku C.Y., Hsu S.K., 2009. Crustal structure and deformation at the northern Manila Trench between Taiwan and Luzon islands, Tectonophysics , 466, 229– 240. Google Scholar CrossRef Search ADS   Kuo-Chen H., Wu F., Chang W.-L, Chang C.-Y., Cheng C.Y., Hirata N., 2015. Is the Lishan fault of Taiwan active?, Tectonophysics , 661, 210– 214. Google Scholar CrossRef Search ADS   Lai K.-Y., Chen Y.-G., Wu Y.-M., Avouac J., Kuo Y.-T., Wang Y., Chang C.-H., Lin K.-C., 2009. The 2005 Ilan earthquake doublet and seismic crisis in northeastern Taiwan: evidence for dyke intrusion associated with on-land propagation of the Okinawa Trough, Geophys. J. Int. , 179, 678– 686. Google Scholar CrossRef Search ADS   Lee J.-C., Angelier J., Chu H.-T., 1997. Polyphase history and kinematics of a complex major fault zone in the northern Taiwan mountain belt: the Lishan Fault, Tectonophysics , 274, 97– 115. Google Scholar CrossRef Search ADS   Lee Y.-H., Byrne T., Wang W.-H., Lo W., Rau R.-J., Lu H.-Y., 2015. Simultaneous mountainbuilding in the Taiwan orogenic belt, Geology , 43, 451– 454. Google Scholar CrossRef Search ADS   Li Z., Fielding E.J., Cross P., 2009. Integration of InSAR time-series analysis and water-vapor correction for mapping postseismic motion after the 2003 Bam (Iran) earthquake, IEEE Trans. Geosci. Remote Sens. , 47, 3220– 3230. Google Scholar CrossRef Search ADS   Liang W.-T., Lee J.-C., Kuo B.-Y., 2005. Left-lateral strike-slip faulting in Ilan: lateral extrusion at the transition between Taiwan mountain range and the Okinawa trough, in Geodynamics and Environment in East Asia International Conference and 5th Taiwan-France Earth Science Symposium , pp. 104– 108, Taitung, ROC. Lin K.-C. et al., 2010. GPS crustal deformation, strain rate, and seismic activity after the 1999 Chi-Chi earthquake in Taiwan, J. geophys. Res ., 115, B07404, doi:10.1029/2009JB006417. Lin S.F., Liu P.M., Lai T.H., 2004. Late Holocene pollen sequence of the Ilan Plain, northeastern Taiwan, and its environmental and climatic implications, Terr. Atmos. Ocean. Sci ., 15, 221– 238. Google Scholar CrossRef Search ADS   Liu C.C., 1995. The Ilan plain and the southwestward extending Okinawa Trough, J. geol. Soc. (China) , 38, 229– 242. Lo C.-L., Chang E.T.-Y., Chao B.F., 2013. Relocating the historical 1951 Hualien earthquake in eastern Taiwan based on tide gauge record, Geophys. J. Int. , 192, 854– 860. Google Scholar CrossRef Search ADS   Nakamura M., 2004. Crustal deformation in the central and southern Ryukyu arc estimated from GPS data, Earth planet. Sci. Lett. , 217, 389– 398. Google Scholar CrossRef Search ADS   Okada Y., 1992. Surface deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 92, 1018– 1040. Prescott W., 1976. An extension of Frank's method for obtaining crustal shear strains from survey data, Bull. seism. Soc. Am. , 66, 1847– 1853. Rau R.-J., Ching K.-E., Hu J.-C., Lee J.-C., 2008. Crustal deformation and block kinematics in transition from collision to subduction: global positioning system measurements in northern Taiwan, 1995–2005, J. geophys. Res ., 113, B09404, doi:10.1029/2007JB005414. Remy D., Bonvalot S., Briole P., Murakami M., 2003. Accurate measurements of tropospheric effects in volcanic areas from SAR interferometry data: application to Sakurajima volcano (Japan), Earth planet. Sci. Lett. , 213, 299– 310. Google Scholar CrossRef Search ADS   Rosen P.A., Hensley S., Peltzer G., Simons M., 2004. Updated repeat orbit interferometry package released, EOS, Trans. Am. geophys. Un. , 85(5), 47, doi:10.1029/2004EO050004. Seno T., 1977. The instantaneous rotation vector of the Philippine Sea plate relative to the Eurasian plate, Tectonophysics , 42, 209– 226. Google Scholar CrossRef Search ADS   Shen Z.K., Jackson D.D., Ge B.X., 1996. Crustal deformation across and beyond the Los Angeles basin from geodetic measurements, J. geophys. Res ., 101, 27 957–27 980. Google Scholar CrossRef Search ADS   Shinjo R., Kato Y., 2000. Geochemical constraints on the origin of bimodal magmatism at the Okinawa Trough, an incipient back-arc basin, Lithos , 54, 117– 137. Google Scholar CrossRef Search ADS   Shyu J.B.H., Sieh K., Chen Y.G., Liu C.S., 2005. Neotectonic architecture of Taiwan and its implications for future large earthquakes, J. geophys. Res. , 110, B08402, doi:10.1029/2004JB003251. Sibuet J., Hsu S.K., 1997. Geodynamics of the Taiwan arc-arc collision, Tectonophysics , 274, 221– 251. Google Scholar CrossRef Search ADS   Sibuet J., Hsu S.K., 2004. How was Taiwan created?, Tectonophysics , 379, 159– 181. Google Scholar CrossRef Search ADS   Sibuet J. et al., 1987. Back arc extension in the Okinawa Trough, J. geophys. Res. , 92, 14 041–14 063. Google Scholar CrossRef Search ADS   Suppe J., 1984. Kinematics of arc-continent collision, flipping of subduction, and back-arc spreading near Taiwan, Geol. Soc. China Mem. , 6, 21– 33. Suppe J., 1987. The active Taiwan mountain belt, in The Anatomy of Mountain Ranges , pp. 277– 293, eds Schaer J.-P., Rodgers J., Princeton Univ. Press. Teng L.S., 1990. Geotectonic evolution of Late Cenozoic arc-continent collision in Taiwan, Tectonophysics , 183, 57– 76. Google Scholar CrossRef Search ADS   Teng L.S., Lee C.T., Peng C.H., Chen W.F., Chu C.J., 2001. Origin and geological evolution of the Taipei basin, northern Taiwan, West. Pac. Earth Sci ., 1, 115– 142. Teng L.S., Lin A.T., 2004. Cenozoic tectonics of the China continental margin: insights from Taiwan, Geol. Soc. Lond., Spec. Pub. , 226, 313– 332. Google Scholar CrossRef Search ADS   Teng L.S., Lee C.T., Tsai Y.B., Hsiao L.Y., 2000. Slab breakoff as a mechanism for flipping of subduction polarity in Taiwan, Geology , 28, 155– 158. Google Scholar CrossRef Search ADS   Tsai M.-C., Yu S.-B., Shin T.-C., Kuo K.-W., Leu P.-L., Chang C.-H., Mei-Yi H., 2015. Velocity field derived from Taiwan continuous GPS array (2007–2013), Terr. Atmos. Ocean. Sci ., 26, 527– 556. Google Scholar CrossRef Search ADS   Ustaszewski K., Wu Y.M., Suppe J., Huang H.H., Chang C.H., Carena S., 2012. Crust–mantle boundaries in the Taiwan–Luzon arc-continent collision system determined from local earthquake tomography and 1D models: implications for the mode of subduction polarity reversal, Tectonophysics , 578, 31– 49. Google Scholar CrossRef Search ADS   Wang K.L., Bilek S.L., 2014. Invited review paper: fault creep caused by subduction of rough seafloor relief, Tectonophysics , 610, 1– 24. Google Scholar CrossRef Search ADS   Wang K.L., Chung S.L., Chen C.H., Shinjo R., Yang T.F., Chen C.H., 1999. Post-collisional magmatism around northern Taiwan and its relation with opening of the Okinawa Trough, Tectonophysics , 308, 363– 376. Google Scholar CrossRef Search ADS   Wang H., Wright T.J., 2012. Satellite geodetic imaging reveals internal deformation of western Tibet, Geophys. Res. Lett. , 39, doi:10.1029/2012GL051222. Wang C.S., Yang M.L., Chou C.P., Chang Y.C., Lee C.S., 2000. Westward Extension of the Okinawa Trough at its western end in the Northern Taiwan area: bathymetric and seismological evidence, Terr. Atmos. Ocean. Sci. , 11, 459– 480. Google Scholar CrossRef Search ADS   Wang H., Wright T.J., Biggs J., 2009. Interseismic slip rate of the northwestern Xianshuihe fault from InSAR data, Geophys. Res. Lett ., 36, L03302, doi:10.1029/2008GL036560. Ward S.N., 1994. A multidisciplinary approach to seismic hazard in Southern California, Bull. seism. Soc. A ., 84, 1293– 1309. Welstead S., 1999. Fractal and Wavelet Image Compression Techniques , SPIE Optical Engineering Press. Wessel P., Smith W.H.F., 1998. New, improved version of generic mapping tools released, EOS, Trans. Am. geophys. Un. , 79(47), 579, doi:10.1029/98EO00426. Wessel P., Smith W.H.F., Scharroo R., Luis J., Wobbe F., 2013. Generic mapping tools: improved version released, EOS, Trans. Am. geophys. Un. , 94, 409– 410. Google Scholar CrossRef Search ADS   Wright T., Parsons B., Fielding E., 2001. Measurement of interseismic strain accumulation across the North Anatolian Fault by satellite radar interferometry, Geophys. Res. Lett ., 28, 2117– 2120. Google Scholar CrossRef Search ADS   Wright T.J., Parsons B.E., Lu Z., 2004. Toward mapping surface deformation in three dimensions using InSAR, Geophys. Res. Lett. , 31, L01607, doi:10.1029/2003GL018827. Wu Y.-M., Zhao L., Chang C.-H., Hsu Y.-J., 2008. Focal-mechanism determination in Taiwan by genetic algorithm, Bull. seism. Soc. Am. , 98, 651– 661. Google Scholar CrossRef Search ADS   Wu Y.-M., Hsu Y.J., Chang C.H., Teng L.S., Nakamura M., 2010. Temporal and spatial variation of stress field in Taiwan from 1991 to 2007: insights from comprehensive first motion focal mechanism catalog, Earth planet. Sci. Lett. , 298, 306– 316. Google Scholar CrossRef Search ADS   Yang M., Chen K.H., Shiao S.W., 2003. A new height reference network in Taiwan, Surv. Rev ., 37, 260– 268. Google Scholar CrossRef Search ADS   Yu S.-B., Kuo L.-C., 2001. Present-day crustal motion along the longitudinal Valley Fault, eastern Taiwan, Tectonophysics , 333, 199– 217. Google Scholar CrossRef Search ADS   Yu S.-B., Chen H.-Y., Kuo L.-C., 1997. Velocity field of GPS stations in the Taiwan area, Tectonophysics , 274, 41– 59. Google Scholar CrossRef Search ADS   Yu S.-B., Kuo L.-C., Punongbayan R.S., Ramos E.G., 1999. GPS observation of crustal deformation in the Taiwan-Luzon region, Geophys. Res. Lett ., 26, 923– 926. Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. APPENDIX S.1 Is the long-term deformation of the Ilan Plain controlled by tectonic or non-tectonic signals? Figure S1. Comparison between long-term (2003–2009) surface deformation from SBAS-PSI and groundwater-level (2005 January–2009 June) changes. Triangles represent a long-term rise in the groundwater level, and inverted triangles represent a long-term decline in the groundwater level. Triangles/inverted triangles outlined in red denote that pumping activity was observed in the aquifers during the monitoring period. Triangles/inverted triangles outlined in black denote that no pumping activity was observed. Figure S2. Long-term groundwater-level changes at the (a) WJ and (b) ZX stations. Each station is equipped with water loggers at four different depths. Red dashed lines are least-squares polynomial fittings to the first-degree polynomial and represent the long-term groundwater-level changes. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Monitoring interseismic activity on the Ilan Plain (NE Taiwan) using Small Baseline PS-InSAR, GPS and leveling measurements: partitioning from arc-continent collision and backarc extension

, Volume 212 (1) – Jan 1, 2018
20 pages

/lp/ou_press/monitoring-interseismic-activity-on-the-ilan-plain-ne-taiwan-using-UztpWa54w5
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx394
Publisher site
See Article on Publisher Site

### Abstract

Abstract The Ilan Plain, located in Northeast Taiwan, represents a transition zone between oblique collision (between the Luzon Arc and the Eurasian Plate) and backarc extension (the Okinawa Trough). The mechanism for this abrupt transition from arc-continent collision to backarc extension remains uncertain. We used Global Positioning System (GPS), leveling and multi-interferogram Small Baseline Persistent Scatterer Interferometry (SBAS-PSI) data to monitor the interseismic activity in the basin. A common reference site was selected for the data sets. The horizontal component of GPS and the vertical measurements of the leveling data were converted to line-of-sight (LOS) data and compared with the SBAS-PSI data. The comparison shows that the entire Ilan Plain is undergoing rapid subsidence at a maximum rate of −11 ± 2 mm yr−1 in the LOS direction. We speculate that vertical deformation and anthropogenic activity may play important roles in this deformation. We also performed a joint inversion modeling that combined both the DInSAR and strong motion data to constrain the source model of the 2005 Ilan earthquake. The best-fitting model predicts that the Sansing fault caused the 2005 Ilan earthquake. The observed transtensional deformation is dominated by the normal faulting with a minor left-lateral strike-slip motion. We compared our SBAS-PSI results with the short-term (2005–2009) groundwater level changes. The results indicate that although pumping-induced surface subsidence cannot be excluded, tectonic deformation, including rapid southward movement of the Ryukyu arc and backarc extension of the Okinawa Trough, characterizes the opening of the Ilan Plain. Furthermore, a series of normal and left-lateral strike-slip transtensional faults, including the Choshui and Sansing faults, form a bookshelf-like structure that accommodates the extension of the plain. Although situated in a region of complex structural interactions, the Ilan Plain is primarily controlled by extension rather than by shortening. As the massive, pre-existing Philippines–Ryukyu island arc was pierced by the Philippine Sea Plate, the Ilan Plain formed as a remnant backarc basin on the northeastern corner of Taiwan. Radar interferometry, Asia, Seismicity and tectonics 1 INTRODUCTION Arc-continent collision and backarc extension are entirely opposite kinematic phenomena in classical orogenic theory. Arc-continent collisions cause mountain building via compressional structures, which typically occur within suture belts that form between converging plates. For example, Taiwan formed from the collision between the Philippine Sea Plate and the Chinese continental shelf (e.g. Huang et al.1997; Lee et al.2015). In contrast, backarc extension represents tectonic collapse that results from the relaxation of tension, the opposite of arc-continent collision. Backarc extension frequently occurs in interior continental areas behind collisional margins. For example, backarc basins are present in the central and western Mediterranean (e.g. Carminati et al.1998, 2012). Northeast Taiwan is one of the few areas on Earth that is currently experiencing two different stress regimes (Hu et al.1996, 2002). On the one hand, the Philippine Sea Plate is obliquely colliding with the southeastern margin of the Chinese continental shelf at a rate of 82 mm yr−1 with respect to the Eurasia Plate (Fig. 1, Yu et al.1997, 1999; Lin et al.2010; Tsai et al.2015). This convergence not only causes crustal thickening along several main mountain ranges in Taiwan (e.g. the Hsüehshan Range, Clark et al.1993; the Central Range, Hovius et al.2000) but also pushes the entire island northwestward (GPS monitoring, e.g. Yu et al.1997; Lin et al.2010) toward the interior of Eurasia (Suppe 1984, 1987; Teng 1990; Sibuet & Hsu 2004). Northeast Taiwan represents an exception to these trends in that this region is experiencing rapid southeastward extrusion (Angelier et al.2009; Hou et al.2009) along with the Ryukyu Islands (Fig. 1; Yu et al.1999; Nakamura 2004). Figure 1. View largeDownload slide Tectonic framework of Taiwan and the surrounding area. The Philippine Sea Plate is moving northwestward (directed N304.1°W) at a rate of 78 mm yr−1 in the MORVEL plate model (DeMets et al.2010). The Ryukyu arc-trench system is retreating to the south with a migration velocity of ∼30–40 mm yr−1. The GPS vectors are referenced to the fixed Chinese continental margin (permanent fixed GPS site SHAO in Shanghai, China; GPS data from Yu et al.1999; Nakamura 2004). LS: Lishan fault and LV: Longitudinal Valley fault. Figure 1. View largeDownload slide Tectonic framework of Taiwan and the surrounding area. The Philippine Sea Plate is moving northwestward (directed N304.1°W) at a rate of 78 mm yr−1 in the MORVEL plate model (DeMets et al.2010). The Ryukyu arc-trench system is retreating to the south with a migration velocity of ∼30–40 mm yr−1. The GPS vectors are referenced to the fixed Chinese continental margin (permanent fixed GPS site SHAO in Shanghai, China; GPS data from Yu et al.1999; Nakamura 2004). LS: Lishan fault and LV: Longitudinal Valley fault. One possible cause of this transition in Northeast Taiwan is a subduction polarity reversal of the Philippine Sea Plate (e.g. Kim et al.2005; Huang et al.2012; Ustaszewski et al.2012). The northern margin of the Philippine Sea Plate is subducting northward beneath the oceanic crust of the Okinawa Trough (Kuo-Chen et al.2015). Backarc basins typically form continent ward of plate collisions, for example, the Okinawa backarc basin (e.g. Sibuet et al.1987; Sibuet & Hsu 1997; Wang et al.1999; Shinjo & Kato 2000). To the south, however, the southwestern margin of the Philippine Sea Plate is obducting onto and overriding the South China Sea basin along the Manila trench (e.g. Seno 1977; Kato et al.2003; Ku & Hsu 2009; Wang & Bilek 2014; Hsiung et al.2015). Consequently, two entirely different tectonic domains have formed within a short distance (<100 km) of this margin in the vicinity of Northeast Taiwan. To the east and south, these domains consist of an extensional regime associated with the westernmost edge of the Okinawa backarc extension and coupled to the rapid southward rollback of the Ryukyu trench. A compressional domain to the west is represented by crustal shortening and thrusting along a major plate boundary and/or reverse faults, represented by the Lishan (LS) fault (Lee et al.1997), which exhibits overthrusting deformation behaviour (Kuo-Chen et al.2015). Compared to other deformation phenomena, interseismic strain accumulation is small in magnitude and is distributed over large spatial distances and temporal spans. The InSAR technique may have the potential to produce regional-scale deformation maps that can identify all active structures and measure the rate of strain accumulation across these active structures. However, short-term monitoring of these tiny but seismically active faults is not easy, especially for active intracontinental structures. Additionally, the expected low-amplitude deformations are comparable in scale to contamination by other errors, such as orbital ramps, temporal variations in the stratified troposphere (e.g. Li et al.2009; Jolivet et al.2011) and topographic effects. Quantitative analysis of interseismic deformation is imperative to understanding the tectonics of an area. The complementary strengths of InSAR, GPS and leveling measurements allow interseismic displacements to be inferred by appropriately combining measurements. Although the displacement fields of the Ilan Plain have been constructed from Permanent Scatterers SAR interferometry (PS-InSAR, e.g. Kang et al.2015) or GPS data (e.g. Angelier et al.2009; Hou et al.2009), a complete and spatially continuous interseismic deformation field associated with the neotectonic activity of the basin is not yet available. InSAR measurements are spatially continuous but are sensitive to only the surface displacement in the line-of-sight (LOS) direction of the radar (e.g. Wright et al.2004). GPS observations provide 3-D velocity vectors but are available only in limited quantities. Leveling measurements face the same problem. Additionally, these three data sets have different reference frames; for instance, leveling (Yang et al.2003; Chen et al.2011) and GPS measurements are relative to a tidal gauge in Keelung along the coast of North Taiwan and the S01R station on Paisha Island (e.g. Angelier et al.2009; Hou et al.2009) and the LOS mean rate is referenced to a ground-to-satellite radar system (Kang et al.2015). In this paper, we monitor interseismic deformation in the Ilan Plain by appropriately integrating InSAR, leveling and GPS measurements in the area. The different reference frames of these techniques are united into a common reference frame. A newly developed method called the Full Rank Matrix Small Baseline Subset InSAR processor (FRAM-SBAS, Biggs et al.2007; Li et al.2009) is applied to produce Small Baseline Persistent Scatterer Interferometry (SBAS-PSI). We compare the horizontal movement based on the GPS data and the vertical movement from the leveling with the movement tendencies of Persistent Scatterers (PS). We thereby characterize which patterns of deformation control the opening of the Ilan Plain, especially focusing on possible subsurface structures and the associated accommodation mechanism between arc-continent collision and backarc extension. 2 SEISMOTECTONIC SETTINGS The Ilan Plain of Northeast Taiwan is situated between two main mountain ranges on Taiwan Island: the Hsüehshan Range to the north and the Backbone Range to the south. Two seismic zones surround the plain. An area known as the Okinawa seismic belt occupies the area between the eastern coast of the Ilan Plain and the Okinawa Trough (seismic zone I; Fig. 2, Kao & Rau 1999). This seismic zone includes a series of normal faults and grabens that generally strike in an east–west direction (Kao & Rau 1999; Wang et al.2000; Huang et al.2012). The second seismic zone, known as the ‘Collision seismic zone’, lies just northeast offshore of Hualien County, Taiwan (Lo et al.2013), and forms an elbow-shaped concentration of earthquakes that occur at depths of approximately 20–30 and 30–50 km (yellow and red points in Fig. 2). The colliding seismic zone represents the location where the reversal of subduction polarity begins. The subduction of the Philippine Sea Plate shifts from northward to northwestward at this hinge point (Huang et al.2012). Figure 2. View largeDownload slide Seismic activity in Northeast Taiwan. The beach balls represent the focal mechanism solutions for earthquake events (red). All the focal mechanism solutions are calculated from the Broadband Array in Taiwan for Seismology (BATS) catalogue. The background seismicity is colour coded by focal depth (the events are from the BATS catalogue). The light grey lines represent the boundary of lithology. The black dotted lines ‘w1’ (Shyu et al.2005) and ‘w2’ (Lai et al.2009; Kang et al.2015) represent two possible northward projections of the Lishan fault. (I) Okinawa Seismic Zone and (II) Collision Seismic Zone. Figure 2. View largeDownload slide Seismic activity in Northeast Taiwan. The beach balls represent the focal mechanism solutions for earthquake events (red). All the focal mechanism solutions are calculated from the Broadband Array in Taiwan for Seismology (BATS) catalogue. The background seismicity is colour coded by focal depth (the events are from the BATS catalogue). The light grey lines represent the boundary of lithology. The black dotted lines ‘w1’ (Shyu et al.2005) and ‘w2’ (Lai et al.2009; Kang et al.2015) represent two possible northward projections of the Lishan fault. (I) Okinawa Seismic Zone and (II) Collision Seismic Zone. In addition to the Longitudinal Valley fault, which acts as a plate boundary between the Philippine Sea Plate and the Eurasian Plate (e.g. Yu & Kuo 2001; Chen et al.2007), another large-scale boundary fault, the LS fault, separates the Hsüehshan Range from the Backbone Range to the west of the Ilan Plain. The LS fault is a steeply dipping reverse fault with dip directions that change from eastward to westward along strike (Lee et al.1997). This structure separates two major geological structural units within the Taiwan orogeny: the Mio-Oligocene Daliao Group, which mainly occurs in the Hsüehshan Range, and the Oligocene Chiayang Group, which mainly occurs in the Backbone Range (Fig. 3, Teng 1990; Lee et al.1997). How far the LS fault extends to the north is a matter of debate. Some researchers have argued that the LS fault extends northwest along the northwestern margin of the Ilan basin (route ‘w1’ in Fig. 2; Shyu et al.2005). The fault becomes a normal fault upon entering the Ilan Plain (Shyu et al.2005). Other models have stated that the LS fault bends to the south and may connect to the east–west-striking Choshui fault at depth when it enters the Ilan Plain (route ‘w2’ in Fig. 2; Kang et al.2015). Figure 3. View largeDownload slide (a) Geological map (Ho 1986, 1988) and cross-sections of the Ilan Plain. The thin black contour lines represent the thickness of the fluvial sediments. (b) Lines 1 and 2 show transects based on exploratory seismic profiles, as reported in Chiang (1976). CS: Choshui fault; HZ: Huzing fault; IL: Ilan fault; JX: Jiaoxi fault; KS: Kengssu fault; NT: Niuto fault; TC: Toucheng fault; SS: Sansing fault and YC: Yaichein fault. Figure 3. View largeDownload slide (a) Geological map (Ho 1986, 1988) and cross-sections of the Ilan Plain. The thin black contour lines represent the thickness of the fluvial sediments. (b) Lines 1 and 2 show transects based on exploratory seismic profiles, as reported in Chiang (1976). CS: Choshui fault; HZ: Huzing fault; IL: Ilan fault; JX: Jiaoxi fault; KS: Kengssu fault; NT: Niuto fault; TC: Toucheng fault; SS: Sansing fault and YC: Yaichein fault. The Ilan Plain is a typical equilateral triangular basin, with its eastern coast facing the West Pacific Ocean. The basin is filled with 500–2000 m of fluvial sediments (Liu 1995; Dadson et al.2003). Consequently, no active structures can be observed on the surface. In terms of sedimentation rate, a previous study (Lin et al.2004) based on one well (the Wuyuan site) south of the Nanyang river suggested that the region featured a very rapid sedimentation rate of 21.7 mm yr−1 between 3400 and 3100 yr B.P., a slow sedimentation rate of ∼1.7 mm yr−1 during the next 2000 yr and an intermediate sedimentation rate of 6.9 mm yr−1 during the last 1000 yr. Two previous studies (Rau et al.2008; Ching et al.2011) proposed high slip rates (>20 mm yr−1) along the two boundary faults (Toucheng and Choshui) based on block models using campaign GPS data. However, possible seismogenic and buried structures were not taken into account in their models. Thus, these results may overestimate the seismic potential of these two boundary faults. The basement of the plain consists of consolidated basal layers that are overlain by shallow, uncemented fluvial sediments (Liu 1995). Previous seismic exploration transects revealed that the basement has an asymmetric washbowl-like shape (Fig. 3a, Chiang 1976). The thickest sediment occurs in an inlet where a river flows into the western Pacific (the 1600 m basement contour in Fig. 3a). However, the current subsidence centre seems to have shifted to the southern bank of the Lanyang River (Kang et al.2015). Chiang (1976) identified a series of blind faults in the basement. From north to south, these faults are the Toucheng (TC), Jiaoxi (JX), Niuto (NT), Huzing (HZ), Yaichein (YC), Ilan (IL), Kengssu (KS), Choshui (CS) and Sansing (SS) faults (fault names are based on Kang et al.2015). Chiang (1976) interpreted these faults as high-angle thrust faults (Fig. 3b) according to the stratal contact relationships but could not evaluate their recent activity. Kang et al. (2015) applied PS-InSAR technology to analyse the relationships with the basement faults and regional geodynamic setting. These authors revealed that these reverse faults could change into normal faults because the plain is experiencing subsidence. Lai et al. (2009) inverted the coseismic rupture process of the 2005 Ilan earthquake doublet with GPS and seismicity data sets. One 15-km-long fault was recognized between the Choshui and Sansing faults and was considered to be the seismogenic structure responsible for the 2005 Ilan earthquake doublet. This active structure exhibits normal faulting with a left-lateral sense of shear. 3 GPS AND LEVELING MEASUREMENTS In terms of GPS and leveling measurements, Taiwan is one of the most monitored places on Earth. The Ilan Plain contains eight permanent GPS stations (provided by Lin et al.2010), 29 temporary GPS monitoring points (provided by the Central Geological Survey) and 99 leveling points from four leveling lines in the leveling measurement project known as the 2001 Taiwan Vertical Datum (Yang et al.2003; Chen et al.2011). The time spans of the GPS and leveling measurements that we used in this study differ: the GPS data were collected from 2004 to 2013, and the leveling data were collected from 2000 to 2008 (Chen et al.2011). However, both types of data were converted to average velocities to determine the relatively stable long-term movement trends, for example, on a 10-yr scale. 3.1 Horizontal and vertical displacement fields We used GPS data collected from 2004 to 2013. Some workers have suggested that the Ilan Plain experiences an escape movement to the southeast (e.g. Angelier et al.2009; Hou et al.2009). The plain generally extrudes towards S139°E (Fig. 4a). A maximum velocity of ∼44 mm yr−1 occurs around the Su’ao coast in the southeastern corner of the plain. The velocity vectors in the northern half of the plain are generally slower than those in the southern half. The northern velocities range from 0.6 to 5.9 mm yr−1 towards 15°–160°, and the southern velocities range from 2.7 to 44 mm yr−1 towards 60°–139° (Hou et al.2009). Figure 4. View largeDownload slide (a) Horizontal velocity field relative to the Paisha station, S01R, which is located on Penghu Island along the Chinese continental margin (2004–2013). The white squares mark permanent GPS stations, and the black triangles markcampaign-mode GPS arrays. The velocity vectors are scaled according to a 95 per cent confidence interval. (b) Vertical velocity field (2000–2008) for the Ilan Plain according to leveling data. The elevation was calculated relative to the Keelung tidal gauge on Taiwan's northern coast (star in inset). The elevation estimates have a mean standard deviation of ± 1.64 mm yr−1 (Chen et al.2011). Figure 4. View largeDownload slide (a) Horizontal velocity field relative to the Paisha station, S01R, which is located on Penghu Island along the Chinese continental margin (2004–2013). The white squares mark permanent GPS stations, and the black triangles markcampaign-mode GPS arrays. The velocity vectors are scaled according to a 95 per cent confidence interval. (b) Vertical velocity field (2000–2008) for the Ilan Plain according to leveling data. The elevation was calculated relative to the Keelung tidal gauge on Taiwan's northern coast (star in inset). The elevation estimates have a mean standard deviation of ± 1.64 mm yr−1 (Chen et al.2011). The leveling results show that the Ilan Plain is presently experiencing rapid subsidence. According to the two E-W-trending profiles that transect the plain and are roughly divided by the 121°36΄E meridian, an obvious velocity shift appears between mountainous areas and the relatively flat plain. The vertical vectors increase from north to south along the N–S-striking survey transects. Several urban areas (e.g. Ilan County and the town of Luotung) show rapid subsidence. Near the Backbone Range, both the southern boundary of the basin and the Su’ao coast exhibit uplift. 3.2 Principal dilation and shearing strain rate In contrast to velocity data (e.g. Fig. 4), the strain rate tensor is independent of the reference frame. However, the maximum shearing and dilation strain rate from GPS data sets could reveal crustal deformation for an active fault during interseismic periods and its possible connection to seismic hazard potential (Ward 1994). We derived strain rates from the velocity estimates according to the procedure described by Shen et al. (1996). This method, which follows previous methods (e.g. Frank 1966; Prescott 1976), interpolates the strain rates from local geodetic measurements through continuous functions within the entire network (for details, see Lin et al.2010). The density of the GPS sites in this study was insufficient; therefore, we created triangles and calculated their centre strain rate. Although this method also applies interpolation method, however, the calculated strain rate is relatively more accurate, because the strain rate in the centre is determined based only on neighbouring real GPS sites. Our results indicate that NW–SE extension is occurring widely in the Ilan Plain (black arrows in Fig. 5a). The largest extensional rate is located on the southern bank of the Langyang River, with a strain rate of 2.66 μstrain yr−1 (μ stands for 10e−6; numbered (1) in Fig. 5a) in the NW–SE direction and a shortening rate of 1.28 μstrain yr−1 (numbered (2) in Fig. 5a) in the NE–SW direction (azimuths 132°–150°; Lin et al.2010). This extensional dilation strain rate pattern is consistent with the description of the coseismic displacement field of the 2005 Ilan earthquake doublet (Lai et al.2009) and the opening direction range of the plain (e.g. Yu et al.1997). Furthermore, this pattern indicates clockwise rotation, as shown by the variation in the direction of the extensional axis from north to south (Fig. 5a), which is also consistent with the opening direction of the Ilan Plain. Figure 5. View largeDownload slide (a) Horizontal principal axes of the strain rates. The amount and direction of the principal strain rates are shown by red (shortening) and black (extension) arrows. The green triangles are the locations of the GPS sites. The Ilan, Choshui and Sansing faults are three major active normal faults inferred from the maximum principal strain rate. Their locations are shifted to the right (dotted lines). From north to south, the maximum extensional axis seems to rotate in a clockwise manner. (b) Left-lateral shearing axes of the strain rates. The amount and direction of the shear strain rates are represented by the gold-coloured bold lines. The above three normal faults exhibit left-lateral senses of shear Figure 5. View largeDownload slide (a) Horizontal principal axes of the strain rates. The amount and direction of the principal strain rates are shown by red (shortening) and black (extension) arrows. The green triangles are the locations of the GPS sites. The Ilan, Choshui and Sansing faults are three major active normal faults inferred from the maximum principal strain rate. Their locations are shifted to the right (dotted lines). From north to south, the maximum extensional axis seems to rotate in a clockwise manner. (b) Left-lateral shearing axes of the strain rates. The amount and direction of the shear strain rates are represented by the gold-coloured bold lines. The above three normal faults exhibit left-lateral senses of shear Fig. 5(b) shows the distribution of the left-lateral shear strain rate in the Ilan Plain. We identified a series of E–W-striking discontinuous shear strain rates in the plain, indicating that some faults that underlie the plain may exhibit a left-lateral sense of shear. Considering the fault identified in seismic reflection profiles (Chiang 1976), we propose that lines (1)–(3) may correspond with the Ilan, Choshui and Sansing faults, respectively. A maximum sinistral strain rate of ∼0.7 μstrain yr−1 appears along the eastern segment of the Sansing fault (number (3) in Fig. 5b). Likewise, the epicentre locations of the 2005 Ilan earthquake doublet (Lai et al.2009) are consistent with the Sansing fault trace. 4 SBAS-PSI MONITORING 4.1 Methodology GPS and leveling data are highly accurate in terms of their respective horizontal and vertical components. However, these data are inherently point based. Current SAR interferometry analysis techniques can provide high-density measurements for monitoring land surface deformation. Among these techniques, such as stacking (e.g. Wright et al.2001; π-RATE, Biggs et al.2007; Elliott et al.2008; Wang et al.2009; Wang & Wright 2012), PS-InSAR technology (Ferretti et al.2001; Hooper et al.2012) and SBASs (Berardino et al.2002), PS-InSAR uses the movement of pixels or points that are coherent throughout the data set to determine the deformation in a particular area. StaMPS is a relatively mature platform that integrates commonly used approaches, such as PS-InSAR and SBAS algorithms (e.g. Hooper et al.2012). Indeed, StaMPS makes no a priori assumptions. However, this method does not have the ability to reduce noise from atmospheric delay, orbital ramps and topography. Furthermore, this method must resolve phase ambiguities in interferogram stacking either by searching a pre-defined solution space or by sparse phase-unwrapping methods. The efficiency and success of phase unwrapping cannot be guaranteed. Finally, the StaMPS software suite works well only in dry, rocky and urban areas. The rainy, wet and cloudy Taiwan Island would definitely introduce more atmosphere noise. In this study, we used a newly developed multitemporal SBAS-PSI method to detect the interseismic deformation in the Ilan Plain. This method is an effective tool for reducing atmospheric artefacts (e.g. Remy et al.2003; Cavalié et al.2008; Elliott et al.2008; Li et al.2009; Jolivet et al.2011, 2012, 2013, 2014) and provides more precise deformation signals (e.g. Biggs et al.2007). We applied the FRAM-SBAS method. FRAM-SBAS is a software suite that combines the open-source code of ROI_PAC (v3.1; Rosen et al.2004), The Generic Mapping Tools (Wessel & Smith 1998; Wessel et al.2013) and an inversion algorithm from Li et al. (2009). We constructed a master–slave network or small baseline image pair under the following conditions to minimize the spatial and temporal decorrelation. First, the perpendicular baselines of master–slave pairs were set to be smaller than 500 m. Second, each acquisition node in the network was given at least four link pairs, meaning that each node has a minimum of four connections with other nodes. Finally, we achieved 170 master–slave pairs (Fig. 6c). Coherent points were extracted based on coherent images; the points that were selected as PS shad coherent values that were greater than 0.38 with a standard deviation of <2.0 mm yr−1. We used the 3-D spatio-temporal approach that was described by Hooper (2008) to unwrap the phase measurements. Figure 6. View largeDownload slide (a) Coverage area of Track 461 for the Ilan Plain (yellow colour in map). (b) A total of 170 pairs were produced under the condition of pre-defined small perpendicular baselines. Figure 6. View largeDownload slide (a) Coverage area of Track 461 for the Ilan Plain (yellow colour in map). (b) A total of 170 pairs were produced under the condition of pre-defined small perpendicular baselines. 4.2 Envisat/ASAR images We used 32 C-band Advanced Synthetic Aperture Radar (ASAR) radar images that were acquired between 2003 November 22, and 2009 June 13, from Environmental Satellite (Envisat) on descending Track 461 (Fig. 6a). All the radar images covered most of the study area (Fig. 6a). The time intervals between the images from 2003 and 2004 and the images from 2008 and 2009 exceeded one year (Fig. 6b). The topographic phase was removed using the Shuttle Radar Topography Mission 3-arcsec Digital Elevation Model (DEM; Farr et al.2007). 5 INTERSEISMIC SURFACE DISPLACEMENT FIELDS IN THE ILAN PLAIN 5.1 Small baseline PSI results We obtained a total of 39 881 high-quality SBAS-PSI pixels from 170 interferogram pairs. All the standard deviations of the PSs were <2.0 mm yr−1. Because the study area includes many densely populated areas, the average data density reached approximately 100 PSs km−2. The LOS mean displacement rate in the Ilan Plain was limited to a range of −11 to 7.5 mm yr−1 (Fig. 7a). A more detailed analysis is provided in Section 5.3. The uncertainties were larger in mountainous regions and in the mountain–plain transition zones (Fig. 7b). We speculate that the large transition from uplift to subsidence and the associated unwrapping error are responsible. Figure 7. View largeDownload slide (a) SBAS-PSI velocity field along the LOS for the Ilan Plain relative to the common reference YILN (indicated by the asterisk). The squares (permanent GPS stations) and diamonds (campaign GPS sites) display the comparison results between the GPS and SBAS-PSI data. The search radius for the PSs was set to 2 km surrounding each GPS station. (b) Estimated standard deviation of the LOS velocity. Figure 7. View largeDownload slide (a) SBAS-PSI velocity field along the LOS for the Ilan Plain relative to the common reference YILN (indicated by the asterisk). The squares (permanent GPS stations) and diamonds (campaign GPS sites) display the comparison results between the GPS and SBAS-PSI data. The search radius for the PSs was set to 2 km surrounding each GPS station. (b) Estimated standard deviation of the LOS velocity. 5.2 Comparison with GPS data GPS, leveling and SBAS-PSI monitoring occurs in different reference frames. The GPS data were relative to the Paisha station (S01R), a permanent GPS site on Penghu Island. The leveling data were referenced to the Keelung tidal gauge on the northern coast of Taiwan. The SBAS-PSI velocity field was referenced to a ground-satellite system. We selected a common reference site to compare these three data sets. Generally, the vertical motion of the GPS measurement is not precise than that of the leveling. Thereafter, when we do the conversion from the eastward, northward, vertical velocity onto the LOS direction, we prefer to select the vertical motion provided by a leveling site that is much closed to a permanent/campaign GPS or under our searched radius. If there will not exist such a leveling site or exceed our searching radius, we have to use the vertical measurement of the GPS site. A permanent GPS station situated to the south of the Ilan County, the YILN was selected (snow sign in Fig. 7a). We achieved a total of 3069 PSs that could fall into a searching circle (2 km as the searching radius). Furthermore, the comparison result is improved by using the new common reference. This site also contains a leveling surveying site, numbered 9039. We unified the references for all three data sets. The eastward, northward and vertical motions were converted into the LOS direction by multiplying them by a full unit vector (Fialko et al.2001). Generally, the vertical motion of the GPS measurement is as not precise as that of the leveling measurements. Hence, during the conversion of the eastward, northward, and vertical velocities into the LOS direction, we used the vertical motion from a leveling site if one was close to the permanent/campaign GPS. If no such leveling site existed or exceeded our search radius, we used the vertical measurement from the GPS site. The descending Envisat/ASAR satellite flies with a mean azimuth angle of 12.41° and with a mean incidence angle of 23°, resulting in an estimated unit vector of (0.3724, −0.0880, 0.9239). The SBAS-PSI density is much higher than that of the GPS stations; hence, the average SBAS-PSI LOS velocity was estimated around each GPS station. We searched all the PSs within a certain distance from the common reference site. Then, we calculated their mean LOS displacement velocities. Here, the search distance is 2 km and the average LOS velocity of the YILN is −0.73 mm yr−1. At the same time, the converted eastward, northward and vertical component of the GPS data into the LOS direction is 0.01 mm yr−1. Thus, all the LOS velocities could be shifted into the common reference frame based on the YILN site by adding 0.74 mm yr−1. We compared all the campaign (displayed by diamonds) and permanent (represented by squares) GPS rate conversion results (triangles) with the PS pixels (coloured circles). This process rested on the assumption that the precision of the vertical GPS measurements was acceptable. The comparison results were satisfactory in the centre of the basin but became worse when approaching the boundary of the Ilan Plain, which may be related to the fact that the campaign and permanent GPS measurements were calculated from two surveying approaches with different measuring cycle periods, precisions, consistencies and stabilities. 5.3 Interseismic activity of the Ilan Plain The 1-D LOS measurement is a composite vector because both horizontal and vertical motions contribute to it. Therefore, the method that combines InSAR data with spatially sparse GPS observations has been widely applied in distinguishing the individual contributions from the horizontal and vertical components, for example, isolating the vertical rates associated with non-tectonic processes (e.g. groundwater and oil pumping, Bawden et al.2001) and/or extracting possible tectonic signals (e.g. significant uplift in the San Francisco Bay Area, Bürgmann et al.2006), by removing the contributions of the GPS-derived horizontal velocity field from the InSAR range-change rates. The Ilan Plain is experiencing not only horizontal extrusion but also rapid subsidence. Vertical deformation is always the result of both non-tectonic (including landslides, and seasonal and long-term groundwater level changes) and tectonic deformations (sediment settling, normal faulting and time-dependent earthquake cycle deformation). Thus, when the horizontal components of the GPS data (Ve and Vn) and the vertical measurements of the leveling data (Vu) are converted into the LOS direction (LOSconv) for comparison with the SBAS-PSI measurements, we are able to constrain the relative contributions of the southeastward extrusion and vertical motions in the opening process of the Ilan Plain. The horizontal components should represent purely tectonic deformation (e.g. the lateral extrusion of the block). Unfortunately, the vertical signals could not be easily isolated because they represent a combination of tectonic movement and the effects of anthropogenic activities, such as water pumping. We use the groundwater data from three major wells (ZX, WJ and LZ) to characterize how the long-term deformation of the Ilan Plain controlled by tectonic or non-tectonic signals (Fig. S1, Supporting Information). Our analysis suggests that surface elevation change is the combined effect of multiple aquifers at the same location; thus, pumping-induced surface subsidence cannot be excluded (Fig. S2, Supporting Information). This result is different from the suggestion proposed by Kang et al. (2015), the detailed discussion of pumping-induced deformation in the Ilan Plain is listed in Supporting Information Appendix S1. We multiplied the horizontal (Ve and Vn) and vertical (Vu) measurements by unit vectors of (0.3724, −0.0880, 0) and (0, 0, 0.9239), respectively. The sum of the vertical and horizontal conversion results should theoretically be equal to the average displacement velocity of the PSs (eq. 1).   $${\rm{LO}}{{\rm{S}}_{{\rm{conv}}}} = \Delta ({\rm{V}}{{\rm{u}}_{{\rm{conv}}}}) + \Delta ({\rm{V}}{{\rm{h}}_{{\rm{conv}}}}) = {\rm{LO}}{{\rm{S}}_{{\rm{SBAS}}}}$$ (1)  $${\rm{RMS}} = {\rm{LO}}{{\rm{S}}_{{\rm{SABS}}}} - {\rm{LO}}{{\rm{S}}_{{\rm{conv}}}}$$ (2) where Δ(Vuconv) and Δ(Vhconv) are the LOS conversion results from the horizontal and vertical components, respectively; LOSconv and LOSSBAS are the sum of the LOS conversion results from the horizontal and vertical components and the SBAS-PSI mean rate change, respectively and RMS is the residual error of LOSSBAS and LOSconv. 5.3.1 Boundary normal faults—Zones (1) and (6) The Ilan triangular basin is bounded by two boundary faults located along the southern slope of the Hsüehshan Range and the southern margin of the plain. These faults are typical normal faults, as inferred from a series of fault scarps observed in the field (e.g. Shyu et al.2005; Kang et al.2015). In our LOS mean displacement rate, these two flanking faults display decreasing LOS rates with respect to the YILN site (Fig. 8a). Warm-toned PS points appear in both the hanging and footwalls of the Hsüehshan Range Front (HSRF) and Backbone Range Front faults (BBRF), which are coloured orange and yellow, respectively (see Fig. 7a). Figure 8. View largeDownload slide (a) SBAS-PSI LOS velocity profile that crosses the Ilan Plain from north to south. The projected SBAS-PSI LOS velocities are represented by blue circles, yellow (campaign-mode GPS array) and blue diamonds (permanent GPS array), and magenta triangles (leveling stations). The blue line represents the topographic relief along a profile based on 40-m DEM data in Taiwan. The GPS and leveling data were collected from within a 5-km-wide swath along the cross-section. The SBAS PS pixels were collected from within a 2-km-wide swath. HSRF: Hsüehshan Range Front fault and BBRF: Backbone Range Front fault. (b) Velocity variation trend. The vertical thin black lines represent a set of possible faults. IL: Ilan fault; CS: Choshui fault and SS: Sansing fault. Figure 8. View largeDownload slide (a) SBAS-PSI LOS velocity profile that crosses the Ilan Plain from north to south. The projected SBAS-PSI LOS velocities are represented by blue circles, yellow (campaign-mode GPS array) and blue diamonds (permanent GPS array), and magenta triangles (leveling stations). The blue line represents the topographic relief along a profile based on 40-m DEM data in Taiwan. The GPS and leveling data were collected from within a 5-km-wide swath along the cross-section. The SBAS PS pixels were collected from within a 2-km-wide swath. HSRF: Hsüehshan Range Front fault and BBRF: Backbone Range Front fault. (b) Velocity variation trend. The vertical thin black lines represent a set of possible faults. IL: Ilan fault; CS: Choshui fault and SS: Sansing fault. Zone (1) exhibits similar movement trends in the leveling (e.g. 9033–9035) and PSI data sets (Fig. 8a). However, the density of GPS data is quite low in this region. We could not determine the deformation from the horizontal components. The vertical velocity decreased by ∼2 mm yr−1 from the Hsüehshan Range to the Ilan Plain (Fig. 8b). We speculate that the HSRF is a normal fault that dips to the south. Relative tothe YILN site, the region to the south of the BBRF appears to have an LOS rate that is ∼8.5 mm yr−1 faster than the area to the north (Fig. 8b). However, this tendency could not be verified by any other geodetic data set because no GPS or leveling survey points exist in the Backbone Range. Although the BBRF is generally considered to be a normal fault, we attribute this shortening tendency to the horizontal components because the GPS velocities indicate that the southern corner of the Ilan Plain is moving laterally to the southeast. Meanwhile, the radar satellites monitor the land deformation with a right-side-looking pattern, facing to the west, which shortens the baseline. 5.3.2 Interseismic activity in the Ilan Plain—Zones (2)–(5) Because the Ilan basin contains many large population centres, the tectonic deformation signals are likely juxtaposed with anthropogenic ground subsidence, for example, the massive pumping of groundwater. Therefore, we could not directly correlate every appearance of rate gradients to an active structure. Below, we describe several zones of movement. In Zone (2), velocity gradients of ∼5 and ∼6 mm yr−1 were measured in the SBAS-PSI and leveling data sets relative to the YILN site (Fig. 8b). The horizontal velocity is separated based on the permanent GPS station ILAN and decreases to the north and south. However, the decreasing rates of −3 and 3 mm yr−1 in these two directions are canceled out (Fig. 8b). Thus, the horizontal motion does not contribute to the LOS mean rate. Across the area north of the ILAN station, only the vertical movement contributes to the LOS displacement mean rate. In Zone (3), the GPS, leveling and SBAS data sets show similarities: the values of each decrease from south to north (Fig. 8a). The horizontal conversion rate change decreases by 3 mm yr−1 from the YILN station to the G250 GPS sites (Fig. 8b). However, no dramatic rate changes were observed in the leveling survey data from 7062 to 7063 (Fig. 8a). Thus, the horizontal subcomponents contributed more to the movement to the north of the YILN GPS station. Zone (4) contains the Lanyang River; thus, no PS pixels exist. Additionally, the GPS and leveling data are sparse. In Zone (5), an area of subsidence extends from the town of Luotung to the eastern coastline (blue zone in Fig. 7a). The maximum subsidence rate reaches −11 mm yr−1 relative to the YILN site along the LOS direction (Fig. 7a). Additionally, the location of this subsidence is consistent with the subsidence observed in the StaMPS data (Kang et al.2015). First, the northern and southern boundaries of the subsidence centre correspond to the Choshui and Sansing faults, respectively (Fig. 8b). The NW–SE extension along these faults is confirmed by the coseismic displacement field of the 2005 Ilan earthquake doublet (Lai et al.2009). On the other hand, the traces of these faults are not obvious in map view (Fig. 7a). The interseismic displacement fields throughout the Ilan Plain have been summarized in Table 1. Most of the regions are controlled by vertical activity (V). Based on the above analysis, we speculate that the subsidence might be controlled by a series of normal faults in the subsurface basement, although their deformation signals are subtle and mixed with human activity. On the other hand, the horizontal component begins to contribute more to the LOS composite velocity in Zone (3) farther to the south. According to the YILN reference site, Zone (3) is also the location where the horizontal velocity vectors begin to change their movement, changing from directed toward the YILN site to a clockwise rotation. The rotational pattern is consistent with that of the backarc extension of the Okinawa Trough. Table 1. Movement tendency comparison between the converted results (LOS) of the horizontal (GPS) and vertical deformation (leveling) and the PSI data. Zone number  LOSSBAS  LOSconv_V  LOSconv_H  LOSconv_SUM  RMS  MD  DS  (1)  3  3  NaN  3  0  V.  NF  (2)  5  6  0  6  −1  V.  NF  (3)  −2.5  0  −3  −3  −0.5  H.  Ext.  (4)  NaN  NaN  NaN  NaN  NaN  NaN  NaN  (5)  4.5  4  −2  2  2.5  V. + H.  NF  (6)  8.5  NaN  NaN  NaN  NaN  V. + H.  NF + Ext.  Zone number  LOSSBAS  LOSconv_V  LOSconv_H  LOSconv_SUM  RMS  MD  DS  (1)  3  3  NaN  3  0  V.  NF  (2)  5  6  0  6  −1  V.  NF  (3)  −2.5  0  −3  −3  −0.5  H.  Ext.  (4)  NaN  NaN  NaN  NaN  NaN  NaN  NaN  (5)  4.5  4  −2  2  2.5  V. + H.  NF  (6)  8.5  NaN  NaN  NaN  NaN  V. + H.  NF + Ext.  Notes: LOSconv_V: rate change in the vertical velocity conversion into the LOS direction; LOSconv_H: rate change in the horizontal velocity conversion into the LOS direction; LOSconv_SUM: sum of LOSconv_V and LOSconv_H; RMS: difference between LOSSBAS and LOSconv_SUM; MD: main deformation; V.: vertical deformation; H.: horizontal deformation; DS: deformation sources; NF: normal faulting; Ext.: lateral extrusion and NaN: no data. View Large 6 INFLUENCE OF THE 2005 ILAN DOUBLET COSEISMIC DEFORMATION The SBAS-PSI mean displacement rate results include one earthquake, known as the 2005 Ilan earthquake doublet that occurred during the studied interseismic period. Lai et al. (2009) simulated the coseismic displacement and suggested that the maximum coseismic slip was less than 50 mm. Although this amount of coseismic displacement could not seemingly contaminate our 6-yr mean interseismic displacement velocity, some scientific questions remain unresolved. For instance, what is the difference between the coseismic and interseismic displacement patterns? Did any strike-slip shear movement associated with the westward propagation of seismic events occur during the coseismic period? Furthermore, no study has applied the DInSAR method to study on the above-mentioned questions. Here, we use a pair of two radar images (captured before and after the 2005 Ilan earthquake doublet) to observe the coseismic displacement associated with the 2005 Ilan earthquake doublet. We also use the strong motion data provided by Lai et al. (2009). 6.1 Coseismic measurements from DInSAR The coseismic deformation due to the 2005 Ilan earthquake was measured using two descending Envisat/ASAR radar images. These two images were obtained three months before and after the 2005 Ilan earthquake doublet (Table 2). One was recorded on 2004 December 11, and the other was recorded on 2005 June 4. Table 2. Envisat/ASAR radar image pair used in this study. Number  Acquisition time  Path  Mode  Orbit  B⊥  Incidence  Ifg1  2004 December 11–2005 June 4  461  Envisat  Desc.  204 m  23.6°  Number  Acquisition time  Path  Mode  Orbit  B⊥  Incidence  Ifg1  2004 December 11–2005 June 4  461  Envisat  Desc.  204 m  23.6°  View Large We used the JPL/Caltech ROI_PAC code (Rosen et al.2004) to process these radar images to produce an interferogram with the shorter perpendicular baselines (∼150 m, Table 2). The coherence was further improved using multilooking methods, typically with 16 ranges. The topographic phase was removed using a simulated phase from the 3 arcsec (∼90 m) SRTM DEM (Farr et al.2007). The phase filtering (Goldstein & Werner 1998), phase unwrapping (e.g. using a minimum cost-flow approach, Chen 2001) and phase-to-displacement conversion and coherence estimation are the same as those of the conventional strip map interferometry (e.g. Fialko 2004; Wright et al.2004). Finally, we obtained a coseismic surface deformation map (Fig. 9). Figure 9. View largeDownload slide Interference fringe map for the coseismic deformation maps of the 2005 Ilan earthquake doublet along the line of sight (LOS). The usage of the ROI_PAC combines the two earthquake events in the Ilan earthquake doublet into one event. The two black beach balls stand for the focal mechanism solutions of the GCMT. Figure 9. View largeDownload slide Interference fringe map for the coseismic deformation maps of the 2005 Ilan earthquake doublet along the line of sight (LOS). The usage of the ROI_PAC combines the two earthquake events in the Ilan earthquake doublet into one event. The two black beach balls stand for the focal mechanism solutions of the GCMT. The coseismic deformation is mainly concentrated on the northern wall of the Sansing fault. Three interference fringes appear to surround its epicentre (Fig. 9). The coherence fringes only appear on the landward side. We therefore speculate that the Sansing fault may extend eastward to the oceanic bottom. Previous research indicated that no surface rupture occurred during the Ilan earthquake doublet (e.g. Liang et al.2005). 6.2 Joint inversion for the 2005 Ilan earthquake To obtain a better first-order description of the 2005 Ilan doublet earthquake source model, we use the SAR interferometry data together with the 13 static strong motion data published by Lai et al. (2009) to jointly inverse the coseismic slip distribution. 6.2.1 Finite-source inversion model We apply a finite fault model developed by Furuya et al. (2010) that can reasonably account for the observed deformation by DInSAR. We assume that the ground displacements are caused by distributed slip on multiple rectangular fault planes in a homogeneous elastic half-space. Using Okada's analytical solutions for the ground displacements due to dislocation sources (Okada 1992), we infer the location, size (length and width), dip angle, strike and slip distribution of the modeled faults. As already noted, the detected ground displacements from the DInSAR data show that the seismogenic structure for the 2005 Ilan earthquake can be simply assigning a single fault plane. Therefore, the fault plane is discretized into many rectangular subfault patches, each with a size of 1 × 1 km along the strike and dip directions, respectively. In order to improve the calculation efficiency, we downsampled the DInSAR data using a 2-D quantization algorithm known as quadtree (Welstead 1999). The number of pixels is reduced from 400 759 to 421. The initial inversions are carried out assuming that there are only left-lateral strike slip with normal slip components of the main shock according to the focal mechanism solutions (e.g. Broadband Array in Taiwan for Seismology, BATS; Global Centroid Moment Tensor Catalogue, GCMT). Therefore, the rake angle is set to vary in the range of [−65°, 65°] and inverse for the fault slip based on assumed fault geometry model. Laplace smoothing constraints of the adjacent fault patches are applied to avoid the abrupt variations of fault slip. The constraint of no significant slip at the fault edges except at the ground surface is applied to further regularize the inversion problem of fault parameters. The best-fitting model predicts that the coseismic rupture occurs along a fault plane with a strike direction of N265.9°W, a dip angle of 36.2°E dipping to the NW and the mean rake angle is −54.1°. The striking direction we searched is consistent with that from the focal mechanism solutions of the GCMT and BATS (see table 1 in Lai et al.2009). And the dip angle is a little smaller than the above-mentioned earthquake catalogues, but similar with the model proposed by Lai et al. (2009). Besides, the rake angle of −54.1° indicates that the seismogenic fault has much more dip-slip component dominated by normal faulting than we previously expected. Because better spatial coverage of the DInSAR data than those form GPS or strong motion, we expect could improve the previous model proposed by Lai et al. (2009). Significant strike-slip and dip-slip distributions are concentrated on the fault patches at depth of 2–5.5 and 3–8 km, respectively. And the main slip zone was located 15 km at length (Fig. 10). The maximum normal faulting and left-lateral shearing are 24 and 31 cm, respectively. Besides, the fault tip does not propagate to the ground surface. It is consistent with the truth that no surface rupture was reported after the earthquake. The seismic moment release of our preferred fault model is 1.80 × 1018 N·m, equivalent to an Mw 6.1 seismic event. This magnitude is larger than those from the of earthquake catalogues such as the GCMT (Mw 5.7) and the BATS (Mw 5.5). The difference could be attributed to the truth that we processed the Ilan doublet as only one seismic event. Figure 10. View largeDownload slide Coseismic slip distribution along the fault plane of the 2005 Ilan earthquake. (a) Strike-slip; (b) dip-slip and (c) slip magnitude. The black arrows indicate the slip vectors. Figure 10. View largeDownload slide Coseismic slip distribution along the fault plane of the 2005 Ilan earthquake. (a) Strike-slip; (b) dip-slip and (c) slip magnitude. The black arrows indicate the slip vectors. 6.2.2 Forward modeling We also used the above optimal fault geometry parameters to inverse the surface displacement. Then, the observed DInSAR results are compared with the simulated ones and the residuals between them are calculated. The horizontal and vertical slip vectors of the strong motion stations are also compared separately. Our simulated DInSAR pattern is roughly consistent with the observed result (see Figs 11a and b). The maximum coseismic displacement appears at the same zone where our interseismic SBAS-PSI measurements have shown (Fig. 8a). Both the left-lateral strike slip and the normal faulting along the Sansing fault determine a faring-away moving tendency along the LOS direction. And the residuals (Fig. 11c) between the observed and the simulated are less than 1.0 cm, especially in the maximum deformation zone (Fig. 11c). Figure 11. View largeDownload slide Forward modeling results of the coseismic displacement field of the 2005 Ilan doublet. (a) Observed; (b) Simulated and (c) Residual of the DInSAR. Comparison with strong motion data. (d) Horizontal and (e) Vertical displacements. The blue and red arrows are used to represent for the observed and simulated slip vectors. The stars stand for the epicentre locations of the 2005 Ilan doublet (Lai et al.2009). The black line represents the Sansing fault proposed seismogenic structure for the 2005 Ilan earthquake. Figure 11. View largeDownload slide Forward modeling results of the coseismic displacement field of the 2005 Ilan doublet. (a) Observed; (b) Simulated and (c) Residual of the DInSAR. Comparison with strong motion data. (d) Horizontal and (e) Vertical displacements. The blue and red arrows are used to represent for the observed and simulated slip vectors. The stars stand for the epicentre locations of the 2005 Ilan doublet (Lai et al.2009). The black line represents the Sansing fault proposed seismogenic structure for the 2005 Ilan earthquake. The horizontal coseismic displacement also shows an NW–SE extension pattern (Fig. 11d). The northern wall of the Sansing fault was moving northwestward; while the southern wall was moving southeastward (Fig. 11d). This pattern is consistent with the model proposed by Lai et al. (2009) model. Except the ILA008 station located at the Ilan Plain was displaying an uplift motion, almost all the other strong motion stations were experiencing a subsidence movement during the 2005 Ilan earthquake (Fig. 11e). Either the horizontal or the vertical coseismic deformation pattern is also consistent with our interseismic SBAS-PSI results (e.g. Fig. 8a) for the Ilan Plain. The modeling results have verified our knowledge achieved through the interseismic monitoring like the GPS, leveling and SBAS-PSI. There exists a set of parallel transtensional faults distributed underlying the Ilan Plain, for example, the Ilan, the Choshui and Sansing faults. They are arranged like a bookshelf-like structure to accommodate the clockwise rotation and the NW–SE opening in the basin. Although the subsidence tendency indicated by the interseismic velocity field by SBAS-PSI could not be only attributed to view tectonic deformation (e.g. anthropogenic activity can also cause land subsidence, see Supporting Information Appendix S1), however, both the co- and interseismic measurements can give us the reassessment of seismic hazards according to the seismogenic structures and strain accumulation across the major faults in the Ilan Plain. 7 DISCUSSION 7.1 Formation of the Ilan Plain and its associated transition mechanism In this paper, we analysed the possible contributionsof horizontal or vertical deformation to the interseismic activity of the Ilan Plain. The entire basin is experiencing fast-lateral extrusion to the southeast and rapid subsidence. Our results from integrated GPS, leveling and SBAS-PSI monitoring indicate that the extension of the basin is mainly controlled by two normal boundary faults on the northwestern and southern margins, the HSRF and the BBRF, respectively. We speculate that these two normal faults differ in their deformation behaviour. The BBRF is likely a passive boundary fault responsible for tearing the Backbone Range away from the Hsüehshan Range and reorienting the front of this range into an E–W direction (Liang et al.2005; Huang et al.2012). Therefore, the Ilan Plain is a typical intermountain basin bounded by the Backbone and Hsüehshan Ranges. Furthermore, the Ilan Plain is situated in the westernmost extension of the Okinawa Trough. The backarc extension of the Okinawa Trough might extend to Taiwan, thereby controlling the opening of the Ilan Plain (Sibuet et al.1987). The Ilan Plain appears to be undergoing decoupled extrusion from the northwest to the southeast at the surface and from shallow depths to the basement because the basin is underlain by basement rocks and an overlying sedimentary cover with two different compositions. The GPS monitoring (e.g. Angelier et al.2009; Hou et al.2009) has shown that the unconsolidated shallow sediments are being extruded laterally and exhibit ductile flow behaviour. At greater depths, however, the lower rigid basal units have been deformed in a typical continuous style and have behaved in an elastic manner. This configuration may explain why there is clockwise rotation around various Euler poles normal to the ground surface (Angelier et al.2009). GPS measurements have shown that an irregular clockwise rotation is observed throughout the Ilan Plain relative to a fixed Eurasia plate (‘irregular’ meaning that rotation occurs around different Euler Poles, Angelier et al.2009; Hou et al.2009). Additionally, the principal shear rate results (Fig. 5b) and the coseismic stress transfer phenomena along the Sansing fault (e.g. Liang et al.2005; Lai et al.2009) indicate that a set of left-lateral structures exist at depth in the basin. Our SBAS-PSI results also indicate that a series of normal faults are present within the basement (see Figs 7a and 8). The basement below the plain is therefore cut by a series of left-lateral normal faults, for example, the Ilan, the Choshui and the Sansing faults (Fig. 12b). These transtensional faults together produce a bookshelf-like structure (Fig. 12a) associated with differential rotation among several second-order blocks and surface subsidence. This movement pattern plays an important role in accommodating the asymmetric transition from arc-continent collision represented by the LS boundary fault and the backarc extension occurring in the Okinawa Trough (Fig. 12b). Figure 12. View largeDownload slide (a) Cartoon map explaining the left-lateral shear senseof the subsurface faults (modified from Liang et al.2005). We propose a bookshelf-like structure. (b) 3-D topographic map of the possible underground structures in the Ilan Plain (depths are not to scale). The opening of the basin is controlled by the bending of the main mountain ranges, the southward retreat of the Ryukyu arc and backarc extension in the Okinawa Trough. A series of transtensional faults is present in the basement. The black arrows indicate clockwise rotation at shallower depths and the extension of the basement. HSRF: Hsüehshan Range Front fault and BBRF: Backbone Range Front fault. Figure 12. View largeDownload slide (a) Cartoon map explaining the left-lateral shear senseof the subsurface faults (modified from Liang et al.2005). We propose a bookshelf-like structure. (b) 3-D topographic map of the possible underground structures in the Ilan Plain (depths are not to scale). The opening of the basin is controlled by the bending of the main mountain ranges, the southward retreat of the Ryukyu arc and backarc extension in the Okinawa Trough. A series of transtensional faults is present in the basement. The black arrows indicate clockwise rotation at shallower depths and the extension of the basement. HSRF: Hsüehshan Range Front fault and BBRF: Backbone Range Front fault. 7.2 Is the Ilan Plain a remnant backarc basin? Several models have been proposed to explain why and how the sudden shift from shortening to extension occurred in Northeast Taiwan. Hsu & Sibuet (1995) and Sibuet & Hsu (1997, 2004) suggested that oblique collision between the Philippine Sea Plate and the Eurasian Plate pushed Taiwan to the northwest (Fig. 13a). The island then acted as wedge, piercing and splitting the pre-existing large-scale Philippine–Ryukyu arc chain (Fig. 13a), which is now expressed in the northeastern and southwestern corners of Taiwan. The Okinawa Trough survives as a backarc basin that is connected with the northeastern corner of the island and includes the Ilan Plain (Figs 13b and c). Teng et al. (2000, 2001) and Shyu et al. (2005) regarded Northeast Taiwan as a post-collisional area. This region appears to no longer experience compression associated with the plate collision and has transitioned into a stage of crustal extension and rifting. To the east of Taiwan's eastern coast, the south-facing Ryukyu arc system is experiencing rapid southerly retreat. This motion has resulted in the opening and translation of the Okinawa Trough, including the westernmost region of extension, that is, the Ilan Plain, to the south. Huang et al. (2012) studied the seismicity and stress patterns in the transitional domain between Ilan and Hualien and found that the fault types change from thrust faults to strike-slip faults to normal faults under the Ilan Plain from south to north. This change is also considered to be caused by the bending of the Central Range belt (Huang et al.2012). Figure 13. View largeDownload slide Cartoon figure illustrating the possible formation history of the Ilan Plain. (a) A pre-existing Ryukyu-Philippine arc chain (RPAC) and the location of Taiwan prior to collision are shown with dotted and solid lines, respectively. TW: Taiwan Island and PSP: Philippine Sea Plate. (b) After arc-continent collision. The Ilan Plain and the Okinawa Trough are located behind the collision zone. LV: Longitudinal Valley fault and LS: Lishan fault. (c) 3-D topographic diagram of Taiwan and its surroundings. The southwestern Okinawa Trough is connected to the Ilan Plain offshore. The LS fault is interpreted to be a suture zone within the former Okinawa-Ilan Plain backarc basin system. The LV fault is the suture zone between the former Ryukyu and Luzon arcs.The blue dotted lines represent the possible evolution history (stages I–III) of the large-scale Ryukyu-Philippine arc chain (modified from Hsu & Sibuet 1995) as it was pierced and ruptured by Taiwan. Figure 13. View largeDownload slide Cartoon figure illustrating the possible formation history of the Ilan Plain. (a) A pre-existing Ryukyu-Philippine arc chain (RPAC) and the location of Taiwan prior to collision are shown with dotted and solid lines, respectively. TW: Taiwan Island and PSP: Philippine Sea Plate. (b) After arc-continent collision. The Ilan Plain and the Okinawa Trough are located behind the collision zone. LV: Longitudinal Valley fault and LS: Lishan fault. (c) 3-D topographic diagram of Taiwan and its surroundings. The southwestern Okinawa Trough is connected to the Ilan Plain offshore. The LS fault is interpreted to be a suture zone within the former Okinawa-Ilan Plain backarc basin system. The LV fault is the suture zone between the former Ryukyu and Luzon arcs.The blue dotted lines represent the possible evolution history (stages I–III) of the large-scale Ryukyu-Philippine arc chain (modified from Hsu & Sibuet 1995) as it was pierced and ruptured by Taiwan. Although Northeast Taiwan experienced uplift during the early stages of the orogenic episode of the island (e.g. Teng & Lin 2004), studies using GPS measurements (e.g. Rau et al.2008; Lin et al.2010), seismotectonic research (e.g. Wu et al.2010; Huang et al.2012) and lithospheric tomography (Kim et al.2005; Wu et al.2010) have all indicated that post-collisional extension rather than compression has dominantly affected the upper crust in Northeast Taiwan. Northeast Taiwan also experiences a considerable number of earthquakes that are associated with normal faults at intermediate depths (Wu et al.2008, 2010). A combination of LOS mean velocity data, leveling measurements and GPS data from the Ilan Plain allows us to integrate interpretations of structures at the surface and at depth and to further clarify the stress regime in Northeast Taiwan. 8 CONCLUSIONS The Ilan Plain experiences a greater degree of extensional stress than compressional stress. This extensional stress arises from both the backarc extension of the Okinawa Trough and the rapid rollback retreat of the Ryukyu Islands. The opening of the basin is controlled by two large-scale normal faults: the HSRF and BBRF normal faults on the northwestern and southern flanks, respectively. Additionally, several E–W-striking transtensional structures, the Choshui and Sansing faults, are distributed in the plain and form a bookshelf-like structure. The LS fault has transformed from a thrust fault into a normal fault. The Ilan Plain may also be a remnant backarc basin, similar to the Okinawa backarc basin. The continuous collision of the Philippine Sea Plate and the Chinese continental shelf plowed Taiwan into the margin, splitting the large-scale Ryukyu–Philippine arc chain and leaving the Ilan Plain and Okinawa Trough as two remnant backarc basins that are connected at the northeastern corner of Taiwan. Acknowledgements The authors are grateful to the editor and reviewers for constructive comments and suggestions. This work was supported by the National Natural Science Foundation (grant nos 41461164002, 41502219); Research Grant from Institute of Crustal Dynamics, China Earthquake Administration (grant nos ZDJ2017-15, ZDJ2017-29). This study also benefits the supports from the Silk Road Project (II): Tibet/Himalaya vs. Caucasus/Iran Orogenic Belts form Excellent Research Projects of National Taiwan University (grant no. 101-105R891204) and Integration Tectonic Research in East Asia from Ministry of Science and Technology in Taiwan under the grant 102-2116-M-002-021. We thank the European Space Agency (ESA) for providing the 32 Envisat/ASAR radar images. We are grateful to Professor Jiankun He and Dr Wenhao Shen, Dr. Xin Zhou for the valuable discussions. REFERENCES Angelier J. et al., 2009. Does extrusion occur at both tips of the Taiwan collision belt? Insights from active deformation studies in the Ilan Plain and Pingtung Plain regions, Tectonophysics , 466, 356– 376. Google Scholar CrossRef Search ADS   Bawden G.W., Thatcher W., Stein R.S., Hudnut K.W., Peltzer G., 2001. Tectonic contraction across Los Angeles after removal of groundwater pumping effects, Nature , 412, 812– 815. Google Scholar CrossRef Search ADS PubMed  Berardino P., Fornaro G., Lanari R., Sansosti E., 2002. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms, IEEE Trans. Geosci. Remote Sens. , 40, 2375– 2383. Google Scholar CrossRef Search ADS   Biggs J., Wright T., Lu Z., Parsons B., 2007. Multi-interferogram method for measuring interseismic deformation: Denali Fault, Alaska, Geophys. J. Int. , 170, 1165– 1179. Google Scholar CrossRef Search ADS   Bürgmann R., Hilley G., Ferretti A., Novali F., 2006. Resolving vertical tectonics in the San Francisco Bay Area from permanent scattererInSAR and GPS analysis, Geology , 34, 221– 224. Google Scholar CrossRef Search ADS   Carminati E., Wortel M.J.R., Meijer P.T., Sabadini R., 1998. The two-stage opening of the western–central Mediterranean basins: a forward modeling test to a new evolutionary model, Earth planet. Sci. Lett. , 160, 667– 679. Google Scholar CrossRef Search ADS   Carminati E., Lustrino M., Doglioni C., 2012. Geodynamic evolution of the central and western Mediterranean: tectonics vs igneous petrology constraints, Tectonophysics , 579, 173– 192. Google Scholar CrossRef Search ADS   Cavalié O., Lasserre C., Doin M.-P., Peltzer G., Sun J., Xu X., Shen Z.-K., 2008. Measurement of interseismic strain across the Haiyuan fault (Gansu, China), by InSAR, Earth planet. Sci. Lett ., 275, 246– 257. Google Scholar CrossRef Search ADS   Chen C.W., 2001. Statistical-Cost Network-Flow approaches to two-dimensional phase unwrapping for radar interferometry, PhD thesis , Stanford University, Stanford, CA, USA. Chen W.-S. et al., 2007. Late Holocene paleoearthquake activity in the middle part of the Longitudinal Valley fault, eastern Taiwan, Earth planet. Sci. Lett ., 264, 420– 437. Google Scholar CrossRef Search ADS   Chen K., Yang M., Huang Y.-T., Ching K.-E., Rau R., 2011. Vertical displacement rate field of Taiwan from geodetic levelling data 2000–2008, Surv. Rev. , 43, 296– 302. Google Scholar CrossRef Search ADS   Chiang S.-C., 1976. Seismic prospecting in the Ilan Plain, Min. Tech. , 14, 215– 221. Ching K.-E., Hsieh M.-L., Johnson K.-M., Chen K.-H., Rau R.-J., Yang M., 2011. Modern vertical deformation rates and mountain building in Taiwan from precise leveling and continuous GPS observations, 2000–2008, J. geophys. Res ., 116, B08406, doi:10.1029/2011JB008242. Clark M.B., Fisher D.M., Lu C.-Y., Chen C.-H., 1993. Kinematic analyses of the Hsüehshan Range, Taiwan: a large-scale pop-up structure, Tectonics , 12, 205– 217. Google Scholar CrossRef Search ADS   Dadson S.J. et al., 2003. Links between erosion, runoff variability and seismicity in the Taiwan orogen, Nature , 426, 648– 651. Google Scholar CrossRef Search ADS PubMed  DeMets C., Gordon R.G., Argus D.F., 2010. Geologically current plate motions, Geophys. J. Int. , 181, 1– 80. Google Scholar CrossRef Search ADS   Elliott J.R., Biggs J., Parsons B., Wright T.J., 2008. InSAR slip rate determination on the Altyn Tagh Fault, northern Tibet, in the presence of topographically correlated atmospheric delays, Geophys. Res. Lett. , 35, L12309, doi:10.1029/2008GL033659. Farr T.G. et al., 2007. The shuttle radar topography mission, Rev. Geophys. , 45, RG2004, doi:10.1029/2005RG000183. Ferretti A., Prati C., Rocca F., 2001. Permanent scatterers in SAR interferometry, IEEE Trans. Geosci. Remote Sens. , 39, 8– 20. Google Scholar CrossRef Search ADS   Fialko Y., Simons M., Agnew D., 2001. The complete (3-D) surface displacement field in the epicentral area of the 1999 MW7.1 Hector Mine Earthquake, California, from space geodetic observations, Geophys. Res. Lett . 28, 3063– 3066. Google Scholar CrossRef Search ADS   Fialko Y., 2004. Probing the mechanical properties of seismically active crust with space geodesy: study of the coseismic deformation due to the 1992 Mw7. 3 Landers (southern California) earthquake, J. geophys. Res ., 109, B03307, doi:10.1029/2003JB002756. Frank F., 1966. Deduction of earth strains from survey data, Bull. seism. Soc. Am. , 56, 35– 42. Furuya M., Kobayashi T., Takada Y., Murakami M., 2010. Fault source modeling of the 2008 wenchuan earthquake based on ALOS/PALSAR data, Bull. seism. Soc. Am. , 100, 2750– 2766. Google Scholar CrossRef Search ADS   Goldstein R.M., Werner C.L., 1998. Radar interferogram filtering for geophysical applications, Geophys. Res. Lett. , 25, 4035– 4038. Google Scholar CrossRef Search ADS   Ho C.-S., 1986. Geological Map of Taiwan. 1: 500,000. Cent. Geol. Surv. , MOEA, Taipei. Ho C.-S., 1988. An Introduction to the Geology of Taiwan: Explanatory Text of the Geology Map of Taiwan , 2nd edn, Ministry of Economic Affairs, Taipei, 192 pp. Hooper A., 2008. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches, Geophys. Res. Lett. , 35, L16302, doi:10.1029/2008GL034654. Hooper A., Bekaert D., Spaans K., Arıkan M., 2012. Recent advances in SAR interferometry time series analysis for measuring crustal deformation, Tectonophysics , 514–517, 1– 13, doi:10.1016/j.tecto.2011.10.013. Google Scholar CrossRef Search ADS   Hou C.-S. et al., 2009. The crustal deformation of the Ilan Plain acted as a westernmost extension of the Okinawa Trough, Tectonophysics , 466, 344– 355. Google Scholar CrossRef Search ADS   Hovius N., Stark C.P., Chu H.-T., Lin J.-C., 2000. Supply and removal of sediment in a landslide-dominated mountain belt: Central Range, Taiwan, J. Geol. , 108, 73– 89. Google Scholar CrossRef Search ADS PubMed  Hsiung K.H., Su C.C., Yu H.S., Chang J.H., 2015. Morphology, seismic characteristics and development of the sediment dispersal system along the Taiwan–Luzon convergent margin, Mar. Geophys. Res ., 36, 293– 308. Google Scholar CrossRef Search ADS   Hsu S.K., Sibuet J., 1995. Is Taiwan the result of arc-continent or arc-arc collision?, Earth planet. Sci. Lett ., 136, 315– 324. Google Scholar CrossRef Search ADS   Hu J.-C., Angelier J., Lee J.-C., Chu H.-T., Byrne D., 1996. Kinematics of convergence, deformation and stress distribution in the Taiwan collision area: 2-D finite-element numerical modelling, Tectonophysics , 255, 243– 268. Google Scholar CrossRef Search ADS   Hu J.-C., Yu S.-B., Chu H.-T., Angelier J., 2002. Transition tectonics of northern Taiwan induced by convergence and trench retreat, Geol. Soc. Am. Spec. Pap. , 358, 149– 162. Huang C.Y., Wu W.Y., Chang C.P., Tsao S., Yuan P.B., Lin C.W., Xia K.Y., 1997. Tectonic evolution of accretionary prism in the arc-continent collision terrane of Taiwan, Tectonophysics , 281, 31– 51. Google Scholar CrossRef Search ADS   Huang H.H., Shyu J.B.H., Wu Y.M., Chang C.H., Chen Y.G., 2012. Seismotectonics of northeastern Taiwan: kinematics of the transition from waning collision to subduction and postcollisional extension, J. geophys. Res ., 117, B01313, doi:10.1029/2011JB008852. Jolivet R., Grandin R., Lasserre C., Doin M.-P., Peltzer G., 2011. Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data, Geophys. Res. Lett. , 38, L17311, doi:10.1029/2011GL048757. Jolivet R. et al., 2012. Shallow creep on the Haiyuan Fault (Gansu, China) revealed by SAR Interferometry, J. geophys. Res ., 117, B06401, doi:10.1029/2011JB008732. Jolivet R., Agram P.S., Lin N.Y., Simons M., Doin M.-P., Peltzer G., Li Z., 2014. Improving InSAR geodesy using global atmospheric models, J. geophys. Res. , 119, 2324– 2341. Google Scholar CrossRef Search ADS   Jolivet R., Lasserre C., Doin M.-P., Peltzer G., Sun J., Rong D., Shen Z., Xu X., 2013. Creep spatial and temporal evolution revealed by InSAR time series analysis, ESA Spec. Publ. , 704, 90. Kang C.-C., Chang C.-P., Siame L., Lee J.-C., 2015. Present-day surface deformation and tectonic insights of the extensional Ilan Plain, NE Taiwan, J. Asian. Earth. Sci. , 105, 408– 417. Google Scholar CrossRef Search ADS   Kao H., Rau R.-J., 1999. Detailed structures of the subducted Philippine Sea plate beneath northeast Taiwan: a new type of double seismic zone, J. geophys. Res. , 104, 1015– 1033. Google Scholar CrossRef Search ADS   Kato T., Beavan J., Matsushima T., Kotake Y., Camacho J.T., Nakao S., 2003. Geodetic evidence of back-arc spreading in the Mariana Trough, Geophys. Res. Lett. , 30(12), 1625, doi:10.1029/2002GL016757. Kim K.H., Chiu J.M., Pujol J., Chen K.C., Huang B.S., Yeh Y.H., Shen P., 2005. Three-dimensional VP and VS structural models associated with the active subduction and collision tectonics in the Taiwan region, Geophys. J. Int. , 162, 204– 220. Google Scholar CrossRef Search ADS   Ku C.Y., Hsu S.K., 2009. Crustal structure and deformation at the northern Manila Trench between Taiwan and Luzon islands, Tectonophysics , 466, 229– 240. Google Scholar CrossRef Search ADS   Kuo-Chen H., Wu F., Chang W.-L, Chang C.-Y., Cheng C.Y., Hirata N., 2015. Is the Lishan fault of Taiwan active?, Tectonophysics , 661, 210– 214. Google Scholar CrossRef Search ADS   Lai K.-Y., Chen Y.-G., Wu Y.-M., Avouac J., Kuo Y.-T., Wang Y., Chang C.-H., Lin K.-C., 2009. The 2005 Ilan earthquake doublet and seismic crisis in northeastern Taiwan: evidence for dyke intrusion associated with on-land propagation of the Okinawa Trough, Geophys. J. Int. , 179, 678– 686. Google Scholar CrossRef Search ADS   Lee J.-C., Angelier J., Chu H.-T., 1997. Polyphase history and kinematics of a complex major fault zone in the northern Taiwan mountain belt: the Lishan Fault, Tectonophysics , 274, 97– 115. Google Scholar CrossRef Search ADS   Lee Y.-H., Byrne T., Wang W.-H., Lo W., Rau R.-J., Lu H.-Y., 2015. Simultaneous mountainbuilding in the Taiwan orogenic belt, Geology , 43, 451– 454. Google Scholar CrossRef Search ADS   Li Z., Fielding E.J., Cross P., 2009. Integration of InSAR time-series analysis and water-vapor correction for mapping postseismic motion after the 2003 Bam (Iran) earthquake, IEEE Trans. Geosci. Remote Sens. , 47, 3220– 3230. Google Scholar CrossRef Search ADS   Liang W.-T., Lee J.-C., Kuo B.-Y., 2005. Left-lateral strike-slip faulting in Ilan: lateral extrusion at the transition between Taiwan mountain range and the Okinawa trough, in Geodynamics and Environment in East Asia International Conference and 5th Taiwan-France Earth Science Symposium , pp. 104– 108, Taitung, ROC. Lin K.-C. et al., 2010. GPS crustal deformation, strain rate, and seismic activity after the 1999 Chi-Chi earthquake in Taiwan, J. geophys. Res ., 115, B07404, doi:10.1029/2009JB006417. Lin S.F., Liu P.M., Lai T.H., 2004. Late Holocene pollen sequence of the Ilan Plain, northeastern Taiwan, and its environmental and climatic implications, Terr. Atmos. Ocean. Sci ., 15, 221– 238. Google Scholar CrossRef Search ADS   Liu C.C., 1995. The Ilan plain and the southwestward extending Okinawa Trough, J. geol. Soc. (China) , 38, 229– 242. Lo C.-L., Chang E.T.-Y., Chao B.F., 2013. Relocating the historical 1951 Hualien earthquake in eastern Taiwan based on tide gauge record, Geophys. J. Int. , 192, 854– 860. Google Scholar CrossRef Search ADS   Nakamura M., 2004. Crustal deformation in the central and southern Ryukyu arc estimated from GPS data, Earth planet. Sci. Lett. , 217, 389– 398. Google Scholar CrossRef Search ADS   Okada Y., 1992. Surface deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 92, 1018– 1040. Prescott W., 1976. An extension of Frank's method for obtaining crustal shear strains from survey data, Bull. seism. Soc. Am. , 66, 1847– 1853. Rau R.-J., Ching K.-E., Hu J.-C., Lee J.-C., 2008. Crustal deformation and block kinematics in transition from collision to subduction: global positioning system measurements in northern Taiwan, 1995–2005, J. geophys. Res ., 113, B09404, doi:10.1029/2007JB005414. Remy D., Bonvalot S., Briole P., Murakami M., 2003. Accurate measurements of tropospheric effects in volcanic areas from SAR interferometry data: application to Sakurajima volcano (Japan), Earth planet. Sci. Lett. , 213, 299– 310. Google Scholar CrossRef Search ADS   Rosen P.A., Hensley S., Peltzer G., Simons M., 2004. Updated repeat orbit interferometry package released, EOS, Trans. Am. geophys. Un. , 85(5), 47, doi:10.1029/2004EO050004. Seno T., 1977. The instantaneous rotation vector of the Philippine Sea plate relative to the Eurasian plate, Tectonophysics , 42, 209– 226. Google Scholar CrossRef Search ADS   Shen Z.K., Jackson D.D., Ge B.X., 1996. Crustal deformation across and beyond the Los Angeles basin from geodetic measurements, J. geophys. Res ., 101, 27 957–27 980. Google Scholar CrossRef Search ADS   Shinjo R., Kato Y., 2000. Geochemical constraints on the origin of bimodal magmatism at the Okinawa Trough, an incipient back-arc basin, Lithos , 54, 117– 137. Google Scholar CrossRef Search ADS   Shyu J.B.H., Sieh K., Chen Y.G., Liu C.S., 2005. Neotectonic architecture of Taiwan and its implications for future large earthquakes, J. geophys. Res. , 110, B08402, doi:10.1029/2004JB003251. Sibuet J., Hsu S.K., 1997. Geodynamics of the Taiwan arc-arc collision, Tectonophysics , 274, 221– 251. Google Scholar CrossRef Search ADS   Sibuet J., Hsu S.K., 2004. How was Taiwan created?, Tectonophysics , 379, 159– 181. Google Scholar CrossRef Search ADS   Sibuet J. et al., 1987. Back arc extension in the Okinawa Trough, J. geophys. Res. , 92, 14 041–14 063. Google Scholar CrossRef Search ADS   Suppe J., 1984. Kinematics of arc-continent collision, flipping of subduction, and back-arc spreading near Taiwan, Geol. Soc. China Mem. , 6, 21– 33. Suppe J., 1987. The active Taiwan mountain belt, in The Anatomy of Mountain Ranges , pp. 277– 293, eds Schaer J.-P., Rodgers J., Princeton Univ. Press. Teng L.S., 1990. Geotectonic evolution of Late Cenozoic arc-continent collision in Taiwan, Tectonophysics , 183, 57– 76. Google Scholar CrossRef Search ADS   Teng L.S., Lee C.T., Peng C.H., Chen W.F., Chu C.J., 2001. Origin and geological evolution of the Taipei basin, northern Taiwan, West. Pac. Earth Sci ., 1, 115– 142. Teng L.S., Lin A.T., 2004. Cenozoic tectonics of the China continental margin: insights from Taiwan, Geol. Soc. Lond., Spec. Pub. , 226, 313– 332. Google Scholar CrossRef Search ADS   Teng L.S., Lee C.T., Tsai Y.B., Hsiao L.Y., 2000. Slab breakoff as a mechanism for flipping of subduction polarity in Taiwan, Geology , 28, 155– 158. Google Scholar CrossRef Search ADS   Tsai M.-C., Yu S.-B., Shin T.-C., Kuo K.-W., Leu P.-L., Chang C.-H., Mei-Yi H., 2015. Velocity field derived from Taiwan continuous GPS array (2007–2013), Terr. Atmos. Ocean. Sci ., 26, 527– 556. Google Scholar CrossRef Search ADS   Ustaszewski K., Wu Y.M., Suppe J., Huang H.H., Chang C.H., Carena S., 2012. Crust–mantle boundaries in the Taiwan–Luzon arc-continent collision system determined from local earthquake tomography and 1D models: implications for the mode of subduction polarity reversal, Tectonophysics , 578, 31– 49. Google Scholar CrossRef Search ADS   Wang K.L., Bilek S.L., 2014. Invited review paper: fault creep caused by subduction of rough seafloor relief, Tectonophysics , 610, 1– 24. Google Scholar CrossRef Search ADS   Wang K.L., Chung S.L., Chen C.H., Shinjo R., Yang T.F., Chen C.H., 1999. Post-collisional magmatism around northern Taiwan and its relation with opening of the Okinawa Trough, Tectonophysics , 308, 363– 376. Google Scholar CrossRef Search ADS   Wang H., Wright T.J., 2012. Satellite geodetic imaging reveals internal deformation of western Tibet, Geophys. Res. Lett. , 39, doi:10.1029/2012GL051222. Wang C.S., Yang M.L., Chou C.P., Chang Y.C., Lee C.S., 2000. Westward Extension of the Okinawa Trough at its western end in the Northern Taiwan area: bathymetric and seismological evidence, Terr. Atmos. Ocean. Sci. , 11, 459– 480. Google Scholar CrossRef Search ADS   Wang H., Wright T.J., Biggs J., 2009. Interseismic slip rate of the northwestern Xianshuihe fault from InSAR data, Geophys. Res. Lett ., 36, L03302, doi:10.1029/2008GL036560. Ward S.N., 1994. A multidisciplinary approach to seismic hazard in Southern California, Bull. seism. Soc. A ., 84, 1293– 1309. Welstead S., 1999. Fractal and Wavelet Image Compression Techniques , SPIE Optical Engineering Press. Wessel P., Smith W.H.F., 1998. New, improved version of generic mapping tools released, EOS, Trans. Am. geophys. Un. , 79(47), 579, doi:10.1029/98EO00426. Wessel P., Smith W.H.F., Scharroo R., Luis J., Wobbe F., 2013. Generic mapping tools: improved version released, EOS, Trans. Am. geophys. Un. , 94, 409– 410. Google Scholar CrossRef Search ADS   Wright T., Parsons B., Fielding E., 2001. Measurement of interseismic strain accumulation across the North Anatolian Fault by satellite radar interferometry, Geophys. Res. Lett ., 28, 2117– 2120. Google Scholar CrossRef Search ADS   Wright T.J., Parsons B.E., Lu Z., 2004. Toward mapping surface deformation in three dimensions using InSAR, Geophys. Res. Lett. , 31, L01607, doi:10.1029/2003GL018827. Wu Y.-M., Zhao L., Chang C.-H., Hsu Y.-J., 2008. Focal-mechanism determination in Taiwan by genetic algorithm, Bull. seism. Soc. Am. , 98, 651– 661. Google Scholar CrossRef Search ADS   Wu Y.-M., Hsu Y.J., Chang C.H., Teng L.S., Nakamura M., 2010. Temporal and spatial variation of stress field in Taiwan from 1991 to 2007: insights from comprehensive first motion focal mechanism catalog, Earth planet. Sci. Lett. , 298, 306– 316. Google Scholar CrossRef Search ADS   Yang M., Chen K.H., Shiao S.W., 2003. A new height reference network in Taiwan, Surv. Rev ., 37, 260– 268. Google Scholar CrossRef Search ADS   Yu S.-B., Kuo L.-C., 2001. Present-day crustal motion along the longitudinal Valley Fault, eastern Taiwan, Tectonophysics , 333, 199– 217. Google Scholar CrossRef Search ADS   Yu S.-B., Chen H.-Y., Kuo L.-C., 1997. Velocity field of GPS stations in the Taiwan area, Tectonophysics , 274, 41– 59. Google Scholar CrossRef Search ADS   Yu S.-B., Kuo L.-C., Punongbayan R.S., Ramos E.G., 1999. GPS observation of crustal deformation in the Taiwan-Luzon region, Geophys. Res. Lett ., 26, 923– 926. Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. APPENDIX S.1 Is the long-term deformation of the Ilan Plain controlled by tectonic or non-tectonic signals? Figure S1. Comparison between long-term (2003–2009) surface deformation from SBAS-PSI and groundwater-level (2005 January–2009 June) changes. Triangles represent a long-term rise in the groundwater level, and inverted triangles represent a long-term decline in the groundwater level. Triangles/inverted triangles outlined in red denote that pumping activity was observed in the aquifers during the monitoring period. Triangles/inverted triangles outlined in black denote that no pumping activity was observed. Figure S2. Long-term groundwater-level changes at the (a) WJ and (b) ZX stations. Each station is equipped with water loggers at four different depths. Red dashed lines are least-squares polynomial fittings to the first-degree polynomial and represent the long-term groundwater-level changes. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations