# Geophysical models integration using principal component analysis: application to unconventional geothermal exploration

Jean-Michel Ars, Pascal Tarits, Sophie Hautot, Mathieu Bellanger

► To cite this version:

Jean-Michel Ars, Pascal Tarits, Sophie Hautot, Mathieu Bellanger. Geophysical models integration using principal component analysis: application to unconventional geothermal exploration. Geophysical Journal International, 2024, 239 (3), pp.1789-1798. 10.1093/gji/ggae357. hal-05323381

HAL Id: hal-05323381

https://brgm.hal.science/hal-05323381v1

Submitted on 21 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geophysical models integration using principal component analysis: application to unconventional geothermal exploration

Jean-Michel Ars,<sup>1,*</sup> Pascal Tarits,<sup>1,2</sup> Sophie Hautot<sup>2</sup> and Mathieu Bellanger<sup>3</sup>

Accepted 2024 September 28. Received 2024 September 24; in original form 2024 April 16

## SUMMARY

Geophysical data integration covers a wide range of approaches, from visual interpretation of model presented side by side to sophistical statistical analyses such as automatic clustering. We present here a geophysical model integration based on principal component analysis (PCA), which allows to gain insight in a multivariable system. PCA for geophysical models integration define a new set of principal component (PC) models, distributed along new orthogonal axes by solving a eigenvalue problem. We show that PC component models patterns reflect different processes sensed by the geophysical observations. We applied this integration method to models obtained from constrained and joint inversion of gravity, ambient noise and MT data in the framework of unconventional geothermal exploration in Massif Central, France. PCA of the log-resistivity, the density contrast and the Vs velocity model has three independent components. The first one (PC1) representing 69 per cent of the total variance of the system is highly influenced by the parameter coupling enforced in the joint inversion process. PC1 allows to point to geophysical structures that may be related to the geothermal anomaly. The second component (PC2) represents 22 per cent of the total variance and is strongly correlated to the resistivity distribution. PC2 correlation to shallow fault structures suggests that it may be a marker of fracturing. The third component, PC3, accounts for 9 per cent of the total variance and is correlated with velocity structures and anticorrelated with density structures, respectively. The contribution of geophysical properties, mainly sensitive to the elastic properties of rock units, and the good agreement with shallow geological features make PC3 useful for the 3-D description of geological units. This statistical approach helps the interpretation of geophysical models into a limited number of geological processes, one possibly geothermal. Lithology may be derived from PC3, fracturing from PC2 and hydrothermal anomaly from PC1. Nevertheless PCA as a geophysical integration methodology is site-dependent and interpretation relies on a priori knowledge of the local geological mechanisms.

**Key words:** Gravity anomalies and Earth structure; Magnetotellurics; Joint inversion; Statistical methods; Hydrothermal systems; Seismic tomography.

## 1 INTRODUCTION

Geophysical exploration of natural resources relies mainly on imaging characteristic of the geological target in order to assess for its depth, geometry, etc. Natural resources can be found in very different geological contexts, and the targets can therefore present a variety of geophysical signatures depending on the imaging methods used. It is the combination of complementary information from

several geophysical techniques that allows specificity of each site. Parallel exploration using several geophysical imaging method is a successful approach that has been applied in many different geophysical exploration contexts (Power et al. 2004; Gola et al. 2017). Such multiphysics strategy may be improved by processing the different data sets together using a joint inversion approach. One can adopt a joint inversion strategy based on geometrical coupling between models which implies the least assumption on the local geology. One classical structural approach is the cross-gradient coupling term which enforces models to exhibit similar structures (Gallardo & Meju 2007). Another solution is the use of empirical

<footnote n="1"><sup>1</sup>Intitut Universitaire European de la Mer, UMR 6538 Laboratoire Geo-Ocean, Plouzané, F29280, France. E-mail: jean-michel.ars@geolinkservice.com</footnote>

<footnote n="2"><sup>2</sup>IMAGIR SARL, Saint-Renan, F29290, France</footnote>

<footnote n="3"><sup>3</sup>TLS GEOTHERMICS SAS, Toulouse, F31200, France</footnote>

<footnote>*Now at: GEOLINKS SERVICE SAS, Paris 75000, France.</footnote>
petrophysics relationship (Jegen et al. 2009; Coutant et al. 2012) to a specific exploration scenario allowing a robust constraint on models for a specific range of rocks. Joint inversion allows to image the underground with several geophysical models whose consistency would be enforced throughout the joint inversion process using coupling terms. Coupled models get closer to ‘True Earth’ model avoiding biased results due to inversion issues (non-unicity, local minima...; Meju & Gallardo 2016).

Whether with a geometrical or petrophysical coupling, most of the similarity measurements between parameters used for joint inversion make the assumption that models display similar interfaces, but it is likely to be wrong. First of all, coupled inversion cannot completely overcome the difference of resolution, the difference of parametrization and the irregular data distribution between methods. Geophysical imaging methods have different sensitivity so that geophysical signatures of the same resource can shape very different patterns (Coppo et al. 2015; Ars et al. 2019) especially around fractured bedrock (Schill et al. 2010). To better understand the geological target, it is important to be aware of the mutual information shared between geophysical signatures from different models. Jointly inverted or not, geophysical models from multiphysics data sets need to be integrated to interpret them properly into a reliable conceptual geological model. The difficulty in the multiphysics interpretation of integrated models is that each data obtained in a complex geological environment with several geological processes is the superimposition of the geophysical signature of each process. A step further is therefore to extract, from jointly inverted multiphysics models, the signature of each process.

