Mapping of Fault and Hydrothermal System beneath the Seulawah Volcano Inferred from a Magnetotellurics Structure

: Magnetotellurics (MT) is an important geophysical method for exploring geothermal systems, with the Earth resistivity obtained from the MT method proving to be useful for the hydrothermal imaging changes of the system. In this research, we applied the MT method to map the geothermal system of the Seulawah Agam volcano in northern Sumatra, a site intended for the construction of a geothermal power plant with an estimated energy of 230 Mwe. Herein, 3D MT measurements were carried out, covering the entire area of the volcano and the various intersecting local faults from the Seulimeum segment in the NW–SE direction. Based on Occam 2D inversion, a conductive anomaly (<10 ohm · m) near the surface was identiﬁed in response to speciﬁc manifestation areas, including the Heutsz crater on the northern side and the Cempaga crater on the southern side. A further conductive anomaly was also found at a depth of 1 km, which was presumably due to a clay cap layer covering the ﬂuid in the reservoir layer below the surface, where the manifestation areas are formed at various locations (where faults and fractures are found) owing to the ﬂuid in the reservoir rising to the surface. The MT modeling also revealed that the reservoir layer in Seulawah Agam lies at a depth of 2 km with a higher resistivity of 40–150 ohm · m, which is the main target of geothermal energy exploration. At the same time, the heat source zone where magma is located was estimated to lie in two locations, namely, on the northern side centering on the Heutsz crater area and the southern side in the Cempaga crater area. A clear 3D structure obtained via Occam inversion was also used to visualize the hydrothermal ﬂow in the Seulawah Agam volcano that originates from two heat source zones, where one structure that was consistent across all models is the conductive zone that reaches a depth of 5 km in the south in response to the regional faulting of the Seulimeum segment. Based on the MT research, we concluded that the volcano has the geothermal potential to be tapped into power plant energy in the future.


Introduction
Geothermal energy presents alternative energy generated from the heat sources deep within the Earth's crust and is generally associated with the volcanic and tectonic activities the volcano is located in Aceh Besar, 85 km from Banda Aceh, while (c) is the topography of the volcano that located in the area close to the Great Sumatran Fault (GSF). P1 is present the Profile 1 of 2D MT measurement evaluated in 2018 [36], while P2 is Profile 2 that was published in 2019 [21].

Seulawah Agam Geothermal Field
Mt. Seulawah Agam is a stratovolcano located in Aceh Besar District, Aceh Province, Indonesia. Currently, a plan to construct a 230-MWe geothermal power plant is underway. In general, the formation system of all volcanoes on the island of Sumatra is influenced by the Great Sumatra Fault's (GSF) activity that stretches some 1700 km from Lampung to the Andaman Islands of India [37][38][39]. In fact, the GSF is not an intact fault block but a collection of 20 segments dispersed across the island of Sumatra [40,41], where the Seulawah Agam volcano lies close to the Seulimeum segment leading to Pulau Weh [37,41]. In addition, several local faults with an NW-SE direction were identified as the main controllers of the volcano [29,42,43] and the sites where the fluid flows above the surface, which means that the fault location is closely related to the presence of various manifestations, as shown in Figure 2.
Various geological maps [42,43] have revealed that Seulawah Agam volcano is dominated by Lam Teuba volcanic rocks consisting of andesite to dacite rocks, pumice breccias, tuff, agglomerates, and dust flows of the Pleistocene age, while during the Holocene, lahar deposits settled on the older Lam Teuba volcanic rocks. The volcano also features various other formations, including the Seulimeum, Kota Bakti, and Idi formations, which are located on the western and southern sides. Meanwhile, the volcanic activity of Seulawah Agam can be identified from various mapped surface manifestations, including one in the south that is almost adjacent to the volcanic peak, namely, the Cempaga crater, and another in the north, namely, the Heutsz crater. In addition, this section also features various hot springs, including those in the Ie Suum and Ie Jue areas. The manifestations in the Ie Jue area are characterized by the presence of warm ground, steaming ground, hot springs, fumaroles, and mud puddles, which are distributed on the northern and southern sides of the volcano. A geothermometer study demonstrated that the surface water temperature of the Ie Jue manifestation is 97.61 • C ± 0.058, while that of the Heutz crater area is 76.6 • C ± 0.08, as indicated by inductively coupled plasma (ICP) analysis [28].  [42]; (B) photo documentation of the Heutsz crater emitting sulfur oxide gas in an area close to the top of volcanic with an altitude of 700 m asl; and (C) documentation of a mud pool with a temperature of 89.05 • C [28]. At the Heutz crater area, some neo-formational minerals were also found, such as sulfur, sulfatara, mudpole, and fumarole, indicating the volcanic activities of Seulawah Agam Mountain.

