Magnetotelluric Imaging of the Zhangzhou Basin Geothermal Zone, Southeastern China

Abstract

The geothermal zone of southeast China, which is one of the country’s known geothermal zones, contains significant natural geothermal resources. To understand the formation of geothermal resources, a magnetotelluric (MT) investigation with a site spacing of 1–2 km was carried out around the Zhangzhou Basin. The recorded MT data were processed by robust time series and remote reference processing techniques. The data analysis results revealed that two-dimensional (2-D) modeling can be used to approximately determine the electrical structure. The joint inversions of TE and TM modes have been performed after distortion decomposition. In the inversion models, a low resistivity cap of 200–800 m thickness was observed, which represented the blanketing sediments composed of Quaternary and volcanic rocks of the late Jurassic period. The presence of high resistivity above a depth of 20 km indicates the granites are widely developed in the upper and middle crust. MT measurements have revealed some deep-seated high conductive zones, which were inferred to be partially melting at depth of 8–17 km, which is likely to be reason behind the formation of higher-temperature hot springs. The results also show that there is a shallower Moho, which indicates that the heat from the upper mantle may have a big contribution to the surface heat flow. Fractures-controlled meteoric fluid circulation is the most likely explanation for the hot springs.

1. Introduction

Six main geothermal zones can be recognized in China (Figure 1), including the Taiwan geothermal area, Southeast coast geothermal zone, Sichuan-Yunnan N-S geothermal area, Tancheng-Lujiang geothermal area, Xizang-Yunnan geothermal area and Qilian-Luliang arc-shaped geothermal area [1]. The southeast coast geothermal zone mainly includes the provinces of Guangdong, Fujian, southern Hunan and eastern Jiangxi. The thermal waters in this zone are mainly low-medium temperature convective geothermal systems. Zhangzhou Basin, which is well known for its numerous hot springs with temperature over 70 °C (Figure 2), represents one of the most extensive hot spring activities in China. The geothermal field of Zhangzhou (red star in Figure 2) has a temperature of 106 °C at the orifice and 121.5 °C at the bottom of the borehole with the depth of 90 m. This is the highest temperature field founded in the southeast coast [2].

It is clear that the occurrence of the hot springs is related to granite and volcanic intrusion, which is controlled by a series of tectonic fractures. To better understand of the thermal activities, numerous geophysical studies have been conducted in Zhangzhou Basin since 1980s. A low-velocity layer is found at the depths of 12–16 km by deep seismic wide-angle reflection and refraction-profile. This layer was inferred to be melting or partial melting, and is likely to be the crustal internal heat source of Zhangzhou geothermal zone [3,4,5]. Seismic tomography shows that the low-velocity anomalies extend to deep depths near Zhangzhou Basin with high convective heat flow [6]. Previous magnetotelluric (MT) suggests that Zhangzhou Basin is located in the uplift of the mantle, and the depth of the asthenosphere is less than 90 km [7]. In contrast to the mantle uplift, a high conductive thin layer was found at depths of 10–13 km [8]. However, the previous research results are mainly regional, there is no detailed study aimed at Zhangzhou Basin. Thus far, there is no unified understanding of the heat sources which mainly focus on three views: the residual heat from magma, the heat generated by radioactive decay of granites and the heat generated by deep molten bodies.

The electrical conductivity of rocks mainly depends on the saline water content and many other interrelated parameters such as temperature, porosity, permeability and mineral alteration [9]. Of all the electromagnetic methods, MT seems to be the most appropriate method for geothermal areas since the conductive anomalies caused by both water-reservoirs and magmatic bodies can be detected by the MT method. Moreover, MT can probe both shallow and deep areas up to several tens of kilometers. Because of the electrical conductivity sensitivity and the investigation depth, MT has been widely used in geothermal explorations [10,11,12]. MT data inversion using geological and geophysical constraints enables us not only to determine the general conductivity distribution in the studied area but also to delineate the geothermal reservoir and detect the cap layer [9]. In this paper, we describe the electrical structures in Zhangzhou Basin geothermal zone with the purpose of showing the subsurface anomalous structures related to geothermal resources. Our results will help to determine the heat source and underlying magmatic system supplying the geothermal fields.

2. Geological Setting

