Image of the five elements and prediction of the geothermal field based on gravity, magnetic and magnetotelluric data in the PanZ area

There are two problems in the prediction of the geothermal field in the PanZ area: (1) the plane scopes have some debates, and (2) the vertical scopes need to be further ascertained. Faced with these two problems, a complete set of methods was developed and summarized, and the details are as follows: a geothermal field can be divided into five elements, i.e., heat source, fault channel, thermal reservoir, cap rock and water; then, they are interpreted and imaged with the help of gravity, magnetic and magnetotelluric (MT) data; and finally, according to the integrity of five elements and the correlation between them, geothermal fields are predicted. In the PanZ area, (1) the normalized vertical derivative of the total horizontal derivative of the Bouguer gravity anomaly was applied to identify the fault channels; (2) the water was recognized using the joint interpretation results from an integrated geophysical profile with gravity and MT data instead of a single MT result; (3) the cap rock was inverted with the Bouguer gravity anomaly, using the Parker–Oldenburg inversion method, and with the help of the MT anomaly in the integrated geophysical profile, the vertical distribution of the geothermal reservoir was further ascertained; and (4) the intermediate acid magmatic rock with radioactivity, i.e., a heat source, was identified with the residual magnetic anomaly, imaged using the magnetic forward formula of the cuboid. Finally, the two geothermal fields were predicted and outlined using the above methods. A comparison of the distributions of the geothermal gradient and the outlet water temperatures of the drill holes indicated that the predicted results are credible. To better understand the effect of the method of predicting the geo-thermal field, a 3D geological model was constructed from the inverted results using GOCAD software, and the operating mechanism of geothermal system was analyzed based on the migration, storage, heating and insulation of the water element in the other four elements. To determine the reason for the formation of the geothermal field, the geological evolution of four elements was discussed, except the water element.


