Subsurface temperature prediction by means of the coefficient correction method of the optimal temperature: A case study in the Xiong’an New Area, China

Temperature is one of the key characteristics of the earth’s interior, the knowledge of which determines our ability to study fundamental scientific and applied geothermal problems (Spichak, 2015). Thus, it is crucial to estimate the subsurface temperature accurately.

At present, the common methods of estimating the earth’s interior temperature can be classified into three categories: (1) borehole temperature logs, (2) heat transfer modeling (Björnsson, 2008; Fuchs and Balling, 2016a, 2016b; Poulsen et al., 2017), and (3) indirect geothermometers (Arnórsson and Gunnlaugsson, 1985; Polyak and Tolstikhin, 1985; Kharaka and Mariner, 1989; Harvey and Browne, 1991). In the first case, the subsurface temperature distribution is usually obtained through spatial interpolation based on the actual temperature logging. However, this method often leads to large errors, especially in areas with complex geologic structures. In the second case, the construction of the temperature models needs to (1) determine the heat transfer mechanism (i.e., determine whether the heat transfer is via conduction, convection, or a combination of both), (2) define the boundary conditions, and (3) accurately define the thermophysical parameters of each medium. This method usually involves considerable error because the thermophysical parameters and boundary conditions can only be roughly estimated (especially the lower boundary conditions). In addition, another alternative method uses indirect geothermometers based on geologic (Harvey and Browne, 1991), geochemical (Kharaka and Mariner, 1989), isotopic (Polyak and Tolstikhin, 1985), or gas composition (Arnórsson and Gunnlaugsson, 1985) data collected on the surface. Although these geothermometers can serve as valuable tools for temperature assessment at some characteristic depths, they cannot estimate the temperature distribution in an area, nor can they be used to interpolate/extrapolate temperature well logs (Spichak and Zakharova, 2012).

Existing methods of temperature estimation are incapable of providing reliable accuracy at depths below the reach of the boreholes and in the cross-borehole space. This difficulty can be overcome by using the temperature-dependent (proxy) petrophysical parameters. In addition, previous results have shown that, among all the types of petrophysical parameters, the resistivity of the rock is the most sensitive to temperature changes, and the temperature dependence of the electrical conductivity of rocks justifies its application to geothermometry (Waff, 1974; Shankland and Waff, 1977; Zhang et al., 2011), i.e., the conversion of the electrical conductivity to temperature values. There are two types of methods (Hermance and Grillot, 1970; Spichak et al., 2007, 2011; Spichak and Manzella, 2009; Meng et al., 2010; Zhang et al., 2011) for estimating the subsurface temperature using a proxy parameter: (1) artificial neural networks (ANNs) and (2) purely empirical formulas. For the first type, although the ANN model has a high subsurface temperature prediction accuracy in the depth range covered by the data samples, the temperature prediction may be inaccurate in the depth ranges without data samples. In addition, a large amount of sample data is required to ensure that the temperature prediction has a reliable accuracy. In the latter case, the validity of the empirical formulas is unjustifiably postulated to be invariant with respect to the spatial coordinates; that is, each parameter in the empirical formula is assumed to be constant in any geologic setting and at all depths.

Based on the disadvantages of the existing methods for predicting the temperature of the earth’s interior, we make full use of the advantages of resistivity in temperature prediction (Waff, 1974; Shankland and Waff, 1977) and the effectiveness of electromagnetic (EM, including magnetotelluric [MT], audio-frequency MT [AMT], and controlled-source AMT [CSAMT]) method in detecting the geoelectrical structure of geothermal fields (Munoz, 2014; Spichak, 2015). An EM geothermometer based on the coefficient correction method of the optimal temperature (CCMOT) is proposed, and the key parameters of the model are improved so that they change with the geologic settings and spatial coordinates. First, the optimal correction temperature (OCT), intrinsic normalized resistivity (INR), and temperature correction coefficient (TCC) for the different depths are derived based on the resistivity-temperature logging data of the study area or adjacent area, and regression analysis is performed on these three key parameters (i.e., $T_0(z)$⁠, $R_{NT0}(z)$⁠, and ${\alpha }_{T0}(z)$⁠). Then, resistivity inversion of the EM data (MT data are used in this research) from the study area is carried out, and the resistivity profile of the EM inversion also is converted to normalized resistivity $R_{N\mathrm{inv}}(x,z)$⁠. Finally, combined with $T_0(z)$⁠, $R_{NT0}(z)$⁠, ${\alpha }_{T0}(z)$⁠, and $R_{N\mathrm{inv}}(x,z)$⁠, the precise temperature prediction model at different spatial coordinates in the study area is constructed, and the subsurface temperature field is calculated layer by layer. This method has strong applicability and practicality, as well as a wide/deep prediction range (corresponding to the coverage and depth of the EM data). Thus, the CCMOT is of great significance to geothermal resource exploration and evaluation, and it also may enable monitoring of the variations in subsurface temperatures based on ground EM monitoring data and the forecasting of future trends during the exploitation of geothermal reservoirs.