The Zhangzhou Basin geothermal zone, which is located in the southeast coast of China (Figure 2), is geologically special as it was created from the interactions of Eurasia Plate, Philippine Plate and Pacific Plate [7]. Multistage tectonic events have taken place in this region since the Paleozoic time, such as the Caledonian fold movement, the Indosinian orogeny and the Yanshanian volcanic magmatic activity [13]. Of these tectonic events, the Yanshanian movement was the period of peak of the tectonic movement and magmatic activity. This has resulted in extensive Yanshanian granites being widely distributed in the study area, while the covering layer is mainly composed of Upper Jurassic volcanic rocks. The thin Quaternary sediments mainly distribute in the middle of the basin and the sides of rivers. Under the influence of multiple tectonic movements, fault activities have become well developed, forming NE and NW fracture systems under the surface. These fracture systems controlling the migration and storage of groundwater are important structure conditions promoting the formation of hot springs. The hot springs are mainly distributed in the region dominated by Yanshanian magmatic activity, while most of them are located in the cross positions of NE and NW direction faults. Most notably, three major fault zones, Changle-Zhaoan fault zone (strike NE), Nanjing-Xiamen fault zone (strike NW) and Jiulongjiang fault zone (strike NW) converge at these points [14]. These deep faults and accompanying faults play an important role in the formation of geothermal fields.

3. Methods, Data Acquisition and Analysis

MT is a passive source electromagnetic method that uses natural magnetotelluric fields to investigate the resistivity structure of the Earth [15]. The varying electric and magnetic fields are coupled through Maxwell’s equations. Subsurface resistivity can be inferred by simultaneously measuring both fields. The penetration depth of electromagnetic waves depends on the period and resistivity. The short period waves that contain shallow depth information are attenuated rapidly. While the long period waves can probe into the deep earth. This phenomenon of the penetration depth of the MT fields increasing with period and resistivity is well-known and is called the skin effect [16]. Generally, the natural MT field above 8 Hz is generated by electrical storms, while this field is due to solar activity at frequencies below 0.125 Hz. In the MT method, two components of the electrical field (Ex, Ey) and three components of the magnetic field (Hx, Hy, Hz) are measured on the surface of the Earth. The x-axis and the y-axis are normally taken to be in magnetic N-S and magnetic E-W respectively. The horizontal electric (E) and magnetic (H) fields are measured in the time domain and are transferred into the frequency domain using Fourier transformation in order to receive the MT response function.

A total of 157 broadband (320–0.00055 Hz) MT sites have been collected using Phoenix MTU instruments. The horizontal electric (Ex, Ey) and magnetic (Hx, Hy) field components were measured in magnetic N-S and magnetic E-W with a site spacing of 1–2 km at each site. The data were processed from time series to response functions using robust remote reference cascade decimation techniques from the study by Jones et al. [17], as implemented by the Phoenix Geophysics software package SSMT2000. The smoothing process created a generated curve of the apparent resistivity and apparent phase, which was conducted with MTEditor program.

In magnetoetelluric analyses, the response of regional structures can be masked by the galvanic distortions produced by small-scale inhomogeneities that are near the surface [18,19]. In this study, the multi-site multi-frequency distortion decomposition developed by McNeice and Jones [18], which is based on the Goom-Bailey (GB) [20], has been used to estimate the galvanic distortion of the MT response as well as the regional strike and dimensionality of the MT data prior to 2-D inversion. Rose diagrams are used to depict GB strike directions of all sites for whole frequencies (Figure 3).

We obtained the tensor principal direction as being almost N45° E (or N45° W) with 90° ambiguity. Furthermore, the strike directions along profile B and C are more consistent than profile A and D. Since the regional tectonic and geological trend supports a dominant NE direction, we chose N40° E for profile A, N45° E for profile B, N40° E for profile C and N50° E for profile D, respectively. Furthermore, the whole frequencies of the MT data are rotated to the respective strike directions before inversion. The absolute values of the twist and shear values for unconstrained data of all the sites obtained by GB decomposition are illustrated in Figure 4 which indicate the level of distortion. Over 80% of the twist angles were distributed in 0°–15° moderately. And the values of the shear angles show a larger relative variation with over 75% in the range of 5°–25°. The values of the twist and shear angles indicate that the level of distortion is moderate in the study area. In addition, the impedance data were recalculated at a constrained geo-electrical strike direction for each profile. And the distribution of average data RMS values for the whole frequency range from the GB decomposition analysis along the four profiles is shown in Figure 5. The average misfit is considered to be an acceptable value when it is less than 2.0 [21]. Most of the sites can meet the requirements.

According to the above analysis, the deep structural trend of the region has NE orientation. That means that two profiles (B and C) are perpendicular to the main structure, while two profiles (A and D) are parallel to the main structure. Since the MT strike is controlled by the deep basin structure, it is not suitable to use profile A and D to obtain the information at a deeper depth in the 2-D model. Furthermore, the results from GB decomposition also portend that the geo-electric structure of profile A and D is more complex. Therefore, we used only profile B and C to obtain information at a deeper depth with the frequency range 320–0.001 Hz, and limited the use of the two other profiles A and D at depths of less than 5 km with the frequency range 320–0.1 Hz.

