The Deep Electrical Structure of Southern Taiwan and Its Tectonic Implications

The Taiwan orogen has formed as a result of the arc-continent collision between the Eurasian continental margin and the Luzon volcanic arc over the last 5 million years and is the type example of an arc-continent collision. The tectonic processes at work beneath Taiwan are still debated; the available data have been interpreted with both thin-skinned and lithospheric collision models. In 2004, the Taiwan Integrated Geodynamical Research (TAIGER) project began a systematic investigation of the crustal and upper mantle structure beneath Taiwan. TAIGER magnetotelluric (MT) data from central Taiwan favor a thick-skinned model for that region. The Taiwan orogen becomes younger to the south, so the earlier stages of collision were investigated with a 100-km-long MT profile in southern Taiwan at latitude of 23.3°N. Data were recorded at 15 MT sites and tensor decomposition and two-dimensional inversion were applied to the MT data. The shallow electrical resistivity structure is in good agreement with surface geology. The deeper structure shows a major conductor in the mid-crust that can be explained by fluid content of 0.4 1.4%. A similar feature was observed in central Taiwan, but with a higher fluid content. The conductor in southern Taiwan extends to lower crustal depths and is likely caused by fluids generated by metamorphic reactions in a thickened crust. Together the central and southern Taiwan MT profiles show a crustal root beneath the Central Range.


InTroDucTIon
The island of Taiwan was formed by the arc-continent collision between the Luzon volcanic arc on the western margin of the Philippine Sea Plate (PSP) and the passive continental margin of southeastern China, beginning some 5 Ma ago (Suppe 1981;Ho 1986).The polarity of subduction changes across the island; in the northeast the PSP is subducting to the northwest beneath the Eurasian Plate (EUP), while in the south a portion of the continental EUP subducts to the east beneath the PSP (Tsai 1986;Seno 1993;Kao 2000).Convergence in southern Taiwan occurs at 80 mm yr -1 in the direction N54°W (Yu et al. 1999).Taiwan represents a young active orogen and is a unique natural laboratory for studying arc-continent collisions.
The major geological structures in Taiwan are aligned from N20°E to N30°E and the major tectonic provinces are separated by narrow fault zones and geologic contacts (Fig. 1).The major geological units in Taiwan are (from west to east): (1) the Coastal Plain (CP) comprised of Quaternary sedimentary rocks; (2) the Western Foothills (WF) which are a classic fold and thrust belt; (3) the Hsuehshan Range (HR) consisting of shallow marine quartizite and slate; (4) the Backbone Range (BR) of metapelite and metagraywackes; (5) the Eastern Central Range (ECR) which includes the pre-Tertiary Tanano schist, phyllite, metasandstone, and metabasite; (6) the Longitudinal Valley (LV) which overlies the suture zone between the EUP and PSP; (7) the Coastal Range (CR) consisting of deformed forearc and remnant volcanic arc materials (Fig. 1).The WF and HR are separated by the Chuchih fault, which merges with the Lishan fault at the southern extent of the HR to form the Laonungchi or Chaochou fault.The HR and BR are mainly uplifted from a thick sequence of Tertiary marine argillaceous sediments.These sediments were altered to slate and phyllite after a major episode of Plio-Pleistocene deformation associated with the arc-continent collision.In general, metamorphic grade increases across Taiwan from west to east (Ho 1999).
A range of tectonic models have been developed for the Taiwan orogen (Suppe 1981;Hsu et al. 1995;Wu et al. 1997;Lin 2002).These models encompass a broad spectrum from the thin-skinned (Suppe 1981) tectonic model to the lithospheric collision model (Wu et al. 1997).The thinskinned model features a shallow eastward dipping décollement that was inferred from geological data, shallow (1 -15 km) exploration in western Taiwan (Wang et al. 2002;Ma et al. 2006;Chiang et al. 2008).In this model, all deformation occurs above the décollement.In contrast, the lithospheric collision model was based on geophysical data (Rau and Wu 1995;Wu et al. 1997Wu et al. , 2007) ) that gave evidence for a crustal root beneath the Central Ranges.
Existing seismic reflection data and geological mapping confirm that thin-skinned thrusting occurs in the western third of Taiwan (Suppe et al. 1981;Wang et al. 2002).However limited geophysical data has been collected in the Central Ranges and studies in these areas are needed to provide further constraints on deep structure and distinguish between the various tectonic models.The TAIGER (Taiwan Integrated Geodynamics Research) project has been investigating the tectonics of Taiwan since 2004 with the goal of addressing these issues.The TAIGER project (Wu and US/Taiwan TAIGER teams 2007;Okaya et al. 2009) aims to provide a detailed study of the Taiwan orogen through an integrated program of petrophysics, geodynamic modeling, seismology and magnetotellurics (MT).The initial results of the TAIGER MT studies in central Taiwan were described by Bertrand et al. (2009) and the profile location is shown in Fig. 2.This study showed that the deformation at that location was dominated by thick-skinned tectonics.The Taiwan orogen becomes younger to the south (Teng 1990;Sibuet and Hsu 2004) and provides an opportunity to study the temporal evolution.This paper describes MT data that were collected in southern Taiwan with the goal of studying an earlier stage in the temporal evolution of the Taiwan orogen.