Geophysical data integration spreads across a wide range of techniques, from qualitative approaches mainly relying on visual interpretation (Kvamme 2006) to dedicated strategy such as Play and Fairway Analysis (Craig 2018) and more complex statistical analysis using either cluster analysis (Paasche et al. 2006), probability density function (Muñoz et al. 2010) or self-organized map (Bauer et al. 2012). We present here another statistical approach for multivariate analysis called principal component analysis (PCA) or Empirical Orthogonal Function (EOF) (Hannachi et al. 2007). This approach provides enlightenment on relationships and variance between variables by recombining the total variance of the system among new principal components. It can not only be used for dimensionality reduction but also for pattern extraction that maybe related to different processes originally mixed into each independent geophysical data set (Hannachi et al. 2007). It has been widely used in geochemistry (Albarede & Pyle 1996), atmospheric science (Hannachi et al. 2007), sociology (Widaman 1993) and other sciences based on multivariate integration. While a large range of statistical geophysical integration has been investigated, PCA for geophysics is quite unknown. It has been successfully applied for geological unit characterization from borehole multiphysic data sets (Kassenaar 1991; Barrash & Morin 1997) or airborne geophysics (Ranjbar et al. 2001) but, to our knowledge, not for geophysical model integration.

Ars et al. (2019) developed a joint inversion technique of 3-D gravity and ambient noise tomography (ANT) constrained by a resistivity model dedicated to deep geothermal exploration. The three geophysical methods provide complementary information on the geology. Magnetotellurics (MT), which resolves resistivity distribution, is sensitive to fluid and mineral phases. ANT, which resolves shear wave velocity, is sensitive to structural properties and fluid content. Finally, gravity is used to map density distribution, mainly related to lithology. Resistivity, density and $V_s$ velocity parameters

are coupled through linear correlation laws. Instead of assuming petrophysical relationships between geophysical parameters, we try to enforce their correlation at model scale. This approach was first applied for the geothermal exploration along a Westward-Eastward strip (red dashed rectangle in Fig. 2) in a crustal fault zone located in Massif Central, France (Ars et al. 2019). Here we focus on the potential of the joint inversion approach for regional 3-D geophysical geothermal exploration using dense multiphysics data sets. We applied the joint inversion methodology in the same survey area located in the French Massif-Central (Fig. 2) but with a larger number of data, enhancing the resolution and the model extension. Geophysical imaging with high resolution data set are increasingly used for non-conventional geothermal exploration, but only with a single imaging method (Wannamaker et al. 2007; Lehujeur et al. 2018). Such resources require broad data coverage to understand the regional geology with enough resolution to image narrow structures such as faults which are not always obvious to characterize with a single imaging method. Using dense multiphysic data sets, we aim to gain in resolution of the 3-D geometry of a major fault zone (MFZ), believed to be permeable enough for deep fluid circulation (Roche et al. 2017), and the 3-D geometry of the local geology mainly composed of a granitic body and a metamorphic unit facing each other on both sides of the fault (see Fig. 2, Bellanger et al. 2019; Duwiquet et al. 2019).

We first describe PCA methodology, then the geophysical data set used and the imaging methodology using joint inversion approach. Following, we focus on geophysical model integration using PCA and how it allows to highlight local geology and potential geothermal features.

## 2 PRINCIPAL COMPONENT ANALYSIS FORMULATION

Principal component analysis aims at describing the covariance of several sets of multivariate parameters with a smaller number of parameters. The operation consists of determining the orientation of a new set of orthogonal axes, or components, in an $N$-dimensional space in order to optimize the variance of a given system of $N$ variables along those new components. For geophysical model interpretation, the given $N$-dimensional system is composed of $N$ geophysical models. Since geophysical parameters maybe of different nature, it is convenient to define reduced parameters (zero mean and normalized by their standard deviations) to be able to compare the combination of them. Let assume a matrix $X$ composed of $N$ reduced geophysical models of $M$ parameters each, so $X$ shape is $(M \times N)$. PCA is thus an eigenvalue problem where the eigenvalues are the variance of the covariance matrix of $X$ along the eigenvectors controlling the $N$ new components (see example of PCA in a 2-D space in Fig 1). Following Albarede & Pyle (1996), the covariance matrix $S$ of the reduced model matrix $X$ is the correlation matrix of the original set of reduced models. The eigenvalues problem can be written as:

$$S = U\Lambda U^T \quad (1)$$