Overview of the Magnetotellurics Method
The MT method is a passive electromagnetic method that involves measuring fluctuations in the natural electric (E) and magnetic (H) fields in orthogonal directions on the earth's surface to determine the distribution of conductivity at depths ranging from a few meters to several hundred kilometers [44,45]. The natural source of the MT field around 1 Hz is thunderstorms, while most of the signal below 1 Hz is caused by magnetospheric currents created by solar activity [46]. In the frequency domain, the measured horizontal components of the electric field E = E x E y T and the magnetic field H = H x H y T are related in terms of an impedance tensor, expressed as E = ZH, where Z is the complex number 2 × 2 of the impedance tensor, which can be mathematically expressed as follows: The relationship between two components of the magnetic field in the horizontal direction, H x . and H y , and the vertical component (H z ) can be written as H z = T T H, where T = T x T y is a vertical magnetic sensor known as a tipper, which is very sensitive to 3D conductivity and can provide additional information, especially regarding lateral resistivity variations [47]. This analysis can be carried out by mapping the direction of the tipper-derived induction arrow to determine the position of the concentrated conductor below the surface [48,49]. The tensor element Z ij is generally represented in the form of apparent resistivity (ρ a ) and phase (ϕ), which can be expressed mathematically as follows: To determine the model, we can formulate the point where the amplitude of the skin depth signal will decrease to 1/e of the surface current as follows [50]: In addition, the Z tensor also contains information regarding the dimensions and directions. In a 1D Earth study, the impedances Z xy and Z yy (where the electric and magnetic fields are parallel) will be zero, while the off-diagonal components that combine the orthogonal electric and magnetic fields will have an equal value but opposite signs. Meanwhile, in a 2D model, the impedance components Z xx and Z yy will be equal in value but will have opposite signs, while Z xx and Z yy will have different values and opposite signs. This can be mathematically expressed as follows in Equation (4a,b): Using coordinate system rotation and distortion effect correction, Z xy can be calculated from E x and H y as well as the three field components E x , H y , H z in the electric transfer (TE) mode, while the Z yx component can be calculated from E y and H x for the transversal magnetic (TM) mode.

Magnetotellurics Data Acquisition
This study collected MT data from 26 sites using a Phoenix MTU-5 instrument in the broadband range (1-10 kHz). The sites' distribution was designed to cover the entire anomaly, which is the target of Seulawah Agam volcano, starting from that of Alue Pu manifestation at the Southern part up to Ie Sueem at the Northern part of the volcano. The MT stations were placed to cover the volcanic area and cross several local faults in view of studying the entire geothermal system in Seulawah Agam. The measurements were carried out in 3D with distances between points of 500-1000 m for x and y, respectively. The entire distribution of this MT data could be divided into five profiles that intersect the faults and volcanic manifestations, as shown in Figure 3. The distances between the points have represented the covering area of the entire volcano, where several MT studies for geothermal activity used the measurement distance more than 1 km. The measurement was conducted in a regional scale, including investigation of the tectonics, where the applied distances between sites were varied between 5 to 10 km. This is similar to the volcano-tectonic mapping of the Southern Andes, Chile, that used distances between stations of ±5 km [51] and the studies on the structure of fracture in Central Tibet with the distances between stations ±10 km [52]. Meanwhile, in the small study area, the MT method was measured with the distance between sites varying from 1 to 5 km. For example, the measurement of the 3D MT at the geothermal in Southwest Iceland, in Grand Canaria Island, and Ayrobera geothermal Northeast Ethiopia with the distances between stations of ±2 km, ±4 km, and ± 1.5 km, respectively [31,32,53]. As a result, the acquired MT site distribution in this study could provide geometrical details from the hydrothermal system of the Seulawah Agam volcano. In addition, volcanoes in the subduction areas such as Indonesia are generally different from other volcanoes' conditions from different areas, where the volcanoes are at high terrain conditions, which allows difficult on-site data collection. The coil arrangement was created according to the standard method of MT measurement [54], where the electric field x (E x ) is measured in the north-south direction and the y field (E y ) in the east-west direction, with a distance between the electrodes of ±50 m. The magnetic field components H x were placed in a north-south direction, H y in an east-west direction, and H z vertically to the ground surface. In each site station, the five components of the time-varying electromagnetic field were recorded for approximately 5 h to collect the required data. In the post-processing stage, the field data was calculated using standard robust processing software (Phoenix Geophysics SSMT2000) to ensure that the data in the frequency domain was obtained in the form of apparent resistivity (ρ a ) and phase impedance (ϕ) for each station. Finally, the *EDI files obtained from the process can be used as input parameters for further data analysis, such as in terms of the phase tensor and the inversion process.