MAgnEToTEllurIc (MT) DATA collEcTIon AnD AnAlySIS
Two types of instruments were used to record MT data on the southern Taiwan transect.Long-period (10 -10000 s) NIMS recording systems designed by Narod Geophysics Ltd. (Canada) were used for investigating the deep crustal structures at all 15 stations.In addition, broad-band (0.001  -1000 s) Phoenix V5-2000A systems were used for enhancing the shallow resolution at 7 out of 15 stations.The broadband MT systems were not deployed at all sites due to the extreme topography in the study region.
The MT time series data acquired at each site were processed using the statistically robust algorithms from Egbert (1997) for the long-period MT data and Jones et al. (1989) for the broadband MT data.A remote-reference station on Penghu Island was used to remove bias for both types of MT data (Gamble et al. 1979).For sites with coincident longperiod and broad-band data, the MT transfer function estimates were merged to form a continuous sounding between periods of 0.001 and 1000 s.The apparent resistivity and phase curves for all stations are shown in Fig. 3.
Tensor decomposition was used to investigate the dimensionality of the MT data.An overall 2D geoelectric strike direction of N29°E was determined using the multistation decomposition from McNeice and Jones (2001).This strike direction is essentially parallel to the geological strike.Single site tensor decomposition also produced local strike directions that were nearly parallel to the major geological boundaries for the period range from 13 -7500 s (Fig. 4b).The agreement between the geoelectric strike and the strike of the geological provinces suggest that a 2D interpretation may be possible.Values of the shear and twist angles, and the corresponding misfit of the tensor decomposition model are shown in Fig. 4a.Large shear and twist angles indicate strong three-dimensional distortion occurs at sites located in the Backbone ranges .
Induction vectors graphically represent the transfer function between the horizontal and vertical magnetic field components, and are widely used for studying the dimensionality of MT data (Santos et al. 2001;Bedrosian et al. 2004;Brasse et al. 2009).The convention of Parkinson (1962) is used throughout this paper in which the in-phase (real) induction vectors point toward conductors.At a period of 10 s the MT signals sample upper crustal structure and the induction vectors do not point in a uniform direction (Fig. 5), and indicate presence of some shallow conductors.At a period of 1000 s the MT signals are sampling the lower crust, most induction vectors point toward to southeast direction reflecting the presence of the Pacific Ocean.If the subsurface resistivity structure is 2D then the induction vectors should be perpendicular to the geoelectric strike direction.This can be seen to be the case at a period of 100 s in Fig. 5 at the east and west ends of the profile.Over the BR and ECR the induction vectors are directed to the south, which could be related to an along-strike conductive anomaly in the crust.Alternatively, the south pointing induction vectors could be related to the coast-effect, where the conductive anomaly is the ocean.These indicators of 3D geoelectric structure suggest that a 2D analysis must be applied with caution.In this situation, features observed in a model derived with a 2D data analysis should be rigorously tested.The analysis must ensure that features in the resistivity model are required by the MT data and free of artifacts from the coast-effect.Details of this analysis are described in section 4.

