Validation and Evaluation of an Estimation Method for Deep Thermal Structures Using an Activity Index in Major Geothermal Fields in Northeastern Japan

: A considerable number of rock bodies with varying percentages of supercritical ﬂuid exist around the brittle–ductile transition (BDT) zone at a depth of several kilometers from the surface of the Earth, in northeastern Japan. As the BDT zone in the granitic basement of the continental crust is estimated to occur at about 380 ◦ C, the identiﬁcation of the depth corresponding to 380 ◦ C is important to utilize the thermal energy inside the “supercritical geothermal systems”. In this study, we focused on an estimation method to determine the depth of the isothermal layer corresponding to 380 ◦ C, using the activity index (AI) obtained from the maximum-temperature data of the geothermal wells and hot springs. The thermal proﬁles of deep and hot exploration boreholes and the hypocentral distribution of natural earthquakes were used to evaluate the characteristics and accuracy of the deep thermal structure, using the activity index. The estimated depth corresponding to 380 ◦ C tended to be higher than the actual depth, with a maximum possible estimation error of approximately 9.7 km. Distribution maps showing the depth of the isothermal layer corresponding to 380 ◦ C were created for six major geothermal ﬁelds in northeastern Japan, using the results from this study.


Introduction
It has been established that supercritical geothermal systems, which are located at depths near or below the brittle-ductile transition (BDT) zone in the granitic basement and that contain supercritical fluid (i.e., at temperatures and pressures >374 • C and >22.1 MPa for pure water and >406 • C and >29.8 MPa for seawater, respectively), have the potential for a huge amount of power generation and subsequent reduction of CO 2 emissions in many countries with high-temperature geothermal resources [1]. In Japan, researchers initiated a project to investigate the feasibility of geothermal development around and beyond the BDT zone in 2010 (Japan Beyond-Brittle Project, JBBP) [2], and concluded that the development of supercritical geothermal resources can solve most of the problems impeding the development of hydrothermal systems and hot dry rock enhanced geothermal systems (EGS) in the country. After drastic changes to the energy policy induced by the Great East Japan Earthquake on March 11, 2011, the development of supercritical geothermal resources was identified as one of the most promising means to change Japanese energy resources to renewable ones. Currently, NEDO (New Energy and Industrial Technology Development Organization), one of the funding agencies of the Japanese government, is funding studies on the development of supercritical geothermal development, feasibility studies, as well as the drilling of boreholes to provide proof-of-concept.
Researchers in the supercritical geothermal project observed that a considerable number of magmatic intrusions containing a supercritical fluid (supercritical geothermal systems), which have their origins in the subduction of oceanic plates [3], exist in the BDT zone at shallower depths in northeastern Japan. This was done by interpreting resistivity images from magnetotelluric (MT) surveys and the hypocentral distribution of inland earthquakes [4]. The depth of the BDT zone in northeast Japan was identified by drilling a deep exploration borehole WD-1a at Kakkonda geothermal field, Iwate, Japan, in 1998 [5,6].
The temperature profile of WD-1a changed from convectional to conductive in the granitic basement at a depth of 3100 m where the borehole temperature was around 380 • C, associated with penetration into the BDT zone. The depth of the BDT zone depends on the physical characteristics of magmatic intrusion and overburden, and the depth of the BDT zone in northeastern Japan is much less than in other volcanic zones worldwide. This shallow depth of the BDT zone is advantageous for the development of supercritical geothermal systems as it can reduce the cost and risks of drilling to acceptable ranges. Hence, depth estimation of the BDT zone has been identified as one of the most important tasks along with resistivity imaging with MT surveys in the NEDO-funded Japanese supercritical geothermal project.
As thermal profiling using temperature logging data is based on the measured value inside a borehole, it is the most reliable method for estimating the deep thermal structure. However, the temperature logging data cannot be directly used for the identification of the depth of the isotherm corresponding to 380 • C because only a few boreholes have been drilled near and in the BDT zone. Thus, extrapolation and interpolation were required to obtain the thermal profile.
An extrapolation method for the thermal structure using the activity index (AI) was developed [7]. Although the activity index is a technique for evaluating temperature ranks of geothermal fields, it is possible to estimate the temperature profile in deeper areas than the measured values in hydrothermal convection areas. In addition to AI, there are several methods to estimate the subsurface thermal structure in the hydrothermal convection area. For example, Bredehoeft and Papadopulos [8] estimated the subsurface thermal structure extrapolating the model curve from the thermal profile fitting using a one-dimensional heat transport model. Ishitsuka et al. [9] proposed a neural kriging method that estimates the thermal structure based on resistivity data and showed that the thermal structure can be estimated with relatively high accuracy.
However, these methods require the acquisition of detailed temperature profiles and resistivity data. On the other hand, vertical extrapolation of the temperature using the AI requires only the drilling depth and corresponding maximum temperature as inputs, and discharge temperature data for hot springs with no requirement for temperature measurements. Fortunately, 27,261 hot springs are distributed all over Japan [9], and we expected that it will be possible to estimate the subsurface thermal structure based on the appropriate spatial interpolation.
The thermal structure in six major geothermal fields in northeastern Japan was estimated and compared with the hypocentral distribution of natural earthquakes. Suzuki et al. [10] evaluated the validity of the thermal structure inferred using the activity index by comparing the hypocentral distribution of natural earthquakes in two geothermal fields in northeastern Japan. The relationship between the thermal structure estimated from the activity index and the hypocenter distribution is very consistent; however, there are some areas with inconsistencies. The reason is likely related to the spatial interpolation accuracy when estimating the thermal structure. As it is difficult to obtain geothermal data, such as geothermal wells and hot springs at even intervals, certain areas have low spatial interpolation accuracy [11]. Masking such areas is considered necessary.
In this study, we quantitatively evaluated the spatial interpolation accuracy of the thermal structure estimated from the activity index using the geostatistical method. The estimated thermal structure was divided into 1 km grids, and we performed a detailed evaluation by comparing with the bottom of seismogenic layer in each grid. The rest of paper is organized as follows. In Section 2, we provide information regarding the used data. Our methodology is described in Section 3. The results and discussion are presented in Section 4, and Section 5 concludes the paper.

