South Ilan Plain High-Re s olution 3-D S-Wave Velocity from Ambient Noise Tomography

The Ilan Plain in northeastern Taiwan is located at a pivotal point where the Ryukyu trench subduction zone, the northern Taiwan crustal stretching zone, and the ongoing arc-continent collision zone converge. In contrast to the North Ilan Plain, the South Ilan Plain exhibits a thin unconsolidated sedimentary layer with depths ranging from 0 1 km, high on-land seismicity and significant SE movements relative to Penghu island. We deployed a dense network of 43 short-period vertical component Texan instruments from June to November 2013 in this study, covering most of the South Ilan Plain and its vicinity. We then used the ambient noise tomography method for simultaneous phase and group Rayleigh wave velocity measurements to invert a high-resolution 3-D S-wave for shallow structures (up to a depth of 2.5 km) in the South Ilan Plain. We used the fast marching method for ray tracing to deal with ray bending in an inhomogeneous medium. The resulting rays gradually bend toward high velocity zones with increasing number of iterations. The high velocity zone results are modified by more iterations and the resolutions become higher because ray crossings are proportional to ray densities for evenly distributed stations. The final results agreed well with known sedimentary basement thickness patterns. We observed nearly EW trending fast anomalies beneath the mountainous terrain abutting to the South Ilan Plain. The Chingshui location consistently exhibited a low S-wave velocity zone to a depth of 1.5 km.


INTRoDucTIoN
Ambient noise refers to noise signals that originate from sources that are randomly distributed in space.About a decade ago, theoretical derivations demonstrated that the cross-correlation of continuous ambient noise between two stations predominantly yields fundamental-mode Surface wave Green's functions (e.g., Lobkis and Weaver 2001;Weaver andLobkis 2001, 2004;Derode et al. 2003;Snieder 2004;Wapenaar 2004;Larose et al. 2005;Snieder and Wapenaar 2010).As a result, the ambient noise tomography method has been widely used to investigate crust and mantle velocity structures (e.g., Lin et al. 2007;Huang et al. 2012Huang et al. , 2015)).In comparison with the conventional approach using earthquake sources (e.g., Yeh et al. 2016), the merits of this application include that it is free from earthquake source distribution limitations and the resolution depends solely on the station distribution.If dense array data are available, the ambient noise tomography resolutions can be sufficiently high to reveal local geological structures such as sedimentary basins and underground heat resources.We aim to resolve high-resolution S-wave crust velocity structures beneath the South Ilan Plain using the ambient noise tomography method applied to local dense array data.Where Taiwan is located on the collision boundary between the Eurasian and Philippine Sea plates, the northeast Ilan Plain forms a triangular alluvial delta that most likely represents the southwestern end of the Okinawa Trough rifting zone.The Ilan Plain is roughly divided into northern and southern areas by the Lanyang River, which originates from the western vertex and continuously accumulates sediments on the plain.The tectonic characteristics are distinctive between the two areas, i.e., the southern area exhibits a relatively thin top sedimentary layer, high on-land seismicity, and significant clockwise block rotation.The South Ilan Plain has been targeted for past and current thermal explorations.Therefore, we focused on the South Ilan Plain and deployed a dense array of 43 Texan instruments to study ambient noise tomography.We used high resolution ambient noise tomography to better understand the seismogenic and geothermal structures of the South Ilan Plain.

