GEOTHERMAL RESERVOIR DEPTH OF SEULAWAH AGAM VOLCANO ESTIMATED FROM 1D MAGNETOTELLURIC

1Geophysical Engineering Department, UniversitasSyiah Kuala, Banda Aceh, Indonesia 2Geological Engineering Department, Bandung Institute of Technology, Indonesia 3Chemistry Department, UniversitasSyiah Kuala, Banda Aceh, Indonesia 4Geological Engineering, UniversitasSyiah Kuala, Banda Aceh, Indonesia 5Energy and Mineral Resources Agency of Aceh Province, Banda Aceh, Indonesia 6Physics Department, University of Brawijaya, Malang, Indonesia 7PT. ElnusaTbk, Jakarta Selatan, Indonesia 8Electrical Engineering Department, UniversitasSyiah Kuala, Banda Aceh, Indonesia 9Geophysical Engineering, Sumatera Institut of Technology, Lampung Selatan, Indonesia * marwan.geo@unsyiah.ac.id


INTRODUCTION
Geothermal energy is a renewable energy source that is also dependable and environmentally friendly, and its use may reduce our reliance on horrifying fossil fuels [1]. Geothermal energy could be generated and stored inside volcanoes as hot water or steam [2], [3]. In some places, volcanoes are used for power plants andgeotourism, which might improve the community's economy [4], [5]. Recently, Indonesia has targeted 23% of its total consumed energy from renewable resources by 2025, increasing to 31% by 2050 [6]. The target relates to Indonesia's geothermal potential of 28,000 MW [7]. Currently, 12 power plants have been built in Indonesia, with Kamojang developing 235 MW, Sarulla developing 330 MW [8], and WayangWindu developing 227 MW [8], [9]. While there are 20 geothermal sources in the northernmost region of Sumatra, Aceh Province, there are no geothermal power plants [10]. Four of the geothermal sources are activestatus, i.e., Jaboi volcano which has an estimated energy capacity of 50 MW [11], PeutSagou volcano with 100 MW [12], [13], geuredong volcano [14], and the largest volcano SeulawahAgam with an energy of 275 MW [15], [16]. Figure 1 shows the topography of SeulawahAgam which is visualized by the Shuttle Radar Topography Mission (SRTM) with a resolution of 30 m/px [17]. Various explorations have been carried out in SeulawahAgam to study a geothermal system, such as those using remote sensing, which can provide an initial picture of the manifestation distribution based on the vegetation and temperature obtained from Operational Land Imager (OLI) and Thermal Infra-Red Scanner (TIRS) sensors [18], [19]. Furthermore, the studies were also carried out in gas and geothermal water surveys [16], [20] and integrated geophysical surveys, i.e., magnetotelluric [21], [22]. In addition, an electromagnetic transient was also conducted to measure the crossing of SelawahAgam volcano to obtain geothermal system information, such as heat sources and reservoirs [15]. These surveys, as previously mentioned, were still limited to certain areas and did not cross the volcano's regional fault, so they do not provide clear information on hydrothermal and fault distribution as the primary regulators of volcanic fluids. As a result, relatively broad data distribution is required to cover the entire volcano. Geothermal power sources are ideal for measuring targets using electromagnetic methods, as they produce substantial variations in subsurface electrical resistivity [23], [24], wherein thermal areas, electrical resistivity is substantially different and generally lower than in areas with cooler subsurface temperatures. In terms of temperature difference, anomaly targets could be easily distinguished from the surrounding area that contains rocks [25], [26]. Because of the ability to map the deep conduction characteristics that play an important role in geothermal energy, magnetotellurics is particularly effective in defining subsurface geology among various electromagnetic methods [23], [27]. The electrical resistivity of geothermal systems depends on factors, such as temperature, porosity and permeability, fluid salinity, and mineralogy alteration [28]. The magnetotelluric (MT) survey provides an understanding of shallow to deep geological features. Because of the sensitivity of electrical conductivity and its depth investigation, MT is an important tool for geothermal exploration [27], [29]. Furthermore, the MT method has a much deeper investigation depth than other electromagnetic (EM) methods, which are typically incapable of determining geological features and detecting geothermal reservoirs more significant than 1-2 km. Therefore, to add more geophysical surveys on the SeulawahAgam volcano, we analyzed magnetotelluric data with a wide volcano distribution to investigate the hydrothermal system and faults.