Geothermal Data for Estimation of Subsurface Thermal Structure in Northeastern Japan
The subsurface thermal structure was estimated using the "Geothermal Potential Map in Japan" [12], the "Database by Suzuki et al. [10]", and the thermal profile data of 19 geothermal wells in the Kakkonda geothermal field, provided by Tohoku Sustainable and Renewable Energy Co. Inc., Tohoku, Japan (TOUSEC). The "Geothermal Potential Map in Japan" includes data on 3066 wells and 7203 hot springs distributed throughout Japan. The data on 3066 wells in this database originate from three main data sources; the "Geothermal wells data by the New Energy and Industrial Technology Development Organization (NEDO)" for 572 sampling points, the "Geothermal gradient database of Japan" for 1809 points [13], and the "Database on the Temperature profiles of boreholes in Japan" for 685 points [14]. However, certain data in the "Geothermal wells data by NEDO" and "Geothermal gradient database of Japan" overlap [15], and, for these data, the position information provided by different geodetic systems is incorrect. Therefore, we corrected the position information for these data using the latest geodetic systems. Duplicate data on 115 hot springs where there was overlap of the "Geothermal Potential Map in Japan" and "Database by Suzuki et al. (2014) [10]" were deleted. In the final analysis, we used 1627 data points distributed mainly over northeastern Japan.

