remote sensing Three-Dimensional Magnetotelluric Characterization of the Travale Geothermal Field (Italy)

: The geoelectrical features of the Travale geothermal ﬁeld (Italy), one of the most productive geothermal ﬁelds in the world, have been investigated by means of three-dimensional (3D) magnetotelluric (MT) data inversion. This study presents the ﬁrst resistivity model of the Travale geothermal ﬁeld derived from derivative-based 3D MT inversion. We analyzed MT data that have been acquired in Travale over the past decades in order to determine its geoelectrical dimensionality, directionality, and phase tensor properties. We selected data from 51 MT sites for 3D inversion. We carried out a number of 3D MT inversion tests by changing the type of data to be inverted, the inclusion of static-shift correction at some sites where new time-domain electromagnetic soundings (TDEM) were acquired, the grid rotation, as well as the starting model in order to assess the connection between the inversion model and the geology. The ﬁnal 3D model herein presents deep elongated resistive bodies between the depths of 1.5 and 8 km. They are transverse to the Apennine structures and suggest a correlation with the strike-slip tectonics. Comparison with a seismic velocity model and well log data suggests a highly-fractured volume of rocks with vapor-dominated circulation. The outcome of this study provides new insights into the complex geothermal system of Travale.

The Larderello-Travale geothermal area (LTGA) is located in Italy (southern Tuscany). It is here where electrical power production and delivery from geothermal resources first began, in 1913, with the first 250-kW geo-thermoelectric unit in the world [24]. The LTGA has been extensively studied over the years and has been classified as a convective (vapordominated) and intrusive (magmatic heat source) geothermal play type [25,26]. In fact, in southern Tuscany, the magmatic activity and geodynamic setting gave rise to a vast Even though several differences occur, the Larderello and Travale fields belong to a common deep geothermal reservoir [42]. The Travale geothermal area is located to the south-east of the town of Larderello and covers a surface of around 50 km 2 [43]. Here, in 2018, 38 wells and 8 plants were in production with an installed capacity of 200 MWe. Two main reservoirs comprise the Travale geothermal field [42,46,47]. The shallow reservoir is a conventional hydrothermal system hosted in the evaporite-carbonate units with an average temperature of around 200 °C. The deep reservoir is hosted by the metamorphic basement and Neogene granitoids at a depth of about 2500-4000 m. The recent deep exploration program started in 1992, registering superheated steam at 350 °C with an initial pressure of 7 MPa at a depth of around 4 km [42,48]. Temperature logs from wells measured 300-350 °C at a depth of 3 km [43,47].
The Travale area is located in the inner part of the Northern Apennines, a sector o the Apennine orogenic belt that developed as a consequence of the Cenozoic collision between the European (Corso-Sardinian block) and the Adria plates [49]. The tectonic evolutional model of the Northern Apennines, proposed by several authors [50][51][52], and references therein, implies two main deformational processes: (i) an initial one related to eastward-migrating compressional tectonics and (ii) a subsequent extensional tectonic process that migrated eastward, and which has been affecting the inner part of the orogenic belt since at least the early Miocene. Alternative models have also been proposed to describe the tectonic evolution of the inner Northern Apennines [49,53]. These models Even though several differences occur, the Larderello and Travale fields belong to a common deep geothermal reservoir [42]. The Travale geothermal area is located to the south-east of the town of Larderello and covers a surface of around 50 km 2 [43]. Here, in 2018, 38 wells and 8 plants were in production with an installed capacity of 200 MW e . Two main reservoirs comprise the Travale geothermal field [42,46,47]. The shallow reservoir is a conventional hydrothermal system hosted in the evaporite-carbonate units with an average temperature of around 200 • C. The deep reservoir is hosted by the metamorphic basement and Neogene granitoids at a depth of about 2500-4000 m. The recent deep exploration program started in 1992, registering superheated steam at 350 • C with an initial pressure of 7 MPa at a depth of around 4 km [42,48]. Temperature logs from wells measured 300-350 • C at a depth of 3 km [43,47].
The Travale area is located in the inner part of the Northern Apennines, a sector of the Apennine orogenic belt that developed as a consequence of the Cenozoic collision between the European (Corso-Sardinian block) and the Adria plates [49]. The tectonic evolutional model of the Northern Apennines, proposed by several authors [50][51][52], and references therein, implies two main deformational processes: (i) an initial one related to eastwardmigrating compressional tectonics and (ii) a subsequent extensional tectonic process that migrated eastward, and which has been affecting the inner part of the orogenic belt since at least the early Miocene. Alternative models have also been proposed to describe the tectonic evolution of the inner Northern Apennines [49,53]. These models have revealed a complex tectonic evolution during the Miocene-Pleistocene with alternating compressive and extensional tectonic events, which contrast with the previously cited model supporting uninterrupted regional extensional tectonics since the early Miocene. The geological framework of the Travale area is well documented [33,47,54,55]. Figure 2a reports a simplified geological map of the Travale geothermal area based on data from the Tuscany regional website (Geoportale Geoscopio website [56]). Figure 2b shows a NW-SE-oriented schematic structural model of the Travale area [46]. The outcropping carbonate formations represent the recharge area of the shallow reservoir, which is cooled down and is therefore less exploited than the deep reservoir. The structural-stratigraphic conceptual model is depicted in Figure 2c and outlines the following geological units: (1) Neogene sediments, (2) Ligurian and sub-Ligurian complex (Jurassic-Eocene), (3) Tuscan Nappe (Triassic-Lower Miocene), (4) tectonic wedge complex (Paleozoic-Triassic), (5) phyllitic and quartzitic complex, (6) mica-schist complex, (7) gneiss complex, (8) Pliocene granite and Quaternary granite. The shallow geothermal reservoir is hosted in units three and four, while the deep metamorphic reservoir is found in units five, six, seven, and eight. The intrusive complex represents the heat source of this long-lived geothermal system [36].
From seismic reflection data, two seismic markers have been detected in the deep structure of the Travale field [33,55,57]. The shallow marker is referred to as the "H-horizon" and appears as a discontinuous high-amplitude reflector. It is located in the metamorphic reservoir at a depth of 2.5-4 km, above the Pliocene granitoids, and represents the actual reservoir target (see Figure 2b). In fact, most of the drilled wells encountering the H-horizon are productive (e.g., in the sectors of Sesta, Radicondoli, and Montieri depicted in Figure 2b) [33]. The deep marker "K-horizon" is represented by a high-amplitude discontinuous reflector showing local bright spot features at a range of depths between 3 km (in Larderello) and 8 km (in Travale). Its origin and nature are still under debate [36,38,[57][58][59]. The K-horizon has often been associated to the 400-450 • C isotherm and the occurrence of supercritical fluids, but recent deep drilling in the proximity of Larderello measured a temperature above the K-horizon of about 510 • C [45]. The most recent seismic study in the Travale area was a local earthquake tomography analysis derived by travel-time inversion [39]. The 3D model of P-wave velocity (v p ) highlighted a deep-rooted low-velocity body (v p = 5.7 km/s) below the Travale area, between the H-and K-horizons (see Figure 2b).
The Travale geothermal area is characterized by a low gravity anomaly of about 10-20 mGal [34]. This anomaly is spatially coherent with the highest heat flow of about 200 mW/m 2 measured at the surface. Thermal numerical-modeling studies were performed in order to investigate the nature of the heat source of the reservoir as well as to predict the future evolution of the geothermal system [48].
Previous MT studies of the Travale geothermal area resulted in 2D resistivity models [29,40]. The 2D inversion followed the algorithm of Rodi and Mackie [60], using a priori geological sections of the profiles as their starting model. The inversion results showed, in essence, a highly conductive shallow subsurface and a large resistive body (around 800 Ω·m) overlying the K-horizon from a depth of 3 to 6 km. However, since previous studies are limited to 2D interpretation, this MT dataset needs to be re-examined in the light of the most recent 3D MT inversion techniques. Tuscan Nappe sediments (Late Triassic-Early Miocene), (5) Tuscan Nappe: basal evaporites (Late Triassic), (6) Metamorphic units (Triassic-Paleozoic). The red squares are the 51 sites selected for MT inversion. Radicondoli7bis is the well shown in Figure 13. The black curves are the main faults [56]; (b) Schematic conceptual model of the Travale area with NW-SE orientation where seismic horizons H and K are indicated (modified from [42]); (c) Schematic sketch of the tectono-stratigraphic and hydrogeological complexes of the LTGA (modified from [24,47]).

Acquisition and Processing
The Travale geothermal field was explored by means of three MT campaigns conducted in 1992, 2004, and 2006-07 [29,40]. The sites are depicted in Figure 1 with different colors according to the year of acquisition. The investigated area covers a surface of 48 km 2 . The average distance from the coastline is 60 km. The minimum distance between the sites is 191 m (sites d8-d9) and the average distance is 3.8 km. The minimum and the maximum elevations are 314 m (site c9) and 652 m (site k5) a.s.l., respectively.
The acquisition settings of the three datasets were slightly different. The 1992 dataset is composed of two long E-W profiles and was acquired as part of an exploration project using Phoenix V-5 systems [28,61]. A remote reference processing technique was applied using the remote site located on the Island of Capraia (Tuscan Archipelago).
Regarding the 2004 dataset, four-component (and five-component for some sites) data were collected using the Phoenix MTU system in the frame of the INTAS (EU) Project. At each site, the MT signal was recorded overnight for (at least) 12 h in the range of 0.003-993 s. Quality controls were carried out on time series, spectra, and transfer functions and a strong noise signature was observed for some sites, as reported by [29]. The source of the noise for short-period MT data was related to the local power plants and geothermal exploiting activities. To solve this, the remote-reference processing technique was applied using simultaneously measured data at the local sites. Long-period MT data were contaminated by noise arising from nearby electrified railways. To solve this, a remote site was placed in Sardinia, around 500 km southwest of Travale. The whole dataset was processed using the Phoenix Geophysics software [62], based on the remote reference robust processing method [63].
The 2006-07 dataset was acquired in the context of the I-GET (EU) Project. Longperiod MT data (period range: 0.2-1000 s) were measured using the MT systems "NIMS" or "LEMI," while audio-MT data (period range: 0.001-10 s) used a Stratagem system [40].
In this study, the whole dataset depicted in Figure 1 (86 sites) was adopted for the dimensionality and directionality analyses. All the sites shown in Figure 1 are meant to provide a complete overview of the MT data previously acquired in Travale and to eventually foster new MT campaigns by ensuring high-quality data and regular site distribution on a 3D grid. The availability of a large MT dataset from different campaigns allowed for comprehensive dimensionality analysis. For 3D inversion modeling, we selected the subset of 51 MT soundings acquired in 2004. The main reasons for this choice were based on the low level of noise, the similarity of the period range covered by the responses and the regular spatial distribution of the site locations. Some MT sites were discarded from 3D inversion due to either their poor quality (the 2006-07 dataset) or to their inadequate spatial coverage (the 1992 dataset is composed of two isolated profiles around Larderello). These two factors would have corrupted the 3D modeling and the meshing, respectively.
The data quality of the sites selected for inversion is shown in Figure 3. The observed off-diagonal apparent resistivity is not corrected for static shift, which can explain some scattering of the curves.

Dimensionality Analysis
The first evaluation of the raw dataset (before considering any static shift correcti or applying any rotation) involved a thorough analysis of the geoelectrical dimensiona by means of the rotational invariants of the impedance tensor [64,65], which gives initial indication of the spatial distribution (1D, 2D, or 3D) of the electrical resistivit depth and identifies the presence of galvanic distortion (3D/2D). Details are given in Supplementary Material (part A). The analysis results show a discrete number of 1D 2D tensors for short periods and a predominant number of 3D cases, which cle outlined the 3D nature of the subsurface electrical structures and thus the priority gi to 3D inversion.
Further insight into the data came from the phase tensor analysis [66][67][68]. The ph tensor (Φ) is based on the observed impedance tensor and returns a unique mathemat solution that, contrary to the invariants, is not affected by galvanic distortion [ Therefore, its analysis provides a distortion-free characterization of the strike dimensionality as well as offering additional information for the interpretation of underlying conductivity distribution.
The direction of the tensor and the three coordinate invariants (Φmin, Φmax, β) graphically represented by means of the tensor ellipse, which discloses the orientatio well as the lateral variations of the underlying electrical structures. The phase ten ellipses and their analysis for different periods are described in the Supplemen Material (part A). The result confirms a general 3D dimensionality with scarce 1D and

Dimensionality Analysis
The first evaluation of the raw dataset (before considering any static shift corrections or applying any rotation) involved a thorough analysis of the geoelectrical dimensionality by means of the rotational invariants of the impedance tensor [64,65], which gives an initial indication of the spatial distribution (1D, 2D, or 3D) of the electrical resistivity at depth and identifies the presence of galvanic distortion (3D/2D). Details are given in the Supplementary Material (Part A). The analysis results show a discrete number of 1D and 2D tensors for short periods and a predominant number of 3D cases, which clearly outlined the 3D nature of the subsurface electrical structures and thus the priority given to 3D inversion.
Further insight into the data came from the phase tensor analysis [66][67][68]. The phase tensor (Φ) is based on the observed impedance tensor and returns a unique mathematical solution that, contrary to the invariants, is not affected by galvanic distortion [66]. Therefore, its analysis provides a distortion-free characterization of the strike and dimensionality as well as offering additional information for the interpretation of the underlying conductivity distribution.
The direction of the tensor and the three coordinate invariants (Φ min , Φ max , β) are graphically represented by means of the tensor ellipse, which discloses the orientation as well as the lateral variations of the underlying electrical structures. The phase tensor ellipses and their analysis for different periods are described in the Supplementary Material (Part A). The result confirms a general 3D dimensionality with scarce 1D and 2D cases at short periods. The strike direction was calculated for the whole period range from Remote Sens. 2022, 14, 542 8 of 25 the impedance tensor Z (Figure 4a) [65], the phase tensor azimuth (Figure 4b) [66], and the induction arrows ( Figure 4c). The three rose-diagrams of the strike direction are in good agreement and define the direction of N130 • E for the geoelectrical strike. This strike direction was also estimated at each period decade (not reported here) and was consistent with that calculated considering all the periods. In detail, there was a clear strike direction in the range of 0.1-100 s because all the rose-diagram bins were aligned toward the preferred orientation of N130 • E. At periods higher than 100 s, a leading direction of around N150 • E was observed, but a small number of bins were oriented around N45 • E. good agreement and define the direction of N130°E for the geoelectrical strike. This strike direction was also estimated at each period decade (not reported here) and was consistent with that calculated considering all the periods. In detail, there was a clear strike direction in the range of 0.1-100 s because all the rose-diagram bins were aligned toward the preferred orientation of N130°E. At periods higher than 100 s, a leading direction of around N150°E was observed, but a small number of bins were oriented around N45°E. The geomagnetic transfer function, also known as the tipper, represents the ratio between the vertical and horizontal components of the magnetic field. Figure 5 shows the tipper vectors (or induction arrows) of the subset of 26 sites (out of 51) whose vertical transfer function has been measured. The color scale of the ellipses ranges from 0° to 90° and refers to the phase-tensor determinant. It deviates from an angle of 45° if there is no 1D regional distribution in the subsurface. For the selected periods of 1 and 10 s ( Figure  5a,b, respectively), the induction arrows present a substantial agreement with the strike direction shown in Figure 4 since they are orthogonal. We adopted the Parkinson criterion to represent the induction arrows, which point toward the region of highest conductance. The vectors suggest a possible NW-SE-oriented conductive region (in line with the phasetensor ellipses shown in Supplementary Material, part A).
The outcome of the dimensionality analysis of the three complete datasets is interestingly related to the geology of the area. Firstly, the direction of the geoelectrical strike corresponds to the NW-SE-oriented major tectonic structures and surface geological units, as shown in the map of Figure 2a. This points out the role of the NNW-SSE-oriented Radicondoli sedimentary basin, which is located to the northeast of the study area (sensu [59] or [69]). Secondly, the ellipses and the tippers at all periods seem to be sensitive to a possible large-scale deep structure oriented about N30-40°E along a certain strike that is transverse to the Apennines (see Supplementary Material, part A). The geomagnetic transfer function, also known as the tipper, represents the ratio between the vertical and horizontal components of the magnetic field. Figure 5 shows the tipper vectors (or induction arrows) of the subset of 26 sites (out of 51) whose vertical transfer function has been measured. The color scale of the ellipses ranges from 0 • to 90 • and refers to the phase-tensor determinant. It deviates from an angle of 45 • if there is no 1D regional distribution in the subsurface. For the selected periods of 1 and 10 s (Figure 5a,b, respectively), the induction arrows present a substantial agreement with the strike direction shown in Figure 4 since they are orthogonal. We adopted the Parkinson criterion to represent the induction arrows, which point toward the region of highest conductance. The vectors suggest a possible NW-SE-oriented conductive region (in line with the phase-tensor ellipses shown in Supplementary Material, Part A).
The outcome of the dimensionality analysis of the three complete datasets is interestingly related to the geology of the area. Firstly, the direction of the geoelectrical strike corresponds to the NW-SE-oriented major tectonic structures and surface geological units, as shown in the map of Figure 2a. This points out the role of the NNW-SSE-oriented Radicondoli sedimentary basin, which is located to the northeast of the study area (sensu [59] or [69]). Secondly, the ellipses and the tippers at all periods seem to be sensitive to a possible large-scale deep structure oriented about N30-40 • E along a certain strike that is transverse to the Apennines (see Supplementary Material, Part A).
The detailed analysis carried out in this section proves that, although a 2D strike direction of around N130 • E could be considered as a "first order approximation," as also suggested by the regional 2D trend of the surface geology, proper interpretation of our data requires 3D MT inversion. The detailed analysis carried out in this section proves that, although a 2D strike direction of around N130°E could be considered as a "first order approximation," as also suggested by the regional 2D trend of the surface geology, proper interpretation of our data requires 3D MT inversion.

Static Shift Correction by Means of TDEM Data
TDEM data have often been combined with MT data from geothermal areas even though their depth of investigation is usually much lower than that of the geothermal target (see for example [70]). The TDEM method is commonly adopted to correct the MT static shift because it provides independent information, which is not affected by telluric distortion (for a detailed description of the correction method see [71]). We carried out a TDEM survey to manage the static shift that could affect the MT apparent-resistivity (ρapp) curves.
We acquired new TDEM soundings in 2019 in order to constrain the MT soundings placed in different geological settings and to ensure a wide spatial coverage. The soundings corresponded to eight MT sites: a1, b2, b6, e1, g1, k1, k4, and k5 (see Figure 1 for their locations). TDEM data were acquired using the TEM-FAST 48 system by the AEMR company. The acquisition configuration was a coincident loop of 100 × 100 m or 75 × 75 m, depending on the accessibility of the site. The injected current was 3 A, the turnoff time was 7-8 μs and a total of 32-40 samples were acquired in the range of 4-4000 μs.
The TDEM data were adopted in order to directly correct the static shift of the correspondent MT site by means of 1D joint optimization and, if needed, to support the correction of the closest MT site only if it was placed on the same outcropping geological formation [29,72]. The static shift correction hence involved six TDEM-MT sites (k1 and k4 are not included in the 3D inversion), plus ten of the closest MT sites. A larger TDEM survey was not possible due to logistical issues.
To correct the static shift, TDEM and MT data were jointly inverted following the method of [73], based on particle swarm optimization (PSO). This metaheuristic approach performs the simultaneous 1D optimization of TDEM and MT data for each MT polarization. The algorithm optimizes both the model parameters, i.e., the resistivity model, and an additional parameter accounting for the static shift. The off-diagonal components of the impedance tensor of the selected MT sites were corrected for static shift in order to perform a 3D inversion of the corrected off-diagonal components to be compared with the inversion of the uncorrected off-diagonal components (see Section 5).

Static Shift Correction by Means of TDEM Data
TDEM data have often been combined with MT data from geothermal areas even though their depth of investigation is usually much lower than that of the geothermal target (see for example [70]). The TDEM method is commonly adopted to correct the MT static shift because it provides independent information, which is not affected by telluric distortion (for a detailed description of the correction method see [71]). We carried out a TDEM survey to manage the static shift that could affect the MT apparent-resistivity (ρ app ) curves.
We acquired new TDEM soundings in 2019 in order to constrain the MT soundings placed in different geological settings and to ensure a wide spatial coverage. The soundings corresponded to eight MT sites: a1, b2, b6, e1, g1, k1, k4, and k5 (see Figure 1 for their locations). TDEM data were acquired using the TEM-FAST 48 system by the AEMR company. The acquisition configuration was a coincident loop of 100 × 100 m or 75 × 75 m, depending on the accessibility of the site. The injected current was 3 A, the turn-off time was 7-8 µs and a total of 32-40 samples were acquired in the range of 4-4000 µs.
The TDEM data were adopted in order to directly correct the static shift of the correspondent MT site by means of 1D joint optimization and, if needed, to support the correction of the closest MT site only if it was placed on the same outcropping geological formation [29,72]. The static shift correction hence involved six TDEM-MT sites (k1 and k4 are not included in the 3D inversion), plus ten of the closest MT sites. A larger TDEM survey was not possible due to logistical issues.
To correct the static shift, TDEM and MT data were jointly inverted following the method of [73], based on particle swarm optimization (PSO). This metaheuristic approach performs the simultaneous 1D optimization of TDEM and MT data for each MT polarization. The algorithm optimizes both the model parameters, i.e., the resistivity model, and an additional parameter accounting for the static shift. The off-diagonal components of the impedance tensor of the selected MT sites were corrected for static shift in order to perform a 3D inversion of the corrected off-diagonal components to be compared with the inversion of the uncorrected off-diagonal components (see Section 5).
As an example, a typical outcome of the static shift correction using PSO is plotted in Figure 6 for xy and yx modes of site a1. Figure 6a shows the MT observed xy and yx apparent resistivity (red and blue dots, respectively) and the calculated response (red and blue crosses, respectively). There is a satisfactory overlap between the MT corrected curves and the TDEM data (black dots) at short periods. In the example, the multiplicative factors calculated for static-shift correction were S xy = 9 and S yx = 2. Figure 6b displays the good agreement between the observed and modeled phases. Figure 6c plots the diagonal components whose low magnitude is evident for the observed xx component (green dots) at long periods.
at long periods.
Regarding the other sites corrected for static shift, the majority of the original MT ρapp curves that, at the shortest periods, laid either above 100 Ω⸱m or below 10 Ω⸱m, were shifted within the range 10-100 Ω⸱m. In particular, the sites over the Neogene sediments (such as a1 and b2) had a corrected ρapp starting from around 100 Ω⸱m, while the sites on the Ligurian unit (such as e1 and d3) had a corrected ρapp starting from around 20-30 Ω⸱m. A detailed overview of the static-shift factors Sxy and Syx is depicted in Figure 7. The histograms show that the range of Sxy and Syx is 0.1-10 and most of them are lower than 2, thus meaning low correction factors. The grey bins in Figure 7 represent the ratio Sxy/Syx, which is generally lower than 1 and represents the level of correction between the two modes.  Regarding the other sites corrected for static shift, the majority of the original MT ρ app curves that, at the shortest periods, laid either above 100 Ω·m or below 10 Ω·m, were shifted within the range 10-100 Ω·m. In particular, the sites over the Neogene sediments (such as a1 and b2) had a corrected ρ app starting from around 100 Ω·m, while the sites on the Ligurian unit (such as e1 and d3) had a corrected ρ app starting from around 20-30 Ω·m. A detailed overview of the static-shift factors S xy and S yx is depicted in Figure 7. The histograms show that the range of S xy and S yx is 0.1-10 and most of them are lower than 2, thus meaning low correction factors. The grey bins in Figure 7 represent the ratio S xy /S yx , which is generally lower than 1 and represents the level of correction between the two modes.

3D MT Modeling and Inversion
The final dataset addressed by the inversion was composed of the 51 sites plotted as red squares on the geological map in Figure 2a. The complete period range was 0.003-

3D MT Modeling and Inversion
The final dataset addressed by the inversion was composed of the 51 sites plotted as red squares on the geological map in Figure 2a. The complete period range was 0.003-1000 s, originally discretized into 75 values, then resampled into 20 values to unburden the timeconsuming 3D computation. The resampling was calculated using a routine in the MTpy toolbox [67].
The 3D inversion was performed using the ModEM software, which is available for the EM research community [9,15]. The inversion scheme of ModEM is based on the nonlinear conjugate gradient (NLCG) and the computation is parallelized. The objective function of the linearized inversion minimizes both the data misfit and model roughness, penalizing the deviations from the starting model. Inversion with different settings were executed using 96 cores (three nodes) of the high-performance computing (HPC) cluster of the Politecnico di Torino. The CPU model of the single node was an Intel Xeon Gold 6130 2.10 GHz with 21 TB (DDR4) of RAM. When the runs were computed, the sustained performance of the cluster was 21 TFLOPS. The computation time was around one hour per NLCG iteration.
The inversion settings, the 3D mesh, and the results were prepared and analyzed using the 3D-GRID Academic software. The domain was 130 km in the horizontal directions away from the model center and around 300 km deep, consistent with the boundary conditions and the electromagnetic skin depth. The mesh included the topography and bathymetry and the resistivity of the sea was fixed to 0.3 Ω·m. One important advantage of including the topography in the 3D model is that the forward computation compensates the static shift caused by topography [70,74,75]. Along the horizontal directions, the model mesh was discretized into 100 × 100 cells, whose size was constant (140 m) for the central 70 × 70 cells, linearly increasing by a factor of 1.37 for the external cells (15 planes for each lateral side towards the boundary of the domain). The vertical direction was discretized into 65 layers, whose thickness was 28 m for the first layer, then increasing by a factor of 1.2.
We performed an extensive number of inversion tests in order to find the model (affected by equivalence) that has the best data fitting and agreement with the complex geology of the Travale geothermal field. The settings of the different inversion tests are summarized here, and fully explained in the Supplementary Material (Part B). These settings regard the starting model, smoothing factors, mesh orientation, tensor(s) to be inverted, and correction of the static shift for some sites. These settings were chosen following the most recent literature findings but also adapted to our specific dataset.
The starting model used for the inversion was a homogeneous half space of 100 Ω·m, which was selected after testing four different resistivity values (10, 100, 300, and 1000 Ω·m). Figure S3 in Supplementary Material shows that the lowest RMSE (1.05) was achieved with the fewest number of iterations by the inversion test using a 10-Ω·m homogeneous starting model. However, the 100-Ω·m resistivity value yielded the second lowest RMSE (1.15) and, since it was coherent with the average apparent resistivity of the sites after static shift correction, it was chosen for the starting model. We also performed a further inversion test with an a priori starting model derived from the 2D inversion result of [29]. It is provided in Supplementary Material (Part B).
The smoothing factor controls the model regularization along the three dimensions. Its choice was as critical as that of the starting model because of the over-parametrization of the inverse problem, which can lead to different outcomes [76]. After some trials, we set the smoothing factor to 0.2 for the horizontal directions and 0.1 for the vertical direction, in order to emphasize the vertical contrasts among the deep structures. When the inversion was initialized by the a priori starting model (test C in Supplementary Material, Part B), the smoothing factors were increased to 0.5 and 0.4, respectively, in order to test, on one hand, if the a priori resistivity contrast could be smoothed-despite the ModEM approach of minimal deviations from the starting model-and, on the other hand, if it was in fact required by the data.
The coordinate system was aligned with the geoelectrical strike. Even though the coherence between the strike direction and the rotation of the model mesh and data is fundamental in 2D modeling, it has turned out to be critical in 3D modeling as well [76][77][78]. The authors of [76] demonstrated that the 3D inversion result is dependent on the coordinate system and recommended a strike-oriented model mesh since it mostly minimizes the coupling among the four components of the impedance tensor (which, in ModEM, are handled independently). The MT data (Z and T) were rotated to N50 • W (which is the same rotation as N130 • E). The location of the MT sites was rotated to N50 • E with respect to the mesh in order to simulate a mesh rotation of N50 • W, because ModEM assumes that the data and the mesh are rotated in the same direction. We also performed two inversion tests (tests A and B in Supplementary Material, Part B) with the coordinate system aligned to the geographic north (N0 • ) and the MT data (Z and T) not rotated.
The uncorrected off-diagonal components of Z were inverted together with the tipper data. Even though in a 3D environment the information about the subsurface resistivity distribution is included in all the tensor elements, which can improve the spatial resolution, we neglected the inversion of the on-diagonal components since they had a relatively high magnitude. The inversion of the tipper T was included to improve the sensitivity at depth and to seek out lateral constraints [10]. After evaluating the quality of our data for the inversion, we selected 24 out of 26 sites with the tipper.
For the sake of completeness, the other inversion tests performed are supplied in Supplementary Material (Part B). They encompass the inversion of: uncorrected Z and T (tests A, C, D), uncorrected Z (test B), uncorrected off-diagonal components of Z with (test E) and without the tipper (test F), static-shift corrected off-diagonal components of Z (test G) considering the sites corrected for static shift (see Section 3.3).
The error floor was set as a portion of the absolute value of the impedance components (|Z ij |) instead of as a portion of the mean of the off-diagonal components (|Z xy · Z xy | 0.5 ), because the mean value could have underestimated one component of the tensor with respect to the other ones [76]. Since the original errors of the data were not negligible, the error associated with the data was set as the maximum between the original error and an error floor of 10% for the off-diagonal components of Z, and of 15% for the on-diagonal components of Z, which generally presented high magnitudes. The error floor of the tipper components was assigned by taking into account the error ratio between the impedance and the geomagnetic transfer function [79], that is, equal to the error floor set for the impedance, i.e., 10% of the absolute value.

Inversion Result
The result we present in this section is the outcome of the inversion performed on the strike-aligned mesh starting with a 100 Ω·m homogeneous model and using the offdiagonal components of Z and T. The error floor was 10% for Z xy , Z yx , and the T. We consider this outcome to be the reference model among the other tests performed for two reasons: it is the most representative model of the geology of the area and it has the best data fitting for the off-diagonals at long periods (the calculated responses for selected sites are shown in Supplementary Material, Part C). In fact, the exclusion of the diagonals improved the data fitting with little effect on the model structures thanks to the inclusion of the tipper data. Furthermore, this inversion did not consider the static shift corrections, but given that we had rotated the data, and that only information from TDEM was available at some sites, it would have added doubts if this correction had been applied to the proper components of the tensor.
The inversion ended when the RMSE did not further decrease (and the minimum damping factor was reached). The total number of iterations was 71 and the final RMSE was 1.12.
The 3D resistivity model obtained is shown in Figures 8-10. (The figures have not been smoothed, but represent the original cell-by-cell visualization.) Figure 8 shows the resistivity distribution in the horizontal plane at six different depths (27 m a.s.l., 222 m, 722 m, 1.4 km, 3.1 km, and 5.3 km b.s.l.) as they reveal the most evident resistivity contrasts of the imaged model. In the shallow structures (Figure 8a,b), we observe a clear difference between the northeastern conductive region roughly oriented NW-SE (<10 Ω·m), hereafter C1, and the southwestern region, hereafter R1, where the resistivity rises by up to 300 Ω·m. From 1 to 2.5 km b.s.l., the resistivity increases up to 180 Ω·m in a resistive body (R2 in Figure 8d). Below a depth of 3 km (Figure 8e,f), a resistive body of maximum 140 Ω·m, hereafter R3, elongates orthogonally to the strike (i.e., in parallel to the y-axis).  [29]. The SW-NE profile is orthogonal to the strike direction and crosses sites from k5 to a8 (see AA′ in Figure 8a). The vertical dotted lines are the traces of the multi-segment section. Figure 9. Vertical cross-section of the model corresponding to the MT profile investigated by [29]. The SW-NE profile is orthogonal to the strike direction and crosses sites from k5 to a8 (see AA in Figure 8a). The vertical dotted lines are the traces of the multi-segment section.
The vertical cross-section depicted in Figure 9 is directly comparable with the 2D model of [29]. The most important feature between the depths of 3 and 8 km is a large 140-Ω·m body (R3). The R3 body extends orthogonally to the strike direction, as imaged also in the horizontal view of Figure 8e-f. Figure 10 shows the vertical cross-sections along the Z Y and Z X planes (the traces of these sections are marked in Figure 8b).   The RMSE of 1.12 can be appreciated from Figure 11a,b, which plots the distribution of the RMSE in the frequency-space domain for Z and T, respectively. The final RMSE normalized for the full period bandwidth was 0.87 for Z and 1.5 for T. The lowest RMSE for Z was measured in site g1 (0.5) and for T in site b8 (0.6). The worst RMSE resulted in site e1 (1.4) for Z and in site c1 for T (3.4).

Discussion
We have explored different inversion setups in this study because 3D inversion in geothermal areas is potentially challenging, mainly due to the complex geology. The strike-aligned approach and the exclusion of the on-diagonals of Z had the advantage of improving the data fitting at long periods.
The xy-plane view of the resistivity distribution at shallow depth in Figure 8a is in agreement with the expected resistivity of the outcropping rocks described in the geological map reported in Figure 2a. The northeastern conductive area (C1) mostly corresponds to the Neogene sediments and the Ligurian complex, and was imaged up to a depth of around 1200 m b.s.l. (see Figures 8a and 9). The central mildly resistive area (10-50 Ω⸱m) between C1 and R1 in the horizontal view of Figure 8b corresponds to the spatial coverage of the Ligurian and sub-Ligurian flysch complex. Finally, the southern region R1, showing resistivities up to 300 Ω⸱m, matches the sediments and carbonates of the Tuscan nappe. It should be noted, however, that the irregular space-covering of the sites may have influenced the shape of the imaged resistive structures, as can be seen in Figure 8. Figure 8e,f displays the loss of resolution with depth, which can be explained by two main reasons. The first reason is related to the dataset, because the spatial aperture of the MT survey, which should be twice the intended depth of investigation, prevented a deep geometrical constraint of the 3D structures [80]. The second reason regards the imaged shallow conductors, such as C1, which reduced the maximum skin depth calculated from the impedance determinant to a few tens of kilometers (instead of several hundreds of kilometers for the sites far from C1).
Interestingly, we found a resistivity contrast in correspondence with the bottom of the Ligurian units (and Neogene sediments) of the Radicondoli basin, separating them from the underlying units piled in the chain (from [33]). Figure 12 shows a north-south section of the model crossing sites from a4 to j3. The geological surface (black-dashed line) was processed in Petrel software from the findings of [33] and represents the base of the cap-rock (Neogene and Ligurian units). To the north (below site b4) the basin-bottom surface descends from a depth of around 600 to 1500 m b.s.l., where the resistivity jumps from 4 Ω⸱m to 40 Ω⸱m. This represents a remarkable agreement between a 3D resistivity contrast and a geological surface from the known geological model.

Discussion
We have explored different inversion setups in this study because 3D inversion in geothermal areas is potentially challenging, mainly due to the complex geology. The strike-aligned approach and the exclusion of the on-diagonals of Z had the advantage of improving the data fitting at long periods.
The xy-plane view of the resistivity distribution at shallow depth in Figure 8a is in agreement with the expected resistivity of the outcropping rocks described in the geological map reported in Figure 2a. The northeastern conductive area (C1) mostly corresponds to the Neogene sediments and the Ligurian complex, and was imaged up to a depth of around 1200 m b.s.l. (see Figures 8a and 9). The central mildly resistive area (10-50 Ω·m) between C1 and R1 in the horizontal view of Figure 8b corresponds to the spatial coverage of the Ligurian and sub-Ligurian flysch complex. Finally, the southern region R1, showing resistivities up to 300 Ω·m, matches the sediments and carbonates of the Tuscan nappe. It should be noted, however, that the irregular space-covering of the sites may have influenced the shape of the imaged resistive structures, as can be seen in Figure 8. Figure 8e,f displays the loss of resolution with depth, which can be explained by two main reasons. The first reason is related to the dataset, because the spatial aperture of the MT survey, which should be twice the intended depth of investigation, prevented a deep geometrical constraint of the 3D structures [80]. The second reason regards the imaged shallow conductors, such as C1, which reduced the maximum skin depth calculated from the impedance determinant to a few tens of kilometers (instead of several hundreds of kilometers for the sites far from C1).
Interestingly, we found a resistivity contrast in correspondence with the bottom of the Ligurian units (and Neogene sediments) of the Radicondoli basin, separating them from the underlying units piled in the chain (from [33]). Figure 12 shows a north-south section of the model crossing sites from a4 to j3. The geological surface (black-dashed line) was processed in Petrel software from the findings of [33] and represents the base of the cap-rock (Neogene and Ligurian units). To the north (below site b4) the basin-bottom surface descends from a depth of around 600 to 1500 m b.s.l., where the resistivity jumps from 4 Ω·m to 40 Ω·m. This represents a remarkable agreement between a 3D resistivity contrast and a geological surface from the known geological model. resistive region R2 (180 Ω⸱m) rather than a conductive anomaly as seen in the 2D model. The resistive body R3 at 3-8 km b.s.l. has a different shape, boundaries, and resistivity values compared to the 2D model. Given all these dissimilarities, our result overcame the main limitations of the 2D model, which could have been strongly biased by its chosen a priori geological cross-section starting model. Moreover, the 2D model might have missed some heterogeneities due to the adoption of a 2D approach for the inversion of 3D data. The volume R2, with more than 180 Ω⸱m at a depth of about 1-2.5 km, is embedded between the low-resistive cap-rock (C1) and the 140 Ω⸱m R3 body. The interpretation of R2 is challenging due to the interplay of multiple processes. At this depth, heterogeneous units (tectonic wedge complex (TWC) and phyllites) mostly occur across a very wide range of electrical resistivities (from 10 −1 Ω⸱m to 10 4 Ω⸱m) measured in the geophysical well logs, such as the wells Radicondoli 7bis and Montieri 4 in the area of interest [29]. Figure  13 shows a comparison between the well log resistivity along the Radicondoli 7bis and the correspondent 1D resistivity profile extracted from our 3D MT model. At the depth of the TWC and metamorphic units (below about 750 m), an averaged MT response is obtained considering the higher order of magnitude of the investigated volume with respect to the resistivity log that recorded a huge variability of the resistivity. At shallow depths (above 750 m), the resistivity from the model was higher than that from the log. In this case, the shallow cap-rock is not properly imaged, probably because the well is located at the The vertical section of the 3D model in Figure 9 crosses the sites belonging to the 2D profile investigated by [29]. The resistivity distribution of the 3D model is partially similar to that of the 2D model, but some features emerged from the 3D inversion. The shallow structure R1 (up to 1 km b.s.l.) beneath sites k5-j2 is far more resistive (300 Ω·m) than in the 2D model. Beneath sites g4-a8, from 1 to around 3 km b.s.l., the 3D model has a resistive region R2 (180 Ω·m) rather than a conductive anomaly as seen in the 2D model. The resistive body R3 at 3-8 km b.s.l. has a different shape, boundaries, and resistivity values compared to the 2D model. Given all these dissimilarities, our result overcame the main limitations of the 2D model, which could have been strongly biased by its chosen a priori geological cross-section starting model. Moreover, the 2D model might have missed some heterogeneities due to the adoption of a 2D approach for the inversion of 3D data.
The volume R2, with more than 180 Ω·m at a depth of about 1-2.5 km, is embedded between the low-resistive cap-rock (C1) and the 140 Ω·m R3 body. The interpretation of R2 is challenging due to the interplay of multiple processes. At this depth, heterogeneous units (tectonic wedge complex (TWC) and phyllites) mostly occur across a very wide range of electrical resistivities (from 10 −1 Ω·m to 10 4 Ω·m) measured in the geophysical well logs, such as the wells Radicondoli 7bis and Montieri 4 in the area of interest [29]. Figure 13 shows a comparison between the well log resistivity along the Radicondoli 7bis and the correspondent 1D resistivity profile extracted from our 3D MT model. At the depth of the TWC and metamorphic units (below about 750 m), an averaged MT response is obtained considering the higher order of magnitude of the investigated volume with respect to the resistivity log that recorded a huge variability of the resistivity. At shallow depths (above 750 m), the resistivity from the model was higher than that from the log. In this case, the shallow cap-rock is not properly imaged, probably because the well is located at the boundary of the area covered by the MT sites (see Figure 2a). There, the model resistivity keeps the starting value of 100 Ω·m, while where the model is better constrained, the cap-rock is properly imaged with a lower resistivity. Even though highly resistive volumes were measured in the heterogeneous units (TWC and phyllites), the embedded extremely low resistivity values (probably owing to graphite) have affected the measured EM signal, which provided an averaged resistivity of this geological unit. This phenomenon is known to occur in the Larderello geothermal system following a case where research was conducted in another sector of the field in the proximity of Lago Boracifero by means of EM measurements and an experimental Surface-Hole Deep ERT [35,44,81]. In the central part of the study area, the 3D model is characterized by a well-defined volume of increased resistivity (R2), which elongates transversely to the Apennines (see Figures 8e,f and 13). A factor that can play a role in the increased resistivity is the occurrence of the vapor-dominated hydrothermal circulation. Indeed, part of the exploited reservoir is located at the same depth as R2. We point out the match between R2 and the current geothermal target represented by the seismic H-horizon along the seismic profile from [39] (see Figure 15). It is possible that in other sectors of the investigated area the H-horizon can lie in correspondence with R3. Further investigations with complete seismic datasets would be beneficial. boundary of the area covered by the MT sites (see Figure 2a). There, the model resistivity keeps the starting value of 100 Ω⸱m, while where the model is better constrained, the caprock is properly imaged with a lower resistivity. Even though highly resistive volumes were measured in the heterogeneous units (TWC and phyllites), the embedded extremely low resistivity values (probably owing to graphite) have affected the measured EM signal, which provided an averaged resistivity of this geological unit. This phenomenon is known to occur in the Larderello geothermal system following a case where research was conducted in another sector of the field in the proximity of Lago Boracifero by means of EM measurements and an experimental Surface-Hole Deep ERT [35,44,81]. In the central part of the study area, the 3D model is characterized by a well-defined volume of increased resistivity (R2), which elongates transversely to the Apennines (see Figures 8ef and 13). A factor that can play a role in the increased resistivity is the occurrence of the vapor-dominated hydrothermal circulation. Indeed, part of the exploited reservoir is located at the same depth as R2. We point out the match between R2 and the current geothermal target represented by the seismic H-horizon along the seismic profile from [39] (see Figure 15). It is possible that in other sectors of the investigated area the Hhorizon can lie in correspondence with R3. Further investigations with complete seismic datasets would be beneficial.  The deep resistive body R3 (see Figure 8e-f, Figures 9 and 10) was imaged in all the inversion tests and showed high resistivities (>140 Ω·m) between the depths of 3 and 8 km. Figure 14 displays the 3D model with a cutaway view on the 160 Ω·m isosurface. Further inspection of the 3D model is supplied in Supplementary Material with a brief animation. For the first time the deep resistive body of the Travale geothermal field, already known from 2D inversion, is revealed with its position and orientation in three dimensions. From Figure 14, the R3 body is oriented orthogonally to the main Apennine strike direction. This was expected because some of the soundings showed a strike direction of N40-50 • E at long periods (see Supplementary Material, Part A). The orientation of R2 and R3 represents a major improvement in the knowledge of the Travale geothermal field. A strict relation between the strike-slip tectonics transversal to the Apennine direction and fluid circulation was recently claimed by various studies [59,[81][82][83].  To assess if R3 was necessary to fit the data, we carried out a sensitivity analysis (see Supplementary Material, Part C, for details). R3 was replaced with a perturbed body of 1, 20, 50, and 100 Ω·m in four different sensitivity tests, respectively. The 1 Ω·m perturbed model most significantly worsened the RMSE. The recalculation of the forward problem led to an overall RMSE increase of 75% and to an RMSE increase from 57% to 101% for the sites directly above the body (see Figure S17 in the Supplementary Material for a representative site). The sensitivity analysis proved that the model is most sensitive to the lowest resistivity value of the perturbed body (1 Ω·m) and that the data are incompatible with a resistivity lower than 50-100 Ω·m for R3, given the limitations of the poor resolution at this depth.
In order to interpret the resistive body R3 below the Travale field, we compared the 3D resistivity model with the 3D local earthquake tomography by [39]. An analogy occurs between the resistive body R3 and the low-velocity body detected in the v p images below the Travale area along an available profile (see Figure 15 from profile CC plotted in Figure 8a). The velocity model reported a 5-km-wide low-velocity body (v p about fragile regime and rock fracturing. It is worth noting from [39] that the largest hypocenter cluster corresponds to a high-velocity anomaly below the K-horizon, a region which is unfortunately not covered by our MT survey. It should be noted that, in [39], the trend of the low-velocity body is different from that of our R3 body, but some tests performed in [39] also showed a transverse NE-SW trend in the area.  A vertical cross-section drawn from the 3D resistivity model is compared with the 3D seismic velocity model from local earthquake tomography by [39]. The R3 body (>140 Ω⸱m) and the Figure 15. A vertical cross-section drawn from the 3D resistivity model is compared with the 3D seismic velocity model from local earthquake tomography by [39]. The R3 body (>140 Ω·m) and the low-velocity body with v p around 5-5.5 km/s both extend between the depths of 3 and 8 km. The trace of the cross-section is the profile CC in Figure 8a and it corresponds to the seismic profile CC of Figure 6 in [39]. The dashed lines represent the absolute seismic P-wave velocity (from 4.5 to 6 km/s). The red contour lines are the isodepths of the H-and K-horizons (from [39] and references therein). The black circles represent the hypocenters and are proportional to the magnitude. The coordinates of the section edges are (latitude north-longitude east): 43.1667 • -11.0240 • (left) and 43.2124 • -11.0732 • (right).
Our interpretation of R3 is that it may be due to the effect of the interplay of a highly fractured volume of crystalline rocks hosting a high-temperature hydrothermal circulation and to the occurrence of a very large (crystallized) granitic intrusion, locally cored by some wells at depths between about 2-4.5 km. The resistive nature of the deep body is explained by the vapor-dominated condition of the reservoir, so that highly resistive fluids (vapors) circulate in highly resistive rocks (e.g., granites). Moreover, a pervasive mineral alteration, which would decrease the bulk resistivity of the rocks, has been excluded by well-cuttings and core-analyses carried out in this area as part of the I-GET project [40,84].
Our 3D model had limited resolution below R3, but showed an interesting relationship with the seismic marker K-horizon. The bottom of R3 is around 8 km deep (see Figures 9,10 and 15) and could be reasonably associated with the depth of the K-horizon in this area (following [47]). We did not recognize the occurrence of melt fractions in the area of study down to 8 km. The occurrence of melted igneous intrusions cannot be excluded at a level deeper than 8 km or in the Travale area, which is not covered by our MT sites.

Conclusions
This work presents the first 3D resistivity model of the geothermal area of Travale (Italy) resulting from derivative-based 3D MT inversion.
The MT dataset was accurately analyzed in terms of geoelectrical dimensionality, strike direction, and phase tensor properties.
Given the equivalence of the solutions of the 3D MT inverse problem, we performed an extensive number of inversion tests by varying the initial parameters. These regarded the starting model, the data to be inverted, the grid rotation, and the correction of the static shift for some selected sites. The static shift was corrected for the off-diagonal components of the impedance tensor by means of TDEM soundings for those MT sites in correspondence with (or in the vicinity of) the TDEM acquisitions. After many tests, we have demonstrated that-if static-shift corrected data of the selected sites are included in the 3D modeling-the shallow structures are more homogeneous and the long period data fit better than for the test with uncorrected data (see Supplementary Material), but there was no appreciable improvement in the final RMSE and model. For these reasons, the 3D model for the final interpretation was selected from the inversion of uncorrected off-diagonal components and the tipper, using a strike-aligned mesh. This model has the best data fitting and agreement with the complex geology of the Travale geothermal system.
The selected model represents an important contribution to the investigation of the Travale geothermal field. The most important outcome is that the final model identifies and characterizes the region's 3D subsurface structures. In this way, our results have expanded the knowledge of the subsurface by adding another spatial dimension, where until now it had only been estimated by means of 2D MT inversion. The validity of our results is supported by geological information, resistivity well logs, and seismic data: -A distinct correlation emerged between the resistivity contrast at shallow depths and the geological surface of the base of the Neogene sedimentary basin; -Two deep resistive bodies were imaged at depths of 1-3 km (R2) and 3-8 km (R3) with resistivities higher than 180 Ω·m and 140 Ω·m, respectively, with a N40-50 • E orientation; - The mildly resistive body R2 lies in a region where the geophysical well logs measured a heterogeneous resistivity value (10 −1 Ω·m-10 4 Ω·m); -A marked analogy was identified between the deep resistive body R3 and a lowvelocity body (v p about 5 km/s) deeply rooted in the crust below the Travale area; - The role played by the vapor-dominated circulation is recognized in these highresistivity bodies (R2 and R3), together with the occurrence of (crystallized) granitic intrusions contributing to R3.
Three-dimensional MT inversion was challenging due also to the time-consuming computation which can require 100 times more unknowns than a 2D inversion. Furthermore, the geology of the area is extremely complex and the framework of the geothermal system is still under debate. Our work offers the first interpretation of the 3D subsurface electrical distribution in the Travale area from a geophysical perspective and could pave the way to further insight about the geothermal system. Future work should broaden the MT characterization of the Travale area by means of new acquisition campaigns that would enlarge the investigated zone, ideally with a regularly spaced coverage of the sites, and would moreover improve the data quality according