BACKGROUND AND GEOLOGICAL SETTING
The distribution of active volcanoes on Sumatra Island is distributed across the volcano. It is generally parallel to the subduction zone, with depths ranging from 100 to 150 kilometers from the subduction plate [30], [31]. The tilted convergence of the Indo-Indo-Australian and Eurasian plate collisions formed a dextrose slip known as the Great Sumatran Fault. Tectonically, the Subduction Zone of Sunda Megathrust in the Western part of Sumatra produces the GSF fault split into 20 segments [31]- [34], wherein the northern part of the GSF is divided into two branches, namely the Aceh segment that leads to Aceh Island and the Seulimeum segment that leads to Weh Island [35]- [38]. Simultaneously, the Seulimeum fault breaks through the SeulawahAgam volcano to the west [39], [40]. The Seulawahvolcano is dominated by Lamteuba volcanic rocks composed of basaltic andesite composed lava and pyroclastic flows, according to the geology map reported by [41], which is dominated by andesite-dacite, sandstone, volcanic breccia, tufa, and agglomerate are the constituents of this formation [41], [42]. The volcanic activity is characterized by several traces of manifestation indicating the presence of geothermal energy in a certain geological condition of subsurface rocks [43], [44]. SeulawahAgam volcano has seven manifestations, two of which are craters, namely Heutz crater on the Northside and Ceumpaga Crater on the South of the volcano [15], [29]. Others, such as IeJue and IeSuum in Krueng Raya, Aceh Besar, appear in the form of hot springs and warm ground. Hot sulfuric hot springs are found in these manifestations, and the water's color is murky due to mixed acidic clay [16].

BASIC THEORY OF MAGNETOTELLURIC
Magnetotellurics is a passive EM method that involves fluctuations in electric fields and natural magnetic fields perpendicular to each other on the earth's surface to determine the value of resistivity in the earth's subsurface up to tens of kilometers away [45], [46]. The frequency used in the MT method is 10 −4 -10 4 Hz, allowing us to map the distribution of rock resistivity values with high depth [47]. The MT field is generated by a broad spectrum of EM waves that occur naturally due to lightning discharge and interaction between the solar wind and the earth's magnetosphere [48]. The source wavefield is a polarized plane wave that propagates vertically to the earth's surface. Mathematically, it could be calculated using Eq.1: Where is the resistivity and is the frequency. The basic principle of electromagnetic theory is Maxwell's equations that explain the properties of electricity and magnetism and their interactions [49]- [51]. The comparison of electric and measurable magnetic fields is commonly referred to as impedance when obtaining information about the resistivity structure of measurements on the earth's surface. The magnetotelluric method has been used to calculate tensor impedance (Z) as a function of frequency, as shown in eq.2 [52].
where (Ex, Ey) is the orthogonal horizontal electric field component and (Hx, Hy) is the magnetic field at a specific frequency. The impedance tensor consists of four sensor elements, each of which is a complex number represented by a matrix. Every component of the impendence tensor has a value of resistivity and phase. While the apparent resistivity (ρa) can be described in the following Eq. 3: Where the phase (Φ) can be described by Eq. 4: Where ω is frequency angular = 2πf. To estimate variations in resistance values to depth, pseudo-resistivity calculations in frequency domains could be used. When the frequency decreases, the magnetotelluric wave's penetration range expands and reaches a greater depth. Depending on the variation in resistance to time, phase changes could be either positive or negative.

