An Investigation of the 3 D Electrical Resistivity Structure in the Chingshui Geothermal Area , NE Taiwan

The Chingshui geothermal area southwest of the Ilan plain is identified as a western extension of the Okinawa Trough in the northern Taiwan subduction system. Numerous geophysical, geological and geochemical investigations have been conducted since the 1970s by the Industrial Technology Research Institute, the Chinese Petroleum Corporation of Taiwan and the National Science Council of Taiwan. These studies indicated that the Chingshui stream is one of the largest geothermal areas for electricity generation in Taiwan. However, the power generation efficiency has not met initial expectations. Magnetotelluric (MT) data analyses show that the Chingshui geothermal region is a geologically complex area. A full three-dimensional (3D) inversion was therefore applied to reprocess the MT data and provide the detailed electrical structure beneath the Chingshui geothermal region. The 3D geoelectrical model displays an improved image that clearly delineates the Chingshui geothermal system geometry. Two conductive anomalies are imaged that possibly indicate high potential areas for geothermal energy in the Chingshui geothermal system. One of the potential areas is located in the eastern part of the Chingshui Fault at shallow depths. A significant conductive anomaly is associated with high heat flow and fluid content situations southwest of the geothermal manifest area at depth. A higher interconnected fluid indicates that this area contains the highest potential for geothermal energy in the Chingshui geothermal system.


INTRODUCTION
The Chingshui geothermal area in northeastern Taiwan is potentially one of the largest high-temperature areas in the country (Fig. 1).The integration of geophysical and geological evidences has confirmed the area to be a hot-waterdominated geothermal system (Su 1978;Tseng 1978;Hsiao and Chiang 1979;Yu and Tsai 1979).The National Science Council of Taiwan established the first geothermal power plant in 1981 to generate an expected maximum 3 MW of electricity.However, the average output was only 1.18 MW in the first year of operation and the capacity declined by more than 50% to 0.52 MW after 2 years.The plant was decommissioned because of low productivity.By 1992 the average electricity output had decreased to 0.18 MW (Fan et al. 2005).
Although the Industrial Technology Research Institute and Chinese Petroleum Corporation of Taiwan had explored the area during the 1970s.The exact location of the geothermal reservoir still remains uncertain.The geophysical investigation was therefore restarted to promote renewable energy in 2006 (Tong et al. 2008).This investigation combines magnetotelluric (MT), magnetic, gravity, and borehole data to reconstruct a conceptual model to provide more details of the complex Chingshui geothermal system.Chang et al. (2014) thought that the surface hot spring locations are inconsistent with the low resistivity region locations proposed in MT inverted images by Tong et al. (2008).They reprocessed the MT data set again and integrated it with direct-current resistivity data and well logs to rebuild a new three-dimensional (3D) visualization model.
Previous MT studies applied one-and/or two-dimensional inversion approaches (Tong et al. 2008;Chang et al. 2014) without dimensionality analysis which may be unable to delineate the geothermal structures if the MT data is mixed with 3D effects, galvanic distortions (Bertrand et al. 2009(Bertrand et al. , 2012;;Chiang et al. 2010Chiang et al. , 2011) ) or topographic effects (Wannamaker et al. 1986;Jiracek 1990).We attempted to carefully reprocess these MT data step by step to evaluate whether the MT data can be appropriately inverted using 1D or 2D inversion.A 3D inversion approach was then applied in this study.

GEOLOGICAL SETTING
The Chingshui geothermal area southwest of the Ilan plain is identified as the western extension of the Okinawa Trough in the northern Taiwan subduction system (Suppe 1984;Yeh 1989;Liu 1995;Lin et al. 2004).The Chingshui Fault traverses this region from the northwest to southeast along the Chingshui stream (Tseng 1978;Hsiao and Chiang 1979), and a northeast-southwest trend anticlinorium axis parallel to the stratum is present in this region, as shown in Fig. 1.The main lithological unit is a fractured Miocene Lushan Slate Formation in the northeast tertiary submetamorphic zone (Fig. 1).The geological strike trends northeast to southwest and dips to the southeast at angles ranging from approximately 35 -90° (Hsiao and Chiang 1979).Lithological analysis shows Jentse, Chingshuihu, and Kulu Members within the Lushan Formation from the surface to deep.The Jentse Member is composed primarily of metasandstone such as albite, chlorite, illite, and quartz (Yui et al. 1993) intercalated with slate, whereas the underlying Chingshuihu and Kulu Member areas are primarily composed of slate (Hsiao and Chiang 1979).
Numerous hot spring outcrops are presented within the slates due to lower porosity and permeability within this unit.Geothermal fluid has difficulty transferring into these outcrops under these conditions.The fractured slate and metasandstone of the Jentse Member (Hsiao and Chiang 1979), joints (Tseng 1978), and other extensive fractures influence the geothermal fluid paths according to drilling wells and geological structures.Chemical fluid measurements also implied that the geothermal fluid may originate from rain or river water flowing into the fracture system and then returning to the surface as indicated by low Na + and Cl - ion concentration measurements in contrast to ancient brines (Tseng 1978).This is insufficiently supported by source water originated from greater depths, because the geochemical samples are contaminated with surface fluids that could potentially dilute the ancient brine concentration.

