Estimation of near-surface temperature in Suwawa Geothermal Prospect, Gorontalo, Sulawesi, Indonesia, based on magnetotelluric and artificial neural network

A geophysical survey using broadband magnetotelluric (MT) technology was carried out in Suwawa Geothermal Prospect Area, Gorontalo Province, Sulawesi Island, Indonesia. The target of that research is to evaluate the geothermal potential hidden below the surface, based on underground resistivity distribution. However, the information about resistivity alone is not enough to get a proper understanding of the geothermal system in this area. Another important subsurface feature that could be useful for the evaluation is temperature. In this study, an attempt to predict the subsurface temperature using resistivity and limited information from a shallow borehole thermogram was carried out. Employing the dependency between resistivity and temperature an indirect temperature estimator was built, thanks to the applicability of artificial neural network (ANN) to learn the pattern connecting both parameters. Comparing some neural network training data shows that the predictive powers of the calibrated neural network highly influenced by the geological difference between the location of borehole and MT station. The best trained ANN was then used to predict the temperature below the other MT stations. The result shows that a proper ANN architecture is important to improve the deeper temperature estimation. The best ANN estimator was obtained from the BT01 and AMT39 data pair, which has the highest correlation as well. This preliminary study gives useful insight into how resistivity could be an alternative tool to delineate the near-surface temperature profile, in order to get a more comprehensive image of the subsurface geothermal system.


Introduction
Geothermal is one of the most promising energy resources preserved in Indonesia, with approximately 23.9 GW of geothermal energy potential [1].The utilization of this energy is highly related to the effort that have been given during the exploration stage. One of the geothermal exploration activities was carried out in Suwawa, Gorontalo, Sulawesi Island, Indonesia by the researchers from the Center of Minerals, Coal, and Geothermal Resources back in 2012 using magnetotelluric (MT) method. The magnetotelluric method is used to image the subsurface condition based on its resistivity value [2][3][4]. However, in order to get a reliable exploration result, one would need other subsurface parameter, such as temperature. Temperature is an important subsurface feature which directly related to  [5,6]. Knowledge about its distribution would be useful to gain a deeper understanding about the geothermal system. Figure 1. AMT stations and shallow borehole location on regional geological map of Suwawa region, Gorontalo province, Indonesia (modified from [7]).
In this study, the attempt to estimate the shallow temperature in Suwawa Geothermal Prospect was carried out. Rather than using direct measurement using temperature logs, temperature profile was predicted from resistivity value below an MT station using artificial neural network (ANN). Resistivity is a temperature-dependent parameter [8][9][10][11], therefore predicting the temperature from resistivity is possible if the relationship between both parameters was recognized [12]. The resistivity data was obtained from MT sounding, while a limited number of temperature data was collected from a shallow borehole thermogram. The degree of correlation between temperature and resistivity information, and the reliability of the ANN technique was firstly evaluated before finally being used to estimate the temperature below MT stations nearby.

Field Data
The geothermal prospect is situated in between 0°28'13.7" N -0°36'54.8'' N and 123°06'00'' E -123°15'00" E. Lithological condition of Suwawa Geothermal Prospect area consists of tertiary to quaternary rocks. The oldest rock is Bilungala Volcanic which includes breccias, tuffs, and andesitic lavas from Early to Late Miocene ages. Diorite Bone, located in the middle of the area, consists of quartz diorite, diorite, granodiorite, and late Miocene granite. In the south-eastern part, there is Pinogoe volcanic group which comprises of tuffs, lapilli tuffs, breccias, and lava from Pliocene to Plistocene ages. The youngest rocks are lake sediment and alluvium [7].
There are two kinds of data sets that were used in this preliminary work. These data were collected from Suwawa geothermal prospect area in Gorontalo Province, Indonesia. One data set needed is shallow borehole thermogram data, and the other is audio-magnetotelluric data. This field has only two borehole thermogram data recorded up to this time: BT01 and BT02. However, the latter does not have any AMT station nearby, so it cannot be used. For BT01, there are 3 surrounding AMT stations: AMT16, AMT17, and AMT39. The location of the borehole point and all AMT stations are shown in figure 1. As seen on the map, these AMT points are situated in different geological conditions and at different distances relative to the exploratory well. The AMT16 station is located on an alluvium layer about 900 m from the borehole. In the same surface geological condition, AMT39 lies 680 m to the west from the borehole. The last station, AMT17 is located at around 700 m from the borehole and sits on an andesitic body. This different condition might affect the correlation between resistivity derived from AMT and borehole temperature value.