The Xiong’an New Area is located in the hinterland of Baoding, Beijing, and Tianjin, China. The administrative division includes Xiong County, Rongcheng County, Anxin County, and some surrounding areas in Hebei Province. The regional geologic structure (Figure 1a) is complex and generally controlled by the tectonic framework of northern China. It is characterized by multiple convex and concave structures developed alternately in an northeast-trending uplift and fault depression pattern (Xiong and Zhang, 1988; Chen et al., 1990; Sun and Niu, 2000). Several studies have suggested that, in the fault depression in northern China, the concave-convex structural pattern and the spatial distribution of the geothermal resources are controlled by a mantle subplume that is mushroom-shaped with an upper surface that reaches a depth of approximately 50 km (Niu et al., 2001).

As shown in Figure 1b, the Xiong’an New Area is tectonically located at the junction between the Rongcheng uplift and the Niutuozhen uplift in the northern part of the Jizhong Depression, Bohai Bay Basin, northern China. The Langfang-Gu’an sag is located to the north, the Xushui sag and Baoding sag are located to the west, the Baxian sag is located to the southeast, and the Raoyang sag and Gaoyang low uplift are located to the south. Several faults control the Xiong’an geothermal field due to the multistage tectonic movement (Cao et al., 2021); the main fault structures include the Niunan fault, Niudong fault, Daxing fault, Xushui fault, and Rongcheng fault. The northeast-trending Rongcheng fault separates the Niutuozhen uplift and the Rongcheng uplift, and the east–west-trending Xushui-Niunan fault separates the two uplifts (Rongcheng uplift and Niutuozhen uplift) from the Gaoyang low uplift. These faults act as heat and water conduction channels and play an important role in the formation of the geothermal system in the Xiong’an New Area (Wu et al., 2018; Liu et al., 2020). Collectively, these details make this a promising area for geothermal energy exploitation. In addition, this area is a typical carbonate rock thermal reservoir (Liu et al., 2020), and there are numerous valuable borehole logging data and MT data for the Xiong’an New Area (Wu et al., 2018). As such, the Xiong’an New Area is selected as the priority study area for this current study.

The preceding behavior for rock resistivity with temperature is primarily for moderate temperature formations because when the temperature is above 200°C, the negative correlation may disappear and the resistivity variation becomes insensitive to the temperature changes (Ucok et al., 1980; Llera et al., 1990; Ussher et al., 2000; Kristinsdóttir et al., 2010). In addition, if the preceding empirical formulas are directly used to predict the subsurface temperature distribution, there are two significant drawbacks.

1. The estimation of the subsurface temperature from untransformed resistivity data (measured resistivity) will result in a large error. For example, due to the differences in lithology, if the resistivity at depth is higher than that of the shallow surface, and the untransformed resistivity is directly used to predict the temperature distribution, then it will be determined that the temperature at depth is lower than that of the shallow surface (in general, rocks have ionic conductivity properties and their temperature is inversely proportional to their resistivity), leading to an erroneous result. In addition, logging resistivity and EM inversion resistivity can effectively reflect the resistivity variation trend of a specific formation, but the resistivity values of both may not be the same.

   To address this, the relationship between the “normalized resistivity” and temperature for different formations is determined, which overcomes the preceding problem and reduces the large deviation in the temperature prediction caused by the disturbance of the untransformed resistivity.

2. The three key parameters, i.e., the base temperature (⁠$T_{25}$ and $T_{18}$⁠), the intrinsic resistivity (⁠$R_{25}$ and $R_{18}$⁠), and the temperature coefficient (⁠${\alpha }_{25}$ and ${\alpha }_{18}$⁠), are fixed values in any geologic setting and depth range. Thus, predicting the temperature of the earth’s interior using these fixed values would undoubtedly lead to large errors.

   Here, based on the temperature-resistivity logging data from the study area or an adjacent area (as for logging resistivity, the resistivity of deep investigation laterolog after the borehole correction is usually used), two sets of constraint equations are selected to derive the OCT, INR, and TCC at different depths in the study area, thereby overcoming the problem that the validity of the empirical formulas is unjustifiably postulated to be invariant with respect to the spatial coordinates (geologic settings).

In a specific stratigraphic segment, the normalized resistivity of the medium at a specific temperature (⁠$T_{\mathrm{test}}$⁠) is used to convert the normalized resistivity of all the measuring points with different temperatures, so that the observed values are closest to the calculated values. This temperature is called the OCT, where $T_0=T_{\mathrm{test}}$⁠. The operation steps of the CCMOT are described next.

• The resistivity-temperature logging data (from $N_L$ [$N_L\ge 1$] boreholes in the study area or its adjacent areas) are divided into N depth segments, and each segment is $L_D$ (m) long. The resistivity of each segment is normalized. It should be emphasized that the borehole correction to the logging resistivity is required; otherwise, the data would be “normalized” to the incorrect quantity.

