Mapping possible subsurface granitic bodies in the northeastern Taiwan mountain belt using the VLF-EM method

Journal of Applied Geophysics 85 (2012) 25–36

Contents lists available at SciVerse ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

Mapping possible subsurface granitic bodies in the northeastern Taiwan mountain belt using the VLF-EM method Yih Jeng a,⁎, Chu-Lin Huang a, Lun-Tao Tong b, Ming-Juin Lin a, c, Chih-Sung Chen a, Hsin-Han Huang a a b c

Department of Earth Sciences, National Taiwan Normal University, 88, Sec. 4, Ting-Chou Road, Taipei, 116, Taiwan, ROC Green Energy and Environment Research Laboratories, Industrial Technology Research Institute, Hsinchu, Taiwan, ROC St. Francis Xavier High School, Taoyuan, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 30 January 2012 Accepted 21 June 2012 Available online 29 June 2012 Keywords: Taiwan orogeny Hoping geological area Granite VLF-EM EEMD

a b s t r a c t Large gneiss bodies have been reported in the metamorphic complex in northern and eastern Taiwan for decades. Some of them are cut or intruded by granitic pegmatite dikes. However, increasing evidence suggests that the gneiss bodies are more likely to be granites or meta-granites. To validate the existence of the granites/meta-granites and propose their potential distribution in the metamorphic complex of northeastern Taiwan, a geological reconnaissance along with a crooked long-distance VLF-EM survey line of 19 km and a 4.4 km controlled experimental line were conducted in the Hoping geological area of the northeastern Taiwan mountain belt. The VLF-EM data were initially processed by using the Fraser linear filter and a nonlinear filtering method based on the ensemble empirical mode decomposition (EEMD) technique to enhance the signal and to evaluate the data quality. A skin-depth added Karous–Hjelt filter was performed to generate the equivalent current density model. With the aid of the 3-D topographic representation, the equivalent current density model clearly indicates that a vast area of granites/meta-granites in the survey area is highly possible. In spite of a large uncertainty of the pseudo-quantitative model, the geological implication of our finding agrees with the tectonic framework that Taiwan and the adjacent southern Ryukyu arc system could be part of the rifted China continental margin before the collision of the Luzon and Ryukyu arcs started in late Cenozoic. © 2012 Elsevier B.V. All rights reserved.

1. Introduction Conventional geological studies of the Taiwan orogeny aim at the metamorphic complex; igneous rocks, particularly the felsic, are scarcely discussed. The geological map published by the Central Geological Survey (CGS) of Taiwan indicates roughly a distribution of gneiss with granitic dikes in the Hoping geological area, northeastern Taiwan mountain belt (Fig. 1). The distribution of igneous bodies in the northeastern Taiwan mountain belt is an important issue for geologists and geophysicists because it should have significant contribution to establish a better model of Taiwan orogeny. Besides, the igneous body such as the granite is valuable for the construction industry in Taiwan due to its scarcity. Some geochemical and tectonic studies have been carried out on Taiwan granitic rocks (Ho, 1988; Jolivet et al., 1990; Lan and Wang Lee, 1992; Qiu et al., 2011); however, the study of granite distribution using geophysical methods in the northeastern Taiwan mountain belt is sufficiently new that little or no related published literature describing specific details is available. A recent geochronologic study

⁎ Corresponding author. Tel.: +886 2 77346416; fax: +886 2 29333315. E-mail address: [email protected] (Y. Jeng). 0926-9851/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2012.06.010

for the rock samples selected from the Tailuko metamorphic belt in northeastern Taiwan revealed the existence of the 191 ± 10 Ma Talun meta-granite, the oldest granitic intrusion ever reported in the Taiwan region and along the eastern coastal area of south China (Yui et al., 2009). This finding may suggest that the northeastern Taiwan mountain belt has close relation with the rock strata of the Eurasian continent margin if the distribution of the granite is significant. Although further studies of this issue are undergoing currently, the results are preliminary and have not been well examined. A major problem hindering the geological investigation of this area is that the rugged terrain and uncultivated virgin forest make any access extremely difficult. Most parts of the northeastern Taiwan mountain belt are investigated only with brief geological reconnaissance; no systematic detailed mappings are available. The main objectives of this investigation are to examine the granite/meta-granite existence and the potential distribution in a selected geological section of northeastern Taiwan, to interpret their possible origin, and to relate their importance to the orogeny of Taiwan. Considering the difficulty of performing geophysical field work in this rugged terrain, we used the very low frequency electromagnetic (VLF-EM) instrument as the main tool for field survey due to its light weight, easy operation, and environmentally friendly attributes. The selected study region is located in the Hoping geologic area where metamorphic complex

26

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

Fig. 1. Tectonic setting of Taiwan (Teng, 1996) and the location of the study area. The survey lines are displayed in a 3-D coordinate system over the plane sketch of the proposed subsurface granite belt. A perspective shaded relief view of the study area is shown in Fig. 2. The map is compiled in Transverse Mercator 2-degree projection (TM2) using TWD97 (Taiwan Datum 1997) in meters, which is the coordination system currently used in Taiwan.