TecToNIc SeTTINgS
Taiwan sits atop the boundary between the Eurasian and Philippine Sea plates, consisting of mountain belts formed by the collision of the China continental shelf and the Luzon arc since the late Cenozoic period (Ho 1988).In northern Taiwan, however, the previous collisional orogeny is currently undergoing crustal stretching as a result of flipped subduction polarity from the northwest vergent Manila trench to the south vergent Ryukyu trench since Plio-Pleistocene times (Teng 1996).The most recent rifting phase (2 MA) of the Okinawa Trough, a back arc basin behind the Ryukyu trench-arc system that developed on attenuated continental lithosphere (Sibuet et al. 1987), is most likely associated with the Philippine Sea plate subduction.
The Ilan Plain located in northeastern Taiwan unfolds toward the Okinawa Trough, suggesting onshore rifting of the southwest end of the Okinawa Trough (Hou et al. 2009).The Ilan Plain is located at the pivotal point where collisional orogeny, to the southwest, and crustal stretching, to the northeast, converge (Fig. 1).
To the northwest the Ilan Plain abuts the Eocene-Oligocene formations of the Hsuehshan range and to the south, the Miocene Lishan formation of the Central Range (Ho 1986).The Hsuehshan range formations have older sediments, higher metamorphic grades, and more symmetrical strain deformations (Lee et al. 1997).The Lishan Fault, which lies between the Hsuehshan and Central ranges, is a major fault zone in the Taiwan collision belt (Lee et al. 1997) and seemingly follows the Lanyang River extending to the Ilan Plain (Sibuet et al. 1987).Intriguingly, the north and south sides of the Ilan Plain, which are divided by the Lanyang River, exhibit distinctive tectonic features.
Based on seismic reflection method results the depths of the unconsolidated Quaternary alluvium on the Ilan Plain display a concave shape dipping shoreward with the deepest area on the North Ilan Plain, i.e., approximately 1500 m, while that on the south side is shallower and reaches approximately 750 m near the Lanyang River (Chiang 1976).The overall shape is consistent with that derived from the Sp wave delay times converted from the basement (Chang et al. 2010).Recent GPS observations revealed that the central and southern flanks of the Ilan Plain are currently undergoing extensional strain rates explained by the lateral extrusion facilitated by the Okinawa Trough rifting, while the strain rates of the northeastern flank are only minor or insignificant (Hou et al. 2009).As a result, on-land seismic activities on the Ilan Plain are considerably higher in the south than those in the north.The most significant activities in the past decade are the M w 6.1 earthquake that occurred on 15 May 2002, and the M w 5.7 doublet earthquakes that occurred on 3 March 2005, with predominantly left lateral strike-slip focal mechanisms on EW trending faults, which were attributed to the lateral extrusion (Fig. 1) (Ku et al. 2009;Lai et al. 2009;Shyu et al. 2016).The South Ilan Plain is also known for previous geothermal explorations (Chingshui) that are expected to continue in the future.It has been proposed that the observed EW trending magnetic anomalies crossing the mid-Ilan Plain may be caused by the acquired magnetisms of the upward magma body when cooled below the Curie temperature.These magma dikes are the likely heat sources (Tong et al. 2008).

DATA AND MeTHoDS
We deployed a dense network of 43 short-period vertical component Texan instruments from June to November 2013, covering most of the South Ilan Plain and its vicinity (Fig. 1).The interval between stations is approximately 2 km.We used the ambient noise tomography method (You et al. 2010;Huang et al. 2012) to obtain a high-resolution 3-D upper crust S-wave velocity model of the South Ilan Plain.There are four main parts to the procedure: (1) preparing data for a single station, (2) extracting the empirical Green's function (EGF) for each station pair by cross-correlation, (3) measuring group and phase velocities as a function of frequency (Rayleigh wave dispersion) from the EGF, and (4) jointly inverting 3-D S-wave velocity models.For the inversion scheme we applied a fast marching method (Rawlinson and Sambridge 2004) for 2-D ray tracing to account for ray bending in an inhomogeneous medium.We provide detailed descriptions of each part of the procedure in the following sections.
The data collected from the field were stored as oneday periods for each file.For each one-day-period data we removed the daily trends, obtained the mean and then deconvolved the instrument response.After preparing the data for each station we calculated the cross-correlation function (CCF) of each station pair to represent their corresponding EGF.The CCF was calculated in the frequency domain and spectral whitening was applied to reduce the resonance effects at certain frequencies.The CCF signal-to-noise ratios were enhanced using Welch's method (Seats et al. 2012), where the final CCF was derived by stacking the daily record results from all CCFs and processing these using a moving window with partial overlap (Fig. 2).The windows with earthquake and/or anomaly signals were omitted for CCF calculations.This was judged by setting a threshold greater than 1.3 times the standard deviation of the amplitude distributions from a particular day.For each station pair we further rejected those daily CCFs that did not resemble most of the others by calculating the cross-correlation coefficient between a particular CCF and the reference CCF determined by stacking all CCFs.Those with a coefficient less than 0.8 were discarded and we stacked the remaining CCFs to represent the EGF for that station pair.In Fig. 3, CCFs of each station pair are stacked, then we stacked once again judging by similar interstation distance despite there are not in the same station pair, arranged with inter-station distance.The stacking CCFs are filtered at the 0.5 -3.7 s period band.The EGF of a station pair is obtained from folding the CCF causal and acausal parts.The EGF represents waveforms excited by a unit impulse source at one station and received at the other and is a function of the velocity structure between the two stations.For vertical component ambient noise observations the EGF is dominated by fundamental Surface waves.We used frequency-time analysis (Levshin et al. 1989) and an image transformation technique (Yao et al. 2005(Yao et al. , 2006) ) to measure the EGF group and phase velocities as a function of frequency.The frequency-time analysis determines group velocities from the arrival time envelopes after applying a Gaussian filter.On the other hand, the image transformation technique constructs a time-period (t-T) image for the surface wave component.Each image column is the normalized amplitude corresponding to the filtered central period (T), and an interval of 0.2 s was used in this study.The phase velocity (c) was converted from time (t) using the far-field representation of the Green's function for the surface-wave fundamental mode (Aki and Richards 2002;Yao et al. 2006).The 2π ambiguity of the phase velocity for one period was calculated from previous group velocity results based on the assumption that velocity increases with the time period.To satisfy far-field conditions the maximum period used was limited by an inter-station distance that was less than 2.5 times the corresponding wavelength.While three (Bensen et al. 2007) and two (Shapiro et al. 2005) times the wavelength were previously used, we chose the 2.5 times wavelength criteria and the choice is validated by the measured velocity smoothly varying with periods.As a result, the dispersion curves were measured over a period ranging from 0.5 -3.7 s.
Following group and phase velocities dispersion determination for each station pair, we inverted the 3-D Swave structures using wavelet-based sparsely-constrained tomography (Fang et al. 2015).In addition to the flexible regularization based on ray densities, the method also deals with ray bending in an inhomogeneous medium, which was traced using the fast marching method (Rawlinson and Sambridge 2004).In the inversion scheme the partial derivatives of the group and phase velocities using S-wave velocity as a function of depth provide a sensitivity kernel for each period and produce velocity constraints on the depth dimension.
We performed ten iterations for the final velocity structures using an updated velocity model, ray paths, and sensitivity kernels after each iteration.Figure 4 shows that the mean residuals reduction is minor after eight iterations with values less than 0.5 s.We first performed a checker board test to assess the inversion resolution by alternatively assigning positive and negative 20% S-wave velocity anomalies to adjacent grids.The grid sizes used are 0.01° × 0.01° × 0.5 km and 0.02° × 0.02° × 0.5 km, respectively (Figs. 5 and 6).Three vertical synthetic and recovered checkerboard profiles are shown in Fig. 7.