MEASUREMENTS AND DATA ANALYSIS
To cover the entire area of SeulawahAgam volcano, MT data measurements were taken into 4 profiles with a total of 26 points crossing the volcano, surrounding geothermal manifestations. Figure 2 shows the measurement design and data distribution of MT.This data distribution was chosen to represent the volcanic area in general. Besides that, the measurement profile is made across the volcanic area, and manifestations that can provide subsurface images of the hydrothermal system of the volcano, the MT profile in Figure 2 is shown by P1 -P4. The distance between the measurement points was made 500 m to acquire a detailed picture of the volcano that could vertically form a cross-section of four profiles.Track one has five measurement points, a total track length of 9 km, and distance variations ranging from 2 to 3 km. Track two has eight measurement points, a total track length of 13 km, and distance variations ranging from 1.5 to 2.5 km. Then there are eight measurement points on track three, with a total track length of 15.4 km and distance variations ranging from 2 to 4 km. There are five measurement points on track four with a total track length of 11 km and variations in distances ranging from 2 -4 km. For each point of measurement, it required 12 hours to record magnetic and electrical field data from the MT method. Magnetotelluric data collected in the field vary in the electric and magnetic field intensities over time. As a result, some standard corrections were required, including time series analysis, robust processing, and data correction with SSMT2000 as basic software from Phoenix Geosystem for data-processing in MT. The data sorting stage was carried out because of interference from vibration and electrical activities such as electricity near the measurement point [53]. The data points from each frequency were averaged resulting in the pseudo resistance curve and impedance phase of each point. Curves with many error bars that are difficult to smoothen could be managed by increasing cross power to produce smooth MT curves in *EDI format that act as MT data standards.

RESULT AND DISCUSSION
Before performing the inversion process, tipper analysis might provide information on subsurface structures [54], where the magnitude of the induction arrow is proportional to the anomalous intensity of the current concentration, which is associated with the magnitude of the conductivity gradient. As a result, the tipper arrow's direction may indicate a conductor location that is distant from the resistive zone [55]. We used a tipper analysis of the entire MT data in this study to determine the location of the conductor's suspected subsurface anomaly, as shown in Figure  3. Arrow symbols suggest the induction magnitude is proportional to the anomaly intensity of current concentration [56]. This is related to the magnitude of the thermal conductivity gradient on the site. In other words, the direction of the tipper heat indicates the presence of a conductor. The tipper analysis was carried out by plotting all measurement points at four different frequencies: high frequencies (86 Hz and 5.6 Hz) were used for shallow anomalies. In contrast, low frequencies (2.3 Hz and 0.05 Hz) were used for deeper anomalies subsurface. At high frequencies (86 Hz) as shown in Fig 3.a, noise on the surface or shallow depth causes a relatively small gradient induction in the entire region. In general, the anomaly leads to a conductor close to the manifestation. However, this anomaly is still random in some places caused by signal interference at frequencies near the surface. Similar results were given by the investigation using a frequency of 5.6 Hz (Fig.3b). At frequencies (2.34 Hz and 0.05 Hz) as shown in Fig.3 c and d, the induction directions tend to be at the point of manifestation, such as the Heutz crater, IeJue, and IeBusuk on the volcano's Northside. A high gradient of induction distinguishes the campaign manifestation on the Southside of SeulawahAgam. It is triggered by conductive zones such as reservoirs, manifestation traces, and subsurface fault layers.

Rho and Phase Data
A comparison of pseudo resistivity and phase against frequency was performed in the data analysis stage, where for transverse electric (TE) mode, data were taken from the pseudo resistivity and phase in xy direction toward the electrical field vector aligned with the measuring track direction and magnetic field (that is perpendicular with the track direction). In the Transverse Magnetic (TM) mode, the data were derived from the pseudo resistivity and phase in the yx direction calculated using a magnetic field vector aligned with the measuring track and an electrical field perpendicular to the measuring track. Data analysis was conducted on the entire measuring points of MT, which specifically, we only showed several examples of rho and phase for xy and yx ( Figure 4).The response obtained for the three locations was relatively equal for rho or phase data, such as at site X1 (Figure 4a), which is located near IeJue and IeBusuk craters and has rho at high xy and yx frequencies, implying resistive layers. A low rho was obtained as a response to the conductive layer at the median frequencies.
Meanwhile, the rho value returned high at low frequency with a relative pattern against the first layer. Other than resistivity data, phase values could indicate the general condition of the subsurface, where values of conductive layers could be indicated by low rho and high phase. The phase values on MT data are classified into three categories: phase 0°-45° indicating the alteration of subsurface resistive anomaly, phase 45°-90°indicating conductive anomaly characteristics, and phase 45°indicating a response against lateral layers. Other locations, such as X2 ( Figure 4b) and X3 (Figure 4c), located close to Hertz crater, also exhibited similar characteristics to location X1. The rho and phase data suggested three subsurface layers: the first was dominated by a resistive zone estimated as a response to rocks at the bottom, the second by a conductive zone suspected as a response to a reservoir anomaly investigated in geothermal studies, and the last layer had resistive characteristics inferred as a response to a heat source anomaly beneath the volcano. Besides, the response of rho and phase data that is relatively equal on station MT could be utilized as quality control. The similarity to the response could indicate that the responding data from the subsurface anomaly was of high quality. In general, the responses shown by MT data from the SeulawahAgam volcano are consistent with resistive characteristics reported by others in other locations [57].