$\Lambda$ is a diagonal matrix ($N \times N$) where the diagonal elements are the $N$ eigenvalues ($\lambda_1, \lambda_2, \dots, \lambda_n$) and $U$ the matrix of the $N$ orthogonal eigenvectors ($u_1, u_2, \dots, u_n$). We call $k^{th}$ PC component the projection of the reduced geophysical models matrix $X$ onto the $k^{th}$ eigenvector. The variance of the $k^{th}$ PC is the $k^{th}$ eigenvalue $\lambda_k$. We may also obtain information on the robustness on the $N$ PC components by expressing their uncertainty with eq. (2) (North
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Principal component analysis applied to a 2 geophysical models system. PC1 explains 78 per cent of the total variance of the system while PC2 explains 22 per cent of the variance. Variance percentages are the 2 eigenvalues of eq. (1), and the 2 eigenvectors yield the orientation of the 2 principal components.</figcaption></figure>

et al. 1982). PC components are classified as degeneracy when variance of the N PCs are not distinguishable within their uncertainties, in such case PCs are not reliable for pattern interpretation.

$$ \Delta\lambda_k^2 = \lambda_k^2 \sqrt{\frac{2}{n^*}} \qquad (2) $$

with $n^*$ the number of independent observations, i.e the number of samples.

Each PC component is a linear combination of all original reduced geophysical models from the projection of the later onto the eigenvectors. As we aim to perform PCA to interpret geophysical models with heterogeneous physical units, the final PC component cannot be expressed with any physical units. Thus to interpret the PCs, it is necessary to know the contribution of each original geophysical model to each PC components. The contribution of each of the N reduced geophysical models to the N PCs can be analysed by computing the matrix C where $c_{ij}$ represents the correlation between the ith geophysical model and the jth PC component (Albarede & Pyle 1996).

$$ C = U \Lambda^{1/2}. \qquad (3) $$

Hence a high correlation between one of the geophysical models and one of the PC component means a high contribution of the model to the loading. When each PC component is well resolved according to eq. (2), following Hannachi et al. (2007), we propose to interpret each PCs component as the signature of a single process affecting all data. The PC component models are, therefore, characteristics patterns of the process. This process may be associated to geological features in the area under study, depending on their different multigeophysical contributions.

# 3 GEOPHYSICAL JOINT INVERSION

## 3.1 Geophysical survey

A dense temporary seismic network of 299 vertical component geophones was deployed in 2017 from November 6th to December 7th (32 days). The geophones were set on a non-regular 30 km × 30

km grid (Fig. 2). The seismic sensors used for the survey are Fairfield Nodal Zland nodes. Processing of the raw data sets was performed by SISPROBE company. We obtained robust dispersion curves with period ranging from 0.9 to 4.7 s using cross-correlation of daily seismic recording (Shapiro & Campillo 2004) and frequency-time analysis (Levshin et al. 1989). Group velocities were inverted into Rayleigh wave group velocity maps in the period range 0.9–4.7 s using Mordret et al. (2013) approach (Fig. 3c). The maps were finally resampled on a regular grid every 700 m in both horizontal directions. The final data set used for the joint inversion consisted of 921 local dispersion curves covering the area limited by the seismic array.

We carried out 45 magnetotelluric (MT) soundings in September 2015. The network was completed with 22 more MT sites in September 2016 and 2017. In this study, we selected 48 MT sites among the 67 soundings, distributed over an area of 23 km × 23 km (Fig. 2). We performed the MT data acquisition with Metronix equipment and pairs of unpolarized Pb-PbCl2 electrodes 100 m apart in average. The magnetic fields were recorded with MF-07 and MF-06 coils. Most of the MT tensors display large diagonal components over 1 s and some tensors are characterized by off-quadrant off-diagonal phases at long periods (Ars et al. 2019, Fig. 3a).

Finally, we carried out 90 gravity measurements using a SCINTREX CG-5 gravimeter (Ars et al. 2019) which completed 537 gravity points from the International Gravity Bureau (IGB) data set. We calculated the complete Bouguer anomaly using the DEM from ASTER GDEM 30 m (NASA, 2009) and a reference density of 2.67 g cc⁻¹. We also calculated a regional Bouguer anomaly with a large data set of values from IGB covering an area four times wider than survey area. The regional anomaly represents the large scale effect of the Massif-Central uplift related to its recent volcanic activity (Granet et al. 1995) and is fairly constant across the study area. Thus, we removed a mean value from the study area gravity map to obtain the residual gravity anomaly (Fig. 3b).

## 3.2 Geophysical inversion

### 3.2.1 Magnetotelluric inversion

Following the workflow and parameters described in (Ars et al. 2019), we performed a new 3-D MT inversion with the 48 sites selected for this study with the inversion software MINIM3-D (Hautot et al. 2000, 2007). We inverted the 48 full impedance tensors for all periods available. The final RMS map is presented in Fig. 4(a). We obtained a resistivity model made of 28 × 26 × 15 cells (north 130° × east 40° × depth) covering an area of 39.1 km × 41.4 km down to 51.2 km depth. A rotation of 40° of the grid proved to be optimum to model accurately the complex phases behaviour. The parametrization led to 3268 unknowns for 11 866 data (real and imaginary part of the four terms of the impedance tensors at every period available). The resistivity model down to the maximum depth resolved by the gravity and ANT data (4 km) is displayed on Fig. 5.

### 3.2.2 Joint inversion