DATA COLLECTION AND ANALYSIS
The MT time series data were measured over 72 h at all stations to improve the data quality.This measurement period would expect data to resolve in periods of 1000 s or at the least 100 s.The MT time series data acquired at each site were pro cessed using statistically robust algorithms from Jones et al. (1989) for the MT data.Two electric field orthogonal components and three magnetic field components were recorded simultaneously using three Phoenix MTU-5A and two MTU-2E systems.The remote reference technique after Gamble et al. (1979) was used to improve the data quality.Impedance data for periods > 10 s were removed due to lower data quality as identified by larger error bars and/or scattered data points, as previous papers (Tong et al. 2008;Chang et al. 2014) did not show the detailed MT responses.We therefore could not compare the difference between old and our reprocessed data.However, previous MT electrical resistivity models (Tong et al. 2008;Chang et al. 2014) (Caldwell et al. 2004), to evaluate whether a traditional 2D inversion is suitable for the MT data.Finally, a 3D electrical model using WSINV3D-MT code (Siripunvaraporn et al. 2005a) was reconstructed to fully represent the Chingshui geothermal structure.

DIMENSIONALITY AND DIRECTIONALITY PARAMETERS
Dimensionality and directionality is the first step of MT data analysis to evaluate whether to apply 2D or 3D modeling (Bahr 1988;Groom and Bailey 1989;Jones and Groom 1993;Smith 1995;Caldwell et al. 2004).These parameters were calculated using a frequency-dependent 2 × 2 complex impedance tensor for two horizontal orthogonal electromagnetic waves.Parameters were derived using the McNeice-Jones scheme (McNeice and Jones 2001) and phase tensor decomposition (Caldwell et al. 2004) for a range of periods at each site, respectively.

McNeice-Jones Scheme
Data were set with an error floor of 3.50% for apparent resistivity, and 1° for the phase.Three-decade wide period bands (0.003 -10 s) were selected for decomposition at all sites.For longer periods the regional strikes in northeastern Taiwan could be approximated by 45° (Bertrand et al. 2012).Because of the inherent geoelectrical strike analysis ambiguity, the strikes were aligned with the geological strike as a reference for periods < 10 s (Fig. 2).The geoelectric strikes appear out of sequence in high-to low-period bands incon-sistent with the dominating regional strike (Fig. 2).These unstable geoelectrical strike orientations demonstrate the complex structure and galvanic distortion probably caused by streams in the region.Additionally, the average shear and twist value histograms (Figs.3a, b) seem higher than those for the Western Foothills in Central Taiwan (Chiang et al. 2008).The shear values were primarily within the ±40° range and the twist values within the ±30° range over whole periods.The marginally smaller twist values indicate significantly distorted sites in the region.Bahr skew values (Fig. 3c) are the potentially significant 3D indicator (Bahr 1988;Simpson and Bahr 2005) that show strong lateral heterogeneities, despite most RMS misfits (Fig. 3d) being a good fit and consistent with the Groom-Baily model.The dimensionality parameters and unstable geoelectric strikes show strong 3D distortion in the region.

