The MT method is a passive source electromagnetic method sensitive to the electrical conductivity contrast of the subsurface. According to the purpose of the study, it provides the opportunity to research from a few hundred meters to tens of kilometers. It is based on the measurement of the time variation of the electric (Ex, Ey) and magnetic (Hx, Hy, Hz) fields perpendicular to each other in the field. Then, the time series data is transferred to the frequency domain by utilizing the Fourier transform (Chave and Jones, 2012). The complex impedance tensor is derived from the relationship between the simultaneously measured electric and magnetic fields.

Where, Z is a frequency dependent 2x2 tensor which called impedance tensor. The impedance tensor also contains information about the electrical resistivity distribution of the underground, the dimensionality and direction of geo-electrical structures. By utilizing the above equation, apperent resistivity (ρa) (Cagniard, 1953) and impedance phase (\(\varnothing\)) can be derived;

## Data and Modeling

The long period MT data which used for this study was collected in August 2018. The data was recorded along 15 hours at 39 sites (Fig. 1b). Data for five components (Ex, Ey, Bx, By, and Bz) were recorded at each measurement site using Metronix induction coil magnetometers, which are sensitive to the AMT band, along with non-polarized Pb-PbCl electrodes. We implemented the remote reference procedures by using Egbert’s code (Egbert and Booker, 1986; Egbert, 1997) to estimate the MT impedance (Z) tensor. After time series analysis, we produced frequency-domain data from 0.001 Hz to 10 kHz for 50 discrete frequencies. In the study area, the vertical coil was often not sufficiently covered, as hard ground was usually encountered at the measurement points. Therefore, the vertical magnetic field component (i.e. tipper data) could not be recorded without noise due to the effect of wind and environmental noise. Because of noisy recording, tipper data is not included in the inversion.

To determine the dimensionality of the study area, the MT data was analyzed using phase tensor analysis (Caldwell, 2004), which is unaffected by galvanic distortion. Depending on frequency the phase tensor is conventionally shown as an ellipse. This ellipse is basically represented by principle axes and skew angle (*β*). The skew angle contains information about the dimensionality of structures. Extreme *β* values (|*β|* > \({3}^{o}\)) indicate the presence of three-dimensional structures. Figure 3 shows the phase tensor maps of the 0.1, 1, 10, and 100 second periods plotted from the observed data. Ellipses are colored according to the skew angle. We generate the phase tensor ellipses by utilizing MTpy software package (Kirkby et al., 2019; Krieger and Peacock, 2014). The fact that the \(\beta\) values are quite high, especially in the 10s and 100s periods, indicates that the structures at those depths are 3-D. Another remarkable issue here is that the maximum axes of the phase tensor ellipses in the 10 s period are oriented to the north. In fact, the minimum axes of some ellipses have very small values. We think that this is caused by conductive sea water (~ 0.3 \(\Omega\)). We did not plot the induction vectors because the tipper data is noisy.

We utilized the 3D Grid-Tools software academic version to create a discrete earth model and to include topography and bathymetry information into the created mesh. The shortest period of our data is 0.0001 s. Thus, the first thickness was calculated as 7 m. As for an initial model of inversion process, we started a homogenous medium with 100 Ωm. The penetration depth was estimated as 159 km, for 1000 s the longest period and with a 100 Ωm. Further, we also added 12 blocks in the x, y, and z directions to avoid negatory effect of the boundary conditions (Ranaganyaki and Madden, 1980; Lindsey and Newman, 2015) with a size increasing factor of 1.15 for each direction. The total size of model comprise of 118x118x103 cells. While the topography variation is represented with high resistive blocks (e.g. 1016 Ωm), sea water is represented with conductive blocks of 0.3 Ωm resistivity value (Fig. 4). Since the whole of the model is very large, we only showed the study area and its vicinity in the figure. As a matter of fact that, nondiagonal components of impedance tensor (Zxy, Zyx) were used with 10% error floor. Due to their low amplitudes (Lindsey and Newman, 2015), we added 30% error floor to diagonal components (Zxx, Zyy).