with igneous dikes (mainly granitic pegmatite and diabase dikes) are most likely to exist (Fig. 1). The field survey was carried out along the Hoping forest trail, traversing roughly from east to west in the Hoping geological area with a total survey length of about 19 km. A controlled survey was conducted along Hoping creek, about 4 km north of the Hoping forest trail where marble/granite contacts are clearly visible. In this study, we analyzed the acquired field data by using three different filtering techniques, and the subsurface earth model is obtained through a 2-D equivalent current density distribution. The results propose a possible detailed spatial distribution of the granitic bodies, marble and schist, and marble/granite contact in the study region. We believe that the subsurface earth model deduced from our findings would provide a better understanding of the northeastern Taiwan orogeny that may draw the attention of geologists for further investigations. 2. Geological settings of the area The island of Taiwan is situated on the boundary between the Philippine Sea plate and the Eurasian plate. The Philippine Sea plate has been moving northwest (Jolivet et al., 1990) and developing two volcanic arcs, the Luzon arc in the south by overriding the Eurasian plate and the Ryukyu arc in the northeast by subducting beneath the Eurasian plate (Fig. 1). Consequently, the Taiwan mountain belt is an active collision zone between the Philippine Sea plate and the Eurasian plate. Teng (1996) proposed a hypothetic model that Taiwan and the neighboring southern Ryukyu used to be part of the rifted Eurasian plate margin before the collision started in late Miocene. Tectonized rock strata of both the China continental margin and the Luzon arc are exposed in the orogenic belt as stacked west-facing folds and thrust sheets trending north-northeast– south-southwest (Ho, 1988). Based on investigations of the CGS of Taiwan, the pre-Tertiary metamorphic complex exposing mainly in the eastern part of the Central Range is the oldest geological element of Taiwan (Central Geological Survey of Taiwan, 2012). From the regional perspective and the tectonic history, the Taiwan orogenic belt should comprise pre-Tertiary metamorphic complex and igneous rocks of the deformed continental crust (Teng, 1996); however, the latter is not well recognized due to the lack of field evidence. In contrast to the compressional tectonism of Taiwan collision orogeny, the Ryukyu arc is developed by extensional tectonism. The Okinawa trough, an active back-arc rifting basin behind the Ryukyu

arc, develops upon the thinning crust of the original orogeny and runs through the Taiwan orogenic belt at the Ilan Plain (Fig. 1). As a consequence, the northern Taiwan orogenic belt is collapsing to give way to the rifting Okinawa trough. Therefore, igneous rocks of the deformed continental crust should be more evident in the northeastern Taiwan mountain belt. Yen (1963, 1981) and historical reports of the CGS of Taiwan indicate that the greenschists, marble, and gneisses are important members of the metamorphic complex in northern and eastern Taiwan, and are exposed extensively throughout the area. The gneiss bodies with sparse plutonic intrusions are scattered in different places in northeastern Taiwan near the eastern coast. In the northeastern Central Range, large gneiss bodies are found in the pre-Tertiary metamorphic complex. Our study region is located in the Hoping geological area of the eastern Central Range, northeastern Taiwan, approximately between 30 km and 35 km south of Ilan Plain (Fig. 1), on the border of Ilan and Hualien Counties with rugged topography (Fig. 2). The geology of this area belongs to the pre-Tertiary metamorphic complex comprising mainly marble, gneisses, and pegmatite. North of our study region, between Tungao and Nanao of Ilan county, the largest gneiss body in northeastern Taiwan was roughly mapped (Central Geological Survey of Taiwan, 2012; Ho, 1988). This gneiss body approximately extends 16 km east to west and 3 km north to south. In our study region, outcrops of large gneiss bodies and some smaller

Fig. 2. Shaded topographic map of the Hoping geological area with the survey lines indicated.

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

ones are also mapped by the CGS of Taiwan. The gneiss bodies usually are found to be cut or intruded by numerous granitic pegmatite dikes, diabase dikes, or quartz veins (Central Geological Survey of Taiwan, 2012). Lan and Wang Lee (1992) conduct mineral chemistry study on the metamorphic rocks in northeastern Taiwan and speculate that the gneiss and schist in the area possibly originate from granite. A petrologic investigation conducted by Qiu et al. (2011) confirms the possibility that granite is the original rock that formed the metamorphic rocks in the Hoping area. Recent geological surveys funded by Tai-Power Company (L.-T. Tong, personal communication, July 21, 2011) and the National Science Council of the ROC in the Hoping area (Qiu et al., 2011) reveal that granites, meta-granites, basalt and meta-basalts are scattered in different places. On the north side of Hoping creek, large xenoliths basalt may exist in the metamorphic basement (L.-T. Tong, personal communication, Oct. 19, 2011). However, the results are tentative and have not been reviewed and published.

27

The VLF-EM instrument measures the two components by reading the magnetic signals detected by a reference coil and a signal coil, the two coils being perpendicular to each other. Theoretically, these two components are derived from the tilt angle θ and the ellipticity (or eccentricity) ε of the polarized elliptic resultant field measured. Two simple equations given below show the relationships (Paterson and Ronka, 1971). The tilt angle θ is the inclined angle of the minor axis of the ellipse to the vertical, and can be expressed as: −1

θ≈ tan


 
 S S sinα cosϕ ≈ sinα cosϕ P P

ð1Þ

where α is the inclined angle of the cylindrical secondary field S to the plane of the primary field. Eq. (1) indicates that tilt angle θ approximates the inphase component (S cos ϕ) of the vertical secondary field. The relationship in Eq. (1) is valid when the value of θ is small. The ellipticity ε is defined as the ratio of the minor axis to the major axis of the polarized elliptic field, and can be represented by:

3. VLF-EM method The VLF-EM technique is an electromagnetic method in which the signal sources are radio waves in the bandwidths of 15–30 kHz from worldwide network radio transmission stations. This technique is a passive electromagnetic method that has been extensively used in geophysical exploration for various purposes such as ore body detection, groundwater survey, shallow fault investigation, geological mapping, and landform studies (Al-Tarazi et al., 2008; Bernard and Valla, 1991; Fischer et al., 1983; Jeng et al., 2004; Lin and Jeng, 2010; Paál, 1965; Paterson and Ronka, 1971; Van Dam, 2012). The signals from VLF transmission stations propagate as spacewave (signal goes directly from transmitter to receiver through space), skywave (signal is refracted and reflected by the ionosphere and return to the ground), and groundwave (wave that travels along the surface of the earth). At long distances from the transmission station, we receive mainly the skywave as the primary wave (Paterson and Ronka, 1971). The VLF primary wave P penetrates the ground surface and generates a phase-shifted secondary field S when coupling with buried conductors. The VLF instrument measures the magnetic field, which is the resultant magnetic field of the primary and the secondary, rotating in space and varying in magnitude as it propagates. Therefore, it describes a polarized elliptic field (Grant and West, 1965). Since the instrument cannot differentiate the primary and secondary fields, the VLF-EM survey actually made in the field is a result of measuring the tilt angle and the ellipticity (or eccentricity) of the polarized elliptic resultant field, then the desired secondary field can be estimated through approximation. The procedure is briefly illustrated by applying the basic electromagnetic (EM) theory as follows. It is well known that the phase difference φ between P and S is the phase lag due to induction plus the phase shift caused by the electric properties of the subsurface substance. For a very good conductor, the phase difference φ equals about π. Conversely, it equals approximately π/2 for a non-conductive body (Telford et al., 1990). In most cases, the phase difference φ ranges between π/2 and π. A common practice is to decompose the secondary field into inphase component (S cos ϕ) and quadrature component (S sin ϕ) to simplify the data processing and interpretation. By the EM convention, the ratio of the quadrature (imaginary) component to the inphase (real) component is called scalar tipper. Some authors call the inphase and quadrature components as real and imaginary components of the tipper, respectively (Beamish, 1994; Monteiro Santos et al., 2006). We use the terms real and imaginary for the convenience of applying complex algebra to describe these two components, and they should not be the same as those defined as the imaginary and real numbers in mathematics.

ε≈

S sinα sinϕ P

ð2Þ

The above two equations reveal that the inphase and quadrature components of the secondary magnetic field measured are ratios of the vertical secondary field to the primary field; therefore, these two components in the case of the VLF-EM method are usually expressed as percentage. With the understanding of the VLF-EM method, some basic considerations must be taken in account: 1. The secondary field may not be unique; in other words, the magnetic field to be measured in the survey area is a resultant field of the primary with one or more secondary fields and may be contaminated by other magnetic fields depending on the conditions of the survey environments. Therefore, appropriate filtering processes are required to avoid false interpretation. 2. Since the tilt angle and ellipticity of the polarized elliptic field are measured to approximate the inphase component and quadrature component of the vertical secondary field, respectively, the magnitude of the secondary field will affect the error of measurements. The error of the inphase measurement increases as the magnitude of the secondary field increases, whereas the error of the quadrature measurement decreases (Jeng et al., 2007). 3. The penetration depth of VLF signals varies with the frequency (f) of the signal and the resistivity (ρ) of the medium. A simple criterion for the penetration of EM waves is the skin depth (δ), the distance in which the amplitude of wave is reduced by a factor of 1/e, or 0.368. This is given by Eq. (3) in the unit of m (meter), 1=2

δ≈500ðρ=f Þ

m

ð3Þ

The skin depth is an important factor to consider when the VLF-EM method is applied in an area with overlying conductive layers. Eq. (3) shows that resistivity is essential in determining the skin depth. In our field measurements, the outcrop resistivities acquired by using a portable resistivity meter vary between 257 and 3367 Ω-m (1204 Ω‐m on average) for granite/meta-granite. Resistivities of marble and schist can be higher than 1000 Ω‐m; however, they depend on the degree of fractures and water contents in the rock. Other resistivity data are available in our previous audio magnetotelluric (AMT) survey along the same VLF survey line. The AMT frequency ranges from 10,000 to 0.1 Hz with the penetration depth of about 100 m to a few kilometers, and the resistivity measured in the area hosted by marble is about 500 Ω‐m. We also

28

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

performed the EM survey by using the HLEM (horizontal loop electro-magnetic) method which penetrates to a depth of about 6 m. The measured resistivity range is between 43 and 5000 Ω‐m (707 Ω‐m on average) with about 820 Ω‐m for areas hosted by marble/schist and 650 Ω‐m for granite/gneiss. Based on the resistivity data obtained by these methods, the minimum skin depth estimated in our VLF-EM survey should be about 100 m and may reach 300 m in the highest resistivity zones. 4. Data processing methods We applied two standard VLF-EM linear filtering methods and one newer nonlinear data analysis technique, ensemble empirical mode decomposition (EEMD), to process the data. The standard filtering results can be improved with an aid of the EEMD technique. For applying the linear inverse filer to calculate the pseudo-current density for geological interpretation, attenuation with depth based on the skin-depth is incorporated into the program with the resistivity obtained in our previous MT survey. 4.1. Linear filtering methods The most commonly employed VLF-EM data processing techniques are the Fraser filter (Fraser, 1969) and the Karous & Hjelt linear filter (Karous and Hjelt, 1983). The Fraser filter first applies a simple difference operator to transform zero-crossings into peaks (shift the phase of the data by 90°), and removes dc bias to increase resolution of local anomalies. The difference operator is similar to a differentiator (quadrature filter) in basic signal processing (Claerbout, 1976), which converts sines into cosines, but those higher frequency amplitudes are amplified with respect to lower frequencies. For example, let the zero-crossing represent a sine function of distance s given by H = c + H0 sin ks with amplitude H0 and frequency (wave-number) k and dc bias c. By applying the differentiator to the sine function we have ΔH dH ≈ ¼ H 0 k cosks Δs ds

ð4Þ