Introduction
It is important to determine the geophysical characteristics of geothermal fields for geothermal exploration, and hence, geophysical methods have been widely applied in geothermal exploration (Baillieux et al. 2013;Tontini et al. 2016).The Parker-Oldenburg method was applied to invert the Moho depths (the main heat sources of geothermal fields) in the Gonghe Basin (Wang et al. 2021), the Nankai subduction zone and neighboring regions (Li 2011).This method has the advantage of fast calculation, and it is suitable for inverting the depths of the thermal reservoir and cap rock when there are density differences between them and their surrounding rock.The MT method with resistivity imaging is a valid geophysical tool for surveying underground fluid because rock with fluid can show a lower resistivity than the surrounding dry rock (Aizawa et al. 2009;Ogawa et al. 2014;Shimojuku et al. 2014;Wannamaker et al. 2009).Hence, this method has been applied to image geothermal fluids to identify potential geothermal resources in the Kakkonda area, Northeast Japan (Ishizu et al. 2021), the Fang area, Chiangmai city, Thailand (Amatyakul et al. 2016), the Gediz Graben, West Anatolia, Turkey (Erdogan and Candansayar 2017) and the Eburru area, Kenya (Maithya and Fujimitsu 2019).When thermal reservoirs and faults are rich in water, which can result in a lower resistivity than dry rock, this method can be utilized to recognize underground water.Magnetic data were applied to estimate the depth of magmatic intrusions (heat sources) in the Wikki warm Spring region, northeastern Nigeria (Abraham et al. 2015).Aeromagnetic data were analyzed to determine the magnetic anomaly characteristics related to basic and ultrabasic rocks and to invert their depths within the Gongola Basin in the Upper Benue Trough, northeastern Nigeria (Abubakar et al. 2010).A ground-based magnetic survey was performed to investigate the magmatic heat source near the Silti Debre Zeyet Fault Zone northwest of the Aluto-Langano geothermal field (Kebede et al. 2022).In addition, magnetotelluric (MT) data were applied to determine magmatic heat sources in the Rotorua and Waimangu geothermal fields, the Taupo Volcanic Zone, New Zealand (Heise et al. 2016), the hydrothermal activity in Long Valley Caldera, California (Hermance et al. 1984) and the Kangding geothermal system along the Xianshuihe fault zone of the eastern Himalayas, eastern Tibetan Plateau (Cheng et al. 2022).In general, intermediate acid magmatic rocks (IAMRs), which play an important role in the formation of geothermal fields, have the characteristics of high magnetic susceptibility and high resistivity, so these two methods can be applied to identify the IAMR and invert its depth.In addition, gravity and magnetic data were jointly applied to invert the Los Humeros geothermal reservoir in Mexico using a correspondence map methodology (Carrillo et al. 2022).Gravity and magnetic data were utilized to calculate the Moho structure and Curie point depth relative to the geothermal gradient in the Sulawesi Sea and periphery of Sulawesi Island (Zhang et al. 2021).Aeromagnetic data and residual gravity data were utilized to prospect the geothermal resources in the Sabalan region in northwestern Iran by constituting Curie point depth, geothermal gradient and heat-flow map and evaluating the overburden thickness (Afshar et al 2017).Furthermore, these data and MT data were applied to identify heat sources, hydrothermal reservoirs, and potential geothermal fluid pathways to investigate the Sabalan geothermal area (Afshar et al. 2023).A joint inversion approach of ambient noise surface waves and gravity, which was constrained by the 3D electrical resistivity distribution, was developed to analyze the potential of geothermal resources in the French Massif Central (Ars et al. 2019).Magnetometry and MT surveys were used to image cap rock, topsoil, bedrock, reservoirs and magmatic rocks to confirm a plausible geothermal reservoir in the Abgarm area, Mahallat, Iran (Hosseini et al. 2021).These technical routes of joint interpretations can be applied to enhance the credibility of the interpretive data.
The PanZ area is located in northern Tianjin city, China (Fig. 1a), in the north-central North China Plain; it has a low and flat terrain, which varies from − 5 to 10 m in elevation.The regional geological structure of the PanZ area is defined by the Bohai Bay rift basin, which is a typical multicycle and multistage basin.The tectonic evolution mainly experienced the formation of a crystalline basement and the development of sedimentary cap rocks (Zhao et al. 2015).According to the classification of the Chinese stratigraphic region (Yao et al. 2016), the area belongs to the North China Plain stratigraphic subregion of the Shanxi-Hebei-Shandong-Henan stratigraphic region of the North China stratigraphic province, spanning three level 3 geological structural units (Fig. 1a), i.e., the Jizhong depression, Cangxian uplift and Huanghua depression, lacking the Upper Proterozoic Sinian, Paleozoic Silurian, and Devonian strata.The land surface is buried by Quaternary cover in the study area, and Neogene strata are widely developed.Paleogene and Mesozoic strata are distributed in the Jizhong depression and Huanghua depression, and Paleozoic and Proterozoic strata underlie the above strata in the whole area; their rock types are shown in Table 2.The faults are mainly oriented along the NE, NNE and NW directions, and secondary fractures develop.They are the boundaries between the geological structural units (Fig. 1b).
The area is rich in geothermal resources, and the middle Proterozoic Jixian System, which is the main thermal reservoir in the bedrock area, is widely distributed (Wang et al. 2020).Owing to the destruction of the North China Craton (NCC), the Western Pacific Plate thrusted under the North China Plate at a high angle in the Mesozoic era; then, the former melted to bring about magmatic underplating to the latter, replacement between them and delamination of the latter, and finally, the crust became thinner in North China (Chen 2010;Zhu et al. 2011Zhu et al. , 2012)).This shortened the distance between Fig. 1 a Map showing the regional location and tectonic divisions of the PanZ area in the Bohai Bay Basin (the blue dashed line marks the coastline) (Chang et al. 2018;Tang et al. 2019;Zuo et al. 2017), b geological map of suboutcrops of uncovered tertiary strata in the PanZ area and c plane scopes of the geothermal field defined using the geothermal gradient (data from Lin 2006) and the MT resistivity from the Paleozoic strata in the PanZ area (data from Yang et al. 2009 andLi et al., 2010) the upper mantle and the subbasement of the uncovered Neogene strata, and as a result, there is a high geothermal gradient generally in the area.In addition, the thermal conductivity of each stratum was different (bedrock (LPP strata) > ambient rock (upper Paleozoic and Mesozoic strata) > cap rock (Cenozoic strata)), and as a result, the geothermal gradient in the uplifts was greater than that in the depressions (Xiong and Zhang 1988).Hence, the raised areas of the bedrock in the Cangxian uplift are considered the primary targets for geothermal exploration and development.At the same time, IAMRs with radioactivity (Zhang et al. 2014), Quaternary and Neogene strata with the function of thermal insulation and faults with the function of water migration (Wang et al. 2020) are widely developed in the area.These conditions provided good conditions for the formation of the geothermal field.
Geothermal development and utilization can be traced back to the 1930s in the PanZ area.In the 1970s, under the guidance of Professor Li Siguang, large-scale geothermal resource surveys were carried out.Zhang et al. (2019) and Jiang and Zhang (2012) inverted the depth of the Mohorovicic discontinuity using teleseismic records and gravity data, respectively, which represent the depths of the main heat sources of geothermal fields, i.e., the upper mantle, but the IAMR, i.e., an auxiliary heat source that plays an important role in the formation of the geothermal field, has not been recognized and interpreted.The depths of the strata related to the geothermal reservoirs and the thermal insulation cap rocks were locally exposed and inverted using drill hole data and geophysical profiles (Chen 2020;Sui et al 2019;Jia 2014, Li et al. 2010;Yang et al 2009), while those in the whole area should be further interpreted.In addition, the faults, which were identified using geophysical data, were also debated (Jiang et al. 2010;Zheng et al. 2018), and the underground water has not been interpreted or inverted.In this case, the sensitivity analysis of the main factors of heat transfer and stress in the area was not completely objective (Chen 2020), which is not conducive to the sustainable development and utilization of the geothermal fields.For the above reasons, the plane scopes of the geothermal fields, which are defined using different methods, have some debates (Lin 2006;Yang et al. 2009;Li et al. 2010) (Fig. 1c), and the vertical scopes need to be further ascertained.Hence, it is necessary to study a method to predict the geothermal field in the PanZ area, which provides basic data for the sustainable development and utilization of thermal resources.
Faced with the existing problem of predicting the geothermal field, and combined with the above geophysical methods, all the elements related to the geothermal field were interpreted and connected together to predict the geothermal fields in the PanZ area, and a complete set of methods for predicting geothermal fields was summarized.First, a geothermal field was divided into five elements in this study, i.e., heat source, including the upper mantle and IAMR, fault channel, thermal reservoir, cap rock and water.Then, with the help of gravity, magnetic and MT data, the five elements were interpreted and imaged in the PanZ area.(1) The normalized vertical derivative of the total horizontal derivative (NVDR-THDR), which has the functions of both edge detection and enhancement techniques (Wang et al. 2009), was applied to process the Bouguer gravity anomaly to identify the NE-, NW-and EW-trending fault channels.(2) Water was recognized in the fault and geothermal reservoir using the joint interpretation results from an integrated geophysical profile with gravity and MT data.The former was processed using an inversion method of the 2.5D prism model (Zeng 2005), and the latter was processed using a nonlinear conjugate gradient method (Rodi and Mackie 2001;Zhang et al. 2017).
(3) The interface between the cap rock and the lower Paleozoic and previous strata in age (LPP strata), including the geothermal reservoir, was inverted to acquire the thickness of the cap rock and the depth of the LPP strata with the Bouguer gravity anomaly, which was processed using a cutting method (Li et al. 2020), the iteration method for downward continuation of the potential field (PFDCIM) (Xu 2007) and the Parker-Oldenburg inversion method (Oldenburg 1974;Parker 1973).With the help of the MT anomaly, the vertical distribution of the geothermal reservoir was further ascertained.(4) The IAMR heat source was identified with the residual magnetic anomaly and imaged using the magnetic forward formula of the cuboid (Kuang et al. 2016), and it was verified by the MT resistivity anomaly in a profile.Finally, two geothermal fields were predicted and outlined according to the residual gravity and resistivity anomaly features.Compared with the distributions of the geothermal gradient and the outlet water temperatures of the drill holes in the study area, it is indicated that the two predicted geothermal fields are credible.To better understand the effect of the method of predicting the geothermal field, a 3D geological model was constructed from the above inverted results using GOCAD software (Zhang et al. 2022), and the operating mechanism of geothermal system was analyzed based on the migration, storage, heating and insulation of the water element in the other four elements.To determine the reason for the formation of the geothermal field, the geological evolution of the four elements was discussed, except the water element.