• The OCT, INR, and TCC of the different segments in the different boreholes are derived, and the regression analysis is performed on the preceding key parameters at corresponding depths. We then obtain $T_0(z)$⁠, $R_{NT0}(z)$⁠, and ${\alpha }_{T0}(z)$⁠.

• Inversion of EM data for the study area is carried out to obtain the resistivity (⁠$R_{\mathrm{inv}}(x,z)$⁠) distribution of any node $P(x,z)$ in the underground space. The resistivity profile is divided into M (integer $M\ge 2$⁠) depth segments, and each segment has $M_D(\mathrm{m})$⁠. The resistivities of the different segments are then normalized to $R_{N\mathrm{inv}}(x,z)$⁠.

• Based on $T_0(z)$⁠, $R_{NT0}(z)$⁠, and ${\alpha }_{T0}(z)$ (obtained in step 2), as well as $R_{N\mathrm{inv}}(x,z)$ (obtained in step 3), the precise temperature prediction model of the underground space in the study area can be obtained, and the subsurface temperature distribution $T(x,z)$ can be calculated layer by layer using $R_{N\mathrm{inv}}(x,z)$⁠.

Figure 2 provides a flow chart of the temperature prediction of an underground space. The details of the derivation are as follows.

In step 2, the detailed process of calculating the OCT, INR, and TCC of the different segments (N segments in total) is performed as described next.

1. Based on the improved empirical formula, the relationship between the normalized resistivity and temperature can be expressed as

   $\frac{R_{NT0}}{R_{\mathrm{Norm}}(i)}=1+{\alpha }_{T0}(T(i)-T_0),$

   (4)

   where $R_{NT0}$ is the INR, $R_{\mathrm{Norm}}(i)$ represents the normalized resistivity logging data sets for all of the measurement points in the ith segment of the borehole, ${\alpha }_{T0}$ is the TCC, $T(i)$ is the logging temperature data sets for all of the measurement points in the ith segment of the borehole, and $T_0$ is the OCT.
