# Clarifying geological structure for coal and marsh gas development using magnetotelluric method.

1. INTRODUCTIONTo clarify the geological structure in an area of coal mining using the resistivity distribution, magnetotelluric (MT) and audio-magnetotelluric (AMT) geophysical prospecting methods are applied in deep underground surveying (Cagniard, 1953). The methods were developed to detect underground resources. The AMT method is more useful than the MT method in the case of shallow underground structures. In this study, resistivity distributions are estimated using both MT and AMT methods. In addition, the method of analyzing MT and AMT data is improved to obtain more accurate results. The study area is southern Kushiro in Hokkaido, Japan. The area has a coal mine for which the recoverable coal reserves are estimated to be about 120 000 000 tons. The coal mine is the only underground casting mine in Japan and produces at least 550,000 tons annually, of which at least one-third goes to the domestic market. The Harutori and Osotsunai faults pass through the area of the coal mine. These faults are related to the distribution of coal and marsh gas. However, the faults have not yet been thoroughly investigated. In the present study, MT and AMT measurements are made on both sides of the faults to clarify the fault structures. The thickness of the Cretaceous layer, which can be a reservoir of marsh gas, is also unknown in this area. One of the objectives of the study is to estimate the thickness of this layer.

2. MT METHOD

2.1. OUTLINE

The MT method employs naturally occurring fluctuations in the Earth's magnetic and electrical fields. The magnetic field fluctuations are due to processes occurring in the Earth's interior. Figure 1 is a schematic view of the generation of electromagnetic waves. The MT method uses ionosphere fluctuations due to the solar wind (about 0.001-5 Hz) and thunder electric discharge (about 1-30 kHz). The frequency range of the MT method is from [10.sup.-4] to [10.sup.3] Hz. The MT method has greater sounding depth than other geophysical prospecting methods because of its broad frequency range. However, it is difficult to judge the signal and noise in MT data. To resolve this problem, a remote reference (Gamble et al., 1979; Takakura et al., 1994) is applied. This method is based on the idea that the local difference in the magnetic field is small and that noise signals in a specific study area are not correlated with noise signals in a remote area. A reference site is thus established separately from the measurement sites. Using magnetic field data at the reference site as a reference signal, data with only slight contamination can be obtained because noncorrelating data between measurement sight and reference sight can be removed as noise data. In this study, the reference site was established 17 km from the measurement sites in a forest without artefacts.

[FIGURE 1 OMITTED]

The MT method cannot clarify details of the shallow underground structure because of its low-frequency range. Therefore, the AMT method is used with the MT method because of its frequency range (i.e., from 101 to 104 Hz). The intensity of the source used in the measurements is thus higher than that in the MT measurements. The AMT method is suitable for estimating a shallow structure.

In the present work, electric and magnetic fields were measured for 15 hours overnight from 18:00 (UT 9:00) to 9:00 (UT 0:00), because the noise effect is slight at this time. The data were amplified and filtered. The measurement period at one site was set to at least three days to obtain data with a high signal-to-noise ratio. In addition, data for which the phase was close to 0[degrees] or 180[degrees] were excluded from the inversion analysis because such data are considered to be affected strongly by artificial noise. Data for which the phase was around 45[degrees] or 135[degrees] were optimal.

2.2. PRINCIPLE

MT and AMT methods are electromagnetic prospecting methods that are based on the assumption of plane electromagnetic waves and a linear relationship between the electric field E and magnetic field H. The complex vectors E and H are related by impedance tensor Z:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (1)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (2)

where the x, y, and z axes are respectively defined along the north-south, east-west, and vertical directions. In the case of a two- or three-dimensional underground structure, variable impedances can be calculated by coordinate rotation. In the case that the structure is two-dimensional, [Z.sub.xx] = [Z.sub.yy] = 0 and [Z.sub.xy] 4 [Z.sub.yx] ([not equal to] 0), which do not depend on the coordinate. [Z.sub.xy] and [Z.sub.yx] are the impedances of the transverse electric (TE; E-polarization) and transverse magnetic (TM; H-polarization) modes respectively. That is, the TE mode can be estimated from the electric field along a strike of the underground structure and the magnetic field at right angles to it. The case for the TM mode is simply the opposite of that for the TE mode. Because the real underground structure is three-dimensional, it is apparent that one-dimensional analysis cannot accurately calculate the resistivity structure, but one-dimensional analysis using the TE mode is known to provide a relatively good approximation (Wannamaker et al., 1984). In this study, the impedances are rotated through 34[degrees] to accord with the strike of the Harutori fault.