Strike and Dimensionality Analysis
In the pre-modeling stage of MT data analysis, phase tensor analysis has been widely applied to calculate the dimensions of the sub-surface anomaly structure, such that it can provide an overview of a reasonable mathematical solution [55,56]. Here, the phase data was obtained from the real and imaginary components of the impedance tensor Z = X + iY, which is generally represented as an ellipse in terms of three values: maximum phase max ( ϕ max ), min ( ϕ min ), and slope angle (β) . Mathematically, this can be described as follows: The above equation describes the phase tensor with the parameter α as the main axis in the coordinate system. On the other hand, β is the inclination angle and has value of zero for symmetrical ellipse. Otherwise, it will describe the rotational deviation ϕ max of the symmetrical ellipse. For the 1D condition, where the ellipse is circular, the values of ϕ max and ϕ min are equal. However, for the 2D resistivity structure (ellipse), the value of β < 3 • is the 2D response of the impedance tensor [57]. Furthermore, the strike direction of the phase tensor can also be calculated using α − β, i.e., the azimuth of the major axis [55]. We performed phase tensor analysis using Profile 1, which intersects the Cempaga crater, and Profile 2, which crosses the Heutsz crater in the SW-NE direction (Figure 3), to ensure that the appropriate inversion solution could be determined for the data, as shown in Figure 4. In this part, we only provided tensor phase analysis toward profile 1 and 2, because they are the closest measurement points to the volcano. In addition, these two profiles cross several local fractures that become the main controller of the fluid elevation onto the surface. Additionally, both tracks also cross two craters which is the outcrop of the biggest manifestation from the volcanic activities of Seulawah Agam, with profile 1 crossing Cempaga crater and profile 2 crossing Heutz crater, so that the analysis carried out on the aforementioned profiles are capable of representing other areas in providing detailed information pertaining to the dimension of subsurface anomaly structure located at various depths. This dimensional analysis is required before the inversion [32], in which algorithms such as 1D and 2D are effective for structure mapping of the subsurface anomaly with the same dimension on the algorithm model [33]. Here, we combined the phase and ellipse data to ascertain the conductivity properties and the dimensional model of the anomaly, while we carried out the phase tensor analysis process using MTpy, a python toolbox for MT [58].
In terms of Profile 1 (Figure 4a), which intersects the area of the Cempaga crater with an area of 30,743 m 2 , the overall phase tensor was obtained with a value of >45 • in the period of 100-10 −4 s in a logarithmic scale of the six MT stations, which corresponded to the conductive anomaly generally obtained in a relatively shallow area [46]. However, at station 03, which is close to the Cempaga crater, the phase value was <45 • , which is in response to the resistive anomaly below the surface. Based on the form of the phase tensor, all of the data in this period range was in the form of an ellipse, indicating that the region is dominated by 1D anomalies. Only a few places were regarded as having 2D anomalies, including the structures at the stations 03 and 01. Meanwhile, in the inner structure with a period range of 10 0 -10 4 s, the tensor phase response was <45 • , indicating a resistive anomaly below the surface. In addition, the phase tensor model in this period range was not in the form of an ellipse representing the dimensions in the 2D direction, which means that a 1D model solution cannot be used to predict anomalies in this period range.
In terms of Profile 2 (Figure 4b), from the measurements at eight MT stations intersecting the Heutsz crater area with an area of 23,702 m 2 , the phase tensor exhibited a similar trend, where, in the log 10 0 -10 −4 range, a phase tensor of >45 • was obtained, indicating that conductive anomalies near the surface dominate the area. However, on the right side of the crater (site 13), a phase tensor value of 45 • was obtained, showing that the anomaly at this specific location was the same as the previous station response. The obtained phase tensor model was also in the form of an ellipse, indicating that the structure of the anomaly is 1D. Therefore, the results suggest that 1D mathematical solutions can be used to predict anomalies in the area. In contrast, in the log 10 0 -10 4 range, as a representative of the deep structure, a phase tensor of <45 • was obtained, which indicated a response to a resistive anomaly with 2D dimensions [46,55]. The 2D regional strike in Profile 1 (Figure 4c) was estimated from the Z-invariants, the phase tensor azimuth, and the tipper strike for all periods of 90 • , where the geoelectric strike analysis revealed a similar direction to the geological strike obtained from local faults found in the volcano area. The induction vectors or arrows represent the vertical magnetic field changes due to the lateral variations in electrical conductivity as a frequency function as shown in Figure 5.
Meanwhile, tipper analysis can also provide information on sub-surface structures [59], where the magnitude of the induction arrow is proportional to the intensity of the current concentration anomaly, which, in turn, is associated with the magnitude of the conductivity gradient. In contrast, the direction of the tipper arrow indicates the presence of the conductor and the distance from the resistive zone [48,60]. The tipper vector data of all measurement points were plotted in the low periods (T = 0.00515 s) to represent the nearsurface area and the high periods (T = 4.65116 s) to demonstrate the conductor in the deep anomaly. In the low periods, the tipper direction indicates a large amount of noise in the near-surface area. For example, on the western side of the volcano, the tipper data points the conductor to the west, indicating the conductor's location in the Alue PU and Ie Masam manifestations. In contrast, the northern side generally exhibited a random tipper direction focusing on the caldera margin area with manifestations such as the Heutsz crater, Ie Suum, and a number of local faults, since the conductive fluid accesses the surface via conduction and convection. On the other hand, in the high periods, all the tipper directions indicated the location of the conductor with a focus on a specific location since the tipper direction is not affected by significant noise, unlike in the area between the Heutsz and Ie Jue craters, where the conductors are likely to be in the form of heat sources. This tipper map of the tensor impedance provided information on the conductive anomaly in the near-surface area converging to the caldera margin close to the Seulawah Agam volcano.