ReSulTS
Although the objective checkerboard velocity structures were highly inhomogeneous in an orderly manner, most of the ray paths between station pairs passed through several alternating velocity anomaly cycles and the inhomogeneity effects were averaged over the entire path.As a result the resolutions are largely constrained by the station distribution.The results showed that the original grids were well recovered up to 1 km depth for 0.01° lateral size (Fig. 5) and up to 1.5 km depth for 0.02° (Fig. 6).The resolutions deteriorated at greater depth because we used short-period (0.5 -3.7 s) signal portions, which are more sensitive to the shallow crust.We thus focus on structures no deeper than 1.5 km in this study (Fig. 7).The ability to resolve highresolution shallow structures is important to studies on sedimentary thickness, seismic wave propagations and thermal explorations.
With respect to the 3-D S-wave velocity structure inversion, because we updated the ray paths after each iteration according to the newly derived velocity model, the final ray paths after ten iterations bend toward the high velocity zones (Fig. 8).Notably in the NE part of the South Ilan Plain, the 0.5-s ray paths bent significantly while the 2.5-s ray path bend was only relatively minor (Fig. 8).This can be attributed to the differences in the depth sensitivity kernelthe longer period over which deeper structures were sampled-and can be explained by the high velocity anomalies that only exist at the top and do not extend to deeper parts where the 2.5-s ray paths were sampled.
We present the inverted 3-D S-wave structures of the South Ilan Plain as six horizontal slices at various depths ranging from 0 -2.5 km at a 0.5-km interval, with overlapping contour lines at the basement depths (Fig. 9).We focused on the depth results from 0 -1.5 km slices 1 to 4 because of the lack of resolution at lower depths (Figs. 5 and 6).Slices 1 to 4 exhibited a northeastward decrease in the S-wave velocity.The NE S-wave velocity of the South Ilan Plain had the lowest value.This is consistent with the northeastward increase in sediment thickness (Fig. 9) suggesting that sediment layers are the main contributor to shallow structure Swave velocity.In slice 2 the S-wave velocities increased from 1.0 -2.0 km s -1 (red to green) and do not correspond exactly to the 0.5-km basement depth contour.This is because depth kernels are not clear cut outside the 0.5 km depth and the low velocity sediment above will still contribute.A nearly EW trending high velocity feature was evident at a depth of 0.5 km (between Sanshin and Suao in slice 2 of Fig. 9), which was sustained to depths greater than 1.5 km.The fast anomalies exhibited distinct features on two sides.At 0.5 km depth the west one is relatively far away from the basement contour while the east one exhibits stark velocity contrast (slice 2 of Fig. 9).Extending to greater depth at 1.5 km the west one displays a northward dipping feature while the east one a vertically dipping feature (Figs. 9 and 10).We speculate the reason being relative warm crustal material above 1.5 km to the west.In retrospect the fast anomalies resulted in ray bending and focusing (Fig. 8).The fact that the 2.5-s rays exhibited only minor bending agrees with the anomalies being smeared at depth (slice 6 of Fig. 9).The Chingshui location consistently exhibited a low S-wave velocity zone to a depth of 1.5 km, which corresponds well with its history of thermal explorations (Tong et al. 2008).Further investigations are required to determine whether the thermal source extends to greater depths.

