Tectonic implications on the 2018 Hualien Earthquake

The 2018 Hualien earthquake occurred at the junction of the deformed conti-nental crust of Eurasian plate and the Heping sea basin, where a northeast trending seismic belt exists. In this study, a 3D seismic tomographic inversion is used to investigate the seismogenic structures of the 2018 Hualien earthquake sequence. An earthquake relocation procedure is performed and the focal mechanisms are analyzed to study the faulting behavior. Our results indicate that the source region of the 2018 Hualien earthquake consists of complex high-angle eastward and westward dipping reverse faulting. We also observe that most earthquakes that have occurred in the area exhibit considerable variation in Vp/Vs ratios. Fluid migration may have played an essential role in causing the Vp/Vs ratio variation. We suggest that the mainshock of the 2018 Hualien earthquake may be associated with a west-dipping fault, which is a blind fault tending toward the Heping sea basin and possibly belonging to the Central Range fault system. (1999) for northeastern Taiwan and Cheng et al. for eastern Taiwan. These studies suggest that, in the basement


INTRODUCTION
As Taiwan is located in the circum-Pacific seismic belt and is frequently affected by earthquakes, a large amount of seismic data is available, and we can take advantage of this wealth of data to gain a better understanding of the tectonic structures in eastern Taiwan. In this region, most of the earthquakes have occurred at shallow depths, so knowledge of the high potential of earthquake hazards in this region is a critical issue. Recently, several disastrous earthquakes have struck eastern Taiwan, for example, a moderate earthquake (M L = 6.4; depth = 15 km) occurred in Ruisui, Hulien County on 31 October 2013, and seven months later, close to the location of this event, another moderate earthquake (M L = 5.9, depth = 18 km) occurred in Fenglin, Hualien County, on 21 May 2014. In early February 2018, a seismic swarm comprising hundreds of small earthquakes occurred off-shore of Heping area, where the collision boundary of eastern Taiwan is located (Fig. 1, yellow dots). The largest magnitude among this swarm is up to M L 5.7. On 6 February 2018, another moderate earthquake (M L = 6.2; depth = 6 km) occurred off-shore of Heping area and caused severe damage, including buildings collapse, bridge destruction, surface rupture and causing casualties in Hualien City, which was previously affected by the 1951 Hulien-Taitung earthquake ). This devastating earthquake raised important questions, including whether the event was triggered by the abovementioned seismic swarm, whether the faulting behavior of this earthquake involved multiple ruptures, and whether the event triggered known fault systems. Therefore, the study of seismogenic structures in this area is a very important issue.
Seismic swarms are generally thought to be caused by fluids invading the seismogenic zone and hence their behaviors pattern differs significantly from that of an earthquake sequence. Hence, the main goal of this study is to investigate the seismogenic structure of the 2018 Hualien earthquake and its tectonic implications beneath the northern part of the Longitudinal Valley in east Taiwan. Several regional-scale seismic tomographic studies in Taiwan have been reported, such as Liaw and Yeh (1983), Yeh et al. (1989), Kao and Rau (1999) for northeastern Taiwan and Cheng et al. (2002) for eastern Taiwan. These studies suggest that, in the basement of the Longitudinal Valley and Coastal Ranges, there exists a narrow stripe-like high-velocity zone that belongs to the Luzon arc and may extend toward 24°N. Although several 3-D models of Vp or Vp/Vs ratio have been proposed based on seismic travel-time data in the Taiwan area (Kim et al. 2005;Wu et al. 2007;Kuo-Chen et al. 2012;Huang et al. 2014), they only provide large-scale velocity structure mappings, and their overall resolution is insufficient to describe velocity variations within a small area. Therefore, we intend to invert a high resolution mapping of 3-D velocity structures beneath the northern Longitudinal Valley, eastern Taiwan, which would provide better insights into the seismogenic structure of the 2018 Hualien earthquake.