Phase Tensor Decomposition
The regional phase information is recovered directly from the observed impedance tensor without assuming regional conductivity structure dimensionality (Caldwell et al. 2004), which is another method to assess the Chingshui MT response dimensionality.Figure 4 shows the ellipse major axis and phase tensor skew angle β.The major axis alleged geoelectric strike can be derived from α -β according to Caldwell et al. (2004).In the ideal 2D case β must be close to zero with a constant major axis direction at a range of periods in the 2D conductivity structure, whereas β is non-zero in the general 3D situation.The phase information (Fig. 4) shows that the majority of β values appear close to zero and the ellipses appear symmetrical at short-period bands (T < 0.1 s).The β angles increase significantly away from zero at longer periods (T > 0.1 s) and the major ellipses axes are not as aligned compared to those at short periods.The misalignment reflects an asymmetrical 3D response, which implies that 2D MT data analysis may produce a poorer result.The phase tensor parameters gradually show a strong 3D effect between sites 20 and 23 (W -E trend) below T = 10 s, which implies that a conductive anomaly may be located at depth in the area.Thus, the 3D inversion process can potentially delineate an improved electrical model for the Chingshui geothermal area compared with previous studies (Tong et al. 2008;Chang et al. 2014).

3D INVERSION
Several 3D inversion algorithms have been developed since 1993 (Mackie and Madden 1993;Newman and Alumbaugh 2000;Sasaki 2004).Most of these algorithms require powerful workstations or parallel computers, which make field application difficult (Siripunvaraporn et al. 2005a, b;Tuncer et al. 2006;Xiao et al. 2010).The 3D inversion practical application has been made more feasible on desktop computers using the data-space approach by Siripunvaraporn et al. (2005a) who significantly decreased the computational costs.This study therefore employs the 3D inversion algorithm to reanalyze the MT data set and reinvestigate the Chingshui region geothermal reservior details.
The initial model is constructed with 44(x) × 42(y) × 50(z) = 92400 cells and 7 air layers.Each horizontal mesh space is set 10 km outside the target area, with the variable spaces set approximately 50 -100 m.The vertical mesh spaces begin from 100 m for shallow seawater depths, before gradually increasing the spaces downward to the floor.The background resistivity is set as a homogeneous halfspace of 100 ohm-m including ba thymetry.The seawater resistivity is set to 0.3 ohm-m and locked in the inversion.Eight impedance elements (real and imaginary parts of Z xx , Z xy , Z yy , and Z yx ) are calculated on a workstation.The RMS rapidly decreases at the first iteration and then gradually converges to 1.4 at the fourth iteration.The RMS is not converged further as a flat curve after the fourth iteration.This study identifies the fourth iteration turning point to be a best-interpretation model.Most of the sites appear to fit the phase data, whereas the apparent resistivities (Fig. 5) show a possible static shift and/or topographic effects at sites 4, 5, 16, 17, 18, and 27.The model could be compensated by these static shift data at shallow depth to incorporate nearsurface structure (Sasaki 2004), which indicates that the inverted model is still fine at deep depth in the inversion (e.g., Sasaki and Meju 2006).The 3D inversion code can only invert impedances and it is thus impossible to assign different error floors to apparent resistivity and phase as in typical 2D inversions.In general 3D cases diagonal components can be as large as the off-diagonal ones.Therefore, we manually set data errors as 30% of Z Z xx yy

# and 20% of Z Z xy yx
# for respective impedances in the inversion.