Two-DIMEnSIonAl (2D) rEgulArIzED In-vErSIon
The regularized 2D MT inversion algorithm of Rodi and Mackie (2001) was applied to convert the MT data into a resistivity model.The apparent resistivity, phase and tipper data were rotated into a co-ordinate system of N29°E determined from tensor decomposition.Inversions used a starting model which included the seawater (0.3 Ωm) in the Taiwan Strait and Pacific Ocean embedded in a 100 Ωm half-space.Many variations of inversion parameters and starting models were tested and the preferred inversion models were generated by inversion of all impedance and the real component of the vertical magnetic field (tipper) data.Error floors of 20%, 4.35° and 0.1 were used for the apparent resistivity, phase and tipper data, respectively.These error floor values were chosen to fit strong distortion in southeastern profile as shown in Fig. 4a and were motivated by the approach used by Booker et al. (2004) and Becken et al. (2008) to down-weight the influence of 3D induction and static-shift effects.Lower error floors were used but resulted in rougher models.Static shift coefficients (Fig. 6) were estimated using two approaches: (a) automatic estimation by the inversion algorithm and (b) down-weighting the apparent resistivity data using commercial 2D inversion software of WinGlink ® .Both these approaches gave similar resistivity models.
Multiple inversions with various values of the trade-off parameter described by Booker et al. (2004) were performed to generate a trade-off curve shown in Fig. 7.The trade-off curve has two branches: (1) an initial steep portion (τ = 1000 -30) in which the NRMS data misfit drops rapidly with only a moderate increase in model roughness and (2) a second portion (τ = 30 -0.1) in which model roughness increases rapidly with minimal decrease in NRMS data misfit.The initial steep portion of the trade-off curve represents models where the smoothing constraint is emphasized; the second portion of the trade-off curve represents models where the data fit is of primary importance.The preferred model represents a compromise between fitting the data (low NRMS data misfit) and generating a smooth (realistic) resistivity model, and was obtained with τ = 30 which gave an NRMS value of 1.6.An ideal value of the NRMS would be 1, and the range of 1 -2 is statistically acceptable.As shown in Fig. 3, the calculated responses from the best-fitting model (Fig. 8) are in good agreement with the measured MT data, except for the long-period TM phases at the westernmost 7 sites where poor fitting indicates strong 3D effects that cannot be accounted for with a 2D response.At these stations, the TM phase data is out of quadrant at period of 100     -1000 s and cannot be fit with a 2D inversion.It may be the result of the galvanic distortion effects or the coast effect as described in the section on 3D modeling.