We run the joint inversion of the gravity and ANT data with the joint inversion algorithm developed by Ars et al. (2019). The inversion is constrained by the resistivity model set fixed. The joint inversion is an iterative process that minimizes an objective function using a non-linear steepest gradient method (Beiner 1970; Press et al. 2007). The objective function is composed of two misfit functions,
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left-hand panel: local geology. Location of the survey area (French Massif-central) is shown on the France map on the lower left-hand corner of the figure. Right-hand panel: geophysical survey of the study. The red dashed rectangle represents the investigated area by Ars *et al.* (2019).</figcaption></figure>

one for each geophysical data (gravity and seismic), that measure
the discrepancy between the model responses and the data. Model
responses are computed using the analytical solution for the gravity
effect of a set of rectangular prisms found in Blakely (1996), and
the dispersion curve computation from Wathelet (2005). It also
contains one regularization term for each geophysical model and
three coupling terms.

The coupling terms are based on the linear correlation of the
three model parameters: resistivity, density and *V*<sub>s</sub> velocity (Ars
*et al.* 2019), computed from the covariance of each pair of models
divided by the product of their standard deviations. High correlation
values between the best models retrieved from the joint inversion
outlines the common features between two parameters (resistiv-
ity/density, resistivity/velocity and velocity/density). The signifi-
cance of the correlation can be measured by the *p*-value. This is
a test of the null hypothesis that the distributions underlying the
samples are uncorrelated and are normally distributed (Kowalski
1972).

Gravity data were inverted for a 3-D density model while the local dispersion curves were inverted for a set of 1-D shear wave velocity models. All 1-D velocity profiles are stitched into a pseudo 3-D velocity model to compute the coupling terms with density and resistivity. To achieve a spatial cells correspondence between models, the resistivity, density and *V*<sub>s</sub> velocity models have the same orientation (N130°) and horizontal extension. The density model is made of 28 × 26 × 15 cells (north 130° × east 40° × depth). The lateral mesh of the density model is identical to the resistivity grid covering

an area of 39.1 × 41.4 km down to 14.7 km. The maximum depth of
the density mesh is chosen to be equal to half of the maximum grav-
ity station separation. The velocity model is made of 38 × 36 × 8
cells. The velocity model covers an area of 39.1 × 41.4 km down to
4.5 km. The depth of the seismic mesh was determined by a stochas-
tic inversion of the mean dispersion curve. As a non-unique solution
problem, ANT is dependent on the parametrization and specifically
the vertical mesh (Renalier *et al.* 2010). We choose to invert the
mean dispersion curve using a global inversion scheme to obtain
the velocities and thicknesses of three layers. We thus obtain an
optimum depth of 4.5 km. Resistivity-velocity and density-velocity
coupling terms are computed in accordance with the seismic mesh,
while the resistivity-density coupling term is computed in accor-
dance with the gravity mesh. Considering that the seismic model
is the shallowest model, the volume constrained by the three cou-
pling terms is restricted to the volume occupied by the velocity
model.

We inverted 627 gravity data and 579 local dispersion curves at all periods available ranging from 0.9 to 4.7 s. Considering the N130° oriented grid and the fact that the horizontal mesh of the seismic model is different from the density and resistivity mesh due to 1-D inverse problem ANT nature, the number of dispersion curves used for the inversion is lower than the total number of available dispersion curves. Thus we run the joint inversion for 10920 density parameters and 5220 velocity parameters for a total of 22050 data (all gravity measurements and group velocities at every period available). The standard deviation used on the data are
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Example of the impedance tensor for a MT site. Apparent resistivity of the four components of the tensor are shown on the upper panel while the lower panel displays the phase. Data are represented with dots. We superimposed the response of the best fitting model with continuous lines. (b) Map of the residual Bouguer anomaly. (c) Group velocity dispersion maps derived from a dense temporary network. Group velocity maps are displayed only for the periods indicated.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Map of the MT normalized error. (b) Map of the gravity error. (c) Map of the seismic error for each dispersion curve.</figcaption></figure>

set, respectively, to 1 mGal and 100 m s⁻¹ on the gravity and the seismic data. The initial density contrast model is homogeneous and equal to 0.0 kg m⁻³ while the initial *Vs* velocity model is a vertical gradient increasing from 3000 m s⁻¹ at the surface to 3600.0 m s⁻¹

at a depth of 4.5 km. The final gravity RMS is 0.82 mGal and the final seismic RMS is 53.38 m s⁻¹. Seismic and gravity models responses RMS are presented in Fig. 4. The gravity fit between the model responses and the data is generally good excepted near
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Layers of the resistivity model. (b) Layers of the density model. (c) Layers of the *Vs* velocity model. Blacked dash line indicates the outline of the granitic body on the surface. Red dashed line indicates location of major fault zones on the surface.</figcaption></figure>

topographic highs. The dispersion curves responses fit well the data, with every local RMS lower than 100 m s⁻¹. The data in the western part of the model are very well fitted (RMS < 60 m s⁻¹). The dispersion curves located in the centre of the survey area displays RMS > 80 m s⁻¹ due to poor model responses at long period.

## 3.3 Geophysical models