Geophysical parameter statistics
In this study, the Quaternary samples were obtained from manual excavation, the sampled rocks of the Cenozoic and Mesozoic strata were drill cores from the down holes in the PanZ area, and those of Paleozoic, Archaean strata and magmatic rocks crop out in the Jixian Mountains near the PanZ area.The density, magnetic and apparent resistivity parameters are measured using a densimeter, proton precession magnetometer and induced polarization instrument, respectively.The statistical results, rock types and numbers are listed in Table 1.The following conclusions can be drawn from Table 1: The density of strata generally presents an increasing trend with geological age from new to old, and the apparent resistivity shows relatively low, low, medium and high resistances.However, the magnetic susceptibility of the Archaean strata is high, while that of the other strata is low.In addition, acidic, intermediate, basic and ultrabasic rocks (the first three are magmatic rocks) are characterized by high susceptibility and high resistance, and their contents of radioactive elements increase significantly with increasing SiO 2 content (Abd El-Naby et al. 2021;Nasr 2021).According to the statistical results of density and resistivity parameters in PanZ and its adjacent areas (Table 1), five density layers are divided, i.e., Quaternary, Neogene Minghuazhen-Paleogene Dongying (NM-PD), Paleogene Shahejie-Cretaceous (PS-C), Jurassic-Carboniferous (J-C), LPP strata, and four resistivity layers are divided, i.e., Quaternary-Neogene Minghuazhen, Neogene Guantao-Triassic, Upper Paleozoic, LPP strata.The density and apparent resistivity of each layer are shown in Table 1.The LPP strata and their overlying strata exhibit differences in both apparent resistivity and density.Hence, the interface between them, which is fitted using the gravity profile data, can be verified by the MT profile with a high vertical resolution.The fluctuations of the interfaces between the density layers generate a

Geophysical field characteristics
The size of the PanZ area is 34 km × 57 km.The gravity and MT data are measured using a CG-5 gravimeter and V5-2000 MT instrument on the ground, respectively, and their station and line spacings are both 500 m.The grid cell size of the Bouguer dataset is also 500 m, and the Bouguer reduction density is 2.67 g/cm 3 .The magnetic data are from an aeromagnetic survey with a scale of 1:200,000.
The anomaly features of the three kinds of data are shown in Fig. 2. As shown in Fig. 2a, the Bouguer gravity anomaly is arranged alternately with high and low values along the SE direction, and its region can be macroscopically divided into the northwestern area with a low gravity anomaly, the middle area with a high gravity anomaly and the southeastern area with a low gravity anomaly.This is the result of overall density differences among geologic bodies in different areas.The Bouguer gravity anomaly zones in the three areas are distributed along the NE direction, which depends on the shape of the structural basement in each area.The gradient zones between the low-value and high-value areas reflect the main controlling faults between the structural units, and their magnitudes are determined by the fault scales.As seen from Fig. 2b, the magnetic anomaly overall is low, but there are three local high-value areas in the SW, central and NE regions.These high anomalies are mainly generated by intrusive magmatic rocks.As shown in Fig. 2c, the resistivity anomaly from the MT data is similar to the Bouguer gravity anomaly.However, some local high anomalies appear in the high-value areas.These are related to the local bulge of the basement in the geological structural unit.

Constructing the initial model of the geothermal system in the PanZ area
In this study, the geothermal field consists of five elements: heat sources, fault channels, thermal reservoirs, thermal insulation cap rocks and water.The heat sources are the upper mantle and the IAMR.The interface between the crust and the upper mantle, i.e., the Mohorovicic discontinuity, has shallow burial depth in the study area due to the destruction of the NCC in the Mesozoic (Wu et al. 2019), which is less than 35 km (Gao and Li 2014), and there is a regional mirror-image relationship between the fluctuations of the Mohorovicic discontinuity and the bedrock top interface (Zhang et al. 2020b).The shallow burial depth shortens the heat transfer distance from the upper mantle to the land surface.Hence, the heat from the upper mantle, the primary heat source, is easily transferred to the land surface.The IAMR may be associated with radioactivity, so its radioactive decay can provide an auxiliary heat source for the geothermal field.The fault channels in the strata provide channels for the migration of water in the geothermal field.The bedrock-fractured geothermal reservoirs in the area include Mesoproterozoic Jixian, Paleozoic Cambrian and Ordovician strata, and the Jixian reservoir is the most stable in terms of stratigraphic distribution and fracture water content.The thermal insulation cap rocks are the sedimentary strata overlying the thermal reservoirs, including Quaternary and Neogene strata.They can play an important role in the thermal insulation of the geothermal fields.Water is a key element in the geothermal field and can connect the other four elements together.In addition, the heat tends to transfer to the raised portion of the bedrock strata due to the higher thermal conductivity, so the water from the geothermal reservoirs in the bulge or fault block of the uplift is easier to heat.Water is easily extracted because of its shallow burial depth characteristics.Hence, this type of reservoir is considered a target for geothermal exploration and development.
According to the above conditions, the initial model of the geothermal system was constructed in the study area (Fig. 3), and its running processes were illustrated as follows.First, the water supply in the thermal reservoirs comes from meteoric water in the Jixian Mountains in the north.Then, the water flows into the thermal reservoirs in the PanZ area through the faults and fractures in the strata.Next, the water from the geothermal reservoirs in the bulge or fault block of the uplift is heated by heat flow from the upper mantle and the IAMR.Due to the thermal insulation of the cap rocks, the water temperature can be maintained.In this case, a geothermal field can be formed.