Obviously, the phase is shifted 90° and the dc bias is removed. Notice that the higher frequency (wave-number) amplitudes are amplified with respect to lower frequencies (wave-numbers) and therefore a low-pass filtering is needed for suppressing the amplified high frequency components. The Fraser filter implements the difference operator by taking the sum of two consecutive data points, and then subtracting it from the sum at the next two consecutive data points. To prevent it from amplifying the high frequency random noise, a smoothing operator is included in the Fraser filter. Based on the VLF-EM theory, the zero-crossing may indicate the center of the subsurface anomalous body; therefore, advantages of using the Fraser filter are enhancing the zero-crossings and making the contouring easier for noisy 2-D data. One advantage that is seldom noticed in applying the Fraser filter is the ability of removing topographic effect because this filter uses the discrete first derivative as one of its components (Fraser, 1969). The Fraser filter is easy to apply; however, the results are qualitative and the subsurface details are unavailable. Other filtering and inversion techniques should be employed to acquire more detailed subsurface structures. The Karous and Hjelt linear filter is an extension of the Fraser filter with a broad-band design, and has a function of inverse filtering because the output of the filtered data is expressed by a subsurface equivalent current density pseudo-section which would generate the measured magnetic field. By assuming the magnetic field measured on the earth surface is caused by the current density underground and by using the Biot–Savart law, Karous and Hjelt (1983) propose that the vertical magnetic field Hz(x) measured on the

earth surface at a distance of x due to the current density Ia(ξ) at horizontal distance ξ and a depth of Δz is H z ðxÞ ¼

∞ h i 1 2 2 ∫ Ia ðξÞ⋅Δz⋅ðx−ξÞ⋅dξ= ðx−ξÞ þ z 2π −∞

ð5Þ

where x and ξ are the distances measured from the origin of the reference coordinates. The current density Ia can be solved by using the linear inverse method as follows. At first, we simplify the conditions that the current density distribution to be calculated is located at points with depth z = Δx. Then the discrete form of Eq. (5) is H zm ðxi Þ ¼

 ∞      2 1 X 2 I a ξj ⋅Δz⋅Δx⋅ xi −ξj = xi −ξj þ Δx 2π j¼−∞

ð6Þ

Note that Hzm(xi) represents the discrete Hz(x) measured at points with equal interval Δx and xi = i ⋅ Δx. Because the location of the current density element is represented as a current line within a square area of Δx ⋅ Δz at depth z; therefore, ξj = (j + j0) ⋅ Δx and 0 ≤ j0 b 1. Let Hi ¼ 2π Δz H zm ðxi Þ and Ia(ξj) = Ij, Eq. (6) can be re-written more compactly with finite number of coefficients as Hi ¼

n X

Ij ⋅K ij with K ij ¼

j¼−n

i−j−j0 : ði−j−j0 Þ2 þ 1

ð7Þ

Finally, the current density Ij can be obtained by using the linear inverse filter method and the solution is Ij ¼

nþ1 X

K i0

−1

⋅H iþj

ð8Þ

i¼−n

In practice, the inverse filtering is compromised by restricting an infinite number of filter coefficients to a finite number. Experiments were conducted by the inventers and others (Pirttijärvi, 2004) with finite filter coefficients to determine an optimum filter used for the inversion. The shortest filter as we used in the study is a six-coefficient filter giving less than 8% error for the inversion of the field data of a single current line (Karous and Hjelt, 1983). The operation of the Karous and Hjelt linear filter approximates to that of the difference operator used in the Fraser filter; therefore, the data should be smoothed to avoid amplifying high frequency noise. By using the data measured at different spacing (e.g., Δx , 2Δx, 3Δx), we can calculate the variation of current densities with depth. The larger spacing used, the deeper current density Ij is obtained, and the total horizontal length of the current density distribution obtained at that depth becomes shorter accordingly. The complete current density pseudo-section is constructed by stacking the results calculated at each level of depth one below the other so as to form a reverse trapezoid pattern if no extrapolated data points are added. The finite-length filter is useful in VLF-EM data processing by giving a general model of the subsurface structure in 2-D cross-section, but the pseudosection does not necessarily correspond to true current distribution. Improved filters can be obtained by adding the attenuation with depth to calculate a more realistic model. Although the ambiguity problems of 2-D inversion may still exist (Beamish, 1998; Loke and Barker, 1996), the Karous and Hjelt linear filter is an inverse filter that is easy to apply and effective in giving an estimated model of the subsurface structure. 4.2. Nonlinear filtering method There are a number of factors that can seriously degrade the VLF-EM signal, particularly the geological noise and background EM waves. The geological noise is the responses from geological structures that are on a

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

scale too small or too detailed compared with our model. It is a signal generated broadband noise, depending on the causative geological structures. The background EM waves indicate all other harmonic and nonlinear noises from remote sources and environments. These noises usually are non-stationary and nonlinear, and are difficult to remove by using conventional linear filtering methods. For removing the aforementioned noises, we applied an adaptive nonlinear and non-stationary data analysis scheme, the EEMD method, to process the data. The EEMD technique is a noise-assisted empirical mode decomposition method that is more effective and accurate to extract the signal and to eliminate the noise of VLF-EM data than its earlier version EMD (empirical mode decomposition), particularly in alleviating the mode mixing problem (Huang and Wu, 2008; Jeng and Chen, 2011; Lin and Jeng, 2010). We only briefly describe this newly developed empirically based adaptive technique here, and more details can be found in the related literature (Huang and Wu, 2008; Huang et al., 1998; Wu and Huang, 2009). The key concept of the EMD data analysis technique is based on decomposing the data into a finite number of orthogonal simple oscillatory modes called intrinsic mode functions (IMFs) (Huang and Wu, 2008; Wu and Huang, 2009) which fulfill two conditions: (1) The number of extrema and the number of zero crossings must be equal, or differ at most by one; (2) The envelope defined by the local maxima and that by the local minima are symmetrical to zero. The IMFs are unlike basis functions of the Fourier analysis which are mono-frequency and stationary, instead they have time-variable amplitudes and frequencies with physical significance, and yield meaningful instantaneous time–frequency attributes of the data after the Hilbert Transform is applied to each IMF component. Mathematically, the IMFs are obtained by carrying out a ‘sifting’ procedure described extensively by the inventors and others (Huang et al., 1998; Jeng et al., 2007). The sifting procedure is first to establish the upper and lower envelopes for the original data, and then to determine the mean of the two envelopes. The second step is to subtract the mean from the data to get the difference d1. The difference d1 is then treated as the input data for the next cycle of sifting. The sifting procedure is performed iteratively n times until the difference d1n satisfies the conditions of IMF. Then d1n is called the first IMF component c1 of the original data. Usually, c1 has the finest scale or the highest frequency component of the data. We then subtract c1 from the data, and the residue is treated as the new input data to obtain the next IMF component. The sifting procedure can be repeated to obtain all the subsequent IMFs until a predetermined stoppage criterion is reached (Huang and Wu, 2008; Jeng and Chen, 2011). Consequently, the original data set x(t) is decomposed into a series of IMFs and a residue as xðt Þ ¼