BACKGROUND GEOLOGICAL SETTING
The terrain distribution in eastern Taiwan can be divided into three geological units: from west to east, the Central Mountain Range (CMR), the Longitudinal Valley (LV), and the Coastal Ranges (CR). The CMR is located in the extreme west of our study area; it has a terrain height of 1000 -2000 m and is composed of metamorphic rock. The LV is located between the CMR and the CR and is a flat valley approximately 150 km long. The LV is segmented mainly by different types of faulting behaviors as well as the accumulation of alluvial deposits from the CMR. The ongoing vigorous orogenic processes in Taiwan result in active tectonic processes and frequent earthquakes. The CR is located above the Philippine Sea plate; it belongs to the Luzon arc and fills up with fore-arc sediments. In the northern segment of the LV, most earthquakes occur on plate boundaries with reverse faulting mechanisms.
The LV is a plate suture zone. As the Philippine Sea Plate moves at a rate of 8 cm yr -1 toward the Eurasian plate in the northwest direction (Yu et al. 1997), violent crustal deformation occurs in the LV, forming the Longitudinal Valley Fault (LVF) system (Shyu et al. 2006a;Thomas et al. 2014). The main LVF is a thrust fault and dips eastward (see Fig. 1). The exposed fault traces are roughly along the foothills of the western part of the CR, which can be divided into several segments from the Hualien to Taitung areas (Shyu et al. 2006a;Thomas et al. 2014). Two major active faults in the Hualien area are the Linding Fault (LDF) (Lin et al. 2009) and the Milun Fault (MLF), which may be considered the northernmost extension of the LVF (see Fig. 1). Shyu et al. (2007) pointed out that the LVF is segmented by different slip behaviors. Another major fault in the LV is the Central Range Fault (CRF) (Shyu et al. 2006b). However, the fault trace of the CRF is not obvious and its position is still controversial. In addition, the CRF is considered to represent the boundary between the CMR and the LV. The geometry of the boundary between the CMR and the LV from previous studies is quite complex and slightly different from south to north; this boundary does not always correspond to the foothills of the eastern edge of the CMR. Geological and geophysical surveys also indicated that the CRF may be a blind fault, especially in the northern part of the LV (Shyu et al. 2006b).
The northern part of the LV experiences high seismicity (Kuo-Chen et al. 2004;Chuang et al. 2014;Thomas et al. 2014) with most of the seismic activities being caused by the LVF. A few seismic activities have occurred in the CRF, such as the 2013 Ruisui earthquake (Chuang et al. 2014;Wen et al. 2016) and the 2014 Fenglin earthquake, whose epicenters were on the eastern edge of the CMR. The distribution of the aftershocks shows that the seismogenic fault is a west-dipping fault. Furthermore, the aftershocks are only distributed to a depth up to about 5 km, which implies that the CRF is still in a locked state above the depth of 5 km (Shyu et al. 2006b).

3D LOCAL TOMOGRAPHY MODEL
In this study, we used travel time data of P and S waves with a local magnitude greater than 2 recorded by the Central Weather Bureau Seismic Network (CWBSN) and the Broadband Array in Taiwan for Seismology (BATS) for the period 2012 -2017 (see Table 1 for the station list). A damped least-square inversion method (Thurber 1983;Eberhart-Phillips and Michael 1998) was used to determine the Vp structures and the Vp/Vs ratios in the crust and the upper mantle of the region. The study area extends from 23.7 -24.4°N in latitude and 120.9 -122°E in longitude. Many seismic events occurred in the region of interest (Fig. 1). In order to improve the resolution of tomographic inversion and to obtain a uniform ray distribution in all directions, we chose events with epicenter location errors of less than 5 km and those with more than six readings for the P and S waves. Under these criteria, we were able to select 10309 events with 105705 P-wave and 94279 S-wave arrival times.
We chose a priori 1-D P-velocity model modified from Chen et al. (1994) and Wu et al. (2007), following the approach of Kissling et al. (1994), to obtain an initial reference model (the minimum 1-D model), which can roughly represents the main features of a known velocity structure. We further parameterized the model by using horizontal layers of constant velocities, which capture the main features of the velocity structure (Table 2). We parameterized the 3D structure by forming even 2-D grids through a trial-and-error pro-cess (see Fig. 2a). The velocity inversion was damped carefully by calculating a range of damping values to determine the tradeoff between the data residuals and model variances (see the Supplement, Fig. S1). Thus, the final damping values chosen for inverting Vp and Vp/Vs ratios were 20 and 10, respectively. We used a horizontal grid size of 6 km × 6 km in this study. It is worth mentioning that, in our inversion, the effect of station elevation on the 3 D tomographic inversion was also taken into account so that the seismic rays can be traced exactly to the true station position.