Fault identification
The normalized vertical derivative of the total horizontal derivative (NVDR-THDR) is a technique for identifying faults using the maximum and dislocation positions of the linear signal of the potential gradient (Wang et al. 2009).The method strengthens the effective information and suppresses the interference information, which can identify the fault structure effectively.The formulas are expressed as follows: (1) Calculate the THDR of the gravity field: where △g(x,y) represents the gravity anomaly.
The Bouguer gravity anomaly was processed via the NVDR-THDR method in the PanZ area, and then, referring to the geological data in Fig. 1b, its NVDR-THDR was used to identify the NE-, NNE-, NW-and nearly EW-trending fault structures in the study area (Fig. 4), including two level 3 faults (F1 and F2), eight level 4 faults (F3-F9), and three level 5 faults (F10-F12), whose names are listed in Table 2.These faults control the formation and evolution of every tectonic unit.The level 3 faults are the boundaries of the level 3 tectonic units and are distributed in the gradient zones between the areas with high and low gravity anomalies.The level 4 faults are the boundaries of the level 4 tectonic units and have important controls on the formation and evolution of the level 4 tectonic units.The level 5 faults are mainly along NW or nearly EW and are formed by the dislocations of the level 3 and 4 faults or the fault blocks of local areas.These faults formed during multistage tectonic events, and their strikes were related to geological activity during every event (Jiang and Zhang 2012).These faults provided channels for water migration in the geothermal field.
To invert the density interface in the following page ("Interpretation of the depths of the geothermal reservoir and the thickness of the thermal insulation cap rock"), according to the fault distribution (Fig. 4), Bouguer gravity anomaly (Fig. 2a) and geological conditions (Fig. 2a and b), the geological tectonic unit divisions were divided (the details are given in Fig. 4 and Table 3).