n X

cj þ r n

ð9Þ

j¼1

where cj, j =1,…,n are the n IMF components from highest to lowest frequency band, and rn is the residue after repeating n times of the

29

decomposing procedure. The EMD is effective in extracting signals from data generated in noisy nonlinear and nonstationary processes (Jeng et al., 2007; Wu and Huang, 2009); however, it still has some difficulties when applied to real geophysical data. One weakness of the EMD technique as stated previously is that the decomposed components may have mode mixing problem, i.e. a signal or signals of similar scale may reside in different IMF components. The mode mixing problem is particularly serious for data consisting of intermitting signal. Since the acquired VLF-EM data are very likely intermittent due to coupling of the transmitting signal and the targets, we used the EEMD rather than EMD method in this study. By reviewing the related publications (Chen and Jeng, 2011; Huang and Wu, 2008; Jeng and Chen, 2011; Lin and Jeng, 2010; Wu and Huang, 2009), the EEMD method is carried out by adding a finite amplitude white noise to the input data set x(t) before sifting process: X ðt Þ ¼ xðt Þ þ nðt Þ  R

ð10Þ

where X(t) is the so called noise added data, n(t) is the added white noise series, and R is the ratio of the standard deviation of the added noise amplitude to that of the data x(t). The EMD procedure is applied to data X(t) rather than x(t). By repeating the procedure (adding the noise to the data and performing EMD) with sufficient times of different white noise series of the same amplitude (i.e. R keeps constant), we obtain an ensemble of corresponding IMF as: Ei ¼

k X

cij

ð11Þ

j¼1

where Ei is the ensemble of the ith IMF component containing k members; cij is the jth noise added trial of the ith IMF component of X(t). The ith EEMD component is obtained by taking the mean of the IMFs in ensemble Ei: ci ¼

k 1X c k j¼1 ij

ð12Þ

where ci denotes the mean of the ith IMF component, and k is the number of trials for that ensemble. Notice that c i is not a true IMF except that k approaches infinity. Theoretically, the white noise added in each trial will be suppressed significantly if the noise is added sufficient times (a few tens to one hundred at least). Therefore the advantage of using EEMD is that we use the added white noise to provide uniformly distributed reference scales in the whole time-frequency space for the signal to reside, but the noise added in each trial will cancel out each other after averaging. As a result, we alleviate the mode mixing problem of EMD. This method can therefore be used for removing nonlinear noises in VLF-EM data before inverse modeling. Another important reason of employing the EEMD filter in this study is to separate disturbing influences from signals on the

Fig. 3. Pictures of outcrops taken in the vicinities of Stations 8 and 31 on the Hoping creek survey line.

30

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

Fig. 4. Original and Fraser filtered data of three different transmitter frequencies acquired on the Hoping creek controlled survey line. Station numbers start from east to west.

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

31

Fig. 5. EEMD analyzed 17.4 kHz data acquired on the Hoping creek controlled survey line. The original data are decomposed into four IMF components and a residue. (a) Inphase data. (b) Quadrature data.

long survey line. Due to the travel and logistic difficulties came upon on the long survey line in the mountain area, the crew members sometimes had to reverse the direction for a certain segment on the survey line or to finish one survey line in two working days. These actions usually cause discontinuities in the data and fortunately, the disturbances can be removed easily by examining the EEMD decomposed data to determine and eliminate the noise components (Jeng et al, 2007; Lin and Jeng, 2010). A more reliable VLF-EM data set can be attained consequently. 5. Field survey and procedure The VLF-EM survey was conducted with a GSM-19 system. To validate the existence of granites, which mostly are pegmatite and meta-granite, and to map the subsurface structures of the study region, we deploy the VLF-EM survey line along the Hoping forest trail that traversed the gneiss belts in which the granite distributions are inferred with our previous geologic reconnaissance. Signals from three VLF transmitter stations of

17.4 kHz (NDT in Yosamai, Japan), 19.6 kHz (GBZ in Oxford, UK), and 22.3 kHz (NWC in NW Cape, Australia) were received with excellent quality in the study region. Based on the Faraday's law in the EM theory, the VLF-EM response is best as the primary magnetic field which is perpendicular to the cross-section of the subsurface targets; therefore, the transmitter should be in the direction of the geological strike. The locations of the three transmitters basically fulfill the requirement as they are located in the far north and south, and the strikes of the Taiwan mountain belt are approximately trending northeast–southwest. The VLF-EM data of the Hoping forest trail were acquired at 190 stations in two days with the station interval of 100 m. To ensure the data quality for the interpretation of geological significance, a total length of about 4.4 km VLF-EM controlled survey line was deployed along the channelbed of Hoping creek from east to west with the same station interval as used on the Hoping forest trail (Fig. 1). Stations 8 and 31 on the Hoping creek survey line exhibit clear marble/granite contacts in close proximity (Fig. 3).