Data Response
In general, MT data in the form of apparent resistivity and phase impedance can provide an initial picture of the characteristics of the sub-surface resistivity without prior knowledge of the depth. Here, the plot of both the apparent resistivity and the phase data from Profile 1 (Figure 6), which consists of six stations cutting across the Cempaga crater in the southern part of the volcano, are presented in logarithmic scale for both electric (TE) and magnetic (TM) transfer modes, with the data first smoothed out to eliminate any noise errors caused by the measurement instruments, magnetic storms, or human error. Overall, the TE and TM data of the apparent resistivity ( Figure 6A) exhibited the same characteristics but were in different phases. For example, in the log periods of 10 −2.5 -10 −1 s, a high resistivity value was obtained, indicating the presence of a resistive anomaly at shallow depths. The same response was also obtained in the log 10 1 -10 2.5 period, which was also caused by the resistive anomaly in the deep depth area. Meanwhile, in the medium period, i.e., log 10 −1 to 10 1 s, a low resistivity response was obtained as a response to the conductive anomaly below the surface. The difference in the response of the two modes on the MT data was due to the sensitivity to different anomalies, with the TE data being more sensitive to changes in the electric field and the TM data being more sensitive to changes in the sub-surface magnetic field. However, the overall anomaly trend indicated by the MT data was highly similar for both the TM and TE data. In addition, the resistivity curve pattern of the six stations was close to the general resistivity structure and the alteration model, which indicates geothermal potential [53]. The phase value from the MT data could also be used to determine the state of the conductor, as shown in Figure 6B. Both phase values on track 2 indicated different responses since they are in different quadrant phases, i.e., TE mode in the range of 0 • -100 • and TM mode in the range of 100 • -200 • . Generally, the phase data is plotted in the same quadrant by reducing the phase TM data by 180 • [44,61]. Here, while the quadrants are different, taken as a whole, the data phase on track 2 indicated a trend response from the same anomaly. When the apparent resistivity value was high, the phase value obtained was <45 • . Conversely, a phase value of >45 • corresponds to the presence of conducting objects below the surface [55].