2. Equation 4 can be converted to (⁠$1/(R_{\mathrm{Norm}}(i))$⁠) $=(({\alpha }_{T0})/(R_{NT0}))(T(i))+((1-{\alpha }_{T0}{T}_0)/(R_{NT0}))$⁠, and the slope and intercept are defined as B and A, respectively. Therefore,

   $\{\begin{array}{c}\frac{{\alpha }_{T0}}{R_{NT0}}=B\\ \frac{1-{\alpha }_{T0}{T}_0}{R_{NT0}}=A\end{array},$

   (5)

   where A and B are the coefficients related to geologic settings and electrical characteristics.
3. Based on a large number of core experiments, the resistivity of rock is logarithmically related to temperature (Qu et al., 2001; Pan et al., 2016):

   $R_{\mathrm{Norm}}(i)=-C\ln (T(i))+D,$

   (6)

   where C and D are the coefficients related to geologic settings and electrical characteristics.
4. The term T_test is substituted into equations 4 and 6 for conversion, and it is combined with equation 5 to obtain equation 7:

   $\frac{1}{A+B{T}_{\mathrm{test}}}=-C\ln (T_{\mathrm{test}})+D,$

   (7)

   where $T_{\mathrm{test}}$ represents a specific temperature, which is used to convert the normalized resistivity (logging data) of all the measuring points with different temperatures.
5. Taking the derivative of equation 7 with respect to $T_{\mathrm{test}}$ gives

   $B^{2}C{T}_{\mathrm{test}}^2+(2ABC-B)T_{\mathrm{test}}+A^{2}C=0.$

   (8)
6. By solving equation 8, we obtain $T_{\mathrm{test}}$⁠:

   $T_{\mathrm{test}}=\frac{-(2ABC-B)+\sqrt{(2ABC-B)-4A^2{B}^2{C}^2}}{2B^{2}C}.$

   (9)

7. Because these coefficients A, B, C, and D can be obtained from the normalized resistivity and temperature data for the boreholes by calculating the OCT of a specific segment, the corresponding INR and the TCC can be calculated using equations 6 and 5, respectively.

Before using the MT data for indirect temperature estimation, dimensional analysis should be carried out. As is known, the lateral geologic irregularities in a medium manifest as frequency dependencies in the dimensionality indicators of the medium. That is, the MT field transforms to zero in a locally homogeneous medium, and the more it differs from zero, the more laterally inhomogeneous the medium is (Spichak et al., 2011).

In this study, the telluric-vector technique is used to distinguish between the local telluric distortion and the regional induction using the phase information from all of the impedance tensor elements (Bahr, 1988, 1991). The regional strike and dimensionality of the MT data prior to the 2D inversion also are estimated.

As shown in Figure 3, most of the dimensional indicators of the three MT profiles are less than 0.3. Therefore, 2D inversion can be used for the three MT profiles. In addition, Figure 4 shows the rose diagrams used to depict the strike directions of all the sites for all of the frequencies. We obtain a principal tensor direction of almost N45°E (or N45°W) with a 90° ambiguity. The strike direction along profile B is consistent with that along profile C. Because the strike of the regional fault plane in the study area supports a predominantly northeast direction (Lei et al., 2020), we select N42°E for profile A, N45°E for profile B, and N45°E for profile C.

Figure 5 shows the skin depths (a transverse electric [TE] mode and a transverse magnetic [TM] mode) of the three MT profiles. The effective penetration depth of MT signals is less than 10 km; thus, we study the electrical structure characteristics within a depth of 5 km in this study.

The 2D inversion analysis of the rotated data is carried out using the nonlinear conjugate gradient algorithm of the WingLink data interpretation package, which was developed by Rodi and Mackie (2001). Furthermore, all of the frequencies of the MT data are rotated to the respective strike directions before the resistivity inversion. A series of inversion parameters need to be optimized to appropriately fit the data before calculating the final model. In addition, TM-mode data are relatively insensitive to the depth extent of a subsurface body (Eberhart-Phillips et al., 1995), and TE-mode data are more sensitive to 3D effects (Wannamaker et al., 1984). Thus, mixed mode (joint TE and TM modes) analysis is used to help clarify the geologic structures beneath the three profiles (Figure 6).

1. The resistivity (Figure 6a) of profile A shows a good correlation between the resistivity distribution and the formation lithology. There is a stratigraphic boundary based on a resistivity difference at a depth of approximately 1000 m, and from the surface to a depth of approximately 1000 m the resistivity is generally less than 10 ohm·m, corresponding to the Quaternary and Neogene strata. Another important stratigraphic boundary (dotted black line in Figure 6a) exists in the depth range of 2000–3000 m. Above the dotted black line, the resistivity of the strata is relatively low, generally less than 100 ohm·m, corresponding to the strata of the Jixian system (Jxw). However, below the dotted black line, the formation’s resistivity increases rapidly to more than 100 ohm·m. At position nos. 3–9, the high resistivity reaches >200 ohm·m, corresponding to the old strata of the Changcheng system (Chg). The thermal reservoir of Jxw is attached to the sand and mudstone strata of Tertiary/Quaternary with poor thermal conductivity, which can be used as an excellent thermal insulation layer (Wang et al., 2012; Cao et al., 2021). At the same time, the thermal reservoir of Jxw has a large thickness (Figure 6a) and well-developed fractures and caves (Cao et al., 2021). Thus, the preceding characteristics of reservoir and caprock can provide favorable conditions for forming a geothermal system.

   In the depth range of 2000–5000 m, below position nos. 3–9, there is a high-resistivity anomaly (i.e., an upward bulge). We conclude that the Rongcheng uplift caused the uplift of the old strata and changed the electric structure. Moreover, the basement crust in the Jizhong Depression has a multilayer structure, which can be generally divided into high-speed and low-conductivity upper crust, low-speed and high-conductivity middle crust, and high-speed and low-conductivity lower crust (Wang et al., 2021). Among them, the high-conductivity layer (low-resistivity layer) in the crust is tectonically active strata, which is associated with the deep faults that cut through the crust and even the lithosphere to transfer heat from the mantle and become a channel for the upwelling of the heat source material of deep mantle, but its buried depth is generally 15–22 km (Wang et al., 2021). Most of the thermal reservoirs in Xiong’an New Area have lower porosity and permeability strata (Dai et al., 2020; Wang et al., 2021). In addition, the fault structures control the development of karsts and the overall pattern of karst palaeogeomorphology in this study area. Therefore, local thermal anomalies in the study area are caused primarily by fracture structure and karst development (Lu et al., 2019; Han et al., 2020), which enhance groundwater circulation heated by deep strata and makes it easier for the deep heat to reach the shallow layers and then promotes the formation of a thermal anomaly region with obvious resistivity change. Thus, according to the distribution characteristics of the relatively low-resistivity layer near the surface and the geologic data, it is inferred that there are faults on both sides of the uplift area (Wu et al., 2018), namely, fault F1 (Rongxi fault) and fault F2 (Xushuinan fault). These faults can be used as good channels for the upward migration of deep geothermal water (Guo et al., 2020).

2. As shown in Figure 6b, on resistivity profile B and above the dotted black line, the low resistivity (<40 ohm·m) reflects the electrical characteristics of the Cenozoic strata (Kz). Below the dotted black line, the high resistivity corresponds to the Proterozoic bedrock (Pt). There is a noticeable discontinuity in the transverse resistivity, which indicates that the resistivity changes greatly under the influence of the geologic structure. Based on the tectonic characteristics of the area and the fact that at a depth of approximately 3500 m below position nos. 8–13 the resistivity changes dramatically with depth (changes from low to high resistivity), it is inferred that the F3 fault separates these resistivity zones. In addition, at a depth of approximately 2800 m below position nos. 18–45 in profile B, which crosses the Niunan fault (Figure 1c), there is an obvious lateral discontinuity in the resistivity profile (below position nos. 24–33), which is presumed to be caused by the Niunan fault. These faults also can serve as a good channel for upward migration of deep geothermal water.

3. Figure 6c shows the resistivity of profile C. The resistivity gradually increases from the near surface to the middle-deep strata. The dotted black line marks the boundary between the Kz strata and the Pt strata. The resistivity above the dotted line is less than 40 ohm·m, reflecting the Kz strata; however, the resistivity below the dotted line is relatively high, corresponding to the bedrock formation (Pt strata). Compared to the resistivity of profile B, the transverse electrical characteristics of profile C have good continuity. There are no obvious electrical structures that indicate the existence of faults, and the strata are relatively flat. Therefore, the heat source of profile C and its vicinity may come from heat conduction and radioactive heat generation of the deep crust (Cao et al., 2021).

The temperature prediction model is obtained by combining equations 9–12. The normalized resistivity of the three MT profiles is obtained using step 3, and the corresponding temperature distribution of the underground space is calculated layer by layer. Figure 8 shows the temperature fields of the three profiles.

1. According to the visualized temperature distribution, within the depth range of 2 km (⁠$z\le 2\mathrm{km}$⁠), the temperature characteristics can be depicted with a higher resolution than those at depths of greater than 2 km. This phenomenon can be explained as follows. The accuracy of the temperature estimation is closely related to the inversion accuracy of the EM data. In the shallow strata, the electrical characteristics can be more finely delineated. In contrast, in the deep strata, the resolution of the electrical structures of the geologic bodies decreases, so the resolution of the temperature prediction also decreases. Thus, it is not easy to characterize the detailed characteristics of the temperature field. In addition, geothermal resources can be divided into three categories based on temperature: high temperature (>150°C), moderate temperature (90°C–150°C), and low temperature (25°C–90°C) resources. As such, the favorable thermal reservoir areas for low-temperature hot water (60°C–90°C) are mainly concentrated in the depth range of 1.5–3.0 km. At a depth of 5 km in the study area, the temperature is approximately 125°C–130°C.

2. Based on the resistivity profiles and temperature fields (Figure 6), we found that the local low-resistance region corresponds to the local high-temperature region. For example, at a depth of approximately 200–500 m below position no. 8 in profile A (Figure 8a) and in the depth range of 500–1000 m below position nos. 6–21 in profile B (Figure 8b) are local high-temperature regions. In the depth range of 1000–2000 m on profile C (Figure 8c), local high-temperature regions of different scales also exist below nos. 3–9 and 21–33. At the same time, it is evident that the high-resistance region corresponds to the local low-temperature area (in the depth range of 4–5 km below position nos. 67–72 in profile C).

   In addition, the temperature field also exhibits a concave pattern between the two fault zones, whereas on the roof of the two fault zones and on the two sides of the fault zone, the temperature exhibits a convex pattern (Figure 8a and 8b), which is consistent with the high conductor temperature propagation properties in heat conduction (Yue et al., 2019; Wang et al., 2020). In addition, the temperature distribution (profile A) is similar to the temperature field (profile I-I′, as is shown in Figure 1) calculated by Wang et al. (2020). That is, the temperature is locally high on both sides of the fault zone and on the roof of the two fault zones, which further indicates that these faults, as heat and water conduction channels, play important roles in the formation of the geothermal system in the Xiong’an New Area. In addition, at a depth of approximately 1000–2500 m (near the faults) below position nos. 5–10 of profile A, there is a very evident local high-temperature region; therefore, we can infer that the faults near profile A have more active thermal convection than the faults near profile B. However, the temperature distribution of profile C may have typical heat conduction characteristics (the heat source may come from the deep crust). Meanwhile, we notice that the geothermal gradients are higher in the shallow formation (<1000 m), further illuminating that the caprock in the research area can be regarded as a suitable medium for heat insulation.

3. To verify the effectiveness of using the CCMOT for temperature prediction, the goodness-of-fit (GOF) is used to determine the accuracy (equation 15):

   $R^2=\frac{\mathrm{SSR}}{\mathrm{SST}},$

   (15)

   where

   $\mathrm{SSR}=\sum _{i=1}^n(T_{\mathrm{pre}}-\overline{T_{\mathrm{act}}}{)}^2,$

   (16)

   $\mathrm{SST}=\mathrm{SSR}+\mathrm{SSE},$

   (17)

   $\mathrm{SSE}=\sum _{i=1}^n(T_{\mathrm{act}}-T_{\mathrm{pre}}{)}^2,$

   (18)

   where R² is the GOF, SSR represents the regression sum of squares, SST represents the sum of squares of total, SSE represents the residual sum of squares, T_pre denotes the predicted temperature, T_act denotes the measured temperature, and $\overline{T_{\mathrm{act}}}$ denotes the mean value of the measured temperature.

The GOF value ranges from 0 to 1, and the closer the GOF value is to 1, the better the fitting degree between the predicted and measured values. The GOF values of profile A (D16, D17), profile B (D35), and profile C (D32) are 0.8669/0.9725, 0.8997, and 0.8804, respectively (Figure 9). Therefore, the accuracies of the temperature predictions of the three profiles range from 86.69% to 97.25%. However, it should be emphasized that although the CCMOT has achieved good results in predicting the subsurface temperature field in the whole Xiong’an New Area, its application range should be determined according to the actual situation in other areas, especially in areas with complex geologic structure. In addition, under the same conditions, we apply an indirect EM geothermometer based on an ANN (Spichak, 2015) to predict the temperatures along the three profiles. Four constraint boreholes (i.e., D16, D17, D32, and D35) are used as training samples. In addition, the optimal training strategy is chosen to train the neural network (Spichak et al., 2011; Spichak, 2015): the input data consist of the electrical resistivity values and the appropriate space coordinates where they are determined, whereas the output data are the temperature values determined for the same locations. The GOF values of profile A (D16/D17), profile B (D35), and profile C (D32) are 0.8136/0.8818, 0.7621, and 0.8443, respectively (Figure 9). Thus, for this ANN model, the accuracies of the temperature predictions of the three profiles range from 81.36% to 88.18%. Therefore, compared to the ANN model, the temperature prediction accuracy obtained using the CCMOT is 5.33%–9.07% better. However, it may not be a universal conclusion that the CCMOT is superior to ANN in predicting subsurface temperature field. In the case of fewer borehole data, we believe CCMOT can provide higher accuracy in temperature prediction; however, the ANN could perform well with a large enough sample size. In addition, the GOF is used to define the accuracy of predicting temperature fields (profile I-I′) from heat transfer modeling (Wang et al., 2020). The GOF value (D16/D17) of profile I-I′ is 0.8591–0.8625, which is lower than the accuracy of temperature prediction using CCMOT. However, this does not mean that the heat transfer modeling is less effective than CCMOT in predicting temperature because it requires a significant amount of prior information (including logging data, core testing, and stratum characteristics) to constrain the boundary conditions and thermal physical parameters. With sufficient prior information, we believe that the heat transfer modeling can better describe the subsurface temperature field because this method can recover the temperature distribution characteristics of the research area from the perspectives of thermal dynamics (the mechanism of geothermal energy accumulation). Based on this, we conclude that the CCMOT has a strong practicality and applicability and a high accuracy for temperature prediction in underground spaces. Thus, CCMOT will help to advance the exploration, evaluation, development, and use of geothermal resources.

To study the stability and sensitivity of the CCMOT for temperature prediction, we study the influences of the number of constraint boreholes (⁠$N_L$⁠), the borehole location, the subdivision spacing (⁠$L_D$⁠) of the logging data, and the spacing (⁠$M_D$⁠) of the grid node of the inversion resistivity profile along the depth direction on the accuracy of the temperature prediction. In addition, to comprehensively compare the degrees of influence of the preceding parameters on the accuracy of the temperature estimation, based on the principles of mathematical statistics, we adopt the coefficient of variation as the index of the sensitivity evaluation (Abdi, 2010; McAuliffe, 2015).

First, the corresponding $T_0(z)$⁠, $R_{NT0}(z)$⁠, and ${\alpha }_{T0}(z)$ of six boreholes (i.e., D11, D12, D16, D17, D32, and D35) with $L_D=200\mathrm{m}$ are calculated. Then, different temperature prediction models are generated based on the permutations and combinations (⁠$N_L=1$⁠, 2, 3, 4, and 5) of the preceding three key parameters. We calculate the corresponding temperature distributions (⁠$M_D=100\mathrm{m}$⁠) of the three profiles and obtain the GOF values (see Table 1). In addition, for the same number of boreholes (⁠$N_L$⁠), the average value and standard deviation of the GOF are calculated (Figure 10).

The mean value of the GOF is positively correlated with $N_L$⁠, which is because the larger the value of $N_L$⁠, the more valuable regional temperature field information can be extracted. Otherwise, only the local temperature field information can be extracted from a small number of boreholes. Thus, the temperature prediction becomes more accurate as $N_L$ increases.

When the constraint boreholes (logging data) are near the temperature profile to be predicted, the GOF will be greater. Taking verification borehole D17 of profile A as an example (⁠$N_L=1$⁠), if D16 is involved in the constraint calculation (⁠$D_{cp}=0.3\mathrm{km}$⁠), the GOF is 0.910. However, if D11 is used as the constraint data, which is far away from (⁠$D_{cp}=12.6\mathrm{km}$⁠) profile A, the GOF is only 0.719 (see Table 1). We found that when N_L = 2, the average GOF value of profile A (D17) is a maximum value. This may be because constraint boreholes D16 and D17 are close to profile A, i.e., 0.3 km and 0.2 km away, respectively (see Table 2). Based on the principle of permutation and combination, among the six constraint boreholes, the probability of the logging data for D16 or D17 being selected as constraint data is 2/3 when $N_L=2$⁠; therefore, the GOF values are relatively large in more than 2/3 of the cases, and the average GOF value is correspondingly greater.

To study the influence of $L_D$ on the temperature prediction accuracy, first we calculated the corresponding $T_0(z)$⁠, $R_{NT0}(z)$⁠, and ${\alpha }_{T0}(z)$ of six boreholes (i.e., D11, D12, D16, D17, D32, and D35) with different $L_D$⁠. We then calculate the temperature distributions of the three profiles with $N=6$ and $M_D=100\mathrm{m}$ under different $L_D$⁠. Figure 11 shows the GOF for the different $L_D$⁠. As $L_D$ increases, the GOF value generally decreases. However, there is an optimal value of $L_D$ (⁠$L_D=L_{\mathrm{best}}$⁠), which results in the maximum GOF, but the value of $L_{\text{best}}$ is different when the GOF value reaches the maximum value for the different verification boreholes. The maximum GOF values of profile A (D16 and D17) and profile B (D35) occur when $L_D=200\mathrm{m}$⁠, but for profile C (D32), the maximum GOF value occurs when $L_D=250\mathrm{m}$⁠. An optimal value (⁠$L_D=L_{\mathrm{best}}$⁠) results in the maximum GOF value, which can be explained by the following reasons. The smaller $L_D$⁠, the more likely it is that local information about the temperature field below the borehole location can be extracted. The larger $L_D$⁠, the less valuable temperature field information can be extracted. Therefore, both cases can lead to inaccurate temperature predictions, and an optimal value (⁠$L_D=L_{\mathrm{best}}$⁠) can make the temperature prediction more accurate. The difference in the temperature prediction accuracy is approximately 32.3%–46.6% when $L_D$ is 50–400 m.

As described previously, the accuracy of the temperature prediction is closely related to the resistivity inversion of the EM data. Thus, it is necessary to study the influence of $M_D$ on temperature estimation. All six boreholes (i.e., D11, D12, D16, D17, D32, and D35) in this study are used as constraint data to calculate the corresponding key parameters (⁠$T_0(z)$⁠, $R_{NT0}(z)$⁠, and ${\alpha }_{T0}(z)$⁠) with $L_D=200\mathrm{m}$⁠. Combined with the normalized resistivity distribution characteristics with different $M_D$⁠, the corresponding normalized resistivity-temperature relationship is obtained, and the temperature distributions are calculated. Figure 12 provides the GOF values for different $M_D$⁠. We find that the GOF is inversely correlated with $M_D$ because the grid nodes of the temperature profile correspond to the grid nodes of the resistivity profile. That is, the larger $M_D$⁠, the larger the spacing of the temperature profile along the depth direction, which leads to a less precise temperature field, so GOF decreases with increasing $M_D$⁠. The difference in the temperature prediction accuracy is approximately 30% when $M_D$ is 100–400 m.

The value of $C_v$ reflects the sensitivity of the GOF to the various influencing factors. The greater the value of $C_v$⁠, the greater the influence of the factor on the GOF.

Figure 13 summarizes the results of the sensitivity comparison of the studied parameters. In the studied parameter ranges, $L_D$ has the greatest degree of sensitivity, followed by $M_D$⁠, and $N_L$ has the least influence. The influence of the parameter $N_L$ on the GOF is relatively small, indicating that the CCMOT used in this study has relatively little dependence on the number of boreholes. Before the temperature prediction, the parameter $L_D$ (the coefficient of variation for $L_D$ is 17.62%–21.46%) must be optimized to ensure the accuracy of the prediction. In addition, the coefficient of variation for $M_D$ is 13.69%–17.47%; thus, it is demonstrated that the accuracy of the resistivity inversion of the EM data is vital to the temperature prediction.

In addition, it should be emphasized that the $N_L$ in the study area (or adjacent area) is significant for the prediction results of subsurface temperature. However, among the three sensitivity parameters in this research (⁠$N_L$⁠, $L_D$⁠, and $M_D$⁠), the $N_L$ has a relatively weak influence on temperature field prediction. Therefore, this also is an indirect indication that CCMOT is superior to ANN with fewer constrained boreholes: more stable and precise subsurface temperature prediction can be obtained. However, it is a disadvantage of the CCMOT that the predicted temperature distribution is strongly dependent on the $L_D$⁠.

This study demonstrates that the proposed CCMOT can effectively transform the resistivity characteristics into a visualized temperature distribution for middle-low temperature geothermal fields (<200°C) with an accuracy of 86.69%–97.25% in the Xiong’an New Area; the accuracy is 5.33%–9.07% better than that of the optimal ANN model. The CCMOT is convenient and practical, and its prediction range is wide and deep (consistent with the coverage range and depth of EM detection).

In addition, the coefficient of variation indicates that, in the studied parameter ranges, $L_D$ has the greatest sensitivity, followed by $M_D$ and then $N_L$⁠. (1) The CCMOT is relatively less dependent on the number of boreholes used as constraint data than ANN, but it is more advantageous if the constraint borehole(s) is (are) close to the profile to be predicted. (2) It is better to optimize $L_D$ before the temperature prediction because there is an optimal subdivision spacing (⁠$L_{\text{best}}$⁠) of the logging data, which results in the maximum GOF. The difference in the temperature prediction accuracy is approximately 32.3%–46.6% when $L_D$ is 50–400 m. (3) High-quality EM data are the key to obtaining a high-precision temperature prediction.

The concept of normalized resistivity of the CCMOT can be used as a reference for other temperature prediction models, such as normalized wave velocity-temperature models. However, the CCMOT focuses primarily on building a one-to-one relationship between the normalized resistivity and temperature. Future studies will establish a new composite model in which the normalized resistivity varies with temperature and other variables (e.g., pressure, fluid composition, and other microstructural parameters) in an underground space.

The authors gratefully appreciate the financial support from the National Key Research and Development Program of China (2018YFC0604303 and 2018YFC1503705) and the National Nature Science Foundation of China (41630317). We also acknowledge the editors and reviewers for their constructive suggestions and comments.

Data associated with this research are confidential and cannot be released.

NOMENCLATURE
 
• A, B, C, D

  coefficients related to geologic settings and electrical characteristics

 
• CCMOT

  coefficient correction method of the optimal temperature

 
• C_v

  coefficient of variation, equation 19

 
• Dcp

  distance from constraint borehole to MT profile

 
• INR

  intrinsic normalized resistivity

 
• L_best

  optimal subdivision spacing of logging data, which results in the maximum R²

 
• $L_D$

  length of each segment after the logging data is divided into N depth segments and step 1 in Figure 2 

 
• $\max [R_{\log }(i)]$

  maximum value of R_log(i)

 
• $\max [R_{\mathrm{inv}}(k)]$

  maximum value of R_inv(k)

 
• $M_D$

  length of each segment after the resistivity profile is divided into M depth segments and step 3 in Figure 2 

 
• n

  the number of sample data, equation 16

 
• N_cv

  the number of sample data, equation 20

 
• $N_L$

  number of constraint boreholes and step 1 in Figure 2 

 
• OCT

  optimal correction temperature

 
• P(x, z)

  coordinates of the grid node of the EM profile (or temperature profile) and x represents the horizontal position, in which z denotes the depth

 
• R²

  goodness-of-fit, equation 15

 
• $R_{18}$

  resistivity of rock at a temperature of 18°C, equation 2

 
• $R_{25}$

  resistivity of rock at a temperature of 25°C, equation 1

 
• $R_{\mathrm{inv}}(k)$

  data sets of resistivity in the kth (k = 1, 2, 3, …, M) segment of the EM profile

 
• $R_{\mathrm{inv}}(k,x,z)$

  resistivity at P(x,z) in the kth (k = 1, 2, 3, …, M) segment of the EM profile

 
• $R_{\log }(i,j)$

  resistivity logging data of the jth measuring point in the ith (i = 1, 2, 3,…, N) segment of the borehole

 
• $R_{\log }(i)$

  logging resistivity data sets for all of the measurement points in the ith (i = 1, 2, 3,…, N) segment of the borehole

 
• $R_{\text{Ninv}}(k,x,z)$

  normalized resistivity at P(x, z) in the kth (k = 1, 2, 3, …, M) segment of the EM profile

 
• $R_{\text{Ninv}}(x,z)$

  normalized resistivity at P(x, z) of the EM profile

 
• $R_{\mathrm{Norm}}(i)$

  normalized resistivity logging data sets for all of the measurement points in the ith segment of the borehole

 
• $R_{\mathrm{Norm}}(i,j)$

  normalized resistivity logging data of the jth measuring point in the ith (i = 1, 2, 3, …, N) segment of the borehole

 
• $R_{NT0}$

  intrinsic normalized resistivity

 
• $R_T$

  resistivity of a rock at a given temperature (T°C), equations 1 and 2

 
• s

  standard deviation of the sample data, equation 19

 
• SSE

  residual sum of squares, equation 18

 
• SSR

  regression sum of squares, equation 16

 
• SST

  sum of squares of total, equation 17

 
• $T_0$

  optimal correction temperature

 
• $T_{18}$

  base measuring temperature, T₁₈ = 18°C, equation 2

 
• $T_{25}$

  base measuring temperature, T₂₅ = 25°C, equation 1

 
• $T_{\mathrm{act}}$

  measured temperature

 
• $\overline{T_{\mathrm{act}}}$

  mean value of the measured temperature

 
• T_test

  a specific temperature, which is used to convert the normalized resistivity (logging data) of all the measuring points with different temperatures

 
• $T(i)$

  logging temperature data sets for all of the measurement points in the ith segment of the borehole

 
• $T_{\mathrm{pre}}$

  predicted temperature

 
• $T(x,z)$

  predicted temperature at P(x,z)

 
• TCC

  temperature correction coefficient

 
• $\bar{y}$

  average of sample data, equation 19

 
• $y_i$

  a data sample, equation 20

 
• ${\alpha }_{18}$

  temperature coefficient generally estimated as 0.018, equation 2

 
• ${\alpha }_{25}$

  temperature coefficient generally estimated as 0.025, equation 1

 
• ${\alpha }_{T0}$

  temperature correction coefficient

Biographies and photographs of the authors are not available.