32

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

Fig. 6. Comparison of EEMD filtered data with the Fraser and EEMD filtered results. The EEMD filtered data are obtained by eliminating IMF component c1 and the residue. The Fraser and EEMD filtered results indicate that the data are EEMD filtered prior to the Fraser filtering.

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

33

Fig. 7. EEMD filtered inphase current density pseudo-sections of Hoping creek with topographic variation. (a) 2-D profile. (b) 3-D demonstration in TM2 coordinates with elevation. The skin depth effect is not considered in the inverse filtering, and the sample interval is 100 m as laid out in the field. (c) Same as (b) but the skin depth effect is applied.

Therefore, the measured data at these two stations are used as the control points for geological interpretation. 6. Results The data are interpreted with the aid of inspecting the single profiles and inverted current density models. Three modes of single profiles are demonstrated: the single trace of raw data, the Fraser filtered data, and the EEMD processed data. A preliminary idea of the distribution of the subsurface lithological contacts can be acquired by looking over the zero-crossings in the raw data. Then, the Fraser filter enhances the signal by removing the topographic effect and bias in addition to converting the zero-crossings into peaks. For data intervened with nonlinear or non-stationary noises, the EEMD filter should be employed. The current density pseudo-section is the final presentation of the data. It gives the possible subsurface distribution of the rocks. The high current density zones may correspond to the marble or fracture zone with water, and the places of low current densities should indicate the granites or meta-granites. 6.1. Hoping creek controlled survey line The original and Fraser filtered data of this survey line with three different transmitter frequencies are shown in Fig. 4. By examining the single profiles of raw data, the 17.4 kHz and 19.6 kHz inphase data are seriously biased and the 22.3 kHz data are relatively noisy, so that we should look at the Fraser or EEMD filtered data which the bias could be removed and signal enhanced. The inphase data of the first two

profiles, especially the one using the 17.4 kHz signal, show abnormally negative values, and the quadrature response is excessively flat. We suspect that the data are strongly affected by the background field and local clutter noise. This inference is verified by examining the EEMD analyzed data (Fig. 5) which show the residue component drifts with high background energy of long wavelength, and the first component c1 contains abundant local clutter energy of short wavelength. The reason is that the filter bank derived from the EMD or EEMD analysis is dyadic, and the components start from the highest (c1) to the lowest frequency band (residue) in sequence (Huang et al, 1998; Wu and Huang, 2009). The component of the highest frequency band is probably not physical or not the response that we are interested in (Huang et al., 1998; Jeng et al., 2007). A complete signal-noise analysis for the EEMD filter bank can also be found in Chen and Jeng (2011). To recover the usable data, the noisy component c1 and the residue should be eliminated to reconstruct the profiles for interpretation. In some cases, the Fraser filter also removes the background bias and clutter noise if the signal and noise are linearly combined. For this reason, we first tried the Fraser filter on the raw data and found that the 17.4 kHz data are greatly improved while the other two are still noisy. Therefore more delicate filtering process may be required. Fig. 6 shows the EEMD filtered results, and the results obtained by using Fraser filter after EEMD reconstruction. The 17.4 kHz data still demonstrate the best quality after further filtering process. We also found the consistency of these two filters by comparing Fig. 4 with Fig. 6. The two prominent events with geological significance in the Fraser filtered data (Fig. 4) are also revealed by EEMD filtering, and are further improved by applying Fraser filtering right after EEMD reconstruction.

34

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

Fig. 8. Resampled current density pseudo-sections of Hoping creek inphase data in 3-D display. The data were up-sampled to 25 m intervals and EEMD filtered. (a) The skin depth effect not considered; (b) the skin depth effect considered.

According to the basic theory we stated previously, the zerocrossings of raw data or the peaks in the Fraser filtered data are related directly with centers of the subsurface anomalous bodies. More information can be inferred by reading the raw data if the data quality is good enough. For instance, the distinct rise of the inphase amplitude may indicate a large anomalous body. The asymmetry of the amplitude in the in-phase and out-of-phase components is caused by the dipping effect of the structure, and the larger anomaly peak gives the downdip side while the reduced broader peak is on the updip side (Coney, 1977; Telford et al., 1990). The horizontal distance between the positive and negative peaks of the inphase profile can be used to estimate the depth of the anomalous body (Paterson and Ronka, 1971). Inspecting the two significant events in the filtered data of 17.4 kHz, one is located in the vicinity of Station 8, and the other is near Station 31. Referring to the outcrops shown in Fig. 3, the two

events should indicate the subsurface contacts of marble and granite. Near Station 31, the larger positive anomaly peak on the west side as compared with the reduced broader negative peak on the east side indicating the contact may dip to the west which is coincident with the outcropped. To acquire necessary geological information for interpreting the VLF-EM data, we referred to the CGS geological map of the Hoping area (Central Geological Survey of Taiwan, 2012) and conducted a brief geological mapping on the controlled survey line. The collected information indicates that the dominant rocks observed between Stations 1 and 8 are mostly marbles, Stations 8 and 31 are granite and meta-granite, and Stations 31 and 44 are marbles and schist. The event near Station 8 shows a close correlation with the observed granite/marble contact which may indicate the eastern margin of the granite, whereas the peak near station 31 thus corresponds to the marble/granite contact which may signify the western margin of the granite. This conjecture is confirmed by comparing the Karous &

Fig. 9. EEMD filtered 2-D profiles of the inphase current density pseudo-sections of the Hoping forest trail with topographic variation. Skin depth effect is not considered. (a) Data acquired from transmitter station with 17.4 kHz. (b) Data acquired from transmitter station with 19.6 kHz.

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

35