4. Two-Dimensional Model

The rotated data were subjected to 2-D inversion analyses, using the non-linear conjugate gradient algorithm developed by Rodi and Mackie [22], which was implemented in data interpretation package WingLink. Before computing the final models, several inversion parameters should be set to fit the data appropriately. TM-mode data is relatively insensitive to the depth extent of a subsurface body [23], and TE-mode data is more sensitive to the three-dimensional effects [24]. Hence, mixed mode analysis was used to help clarify the geologic structures beneath the profile. Static shift, which is caused by local near-surface inhomogeneity, can result in the parallel movement of the apparent resistivity curves, while the phase curves are basically consistent. During the last few decades, many techniques have been developed for static shift correction [25,26,27]. Setting a large error for the apparent resistivity compared with those for the phase during inversion is another effective way to accommodate the static shift [28]. In our models, the error floors for apparent resistivity and phase were set to 10% and 5% for the TM mode and 30% and 10% for the TE mode. Additionally, the initial model was given a homogeneous half space with 100 $\mathsf{\Omega }·\mathrm{m}$ in the inversion of the data.

The study area is located on the southeast coast of China, while the easternmost sounding site is located approximately at a distance of 10–20 km from the East China Sea coastline. The large contrast in electrical conductivity between the seawater and land can severely distort the observed MT responses, especially when the separation distance from the coast is smaller than the skin depth of the frequency [29]. This is generally known as the “sea effect”. One possible approach to overcome this influence is to correct the sea effect by an iterative method [29] or to directly incorporate the sea into the model for inversion [30,31]. As profile D and the SE limits of profiles B and C are very close to the sea, the sea effect must be expected. Therefore, in order to take the sea effect into account, the sea has been inserted in the models of profile B and C during inversion. Because of the limited inversion depth, the influence of sea effect will be small for profile D.

The validity of the inversion models should be estimated before looking at the final models. The quality of the fitting between the observed and the calculated data of the model can be evaluated by the RMS misfit. A low RMS misfit indicates a likely geo-electric model. One typical apparent resistivity and phase curve for each profile is shown in Figure 6. It can be seen that the calculated data fits well with the observed data at all periods. The final MT resistivity models obtained along the four profiles have an average RMS misfit of 2.05, 2.21, 2.12 and 2.38 for profile A, B, C and D, respectively. The RMS values show good consistency between the observed data and the model calculated. Furthermore, these values are considered to be acceptable for the vast majority of the sites along the profiles. The comparisons between the observed and the calculated pseudo-sections for the models are shown in Supplementary Material. The calculated and observed apparent resistivity and phase of both modes fits well, while the fit of TM mode is better than TE mode.

5. Results and Discussion

The MT results have provided a valid electric structure of Zhangzhou Basin geothermal zone. Of these, the results of profiles A and D are mainly used for revealing details of the shallow portion of the study area, while the results of profile B and C are mainly used for revealing the deep portion. For profile B and C, a version of the model restricted to 5 km is provided to understand the shallow details. In Figure 7, a low resistivity layer with a thickness of 200–800 m was inferred along the four profiles. With local and regional geology, we inferred that the low resistivity layer could contain the sediments of Quaternary and volcanic rocks of the late Jurassic. Rybach et al. [32] regarded these low resistivity sediments as blanketing sediments. It can effectively prevent the loss of deep heat and is one of the important conditions encouraging the formation of a geothermal field. Furthermore, the thickness of the low resistivity cap in the southwest and the east of the study area is thicker, especially in the estuary of Jiulongjiang masked in profile D. The results show that the relatively low resistivity extends deeper, which is possibly due to the reservoir.