The Fig. 5 displays the resistivity model from the 3-D MT inversion (a), the density and the velocity models (b, c) from the 3-D joint inversion corresponding to the volume occupied by the velocity model, down to 4.5 km. The models are shown for four different layers at the same depth, ~0.5, ~1.5, ~2.5 and ~4 km.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Eigenvalues spectrum of the PCA, uncertainties on eigenvalues are represented with orange margins. b: Geophysical contribution from each geophysical model to each PCs model quantified by the correlation value from eq. (3).</figcaption></figure>

The final correlation between density and resistivity (DR), and Vs velocity and resistivity (VR) are respectively 0.52 and 0.42 while the final correlation between density and velocity (DV) is 0.73. The imbalance between the 3 correlation values (DR and VR correlation is lower than DV correlation) is due to the structure of the joint inversion. On the one hand since the resistivity remains unchanged during the joint inversion procedure, it tends to penalize the optimization of the DR and VR correlations as only the density and the velocity models can vary to correlate as much as possible to the resistivity model. On the other hand, the DV correlation reaches a higher value than the two other correlation DR and VR because both density and velocity models can vary to increase the correlation between them. The statistically significant correlation values, although small for DR and VR, suggest that some regions of the models are better correlated than others. The constraint of a statistical linear coupling between models allows sufficient freedom within models so that some regions of the models may not be (or weakly) linearly correlated while others may be strongly correlated. Since the linear relationship is not enforced locally, it does not preclude local non-linear relationships between the models so that they may reproduce the data as well as possible (Ars *et al.* 2019).

Thus, the model-scaled correlation between density and velocity locally increases in accordance with local geological unit such as low velocity structure correlated to a low density structure on the west in the shallow layers (Figs 5b and c) near the location of the granitic unit (dashed black circle in Fig. 2). Also, small low scattered velocity structures seem to be correlated with small conductive bodies below 2 km depth near the location of major fault zone (dashed red area in Figs 5a and c). Local correlation can be highlighted by interpreting models in parallel, but it remains very hard to properly distinguish geological units, faults, fluid saturation or alteration within the geophysical models. Indeed geological features may display various geophysical signatures depending on imaging methods. In order to interpret geological structures from the three geophysical models, we decide to use PCA to gain statistical insight into the multivariable system made of the resistivity, density and velocity models.

# 4 MODELS INTEGRATION

## 4.1 PCA application

This section presents the PCA statistical integration of the jointly inverted models in order to clarify the distribution of the parameters correlations. For this purpose, we tested the PCA with the log-resistivity, density contrast and Vs velocity parameters within

the common part of the models obtained from the joint inversion constrained by the resistivity distribution. We first centre the distribution of each parameter set by removing their mean values. Each set is then normalized with its own standard deviation. The eigenvectors and eigenvalues are then computed from the covariance matrix S (eq. 1). The three components PC1, PC2, PC3 are sorted from 1 to 3 according to their decreasing eigenvalues. They explain, respectively 69, 22 and 9 per cent of the total variance. Their uncertainties are, respectively 7.5, 2.4 and 1 per cent of the total variance (Fig. 6a). We observed that none of the three eigenvalues including their error-bars are overlapping each other (Fig. 6a). This indicates than all three loading are robust components (Hannachi et al. 2007). We then computed the correlation matrix (eq. 3) in order to quantify the contribution of each geophysical parameter to each PC (Fig. 6b). PC1 shows an almost identical negative correlation with the three parameter sets (see Fig. 6b). This component is strongly influenced by the correlation between all model parameters enforced by the coupling terms of the joint inversion. PC2 is mostly driven by an anticorrelation with the resistivity set (see Fig. 6b). While PC3 is moderately correlated to the velocity set and moderately anticorrelated to the density set and with no contribution from the resistivity (see Fig. 6b).

Thus, we plotted the three normalized PC component models (see Fig. 6) by attributing to each cell of the mesh common to the velocity, density, resistivity models the loading values for each PCs computed as the projection of the reduced models onto the eigenvectors. Since the PC component are a linear combination of the normalized geophysical parameters (log-resistivity, density and velocity), it is dimensionless and only their variations across the mesh can be interpreted. PCs models are shown for four different layers of approximately the same depth for each model: ~0.1, ~1, ~2 and ~4 km (see Fig. 7).

## 4.2 PCA highlighted geophysical integration

PC models may be interpreted in terms of an underlying physical process that affects differently the resistivity, density and velocity and their distribution can be interpreted in terms of the geological structures where this physical process dominates. Such geophysical contributions to each PCs models are displayed in Fig. 6 as correlation value. The following discussion is based on comparison between PCs model patterns, prior fault and geological models (see Fig. 7) considering their geophysical contribution.

PC1 model display a positive loading above 250 m depth and a negative loading below 3.0 km depth. Considering the global anticorrelation of the 3 parameters to PC1, such loading may be
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Layers of the normalized PC1 model overlaid with location of the local faults on the surface. (b) Layers of the normalized PC2 model overlaid with location of the local faults on the surface. c: Layers of the normalized PC3 model overlaid with the prior geological model. Granitic units are represented with tilted cross hatches, metamorphic units with vertical cross hatches and volcanic units with dot hatches. Green circles indicate A, B, C and D anomalies and green dashed lines indicate SW-NE corridor discussed in Section 4.2.</figcaption></figure>