DIScuSSIoN
The application of ambient noise tomography, together with dense array data from the South Ilan Plain, enabled us to derive a high-resolution 3-D S-wave velocity model for near surface structures.According to the checker board test results (Figs. 5 and 6), the resolution can reach as high as 1.0 km up to a depth of 1.0 km and as high as 2.0 km up to a depth of 1.5 km.The high-resolution results for the S-wave velocity model exhibit both low and high velocity anomalies.The low velocity anomalies above 1.0 km on NE of the South Ilan Plain are in agreement with a previously established sedimentary layer basement.The EW trending fast anomalies below 0.5 km are distributed at Central Mountain Range locations abutting the South Ilan Plain.
Fast marching method adoption to deal with ray bending in an inhomogeneous medium produced an uneven ray distribution and thus uneven resolution (Fig. 8).This does not necessarily mean a lack of constraints on regions within the array without ray passing, since the initial velocity model was 1-D and the rays were evenly distributed before any iteration.With increasing number of iterations, the rays gradually bend toward the high velocity zones with increasing ray densities therein.As a result, the high-velocity region results are modified by more iterations and the resolutions increase because ray crossings are proportional to ray densities for evenly distributed stations.
The EW trending fast anomaly might be caused by metamorphic rocks existing in the Central Mountain Range.It is interesting to note that the east side of the fast anomaly exhibited a vertical downward feature suggesting overall cold material in this part of the mountain.The mountain range tightly abutting sedimentary layers results in a northward stark velocity contrast a 0.5-km depth (near Suao in slice 2 of Fig. 9).Its potential association with Sanshin Fault activity and the in situ high seismic activity is worth further exploration.The tomography results show that the Chingshui location, at least the shallow part, is a relative low velocity zone.However, we need to combine these results with a P-wave velocity model to determine the cause, i.e., high temperature and/or the presence of fluids.In addition, we need to combine these data with body wave tomography to determine the depth extension.
Based on a previous geomagnetic survey, Yu and Tsai (1979) proposed an intrusive igneous dyke with high magnetic susceptibility to explain the high EW trending magnetic anomaly located between Ilan and Luodong (Fig. 10), or just north of our fast anomaly zone (Tong et al. 2008).Since our tomography results only have resolutions up to a depth of 1.5 km and the igneous dyke may not intrude to the surface, the sedimentary basement data remain the most prominent features in our study.Tong et al. (2008) conducted a magnetotelluric survey in Ilan and presented data indicating a geothermal reservoir that surfaced at Chingshui.To the north, however, the reservoir was covered by low resistivity cap rock.This cap rock may be as deep as 1 km and corresponds well to the EW trending fast anomalies in our results (the cap rock in Fig. 10).

coNcluSIoNS
A high-resolution 3-D S-wave velocity model of the South Ilan Plain was inverted using the ambient noise tomography method using local dense array data.The fast marching method adoption to deal with ray bending in an inhomogeneous medium could map the travel time anomalies more accurately to the true locations of heterogeneous structures, e.g., high velocity zones in particular.According to checker board test results, resolutions could be reached as high as 1.0 km at a depth of 1.0 km and as high as 2.0 km at a depth of 1.5 km.The resulting patterns were predominantly constrained by the sedimentary basement thickness.In addition, nearly EW trending fast anomalies were observed beneath the mountainous terrain abutting to the South Ilan Plain, which may be attributed to existing metamorphic rock and may correspond to the cap rock of a proposed thermal reservoir.The Chingshui location consistently exhibited a low S-wave velocity zone to a depth of 1.5 km, corresponding well with its history of thermal explorations.

Fig. 1 .
Fig. 1.Seismicity (dots), faults (dashed lines), relative plate motions (red arrows), and sedimentary basement contours (color lines) of the Ilan Plain and its vicinity.The inset shows the Ilan Plain location.The white pentagons are cities for reference.Focal mechanisms of the 2005 doublet earthquakes are shown and yellow triangles indicate the distributions of local dense array in this study.(Color online only)

Fig. 3 .
Fig. 3. Stacking EGFs for station pairs with similar interstation distances and arranged by interstation distances.The final EGF for each station pair for stacking is in term stacking of all qualified Julian day cross correlation functions.The stacking CCFs are filtered at the period band 0.5 -3.7 s. (Color online only)

Fig. 4 .
Fig. 4. Mean residuals and associated standard deviations as a function of iteration numbers.(Color online only)