The deep resistivity structures are shown in Figure 8. The two models present the high resistivity at depths of more than 20 km. One possible interpretation is that the intrusive rocks are widely developed in the deep Earth. It can be noted that the study area is a part of South China Block and is located between the Philippine plate and Eurasia plate (Figure 2). The tectonic activity of the area has been intermittent since the Caledonian stage due to the subduction of the Pacific plate. There are several volcanic magmatic activities, and Yanshanian movement is the period of peak magmatism. This has led to extensive granites of late Yanshanian being distributed widely in the study area (Figure 2). There are several view related to the tectonic background of the formation of magmatic rocks, which mainly include the flat-slab subduction model of the Pacific plate [33], changing subduction angles for J-K magmatism [34], mantle plume model [35] and lower crustal delamination model [36]. Among these views, the model of changing subduction angles, involving the interactions between subduction and underplating, appears to be consistent with most of the factual evidence [34]. Since most of the hot springs is exposed in igneous rocks, and the granites are widely distributed in the study area, the residual heat of Yanshanian intrusive rocks was once considered as the heat source of Zhangzhou Basin [14]. As we know, the role of magmatic activity in the formation of geothermal anomalies is not only related to the geologic time but is also related to the scale and depth of the intrusive rocks and the thickness and thermal insulation performance of the cap rock. Through previous analyses, we know that the cap rock is not very thick, which means that there is no great advantage in thermal insulation performance. In addition, Smith and Shaw once provided an estimate of the amount of the thermal energy that still remains in the magma for the conduction mode, which is referred to as dry rock [37]. The theory proved that the magma, which invaded 100 million years ago, is unlikely to have formed any hot geothermal system. Therefore, the heat source of Zhangzhou Basin is unlikely to be the residual heat of Yanshanian intrusive rocks. However, the hot spring distributions are indeed closely related to Yanshanian granites. In fact, the heat in the earth originates mainly from the heat generated by the decay of radioactive isotopes, which are predominantly ²³²Th, ²³⁸U and ⁴⁰K. The average abundances of Th, U and K for local outcropped granites are 23.6 ppm, 6.6 ppm and 4.28%, respectively, while the average heat production is 3.7 $\mathsf{\mu }\mathrm{W}/{\mathrm{m}}^3$ [38], which is 1.48 times that of the global average granite (2.5 $\mathsf{\mu }\mathrm{W}/{\mathrm{m}}^3$: [39]). The results of seismic tomography showed the presence of a wide range of low-velocity zone range from of 3–20 km, which is the reflection of high geothermal anomalies [6]. This prompts us to believe that the heat generated by the radioactive decay of heat-producing elements in granites contributes to the formation of hot springs.

The deep resistivity structure of profile B shows two low resistivity zones (C1 and C2) present at depth of 8–17 km with resistivity of 30 $\mathsf{\Omega }·\mathrm{m}$. Furthermore, there is an upwelling zone at the end of profile C. The surface geothermal activity in the area is confined within two NE trending fault zones, a NW trending fault zone and an approximate EW trending fault zone, which are namely the Changle-Zhao’an fault zone, Fu’an-Nanjing fault zone, Jiulongjiang fault zone and Nanjing-Xiamen fault zone, respectively. Furthermore, the depth extensions of these fault zones may be controlling the lateral extension of the conductive zones. The upwelling zone at a depth 15 km in profile C may indicate the Changle-Zhao’an fault. Furthermore, the surface positions of C1 and C2 coincide geographically with the major fault zone named the Fu’an-Nanjing fault zone and Changle-Zhao’an fault zone on the geological map. The high-conductive zone in the crust has been considered to play an important role in crustal dynamics, and it is a window to reveal the magmatic activities and earthquakes in the plate. About the formation, many geologists have proposed different opinions, which mainly includes partial melting [40,41], ductile shear zone [42] and brine layer caused by dehydration [43]. Because the high-conductivity zone is complex in both structure and distribution, these views are only established under certain tectonic conditions. However, in volcanic and geothermal areas, the cause of high conductivity at a depth of approximately 10 km may be due to the presence of melting or partial melting materials. The depth of the high-conductive zones is often similar to the low-velocity layer discovered by seismic data [3,4,5]. Xiong et al. considered that three situations can result in the formation of the low-velocity layer: the existence of a crushed zone; the change of pressure, temperature and character of water; and the existence of partial melting materials. Of these situations, they placed more emphasis on the third scenario [4], which was consistent with the ideas of Liao et al. [5]. Hu et al. demonstrated that the conductive layer of the southeast coast may be composed of the partial melting materials caused by the upwelling asthenosphere and the basalt bottom invasion through the bouguer gravity anomaly results and the amplitude of the three-dimensional analytic signal of magnetic anomaly [44]. However, the molten materials cannot be formed in Yanshanian period although it was the period of peak magmatic activity, because the heat in magma has been lost after 100 million years. We think that it is likely that these molten materials were formed in the Cenozoic. In fact, the Cenozoic volcanic belt along the southeast coast of China is an important part of the Pacific tectonic-magmatic belt, and the Cenozoic volcanic rocks are dispersed and widely distributed in the southeast coastal areas. The rock types of the volcanic rocks are mainly tholeiite basalt and alkaline basalt which indicate that the magma originates from the asthenosphere [35]. Therefore, the high-conductive zones are likely to be formed by the upwelling of the asthenosphere mantle along the faults in the extension environment. Usually, high-temperature geothermal systems occur where the magma intrudes into high crustal levels (<10 km) and hydrothermal convection can take place above the intrusive body [9]. The Ar-Ar isotope shows that the volcanic rocks of Niutou Mountain in Zhangzhou were formed 15–17 million years ago [35], and the heat flow caused by molten magma at the depths of 12 km can be reflected in the underground that is at a depth of 3 km in just 95,000 years [45]. Therefore, the high-conductive zones may be the deep heat sources of the hot springs. At least they contribute to raising the temperature of the hot springs.