Fig. 10. EEMD filtered current density pseudo-sections of the Hoping forest trail 17.4 kHz inphase and quadrature data in 3-D display. (a) The skin depth effect not considered; (b) the skin depth effect considered. The TM2 coordinates of Stations 110 and 120 are (319814, 2688625) and (319667, 2688226), respectively.

Hjelt linear filtering result in the latter with the EEMD or the Fraser and EEMD filtered data shown in Fig. 6. A more detailed subsurface model of the survey area can be obtained by utilizing the Karous and Hjelt linear filter. The output of this filter is expressed by relative current density pseudo-sections. The representation of the filtered output data is in an opposite sense to the resistivity commonly used in the EM data processing, i.e., lower values of the current density indicate higher values of resistivity. Fig. 7 shows the subsurface relative current density distribution of the Hoping creek controlled survey line with topographic variation both in 2-D profile and in 3-D coordinates. Because the current density is inversely proportional to the resistivity, lower current density may indicate high resistivity and vice versa. With respect to the rock resistivity, granites are much higher than the schist and marble by about one to four orders of magnitude if no water bearing in the rocks (Telford et al., 1990);

therefore, the lower current density zones are inferred to be granites/ meta-granites, and higher values could be associated with marble or schist as the rock types in the study area are granites, greenschists, marble, and gneisses (Ho, 1988; Qiu et al., 2011; Yen, 1981). The shallow (within roughly 50 m depth) distinct subsurface current density contrast in the vicinity of Station 8 and a contact dipping to the west near Station 31 are observed. Since the depth of VLF-EM investigation is typically less than 100 m for conducting bodies which is the minimum skin depth in our case, we considered the skin depth effect and up-sampled the data to 25 m intervals for comparison. The result is shown in Fig. 8. The current density distributions generated by using original data and resampled data (Figs 7 and 8) are basically in good agreement in shallow depth (within 100 m). More detailed structures are recovered in the resampled section; however, it must be interpreted carefully because artifacts may exist.

36

Y. Jeng et al. / Journal of Applied Geophysics 85 (2012) 25–36

6.2. Hoping forest trail survey line With the success of the controlled survey, the data acquired from the Hoping forest trail are processed with the similar procedure as applied on the controlled site. Fig. 9 demonstrates the 2-D profiles of current density pseudo-sections of the Hoping forest trail with elevation. The c1 component with high frequency clutter noise and the low frequency residue were removed in the EEMD data reconstruction to enhance the signal of geological interest. The section derived from 17.4 kHz data correlates well with the one derived from 19.6 kHz data. In general, the current densities inverted in this area are relatively low except a relatively high current density zone is found at around between 11 km (Station 110) and 12 km (Station 120); this could be a water bearing fracture zone because landslides are observed at 11.5 km (Station 115), 11.6 km (Station 116) and 12 km (Station 120). Landslides are also observed at Stations 175 (17.5 km) and 176 (17.6 km); however, the current densities inverted are not as high as at previous three stations. This may indicate that the water saturation rate is lower at these two stations. The overall relative low current density found in the inverted section may be a sign of a vast distribution of granite/meta-granite in this area. However, by referring to the result obtained from the controlled site, the current density distributions below 100 m depth could be impractical. Fig. 10 shows the current density pseudo-sections of the Hoping forest trail with topographic variation in 3-D coordinates. The quadrature profiles show low current density distribution and insignificant contrast throughout the line with or without taking the skin depth effect, which may indicate the poor conductivity of the subsurface. With this evidence and the previous discussion regarding the inphase profiles, the granitic bodies could be the dominant rock type in this area. 7. Conclusions As located in an area of subduction zone, granitic bodies were rarely discovered in the Taiwan orogenic belt and seldom discussed in the literature. With more evidence found by recent geological investigations, it is almost certain that the main rocks, gneiss and schist in this area originate from granite. In this study, we provide the geophysical evidence showing that the granite distribution in the Hoping geological area may not be limited to the spots with obvious outcrops. We suspect that granite and meta-granite are the main rock types in the northeastern Taiwan mountain belt. However, the origin of the granite in this area is still an interesting issue. From the tectonic point of view, we may infer that northern part of the Taiwan collision orogeny is transferred to subduction and results in losing the compressional stress. Consequently, the orogeny of northeastern Taiwan is subjected to dilatational stress. The depression of Ryukyu trench stretches the overriding lithosphere and causes the collapse of the collision orogenic belt in northeastern Taiwan. Recent geological studies of the Hoping geological area indicate that basalt and meta-basalt are the dominant rocks north of Hoping creek while granite and meta-granite are found in the south. With the evidence of the granite outcrops we found in the study region and the near surface relatively high resisitivity distribution inferred from our VLF-EM data, we may arrive at a conclusion that the northeastern Taiwan mountain belt could be part of the deformed Eurasian plate margin, and the basalt or meta-basalt could be the rocks of the deformed oceanic crust exposed due to the collapse of the overriding orogenic wedge. The details should be taken into consideration in future studies. Acknowledgments This research was financially supported partly by the National Science Council of Taiwan, ROC, Grant No. NSC 100-2116-M-003 –005–. The authors appreciate the constructive comments of two anonymous reviewers. We are also grateful to the field team of Green Energy and