3D MODEL OF ELECTRICAL RESISTIVITY
Unstable tensor parameters (e.g., strikes and phase tensor skew angles) show that the Chingshui geothermal region is a geologically complex area.This study performed a full 3D inversion using the WSINV3DMT code (Siripunvaraporn et al. 2005a) to reanalyze the MT data acquired in the Chingshui geothermal area.Figure 6 shows horizontal slices of the final electrical resistivity structure.A significant electrical resistivity boundary appears along the synclinorium axis in the south (Fig. 6) that implies an important electrical resistivity discontinuity probably related to the geothermal reservoir along the structure.The regions far away from MT stations can be identified as compensation effect at near-surface.Two conductivity anomalies appear in range of depth 0.1 -0.5 km (C1 in Figs.6a -c) and depth 1.1 -1.7 km (C2 in Figs.6f -i), respectively.The shallow conductive anomaly C1 and southeastward dipping of the Lushan Formation (Hsiao and Chiang 1979) could be considered as fluid sources supporting the conductive anomaly C2.While C2 is located in a high heat flow region at the northern synclinorium axis (discontinuity of electrical resistivity), ideal geothermal system circulation is implied.
Two profiles near the MT stations are shown in Fig. 7 to represent the electrical structures of the Chingshui geothermal system.The northern conductive anomaly C1 ( A -A' profile in Fig. 7a) in horizontal distance from 0 -1.3 (Chingshuihu Member, CM) can be correlated with alluvial deposits consisting of hydrous mud, sand, and gravel, which gradually vanish in the headstream of the Chingshui stream as the slates or metasandstones appear to the south.A ~50 ohm-m of conductive anomaly (2.2 -3.1 km in horizontal distance) is located in a geothermal well and hot spring outcrop area (Fig. 7a).The conductive anomaly location corresponds to Chang et al. (2014), inferred as a possible fracture zone for geothermal reservoirs from depths of 0 -1 km.Tong et al. (2008) identified this as cap rock, inconsistent with the Chang et al. ( 2014) results below depths of 1 km.The production casing top depths are located in the 159 -1067 m range around the resistive anomaly R1 (Fig. 7a), whereas R1 is formed by argillites and slates (Yui et al. 1993).Additionally, the B -B' profile (Fig. 7b) shows significant conductivity anomalies located at depths of 1.1 -1.7 km (C2).A ~50 ohm-m of conductive anomaly beneath the geothermal area shows similar features with a previous study that was recommended as a candidate area for drilling production wells by Tong et al. (2008).The production casing top is located mostly in the electrical resistivity discontinuity.According to the electrical resistivity model we infer that surface cold water (blue arrows shown in Fig. 7a) could flow through R1 and/or be stored in the large electrical resistivity contrast C2 (Figs. 6g, h, and 7b).The geothermal fluids (red arrows shown in Fig. 7a) are recycled to the surface outcrops along the Chingshui Fault or fracture zones under higher geothermal gradient and pressure conditions.The conductive anomaly C2 is relatively different from those in previous studies (Su 1978;Tong et al. 2008;Chang et al. 2014).Therefore, the C2 anomaly needs to be fully verified by 3D forward modeling To comfirm the C2 anomaly we replaced it with two resistivities of 10 and 1000 ohm-m with four different sizes (Table 1).The difference in response between the sizes and resistivities is shown in Fig. 8.The difference in apparent resistivity and phase curves appear to be significantly affected at four stations, 20 -23 (Figs.8a, b).The maximum difference in phase responses appears around 0.1 s periods, whereas the replaced block in model III and IV (Table 1) shows less sensitivity at all stations (Figs.8c, d).The conductive anomaly C2 size is demonstrated larger than 1.35 km 3 and would be detected by the MT data.This result indicates that C2 is well constrained and absolutely required at depths of 1.1 -1.7 km.

DISCUSSION
The conductive anomaly C2 has not yet been observed.Its presence probably indicates a degree of increased fluid interconnectivity within the slates, which indicates a higher interconnected fluid zone within the slates or metasandstones.The percentage of fluids can be deduced using Archie's law (Archie 1942;Bertrand et al. 2009;Chiang et al. 2010) as Eq. ( 1): where z denotes the effective porosity which correlates to the interconnected fluid fraction if saturated, R m is the electrical resistivity of the fluid-saturated rock, R w represents the electrical conductivity of the brine, s w is the rock saturation, "m" is the rock cementation exponent.For most sedimentary rocks it is found to be in the 1.1 -2.5 range (e.g., Doveton 1986;Ruffet et al. 1991;Guéguen and Palciauskas 1994), "n" is the saturation exponent, and "a" is the tortuosity factor.The average Na + concentration is 1122 ppm (0.1 wt%) measured from geothermal fluids produced onsite at 225°C (Fan et al. 2005) which corresponds to a solution resistivity of about 0.80 ohm-m.This unusual higher converted resistivity is inconsistent with brines under high temperature situations.This fluid could be mixed with surface water or other materials.This study set parameters of s w n -= 1, a = 1 in fluid saturation; R w is set 0.80 ohm-m measured form produced geothermal fluids; and the resistivity of 200°C at depth of 1500 m is assumed to be 0.06 ohm-m derived from the Dakhnov (1962) experiment.The degree of fluid interconnection also depends on the lower and upper bound of rock resistivity values as defined by the 3D model from 6 -24 ohm-m.This unusual low resistivity range appears in 1.5 km depth that could be caused by enhanced temperature, effective porosity or the presence of saline aqueous fluids.Based on measured data from wells the temperature at a depth of 1 -500 m is recorded to be approximately 200°C (Tong et al. 2008).This temperature may explain the observed low resistivity.The effective porosity typically decreases exponentially with depth (Rubey and Hubbert 1959) without a damage zone associated with a fault zone.Brines are widely distributed within sediment in depth based on the Taiwan tectonic setting (Chiang et al. 2008(Chiang et al. , 2010;;Bertrand et al. 2009).Saline aqueous fluids are good electrical conductors and laboratory studies have shown that brine electrical resistivity depends on the salinity, pressure and temperature (Sourirajan and Kennedy 1962).Below 300°C the resistivity of most saline fluids is relatively independent of pressure, while at higher temperatures the resistivity increases with pressure (Quist and Marshall 1968).Therefore, the low resistivity may be explained by the presence of saline aqueous fluids with enhanced heat flow below the Chingshui geothermal area.
The percentages of interconnected fluids were evaluated from the conductive anomaly C2 for different value settings and "m" values according to Doveton (1986) for typical sediment rocks.The differences between the minimum and maximum values are 1.1, 2.6, and 5.7%, respectively, derived from m = 1.1, 1.5, and 2.2 using R w = 0.06 ohm-m, whereas the differences are 11.5, 15.7, and up to 18.7% using R w = 0.80 ohm-m (Table 2).The fluid content curve percentages are shown in Fig. 9.The result indicates that Archie's equation is dominated primarily by R w and "m" values.Thus, we define m = 1.5 in agreement with previous MT research (Chiang et al. 2008)  1.8 -4.6% range (blue dashed line shown in Fig. 9).The effective porosities are consistent with typical porous sediment.Moreover, fluid resistivity of 0.80 ohm-m was measured from produced geothermal wells at 225°C (Fan et al. 2005) and used to roughly estimate the percentage of interconnected fluids from the C2 conductive anomaly concluded from 10.4 -26.1% (blue solid line shown in Fig. 9).A high heat flow zone (Lee and Cheng 1986) appears in the Chingshui geothermal area.These features imply conditions suitable for drilling geothermal energy production wells.

CONCLUSION
The MT data collected in the Chingshui geothermal region demonstrated a 3D electrical feature according to the MT tensor parameters.Therefore, 3D inversion is a challenging yet useful tool to enhance MT data interpretation in the geologically complex region.A novel 3D electrical model clearly delineated Chingshui geothermal system geometry, which provided an improved image for geological interpretations.The C2 conductive anomaly implies a deeper geothermal energy stored with high heat flow and fluid content conditions.A high heat flow zone appears in the Chingshui geothermal manifest area.The effective porosity of geological rocks is typically smaller than 10%, whereas the higher effective porosity estimations indicate influence from lower aqueous saline concentrations.This study shows that the C2 conductive anomaly could be a high potential area for geothermal energy production in the Chingshui geothermal area.Hydrothermal alteration and high temperature enhances the surface conductivity of the rock matrix, which contributes to decreasing resistivity.This is one of the important factors for precise effective porosity estimation.Geochemical evidence (e.g., pore fluid resistivity) and rock samples are required for precise fluid content estimation in the geothermal area.Table 2. Statement of interconnection fluids measured varying from "m" and R w , where "m" is cementation, and R w is the water resistivity factor.
Fig. 1.(a) Geological map of the Chingshui area overlaid with MT stations, CF: Chingshui Fault.(b) Location of the study area in northeastern Taiwan.(c) Topography map of the Chingshui area overlaid with MT stations.

Fig. 4 .
Fig. 4. Measured phase tensor ellipses (major axis) at all periods, filled with color that shows the β invariant.

Fig. 5 .
Fig. 5. MT data and its calculated Rxy and Ryx responses using a 3D inversion algorithm, where x, y indicates N -S and E -W directions.

Fig. 8 .
Fig. 8. Apparent resistivity and phase differences between forward and calculated responses.(a) Block I, (b) block II, (c) block III, (d) block IV as shown in Table1.

Fig. 9 .
Fig. 9.The curves of the interconnected fluids assessed from the 3D electrical model.The black, red, and blue lines respectively indicated m = 1.1, 1.5, and 2.2, solid and dashed lines respectively displayed in R w = 0.06 and 0.80 ohm-m.

Table 1 .
Model dimensions and resistivities.