From the inversion results (Figure 8), we can also see that both models invariantly present an electrical boundary at a depth of approximately 30 km, which we may interpret it as electric Moho. This depth of electric Moho is similar to the depth (29–32 km) found by seismic data [3,4,5]. The shallower Moho may indicate that the heat flow from the upper mantle has a big contribution to the surface heat flow. Heat generated by radiogenic heat production in the crust and heat flow from the upper mantle are two primary components of heat flow measured in stable continental areas [46]. The mean average value of surface heat flow in Zhangzhou is 90 mW/m² which is higher than the global mean heat flow of 65 mW/m² for the continental crust [47]. According to the average heat production of 3.7 $\mathsf{\mu }\mathrm{W}/{\mathrm{m}}^3$, a previous study [38] found that the contribution of heat generated by the heat-producing elements to the surface heat flow was 37–45 mW/m², which is less than the contribution of heat from the upper mantle. In fact, Zhangzhou has a strong crust and mantle mixture in terms of both minerals and energy. The changing subduction angles result in vigorous basaltic underplating and widespread magma mixing between basaltic and felsic underplates. In Zhangzhou granite, there is a significant number of dark inclusions and mafic dike. The ${ϵ}_{ND}\left(T\right)$ walues of granitoids and basalts for the calc-alkaline complexes of Zhangzhou are −4.2 and −3.12 [34]. The negative ${ϵ}_{ND}\left(T\right)$ values of mafic rocks can be generally explained by the mafic rocks originating from enriched mantle wedge and/or the lithospheric mantle [34]. A large amount of mantle upwelling inevitably brings a considerable amount of mantle heat.

Through the geophysical and geochemical analysis evidence provided above, we have reason to believe that the heat source of the basin is mainly composed of three primary components: heat generated by heat-producing elements in the crust, heat from the upper mantle and heat from the partial melting materials in the mid-upper crust. However, not all these three sources have an effect on all the hot springs. In fact, the formation of hot springs is a complicated process. It is not only related to the source of heat, and the blanketing sediments, but also to the source and channel of water and so on. Deuterium-oxygen isotope analyses have shown that the geothermal fluid was derived from the meteoric water [48]. Approximately 90% of the hot springs is obviously related to one or two groups of faults, while 82% of the hot springs are exposed to the magmatic rock area [49]. Therefore, it is clear that the outcrop of the hot springs is related to granite and volcanic intrusion, and is controlled by a series of tectonic fractures. In view of the above mentioned facts, we consider the formation of the hot springs as shown in Figure 9. The cold atmospheric precipitation, which permeated and converged in the faults, became warmed by the bedrock. The heat of the bedrock mainly comes from heat generated by the radioactive decay of granites and upper mantle, forming medium and low-temperature hot springs in the favorable parts due to the effect of pressure difference. In this process, the faults supplied the space and channels for groundwater transportation. In some places, due to the presence of anomalous low-resistivity zones in the middle crust, the temperature of the hot springs will be higher.

6. Conclusions

The MT method has been used to study the geo-electric structure of the Zhangzhou Basin geothermal zone. The 2-D resistivity inversion models of MT data suggest a cap rock composed of Quaternary and volcanic rocks of the late Jurassic period, which are mainly distributed in the west and east of the basin. However, the blanketing sediment layer is not very thick. The reservoir of the study area is a dense fracture network, which consists of host rock and the natural fluids contained in its fractures and pores. The fracture network plays a dominant role in water circulation and aggregation, while the faults provide channels for water movements. In view of the widely-distributed granite and the shallower Moho, the heat source of medium and low temperature hot springs is mainly composed of two parts: the heat generated by granite in the crust and the heat from the upper mantle. In addition, our results also show the presence of anomalous low-resistivity zones in the middle crust which may influence the locations of high-temperature hot springs.