Model Examination
We examined our final inversion results for various model parameters by using checkerboard test and derivative weight sum test (Toomey and Foulger 1989), and present the outcomes in this section.  Table 2. P-wave velocity and Vp/Vs ratio of 1D initial model.

Checkerboard Test
We used a checkerboard test to illustrate the resolution of the 3D tomographic model. First, the 1D velocity model (Table 2) was perturbed by ±4% as the initial velocity model and the theoretical traveling time data was calculated. Then, the synthetic data was inverted to obtain the checkerboard structure. Because only a few stations were located in the oceanic area, the density of the rays across the oceanic area was relatively low resulting in a low spatial resolution of the tomographic inversion for the oceanic area. As the distributions of the seismic events on land were clustered, high-resolution inversion results were obtained at the LVF and in a ~15 km-wide band along the sides of the LVF (see Figs. 2a and 3a). However, the resolution at depth 0 -2 and 25 -35 km was relatively low. The low resolution for depths of 0 -2 km may be because the incident angles were almost perpendicular to the surface, which lowered the lateral resolutions of the velocity. For depths greater than 25 km, the majority of the seismic events occurred in oceanic region and only a few were selected under our criteria, leading to a low resolution. For the deeper part of the crust or in the upper mantle, low seismicity also resulted in a low ray density, thus producing a poor resolution. Nevertheless, the overall resolution in each layer was good. The high-resolution area indicates that the seismic rays crossed most of the crustal region and can be reliable used to determine the geological structures beneath the northern LVF area more precisely.

Derivative Weighted Sum
In addition to the checkerboard test, in order to ensure the confidence level of our model, we used the derivative weighted sum (DWS) was first proposed by Toomey and Foulger (1989) as a parameter to estimate the seismic ray density. Small DWS values indicate poorly sampled nodes (Toomey and Foulger 1989). In Figs. 2b and 3b, the areas enclosed by white curves are those with a ray density greater than 50. At the depths 5 -25 km, the DWS is high, indicating a high confidence level for this area. This phenomenon is obvious in the vicinity of the MLF and LDF and especially in the northern part of the LV.