Borehole thermogram data
Borehole BT01 is located on an alluvial overburden. The well log of BT01 is shown in figure 2. From the drilling results, it is known that lithology below the well site consists of alluvial deposits (Al) from surface to 34 m depth, unaltered to weakly altered polymictic breccias (BX) at a depth of 34 -120 m, and a moderate to very strong altered polymictic breccia layer (BXT) at a depth of 120 -250 m. Subsurface temperature was measured up to 250 m deep at 2 m intervals. The profile shows that the temperature is mainly increasing with the depth. In parallel with the heat measurement, investigation of lithology, as well as alteration intensity, was also carried out. The existence of the altered breccia also shows a high alteration at the respective depth, which is caused by mineralization of geothermal fluids. This argillic type of alteration rock functions as the clay cap of the geothermal system below. The increasing temperature verifies that this borehole penetrates the up-flow zone [13].

Audio-magnetotelluric data
Three AMT impedance data sets from the nearest sites from borehole BT01 were used for calculating the relation between resistivity and subsurface temperature. These data sets were actually derived from a combined AMT-MT survey. The frequency range of the signal is recorded from a high of 8,192 Hz to a low of 0.117 Hz in average. These data sets were collected by using Zonge Geophysical Receiver GDP-32II. For balancing the comparison, the electromagnetic data that was used for correlation was only the high-frequency part (AMT), which also has better resolution for the shallower part [14,15]. The AMT responses (apparent resistivity and impedance phase) at each station are shown in figure 3.  Figure 3. AMT responses, apparent resistivity (upper part) and phase (lower part) at AMT16, AMT17, and AMT39.
From both measurement modes (TE and TM), the average one-dimensional sounding profile is calculated as the determinant. In this work, the audio-magnetotelluric sounding responses, resistivityphase versus frequency, were inverted by using Bostick transformation for a quick calculation. This simple inversion scheme produces an approximate curve of resistivity versus depth [16,17]. Resistivity values at every 10-meter depth were then calculated by employing the linear interpolation technique. Tens of data pairs of resistivity and temperature at the same depth were examined before finally being used as the training basis.

Dependency analysis
An attempt was made to determine the relationship between resistivity and temperature at the same depth below the surface. Bostick resistivity values from all AMT stations were first plotted against temperature at the same depth, then the dependency between both variables was approximated by using several types of regression formulae to see which one was best for this case. The quality of fitness was analysed by the coefficient of determination value (R-squared). This coefficient gives information about how well a regression equation represents the data distribution. The R-squared value for each regression equation is listed in Table 1. The value itself is the ratio between explainable and unexplainable data values in a set of data pairs. Based on the regression result, it can be roughly said that temperature was more strongly correlated to resistivity in an exponential equation than in other regression types. Due to low R-squared value, the result is not so satisfying. Nonetheless, this kind of analysis remains beneficial to explain how the variables are related one to another.

Effect of Resistivity-Temperature Dependency
By using some training data, which was extracted from the shallower part of the data pairs, we tried to estimate the rest. This analysis was carried out to evaluate the effect of correlation value obtained in the previous section, to the prediction accuracy. The result is shown in figure 4 along with the observed temperature value for comparison, with the root mean square error (RMSE) value.
From this result, the calculated temperature profile produced from the BT01 and AMT39 pair has the best fit among all the results, as it has the smallest RMSE of the samples. This result also gives information about how strongly the geological condition and distance affects the relation between resistivity and temperature. Having the least error, the BT01 and AMT39 data pair can be used in the network training phase, before being employed to estimate the deeper temperature values below the surrounding AMT stations. However, the predicted values from the BT01 and AMT16 data pair does not give poor result. It could still be utilized, but with a lower degree of confidence. On the other hand, the BT01 and AMT17 data pair gives quite rough predictions even visually. Therefore, it is unwise to choose the pair as the neural network basis. . Measured (dashed line) and calculated (solid line) temperature profiles from three different data pairs. The result of BT01-AMT39 has the better prediction than the others, as this pair has the highest R-squared value as well.