One-Dimensional Inversion
A 1D earth model was created using the Occam inversion code developed by Constable et al. [62] and Bostick [63] to produce an appropriate solution for the sub-surface 1D structure. The MT sounding curve and the 1D inversion data can provide information regarding the general structure expected from the subsequent 2D and 3D modeling. Here, we only conducted 1D inversion at three stations close to the volcanic manifestations, as shown in Figure 7. The apparent resistivity and phase data used as input parameters in the inversion process were corrected in advance, especially in terms of the high frequencies caused by a large amount of noise near the surface [64]. The 1D inversion process was carried out with the same conditions for the three MT stations. Five specific layers were used to classify the resistivity data based on a standard initial model [65], with a homogeneous half -space of 100 ohm·m. Meanwhile, the Occam model was also given five layers with a solution depth of 100-5000 m, with 10 iterations used to produce a model suitable for sub-surface conditions, as shown in Figure 7. The data at station 02, which is close to the Cempaga crater (Figure 7a), returned 1D Occam inversion results that were similar to the resistivity layer at a depth of 5 km with an RMS of 0.02%, where a depth of 0-300 m is indicated by a moderate resistivity (60-250 ohm·m) in response to the aquifer layer near the surface. Furthermore, a depth of 300-1000 m is a conductive area (<10 ohm·m), presumed to be a layer of clay cap that prevents fluid from rising to the surface. At a 1-5 km depth, the resistivity value was moderate and equal to the aquifer layer. It was deemed to be a response to the conductive anomaly of the sub-surface reservoir and caprock. In terms of station 03, which was located near the Cempaga crater, 1D modeling results were obtained with the same model and layer response as those obtained with station 02, especially at deep depths. At the same time, they were comparatively different at shallow depths with an RMS of 0.02%. One of the factors that led to the difference in response is the distance from the measurement point to the crater location, meaning station 03 is more conductive than station 02 in shallow areas since it is closer to the crater where the fluid activity rises to the surface.
In addition, at site 19, which is located close to the Heutz crater on the northern side of the volcano, a resistivity model with an RMS of 0.5% was obtained, with the same characteristics as the Cempaga crater (site 02), especially in terms of the deep anomaly. Here, a depth of 300 m to 1 km resulted in a low resistivity value, indicating that the area is dominated by conductive rock and is deemed a clay cap layer in the form of liquid fluid that requires fractures to rise the surface. This conductive area indicates a potential geothermal zone that could be further modeled using 2D inversion. In addition, at a depth of 1-2.3 km, a high resistivity value was obtained, likely due to a hydrothermal reservoir, while at this point, a different response was obtained at a depth of 2.3-5 km, with a relatively low resistivity compared with the reservoir layer. This layer is presumed to be a heat source area where a geothermal zone is located, while it could not be identified at the other sites since 1D solutions are not sensitive to 2D anomalies. This can also be verified by phase tensor data analysis, where the layer is dominated by 2D anomalies, meaning 1D mathematical solutions are not applicable.