associated to very altered rock with low velocity, low density and high conductivity near the surface and undamaged rock with strong velocity, strong density and strong resistivity at depth. In between, PC1 layers display contrasted loading with a trend towards negative loading in a SW-NE corridor (green dashed line on Fig. 7a) bounded by faults and positive loading on the north of the survey. This may also indicate a corridor of unweathered rock material surrounded by a highly fractured environment. At the surface, the negative small feature A (Fig. 7a) coincides with a correlated high resistivity, high

velocity and high density anomaly located at fault intersections. The presence of volcanic features at the surface (dot hatches on Fig. 7c) and the contrasting nature of the A anomaly may lead to associate it to intact volcanic intrusion within strongly altered layers. In contrast, the deep positive features (B, C and D in Fig. 7) coincide with low resistivity, low density and low velocity anomalies correlated in a supposedly unweathered medium and close to the fault intersection. These features contrast sharply with the surrounding bedrock and may be associated with a deep-seated weathered medium located in
a permeable fault zone. Such interpretations based on PC1 loading behaviour lead to understand PC1 as an indicator of rock integrity allowing to identify altered material that could be related to geothermal activity.

As PC2 is mostly anticorrelated to resistivity model (Fig. 6), PC2 positive features may be seen as conductive features. PC2 positive structures extend continuously from surface to depth and superimpose with the intersection of the main local faults (Duwiquet et al. 2019) (see Fig. 7b). The superposition with the prior fault model and the anticorrelation to resistivity indicate that positive PC2 pattern may be interpreted an indicator of fluid-saturated fracture zones and highlighting the presence of fluids of interest for geothermal exploration. Moreover the B, C and D anomalies display strong positive PC2 loading (Fig. 7b) which could be associated to saturated fault zones, coherent with deep-rooted altered zones described from PC1. Contrary, the previously described unaltered SW-NE corridor(green dashed line on Fig. 7b) displays negative loading which is consistent as an unweather dry material. Furthermore the slight correlation of PC2 to density and velocity may lead to interpret positive PC2 anomalies as indication of local process that slightly increase velocity and density in faulting zone, like cementation (Tripetta et al. 2017) due to geothermal fluid circulation in permeable fault.

The geological model presented in Fig. 2 (Bellanger et al. 2019) superimposed on the PC3 model (see Fig. 7c) display a good agreement between the lithology and PC3 patterns down to 1 km depth. Granitic body (tilted mesh hatches on Fig. 7c) are correlated to positive PC3 feature while metamorphic rock (vertical mesh hatches on Fig. 7c) correspond to negative feature. Such differentiation may be consistent with the velocity correlation and density anticorrelation to PC3. Indeed, granite tends to be lighter and faster than metamorphic units, due to their stratified nature. (Abers & Hacker 2016). Beyond 1 km depth, both geological and PC3 component models display the same trend to extend the granitic body eastward but at different depth. Superposition of the geological model and the PC3 model is lost at such a depth, mainly because of the lack of constrains with depth on the geological model. However SW-NE corridor mentioned in PC1 and PC2 seems to lay in depth within a null-positive PC3 (green dashed line on Fig. 7c) that has been associated in the shallow layer to granite lithology. Such rock material would be coherent with the unaltered, unsaturated nature described with PC1 and PC2 loading. Finally, positive PC3 loading associated to granite seems to be encircled by main PC2 positive features. Such symbiosis between PC2 and PC3 loading may indicate that PC3 seems to be a component controlled by the density and the velocity and its loading may vary in accordance with the lithology while PC2 highlight saturated fractures that lay at the geological interfaces which may act as natural drain for geothermal circulation.

## 5 CONCLUSION

We, in this paper, the application of a PCA methodology for geological interpretation in a context of multigeophysical imaging. PCA allows to retrieve a set of *N* PC models by solving an eigenvalue problem for a matrix composed of *N* geophysical models. PC components are linear combinations of the original geophysical models. We show that geophysical signatures within each PCs models allow site-dependent geological feature interpretation. As an application, we used models derived from the joint inversion of MT, gravity and ambient noise tomography data from the French Massif-Central

applied to geophysical exploration for unconventional geothermal system. PC1 model is mostly influenced by the coupling enforced during the joint inversion and seems to act as a rock material integrity indicator and allows the identification of altered regions that may be linked to geothermal activity. PC2, which is mostly anticorrelated with resistivity, plays the role of fluid-saturated zone highlighter, distinguishing faults of interest for geothermal exploration. Finally, PC3 provides information about the lithology and allows the reconstruction of a conceptual model of 3-D geological units down to 4.5 km. PC component models also suggest the presence of faults connected to deep conductor structures, that may be related to geothermal circulation and/or alteration. PCA demonstrates its ability to bring out complex interactions between geophysical models allowing their interpretation in terms of geological units, presence of fluid or integrity of rock material within the scope of geothermal exploration. In general, multivariate PCA seems to be efficient to help identify geological units and/or geological processes by highlighting the contribution of each geophysical parameter.

## ACKNOWLEDGMENTS