Result and Discussion
From the comparison result, the deeper portion of temperature values was predicted using the transfer function based on the ANN training of BT01 and AMT39 data pair. With this transfer function, the temperature below the additional AMT station AMT15, AMT16 and below the AMT39 itself were then tried to estimate. The AMT15 station lies 1 km northeast from AMT16 as shown earlier in the map earlier in figure 1. Although it is unfair to forecast the temperature below the other AMT stations with only one training data using the proposed ANN model (which roughly means that the effect of relative distance from the AMT station to the heat source is neglected), the estimation was still carried out in order to test the impact of resistivity differences between the AMT stations. The results are shown in figure 5.  Figure 5. Predicted temperature profile after temperature forecasting under AMT sites nearby, and its updated correlation coefficient with the Bostick resistivity value.
Due to the limited training data, this prediction can only be performed for a limited depth, up to 600 m, to avoid over-estimated values. As shown in the results, the estimated temperature values are different for each resistivity profile, even though the differences are small compared with the resistivity changes. The correlation coefficient was also recalculated. The numbers mainly show that the resistivity and the estimated temperature are associated. In this result, the temperature below the AMT39 station is assumed to be the same as the temperature measured from borehole BT01, because the resistivity from this point is calibrated with the temperature from the borehole. Meanwhile, the temperature values in degree Celsius below the AMT15 and AMT16 stations are calculated from the transfer function of the calibrated ANN, assuming that they are lying at the nearly same geological condition and at about equal relative distance to the heat source with AMT39. Predictions below other points that lie in a different geological condition were not carried out because of the probability that the transfer function is only valid for one certain environment, although this behaviour needs to be further studied.
From the estimated temperature value, regardless of the large difference in the AMT resistivity responses between all stations, the calculated temperature values have only slight variation. From this result, the estimated temperature is increasing with depth at a rate of approximately 70 °C/km of geothermal gradient. However, the temperature profile is not in linear form like the constant geothermal gradient due to the effect of subsurface resistivity values. This estimated gradient is reasonable for a geothermal gradient value in a geothermal field for near-surface temperature.
Then, a minimum curvature gridding technique was used to interpolate the grids between AMT15, AMT16, and AMT39 in order to test the appropriateness of representing subsurface temperature distribution in a two-dimensional vertical cross section in this introductory study. Due to the limited number of data, the horizontal grid is minimized to prevent error and over-estimating value. The result is quite satisfying and plausible as shown in figure 6.
In the shallow and middle parts, subsurface temperature is distributed almost uniformly. While, in the deeper zone, the contour starts bending with the temperature below AMT16, and is slightly higher than its surrounding. The minimum curvature method minimizes the bends of the interpolated grid. So, the curve seen in the graph mainly comes from the actual data. Despite the limited displayed data in this study, it is possible and promising to create a larger temperature distribution map, if some additional borehole data are at hand in the future.  Figure 6. Predicted temperature distribution below AMT15, AMT16, and AMT39 up to 600 m based on the result of artificial neural network using BT01 and AMT39.

Conclusions
In the subsurface, resistivity and temperature at certain depths are most likely correlated exponentially from the coefficient of determination analysis, despite the unsatisfactory value due to the far distance between borehole and AMT points. In order to estimate subsurface temperature, the ANN technique was applied using resistivity data from the AMT method and temperature data from a borehole. A reasonable ANN model is needed to sharpen the temperature estimation for the deeper subsurface. In this case, prior knowledge about the factors that affect the subsurface resistivity value is needed. As a result, deeper temperatures were able to be estimated by resistivity and depth information. From the comparison of ANN test results, the best temperature estimation was obtained from the BT01 and AMT39 data pair, which has the highest correlation coefficient as well. This shows that the geological condition, as well as distance, should be wisely considered. In other word, AMT-Temperature data pairs that come from the same geological location and/or close in distance would be a better ANN estimator training basis. Therefore, the data pair that will be used in this indirect geothermometer method should be chosen carefully.
The predicted temperature profiles below each AMT stations are slightly different, regardless of the large difference in resistivity values. It shows that the influence of resistivity in calculating the temperature value is small compared with the effect of depth or other influential factors. A twodimensional cross section of subsurface temperature is possible and could clearly describe the temperature distribution below the AMT stations. Nevertheless, more training data necessary for obtaining a better training result. So, deeper and larger sets of borehole data will surely provide a better result.