Resistivity [rho] is defined by

[rho] = (k/f) x [[absolute value of E].sup.2]/[[absolute value of H].sup.2], (3)

where f and k denote the frequency and a constant, respectively. This equation can only be applied to homogeneous resistivity distributions, whereas real resistivity structures are highly heterogeneous. Therefore, [rho] is apparent resistivity [[rho].sub.a] when obtained by measurement. Therefore, inversion analysis is indispensable in estimating true resistivity distributions (Vozoff, 1986).

The phase of impedance is defined by

[phi] = arg (Z), (4)

[phi] is a parameter related to the resistivity change underground. It takes a value of 45[degrees] for a homogeneous underground structure and a value less (greater) than 45[degrees] in the case that apparent resistivity increases (decreases) with decreasing frequency. The main axis of impedance is the direction in which [[absolute value of [Z.sub.xx]].sup.2] = [[absolute value of [Z.sub.yy]].sup.2] takes a minimum value. Two directions at right angles are thus defined. If the underground is a two-dimensional structure, one of main axes aligns with the strike of the underground structure. Therefore, the main axis of impedance is used to estimate the strike of the underground structure. As a parameter that weakens the electromagnetic field, skin depth 5 is given as

[delta] = 500 x [([[rho].sub.a]/f).sup.1/2], (5)

Here [delta] is the possible exploration depth when using MT and AMT methods. When [delta] and f are respectively 5000 m and 0.3 Hz in (5), [[rho].sub.a] is 30 [ohm] x m. This [[rho].sub.a] value is reasonable in this study. Hence, [[rho].sub.a] data over 0.3 Hz are used in the analysis carried out in this study.

2.3. REMOTE REFERENCE METHOD

MT data are readily affected by artificial noise in the neighborhood. The remote reference method is used to reduce noise in the data. The method is based on the idea that the local difference in the magnetic data is small in a relatively small area and that there is no correlation in noise recorded at a measurement site and a reference site area short distance away (Gamble at al., 1979). Therefore, a remote reference site is established apart from the measurement sites.

[FIGURE 2 OMITTED]

Using magnetic field data recorded at the reference site as a reference signal, measurement data that do not correlate with reference data are removed as noise. Considering impedance Zxx, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (6)

where [H.sub.xr] and [H.sub.yr] are respectively the magnetic fields in the x and y directions at the reference site, < > denotes the cross power, and * indicates the complex conjugate number. A measurement F can be divided into the signal (FS) and noise (FN); F = FS + FN. The mutual power spectrum is then

[F.sub.1][F.sup.*.sub.2] = <F[S.sub.1] x F[S.sub.*.sub.2] + F[N.sub.1] x F[N.sup.*.sub.2]>, (7)

where the subscripts 1 and 2 denote different sites. If F[S.sub.1] and F[S.sub.2] are correlated, and F[N.sub.1] and F[N.sub.2] are not, and F[N.sub.1] x F[N.sup.*.sub.2] is zero. Noise can thus be removed from the measurement data.

3. STUDY AREA AND MEASUREMENT

3.1. STUDY AREA AND SITES

The study area was southern Kushiro in Hokkaido, Japan. The area includes a coalfield containing Japan's only coal pit. The pit accounts for more than one-third of domestic coal production. The coal deposit is estimated to be about 120 000 000 tons. Drilling investigations have been carried out at several locations in the study area, but the thickness of the Cretaceous layer, which may provide marsh gas, is not clear near the Harutori and Osotsunai faults. MT and AMT measurement sites were established on both sides of the faults in the present work. The site north of the faults is denoted S1 and the site south of the faults is denoted S2. Another measurement site S3 was established about 1.5 km from S2 on the southern side of the faults. In addition, a reference site for the MT measurements was established about 30 km from the measurement sites. The AMT method does not need a reference site because it is robust against artificial noise. Measurement and reference sites and estimated fault positions are shown in Figure 2.

3.2. MEASUREMENT EQUIPMENT