Two-Dimensional Occam Inversion
For the 2D MT inversion process for the Seulawah volcano, we used the Occam 2D inversion code developed by [66], which is based on the minimization of the following unconstrained function: where ∂ y m 2 + ∂ z m 2 is an equation model for the roughness model for data smoothing, and µ −1 is the Lagrange value for ascertaining the minimum and maximum values. The field data will fit the calculated data, W is the diagonal weighting matrix M × M, d is the observation vector, and F(m) is the model response. The 2D modeling was carried out for all the MT profiles covering the local faults and volcanic manifestations. However, to explain the geothermal formation system in 2D, only two paths were necessary, i.e., Profile 1, which covers the Seulawah Agam volcano area and various manifestations, such as the Cempaga crater and Alue PU, and Profile 2, which covers the Heutsz crater on the northern side and various other local faults, as shown in Figure 8. Based on the 2D Occam model, the resistivity value was obtained from log 0.5-3.0 ohm·m to a penetration depth of 5 km below the surface. In terms of Profile 1, there were three anomaly conductors, namely, C1, C2, and C3, whereas the anomaly resistors were identified in two locations, namely, R1 and R2. At shallow depths (<1 km), a number of conductive anomalies with a value of 10 ohm·m were obtained in response to the manifestation area, such as C1 on the western side, where there is a manifestation of Alue PU in the form of hot springs, and location C2, which includes the Cempaga crater from a distance of 2-8 km and is connected to the conductor on C3. Furthermore, the conductive anomaly from C1 to C3 was presumed to be a clay cap area that becomes a hydrothermal flow that delivers fluid to the surface. A report suggests that several volcanoes have an average temperature of the hydrothermal reservoir at 220 • C [67]. However, in [28], using geothermometer gas with the Giggenbach method, the reservoir temperature of Seulawah Agam volcano was predicted to be >225 • C. According to Houcstein [9], this suggests the volcano is suitable for power plant development. In addition, the 2D model also indicated various conductor zones connected to local faults in the measurement area, such as C1 in the S-NW direction and conductors C2 and C3 in the NW-SE direction. These faults are associated with low-resistivity anomalies that can be interpreted as channels transporting thermal energy from the heat source to shallow depths [68].
In addition, there is an accumulation of volcanic deposits on the surface area due to topographic differences with high resistivity (350 ohm·m). This anomaly can be found on the western and eastern sides, where it prevents fluid from rising to the surface, suggesting no formation of hot springs. Resistors R1 and R2, with a resistivity of >600 ohm·m, are separated by a conductor at a measurement distance of 4-6 km. These two resistors extend from the surface almost to a depth of 2 km, which is thought to be in response to the heat source area. However, one of the unique features of Seulawah Agam compared with other volcanoes is the presence of heat source areas in two locations, that is, the northern side centered around the Heutsz crater area and the southern side around the Cempaga area. Based on this information, we attempted to invert Profile 2, which only covers the Heutsz crater, as shown in Figure 9. The results indicated that four-conductor anomalies and two-resistor anomalies could be found at a depth of 1-5 km below the surface. Conductor anomaly C1 was connected to C2-C4, which is located at 0-3 km below the surface. This conductive area with a resistivity of <10 ohm·m is thought to be a clay cap layer with highly conductive characteristics that acts as a cover layer that stores heat and hot liquid in the volcanic system [69,70]. The formation process of clay cap layers in geothermal systems relates to liquids at lower temperatures or to steam condensing at shallow levels to create a highly conductive condensate layer [2]. In addition, the three conductors are separated by various other resistive anomalies in the area near the surface, ranging from 31-110 ohm·m. For example, at 3 km (R2) from the measurement profile, a volcanic deposit was deemed to exist, as was the case at 6-9 km, which is believed to be a response from volcanic deposits due to differences in the topography of the volcano. In addition, the resistive R2 anomaly is closely related to the contrast response to two local faults in the almost NW-SE direction. At this location, the fault is divided into two branches: one extending from the Ie Asam manifestation to Alue Pu, and the other the Seulimeum segment regional fault that is part of the GSF on the island of Weh Sabang [39,71].
In the manifestation area of the Heutsz crater, a low resistivity value (<10 ohm·m) was obtained in response to the conductive zone, where, on the western and eastern sides of the crater, fractures are crossing through the clay cap layer, which helps the fluid to rise to the surface, forming a hot spring with a temperature of 98 • C [28]. Liquids circulating at deeper levels tend to flow to the surface through the path with the highest permeability, such as through fracture, which provides conductive heat transfer between the heat source and the reservoir. The 2D Occam results also indicated that the reservoir layer is >2 km below the surface with a resistivity of 40-130 ohm·m, which is desirable for drilling. In addition, at 5-12 km from the measurement path, a resistive anomaly (>600 ohm·m) was identified, which was deemed to be the heat source of the Seulawah Agam geothermal zone. However, the heat source area is focused on one location only, as the second path only passes through the northern part of the volcano. This indicates that the Seulawah Agam heat source area is not centered under the top of the volcano but is located on its northern and southern sides.