ThrEE-DIMEnSIonAl ForwArD MoDElIng AnD ThE coAST EFFEcT
To ensure that conductive features in the MT model (Fig. 9c) are not contaminated by the coast-effect, data analysis was performed to determine the sensitivity of the MT data collected in southern Taiwan to the presence of the surrounding ocean.A 3D resistivity model for Taiwan was generated (Fig. 10a) which included the regional bathymetry (0.3 Ωm) and the conductors C1 -C4 imaged by the 2D inversion (Fig. 9c).The horizontal mesh size was 10 km in both east-west and north-south directions.The vertical mesh used 57 vertical layers with 10 air layers.The 3D model had a horizontal extent of 1100 km and a vertical extent of 410 km.The background resistivity is assumed to be 1000 Ωm.The near-surface conductors C1, C2 and C3 were represented by blocks with resistivity of 7.3 Ωm extending from the surface to a depth of 5 km.The mid-crustal conductor C4 was modeled as a 34 Ωm block to a depth of 30 km.Also, the NNE-SSW direction for C1 -C4 is assumed on the basis of geological constraints.These properties were chosen to best represent the conductors imaged in the inversion model in Fig. 9c.The boundaries of these starting features are denoted by the dashed rectangular blocks in Figs.10b and c.Synthetic MT data were then computed with a commercial 3D modeling software of WinGlink ® (Mackie et al. 1993) for the period range 0.001 to 10000 s at the locations where data was collected in southern Taiwan.These synthetic 3D MT data were then inverted using the same 2D inversion that produced the model in Fig. 10b.The 2D synthetic inversion model (Fig. 10b) was then compared with the 3D model used to generate the synthetic MT data and the percent difference is displayed in Fig. 10c.The plot shows that the smallest difference between the synthetic and actual models occurs in the regions where the conductors C1 -C4 are located.This suggests that the conductive features C1 -C4 are well resolved, and not contaminated by the coast-effect from the surrounding ocean.The difference becomes significant below a depth of 50 km.The synthetic 3D data are illustrated in Fig. 10d.Note that the presence of the 3D ocean does not produce the out-of-quadrant phases at the western end of the profile.These phases are likely due to distortion by a 3D conductor onshore that is somehow coupled to the ocean (Lezaeta and Haak 2003;Ichihara and Mogi 2009).
Figure 11 shows other quantities that can be used to investigate the dimensionality of MT data.These include the Bahr skew η defined by Bahr (1991) and the phase tensor term β defined by Caldwell et al. (2004).Together these quantities indicate that the MT data along our south MT pro-file are not significantly affected by the coast effect within the period range 10 -1000 s.The MT data are influenced by surrounding ocean at long periods (1000 -10000 s), but with different scales from these two decomposition analysis.All the parameters indicate that the coast effect is less on the southern profile than elsewhere in Taiwan.These results in addition to the well defined geoelectric strike described in section 2 support a 2D interpretation of the MT data collected in southern Taiwan up to a period of 1000 s.