MT and AMT measuring instruments, denoted MTU-5U and MTU-5A respectively, were made by Phoenix Corporation (Canada). These instruments obtain data of the two horizontal components of the electric fields ([E.sub.x], [E.sub.y]), and the two horizontal components and one vertical component of the magnetic field ([H.sub.x], [H.sub.y], [H.sub.z]). Three coils and four electric sensors are used to measure the magnetic field. Coils are set in north-south, east-west, and vertical directions. Electric sensors are set 50 m to the north, south, east, and west. To prevent damage to the equipment by rain and animals, the main instrument box and battery are packaged in a plastic box. Measurement instruments and the layout are shown in Figures 3 and 4. Measurements were carried out from 22nd August to 2nd September 2011. The sites are shown in Figure 5.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

3.3. CLEAR OF STATIC SHIFT

The low-frequency band of AMT measurements does not continue to the high-frequency band of MT measurements for a static shift. The static shift is induced by the local resistance abnormality of the superficial part of the ground that corresponds to the higher-frequency band of the MT measurements (Uchida and Murakami, 1989). Figure 7 is a conceptual diagram of the static shift. The static shift can be adjusted in AMT measurements because the AMT method uses higher-frequency data. The formula used to correct the static shift in the MT data is (Ben et al., 1988; Alan, 1988):