Activity Index
The activity index is an index for quantitative evaluation of the potential of a geothermal system or a geothermal well [7]. The activity index is defined by Equation (1).
Tm is the maximum temperature at the observed depth, Tb is the temperature on the boiling point curve for pure water at hydrostatic pressure corresponding to the observed depth [16], and Tg is the temperature corresponding to an average geothermal gradient (30 • C/km) at the observed depth. In this study, to calculate the activity index of hot springs without the temperature profile data, Tm (the maximum temperature at the observed depth) was set to their discharge temperatures, as suggested by Hayashi [17]. In general, the subsurface temperature increases with depth; therefore, Tm is higher than the discharge temperature. Nevertheless, Tm was set equal to the discharge temperature to avoid the activity index overestimating the temperatures in the thermal structure.
The concept of temperature curve estimation using the activity index is shown in Figure 1; the estimated temperature curve is shown by the red dotted line. The activity index was determined as a percentage of the difference between the temperature profiles for AI = 100 and AI = 0. Thus, the temperature at any point can be automatically estimated from the point data for the maximum temperature at the observed depth. Pure water reaches a critical point temperature near a depth of 3.5 km. However, the critical point of water dramatically shifts to higher temperature values with increasing salinity in the natural brine system, and so the AI = 100 curves were simply extrapolated to all the depths for this study, according to Muraoka et al. [15]. The activity index was calculated for all the data available for northeastern Japan, and the depth corresponding to 380 • C was estimated using the activity index. Figure 2 shows the topographic map, the activity index distribution map, and the distribution map showing the isothermal depth for 380 • C estimated using the activity index, for northeastern Japan. Generic mapping tools 5.4.3 [18] was used for drawing with the tension factor at 0.3 and the value smoothed at 2.5 min intervals. Shallow-depth anomalies of the activity index correspond to the presence of quaternary volcanoes and their associated geothermal fields. In this study, therefore, a depth corresponding to 380 • C for northeastern Japan where quaternary volcanoes are present, was estimated to be shallow. A detailed study of the subsurface thermal structure estimated by the activity index was undertaken in six geothermal fields, namely, Mutsu, Hakkoda, Tamagawa, Kakkonda, Kurikoma, and Onikobe.

Kriging
It is not always easy to acquire data for locations that are crucial or ideal for a study, and it can become necessary to estimate the values at points for which no data exists. In such cases, we can quantitatively evaluate the accuracy of the interpolated values using geostatistical methods; kriging is considered one of the most effective methods for this purpose [19]. Kriging is a method of estimating the data value at a point for which data are not available, as a weighted average of the measured values at surrounding points.
First, semivariogram models for the six geothermal fields were created from the data estimated for points at depths corresponding to 380 • C, using the activity index. The model parameters, nugget effect, sill value, and range are shown in Figure 3. The model parameters were obtained from the scatter plot of the experimental variogram using the least-squares method [20]. Anisotropy was not considered when creating the semivariogram models. The depth distribution maps corresponding to 380 • C of six geothermal fields were created by the ordinary kriging method, using the semivariogram models.
The error in the kriging method (standard error) can be quantitatively calculated using only the distance between the points of data estimation and collection. Figure 4 shows the kriging standard error for each field; the estimation accuracy at any location in these fields can be quantitatively evaluated using this.

Kriging
It is not always easy to acquire data for locations that are crucial or ideal for a study, and it can become necessary to estimate the values at points for which no data exists. In such cases, we can quantitatively evaluate the accuracy of the interpolated values using geostatistical methods; kriging is considered one of the most effective methods for this purpose [19]. Kriging is a method of estimating the data value at a point for which data are not available, as a weighted average of the measured values at surrounding points.
First, semivariogram models for the six geothermal fields were created from the data estimated for points at depths corresponding to 380 °C , using the activity index. The model parameters, nugget effect, sill value, and range are shown in Figure 3. The model parameters were obtained from the scatter plot of the experimental variogram using the least-squares method [20]. Anisotropy was not considered when creating the semivariogram models. The depth distribution maps corresponding to 380 °C of six geothermal fields were created by the ordinary kriging method, using the semivariogram models. The error in the kriging method (standard error) can be quantitatively calculated using only the distance between the points of data estimation and collection. Figure 4 shows the kriging standard error for each field; the estimation accuracy at any location in these fields can be quantitatively evaluated using this.   Figure 5 shows the results of a comparison between the depth distribution of the 380 °C isothermal layer and the earthquake hypocenter distribution. The depth distribution maps for the 380 °C isothermal layer were created by applying the kriging method to the depth data corresponding to 380 °C estimated using the activity index. The cross-sectional view shows the 380 °C isothermal layer depth (solid red line) along the survey line (solid black line), with all the earthquakes that occurred within the plan view plotted using the earthquake catalog provided by the Japan Meteorological

Verification of Thermal Structure by Comparison with Seismicity
A strong correlation is known to exist between the bottom of the seismogenic zone and the subsurface thermal structure, due to the transition to the ductile zone [23][24][25]. The BDT is considered to indicate the bottom of the seismogenic zone because the brittle fracture zone causes for an earthquake are controlled by temperature. Hence, to evaluate the validity of the estimated depth of the 380 • C isothermal layer using the activity index, the estimated 380 • C isothermal depth was compared with the bottom of the seismogenic zone, within a column of width of 1 km and depth of 15 km ( Figure 6). Within this column, the number of earthquakes was integrated, starting from the ground surface, and the depth when the number of events reached 90% (D90) was established as a measure of the bottom of the seismogenic zone.