We thank KIT Karlsruhe and INSU, EOPGS Strasbourg to let us use the magnetotelluric equipment. We would like to acknowledge the International Gravity Bureau for providing us gravity data. We also thank all the local land owners for their permission to install the equipment, Christophe Voisin for the proofreading and advice to make this paper more understandable, and Luis Gallardo and an anonymous reviewer for their fruitful comments which enabled us to greatly improve the manuscript. The financial support was provided by TLS-GEOTHERMICS SAS and Project ANR-19-CE05-0043 GERESFAULT.

## DATA AVAILABILITY

Data are the property of the private company TLS-Geothermics. They are not freely available.

## REFERENCES

* Abers, G.A. & Hacker, B.R., 2016. A matlab toolbox and e xcel workbook for calculating the densities, seismic wave speeds, and major element composition of minerals and rocks at pressure and temperature, *Geochem. Geophys. Geosyst.*, **17**(2), 616–624.

* Albarede, F. & Pyle, D., 1996. Introduction to geochemical modeling., *Geol. Mag.*, **133**(1), 121–121.

* Ars, J.-M., Tarits, P., Hautot, S., Bellanger, M., Coutant, O. & Maia, M., 2019. Joint inversion of gravity and surface wave data constrained by magnetotelluric: Application to deep geothermal exploration of crustal fault zone in felsic basement, *Geothermics*, **80**, 56–68.

* Barrash, W. & Morin, R.H., 1997. Recognition of units in coarse, unconsolidated braided-stream deposits from geophysical log data with principal components analysis, *Geology*, **25**(8), 687–690.

* Bauer, K., Muñoz, G. & Moeck, I., 2012. Pattern recognition and lithological interpretation of collocated seismic and magnetotelluric models using self-organizing maps, *Geophys. J. Int.*, **189**(2), 984–998.

* Beiner, J., 1970. *Fortran Routine Mindef for Function Minimization*, pp. 1–13, Institut de Physique, Univ. of Neuchatel.

* Bellanger, M., Hermant, B., Galibert, S. & Auxiètre, J., 2019. Fault-controlled hydrothermal system associated with major crustal fault zone: future drilling target to assess the deep geothermal potential—the Sioule license project, in *Proceedings of the European Geothermal Congress 2019*, Den Haag, The Netherlands, 11–14 June 2019.
Blakely, R.J., 1996. *Potential Theory in Gravity and Magnetic Applications*. Cambridge Univ. Press.

Coppo, N., Baltassat, J., Bouchot, V., Gadalia, A., Girard, J., Hautot, S. & Wawrzyniak, P., 2015. 3-D Magnetotelluric Investigations in Martinique (Lesser Antilles). Characteristic Deep Resistivity Structures and Shallow Resistivity Distribution Matching Heliborne TEM results, in *Proceedings World Geothermal Congress*.

Coutant, O., Bernard, M., Beauducel, F., Nicollin, F., Bouin, M. & Roussel, S., 2012. Joint inversion of p-wave velocity and density, application to la soufrière of guadeloupe hydrothermal system, *Geophys. J. Int.*, **191**(2), 723-742.

Craig, J.W., Faulds, J.E., Hinz, N., Earney, T.E., Schermerhorn, W.D., Siler, D. & Deoreo, S., 2021. Discovery and analysis of a blind geothermal system in southeastern Gabbs Valley, western Nevada, USA. *Geothermics*, **97**, 102-177.

Duwiquet, H., Arbaret, L., Guillou-Frottier, L., Heap, M.J. & Bellanger, M., 2019. On the geothermal potential of crustal fault zones: a case study from the pontgibaud area (French Massif-Central, France), *Geothermal Energy*, **7**(1), 33, doi:10.1186/s40517-019-0150-7.

Gallardo, L.A. & Meju, M.A., 2007. Joint two-dimensional cross-gradient imaging of magnetotelluric and seismic traveltime data for structural and lithological classification, *Geophys. J. Int.*, **169**(3), 1261-1272.

Gola, G. et al., 2017. Data integration and conceptual modelling of the Larderello Geothermal Area, Italy, *Ener. Proc.*, **125**, 300-309.

Granet, M., Wilson, M. & Achauer, U., 1995. Imaging a mantle plume beneath the French Massif Central, *Earth planet. Sci. Lett.*, **136**(3-4), 281-296.

Hannachi, A., Jolliffe, I. & Stephenson, D., 2007. Empirical orthogonal functions and related techniques in atmospheric science: a review, *Int. Climatol.*, **27**(9), 1119-1152.

Hautot, S., Tarits, P., Whaler, K., Le Gall, B., Tiercelin, J.-J. & Le Turdu, C., 2000. Deep structure of the Baringo Rift Basin (central Kenya) from three-dimensional magnetotelluric imaging: implications for rift evolution, *J. geophys. Res.*, **105**(B10), 23 493-23 518.

Hautot, S., Single, R., Watson, J., Harrop, N., Jerram, D., Tarits, P., Whaler, K. & Dawes, D., 2007. 3-D magnetotelluric inversion and model validation with gravity data for the investigation of flood basalts and associated volcanic rifted margins, *Geophys. J. Int.*, **170**(3), 1418-1430.