[[rho].sub.a] = [[rho]'.sub.a][D.sup.-2], (8)

where [[rho]'.sub.a] is the true apparent resistivity (AMT data) and [D.sup.-2] is the static factor. Data for the apparent resistivity without a static shift are shown in Figure 8. As a result, the effect of the static shift on MT apparent resistivity data can be removed using AMT data of the apparent resistivity.

4. INVERSION ANALYSIS

To estimate the structure of real resistivity [rho], measurements of apparent resistivity [[rho].sub.a] have to be applied to inversion analysis. In inversion analysis, the resistivity model is considered the real resistivity structure when the apparent resistivity calculated with a resistivity model is similar to the measured resistivity. In this study, we assumed a horizontally layered structure model of homogeneous media consisting of N layers and a plane electromagnetic wave propagating in the vertical direction from above at an angular frequency [omega]. This model is often used with TE mode impedance in inversion analysis. Under these assumptions, the impedance, Z, which consists of the magnitudes of the electric field, E, and magnetic field, H, can be defined by the recurrence formula

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (9)

where j, [u.sub.[i.bar]], [h.sub.i] and [[sigma].sub.i] are the imaginary unit and magnetic permeability, thickness and conductivity of the -th layer, respectively. This equation is used to determine the apparent resistivity at the surface, [[rho].sup.a.sub.m]

[[rho].sup.m.sub.a] = (1/[omega][mu])[[absolute value of [Z.sub.1]].sup.2], (10)

where [mu] is magnetic permeability in a vacuum.

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

The objective function of the inversion is the root mean square (RMS) error that expresses the difference between [[rho].sup.m.sub.a] calculated using Eq. (10) and [[rho].sub.a] obtained from the TE mode of the measurement

RMS = [1/n] [summation] [square root of [([[rho].sub.a] - [[rho].sup.m.sub.a]/[[rho].sub.a]).sup.2]], (11)

where n is the number of frequencies used in the MT survey. The error is normalized by [[rho].sub.a] at each frequency, because the values of [[rho].sub.a] differ largely with the frequency. The appropriateness of the resistivity model can be judged by the smallness of the RMS error. Even if this RMS error is almost zero, the model may differ from the true structure owing to overfitting. Simulated annealing (SA) (Aarts et al., 1990) was adopted to overcome this problem because it is well known to produce reasonable results that are independent of initial models and escape from local minima of the error (e.g., Mrinal et al., 1993). The details of the method are described by Asaue et al. (2006).

5. DETERMINING THE SMOOTHING TERM USING AKAIKE'S BAYESIAN INFORMATION CRITERION

The inversion method can be classified as one-, two-, or three-dimensional according to the model used. One-dimensional inversion assuming a horizontally layered structure with SA (Aarts and Korst, 1990) is applied to the estimation of the resistivity distribution in the study. In addition, a smoothing term is applied to inversion analysis to avoiding overfitting. A suitable result is obtained by minimizing U.

[FIGURE 8 OMITTED]

U = [[RMS].sup.2] + [[alpha].sup.2][[C x m].sup.2], (12)

where C is a Laplacian filter, m is the resistivity model of the horizontally layered structure, and [alpha] is a weight coefficient for the smoothing term. The Laplacian is applied to the smoothing term. Akaike's Bayesian information criterion (ABIC) is also applied to decide a suitable value of [alpha] (Uchida, 1993). The value [alpha] is the statistical likelihood that the estimated model for the phenomenon expresses the true model. The idea is to judge the quality of the estimated model by maximizing entropy. In addition, it is important that likelihood L from the Bayesian model is applied. L is defined as

L (m|d) = [integral] p (m|d)[pi](m)dm (13)

where p(m|d) and [pi](m) are the probability distribution and prior distribution respectively. In the Bayesian model, model m is expressed by the prior distribution according to transcendental information, and the occurrence probability of the phenomenon is calculated with a simultaneous density function using probability and prior distributions. This is equivalent to minimizing U to maximize Bayesian model L.

The ABIC is defined by

ABIC = N log (2[pi] [U/N]) - log [absolute value of [[alpha].sup.2][C.sup.T]C] + B, (14)

where B is log[absolute value of [Z.sup.T]Z] + N + 2; N is the number of apparent resistivities, which depends on the number of frequencies: and C is a Laplacian filter. When ABIC is minimized, a suitable a of U is obtained.

6. RESULTS AND DISCUSSION

6.1. COMPARISON WITH GEOLOGY

Measurements and calculations of apparent resistivity are compared in Figure 8. The frequency ranges of the MT and AMT methods are approximately [10.sup.4] to [10.sup.2] Hz and [10.sub.2.5] to [10.sup.0.5] Hz respectively. The results are satisfactory since both curves fit the data well.

Resistivity distributions for sites S1 and S2 are shown in Figure 9. Only one drilling investigation was carried out near the measurement sites. Geology and fault positions were estimated from these drilling data. Hence, these fault structures are not well clarified. The resistivity distribution for site S3 and the nearest geological column, about 1500 m from S3, are shown in Figure 10. Water-selected hard rocks are deposited on the surface layer. It is thus estimated that the low-resistivity zone in the shallow part of S1 and S2 corresponds with groundwater. The low-resistivity zone at depths from 100 to 120 m at S1 is estimated to be a coal or clay layer of the Yubetsu coal layer, but there is no low-resistivity zone in the geological column at a depth of 330 m corresponding to the coal layer. It is thus estimated that the Osotsunai fault exists at depths from 200 to 330 m and the coal layer is displaced by the fault. Therefore, the Osotsunai fault lies between the bottom of the Yubetsu layer and Tennei gravel, and it is connected with the Cretaceous layer. The zone deeper than 350 m is estimated to be a Cretaceous land block that is divided into Osotsunai and Harutori faults. In addition, the low-resistivity zone at a depth of -10 m at S2 corresponds to the aquifer under the Oboro gravel. It is thought that the Harutori and Osotsunai faults are on the north side of S2 and the two dips; it is estimated that there is a northern dip of 210 m for the Harutori fault and an eastern dip of 180 m for the Osotsunai fault, which almost coincide with the dip in the resistivity distribution of 400 m. The low-resistivity zone at a depth of -30 m is due to groundwater in the shallow parts of S1 and S2. The high-resistivity zone at a depth of -100 m corresponds to the Harutori coal layer because this layer has high resistivity at S2.

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

6.2. ESTIMATION OF THE THICKNESS OF THE CRETACEOUS LAYER

Resistivity distributions to a depth of 5000 m are shown for all sites in Figure 11. The typical geologic column in Kushiro is shown in Figure 10. The Cretaceous layer consists mainly of sandstone, mudstone, tuff, and lave in order of depth. These geologic materials have comparatively high resistivity. Therefore, layers with resistivity exceeding 60 [ohm] x m correspond to the Cretaceous in Figure 12. According to the results, layers below a depth of 500 m can be regarded as the Cretaceous layer.

The thickness of the Cretaceous layer is estimated as at least 4000 m in this area. However, the bottom of the Cretaceous layer was not detected at any site. This finding is in agreement with previous results indicating that the thickness of the Cretaceous layer exceeds 3000 m. Clarification will require drilling to depths greater than 800 m.

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

7. CONCLUSIONS

MT and AMT methods were applied to clarify fault structures in Kushiro, Hokkaido, Japan and demonstrate the exploration of deposits of coal and marsh gas. One-dimensional inversion with a smoothing term was carried out to improve the precision of analysis. The main results obtained in the study were as follows.

1. MT and AMT methods were implemented in a coal-producing area and inversion analysis combining these methods was applied to estimate the geological structure. In particular, the groundwater level and coal layers were clarified.

2. The depth of the Osotsunai fault was estimated to range from 200 to 330 m, which differs from what is seen in the geological column.

3. The thickness of the Cretaceous layer was estimated to be at least 4000 m because the resistivity distribution at all sites was almost the same at depths from 500 to 5000 m.

ACKNOWLEDGMENTS

The authors wish to express their sincere thanks to Kushiro Coal Mine Co., Ltd. for assistance with measurements. Sincere thanks are extended to two anonymous reviewers for their valuable and constructive comments that helped improve the clarity of the manuscript.

REFERENCES

Aarts, E.H. L. and Korst, J.H.M.: 1990, Simulated annealing and Boltzmann machines. John Wiley & Sons, New York, 267 pp.

Aarts, E.H.L., and Korst, J.H.M.: 1990, Simulated annealing and Boltzmann machines. John Wiley & Sons, New York, 267 pp.

Alan, G.J.: 1988, Static shift of magnetotelluric data and its removal in a sedimentary basin environment, Geophysics, 53, No. 7, 967-978.

Asaue, H., Koike, K., Yoshinaga, T. and Takakura, S.: 2006, Magnetotelluric resistivity modeling for 3D characterizing the geothermal reservoir in the western side of Mt. Aso, SW Japan. Journal of Applied Geophysics, 58, 296-312.

Cagniard, L.: 1953, Basic theory of the magnetotelluric method of geophysical prospecting. Geophysics, 18, 605-635.

Gamble, T.D., Goubau, W.M., and Clarke, J.: 1979, Magnetetellurics with a remote magnetic reference. Geophysics, 44, 53-68.

Mrinal, K.S., Bimalendu, B.B. and Paul, L.S.: 1993, Nonlinear inversion of resistivity sounding data. Geophysics, 58, 59-74.

Sternberg, B.K., Washburne, J.C. and Pellerin, L.:. 1988, Correction for the static shift in magnetotelluric using transient electromagnetic soundings. Geophysics, 53, No. 11, 1459-1468.

Takakura, S., Takeda, M. and Matsuo, K.: 1994, Effects of regional noise on magnetotellurics and their removal by far remote reference method. BUTSURI-TANSA (Geophys. Explor.) 47, 24-35, (in Japanese with English abstract).

Uchida, T. and Murakami, Y.: 1989, Geothermal reservoir and its resistivity structure. BUTSURI-TANSA (Geophys. Explor.) 42, 458-468, (in Japanese with English abstract).

Uchida, T.: 1993, Smooth 2-D inversion for magnetotelluric data based on statistical criterion ABIC. J. Geomag. Geoelectr. 45, 841-858.

Vozoff, K., ed.: 1986, Magnetotelluric method. Geophysics reprint series 5, SEG, Tulsa.

Wannamaker, P.E., Hohmann, G.W. and Ward, S.H.: 1984, Magnetotelluric responses of three-dimensional bodies in layered earths. Geophysics, 49, 1517-1533.

Hisafumi ASAUE (1) *, Masahito SASAHARA (2), Toru YOSHINAGA (1) *, Yuzo OBARA (1), Kagemi UCHIDA (3) and Hiroyuki MATSUMOTO (3)

(1) Kumamoto University, 2-39-1 Kurokami, Chuo-ku, Kumamoto, Japan

(2) Toda Corporation, 1-7-1 Kyoubashi, Chuo-ku, Tokyo, Japan

(3) Kushiro Coal Mine Co., Ltd., 5-2-23 Okitsu, Kushiro, Hokkaido, Japan

* Corresponding author's e-mail: asaue@kumamoto-u.ac.jp

(Received July 2012, accepted May 2013)

Printer friendly Cite/link Email Feedback | |

Author: | Asaue, Hisafumi; Sasahara, Masahito; Yoshinaga, Toru; Obara, Yuzo; Uchida, Kagemi; Matsumoto, Hiroyu |
---|---|

Publication: | Acta Geodynamica et Geromaterialia |

Article Type: | Report |

Geographic Code: | 9JAPA |

Date: | Apr 1, 2013 |

Words: | 3782 |

Previous Article: | Registrations of seismic foci under mining coal seam --the experience of polish mines. |

Next Article: | Permeability and volumetric changes in coal under different test environment. |

Topics: |