Hydrothermal System from Three-Dimensional Resistivity
The hydrothermal system of the volcano is controlled by the activity of faults and fractures where fluids pass through [72,73]. The visual representation of the 3D model was carried out by combining all of the 2D Occam inversions sliced at some depth, as shown in Figure 10. The resistivity distribution of the preferred 3D structure is now correlated to the geology of the survey area. In the near-surface area of 0-500 m in depth (Figure 10a), conductive anomalies are dominantly distributed irregularly in several locations, such as on the northwest and southeast sides, where the manifestations are generally found. Conductive anomalies can also be found around the volcano area since, in the volcanic area, there are several conductive craters, such as the Cempaga and Heutsz craters. Meanwhile, resistive anomalies are present in various locations, such as at the western and central sides, which are far removed from the manifestation locations. In addition, the contrast between conductor and resistor anomalies is closely related to the presence of various local faults in the NW-SE direction, which include the fault close to the Ie Suum manifestation, the fault from the direction of the Cempaga crater to Ie Jue, and the fault near Alue PU extending to Ie Suum, which is part of the Seulimeum Sumatra Fault segment. In terms of this area near the surface, the Occam inversion results indicated that a number of the conductors were caused by conductive pyroclastic manifestations and materials, meaning the conductive anomaly is more than the resistive one. At a depth of 500-1000 m (Figure 10b), where the 2D cross-section revealed the depth to be a clay cap area covering the reservoir zone, the conductive anomalies are more concentrated in a number of manifestation areas, which are clearer than the previous resistivity depth. For example, in the Ie Jue and Ie Busuk manifestation areas, there is a conductive anomaly in the SW-NE direction that is presumed to be hydrothermal flow caused by local faults. A similar phenomenon can be observed in the Alue Pu, Cempaga, and Ie Masam areas, which exhibit a conductive anomaly with the same patterns as the fault direction, i.e., an almost NW-SE direction. Meanwhile, on the eastern side of the volcano, an anomaly was identified that is more resistive than the clay cap layer and the hot spring over the geothermal system. This is believed to present a response to the reservoir zone that carries fluid to the clay cap layer. At a depth of 1500-2000 m (Figure 10c), a conductive anomaly was found in various manifestations. On the northern side of the volcano, the anomaly began to separate between the manifestations and the other hot springs, while on the western side, the conductive anomaly was still connected between specific manifestations. This further demonstrates the uniqueness of Seulawah Agam, where the heat source zones lie on the northern and southern sides, as observed in the 2D Occam cross section. The same phenomenon was also evident at a depth of 2000-2500 m (Figure 10d), where the conductive anomaly was only found in the manifestation of the heat source on the southern side, namely, Alue Pu and Ie Masam, but the anomaly begins to separate, which shows the discontinuity of the fault in that area, indicating the fault's discontinuity. In addition, at a depth of 1500-2500 m, various resistive anomalies can be found, which are widely distributed around the volcanic area. A large number of resistive anomalies beneath the highly conductive clay cap layer is generally interpreted to indicate a zone of high-temperature water circulation. Therefore, these resistive anomalies can be interpreted as forming the center of the geothermal reservoir system, which is a potential target for future development.
At large depths of 3000-3500 m (Figure 10e), there was deemed to be a heat source layer centered on the northern and southern sides of the volcano, where the presence of a number of anomalies is evidenced through the high resistivity. Meanwhile, on the central side of the volcano, there exists an anomaly with lower resistivity, which is believed to be within the hydrothermal reservoir area that limits the two heat source areas. The conductive anomaly at this depth is a response to the various local and regional faults that are the main mechanism of the Seulawah Agam volcano. A similar phenomenon can be observed at a depth of 4000-5000 m, where the identified conductors are strongly associated with faults, especially the eastern side, which is a regional fault (part of the Seulimeum segment) of the Sumatran fault [21,42]. Overall, the 3D resistivity structure data provided an understanding of the geothermal system of Seulawah Agam, where the low resistivity corresponds to the thermal conduction zone, the fault contrast, and the clay cap layers; the middle resistivity corresponds to the reservoir zone over the systems; and the high resistivity presents a response to the heat source on the northern and southern sides of the volcano.