Recognition of the water in faults and geothermal reservoirs
According to the statistical density results of the sampled rocks (Table 1), there are five density layers and four resistivity layers in the PanZ area, and there is a common interface between the two kinds of layers, i.e., the interface between the LPP strata and their overlying strata, named the LPP top interface.Hence, the interface, which is inverted using gravity data, can be verified by the MT result.In the study area, an integrated geophysical profile with gravity and MT data was obtained (its position is shown in Fig. 1b).
With the help of the identified faults and the geological structure (Fig. 4), the gravity data were used to invert the density layers using an inversion method of the 2.5D prism model, which is expressed as follows (Zeng 2005): where G is the gravitational constant, σ is the residual density, i is the angular point label of the prism, N is the number of edges of the prism, and where x i and z i are the coordinates, as shown in Fig. 5.
The MT data were applied to invert the electrical structure with the depth of 10 km using a nonlinear conjugate gradient method, whose main formula is expressed as follows (Zhang et al. 2017): where ϕ(m, d) is the inverted objective function, d is an N × 1 observed data vector, m is an M × 1 resistivity parameter vector, f(m) is the forward representation of the resistivity parameter vector (m), C −1 d is the covariance matrix of the observed and forward data, m 0 is the prior information of the model, and C −1 m is the covariance matrix of the model parameters.
Both inverted results are shown in Fig. 6.The error between the measured and fitted gravity curves is small in Fig. 6a, and additionally, the error between the interface depths, which were inverted using gravity data and exposed by drill holes in Fig. 6b, is also small.It is illustrated that the inversion result is reliable.Figure 6b shows that (1) the profile crosses three level 3 tectonic units, the Jizhong depression, Cangxian uplift and Huanghua depression, and five level 4 tectonic units, the Wuqing sag, Yangcun slope, Nanwangping bulge, Chituqiao fault sag and Beitang sag, from NW to SE. (2) The variation in the Bouguer gravity curve reflects the fluctuation in the LPP top interface to a certain extent.(3) Combined with the geological map of the bedrock of the uncovered tertiary stratum (Fig. 1b), Mesozoic strata are distributed in only the two depressions and in the Dabaizhuang subsag, upper Paleozoic strata appear in a small part of the Cangxian uplift, and Quaternary, Paleogene, and LPP strata are distributed in all tectonic units.By combining the density layer shown in Fig. 6b with the geological map of the bedrock (Fig. 1b), the density layer distribution in the whole study area can be concluded, which provides fundamental data for the depth inversion of the LPP top interface in the following section.
The LPP top interface and faults, which were interpreted using gravity data (Fig. 6b), were put in the electrical structure, which was inverted using MT data (Fig. 6c).As shown in Fig. 6c, (1) the LPP top interface is mainly located in the gradient zone (5 between the low-value and high-value areas of the apparent resistivity.This further confirms the reliability of the inverted result from the gravity data.(2) In the positions of the Yangliuqing (F1), Cangdong (F2) and F10 faults, the apparent resistivities exhibit obvious low values.These faults are considered to be filled with water.(3) There is an inward homotropous anomaly of the apparent resistivity in the region from 39 to 46 km, i.e., a relatively low anomaly in the background with high resistivity values, and the anomaly is located in the LPP strata, including the geothermal reservoir.It is identified that the anomaly is generated by the water in the geothermal reservoir.According to the electrical structure shown in Fig. 6c, the water in the faults and geothermal reservoir can be recognized.

Interpretation of the depths of the geothermal reservoir and the thickness of the thermal insulation cap rock
The Jixian Formation, the target geothermal reservoir in the study area, is distributed in the LPP strata, whose overlying sedimentary strata are the thermal insulation cap rock.Hence, the depth of the LPP top interface represents the approximate depths of the geothermal reservoir and the thickness of the thermal insulation cap rock.
In the study area, the density layer distribution is different in each level 4 tectonic unit, i.e., the density difference between the LPP strata and its overlying strata is different.According to the density layer distribution of the profile in Fig. 6, the geological map of the bedrock in Fig. 1b and the density statistical results of the sampled rocks in Table 1, the overlying strata of the LPP strata and the density differences are listed in Table 3 in each tectonic unit.Then, the gravity anomaly generated by the fluctuation of the LPP top interface was extracted using the cutting method (Li et al. 2020), it was extended downward to a lower plane using the PFDCIM (Xu 2007), and the LPP top interface depth in each tectonic unit was inverted using the Parker-Oldenburg inversion method (Oldenburg 1974;Parker 1973).The detailed inversion steps were as follows: (1) According to the division of the level 4 tectonic units, the PanZ area was divided into ten parts.(2) The initial LPP top interface depth in each part was determined by the drill hole and the inverted interface depth in Fig. 6b.(3) The gravity anomaly generated by the LPP top interface fluctuation was separated from the Bouguer gravity anomaly using the cutting method (Li et al. 2020) in each tectonic unit.The method is described as follows: where α is the weighting coefficient related to the half second-order difference value, r is the cutting radius equal to the cutting depth, (x, y) are the plane coordinates, g 0 (x, y) is the measured gravity anomaly, and g n (x, y) represents the nth separated regional anomaly. (7) +g n x − r, y + g n x, y + r + g n x, y − r (4) To reduce the oscillation of high-frequency signals, which is caused by the downward continuation factor of the Parker-Oldenburg formula, the separated gravity anomaly is extended downward to a reference plane whose depth is smaller than that of the LPP top interface using the PFDCIM, which is described as follows (Xu 2007): where FFT() is the fast Fourier transform, IDFT() represents the inverse Fourier transform, h is the continuation depth, ∆g 0 (x, y) is the downward-continued gravity anomaly, ∆g h (x, y) is the upward-continued gravity anomaly, and k x and k y are x +k 2 y h FFT g 0 x, y  the wavenumbers in the x and y directions, respectively.Through the repeated comparison of ∆g h (x, y) and the separated gravity anomaly and the correction of ∆g 0 (x, y), the downward continuation can be completely carried out.
(5) With the help of the downward-continued gravity anomaly and the density difference between the LPP strata and its overlying strata, the Parker-Oldenburg inversion method was applied to calculate the LPP top interface depth in each part, which is expressed as follows: where σ is the density difference, z 0 is the average depth of the interface relative to the reference plane, h(x) denotes the fluctuation depth relative to the average depth, G is the universal gravitational constant, and k = k 2 x + k 2 y , k x and k y are the wave- numbers in the x and y directions, respectively.
The calculated depth plus the continuation depth is the actual inverted depth of the LPP top interface.
(6) The full LPP top interface was obtained by combining the interfaces of every part (Fig. 7a).
Eleven drill holes were distributed in the study area, except that drill hole Da did not reach the depth of the LPP top interface.In addition, the geological profile, which was inverted by gravity data and verified by the MT results with a high vertical resolution (Fig. 6), was applied to make ten virtual drill holes.All these drill hole data were compared with the inverted interface depth (Table 4).As shown in Table 4, the absolute errors between the two depths range from 0.00 km to 0.45 km, and most errors are lower than 0.10 km.The relative errors range from 0.00% to 7.07%, and most errors are smaller than 5.00%.It is illustrated that the inversion accuracy is high and the inverted interface depth is reliable.
To demonstrate the effect of the PFDCIM, the cut gravity anomaly in step (3), which was named the original gravity anomaly, was used to invert the LPP top interface depth using the Parker-Oldenburg inversion method in step (5).Then, the depths in the drill ( 9) hole positions were extracted for comparison with those that were inverted by the extended gravity anomaly and exposed by drill holes (Table 4).As shown in Table 4, after the gravity anomaly is extended using the PFDCIM, the absolute values of the inversion accuracies increase by 0.01 km ~ 0.13 km, and the relative values increase by 0.58 ~ 5.76%.It is demonstrated that the PFDCIM plays an important role in improving inversion accuracy.
To further ascertain the vertical distribution of the geothermal reservoir, with the help of the MT anomaly, i.e., the relatively low anomaly in the background with high resistivity values in the region from 39 to 46 km (Fig. 6c), the geothermal reservoir is inferred to be distributed at depths between 3.65 km and 4.43 km and to have a thickness of 0.78 km.These results are consistent with those of Chen ( 2020) and Sui et al. (2019).

Inversion of geothermal sources
The geothermal flow can be divided into two parts in the study: one part is from the upper mantle, and the other part is from the IAMR.In this area, the depths of the Mohorovicic discontinuity, less than 35 km, have little difference (Fig. 8) (Jiang and Zhang 2012;Zhang et al. 2019).Hence, the difference in geothermal flow depends on the distribution of the IAMR in addition to the geologic structure.According to the statistical results of the physical properties (Table 1), the IAMR has a high magnetic susceptibility and a high resistance.The former was used to identify the IAMR and invert its depth, and the latter was applied to verify the results.The detailed steps were as follows: (1) The residual magnetic anomaly was obtained by subtracting the regional field of the magnetic anomaly from itself (Fig. 9a).As shown in Fig. 9a, there was a high residual magnetic anomaly in the north-central part, which could be generated by the magmatic rocks.
(2) By setting the spatial size, plan position, depth, magnetic inclination, declination and susceptibility of the magmatic rock multiple times, the forward formula for the magnetic anomaly of the cuboid was used to calculate the magnetic anomaly of the magmatic rock, which is expressed as follows (Kuang et al. 2016): where μ 0 is the magnetic permeability in vacuo, M is the magnetic intensity, k 1 , k 2 , k 3 , k 4 , k 5 and k 6 are the coefficients related to the direction cosine of the geomagnetic and total magnetic intensity, r is the distance between one point in the cuboid and one point in the ground, (ξ, η, ζ) is the spatial coordinate of one point in the cuboid, (x 0 , y 0 , z 0 ) is the one of the center of the cuboid, (x, y, z) is the one of one point in the ground, and a, b and c represent the length, width and height of the cuboid, respectively.The reduction to the pole (RTP) anomaly was calculated using an RTP method in the frequency domain, whose factor is expressed as follows (Jing et al. 2017): where θ = atan(v/u), u and v are the circular frequencies in the x and y directions, respectively, I 0 is the geomagnetic inclination, D 0 is the geomagnetic declination, and i = √ −1.
Then, the RTP anomaly was compared with the residual magnetic anomaly (Fig. 8b,  c).When the error between the two anomalies is the smallest, the parameters of the magmatic rock are thought to be actual.The size, depth, magnetic inclination, declination and susceptibility are shown in Table 5, and the plan position is shown in Fig. 9c.As seen from these parameters, the magmatic rock includes two parts.The main part was almost located in the central part of the study area and had a size of 4.5 km × 2.4 km × 3.6 km and a depth of 3.9-7.5 km.According to its magnetic susceptibility, the magmatic rock was considered to be a mixed-melting type magmatic rock, in which the IAMR with strong radioactivity can provide auxiliary geothermal flow for the geothermal field.
(1)(3) To confirm the reliability of the calculated results, an MT profile across the main magmatic rock was inverted to acquire the electrical structure using a nonlinear conjugate gradient method (Fig. 9d).There is a remarkably high resistivity zone in the horizontal position of 17.5-20 km and the vertical position of 4.0-7.4km.
The high resistivity zone is consistent with the calculated position of the mag- matic rock from the residual magnetic anomaly.Compared with the statistical results of the resistivity parameters of the magmatic rocks (Table 1), it can be confirmed that the above result is credible.

Prediction of the geothermal field
LPP strata with the Jixian geothermal reservoir are widely distributed in the PanZ area, and the overlying cap rock is thick enough for thermal insulation (Fig. 7).Based on the Residual gravity anomalies and resistivity anomalies can highlight the anomalies generated by local geological structures, and both can reflect geological bodies of different depths and sizes to a certain extent.This is of great significance in the qualitative inference and interpretation of local geological bodies.As seen from the residual gravity anomaly (Fig. 10a) and the resistivity anomaly (Fig. 2c) in the PanZ area, their local high anomaly reflects the local bulge or fault block of the LPP strata in the Cangxian uplift, and with the help of fault, water and IAMR distributions (Fig. 10a), the two geothermal fields, which were named "A1" and "A2", respectively, were predicted.
To prove the credibility of the two predicted geothermal fields, through comparison of the current geothermal gradient in the study area (Wang et al 2020) and the outlet water temperatures of 5 drill holes in the two predicted geothermal fields, it can be seen that both the A1 and A2 geothermal fields are located in the high-value zones of the geothermal gradient, and all the outlet water temperatures are high enough for the development of the geothermal fields (Fig. 10b).These illustrate that the predicted result is reliable.

3D geological modeling and operating mechanism analysis of geothermal system
To better understand the effect of the comprehensive interpretation method based on gravity, magnetic and MT data for predicting the geothermal field, the interpreted contents in "Imaging of five elements of the geothermal field" were integrated into a 3D model of the geothermal system using GOCAD software.Then, five elements of the geothermal field were clearly shown, and the geothermal system was described in detail.The main steps were as follows: (1) The LPP strata with the Jixian geothermal reservoir, faults, thermal insulation cap rock, upper mantle and IAMR were sorted and converted into a 3D format.The LPP strata, cap rock and upper mantle are all viewed as three layers.To have the 3D model displayed more clearly, the layer between the Mohorovicic discontinuity and LPP strata was omitted, and the depth of the upper mantle under the Mohorovicic discontinuity was reduced by 15 km in the model.(2) According to the spatial distributions of the above elements, the initial model was established, mainly including setting the number of strata, the contact relationships between strata (conformable, eroded, baselap and unconformable) and the fault style (normal or reverse faults).
(3) The spatial range of the model was entered: the plane range adopted the boundary of the study area, and the vertical range was set to 0-20 km.(4) The 3D grid spacing of the fault (plane and vertical spacings are 500 m and 250 m, respectively) and the plane grid spacing of the stratigraphic interface (500 m) are set, and then the faults are established.( 5) The 3D grid spacing of the model (plane and vertical spacings are 500 m and 200 m, respectively) is set and the 3D geological model is built.(6) According to their spatial positions, the IAMR and geothermal reservoir were added to the model for property modeling, while the LPP strata, cap rock and upper mantle were assigned different property values.Thus, the 3D model of the geothermal system (3DM-GSM) was finished (Fig. 11).
To display its internal structure, 3DM-GSM was made partially transparent.As shown in Fig. 11, four elements of the geothermal field, i.e., one upper mantle and two IAMR heat sources, twelve fault channels, two Jixian thermal reservoirs and one sedimentary cap rock, are displayed in the 3DM-GSM.The LPP strata underlie the sedimentary cap rock, which is thick in the Jizhong depression and Huanghua depression and thin in the Cangxian uplift.The IAMR heat sources and the thermal reservoirs are located in the interior of the LPP strata, the faults are distributed in both the LPP strata and the sedimentary cap rock, and the upper mantle is located at the bottom of the 3DM-GSM.The four elements are connected together by the water element, which is distributed in the faults and thermal reservoirs.The connection can be described as follows: First, the water is supplied from meteoric water in the Jixian Mountains with elevations of 318-548 m in the north (Zhang et al. 2020a).Then, the water migrates to the study area through karst fissures and faults at some depth (Brandi et al. 1986) and is heated to a certain temperature by the surrounding rock during migration.Next, the water rises to the shallow part under the action of hydrostatic pressure and flows into the A1 thermal reservoir through the F1, F4 and F10 faults and the A2 thermal reservoir through the F2, F6 and F8 faults.Finally, it is heated to a higher temperature by the heat flow from the upper mantle and the IAMR, and the sedimentary cap rock plays a very good role in thermal insulation to make the water maintain a high temperature.Similarly, the systems of the two geothermal fields run continuously.Through the running process of the geothermal system, the effect of every interpreted element can be seen, which indirectly reflects that the comprehensive method plays an important role in the prediction of the geothermal field in PanZ.

Geological evolutions of the four elements of the geothermal field
To determine the reason for the formation of the geothermal field, the geological evolution of the four elements was discussed, except the water element in the PanZ area.The study area was located in the eastern part of the NCC, which was stable before the Late Triassic (Zhai et al. 2021).In the Late Triassic, the Yangtze plate collided against the NCC along the NNE direction (Liu et al. 2022;Zhao et al. 2018).Hence, a series of NWtrending faults (F7, F8 and F9) and folded structures formed, the study area was uplifted, and the Triassic and Paleozoic strata were denuded.In the Early and Middle Jurassic, the paleo-Pacific plate continuously subducted beneath the eastern part of the NCC toward NW direction (Guo et al. 2022).As a result, the study area was uplifted second, a series of NE-and NNE-trending faults with sinistral strike-slip and thrust properties (F1, F2, F3, F5 and F6) and folded structures formed, and the Upper Jurassic strata suffered denudation.The NW-, NE-and NNE-trending faults were the primary channels of water migration in the geothermal field.These folded structures provide the conditions for the formation of geothermal reservoirs.In the Early Cretaceous, the slab rollback of the paleo-Pacific plate changed the primary stress from NW compression to NNW extension in the eastern NCC, resulting in asthenospheric upwelling, lithospheric stretching and instability.Then, continuous rollback induced continuous upwelling, intense stretching, crustal thinning and fault property transformation from thrust to normal (Li et al. 2022;Quan et al. 2022).Changes in the primary stress dislocated some faults to form new faults (F10, F11 and F12), causing the Jizhong depression to decline and the Cangxian uplift and the Huanghua depression to rise.At the same time, intense magmatic activity occurred, and the magma widely intruded into the shallow sedimentary layers.This was an important thermal event, and the main performances were as follows: crustal thinning provided the conditions for the formation of the primary heat source from the upper mantle in the geothermal field, magmatic intrusion provided the conditions for the formation of the auxiliary heat source from the IAMR, fault dislocation provided the conditions for the formation of the auxiliary channels for water migration, and the decline and rise of the structures provided the conditions for thermal accumulation.In the Late Cretaceous, the paleo-Pacific plate subducted rapidly beneath the NCC at a low angle, leading to regional compression, which resulted in the denudation of part of the Lower Cretaceous strata (Meng et al. 2022).During the Paleocene-Oligocene, changes in the subduction angle and direction of the paleo-Pacific plate caused the study area to stretch, resulting in the subsidence of the Huanghua depression and the formation of local sags and bulges (Su et al. 2014).They all promoted thermal accumulation.After the Miocene, the subduction angle and direction of the paleo-Pacific plate changed again.Tectonic movement was weakened, magmatic activity was also very low, and Neogene and Quaternary strata were deposited (Li et al. 2014).The sedimentary layers provided the cap rock for thermal insulation.On the basis of the evolution, the reason for the formation of the geothermal field is easily understood.

Discussion
In this study, five elements of the geothermal field, i.e., fault channel, water, thermal reservoir, cap rock and heat source, including upper mantle and IAMR, were interpreted and imaged with gravity, magnetic and MT data in the PanZ area.Then, with the help of the residual gravity anomalies and resistivity anomalies, the two geothermal fields were predicted.Compared with the previous extent defined using the geothermal gradient (Lin 2006), the A1 geothermal field is completely within the one, one part of the A2 geothermal field is within the one, and another part is beyond the previous extent.Compared with the previous extent defined using the MT resistivity from Paleozoic strata (Yang et al. 2009;Li et al. 2010), the A2 geothermal field is close to the one, and there is no overlap between the A1 geothermal field and the one.Through comparisons with the updated geothermal gradient figure (Fig. 10(b)), the current extent of the geothermal field is closer to the extent of the high-value zones than the previous two extents.This proves that the current extent is more accurate for geothermal development and utilization.
The faults, which provide channels for water migration in the geothermal field, are identified using the NVDR-THDR of the Bouguer gravity anomaly, and the results are consistent with those of Zheng et al. (2018).In previous studies, groundwater was generally recognized using the MT method.However, the recognized water has ambiguities due to the diffusive electromagnetic fields and nonunique features in the inversion (Chave and Jones 2012;Ogawa and Uchida 1996).To solve this problem, recognition of the water in the faults and geothermal reservoirs has two steps in this study.The first step is that the faults and geothermal reservoirs are interpreted using the gravity profile; the second step is that the water in the faults and geothermal reservoirs is recognized using the relatively low resistivity anomaly of the MT profile in the same position as the gravity profile.This can avoid misjudgments from a single MT result.For example, a geological structure with the characteristic of a low resistivity can also cause a low resistivity anomaly.The reason is that the first step limits the extent of recognizing water to the faults and geothermal reservoirs, whose resistivities are reduced by the water.
The depths of the geothermal reservoir and the thickness of the thermal insulation cap rock are inverted after the gravity anomaly is extended downward to a lower plane using the PFDCIM.Through the comparison between the inversion depth and exposure depth from the drill hole, it is illustrated that the inversion accuracy is high.A comparison of the inversion accuracies with and without the PFDCIM illustrates that the PFDCIM plays an important role in reducing the oscillation of the inversion signals (Fig. 7b) and improving the inversion accuracy.The IAMR with strong radioactivity, which can provide auxiliary geothermal flow for the geothermal field, was first explored and inverted in the study area.The mutual confirmation of the residual magnetic anomaly and the MT resistivity anomaly makes the results credible.In addition, the plane extents of the geothermal fields were predicted using the residual gravity anomalies and MT resistivity anomalies, and the vertical distributions of the geothermal reservoirs were ascertained using the latter.This provides the most direct data for the development and utilization of geothermal fields.3D geological modeling integrates all the interpreted results into a 3D model, which can aid in understanding the effect of the above interpretation methods and analyzing the operating mechanism of geothermal system.A discussion of the geological evolutions of the four elements helps in understanding the reason for the formation of the geothermal field.

Conclusions
This study made full use of gravity, magnetic and MT data, interpreted and imaged the five elements of the geothermal field in the PanZ area, predicted two geothermal fields and described the running processes of the geothermal system.This provided a complete set of methods for predicting the geothermal field.The main conclusions are as follows: (1) The five elements of the geothermal field in the PanZ area were interpreted and imaged with gravity, magnetic and MT data, which provided fundamental data for the prediction of the geothermal field.(2) To avoid misjudging the low resistivity due to the geological structure, the joint interpretation results from an integrated geophysical profile with gravity and MT data instead of a single MT result were applied to recognize the water in the fault and geothermal reservoir.(3) To reduce the oscillation of high-frequency signals, the gravity anomaly is extended downward to a plane whose depth is smaller than that of the LPP top interface using the PFDCIM before being calculated using the Parker-Oldenburg formula.(4) The IAMR heat source, which was identified and imaged with the residual magnetic anomaly, was verified by the MT resistivity anomaly in a profile.(5) On the basis that the spatial distributions of the geothermal reservoir, cap rock, fault, water, IAMR and Mohorovicic discontinuity were comprehensively analyzed, the two geothermal fields were predicted according to the local bulge or fault block structures, which were reflected by the residual gravity and resistivity anomalies in the PanZ area.(6) To better understand the effect of the comprehensive interpretation method for predicting the geothermal field in this study, a 3D geological model of the PanZ area was constructed from the inverted five elements of the geothermal field using GOCAD software, and the operating mechanism of the geothermal system was analyzed based on the migration, storage, heating and insulation of the water element in the other four elements.(7) To determine the reason for the formation of the geothermal field, with the help of the main geological events in the eastern NCC, in which the study area was located, the geological evolution of the four elements was discussed, except the water element.

Fig. 2
Fig. 2 Map of the (a) Bouguer gravity anomaly, (b) magnetic anomaly, and (c) MT resistivity anomaly on a cross-section at a depth of 4 km in the PanZ area

Fig. 3
Fig. 3 Sketch of the supplement and runoff of water in the geothermal system in the PanZ area

Fig. 4
Fig. 4 NVDR-THDR of the Bouguer gravity anomaly and fault distribution in the PanZ area

Fig. 5 Fig. 6
Fig. 5 Schematic diagram of a polygon prism model

Fig. 7
Fig. 7 Inverted LPP top interface depth (a) with and (b) without the PFDCIM in the PanZ area

Fig. 9
Fig. 9 Contour map of (a) residual magnetic anomaly in the PanZ area, (b) high residual magnetic anomaly, (c) fitted magnetic anomaly and (d) MT resistivity anomaly in a profile across the area with high residual magnetic anomaly

Fig. 10
Fig. 10 Map of the inferred geothermal field in the PanZ area, whose backgrounds are (a) the residual gravity anomaly, the distribution of faults, IAMRs and water, and (b) the geothermal gradient (data from Wang et al. 2020)

Fig. 11
Fig. 11 3D model of the geothermal system in the PanZ area (× 2 in vertical coordinates)

Table 1
Statistical results of density, magnetic and resistivity parameters in PanZ and its adjacent areas

Table 1 (
continued) AD average density, AGD average group density, AMS average magnetic susceptibility, AAR average apparent resistivity, AGAR average group apparent resistivity

Table 2
Corresponding table of fault numbers and names in the PanZ area

Table 3
Geological tectonic unit divisions in the PanZ area

Table 4
Comparison table of the LPP top interface depths, which were inverted by gravity data and exposed by drill holes or virtual drill holes DHN drill hole number, ED exposure depth, ID inversion depth, AE absolute error, RE relative error, PFDCIM the iteration method for downward continuation of the potential field

Table 5
Calculated parameter list of the magmatic rocks