RELOCATIONS AND FAULT PLANE SOLUTIONS
After tomographic inversion, an earthquake relocation procedure was performed using the velocity structure derived from our inversion results. In Fig. 4, the spatial distribution of relocated events is denoted by circles, with their size scaled according to the earthquake magnitudes and their color representing the focal depth, scaled by a grey scale   (Fig. 4). These relocated earthquakes yielded a significantly lower total travel-time misfit with an average RMS residual of 0.236 s (the original RMS residual is 0.607 s). The average horizontal shifts of their epicenters was 3.53 km and that of their focal depths was 2.37 km. However, the epicenters of relocated events that occurred in the oceanic region tended to shift toward land, which is reasonable as there is a lack of seismic stations in the ocean.
To understand the rupture behavior beneath the fault zones, we further used the recalculated azimuths and takeoff angles to obtain the focal mechanism before 2018. Figure 5a shows the fault plane solutions in the study area; approximately 700 fault plane solutions were obtained, where the normal and strike-slip faulting types were distributed in the CMR. To improve the accuracy of the analysis, we calculated the take-off angles and azimuths from the obtained 3 D tomography structure. With the polarities of P wave arrivals (Reasenberg and Oppenheimer 1985), we were able to obtain more than 700 focal mechanisms with magnitudes mostly between 3 and 5. The results show several earthquake clusters in our study area. The distribution of relocated seismicity shows a distinct strip-shaped seismic zone near the CMR in the LV (see Fig. 4), with a predominantly reverse-faulting rupture mechanism. Furthermore, another obvious stripe-like seismic zone was observed near the CR, where the predominant rupture mechanism is reverse faulting and strike-slip faulting. In particular, in the vicinity of the MLF and LDF, shallow and small earthquakes occur frequently as reported by the Central Weather Bureau (https://www.cwb.gov.tw/V7/earthquake/rtd_eq.htm). Along the fault to the oceanic side and close to off-shore of the Heping area (see the left upper inset in Fig. 1), another seismic cluster was observed, with thrust faulting as the dominant mechanism (see Fig. 5a). In the northwest of the study area, another group of shallow and small earthquakes with a predominant normal faulting rupture mechanism was observed (Fig. 5a).
For the 2018 Hualien earthquake source area, the relocated seismic events can be grouped into two zones (marked by H1 and H2 in Fig. 4) for further discussion. For a more systematic comparison, triangular charts that specify the distribution of focal mechanism (Figs. 5b and c) were analyzed (Frohlich 1992(Frohlich , 2001. The results show that a greater number of events had a thrust-like faulting regime than the strike-slip faulting. The triangular diagram for the H1 zone (Fig. 5b) shows that the seismic events were of the strikeslip (~30%) and thrust (~67%) rupture types. The H2 zone is the main location of the 2018 Hualien earthquake sequence, which the pattern for the H2 zone (see Fig. 5c, strike-slip type: ~27%; thrust type: ~63%) is also quite similar to that of the H1 zone. Therefore, the analysis of the past seismic data is still dominated by the thrust faulting pattern in the two zones, which is consistent with the previous studies (Wu et al. 2008;Hsu et al. 2010).

RESULTS AND DISCUSSION
This study aims to map the seismogenic structures before the 2018 Hualien earthquake. In this section, we will discuss the results in two parts. In the first part, we will discuss the results of deriving the velocity structures and Vp/ Vs ratios. The epicenters were relocated and the velocity structures were derived through an iterative process. In the second part, we will discuss cross-sectional examination of the obtained velocity structures. In order to understand the velocity structures beneath the Hualien area, three cross sections were examined (Fig. 4). Two profiles almost perpendicular to the MLF, and one parallel to the MLF and LDF were observed. Thus, we will outline the relationships between the velocity structures, fault zones and seismicity.
As shown in Fig. 2, although numerous seismic events exist in each layer, most of them are located in a depth range of 2 -25 km. Therefore, in our inversion process, the seismic rays pass through most of the nodes. At depths of 0 -5 km, high and low Vp values show a belted distribution. In addition, low Vp anomalies are distributed in the vicinity of the fault zone (see Fig. 2c). At depths of 5 -15 km, the high Vp anomalies increase with depth and expand in a southeast. This may be related to the existence of oceanic crust structures beneath this area. Several cluster events existing in each layer and mostly show low Vp anomalies. In the vicinity of the fault zones, a high Vp/Vs ratio anomaly is due to the increasing pore pressure, which leads to a decrease in the S wave velocity. We also observed that a low Vp/Vs ratio anomaly broadens from north to south between 2 and 15 km and that the low anomaly zone extends to the CR with increasing depth (see Figs. 2c and 3c). This phenomenon can be attributed to the older and denser rock formations beneath the CR, containing less fluid or SiO 2 , which will in turn result in lower Vp/Vs ratios. Our results indicate that most seismic events are located in areas that have high or greatly varying Vp/Vs ratios. Figure 6a shows high and low anomaly zones in the Vp, Vs and Vp/Vs profiles in profile AA'. It can be observed that the velocity structures at the northern end of the LV are inclined to the north. In addition, a high Vp and Vs anomaly zone can be observed at the southern end of the AA' profile. The aftershock sequence of the Hualien earthquake (the preliminary locations which are indicated as the black dots in Fig. 6a) stops rupturing at the black arrow indicated in Fig. 6a. This feature needs further investigation in our future work. Figures 6b and c show the anomalies in the Vp, Vs and Vp/Vs profiles in profile BB' and CC'. They indicate an eastward leaning fault geometry beneath the LVF zone at a depth of 5 km, which implies that the fault may be an extension of the MLF; this can be considered as a passive rupture plane from the 2018 Hualien earthquake source inversion (Lee et al. 2018;Wen et al. 2019). In particular, in profile CC' (Fig. 6c), a low-velocity anomaly zone is observed in the shallow layer below the Milun Platform, where the earthquake cluster tends eastward. On the west side of the fault and at depths greater than 5 km, a high Vp with high Vp/Vs ratio anomaly area is observed, where an earthquake cluster tends to dip toward the west. The westdipping seismic zone may belong to the CRF and exhibits anomalous seismic activity, possibly owing to the high pore pressure, which can reduce the pressure to rupture (Ougier-Simonin and Zhu 2015). The abovementioned features were also observed from the velocity images of the eastern subduction zone obtained from Huang et al. (2014). Their results indicate an east-dipping seismic zone is in the northern part of the LV, and at depth greater than 10 km, there also exists a seismic zone inclined toward to west [see Huang et al. (2014), profile DD' in Fig. 5]. In addition, both seismic zones are located at the boundary between high-and lowvelocity anomaly. Therefore, this anomaly area can be interpreted as an oversaturated pore pressure zone. These results also imply that the rupture behavior may be caused by the occurrence of localized supra-lithostatic fluid pressure on the steeply inclined thrust fault due to the inefficiency of hydro-fracturing (i.e., low friction coefficient) (Carminati et al. 2004). We conclude that high seismicity exists in zones that exhibit a low P wave velocity and where the gradients of the Vp/Vs ratio vary greatly in each profile. From profiles BB' and CC' (Figs. 6b and c), a west dipping seismic zone is clearly observed. The geometry of the seismic anomaly zone can also be identified as the main rupture plane of the 2018 Hualien earthquake as considered in the finite-fault model (Lee et al. 2018;Wen et al. 2019). A low rupture velocity across the Milun fault (e.g., Lee et al. 2018;Hwang et al. 2019) is possibly related to low velocity structures in the Milun Tableland (Fig. 6). Such low-velocity structures on land might be a factor contributing to the generation of pulse-like velocity waveforms at strong-motion stations located at the southernmost tip of the Milun fault (Kuo et al. 2019;Wu et al. 2019). In addition, low-velocity structures down to a depth of 10 km could restrain the rupture energy from growing into deeper regions. The teleseismic P-wave back-projection proposed by Jian et al. (2019) suggested that the rupture energy of the Hualien earthquake occurred above a depth of 10 km. Nakajima et al. (2001) and Takei (2002) concluded that the variations in Vp/Vs ratios may be due to water-or gas-filled cracks with high aspect ratios. They derived a coefficient, dlnVs/dlnVp, which represents a fractional change in Vp and Vs and can be used to determine whether the variation in Vp/Vs ratios is caused by gas, water, or magma at a specific area (Takei 2002). Figure 7 shows that dlnVs/dlnVp is higher than 2 in the region with high seismicity and the anomaly areas, which indicates the presence of water-filled or high-density cracks. Therefore, results obtained from the tomography inversion at shallower fault segments show low Vs and high Vp/Vs anomalies, which may imply that those zones are filled with fluid at crack openings (Figs. 6 and 7). This phenomenon is clearly observed in the fault zone in the northern part of the LV.

CONCLUSIONS
In this study, we applied a damped least-square inversion method to investigate the 3D Vp structures and Vp/ Vs ratios of the crust and upper mantle beneath the Hualien area, eastern Taiwan, by using the travel-time data of seismic waves. By using the time difference between the observed P and S waves, we inverted the Vp/Vs ratios. From the results, we could not only relocate earthquakes but also deduce the relationship between seismicity and the regional geological structures. In addition, we found that most earthquakes occurred in areas where the Vp/Vs ratio gradients varied greatly. We also suggested that, instead of breaking the MLF, the mainshock of the 2018 Hualien earthquake may be associated with a west-dipping fault that could reach out to the Heping basin. Moreover, we concluded that the 2018 Hualien earthquake sequence occurred along a westdipping fault which is a blind fault and may belong to the CRF. The 3D velocity model derived from this study will help improve the accuracy of earthquake location and is an important factor in estimating strong ground motion. Therefore, the results obtained for the Hualien area using this model are of great significance in the prevention and mitigation of earthquake hazards.