Jegen, M.D., Hobbs, R.W., Tarits, P. & Chave, A., 2009. Joint inversion of marine magnetotelluric and gravity data incorporating seismic constraints: preliminary results of sub-basalt imaging off the Faroe Shelf, *Earth planet. Sci. Lett.*, **282**(1), 47-55.

Kassenaar, J., 1991. An application of principal components analysis to borehole geophysical data, in *Proceedings of the Fourth International Symposium on Borehole Geophysics for Minerals, Geotechnical and Groundwater Applications*, pp. 211-218, Toronto, Ontario.

Kowalski, C.J., 1972. On the effects of non-normality on the distribution of the sample product-moment correlation coefficient, *J. R. Stat. Soc.*, **21**(1), 1-12.

Kvamme, K.L., 2006. Integrating multidimensional geophysical data, *Archaeol. Prospect.*, **13**(1), 57-72.

Lehujeur, M., Vergne, J., Schmittbuhl, J., Zigone, D., Le Chenadec, A. & Team, E., 2018. Reservoir imaging using ambient noise correlation from a dense seismic network, *J. geophys. Res.*, **123**(8), 6671-6686.

Levshin, A.L., Its, E., Lander, A., Bukchin, B., Barmin, M. & Ratnikova, L., 1989. *Seismic Surface Waves in a Laterally Inhomogeneous Earth*, Vol. 9, Springer Science & Business Media.

Meju, M.A. & Gallardo, L.A., 2016. Structural coupling approaches in integrated geophysical imaging, in *Integrated Imaging of the Earth: Theory and Applications*, pp. 49–67, eds Moorkamp, M., Lelièvre, P. G., Linde, N. & Khan, A., Wiley.

Mordret, A., Landès, M., Shapiro, N., Singh, S., Roux, P. & Barkved, O., 2013. Near-surface study at the Valhall Oil Field from ambient noise surface wave tomography, *Geophys. J. Int.*, **193**(3), 1627–1643.

Muñoz, G., Bauer, K., Moeck, I., Schulze, A. & Ritter, O., 2010. Exploring the Groß SchöNebeck (Germany) geothermal site using a statistical joint interpretation of magnetotelluric and seismic tomography models, *Geothermics*, **39**(1), 35–45.

North, G.R., Bell, T.L., Cahalan, R.F. & Moeng, F.J., 1982. Sampling errors in the estimation of empirical orthogonal functions, *Mon. Weather Rev.*, **110**(7), 699–706.

Paasche, H., Tronicke, J., Holliger, K., Green, A.G. & Maurer, H., 2006. Integration of diverse physical-property models: subsurface zonation and petrophysical parameter estimation based on fuzzy c-means cluster analyses, *Geophysics*, **71**(3), H33–H44.

Power, M., Belcourt, G. & Rockel, E., 2004. Geophysical methods for kimberlite exploration in northern Canada, *Leading Edge*, **23**(11), 1124–1129.

Press, W., Teukolsky, S., Vetterling, W. & Flannery, B., 2007. *Numerical Recipies*, 3rd edn, Cambridge Univ. Press.

Ranjbar, H., Hassanzadeh, H., Torabi, M. & Ilaghi, O., 2001. Integration and analysis of airborne geophysical data of the Darrehzar Area, Kerman Province, Iran, using principal component analysis, *J. Appl. Geophys.*, **48**(1), 33–41.

Renalier, F., Jongmans, D., Savvaidis, A., Wathelet, M., Endrun, B. & Cornou, C., 2010. Influence of parameterization on inversion of surface wave dispersion curves and definition of an inversion strategy for sites with a strong vs contrast, *Geophysics*, **75**(6), B197-B209.

Roche, V., Sternai, P., Guillou-Frottier, L., Jolivet, L. & Gerya, T., 2017. Eastern Mediterranean geothermal resources and subduction dynamics, in *EGU General Assembly Conference Abstracts*, Vol. 19, p. 13427.

Schill, E., Geiermann, J. & Kümmritz, J., 2010. 2-D magnetotellurics and gravity at the geothermal site at soultz-sous-forêts, in *Proceedings World Geothermal Congress*, pp. 25-29.

Shapiro, N.M. & Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, *Geophys. Res. Lett.*, **31**(7), 614, doi:10.1029/2004GL019491.

Trippetta, F., Carpenter, B., Mollo, S., Scuderi, M., Scarlato, P. & Collettini, C., 2017. Physical and transport property variations within carbonate-bearing fault zones: Insights from the Monte Maggio Fault (Central Italy), *Geochem. Geophys. Geosyst.*, **18**(11), 4027–4042.

Wannamaker, P.E., Doerner, W.M. & Hasterok, D.P., 2007. Integrated dense array and transect Mt Surveying at Dixie Valley Geothermal Area, Nevada; structural controls, hydrothermal alteration and deep fluid sources, in *Proceedings of the 32th Workshop on Geothermal Reservoir Engineering*, Stanford University.

Wathelet, M., 2005. Array recordings of ambient vibrations: surface-wave inversion, *PhD dissertation*, Faculté des Sciences Appliquées, Liège University, 161pp.

Widaman, K.F., 1993. Common factor analysis versus principal component analysis: differential bias in representing model parameters?, *Multivar. Behav. Res.*, **28**(3), 263-311.