Verification of Thermal Structure by Comparison with Seismicity
A strong correlation is known to exist between the bottom of the seismogenic zone and the subsurface thermal structure, due to the transition to the ductile zone [23][24][25]. The BDT is considered to indicate the bottom of the seismogenic zone because the brittle fracture zone causes for an earthquake are controlled by temperature. Hence, to evaluate the validity of the estimated depth of the 380 °C isothermal layer using the activity index, the estimated 380 °C isothermal depth was compared with the bottom of the seismogenic zone, within a column of width of 1 km and depth of 15 km (Figure 6). Within this column, the number of earthquakes was integrated, starting from the ground surface, and the depth when the number of events reached 90% (D90) was established as a measure of the bottom of the seismogenic zone. However, in this study, there were some points where the number of earthquakes required to calculate D90 were insufficient because the set study range was too narrow. When the number of earthquakes in the column was less than 10 and D90 could not be calculated, the hypocenter depth (D100) at the deepest part of the column was assumed as the lower limit of the seismogenic zone. We defined EDA as the difference between the depth to the bottom of the seismogenic zone thus obtained, and the depth of the 380 °C isothermal layer estimated by the activity index. The EDA was calculated for each of the six study areas. Figure 7 shows the EDA values in each geothermal field. Blank points indicate that no earthquakes occurred in the column. The EDA assumed negative values at many points, being negative when the 380 °C isothermal layer estimated by the activity index was located deeper than the bottom of the seismogenic zone. However, in the Mutsu and Kurikoma geothermal fields, there are fixed points where the EDA showed positive values at the edge of the study area (black open circle in dotted line). These areas had low data density and higher values of the kriging standard error (Figure 4). However, in this study, there were some points where the number of earthquakes required to calculate D90 were insufficient because the set study range was too narrow. When the number of earthquakes in the column was less than 10 and D90 could not be calculated, the hypocenter depth (D100) at the deepest part of the column was assumed as the lower limit of the seismogenic zone. We defined EDA as the difference between the depth to the bottom of the seismogenic zone thus obtained, and the depth of the 380 • C isothermal layer estimated by the activity index. The EDA was calculated for each of the six study areas. Figure 7 shows the EDA values in each geothermal field. Blank points indicate that no earthquakes occurred in the column. The EDA assumed negative values at many points, being negative when the 380 • C isothermal layer estimated by the activity index was located deeper than the bottom of the seismogenic zone. However, in the Mutsu and Kurikoma geothermal fields, there are fixed points where the EDA showed positive values at the edge of the study area (black open circle in dotted line). These areas had low data density and higher values of the kriging standard error (Figure 4). Figure 8 shows the relation between the EDA and the kriging standard error. In all areas, the EDA was within −10 to 4 km. We also found that EDA showed negative values at many points. A negative EDA indicates that the depth corresponding to 380 • C estimated by the activity index was deeper than the seismogenic zone ( Figure 6). This result suggests that the depth corresponding to 380 • C estimated by the activity index was deeper than the actual BDT zone, assuming the bottom of seismogenic zone was equivalent to the actual BDT. However, some points with large kriging standard error at the Mutsu and Kurikoma geothermal fields showed positive values. The data points surrounded by the blue open circles in dotted lines in Figure 8, correspond to those surrounded by black open circles in dotted lines shown in Figure 7. For these points, the interpolation accuracy was low due to the low spatial density of the data points. In addition, the estimation accuracy for the lower limit of the earthquake occurrence layer is related to the number of earthquakes; the occurrence of earthquakes was less than 10 times at all points where EDA showed a positive value.  Figure 8 shows the relation between the EDA and the kriging standard error. In all areas, the EDA was within −10 to 4 km. We also found that EDA showed negative values at many points. A negative EDA indicates that the depth corresponding to 380 °C estimated by the activity index was deeper than the seismogenic zone ( Figure 6). This result suggests that the depth corresponding to 380 °C estimated by the activity index was deeper than the actual BDT zone, assuming the bottom of seismogenic zone was equivalent to the actual BDT. However, some points with large kriging standard error at the Mutsu and Kurikoma geothermal fields showed positive values. The data points surrounded by the blue open circles in dotted lines in Figure 8, correspond to those surrounded by black open circles in dotted lines shown in Figure 7. For these points, the interpolation accuracy was low due to the low spatial density of the data points. In addition, the estimation accuracy for the lower limit of the earthquake occurrence layer is related to the number of earthquakes; the occurrence of earthquakes was less than 10 times at all points where EDA showed a positive value.