ModEM3DMT (Egbert and Kelbert, 2012; Kelbert et al., 2014) parallelized nonlinear conjugate gradient (NLCG) algorithm was used to carry out inversion process. The purpose of the algorithm is to minimize the penalty function. For detail information about ModEM3DMT, we refer to the reader Egbert and Kelbert (2012) and Kelbert et al (2014).

There are several synthetic tests carried out to determine the most proper inversion parameters (Miensopust et al., 2013; Avdeeva et al., 2012; Tietze and Ritter, 2013; Miensopust, 2017; Slezak et al., 2019). Recently, detailed study about inversion parameters was performed by Robertson et al. (2020). They used some part of AusLAMP Project and they tried to specify most reasonable initial damping parameter (\(\lambda\)), initial model resistivity, covariance, cell size etc. We also implemented some synthetic test using several different models and inversion parameters. We used a similar synthetic model (Fig. 5) with Tietze et al., (2015). While creating the synthetic model, we added the blocks representing the sea and the topography to the model as in Fig. 4.

It can be seen that the synthetic model consists of two adjacent blocks which have 50 Ωm and 500 Ωm resistivity values in Fig. 5a. A conductive block with 1 Ωm is embedded into these adjacent blocks at depths from 0.97 km to 3.04 km. The initial model includes both bathymetry and topography layers. Topography and bathymetry information implemented into the model by coordinate information of the MT stations. The X direction represents the geographic north. Our goal in performing the synthetic test is to see how well we can model the anomaly with our measurement distribution. In order to see the measurement distribution clearly, only the study area is focused on in Fig. 5. As with field data, topography and bathymetry were added to the initial model while modeling synthetic test data.

Sea level is located in the x direction of the model and 1.5 km far away from the first station profile approximately. Inversion results are shown in Fig. 6.

Figure 6a shows 1.94 km depth slice of inversion result. The North-west part of the conductive block was not resolved correctly because of sea effect. As we referred before sea level is located north side of the station geometry.

After checking reliability of model parameters, we proceeded to processing of field data. During data processing, we used the model parameters as we mentioned above. Since the cell size of the model to be inverted very big, it requires very powerful computers. For this purpose, we benefited from TRUBA (Turkish Science e-Infrastracture). We used 80 cores throughout the process with 72 iteration steps. Iterations started with an RMS value of 22.08, and reached an RMS value of 2.85 at the end of the 72nd iteration. We specified the range of the resistivity scale as 1-1000 Ωm. Figure 7 shows four cross-sections in the north-south direction. It can be distinguished three distinct layers from the surface to the deep down in Fig. 7–8. The cross-sections from east to west are illustrated in Fig. 8. In the vertical direction at depth slices of the obtained 3D model based on MT inversion result are displayed in Fig. 9.

As can be seen from the sections, basically 3 different units stand out up to a depth of 5 km from the surface. At the top are conductive the Eocene volcanics, just below this there are altered aquifers with relatively higher resistivity. At the bottom is the bedrock intrusion with high resistivity. Since the measurement points are very close to the sea in the north, we cannot get information from deeper zones due to the sea effect. As an example, we show the observed and calculated data fits of the apparent resistivity and phase curves at two MT sites (05S and 09S) in Fig. 11.

We demonstrated resistive (> 150 Ωm) bedrock as volumetrically (Fig. 10). As can be clearly seen from the volumetric representation, the bedrock rises and gets closer to the surface at the middle of the study area. It moves to the northward slope and extends under the young sediments by making undulations.

A geothermal borehole was drilled by the Termal Municipality in the study area, reaching a depth of 804 meters (Fig. 1) (Information are provided by Termal Municipality). During the drilling, hot water at a temperature of 60°C was encountered in the altered serpentine unit at a depth of 680 meters. Based on the 3D subsurface conductivity model, two cross-sections (N-S and E-W) were generated over the borehole (Fig. 12). The red pipe in the figure indicates the location of the drilling. Figure 12 clearly shows a high resistivity intrusion at approximately 700 meters depth, which corresponds to the outlet of the hot water, just below the borehole. We propose that the hot water is emanating from these depths due to the contrast in resistivity between the high resistivity intrusion and the surrounding, relatively more conductive structure.