Discussion
Magnetotelluric imaging has become a common method in investigating volcanoes worldwide [3,20,31]. The method works by utilizing the Earth's magnetic field with a frequency of 10 −4 -10 4 Hz. In this study, we applied the MT method by 16 stations covering the volcano area to investigate the hydrothermal system of Seulawah Agam, which is an important aspect of geothermal exploration. Generally, volcanoes in Indonesia, especially Seulawah Agam, are located at high terrain conditions with an altitude of 1810 mdpl, thus there are only a few regional studies that cover the entire volcanic area. There are currently reports published, but these are limited only to certain manifestations such as Ie Suum and Ie Jue conducted by geophysics and geochemistry [28,34]. Herein, we compared the results obtained based on the 3D structure of resistivities with those obtained in previous reports. Based on the 2D structure (Figures 8 and 9), along with 3D anomaly (Figure 10) of MT resistivities, the main anomalous features found in the resistivity cross-sections which indicate a geothermal system, including clay cap acquired at a depth of 1 km, a reservoir layer at 2 km, and a heat source zone at 3 to 5 km below the surface. The models we had constructed in this study are preferred for geological interpretation [42].
The GSF tectonic activity on the island of Sumatra manifests itself locally through the formation of numerous buried faults striking in the Seulawah Agam volcano in the SE-NW direction. These faults are associated with a low-resistivity anomaly that can be interpreted as a channel transporting thermal energy from a heat source to a shallow depth [35]. Therefore, a number of low-contrast anomalies with various conductors are closely related to the presence of the various local faults and fractures of Seulawah Agam. The flow zone of the volcano's hydrothermal system is presented in Figure 11. The spatial distribution of the anomalies with high conductivity in the northern and southern areas corresponds to volcanic manifestations, where partial dissolution of the conductive material responsible for the conductive anomalies in the hot spring [74]. At deep depths, high-resistive anomalies can be found on the northern and southern sides, representing the heat source area of the Seulawah geothermal system. Based on this information, we simulated the hydrothermal flow of the volcano (as shown by the black lines in Figure 11), where there is a distribution of manifestations and direct contact with faults that are the medium for fluid transfer to the surface. In addition, the hydrothermal model also revealed that the fault on the right does not feature any manifestations since the depth of the fault does not penetrate the clay cap layer as a fluid cover of the reservoir, thus preventing the fluid from rising to the surface. On the contrary, the other two faults penetrate the clay cap area, showing that hydrothermal displacement from the reservoir to the surface takes place and forms various manifestation areas. Geothermometer data predicted that the reservoir temperature of Seulawah Agam reaching >225 • C [28], which is in line with the thermal gradient of 50 • C/km [75]. Although the thermal gradient was close to that of regional estimation, investigating a geothermal requires magmatic depth as the heat source mobilizes the hydrothermal system. The 2D MT model of profile 1 (Figure 8) and profile 2 (Figure 9), which cross the volcano, shows resistivity anomaly with a value of >600 ohm·m. It was suspected to be a heat source zone at the 3-5 km depth, in which the anomaly is located at the west side, directing to the Alue PU manifestation, and at the east side, to the Heutz crater.
Based on the pseudo-3D resistivity, the conductive anomaly was estimated as clay cap with 1 km depth in the volcano is generally directed to tracks, the west, and the east. According to on-site observation on the west side, several manifestations were obtained were obtained from the geochemistry studies, starting from the Heutz crater at 76.6 • C, the Ie Jue crater at 97.61 • C, and the Ie Busuk crater at 47.49 • C [28]. Meanwhile, on the east side, the clay cap area flows from the Cempaga crater to Ie Suum with a temperature of 86.02 • C. The fluid obtained in the clay cap layer could flow up to the surface, forming manifestations resulting from fractures that act as access media. The fracture could be observed in the 3D resistivity structure of the MT method. This interpretation is corroborated by the Transient Electromagnetic data specifically used to map the shallow structure (< 1 km) of the volcano. Altogether, the data revealed the presence of a fault directing to the Heutz crater and Ie Jue manifestations along the clay cap conductive anomaly [21,76]. At the Ie Jue and Ie Busuk manifestations, there were also found some fractures with low resistivities, indicating the fluid movement from reservoir layer, which is in line with the magnetic data having low intensity at the fracture area, as well as high conductivity that was obtained from VLF-R method [34]. The present pseudo-3D model is a combination of 2D model inversion of several profiles, although the analytical data of the tensor phase showing that on the layer characterized by 2D and 3D anomalies. Indeed, many studies suggest that the 2D inversion response could have explained the geothermal system very well [31,32,35], especially for volcanoes located in tropical regions on high topographies. However, for better results, a full 3D inversion should be applied [31,32,35].

Conclusions
Seulawah Agam has been intended for a potential supply source for a geothermal power plant. In geothermal exploration, the mapping of fractures along with fluid flows is a fundamental element because it provides information about the hydrothermal system formation and the flow path connected to the reservoir wells for fluid production. In this paper, a 3D MT-based survey was conducted to study the resistivity structure associated with the hydrothermal activity and the geothermal system of the volcano. The near-surface conductive layer resulting from the Occam inversion indicated the presence of volcanic deposits that can be largely found at the top of the volcano area. Another conductive anomaly was found to be related to the presence of various manifestations on the northern side, namely, the Heutsz, Ie Jue, and Ie Busuk craters, while on the southern side, there exists a number of manifestations related to the Cempaga crater, as well as to Ie suum and Alue PU. In addition, at a depth of 1 km, a low-resistive anomaly was identified as a response to a highly conductive clay cap area that functions as a zone of high-temperature water circulation and covers the fluid rising above the surface. Meanwhile, the resistivity model revealed various fractures and faults that act as a fluid access pathway whereby surface manifestations are formed. The reservoir layer was found to be at a depth of >1.5 km, with a relatively low resistivity compared with the clay cap area. This layer is the center of the geothermal reservoir system and presents a potential drilling target for future geothermal energy development. The 3D resistivity structure also revealed that a heat source area at a depth of 3 km is present in two locations with high resistivity values, i.e., on the northern side centered on the Heutsz crater area and on the southern side centered on the Cempaga crater area, which results in the hydrothermal flow between the manifestations in the north being different from that in the south. Based on data analysis, the MT method was applied to study the geothermal system of the Seulawah Agam volcano. However, the MT data must be 3D-inverted using ModEM [76] or WSINV3DMT [45] software to produce a better model.