One Dimensional Inversion
1Dinversion process was conducted using IPI2Win MT software developed by Moscow State University. This method solved the inversion problem by utilizing variance from the Newton algorithm based on the least number of layers or the Tikhonov approach. IPI2Win toward three modes could perform inversion: mode XY (TE), mode yx (TM), and invariant as a combination of TE and TM. The comparison process of the calculation model using the observation was conducted against the entire MT points, where the 1D model could represent horizontally layered homogenous conditions. The 1D Inversion would show the sounding curve of the earth's resistivity variation as a function of depth. All MT stations were inverted 1D with five rho and thickness layer data differences, and 14 iterations were conducted to obtain a small RMS value ranging from 0.8% to 2.54%. Specifically, Figure 5 shows 1D inversion results for TE and TM modes carried out at locations X1 close to the IeJue manifestation, X2-X4 near Cempaga crater, and X5-X6 in Heutz crater. At location X1 (Figure 5a) in the first profile (12 km from IeJue manifestation), a high response of rho at 0-1 km depth as a representation against the resistive layer was obtained; this anomaly is similar to that between TE and TM. Meanwhile, at 1-2 km depth, a slightly different response was obtained, with TM being more resistive than TE. At the 3 to 5 km depth, the value of rho from TE was higher than that of TM. At other locations, such as X2 ( Figure  5b) and X4 (Figure 5d) that were close to Cempaga crater, rho values were obtained similar, in which at the 0-1 km depth the values were dominated by low resistivity, while at a depth of 1-2 km-dominated by high resistivity. According to the TM data, the rho value from both locations was different, with the TM rho from X2 being relatively similar to the TE rho at X2 and X4. These differences were attributed to response sensitivity from TM mode being more dominant against magnetic parameters of an anomaly. In contrast, TE mode is more dominant against the electrical parameters of an anomaly. At the deep depth (2-5 km) both stations indicated a similar response with a low rho value due to subsurface conductive anomaly. At the Heutz crater location, where points X5 ( Figure 5e) and X6 (Figure 5f) were located, the rho values obtained from TE mode were equal for both sites, in which the resistive zone dominated the first layer at a depth ranging from 0 to 400 m, The characteristics of the conductive zone up to 1 km depth are followed by those of the resistive zone at 5 km depth. Meanwhile, rho data from TM may show a similar pattern for both data sets. Following that, we combined the overall results of the 1D Inversion to develop a cross-section that crosses several geological features such as faults, manifestations, and volcanoes.