Evaluation of the Estimation Method for Subsurface Thermal Structure Using the Activity Index
The possibility that the temperatures in the thermal structure estimated by the activity index can be lower than the actual values can be explained by the conceptual diagram shown in Figure 9. The activity index was regulated by a 30 °C /km average geothermal gradient (AI = 0) and the boiling point curve for pure water (AI = 100). Therefore, the temperature curve estimated for activity index values greater than 0 will show a curved pattern that is convex upward. However, the geothermal gradient in most of the geothermal areas is greater than 30 °C /km, and the boiling point curve shifts to the right with increasing NaCl concentration [16]. Hence, extrapolation of the temperature curve using the activity index is likely to underestimate the values for temperature in active geothermal fields where hydrothermal convection exists. Generally, the geothermal gradient and NaCl concentration of hot water are different in each geothermal area. While these data may improve the estimation accuracy, it is a strength of the activity index that the thermal structure can be estimated from in geothermal areas without such data. The red solid line in Figure 9 shows an example of the temperature profile of a geothermal well drilled in

Evaluation of the Estimation Method for Subsurface Thermal Structure Using the Activity Index
The possibility that the temperatures in the thermal structure estimated by the activity index can be lower than the actual values can be explained by the conceptual diagram shown in Figure 9. The activity index was regulated by a 30 • C/km average geothermal gradient (AI = 0) and the boiling point curve for pure water (AI = 100). Therefore, the temperature curve estimated for activity index values greater than 0 will show a curved pattern that is convex upward. However, the geothermal gradient in most of the geothermal areas is greater than 30 • C/km, and the boiling point curve shifts to the right with increasing NaCl concentration [16]. Hence, extrapolation of the temperature curve using the activity index is likely to underestimate the values for temperature in active geothermal fields where hydrothermal convection exists. Generally, the geothermal gradient and NaCl concentration of hot water are different in each geothermal area. While these data may improve the estimation accuracy, it is a strength of the activity index that the thermal structure can be estimated from in geothermal areas without such data. The red solid line in Figure 9 shows an example of the temperature profile of a geothermal well drilled in the Kakkonda field, showing that the actual temperature curve is more toward the right compared to that estimated by the activity index, and follows an upward convex function.  Generally, the geothermal gradient and NaCl concentration of hot water are different in each geothermal area. While these data may improve the estimation accuracy, it is a strength of the activity index that the thermal structure can be estimated from in geothermal areas without such data. The red solid line in Figure 9 shows an example of the temperature profile of a geothermal well drilled in the Kakkonda field, showing that the actual temperature curve is more toward the right compared to that estimated by the activity index, and follows an upward convex function. Figure 9 suggests that the accuracy for the estimation of the depth corresponding to 380 • C using the activity index is low, with variations of several kilometers. The depth at which the average geothermal gradient (AI = 0) reaches 380 • C is approximately 12.7 km, which is considered the deepest point in an existing geothermal field. The activity index may also exceed 100 due to high concentrations of NaCl in hydrothermal water, but no case where the value exceeds 105 has been reported either in Japan or elsewhere [15,17].
Muraoka et al. [15] assigned a very large activity index to WD-1a, which was calculated based on the highest value in the thermal conduction zone below the BDT. The activity index calculated based on the BDT inflection point was approximately 104. The temperature curve corresponding to an activity index 105 was estimated to have a NaCl concentration of approximately 5% to 10% by weight, as compared with the approximate curve by Haas [16]. When the temperature curve with an activity index of 105 was simply extrapolated to 380 • C, the depth was observed to be about 3 km, which is considered the shallowest depth. Assuming the presence of hydrothermal convection, the range that the actual temperature curve can take is likely to be limited to the range shown in gray in Figure 9. Thus, the deep extrapolation of the temperature curve using the activity index is limited to this range. In this case, the maximum error in the estimation of the depth corresponding to 380 • C using the activity index was estimated to be about 9.7 km.
This result is consistent with the EDA values of -10 to 4 km in all the study ranges ( Figure 8). However, the isothermal layer corresponding to 380 • C inferred using the activity index tended to be deeper than the bottom of the seismogenic zone, so it is estimated that the actual depth for 380 • C is between the depth obtained using the activity index and 3 km, which is the depth corresponding to 380 • C when AI = 105, in most of the geothermal fields. This suggests that the estimation accuracy for depth using the activity index may improve with the increasing value of the activity index.
To evaluate the error of the 380 • C depth estimation method based on the activity index, the depth corresponding to 380 • C of WD-1a in the Kakkonda geothermal field where data are available, and the average value of the depth corresponding to 380 • C, estimated by the activity index of five geothermal wells drilled in the vicinity (<500 m), were compared; the difference was approximately 3.1 km ( Figure 10). Thus, as this value is less than 9.7 km, the estimation is consistent with the above hypothesis ( Figure 9). As the geothermal gradients in actual geothermal areas are often much greater than 30 • C/km, the error in estimating the depth corresponding to 380 • C, using the activity index, is likely to be considerably lower than 9.7 km.  The thermal gradient of WD-1a, as well as the activity index, sharply increased at a depth of around 700 m above the mean sea level (Figure 10). This point is located at the boundary between the shallow reservoir zone and the deep reservoir zone, and the permeability is less in the deep reservoir zone [5]. We observed that this kind of heterogeneity in the convective nature introduces an error in the estimation of the deep thermal structure using the activity index. In the shallow reservoir zone with high permeability, the activity index was low due to mixing, and thus the estimated depths corresponding to 380 • C using data from wells b, c, d, and e, were higher.
We observed that the value of the estimated depth corresponding to 380 • C using the activity index was more likely to be higher than the actual value. However, the EDA was positive in certain areas of the Mutsu, Kakkonda, and Kurikoma geothermal fields. In these areas, the EDA was generally less than 1 km, and the number of earthquakes used to locate the bottom of the seismogenic zone were less than 10; therefore, the estimation of the bottom of the seismogenic zone would be unreliable. In addition, since the accuracy of the seismic depth determination varies within a maximum of 1 km, the depth estimation for 380 • C using the activity index is generally valid for points where the EDA is less than 1 km. However, there are some points where the EDA exceeded 1 km at the edges of the Mutsu and Kurikoma geothermal fields ( Figure 8); these points have large kriging standard errors. Hence, we decided to mask the areas of the map showing depths corresponding to 380 • C ( Figure 11) where the kriging standard error was greater than 1.6 for Mutsu and 2.3 for Kurikoma.   This map can be effectively used to roughly evaluate the depth of magmatic intrusions surrounding existing geothermal fields, although there would be some error in the estimation of the depths corresponding to 380 • C due to the higher geothermal gradients and heterogeneity of structure in the convection zones as discussed in the previous section.

Conclusions
We created activity index distribution maps for northeastern Japan based on existing borehole data. In addition, depth distribution maps corresponding to 380 • C estimated from the activity index were created for six major geothermal fields. The spatial interpolation accuracy between the data of these distribution maps was quantitatively evaluated using the Kriging method. The results showed an anomaly occurring in the relationship between the depth distribution maps corresponding to 380 • C estimated from the activity index and the bottom of the seismogenic layer at the edge of area of Mutsu and Kurikoma geothermal fields. In these areas, the spatial interpolation accuracy was low due to the small amount of data.
In addition, we estimated that the estimation error of the depth of 380 • C determined from the activity index was 9.7 km or less in the areas with higher spatial interpolation accuracy. This estimation was evaluated as valid by comparing the measured value of WD-1a with the 380 • C depth estimated from the activity index. Based on these results, the depth distribution maps of 380 • C estimated from the activity index presented higher accuracy than that of previous research. Thus, we confirmed that the estimation method for deep thermal structures using the activity index was effective for identifying the BDT zone at shallower depths (although with the probability for some error), and for recommending areas suitable for supercritical geothermal exploitation.