gEoElEcTrIc STrucTurE AnD ITS InTEr-prETATIon
Figure 9c shows the preferred inversion model (τ = 30) overlaid with selected earthquake hypocenters (M W ≥ 3.0) from Central Weather Bureau, CWB located within 10 km of the MT profile.At shallow depths (0 -10 km), the resistivity progressively increases from west to east across Taiwan, which is as expected for the change in lithology from sedimentary rocks (CP), to meta-sedimentary rocks (WF and HR), to higher-grade metamorphic rocks (BR and ECR).Three significant conductivity anomalies are observed in the shallow crust (C1, C2 and C3) and a fourth is located in the mid-crust (C4).
C1 is spatially associated with the Quaternary sedimentary rocks of the Coastal Plain in southwestern Taiwan and coincides with a zone of low seismic velocity (Kim et al. 2005;Wu et al. 2007).High porosity sedimentary rocks, saturated with interconnected fluids can account for this conductive C1 anomaly.Higher resistivities are observed beneath the Western Foothills and indicate the lower porosity of Neogene sedimentary rocks with low-grade metamorphism.The geothermal gradient in the Central Range is high (Wang et al. 1994;Lin 2002) and the region has numerous geothermal phenomena (Ho 1999).Relatively high heat flow of 350 mW m -2 is observed close to C2 (Fig. 9b).Thus, C2 may be related to the circulation of geothermal fluids in this area where the fractures associated with the Central Range Fault (CERF) between the BR and ECR could give the necessary permeability.Further east, feature C3 is located within the Coastal Range and may be related to a typical marine mélange with heterogeneous rock materials consisting of a sheared argillaceous matrix mixed thoroughly with native and exotic tectonic fragments (Hsu 1956(Hsu , 1976;;Ho 1999).Feature C3 is also coincident with the low V p and high V p /V s zone reported by Wu et al. (2007).This phenomenon indicates that the conductive anomaly may be related to a fractured formation and fluid-filled seismicity zone (Chen andChen 2000, 2002;Chen et al. 2007).
The most important feature shown in the inversion model is the conductive anomaly C4 located beneath the Central Ranges in the mid-crust with the top at a depth of 15 km (Fig. 9c).The spatial smoothing imposed during MT inversion can smear a conductor to depth, mak-   ing it difficult to interpret the depth extent of a conductive feature.To solve this problem, constrained inversions (Li et al. 2003) were used to verify that the C4 conductor is required by the MT data within the mid-crust.In this approach, the model is forced to be resistive below a given depth.If the MT data cannot sense the bottom of the conductive layer, then the resistive region will not affect the data misfit.The depth of the top of the resistive half-space is moved upwards until the MT data can no longer be fit.This determines the shallowest conductor of 15 km permitted by the data.This conductor is coincident with a zone of low shear-wave velocity observed beneath the Central Ranges (Rau and Wu 1995;Kim et al. 2005;Wu et al. 2007).The resistivity of C4 is an order of magnitude lower than the surrounding crust at the same range of depths.The most likely explanation of the low resistivity might be attributed to saline aqueous fluids.The electrical resistivity of brine depends on the salinity, pressure and temperature (Sourirajan and Kennedy 1962;Quist and Marshall 1969).
The resistivity of C4 is 30 Ωm and this can be explained with fluid content of 0.4 -1.4% if the fluid resistivity is in the range 0.01 -0.05 Ωm (Nesbitt 1993) and the cementation factor value of m = 1.5 in Archie's law (Archie 1942).The existence of deep crustal fluids in southern Taiwan is significant because fluids may play an important role in crustal deformation in active orogens (e.g., Chen et al. 2002;Soyer and Unsworth 2006).The elevated fluid content and high geothermal gradient shown in Fig. 9b, have been related to uplift and thickening of the crust in the Central Range (Lin 2002;Wang and Hung 2002).
Two widely discussed end-member models of Taiwan tectonics are the thin-skinned and lithospheric collision models.As observed in central Taiwan (Bertrand et al. 2009), C4 is located within the mid-crust and is inconsistent with a thin-skinned model in southern Taiwan.The thin-skinned model assumes that all the collisional deformation occurs above the inferred upper crustal decollement (Carena et al. 2002).However, C4 indicates that significant deformation occurs below the decollement.In the lithospheric collision model (Fig. 9d) the presence of a ductile mountain-root (a thickened crust) beneath the Central Range suggests ductile deformation occurs contiguously within the upper crust to upper mantle beneath the Central Range (Wu et al. 1997).Deformation within a developing crustal root would facilitate pro-grade metamorphic reactions that would produce the aqueous fluids imaged as the C4 conductor (Wannamaker et al. 2002).Thus, in the thickened crust, fluids together and possible thermal effects (Wang et al. 1994;Lin 2002) related to the process of ductile deformation could explain the presence of a low-resistivity zone in the mid crust.
The resistivity structure observed in southern Taiwan has some similarity with that reported in central Taiwan by Bertrand et al. (2009).In both areas a mid-crustal conductor was observed beneath the western Central Ranges.However the conductor beneath central Taiwan appears to be stronger than that observed in southern Taiwan.As stated above, the fluid content appears to be higher in central Taiwan (1 -2%), compared to southern Taiwan (0.4 -1.4%) assuming that salinity is the same in each location.The northward increase in fluid content can be understood by considering that the age of the collision zone increases northward.The Luzon arc to the south of Taiwan is a subduction zone where the Eurasian Plate subducts with minimal internal deformation at shallow depths.No MT data is available from the Luzon Volcanic Arc south of Taiwan, but in comparable subduction zones elsewhere, such as the well studied Cascadia subduction zone, a dipping conductor is observed on top of the subducting slab.Once the slab reaches a depth around 100 km, slab dewatering triggers mantle melting and forms the volcanic arc (Wannamaker et al. 1989;Soyer and Unsworth 2006).
The resistivity models observed in Taiwan are quite different to those reported for other subduction zones.As the continent-continent collision matures to the north, the thickness of the continental crust in the orogen increases and produces increasing volumes of metamorphically derived fluids (Wannamaker et al. 2002;Bertrand et al. 2009).
It also appears that the conductor is narrower in central Taiwan than that in southern Taiwan.This may be due to the localization of fluid flow associated with shear zones as the orogen develops.Both porosity analysis and spatial resistivity beneath Central Range differences between two MT profiles indicate that the stages of the Taiwan orogen gradually change from active-(central Taiwan) to post-collision (southern Taiwan).

concluSIonS
This paper describes the analysis and interpretation of joint long-period and broad-band MT data collected as part of the TAIGER project in southern Taiwan.These MT data have been carefully analyzed and exhibit an overall 2D nature.The 3D forward modeling that includes the seawater suggests that the 2D inversions are valid at periods less than 1000 seconds.These conductive anomalies observed in the MT model are interpreted to be associated with interconnected fluids within the crust.The shallow electrical structures correlate well with geological units mapped at the surface.The resistivity gradually changes from lower to higher resistivities from west to east as the lithology changes from sedimentary to metamorphic rocks.A mid-crustal conductive anomaly (C4) is related to the existence of interconnected fluids beneath the Central Range, coupled with a thermal effect.The existence of this conductive anomaly (C4) within the mid-crust reflects the crustal root in southern Taiwan.

Fig. 1 .
Fig. 1.Tectonic map of the Taiwan region showing major geological province boundaries, geological provinces are: CP: Coastal Plain; WF: Western Foothills; HR: Hsuehshan Range; BR: Backbone Range; ECR: Eastern Central Range; CR: Coastal Range; LV: Longitudinal Valley.Location of geothermal measurements are shown in white circles and red dots used in the profile.MT study area indicated by white solid line.EUP: Eurasian Plate; LA: Luzon Volcanic Arc; MT: Manila Trench; OT: Okinawa Trough; PSP: Philippine Sea Plate.

Fig. 2 .
Fig. 2. Study area in southern Taiwan showing MT station locations with dotted symbols, the red dots show NIMS stations, black dots show combined NIMS and V5 stations.The profile of Bertrand et al. (2009) in central Taiwan is shown with inverted triangles.Geological provinces are: CR: Coastal Plain; WF: Western Foothills; HR: Hsuehshan Range; BR: Backbone Range; ECR: Eastern Central Range; CR: Coastal Range; LV: Longitudinal Valley.

Fig. 3 .
Fig. 3.All apparent resistivity, phase and Hz curves for all sites are shown rotated to a geoelectric strike direction.The calculated responses belong to the model in Fig. 9c.

Fig
Fig. 3. (Continued) Fig. 6.Static shift for 2D inversion along the profile.The default static shift coefficient is set 1, static shift values small than 1 indicate that the curve needs to be shifted down and vice versa for values bigger than 1.

Fig. 7 .
Fig. 7. Variation of normalized RMS misfit and model roughness as a function of the trade-off parameter, τ.The vertical-to-horizontal smoothness parameter, α, was set at 3.

Fig. 8 .
Fig.8.Pseudo-sections of apparent resistivity, phase and tipper data.The two MT modes are distinguished by the direction in which the electric current flows.This direction is parallel to strike in the transverse electric mode (TE) and perpendicular to strike in the transverse magnetic (TM) mode.Data at some sites and periods (white regions) were excluded from the inversion analysis due to severe distortion and noise.Response of the inversion model in Fig.9cis also shown.
Fig. 10.3D forward modeling of the coast-effect and its influence on features observed in the 2D inversion (a) Map view of the 3D forward model with mesh size of 10 × 10 km.(b) 2D inversion result of synthetic MT data computed from the 3D model.(c) Percent difference between the actual model Fig. 9c and the synthetic inversion result in (b).(d) Typical responses of apparent resistivity, phase curves from forward modeling are shown.

Fig. 11 .
Fig. 11.Investigation of dimensionality of the MT data in Taiwan.The parameter η is the skew suggested by Bahr (1991), and β was suggested by Caldwell et al. (2004).If these quantities are small, the MT data can be considered 1D or 2D (white areas).Larger values indicate 3D behavior associated with the surrounding ocean.