2D Cross-Section from 1D Inversion
The constructed 2D model is an interpolation result of 1D Inversion at each measuring point in a tract correlated and displayed in distance area and elevation. In the 1D inversion model, the overall track of the chosen tract was in yx mode due to its good sensitivity in vertically mapping the subsurface resistivity value. A geothermal system could be represented as fluidal convection trapped in heat transfer from the heat source to the heat reservoir on a free surface on the earth's crust. The subsurface magma formed the geothermal system due to plate collisions, which causes the subsurface rocks to experience conductive heat. A geothermal system comprises five primary components: a heat source, a reservoir, a clay cap, a geological structure in faults, and a heat-carrying fluid. Specifically, cross-sections 2 D profile 1-4 were obtained from the 1D inversion results against all MT stations ( Figure 6).
In profile 1 (Figure 6a), a low resistivity was obtained within the range of 5.8-20 Ωm up to 2 km depth thickness to the east direction and thinning at the center part as a representation of conductive anomaly. This anomaly is also thought to be a clay cap that has risen to the surface and covered the reservoir zone. Furthermore, this clay cap layer is recommended to preserve heat fluid accumulation in reservoir rocks located beneath it. The clay contents resulting from the alteration of rocks with accumulated geothermal fluids are characterized by a resistivity of less than 10 Ωm. Next, the resistivity zone was obtained within 42-98 Ωm found in the central part of the measuring track. This layer is suspected to be a reservoir layer functioning as a place for the heat fluid layer filling the reservoir rocks. The increase in resistivity beneath the conductive clay cap layer indicates an elevated temperature in the deep depth, which is consistent with a high-temperature geothermal system [44]. The zone of high resistivity observed on the volcano's east side in 1.5-5 km subsurface depth was suspected to respond to the heat source zone becoming a fluid source in reservoir areas.
In profile 2 (Figure 6b), a response model was obtained that is relatively similar, especially with the conductive zone generally observed in the depth of 1-2 km. This shallow layer was observed to have low resistivity from both tracks, which is consistent with the field condition. The area is near IeJue manifestation in fumaroles and hot springs, and it is traversed by several faults that primarily regulate the manifestation system. Meanwhile, high resistivity was obtained at depths ranging from 2 to 5 km and was assigned as the heat source. From tracts P1 and P2, it was possible to see the presence of a resistive anomaly passing through the surface, which was suspected to be caused by magmatic instruction of the volcanic regions. The rocks are estimated to be located at a depth of 3-5.5 km, assigned as the heat source zone of the SeulawahAgam volcano. The heat source is included as an essential component in a geothermal system, characterized by a significantly high resistivity. In profile 3 ( Figure 6c) and profile 4 (Figure 6d), the top layer controlled a moderate resistivity value and was assigned as the topsoil layer affected by volcanic sediment. For shallow depth at P3, the resistivity value obtained was ranged between 52 and 75.7 Ωm at a 12-15-km distance from the measuring point. This anomaly could be due to the presence of local structure as the primary fluid regulator rising to the surface and colliding with Sumatra regional faults. A similar correlation was observed at P4, where high resistivity was observed between 9 and 11 km, indicating the presence of local structure, though no near-surface manifestation was observed. Fluids directed to the surface were suspected of being trapped in an impermeable zone or a mere steaming ground.  Thus, there was no manifestation emerging to the surface. At the following layer with a thickness of 1 km, the conductivity was higher than the above layer. Because of its low resistivity, this layer was possibly assigned as a clay cap layer and stretched from 0 km to 11 km. Tracts three and four revealed that both models were traversed by minor faults associated with the Sumatra fault and followed a similar pattern.This is because a geological faulting region causes significant fracturing, resulting in a permeable subsurface zone, which acts as a channel for geothermal fluids to flow from the reservoir to the surface [11], [24]. The pseudo-2D modeling from 1D input of MT data confirmed the subsurface geological condition through the generated resistivity variation. The resistivity model depicts the subsurface material as well as existing manifestations at a depth of 5 km. Furthermore, the inversion results using data from the 1D inversion model produce a better result as well as the ability to control the subsequent inversion process.

CONCLUSION
To study the geoelectric structure and volcanic system related to the geothermal SeulawahAgam volcano, we applied the MT method spread over 26 points that intersect several manifestations such as the Heutz crater, Cempaga,IeJue and IeSuum hot springs. This data also crosses several local faults that become the central controller of the volcanic manifestation system. Several data analyses such as tipper analysis, 1D inversions of the transverse electric (TE) and transverse magnetic (TM), and pseudo-2D cross-section of 1D models in four profiles have been performed. The results of the cross-section model show that the layer with volcanic deposits has a resistivity value ranging from 57 -98 m and is present in the entire measurement trajectory. In contrast, the clap clay layer that covers the fluid rising to the surface can be found at a depth of 2 km. Furthermore, a reservoir layer with a resistivity value of 94 -188 m is obtained at a 2-3 km depth. The area under the reservoir is estimated to be a heat source generally being at a depth of 3-5 km, and this anomaly is characterized by high resistivity values >1000 m. The results obtained from the 1D model can provide an initial description of the characteristics of the SeulawahAgam geothermal system. Still, the obtained model needs to be verified with 2D and 3D inversions, which are known to map 2D rock anomalies below the surface.