Environment Research Laboratories, Industrial Technology Research Institute, Taiwan for helping in the data acquisition. References Al-Tarazi, E., Abu Rajab, J., Al-Naqa, A., El-Waheidi, M., 2008. Detecting leachate plumes and groundwater pollution at Ruseifa municipal landfill utilizing VLF-EM method. Journal of Applied Geophysics 65, 121–131. Beamish, D., 1994. Two-dimensional, regularised inversion of VLF data. Journal of Applied Geophysics 32, 357–374. Beamish, D., 1998. Three-dimensional modeling of VLF data. Journal of Applied Geophysics 39, 63–76. Bernard, J., Valla, P., 1991. Groundwater exploration in fissured media with electrical and VLF methods. Geoexploration 27, 81–91. Central Geological Survey of Taiwan, MOEA, 2012. Geology of Taiwanaccessed May 7, 2012 at http://www.moeacgs.gov.tw/english/twgeol/twgeol_eastern_1.jsp2012. Chen, C.-S., Jeng, Y., 2011. Nonlinear data processing method for the signal enhancement of GPR data. Journal of Applied Geophysics 75, 113–123. http://dx.doi.org/ 10.1016/j.jappgeo.2011.06.017. Claerbout, J., 1976. Fundamentals of Geophysical Data Processing: with Applications to Petroleum Prospecting. McGraw-Hill, New York. Coney, D.P., 1977. Model studies of the VLF-EM method of geophysical prospecting. Geoexploration 15, 19–35. Fischer, G., Le Quang, B.V., Muller, I., 1983. VLF ground surveys: a powerful tool for the study of shallow two-dimensional structures. Geophysical Prospecting 31, 977–991. Fraser, D.C., 1969. Contouring of VLF-EM data. Geophysics 34, 958–967. Grant, F.S., West, G.F., 1965. Interpretation Theory in Applied Geophysics. McGraw-Hill, New York. Ho, C.S., 1988. An Introduction to the Geology of Taiwan: Explanatory Text of the Geological Map of Taiwan, second edition. Central Geological Survey, Ministry of Economic Affairs, Taipei, Taiwan. Huang, N.E., Wu, Z., 2008. A review on Hilbert–Huang transform: the method and its applications on geophysical studies. Reviews of Geophysics 46. http://dx.doi.org/ 10.1029/2007RG000228 (RG 2006. Huang, N.E., Shen, Z.S., Long, R.M., Wu, C., Shih, H.-H., Zheng, Q., Yen, N.-C., Tung, C.-C., Liu, H.-H., 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis. Proceedings of the Royal Society of London Series A 454, 903–995. Jeng, Y., Chen, C.-S., 2011. A nonlinear method of removing harmonic noise in geophysical data. Nonlinear Processes in Geophysics 18, 1–13. Jeng, Y., Lin, M.-J., Chen, C.-S., 2004. A very low frequency-electromagnetic study of the geo-environmental hazardous areas in Taiwan. Environmental Geology 46, 784–795. Jeng, Y., Lin, M.-J., Chen, C.-S., Wang, Y.-H., 2007. Noise reduction and data recovery for a very low frequency electromagnetic survey using the nonlinear decomposition method. Geophysics 72, F223–F235. Jolivet, L., Huchon, P., Rangin, C., 1990. Tectonic setting of western Pacific marginal basins. Tectonophysics 160, 23–47. Karous, M.R., Hjelt, S.E., 1983. Linear filtering of VLF dip-angle measurements. Geophysical Prospecting 31, 782–794. Lan, C.Y., Wang Lee, C.M., 1992. The mineral chemistry of Fanpaochienshan gneiss and associated amphibolites, northeastern Taiwan. Journal of Geological Society of China 35, 45–76. Lin, M.-J., Jeng, Y., 2010. Application of the VLF-EM method with EEMD to the study of a mud volcano in southern Taiwan. Geomorphology 119, 97–110. http://dx.doi.org/ 10.1016/j.geomorph.2010.02.021. Loke, M.H., Barker, R.D., 1996. Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method. Geophysical Prospecting 44, 131–152. Monteiro Santos, F.A., Mateus, A., Figueiras, J., Gonçalves, M.A., 2006. Mapping groundwater contamination around a landfill facility using the VLF-EM method — a case study. Journal of Applied Geophysics 60, 115–125. Paál, G., 1965. Ore prospecting based on VLF-radio signals. Geoexploration 3, 139–147. Paterson, N.R., Ronka, V., 1971. Five years of surveying with the very low frequency electromagnetic method. Geoexploration 9, 7–26. Pirttijärvi, M., 2004. Karous–Hjelt and Fraser Filtering of VLF Measurements: Manual of the KHFFILT Program. Department of Geosciences, University of Oulu, Finland. Qiu, Y.-P., Wu, S.-M., Cao, J.-X., Lin, Y.-L., Liu, Y.-C., Chen, J.-A., 2011. Petrograph study of the Ho-Ping granite, eastern Taiwan. Proceedings of the 2011 Joint convention of CGS and GST, p. 451 (in Chinese). Telford, W.M., Geldart, L.P., Sheriff, R.E., 1990. Applied Geophysics, second edition. Cambridge University Press, Cambridge, UK. Teng, L.S., 1996. Extensional collapse of the northern Taiwan mountain belt. Geology 24, 949–952. Van Dam, R.L., 2012. Landform characterization using geophysics—Recent advances, applications, and emerging tools. Geomorphology 137, 57–73. http://dx.doi.org/ 10.1016/j.geomorph.2010.09.005. Wu, Z., Huang, N.E., 2009. Ensemble empirical mode decomposition: a noise-assisted data analysis method. Advances in Adaptive Data Analysis 1, 1–41. Yen, T.P., 1963. The metamorphic belts within the Tananao Schist terrain of Taiwan. Proceedings of the Geological Society of China 6, 72–74. Yen, T.P., 1981. Review of the metamorphic geology of the Central Range of Taiwan. Bulletin of Geophysics National Central University 21, 61–67. Yui, T.F., Okamoto, K., Usuki, T., Lan, C.Y., Chu, H.T., Liou, J.G., 2009. Late Triassic-Late Cretaceous accretion/subduction in Taiwan region along the east margin of South China – evidence from zircon SHRIMP dating. International Geology Review 51, 304–328.