Nolwenn, Mérédith, Franck, Daniel Fred Catalina

Tobi Sara

Camille Jestin

M Aquino

G Marquis

J Vergne

surface-wave dispersion curves using correspondence maps

Keywords: Joint inversion, correspondence maps, magnetotelluric, ambient seismic noise, surface-wave dispersion curves Passive methods, magnetotelluric, ambient seismic noise, joint inversion, correspondence maps

In this Thesis, I present a joint inversion framework using correspondence maps (CM) between surface-wave dispersion curves recovered from ambient seismic noise correlations and magnetotelluric data to solve for shear velocity and resistivity with depth, but at the same time for a possible relationship between them. The method is stable using synthetic and field data for 1D models and shows a reduction of the model space that fit all datasets simultaneously.

A sensitivity study of surface-wave dispersion curves using 2D velocity models exposes some limitations on their lateral resolution for small-sized anomalies (<10 km) at the periods of interest (<10 s), and where MT data could be beneficial. Synthetic joint inversion tests in 2D using a first-order relationship show that the velocity model is substantially improved by the inclusion of magnetotelluric data and the resistivity-velocity relationship is better recovered than with separate inversions.

The application of the joint inversion methodology to field data from North Alsace exposes some of its advantages but also some of its difficulties. For an area that shows a 1D model behavior, the inclusion of the correspondence maps produces more coherent models, in both petrophysical and structural aspects. For the case of a more complex subsurface, where several relationships could be present, the CM term is not able to produce a single relationship applicable to all structurally different zones.

Résumé

Dans cette thèse, je présente un cadre d'inversion conjointe entre des courbes de dispersion des ondes de surface récupérées à partir de corrélations de bruit sismique ambiant et des données magnétotelluriques pour imager la vitesse de cisaillement et la résistivité avec la profondeur, en y ajoutant des cartes de correspondance (correspondence maps, CM) pour trouver une éventuelle relation entre ces propriétés physiques. Testée d'abord sur des modéles synthétiques 1D, la méthode est stable et amène bien une réduction de l'espace des modèles qui satisfont tous les jeux des données simultanément.

Une étude de sensibilité des ondes de surface à l'aide de modèles de vitesse 2D montre les limites de leur résolution latérale pour de petites anomalies (<10 km) aux périodes

Research context

Following Tarantola [2005], one of the leaders of inverse methods in geophysics, a physical system can be defined by three elements:

1. Parametrization of the system: set of model parameters whose values fully characterize the system.

2. Forward modeling: definition of physical laws that allow to make predictions i.e. observable data from given values of model parameters.

3. Inverse modeling: consists of using the actual measurements (observable data) to infer the values of the parameters that characterize the system i.e. model parameters.

This Thesis focuses on the third item i.e. inverse modeling. In geophysics, inversion is commonly referred as the problem of using observable data to estimate earth model parameters (seismic velocities, density, resistivity etc.). For example, we use seismic travel times, to estimate seismic velocities with depth. Inverse problems usually regard one type of model parameter. For the previous example, seismic velocity is the model parameter to estimate. In the case of gravity, we invert measurements to estimate the earth's density.

The attempt to combine multiple datasets to solve for several model parameters at the same time is called joint inversion, which is the focus of this Thesis. I use the term joint inversion in this Thesis as defined by [START_REF] Moorkamp | Introduction, chapter 1[END_REF]: "the approach where Chapter 1. Introduction different types of observable data are inverted with a single objective function and where all model parameters are adjusted concurrently throughout the inversion". Contrary to cooperative inversion where one single data is inverted and the results of another data type inversion are used for reference.

Joint inversion of geophysical data has been proposed as means for understanding subsurface structures, since taking advantage of different data sensitivities can reduce the range of acceptable model parameters (i.e: the 'non-uniqueness') [Moorkamp et al., 2007]. One can implement joint inversion in many different ways: 1) joint inversion for a single model parameter from multiple data sets sensitive to it, also called mono-model joint inversion, e.g apparent resistivity from DC resistivity and magnetotelluric data [Jupp and Vozoff, 1975] or shear velocity from receiver functions with surface-wave dispersion curves [Julià et al., 2000], and 2) joint inversion of data sensitive to different physical parameters, multi-modal joint inversion , e.g. resistivity and seismic velocity [Hu et al., 2009] [Gallardo andMeju, 2003]. For the first type, the joint inversion procedure is simpler as both methods share a common subsurface model. The second type however is more challenging since sensing different physical properties involves questions about their mutual coupling (or lack thereof). In this Thesis I focus on the second type of joint inversion (i.e. multi-modal joint inversion) to obtain shear-wave velocity and electrical resistivity models, from surface-wave dispersion curves and magnetotelluric data respectively. The core element of any joint inversion methodology is therefore the hypothesis of a link between the different physical property models [Carrillo and Gallardo, 2018]. In a general way, these links or features can be classified as structural or petrophysical.

Structurally-coupled joint inversion assumes that different geophysical methods sense the sale underlying geology, therefore it is the structure that controls the distribution of petrophysical properties [Gallardo and Meju, 2003]. With this approach, the structural similarity of the different physical property models are considered rather than petrophysical relationships between these properties. Today, the most common structural coupling method is the cross-gradient constraint, introduced by Gallardo and Meju [2003] and Gallardo and Meju [2004] and for which structural similarity is enforced by minimizing the cross-gradient vector. Other structural constraints include curvature based methods e.g. [Zhang and Morgan, 1997], [Haber and Oldenburg, 1997] or Gramian constraints e.g. [Zhdanov et al., 2012], [START_REF] Ogunbo | Gramian constraints in electromagnetic multi-physics joint inversion[END_REF].

Structural coupling does not make any assumption about the petrophysical relationships between the parameters. If there is a priori information about such relationships, joint inversion results can be improved significantly. Hence, structural coupling might overlook the full potential of the joint inversion [Moorkamp, 2017]. Joint inversion methods using parameter relationships are based on the fact that for some specific geological environments, geophysical parameters can be related by some mathematical relationship.

For example, Carcione et al. [2007] showed how in sedimentary environments, porosity and fluids of the rock matrix can provide a physical link that causes a correlation between seismic velocities and electrical conductivities.

So-called petrophysical joint inversion, has generally been approached in two ways:

1.1. Research context 1) when there is a relationship between the parameters which is used to constrain the depending parameters, in a sequential inversion form e.g [START_REF] Panzner | Integrated seismic and electromagnetic model building applied to improve subbasalt depth imaging in the faroe-shetland basin[END_REF] and2) when the parameter relationship is included as part of the joint inversion workflow. From this second group, Tiberi et al. [2008] and Dufréchou et al. [2018] retrieved linear relationships with depth between density and compressional velocity, in a Birch relationship type [START_REF] Birch | The velocity of compressional waves in rocks to 10 kilobars: 1[END_REF]. Sun and Li [2016] use the fuzzy C-means (FCM) clustering technique to handle different shapes of petrophysical relationships. Carrillo and Gallardo [2018] proposed a way to retrieve higher-order (non-linear) relationships for a number of heterogeneous zones while fitting the data using correspondence maps (CM). Generally, joint inversion approaches using petrophysical relationships had the inconvenient of being difficult to generalize since using parameter relationships suited for specific geological settings.

In geophysical exploration, the methods used can be classified in 1) active or 2) passive methods depending on the source from where the data is generated. Active methods tend to be more expensive and involve a higher environmental footprint than passive methods. Among passive methods we can find: magnetotelluric, gravity, magnetics, ambient seismic noise etc. Magnetotelluric studies have been applied in geophysics since several decades, while Ambient Seismic Noise tomography has been expanding rapidly since early 2000s, [Shapiro and Campillo, 2004], [Shapiro et al., 2005] , [START_REF] Sabra | Extracting time-domain green's function estimates from ambient seismic noise[END_REF], [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF] as seen in Fig. 1.1. However, as seen in Fig. 1.2, these two methods represent the lowest contribution of the studies using joint inversion approaches.

We focus in the least studied geophysical methods in joint inversion, i.e. Magnetotelluric data and Ambient Seismic Noise.

Throughout this Thesis I study a joint inversion methodology that includes the CM approach described by Carrillo and Gallardo [2018], implemented as an Occam-type inversion [Constable et al., 1987] to solve for shear velocity and apparent resistivity.

Research questions

The main questions I aim to answer in this Thesis are:

1. Can we implement a joint inversion approach using ambient seismic noise and magnetotelluric data in one and higher-dimension models?

2. Can we obtain meaningful velocity-resistivity relationships from geophysical data? 3. How similar or different models from structural and parameter-relationship joint inversion approaches can be?

Manuscript organization

This Thesis is based on our published work and works in progress. Fig. 1.3 exemplifies the organization of this manuscript.

In Chapter 2, I present the two passive geophysical methods used throughout this Thesis i.e. magnetotelluric and surface-wave dispersion curves from ambient seismic noise correlations.

In Chapter 3, I review the state of art of several joint inversion approaches and introduce the correspondence maps used in this Thesis.

In Chapter 4, based on [START_REF] Aquino | Joint one-dimensional inversion of magnetotelluric data and surface-wave dispersion curves using correspondence maps[END_REF], I use a CM approach to first asses onedimensional problems. I calculate the retrieval of polynomial coefficients given model parameters known and after jointly invert surface-wave dispersion curves and magnetotelluric data for both synthetic and field data. Some work done that fall outside of the main focus of the Thesis i.e. joint inversion is described in the appendix sections:

Appendix A shows the processing carried out for magnetotelluric data used for the field case study in Chapter 6.

In Appendix B, I compare the results of two inversion methods of surface-wave dispersion curves using the EstOF data set in North Alsace, France 1) Least-Squares approach i.e. Occam inversion [Constable et al., 1987] and 2) Monte-Carlo inversion i.e. Metropolis sampling. These two results are compared with a reference model computed by Lehujeur et al. [2018].

Appendix C shows some preliminary work using L-BFGS-B optimization method in 1D.

Contributions

The research presented in this manuscript has been the object of one published article, three oral presentations in international conferences, and works to be submitted for pub- My research makes use of passive geophysical data, namely magnetotellurics and surface-wave dispersion curves obtained from ambient seismic noise correlations. In this Chapter, I describe these two methods and their respective data processing steps.

Chapter 2. Passive geophysical methods

Magnetotelluric method (MT)

Magnetotellurics (MT) is one of the most widely-used geophysical methods to investigate the Earth's subsurface electrical conductivity. It is based on the presence of natural magnetic fields that produce electric currents into the Earth. The interest of the MT method comes from the different natural sources that enable a characterisation of the subsurface over different depths. An example of the spectrum of electromagnetic (EM) sources surrounding the Earth extending across a wide frequency range is shown in Figure 2.1.

The signals measured in MT are the natural electric and magnetic fields in the frequency range 10 -4 to 10 4 Hz, although in general only a part of this range is used in any particular application [Vozoff, 1990]. In this frequency band, the source fields originate either in the Earth's magnetosphere where the complex interaction of the solar wind and the geomagnetic field generates hydromagnetic (plasma) waves [Vozoff, 1990] for frequencies below 1 Hz, or from atmospheric (sferics in 2.1) electrical discharges for frequencies greater than 1 Hz.

The first 100 km of the atmosphere is an electrical insulator, and so the fields generated above this level form plane electromagnetic waves that propagate toward the Earth. At the Earth's surface most of the plane-wave electromagnetic energy is reflected back into the atmosphere and a small amount penetrates into the ground. This source magnetic field induces electrical (telluric) currents in the Earth, which in turn can create secondary magnetic fields in conductors and where electrical conductivity contrasts occur, e.g. [Key, 2003]. As will be shown below, by measuring the electric and magnetic fields at the surface of the Earth, information on the electrical conductivity of the subsurface can be obtained.

Electric and magnetic fields in a homogeneous medium

In the MT method, the physical phenomenon considered is the electromagnetic induction, which consists on the creation of electric currents in a conductor given some time variation of a magnetic source field. The phenomenon can be described from Maxwell's equations, following the approach described in e.g. [Ward and Hohmann, 2012] ∇ × E = -∂B ∂t (2.1)

∇ × H = J + ∂D ∂t (2.2)
where:

2.1. Magnetotelluric method (MT)

D = dielectric displacement (C m -2 )
In an homogeneous isotropic medium the scalar constitutive relations apply

D = ϵE B = µH J = σE
where ϵ is the dielectric permittivity (F m 1 ), µ is the magnetic permeability (H m -1 ) and σ is the electrical conductivity (S m -1 ). Applying the curl to equation 2. Using the vector identity ∇ × (∇ × u) = ∇(∇ • u) -∇ 2 u in an isotropic homogeneous medium, i.e. where no electrical conductivity contrast is present, ∇ • E = 0 which yields the wave equation

∇ 2 E = µ(σ ∂E ∂t + ϵ ∂ 2 E ∂t 2 ) (2.
3)

The Fourier transformation of equation 2.3 yields the Helmholtz equation

(∇ 2 + k 2 )E = 0 (2.4)
where k is the complex wave number in the medium and ω the angular frequency.

k 2 = µω(ϵω -iσ) (2.5)
For the majority of EM methods in geophysics, with the notable exception of groundpenetretating radar, the frequencies used are relatively low, i.e. lower than tens of KHz, therefore

ω 2 µϵ << ωµσ → k 2 ≈ -iωµσ (2.6)
This assumption, known as the quasi-static approximation, means that the second term on the right-hand-side of the wave equation 2.3 can be neglected, leading to a diffusion equation:

∇ 2 E = µσ ∂E ∂t (2.7)
Using the same line of reasoning for the magnetic field intensity,

∇ 2 H = µσ ∂H ∂t (2.8)
These equations describe the behavior of the electric and magnetic fields in the Earth in the frequency band used in MT. They are similar to the heat equation and hence describe diffusive fields. The frequency-domain equivalents of equations 2.7 and 2.8 are

∂ 2 E ∂z 2 + k 2 E = 0 (2.9) ∂ 2 H ∂z 2 + k 2 H = 0 (2.10)
with solutions E = E 0 e -ikz = E 0 e -iαz e -βz (2.11)

Magnetotelluric method (MT)

H = H 0 e -ikz = H 0 e -iαz e -βz (2.12)

where

α = β = σµω 2 (2.
13) E 0 and H 0 are the values of the fields at z = 0. The term e -βz in equations 2.11 and 2.12 describes an exponential attenuation in the vertical direction. The depth at which the field amplitude is reduced by a factor 1/ϵ or about 37% is known as the skin depth

δ δ = 1 β = 2 σµω ≈ 503 1 σf meters (2.14)
It is worth noting here that δ decreases with increasing frequency and increasing conductivity. Using several frequencies allows for different δ and results in a resistivity depth sounding.

Relationship between electric and magnetic fields

Faraday's law ( equation 2.1) in the frequency domain imply

∂E z ∂y - ∂E y ∂z = -iµωH x (2.15) ∂E x ∂z - ∂E z ∂x = -iµωH y (2.16) ∂E y ∂x - ∂E x ∂y = -iµωH z (2.17)
In a uniform plane wave H propagating in the z direction in a half space both E z and H z and their derivatives vanish. Substituting equations 2.11 and 2.12 in 2.16 and 2.17 and substituting derivatives in z by a multiplication by k yield

E y H x = - µω k (2.18) E x H y = µω k (2.19)
The ratio E/H is defined as the electrical impedance Z. Two parameters are defined from this impedance i.e. resistivity ρ (Ω .m) and phase (in radians) [Cagniard, 1953]

ρ = |Z| 2 µω = 1 µω E x H y 2 (2.20) ϕ = arg E x H y = arg µω k = arg ωµ σ e iπ/4 = π 4 (2.21) Chapter 2.

Passive geophysical methods

The phase for a half-space is always equal to π 4 ± π.

Apparent resistivity

Early magnetotelluric studies were based on an interpretation of MT data with a layered medium, i.e. where resistivity only varies vertically as in Figure 2.2. The impedance at the surface of a one-dimensional model can be found by adding the effects of upward and downward traveling electromagnetic fields in each layer.

Following the notation used by Ward and Hohmann [2012], the effective impedance at the top of layer i ( Ẑi ) in a 1D model with n layers can be calculated recursively

Ẑi = Z i Ẑi+1 + Z i tanh(ik i h i ) Z i + Ẑi+1 tanh(ik i h i ) (2.22)
h i is the layer thickness and Z i is the intrinsic impedance of layer i.

Z i = ωµ k i (2.23)
The recursion starts at layer n and migrates up to the surface. It is important to mention that the surface impedance at the top of a given layer is independent of the layers above it. The effective impedance Ẑi can be transformed to apparent resistivity ρ a

ρ a = 1 µω | Ẑ1 | 2 = 1 µω E x H y 2 (2.24)
ρ in equation 2.20 is the true resistivity of the half-space while ρ a due to the diffusive process is called apparent as it contains superimposed effects from each layer. Higher frequencies being more sensitive to shallow layers and lower frequencies to deeper layers imply that ρ a varies with frequency. Apparent phase ϕ a is no longer restricted to π 4 ± π and also varies with frequency.

Magnetotelluric method (MT)

Impedance tensor

In media where the distribution of conductivity shows a higher dimensionality, i.e. 2D or 3D structures, the impedance (and hence the apparent resistivity) is dependent on the source field polarization and on the measurement directions of the electric and magnetic fields. The relationship between the horizontal fields in the frequency domain can be then expressed in a matrix equation form E = ZH.

E x E y = Z xx Z xy Z yx Z yy H x H y (2.25)
Z is now an impedance tensor, in this Thesis, Z is calculated using the Razorback Python library from [START_REF] Smaï | Razorback, an open source python library for robust processing of magnetotelluric data[END_REF] described in Section 2.2.2. The four elements of Z help determine its dimensionality and directionality, both related to those of the subsurface resistivity distribution.

For the 1D case we have

0 Z xy -Z xy 0 (2.26)
For the 2D case, if x or y is along the geoelectrical strike

0 Z xy Z yx 0 (2.27)
The off-diagonal impedances in this case form two independent modes

• The transverse electric (TE) mode, the electric field is oriented along the strike of conductivity structure (e.g. x direction) and the magnetic field is oriented across the strike.

• The transverse magnetic (TM) mode, the magnetic field is oriented along the strike and the electric field is oriented across the strike.

For the 3D case

Z xx Z xy Z yx Z yy (2.28)
The source polarization apparent resistivity (here for E x and H y ) is defined

ρ xy = 1 ωµ |Z xy | 2 (2.29)
and phase

ϕ xy = arg(Z xy ) (2.30)
Chapter 2. Passive geophysical methods

MT data workflow

MT data acquisition

Magnetometers and electrodes are disposed in the field as shown in the schema of Figure

MT data processing

The first step in MT analysis is the estimation of transfer functions (TF) (i.e. impedance tensor Z). Normally MT time series are transformed to evenly spaced frequencies (Fourier transform). For each frequency, the impedance tensor in equation 2.25 can be rewritten

as: e = bz + ϵ, with N observations (i.e. N Fourier transforms of N independent data sections at a given frequency), e and b are the horizontal electric and magnetic field components at a specific site and frequency, and ϵ is an N-vector of random errors. Thereafter the impedance z is a TF to be estimated. If O is the operator yielding those estimations, residual r can be computed as r = e -bO(e, b).

Single site least-squares estimate

Least-squares estimation seeks to find the solution that minimizes the sum of the squared data residuals. The O operator can be computed for a single site data set as

z = (b H b) -1 (b H e) (2.31)
where the super index H denotes the Hermitian (complex conjugate) transpose.

Remote-reference least-squares estimate

The remote-reference (RR) method uses measurements of an additional magnetic field b r at a remote site to reduce the effects of noise. The RR location should be located far away so that the noise between the station and the remote is incoherent, but close enough so that the magnetic source is uniform between the sites. The RR estimate is given by:

z = (b H r b) -1 (b H r e) (2.32)
MT data can still be affected by correlated noise that can be reduced with several RRs.

M-estimator

Remote-reference technique greatly improved MT processing, but it lacked the ability to handle large outliers present in both the station of interest and the RR. M estimator is a robust TF estimator designed to handle large residuals in the regression [Egbert and Booker, 1986a]. M-estimates are a variation of maximum likelihood estimates and seek Chapter 2. Passive geophysical methods to find the model parameters that are most likely given the probability distribution of the data. The solution is similar to the least-squares normal equations, but with an adaptive weight term v added

z = (b H vb) -1 (b H ve) (2.33)
Two weighting functions are proposed, i) Huber and ii) Thomson. The Huber weighting function is defined with the following diagonal elements,

v =    1 if |x| < α α/|x| if |x| ≥ α (2.34)
Typically α is set to 1.5 times the standard deviation. In small residuals i.e. |x| < α the weighted observations are the same as the original observations, but large residuals are weighted down. z is solved iteratively using least squares and then with weights based on the residuals (weighted least squares) until convergence. The worst points are eliminated in the last iterations using Thomson weighting function [START_REF] Egbert | Robust estimation of geomagnetic transfer functions[END_REF],

v = exp(-exp(α(|x| -α)) (2.35)
α=2.8. The Thomson function is more robust, but unlike the Huber function, is does not ensure stability [START_REF] Smaï | Razorback, an open source python library for robust processing of magnetotelluric data[END_REF].

Bounded Influence estimator (BI)

The M estimator handles large residuals, but is highly sensitive to extreme values of the magnetic field (i.e. leverage points). The BI method is a variation of the M-estimator designed by [START_REF] Chave | Bounded influence magnetotelluric response function estimation[END_REF] where the diagonal weighting matrix is enhanced to provide protection against leverage. In the processing done during this Thesis, we use the BI estimator within the Razorback library options. All MT data processed using Razorback during this Thesis can be found in Appendix A.

Razorback run

The Razorback library allows the user to choose classical M-estimator or BI regression with single or multiple combinations of RR. We show in Figure 2.5 an example taken from [START_REF] Smaï | Razorback, an open source python library for robust processing of magnetotelluric data[END_REF] of the processing for one station acquired in Northern Alsace, France. A total of 31 combinations of RR are used in combination with the single site one.

MT forward modeling codes

Two modeling codes for MT were used during this Thesis, i) one-dimensional and ii) two-dimensional modeling. Below we describe the particularities of each these two.

MT forward modeling codes

Two-dimensional modeling: PW2D

For a 2D model, the earth is parameterized by a grid of rectangular prisms; each having a uniform resistivity. To perform the forward calculations, a finite-element code described by Wannamaker et al. [1987] is used, here called PW2D. This program uses a Figure 2.9: Example of response for model in Figure 2.8 obtained using the PW2D code of [Wannamaker et al., 1987]. Four components are output: apparent resistivity and phase for TE and TM modes.

rectangular array of nodes to perform the finite-element calculations, which is called the finite-element mesh. In Figures 2.8 and 2.9 we show an example of a two-dimensional resistivity model and its response using PW2D. The model is a rectangle of dimensions 200 km wide by 10 km depth which includes two conductive bodies (10 Ωm) among a resistive material (1000 Ωm), with receivers located every 20 km. The response includes apparent resistivity and phase for both TE and TM modes. The TE mode is easily identified as the most sensitive to the two conductive bodies.

Ambient seismic noise correlation (ASNC)

Traditionally seismic imaging is based on direct waves emitted by 'artificial' sources (e.g.

air-guns or explosives) or 'natural' sources (e.g. earthquakes). However, with the development of large seismic networks, and access to huge quantities of continuous recordings, a big spectrum of opportunities for imaging without source (passive imaging) have been unlocked. Presently the most popular one is the correlation of long ambient seismic noise records, which allows to reconstruct the impulse response between a pair of receivers, as if a source was placed at one of them.

Here we give a brief introduction to the principles of using seismic ambient noise correlations as virtual seismograms. However, further reading is recommended from seminal articles [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF], [START_REF] Weaver | Diffuse fields in open systems and the emergence of the green's function (l)[END_REF], [Shapiro and Campillo, Chapter 2. Passive geophysical methods 2004], [Shapiro et al., 2005], [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF], [START_REF] Campillo | 1.12 -crust and lithospheric structure -seismic imaging and monitoring with ambient noise correlations[END_REF], among others.

Sources of ambient seismic noise

Ambient seismic noise (ASN) has been studied since early in seismology, and it was soon recognized that the main sources for the permanent agitation of the ground surface are different for different period/frequency bands [START_REF] Gutenberg | Microseisms[END_REF]. Generally, we can separate them in four domains: [2000]. Dashed black lines show the extreme noise models from [START_REF] Peterson | Observations and modeling of seismic background noise[END_REF]. The shaded blocks in green, black, red and blue represent different sources process described in the text.

• For very short periods (below 1 s) human activities (traffic, industries etc) are the primary causes of seismic noise (green shaded area in figure 2.10). They tend to decrease with increasing period and are barely seen at periods > 1 s where a noise local minimum is usually visible in noise spectra. Some natural sources (wind turbulence , rain etc.) can also generate ambient noise in this period range.

• At periods greater than 1 s the noise level increases and peaks at ≈ 7 s (black shaded area in figure 2.10), and is known as the secondary microseismic peak or double frequency peak. It has long been established that it is related to interaction of oceanic waves in deep oceans with similar frequency but travelling in opposite directions [START_REF] Ramirez | An experimental investigation of the nature and origin of microseisms at St. Louis, Missouri: Part two[END_REF].

• The noise peak at ≈ 14 s (red shaded area in figure 2.10) is known as single frequency peak and has a smaller amplitude than the double frequency peak. It also originates from oceans but associated to oceanic waves striking the coast and induc-2.4. Ambient seismic noise correlation (ASNC)

ing a direct transfer of their energy into Rayleigh waves through non linear coupling of waves and bathymetry [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF][START_REF] Stutzmann | GEOSCOPE Station Noise Levels[END_REF]].

• At very long periods (100 -400 s) seismic noise is dominated by infra-gravitational waves linked to the interaction of fluids (atmosphere and oceans) and the Earth [START_REF] Rhie | Excitation of earth's continuous free oscillations by atmosphere-ocean-seafloor coupling[END_REF]. Blue shaded area in figure 2.10

Cross-correlation of noise

The principle of the Ambient Seismic Noise Correlation (ANSC) relies on the demonstration that the impulse response (Green's function) between two points can be retrieved by performing a cross-correlation of the ambient seismic noise records at those two points assuming we can reach an equipartition of the seismic wavefield. This requires a homogeneous distribution of the noise source and/or an efficient diffracting medium. The basis of the approach is now well established by theoretical [START_REF] Weaver | Ultrasonics without a source: Thermal fluctuation correlations at mhz frequencies[END_REF], [START_REF] Weaver | Diffuse fields in open systems and the emergence of the green's function (l)[END_REF], [START_REF] Sanchez | Retrieval of the green's function from cross correlation; the canonical elastic problem[END_REF], [START_REF] Campillo | 1.12 -crust and lithospheric structure -seismic imaging and monitoring with ambient noise correlations[END_REF] and experimental studies [START_REF] Ritzwoller | Eurasian surface wave tomography: Group velocities[END_REF]], [Shapiro and Campillo, 2004], [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF], [START_REF] Lin | Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps[END_REF].

The correlation function recovered at two stations is used as a virtual seismogram, where one station acts as an impulsive source recorded at the other (the receiver). In time domain, for two signals of duration T, the correlation of a signal f at receiver r 1 and g at receiver r 2 is given by

C 1,2 (t) = 1 T T 0 = f (τ )g(t + τ )dτ (2.36)
similarly, in the frequency domain

C 1,2 (w) = f (w) * g(w) (2.37)
Correlation is commonly a way to measure resemblance of two signals. Noise correlation methods are based on the possibility of extracting the coherent part of the signal that contains information about the propagation of the wave between the stations. The crosscorrelation function allows to extract this information without the knowledge of position and time of sources [Lehujeur, 2015].

Practical limitations

Figure 2.11 exemplifies how a uniform distribution of sources lead to a nearly perfect reconstruction of the Green's function of the medium. However, in practical applications, the sources are likely to be distributed non uniformly and the exact response of the Earth is not likely to be recovered. It is therefore important to understand the effect of an non even distribution of sources in the cross-correlations.

Figure 2.12 shows a similar test as that shown previously in Figure 2.11 but now with the sources located at one segment only leading to a poor azimuthal coverage. Looking at Figure 2.12 b), c) and d) we see the impact in the poorly focused cross-correlations and the temporal symmetry observed in the previous experiment is not recovered.

We can point out some effects of an heterogeneous distribution of sources that have been described in the literature:

2.4. Ambient seismic noise correlation (ASNC)

Figure 2.12: Snapshots of cross-correlations as in Figure 2.11, but with a source distribution limited to a line. Reproduced from [START_REF] Paul | Empirical synthesis of time-asymmetrical green functions from the correlation of coda waves[END_REF] • The amplitudes of the correlation for the causal and acausal sides are different [START_REF] Paul | Empirical synthesis of time-asymmetrical green functions from the correlation of coda waves[END_REF]. The asymmetry of the correlation is explained by the effect that the energy contribution is not the same in each side.

• The envelope of the signal can be shifted and induce an error in the estimate of the group velocity [START_REF] Pedersen | Influence of the seismic noise characteristics on noise correlations in the baltic shield[END_REF]] [Lehujeur et al., 2016].

• The phase of the signal can also be affected which directly translates to a bias in the estimation of the phase velocity. Some approaches have been proposed to estimate the bias [START_REF] Weaver | On the correlation of non-isotropically distributed ballistic scalar diffuse waves[END_REF].

This Section aimed to expose the limitations of why in general, cross-correlations of noise can not be considered directly as actual Green's functions for imaging purposes and the noise sources have to be considered as observables such as arrival times and amplitudes [START_REF] Campillo | 1.12 -crust and lithospheric structure -seismic imaging and monitoring with ambient noise correlations[END_REF].

Chapter 2. Passive geophysical methods

ASNC data workflow 2.5.1 Surface wave dispersion curves (SWD) from ASNC: data processing

We have seen that cross-correlations of ambient noise can provide an estimate of the Green's function between the stations. The Green's functions are composed of both body and surface waves. Due to their high amplitude, surface waves represent a high content on the signal, in this Thesis we focus on these type of waves only.

During the scope of this Thesis, we did not cover the processing of ambient seismic noise to recover surface-wave dispersion curves on field data. Therefore the purpose of this section is just to give a brief introduction on how to obtain surface-wave dispersion measurements from ambient seismic noise. Below, we describe the main steps of an standard workflow, for more details the reader should refer to [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF].

Phase 1: Pre-processing

An obvious difficulty when processing ambient seismic noise data is the presence of strong energy transients that can dominate in the correlation. The objective of the pre-processing Phase 1 is both to mitigate the effects of strong directive sources, (large earthquakes, local storms or others) using time domain normalization and to normalize the amplitude spectrum to mimic impulsive sources (spectral whitening).

This process is typically applied to day-long data series.

Phase 2: Cross correlations

The output of Phase 1 are daily pre-processed time series, the next step are the computations of cross-correlations. Cross-correlations are usually performed daily in the frequency domain and then stacked. The purpose of stacking the cross-correlations is to increase the signal-to-noise ratio (SNR) and converge towards the Green functions. Depending on the frequency band of interest and the distribution of noise sources, it may require weeks to years to stabilize the cross-correlations. The output of Phase 2 is an estimate of the Green's function. For a network of n stations, we end up with a total of n(n -1)/2 cross-correlations.

Phase 3: Dispersion measurement, Frequency Time ANalysis (FTAN)

During Phase 3 group and/or phase velocity are measured on the stacked cross-correlations using traditional FTAN (Frequency Time ANalysis) analysis [START_REF] Levshin | Automated Detection, Extraction, and Measurement of Regional Surface Waves[END_REF] [ Bensen et al., 2007] [Keilis-Borok and[START_REF] Keilis-Borok | Seismic surface waves in a laterally inhomogeneous earth[END_REF] ) is computed from the peak of the envelope so that the group velocity is U (ω 0 ) = r/τ (ω 0 ) where r is the offset distance (distance between the stations).

The output of this step are group/phase dispersion curves for each Green's function estimate provided from Phase 2. Analog to MT where phase and apparent resistivity curves are the data used for inversion, for this method surface-wave dispersion curves (group or phase) are the observed data to use in inversion. In this Thesis, we use the disba module implemented by Luu [2019]. It is a Python library for one-dimensional modeling of surface-wave dispersion curves that implements the codes from Computer Programs in Seismology CPS [Herrmann, 2013]. Sensitivity kernels in Figure 2.14 indicate that sensitivity of Rayleigh-wave group velocity for compressional velocity and density is less compared to shear velocity, as has extensively been discussed in literature e.g. [START_REF] Pan | Sensitivity analysis of dispersion curves of Rayleigh waves with fundamental and higher modes[END_REF]. In this Thesis, we invert only for shear velocity and assume a known value or relationship for the other parameters.

Chapter 2. Passive geophysical methods 2.6. SWD forward modeling codes

Three-dimensional modeling: SW4

Modeling and sensitivity analysis of surface waves are generally done based on the onedimensional approach of Herrmann [2013] mentioned previously. As part of this Thesis we were interested to directly evaluate the sensitivity of surface waves in models with lateral changes (e.g. 2D or 3D models). We use SW4 2.01 [Petersson and Sjögreen, 2015] published under the GPL 2 license to produce synthetic seismograms for such 2D/3D models in Chapter 5. SW4 is a program for simulating seismic wave propagation by solving the elastic and visco-elastic wave equations in displacement formulation. It is preferably run on parallel computers, and supports fully 3D heterogeneous models.

We have developed a workflow to produce group velocity dispersion curves maps from SW4 models and is summarized here in three main steps: i) waveform modeling using the elastic code SW4 [Petersson and Sjögreen, 2015], ii) computation of dispersion curves using FTAN and iii) regionalization process. Figure 2.16 illustrates these three main steps.

Waveform modeling with SW4

The cross-correlations of experimental ambient seismic noise records provide an estimation of the impulse response between two stations. In other words, the seismogram recorded at Station 2 if the source was located at Station 1. In order to replicate the effect of having a source in each receiver, we propose to compute n models using SW4 by locating the sources in each of the n receivers each time. That is, if n receivers are deployed, we obtain a total of n 2 seismograms (.sac files). Every time we compute a model response, only one source is active, so that no simultaneous sources occur. Chapter 2. Passive geophysical methods

Dispersion curves computation using FTAN

For each of the receiver and sources combinations FTAN analysis as described in Section 2.5.1.3 is performed in order to compute the group or phase velocity for each period. We use the ftanpy Python module from [START_REF] Ensing | Probing the subsurface of the auckland volcanic field with ambient seismic noise[END_REF] for the generation of dispersion curves. We compute a total of n(n -1) group dispersion curves, where n is the number of receivers deployed.

Regionalization

The inverse problem appears in two steps when imaging the subsurface using surfacewave dispersion curves recovered from cross-correlation of ambient seismic noise (or any travel time data). A first inversion is generally used to estimate local dispersion curves per period, from a distribution of group or phase velocity dispersion curves, this process is sometimes called regionalization, and it is a linear inverse problem. The second step where inversion appears is to solve for the shear velocity with depth from each local dispersion curve, it is a non-linear inversion that is discussed in Chapters 4, 5 and 6.

Generally, dispersion measures are computed in a far-field assumption i.e. where source-receivers distances are greater than a certain number of wavelengths (λ). In the case of ambient seismic noise tomography this condition is generally fixed at 3λ ( [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF], [START_REF] Lin | Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps[END_REF]), in other words, all pairs of stations with offset < 3λ are ignored.

A forward problem is defined as d = Gm [Tarantola, 2005], where d is the observable data, m the model parameter, and G a linear operator. During regionalization, we look for m i.e. the slowness (either relative or absolute) at a given period for each position (pixel) of a pre defined grid. The linear inverse problem is solved as m = G -1 d, with m the slowness to find, the matrix G -1 is the generalized inverse of the path traversed by the wave between the pair of stations and d are the observed slowness (reciprocal of the group velocities) for each pair of stations at the given period. This process is repeated for each period. The result are dispersion curves for each pixel of a grid, that are later to be used for depth inversion (non-linear inverse problem).

Conclusion

In this Chapter, we have introduced the two passive geophysical methods that are implemented in a joint inversion framework in this Thesis i.e. MT and SWD from ASNC. We have reviewed the basis of both methods and described the forward modeling implementation for each of them.

In the next Chapter we will describe several approaches for performing multi-physics joint inversion and how we use these two methods in a joint inversion framework. In summary:

• MT data is used to study the Earth's electrical conductivity (or the inverse i.e.

resistivity)

• ASNC method is used here to obtain surface-wave dispersion maps for building shear velocity models.

• SW4 will be used to compute synthetic seismograms for sensitivity analysis of surface waves using three-dimensional models. Multi parameter joint inversion methods are generally grouped in two classes, structural and petrophysical approaches. In this Chapter we review some of the main methods of these two classes. From the structural approach i) the very well established cross-gradients approach introduced by Gallardo and Meju [2003], ii) the dot product [Molodtsov et al., 2011], [START_REF] Tarits | Mapping density models onto resistivity structure through joint inversion of gravity and mt[END_REF] and iii) gramian constraints [Zhdanov et al., 2012]. On the petrophysical approach we review two of the newest approaches i) the correspondence maps (CM) approach from Carrillo andGallardo [2018] [Carrillo et al., 2021] that is implemented in this Thesis in Chapters 4, 5, and 6, and ii) variation of information

Structural approach

Structural joint inversion methods assume that it is the structure geometry that controls the distribution of petrophysical properties [Gallardo and Meju, 2003]. They require specification of mathematical conditions at structural boundaries in a stable way. Today the most common structural coupling method is the cross-gradients constraint, introduced by Gallardo and Meju [2003].

Chapter 3. Joint inversion methods

Cross-gradients

Generally, properties vary with position and these variations can occur in any direction. At any position, the variation can be characterized by the magnitude and direction of change.

These attributes are represented mathematically by the vector field of the gradients of the geophysical properties. Several joint inversion approaches rely on gradient-based relationships, e.g. the cross-gradient constraint [Gallardo and Meju, 2003]. This approach seeks structural similarity in the images exploited by different geophysical methods. They introduce the cross-gradient function given two model parameters m 1 and m 2 in a 2D model in the XZ plane as

τ (m 1 , m 2 ) = ∇m 1 × ∇m 2 = |∇m 1x ∇m 2z -∇m 1z ∇m 2x | (3.1)
Where ∇m x and ∇m z are the gradients for x and z directions respectively, computed for each model parameter. The cross-gradients criterion requires the problem to satisfy the condition τ = 0, i.e. where any spatial changes occurring in both m 1 and m 2 must point in the same or opposite direction irrespective of the amplitude. In a geological sense this implies that if a boundary exists, then it must be sensed by both methods in a common orientation. An additional flexibility of the technique is that the cross-gradients constraint is also satisfied when either ∇m 1 or ∇m 2 vanishes in some part of the model, thus giving the models the possibility of admitting a geological boundary which has a significant change only in the electrical resistivity or seismic velocity of the adjoining rocks [Gallardo and Meju, 2003]. The cross-gradients approach is flexible and can be applied to several model parameters, not being restricted to two models.

Dot products of image gradients

Noting that the cross-gradient method steers the gradients of the two subsurface models toward co-linearity (either parallel or anti-parallel), Molodtsov et al. [2011] have proposed a variation on the above method. They argue that the correlation between velocity and resistivity is usually positive, i.e. velocity increases with increasing resistivity, and hence parallelism between the gradients should be favoured by the joint inversion. The method is based on the assumption that positive or negative correlation can be a priori distinguished. It differs from the cross-gradient method described above because they use the dot product. [START_REF] Tarits | Mapping density models onto resistivity structure through joint inversion of gravity and mt[END_REF] describe the modified cross-gradients term

τ mod (m 1 , m 2 ) = | sin(θ)| 2 (3.2)
where

θ = cos -1 ( ∇m 1 •∇m 2 |∇m 1 ||∇m 2 |
). The angle θ represents the angle between the vector gradients of m 1 and m 2 , it is zero if the vectors are parallel as in the cross-gradients approach, but it is undefined if either of the gradients is zero.

Petrophysical approach

Gramian constraints

The use of Gramian constraints was first demonstrated in geophysical exploration by Zhdanov et al. [2012]. They introduce the use of a Gramian space of the model parameters with its main property being that the norm of a function in Gramian space provides a measure of correlation between the function and the model parameters.

Gramian constraints are based on the minimization of the determinant (S G in equation 3.3) of the Gram matrix of a system of different model parameters e.g. m 1 and m 2 . The Gram matrix G([m 1 , m 2 ]) is the inner product between the model parameters or gradient of model parameters or some transform parameters [Ogunbo, 2019].

S G = |G([m 1 , m 2 ])| = (m 1 , m 1 ) (m 1 , m 2 ) (m 2 , m 1 ) (m 2 , m 2 ) (3.3)

Petrophysical approach

Joint inversion methods using petrophysical relationships are based in the fact that for some specific geological environments geophysical parameters can be related to some mathematical relationships. Indeed, in many exploration cases, there are down-hole logging data available (e.g. sonic, laterolog) that can be used to obtain a relationship between them, and some degree of proximity to this relationship can be included as an apriori constrain in the objective function of a joint inversion. In the case studies presented here, we assume that subsurface petrophysical relationships are unknown.

Correspondence Maps (CM)

The CM joint inversion was introduced by Carrillo and Gallardo [2018], it is an innovative way to compute joint inversion from a parameter relationship approach, without imposing an specific relationship between the variables. Carrillo and Gallardo [2018] propose a general form of relationship between two model parameters (m 1 and m 2 ), if p and q are the maximum power allowed for each parameter, and a ij are the polynomial coefficients linking m 1 and m 2 , they introduce the function g

g(m 1 , m 2 , a) = q i=0 p j=0 a ij m i 1 m j 2 = 0 (3.4)
From equation 3.4, g is equal to zero in the general case. It is convenient to manipulate it algebraically to avoid the trivial solution, i.e. when all the coefficients describing the polynomial are equal to zero. Carrillo and Gallardo [2018] circumvent this problem by normalizing equation 3.4 by the independent term coefficient a 00 . The normalized equation 3.4 becomes

ĝ(m 1 , m 2 , a) = q i=0 p j=0 âij m i 1 m j 2 = 0 (3.5)

Chapter 3. Joint inversion methods

where, âij = a ij a 00 ; a 00 ̸ = 0; â00 = 1 (3.6) Below we show an example for degree 1 i.e. p and q are equal to 1

ĝ(m 1 , m 2 , a) = â00 + â01 m 0 1 m 1 2 + â10 m 1 1 m 0 2 + â11 m 1 1 m 1 2 = 0 (3.7)
Rearranging and substituting â00 of equation 3.6 in equation

3.7 ĝ(m 1 , m 2 , a) = â01 m 0 1 m 1 2 + â10 m 1 1 m 0 2 + â11 m 1 1 m 1 2 = -1 (3.8)
Before performing a full joint inversion i.e. to solve for m 1 , m 2 and a ij , we assess how the ĝ function behaves if we have as known parameters m 1 and m 2 and we only solve for the coefficients (a ij ) in the inversion. Following Tarantola [2005] notation d = Gm, d is the observed data i.e. a vector of -1, G are the polynomial terms formed by ĝ(m 1 , m 2 , a ij )

and m are the coefficients a ij to be found. We then invert for the coefficients using the standard form m = G -1 d. Two experiments are performed and described below in Subsections 3.2.1.1 and 3.2.1.2.

Separate inversion of the coefficients: circular test

In this first experiment we look at a circular function given by

-m 2 1 -m 2 2 = -1 (3.9)
The explicit form of equation 3.5 when p = 2 and q = 2 becomes

ĝ(m 1 , m 2 , â) = â01 m 2 +â 02 m 2 2 +â 10 m 1 +â 11 m 1 m 2 +â 12 m 1 m 2 2 +â 20 m 2 1 +â 21 m 2 1 m 2 +â 22 m 2 1 m 2 2 = -1 (3.10) Coefficients â01 â02 â10 â11 â12 â20 â21 â22 True 0 -1 0 0 0 -1 0 0 i) 0 % noise -5.5e-17 -1 1.63e-16 0 -5e-16 -1 1.1e-16 1.54e-9 ii) 5 % noise 0.02 -0.92 0.01 -0.008 -0.02 -0.95 -0.04 -0.05
Table 3.1: Coefficients recovered from inversions described in Section 3.2.1.1. The three rows refer to the true coefficients, the first test of inversion i) with no noise added to the parameters m 1 and m 2 and the second test ii) with a 5 % noise added. Highlighted in gray are the coefficients that should be different from zero.

If we compare equation 3.10 with equation 3.9 we can see that all coefficients in 3.10 should be equal to zero except for two, â02 and â20 i.e. the ones multiplying the terms m 2 1 and m 2 2 . We propose to do two tests of inversion to recover the coefficients a ij knowing m 1 and m 2 . The first test i) consists of using the parameters m 1 and m 2 that perfectly describe the circumference function (blue graph a) of Figure 3.1), i.e. with no noise added, and use the inversion to find the coefficients. In the second test ii) we use as input observe that the recovered coefficients slightly differ from the true ones, which translates in larger errors on the data fit, nevertheless the data remains within an error of 5 % i.e -1.05 < g < -0.95. In table 3.1 we summarize the coefficients that were recovered in both tests, compared to the true ones; and we observe that the inversion is stable in the first test i), while in the second test ii) the coefficients show a higher error, as expected.

Separate inversion of the coefficients: well-log test

In the next tests we analyzed the inversion for the coefficients when a relationship is possible to recover given the pairs of m 1 and m 2 . Then, we study what happens when there is no relationship between them, i.e. they behave randomly. We remind the reader, that in both of these tests as for the previous one, m 1 and m 2 are known and the inversion only focuses on recovering the coefficients a ij

i) Original well-log

In this case, we use data downloaded, publicly-available, from the Kansas Geological Survey database (http://www.kgs.ku.edu/Magellan/Logs/). Three parameters were studied, shear velocity, resistivity and density, with the corresponding cross-plots displayed in Figure 3.2 (upper pannels). Therefore three inversions were computed for the different combinations of model pairs i.e. 1) resistivity-velocity, 2) resistivity-density and 3) velocity-density. Since in equation 3.5 a key parameter is the definition of the degree, Chapter 3. Joint inversion methods we performed one set of inversions using degree 1 with the form ĝ = â01 m 2 + â10 m 1 = -1, and a second set using degree 2 with the form ĝ = â01 m 2 + â02 m 2 2 + â10 m 1 = -1. A first glance at the cross-plots in Figure 3.2 shows that there appear to be relationships between the properties and they seem related to depth. The results of the inversion, (three plots in right side of Figure 3.2), show that by using a degree 1 polynomial relationship, it is hard to explain the behavior of the parameters, however by using a degree 2 polynomial, the three combinations show a better fit with the target ĝ = -1. From the three combinations, we can observe that resistivity-velocity shows a higher misfit than the other two combinations that have density as common parameter. One reason could be the lower range of variability of density (1.8-2.8 g/cm 3 ), which translates to 'simpler' relationships, compared to the resistivity-velocity (first cross-plot in Figure 3.2) or a higher coupling between these two parameters. (bottom). Right, the data from inversion for the three combinations of pairs (resistivitydensity, density-resistivity and velocity-density), in orange using p = 1 and q = 1 and in green with p = 2 and q = 2 in equation 3.5, black line is the observed data, i.e. ĝ0 = -1

ii) Random well-log

To verify that CM will not force a relationship when one is not warranted by the 3.2. Petrophysical approach observations, we performed the same experiment as before but assuming a random data set, i.e. no obvious relationship exists (Figure 3.3). The results of the inversions show that for none of the three combinations (resistivity-velocity, density-resistivity, velocitydensity), neither a degree 1 or degree 2 polynomial can provide a relationship between the parameters, resulting in a poor data fit. This result reassured us that indeed, when no relationship is present, CM will show a very poor fit. The risk of producing a meaningless relationship is minimal, at least for polynomials up to degree 2.

Variation of Information (VI)

Variation of Information (VI) is an unsupervised machine-learning method which constructs, where possible, a one-to-one relationship of an a-priori unspecified form between the physical properties under consideration [Moorkamp, 2021a]. It is a theoretical measure of the amount of information contained in property m 1 about property m 2 and very closely related to the concept of Mutual Information (MI) [Haber andGazit, 2013] [Mandolesi andJones, 2014]. Although widely used in medical imaging, it had not been implemented successfully in geophysics until very recently [Moorkamp, 2021b] [Moorkamp, 2021a]. The VI between two properties m 1 and m 2 is defined as

Chapter 3. Joint inversion methods V I(m 1 , m 2 ) = 2H(m 1 , m 2 ) -H(m 1 ) -H(m 2 ) = H(m 1 ) + H(m 2 ) -M I(m 1 , m 2 ) (3.11)
Here H(m) =i p(m i ) log p(m i ) is the Shannon entropy and p(m i ) is the probability density function and MI is the mutual information term M I( m 1 ,m 2 ) = H(m 1 ) + H(m 2 ) + H(m 1 ,m 2 ). The joint inversion tries to minimize the VI term. shows that MI is large (1.47) when knowledge of variable m 1 provides significant information on variable m 2 . On the contrary, when the variables are randomized as in Case 2 in Figure 3.4, the computed MI is low (0.17), meaning that the knowledge of m 1 does not provide significant information about m 2 .

Conclusion

We have reviewed several joint inversion approaches from the two main classes i.e. structural and petrophysical. For this Thesis, we have chosen to implement the correspondence maps due to its adaptability to perform joint inversions with one-dimensional models compared traditional cross-gradients or Gramian constraints applied to 2D/3D models.

In the next Chapter we show the joint inversion methodology implemented within this Thesis using CM.

We have been able to recover polynomial coefficients (e.g. circular test and well-log test) and showed that for random well-log data, no relationship is possible to find, as one would expect. Due to time limitations we did not try the VI approach recently published by Moorkamp [2021a] but it is indeed encouraged to compare these two approaches i.e.

CM and VI joint inversions.

In summary:

• Among several joint inversion methods, we implement the correspondence maps due to its flexibility to be adapted from one-dimensional to higher dimensional models and the advantage to be implemented without prior knowledge about the petrophysics of the models.

Bibliography

Carrillo, J. and Gallardo, L. A. (2018). Joint two-dimensional inversion of magnetotelluric and gravity data using correspondence maps. Geophysical Journal International, 214(2):1061-1071. Carrillo, J., Perez-Flores, M. A., Gallardo, L. A., and Schill, E. (2021). Joint inversion of gravity and magnetic data using correspondence maps with application to geothermal fields. Geophysical Journal International, 228(3):1621-1636. [START_REF] Gallardo | Characterization of heterogeneous near-surface materials by joint 2d inversion of dc resistivity and seismic data[END_REF]. Characterization of heterogeneous near-surface materials by joint 2d inversion of dc resistivity and seismic data. Geophysical Research Letters, 30(13). This Chapter is based on a published paper [START_REF] Aquino | Joint one-dimensional inversion of magnetotelluric data and surface-wave dispersion curves using correspondence maps[END_REF] (appendix D) and an oral presentation [START_REF] Aquino | Joint inversion of magnetotelluric and ambient seismic noise data using cor-respondence maps[END_REF].

Introduction

In the previous Chapter we have introduced the CM term i.e. ĝ(m 1 , m 2 , a) function in equation 3.5. We have made several tests inverting only the coefficients given fixed parameters m 1 and m 2 . In this Chapter, we invert for all three model parameters of the ĝ(m 1 , m 2 , a) function. We present a joint inversion methodology that includes the correspondence map approach described by Carrillo and Gallardo [2018], implemented as an Occam-type inversion [Constable et al., 1987]. Our approach, in addition to inverting the data using geometric constraints, seeks simultaneously for any pair of model parameters Chapter 4. One-dimensional approach a functional relationship. The relationship is treated as a random variable as part of the joint inversion algorithm.

Our work focuses on exploring the use of CM to invert for the 1D shear velocity from dispersion curve data (from ambient seismic noise) and the resistivity (from magnetotelluric data). We first illustrate the method on synthetic data, first and higher-order relationships, and then we apply our methodology to field data. Throughout the Chapter, we use the same notation as Carrillo and Gallardo [2018] and Constable et al. [1987] to help the reader relate our work to these seminal articles.

The method: joint inversion using CM

The joint inverse problem including correspondence maps using Occam's approach is defined as finding the simplest (i.e. smoothest) shear velocity (m 1 ) and resistivity (m 2 ) models that can be correlated by a polynomial relationship and reproduce the observed data. The problem is set to find m 1 , m 2 and the normalized coefficients âij . Therefore we propose the objective function as an expanded version of the matrices described in the main equation of Occam, i.e.

m i = [λ∂ T ∂ + (W J) T W J] -1 (W J) T W (d -F [m i-1 ] + Jm i-1 ) (4.1)
Here m i is the model parameter to find at iteration i and m i-1 the model at previous iteration, λ∂ T ∂ is the smoothing matrix (here the difference in physical property between consecutive layers); λ is the commonly known Lagrange multiplier, a large λ brings a smooth model; J the sensitivity matrix, d and F [m] are the observed and calculated data respectively. W is the data weighing matrix, modified from the one in Constable et al.

[1987] by including a factor F w , to weigh for each of the three data types. This factor ensures that a data type with more elements than the others will not inherently bias the inversion. W is then defined W = diag(F w /σ), where σ is the uncertainty associated to each data value and F w :

F w = A D N N D (4.2)
where A D is the weight given to data type D with 3 D=1 A D = 1. N D is the number of data for each type and N =

X = ||W d -W F [m]|| (4.
3)

The modified model matrix m is defined as a column vector that comprises shear velocity, m 1 , resistivity, m 2 , and the normalized polynomial coefficients, âij

4.3. Synthetic tests m =    m 1 m 2 â    (4.4)
Similarly the observed data d and predicted F [m] are extended to

d =    d 1 d 2 -1    (4.5) F [m] =    F SW D [m 1 ] F M T [m 2 ] ĝ[m 1 , m 2 , â]    (4.6)
J is the general Jacobian matrix which integrates the sensitivity of the given model parameters. More details regarding the Jacobian matrices of the correspondence map function ĝ can be found in Carrillo and Gallardo [2018]. The complete J matrix given the model parameters in equation 4.4 and data in equation 4.5 is

J =    J[m 1 ] 0 0 0 J[m 2 ] 0 J[g m1 ] J[g m2 ] J[g a ]    (4.7)
From equation 4.1, the first term (λ∂ T ∂) relates to the smoothness of the model parameters. In our case of joint inversion using correspondence maps, there is no physical or mathematical reason to smooth the coefficients, and hence a smoothness of zero is imposed to the coefficients term. The smoothing term in equation 4.1 is then,

∂ =    ∂m 1 ∂m 2 ∂a    =    ∂m 1 ∂m 2 0    (4.8)
The experiments presented in the next sections include the retrieval of all parameters described in equation 4.4. The aim is to compare separate versus joint inversions for shear velocity and resistivity. The workflow followed in all experiments is described in Figure 4.1

Synthetic tests

To verify the reliability and stability of the inversion algorithm described in the previous Section, we test several scenarios of correspondence map relationships using synthetic data. 1) first-order relationship, 2) second-order relationship and 3) unknown relationship from well-log data. All experiments start with an homogeneous model that does not have to be close to the true model. For the case of the coefficients, the initial model

m 2 = 6.13m 1 + 1.29 (4.11)
i.e. very close to the true relationship in equation 4.9.

Higher-order relationship: synthetic test

The second experiment consists of a resistivity model (m 2 ) computed from a set of shear velocity values (m 1 ) using a degree 2 polynomial relationship. This type of relationship has been observed in geothermal environments, where shallow formations dominated by sediments show a decrease trend in resistivity and an increment towards deeper granitic

Synthetic tests

formations with a persistent increase in the velocity through the whole section resulting in higher-order parameter relationships [Glass et al., 2018]. Carcione et al. [2007] showed several examples of this type of relationships between velocity and conductivity in oilsaturated shales and sandstones.

The true relationship from where m 2 is computed for this experiment is m 2 = 26.0m 2 1 -23.0m 1 + 6.00 (4.12)

For this experiment, the normalized correspondence map function to search has the form

ĝ(m 1 , m 2 , a) = â20 m 2 1 + â10 m 1 + â01 m 2 = -1 (4.13)
After 15 iterations the recovered relationship is

m 2 = 25.7m 2 1 -21.5m 1 + 5.55 (4.14)
The correspondence map plot in Figure 4.3 shows that the joint inversion (magenta) from equation 4.14 is more consistent with true relationship (black) from equation 4.12 than the results from separate inversions after 25 iterations (blue) that are more scattered, especially for the deeper section of the shear velocity model.

Discussion on first and higher-order relationships from synthetic tests

To describe how accurate the joint and separate inversions results from the synthetic models are, we assume three factors of the inversion algorithm that are inherently correlated:

smoothing, correspondence maps (resulting models) and the data (misfit). In the joint inversion using correspondence maps, we were able to retrieve parameter relationships without compromising the smoothing and the data fit.

In Figure 4.4 we can see how joint inversion results in magenta represent a closer shape to the true relationship of the models than those of separate inversions for both first and second degree synthetic tests. In order to quantify how close the correspondence maps obtained from separate and joint inversions were to the true relationships, we compute the ĝ function using the pairs of models retrieved from each inversion and compare it to the target ĝ = -1 within a certain error or tolerance. We obtained, for the first-order relationship test, that 100% of the joint inversion model-pairs (shear velocity-resistivity) have ĝ values, computed using the true relationship, between -1.05 and -0.95 vs only 15% of the separate inversion pairs. In the second degree relationship synthetic test, it is the case for 80% of the joint inversion pairs but only 45% of the separate inversion pairs. These joint inversion results lead to a meaningful improvement on the modelsrelationships without compromising the data fit. Hence we can say that joint inversion using correspondence maps leads to a reduction in the number of acceptable models that could explain the same observed data ('reduction of non-uniqueness').

Unknown relationship: well-log data

The third test consists of a numerical experiment using subsurface parameters extracted from well-log data. This test is approached by a two-step procedure. The first step refers to the problem of defining a proper degree of the polynomial for joint inversion and the second step to the actual results of the joint inversion. The data from well API 15023214760000 was downloaded from the publicly-available Kansas Geological Survey database (http://www.kgs.ku.edu/Magellan/Logs/). Amongst other logs, it includes deep resistivity (DIL) and compressional slowness (DT) logs. For the purpose of this exercise, shear velocity is computed using a constant V p /V s ratio of 1.7. Contrary to the previous two experiments, the ρ-V s relationship here is unknown, hence closer to a field data case.

Since the correspondence map joint inversion requires the definition of the maximum degree (p and q) in equation 3. After analysing the degree of polynomial to be used, the second step consists of the actual joint inversion using correspondence maps. For modelling purposes, the well-log models are binned into layers 150m thick to make the computations faster. In Figure 4.6 it is noted that RMS misfits for the three data types (MT, SWD curves and g) marginally decrease after the fifth iteration for the joint inversion test. From the parameter side, the shallow section represents a model where velocity increases and resistivity decreases, while in the deeper section both properties increase. The joint inversion approach tries to couple these two behaviors as a parabola form. The recovered relationship for the joint inversion scenario d) after 15 iterations is:

m 2 = 5.32m

Field example

The previous experiments were necessary to asses the reliability of our methodology. The datasets tested in this Section consist of one pair of MT sounding from the ECOGI project [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF] and one fundamental mode Rayleigh-wave group dispersion curve from the EstOF network [Lehujeur et al., 2018] [Lehujeur et al., 2015]. The selected pair of stations (MT-SWD) have a surface separation of at most 150 meters, for this purpose we consider them to be coincident. The processed MT data cover a frequency range from 0.01 to 100 Hz [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF]. The seismological data processed by Lehujeur et al. [2015] consist of fundamental and first overtone Rayleigh-wave dispersion curves in the range of 0.2 o 1.28 Hz (0.7-5 s).

Generally, all components of the MT impedance tensor are non-zero and require a 3D modeling approach to be reproduced. Since the scope of this Chapter is 1D, we choose an MT station for which the phase tensor ellipse [Caldwell et al., 2004] is closest to a circle over the broad frequency range of interest, i.e. suitable for our 1D modeling and inversion approach.

A normalized correspondence map function like equation 4.13 is used for this data set, since the results from separate inversions (blue in Figure 4.8) plot as a parabola in a cross-plot. Models obtained by [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF] and Lehujeur et al. [2018] are displayed as reference in gray. Note that resistivity model from [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF] was obtained from a 2D inversion using MARE2DEM [Key, 2016], displayed in Figure 4.8 is the extracted model at the station selected in this joint inversion test. The shear velocity model displayed as reference, was obtained from Lehujeur et al. [2018], he used a Monte Carlo inversion approach from where the median of the best 2000 models was retained as the solution of the inversion.

The results of the joint inversion (Figure 4.8) confirm the relevance of a parabolic parameter relationship. The introduction of the correspondence map in the joint inversion algorithm results in a shear velocity model below 2.5 km depth closer to the shear velocity model from Lehujeur et al. [2018], compared to that obtained by separate inversion. Note however that in Lehujeur et al. [2018], bounds for shear velocity based on well-log data were used as input due to the nature of the Bayesian approach used in their work. In our case, we did not provide any prior information or bounds for shear velocity and resistivity and half-space initial models were used for both parameters. The recovered relationship after 15 iterations was:

m 2 = 2.49m 2 1 -1.

Effect of the degree of ĝ(m 1 , m 2 , â) in the joint inversion

We have seen that the degree of the polynomial relationship used for the joint inversion is an important parameter. Therefore it is worth to investigate how to choose this parameter and how is the joint inversion affected by it. This study started with joint inversions from synthetic models and responses, where predefined relationships were used, i.e degree 1 and higher. In each of these two synthetic tests, the chosen relationship for the joint inversion was inherently used as the 'correct one', for instance, the higher degree relationship joint inversion used a degree 2 relationship. To study the effect of using different degrees, we took the same synthetic models from Section 4.3.2 and perform joint inversions using purposely a wrong degree for relationship. That is, for Section 4. 3.2, where the known relationship is of degree 2, we perform a joint inversion using degree 1. Figure 4.9 shows the cross-plot ρ-V s and the RMS evolution at each iteration for this test. We observe that the joint inversion using a degree 1 does not converge like when using a degree 2, confirming that the ρ-V s relationship can not be explained by a low order function, and therefore a higher order is required. This suggests that when the relationship is unknown, a higher-order relationship should be initially tested; following with decreasing degrees to simplify the solution. If the simpler relationship can explain the data, the joint inversion will converge, otherwise, the higher-order relationship must be kept.

As logical approach to decide the degree of the joint inversion, one can think of looking at the correspondence maps from separate inversions, and use this as a starting point to estimate the degree to use. If we look at the field example in 4.4, both, the reference models and the models from separate inversions describe a somehow parabolic relation-Figure 4.9: Results of Joint inversion of data in Section 4.3.2 using degree1 and degree2, correspondence maps (left) and RMS misfit evolution (right). Non convergence is observed when performing joint inversion using degree 1.

ship and this was the reason why the degree 2 was chosen for joint inversion. Appendix 4.B shows the result of the joint inversion for this data set but using a degree 1 relationship instead. Figure 4.11 shows how the fit of the degree 1 polynomial relationship joint inversion for the three datasets (dispersion curves, MT: apparent resistivity/phase and ĝ)

is significantly lower than when using a degree 2.

Conclusions

To the best of our knowledge, we have successfully applied a correspondence map joint inversion between surface-wave dispersion curves and MT data for the first time. It is an effective way to find meaningful physical parameter relationships while retrieving the shear velocity and resistivity models. The inclusion of correspondence maps in the joint inverse problem has shown to converge without the definition of prior bounds as input and by using homogeneous initial models.

The results of the first two tests using synthetic data show how the recovered relationships were essentially the same as the true relationships. They also prove, that the models retrieved from joint inversion are importantly closer to the true models than those obtained from separate inversions. We observed when evaluating the ĝ function using model-pairs from joint and separate inversions, a significant higher percentage of modelpairs from joint inversion had a ĝ between -1.05 and -0.95, reducing the number of acceptable models without compromising the data fit ('reduction of non uniqueness'). The first part of the well-log example, the assessment of separate inversions of only the coefficients, demonstrate an effective way to evaluate and select an appropriate polynomial degree for the joint inversion.

The field data experiment illustrates how joint inversion with correspondence maps achieve a very similar velocity model to the one obtained by a Bayesian approach [Lehujeur et al., 2018], without forcing bounds to the model, but rather by imposing the existence of a parameter relationship. It exemplifies how, if a cross-property mathemati-

4.A. Stability of correspondence map joint inversion

cal correspondence exists that can explain observed apparent resistivity and dispersion curves, the resulting models have a higher coupling than when not using the correspondence map term. The addition of a proper correspondence map term has shown not to compromise the data fit in the joint inversion algorithm. It is essential however to remind the reader that the convergence of the joint inversion is inherently related to the appropriate selection of the polynomial degree, as described in Section 4.5.

4.A Stability of correspondence map joint inversion

In order to evaluate the stability of the joint inversion algorithm for a degree 2 using the full expression of ĝ (p = 2 and q = 2), we perform a joint inversion test using the data and model responses from the well-log data in Section 4. The ĝ functions used in these two joint inversion tests are explicitly described below.

a) Constrained degree 2:

ĝ(m 1 , m 2 , â) = â20 m 2 1 + â10 m 1 + â01 m 2 = -1 (4.17) b) Non constrained degree 2: ĝ(m 1 , m 2 , â) = â01 m 2 +â 02 m 2 2 +â 10 m 1 +â 11 m 1 m 2 +â 12 m 1 m 2 2 +â 20 m 2 1 +â 21 m 2 1 m 2 +â 22 m 2 1 m 2 2 = -1 (4.18)
The results in Figure 4.10 show that the joint inversion can be stable using a restricted number of coefficients, i.e. the three coefficients in equation 4.17), while it can take more iterations for the non constrained form to converge or even not converge. For this test, after 15 iterations the non constrained inversion did not converge to a solution like the constrained inversion. Carrillo and Gallardo [2018] ran into the same problem: they describe instability in the joint inversion of two-dimensional models (gravity and MT) when using a large number of coefficients i.e. eight coefficients (p = 2 and q = 2).

4.B First vs higher-order correspondence map joint inversion in field data

We are also interested to evaluate the influence of the chosen degree of the polynomial relationship (p and q values in equation 3.5) on the correspondence map joint inversion. We hence perform two more joint inversion tests on the same field data: i) correspondence Figure 4.10: RMS misfit evolution for constrained (left) and non constrained (right) joint inversion of higher-order relationship test map joint inversion using a maximum degree of 1 (linear relationship) and ii) degree 2 (parabolic relationship). The parabolic relationship is described in Section 4.4. It is expected that using a higher-order relationship would show better results than including a degree 1 relationship in the joint inversion, since the reference models do not show any hint of a first-order type of relationship.

The results are consistent with our expectations. The degree 1 relationship correspondence map is never able to converge into a solution that could simultaneously fit the data while finding a relationship. This indicates that even if the models can be improved via joint inversion using correspondence maps, the convergence of the inversion is restricted to a proper selection of the degree of the relationship. The model parameters of this field data set require a higher polynomial degree in order to fit both observed data sets (d 1 and d 1 ) while explaining a common relationship as observed in Figure 4.11.

4.C Weight analysis of the data in the joint inversion framework

We investigate the effect of the weight, A D in equation 4.2, of each data set in the joint inversion framework. Four scenarios of distribution of weights are tested and described in table 4.1. The observed data sets are based on the synthetic data sets of Section 4.3.2.

Test A 1 = W eight(SWD) A 2 = W eight(MT) A 3 = W eight(Corr. Map) 1 1/3 1/3 1/3 2 1/2 1/4 1/4 3 1/4 1/2 1/4 4 1/2 1/2 0 Table 4.1:
Parameters of the tests for weight analysis of joint inversion using correspondence maps.

The RMS plots in Figure 4.12 a) to c) show how similar the results of the first three tests are in terms of iterations needed to converge. Figure 4.12 d) evidences the effect of not giving any weight to the ĝ term, which results in a high RMS for this term.

4.D What if there is no relationship?

The previous tests in this Chapter are based on the hypothesis that a relationship between the property models m 1 and m 2 exists, however we can think what would the behavior of the joint inversion be in the scenario where no relationship is evident between the two parameters. To test this situation we built a resistivity profile using random values. The pair of models from where synthetic responses were computed for this test are observed in Figure 4.14, they describe a random relationship between one another, and we can say that no evident relationship can be recovered from the cross-plot. We carried out two tests of joint inversion using correspondence maps, (i) using equal distribution of weights for each data set i.e. A = [1/3, 1/3, 1/3] (ii) using zero weight for the correspondence maps term i.e A = [1/2, 1/2, 0]

The results of the joint inversion are shown in Figure 4.15. We can see that for test ii)

where no weight is given to the correspondence maps term, the results are similar to those of separate inversions, which confirms the observation from appendix 4.C regarding the impact of the weights in the joint inversion framework. When we have a pair of models that do not show a clear relationship between them, we can circumvent this issue in the i.e in the worst case where no coupling is observed between the parameters, the results will be similar to those of separate inversions.

In summary:

• First joint inversion using correspondence maps between magnetotelluric data and surface-wave dispersion curves.

• The convergence of the joint inversion using correspondence maps is inherently related to the appropriate selection of the polynomial degree in the

ĝ(m 1 , m 2 , â) function.

Bibliography

Abdelfettah, Y., Sailhac, P., Girard, J., Dalmais, E., Maurer, V., and Genter, A. (2019).

Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric. In European Geothermal Congress 2019. European Geothermal Congress.

Aquino, M., Marquis, G., and Vergne, J. ( 2021). Joint inversion of magnetotelluric and In this Chapter we extend the methodology presented for 1D in Chapter 4. We study the sensitivity of surface-wave dispersion curves to shear velocity using a threedimensional modeling code (SW4) from Petersson and Sjögreen [2015] and implement a methodology to recover these dispersion curves. At the end of the Chapter, we study the feasibility of including MT data from 2D models in a joint inversion framework.

On the use of surface-waves for multi-dimensional models

Introduction

The standard approach for recovering shear wave velocity model from surface-wave dispersion curves involves a two-step procedure: (1) the estimation of group/phase velocity maps at different periods using regionalization and (2) the inversion of these maps to obtain 1D shear velocity models.

The first step starts with the computation of dispersion curves using FTAN for all correlations and next regionalization (linear inversion) to obtain dispersion curves in every pixel of a pre-defined grid. After this, the second step is the inversion of the each Chapter 5. Towards a higher-dimensional approach dispersion curve in the grid to recover 1D shear velocity models (non linear inversion), which are usually interpolated to form pseudo-3D velocity models. Gaudot et al. [2020] presented an alternative approach, which consists on considering the whole dispersion diagram as the data space to produce path average shear velocity 1D profiles and then regionalize at each depth to produce a pseudo-3D velocity model. This two-step procedure is considered valid and equivalent to a 3D inversion, if only lateral heterogeneities of long wavelength are investigated [Montagner, 1986]. Such two-step approach is preferable to a direct 3D inversion because the computations are much simpler [START_REF] Yanovskaya | Smoothness criteria in surface wave tomography[END_REF].

In this Chapter, we intend to asses model responses using higher-dimensional models, instead of the standard approach of 1D velocity profiles, for different structures. In Chapter 2 (Section 2.6.2), we described three steps for obtaining surface-wave dispersion curves maps using SW4. We applied the same methodology for the models built in this Chapter. In Section 5.1.2 we explore first a case without lateral changes, i.e. as a 1D model to directly compare the response to one obtained using sdisp96 [Herrmann, 2013].

In Section 5.1.3, we then explore models with varying lateral changes.

One-dimensional model

For validating the methodology introduced in Chapter 2, Section 2.6.2, we first use a three-layer model (Figure 5.1). It has dimensions of x=60 km y=6 km and z=15 km and cells of 80 m, we show here a section in the XZ plane at y=3 km, which is where all sources and receivers were positioned. Note that Figure 5.1 shows the shear velocity model and in order to model the response two other properties are needed i.e. compressional velocity and density. Compressional velocity was computed using a relationship V p = 2.0V s and a constant density of Ψ=2.4 g/cm 3 was assigned.

Synthetic modelling and dispersion curve computation

Synthetic seismogram modeling

To compute the synthetic seismograms, we deployed 21 receivers spaced every 2 km in the x direction at y=3 km and z=0. A record length of 70 s and a sample rate of 15 ms (i.e. 63 Hz) was used. The source was defined as a triangle of 1 Hz. A total of 21 computations were run, i.e. one computation per source position. In this case the positions of the receivers in order to reproduce the effect of having each receiver as a source. Chapter 5. Towards a higher-dimensional approach tions, leading to n(n -1) seismograms with n the number of receivers (i.e. 21 receivers, so 420 time series in this case). For each of them we calculated surface-wave dispersion curves using standard FTAN as described in Section 2.5.1.3 for a period range extending from 0.2 to 10 s.

Figure 5.3 shows dispersion curves obtained for three source-receiver pairs, source at

x s =18 km and receivers at x r =16 km, x r =36 km and x r =56 km. The estimated group velocities for the pairs with the longer offsets (i.e. e) and f)) are closer to those obtained from sdisp96 [Herrmann, 2013]. The pair with the shortest offset i.e. d) however shows a large difference between the group velocities obtained by the two approaches. Further investigation is encouraged to understand the contribution of the several waves that emerge from the signal. At this stage we can conclude that a sufficient offset separation is needed to properly isolate the surface-wave content in the signal and produce accurate surface-wave dispersion curves.

Regionalization

Regionalization is the computation of group or phase velocity maps at a given period from individual measurements between pairs of stations through a linear inversion scheme as described in Chapter 2. First, we define a grid of cells where group velocity curves are to be calculated. For each period, we define the grid with 11 cells spaced every 4 km i.e. twice the receiver spacing, marked as the cyan crosses in Fig 5 .1. From the previous step we have computed 420 estimates of dispersion curves from 420 Green's functions. The number of selected measurements depends on the period. Note that each stations pairs are used twice i.e. r 1 -r 2 (source at r 1 ) and r 2 -r 1 (source at r 2 ).

The number of selected measurements depends on the period. In Figure 5.4 we see the pairs of stations used for surface-wave tomography at selected periods. These are the pairs that satisfy the condition where the offset is greater than 2λ (wavelengths). The contributing pair distribution is relatively uniform until 2 s, i.e. wavelengths of about 5-6km. Above this period some pairs of stations are discarded, up to >8 s where none of the pairs satisfy the condition (2λ), and therefore data for depth inversion is not used above this period. A linear inversion is performed using these station pairs and a group velocity map of the profile is generated as seen in Figure 5.5 b). We can compare in Figure 5.5 a) the group velocity obtained using sdisp96 code (assuming a locally homogeneous model) versus that obtained from regionalization of SW4 data 5.5 b). The dispersion curves recovered from our workflow using synthetic seismogram from a 3D elastic code is in good agreement with that obtained from using the standard 1D modeling codes. This is expected and acts as a quality control that the methodology produces similar results as a 1D code since the three-layer model is indeed of one-dimensional behavior. This data will then be used for shear velocity inversion.

Depth inversion

Converting group dispersion curves into a depth model requires solving a non-linear inverse problem that has generally multiple solutions. Here we invert group velocities obtained from regionalization (i.e Figure 5.5 b)) at each of the 11 locations of the predefined grid using the Occam approach [Constable et al., 1987] described in the previous

Chapter. The density is assumed to be constant (2.4 g/cm 3 ) and the V p /V s ratio is kept constant at 2.0.

The inversion performed here is using the standard approach where we determine the 1D shear velocities at each point of the area under investigation from corresponding dispersion curves. Inversion starts with an homogeneous initial model with V s =2 km/s. We constrain shear velocity with lower (2.8 km/s) and upper (3.2 km/s) bounds after 5 km depth, where almost no influence is observed from the sensitivity kernels (not displayed here). The inverted models at each position of the grid are very well recovered as can be seen in Figure 5.6. Inversions converge at the first two iterations and the interfaces at depth are easily recovered. Note that this is a very simple model and was mostly used to corroborate the methodology for obtaining dispersion curves using a three-dimensional elastic forward code.

Two-dimensional model

The previous Section exemplifies the methodology using a model with no lateral changes mainly to make direct comparisons with responses obtained from traditional 1D approach as sdisp96 [Herrmann, 2013]. In this Section, we define four models with lateral changes of various dimensions, processed using the same methodology as in the previous Section.

The four models are labeled as i) one block 35 km, ii) one block 20 km wide, iii) one block 10 km wide and iv) one block 5 km wide and are shown in Figure 5.7. They all start with a two-layer background from where the three-layer model of Figure 5.1 is inserted, all with different lateral extents, in order to analyze the impact of anomaly size on sensitivity.

5.1. On the use of surface-waves for multi-dimensional models Chapter 5. Towards a higher-dimensional approach

Computations for these four models use the same number of receivers as the three-layer model of Figure 5.1 (i.e. 21 RX spaced every 2 km).

Synthetic data processing

The four models were processed following the steps described previously. Each model was computed 21 times (one per source location) producing 420 dispersion curves for all station pairs. After this, group velocity maps were generated and can be seen in Figure 5.9. As it was done for the three-layer model before, we computed the group velocities using 1D models at each position of the grid cells using sdisp96 [Herrmann, 2013]. The equivalent response using 1D models for the four models tested is shown in Figure 5.8.

Since the geometry of sources and receivers is the same for all models, the lack of station pairs to contribute to long periods is the same for all cases as in Figure 5.4. A way to overcome this would have been the deployment of even longer offsets (longest is 50 km and shortest 2 km). Data from models i) and ii) (Figure 5.9) seem to recover the main features as the responses using 1D models (Figure 5.8). Data from models iii) and iv) are harder to analyse. Data from model iii) appears to have a lower velocity than 1D responses for long periods (>6 s), this could be simply artifacts from the small number of station pairs contributing at these periods or a difficulty to separate the wave content at these anomaly sizes (i.e. <10 km wide). Clearly, we need to investigate in more detail the wave contents in the signal for models iii) and iv) and how it affects the computation of the dispersion curves. These topics are out of the scope of my Thesis and for now we will focus on the depth inversion part, but we encourage further investigation on this.

Depth inversion

The purpose of the Section is to analyse the sensitivity of surface-wave dispersion curves obtained from synthetic seismograms of 2D models. Before performing the depth inversions we computed the sensitivity kernels using disba module introduced in Section 2.6.1.

To compute the kernels we took two different 1D models, one at x=20 km and the second at x=28 km from models in Figure 5.7. The 1D models at these two positions are the same in all four models of Figure 5.7.

In Figure 5.10 model at x=20 km is labeled two layers, and model at x=28 km is labeled three layers. Figure 5.10 shows the kernels for group velocity relative to shear velocity for the fundamental model of Rayleigh wave. They represent the depth that determines the group velocity for a given period. Normally, longer periods are greater investigation depths. Periods <1.2 s are sensitive to the first kilometer. 2 s period has a slight wider sensitivity below 1 km up to 1.5 km. Wider peaks are distinguished at periods between 2 s and 7 s for depths from 1 to 4 km.

1D inversion

A total of 11 one-dimensional inversions were performed and results are displayed in Figure 5.11. Bounds were set for shear velocity from 1 km/s to 3.2 km/s, nevertheless 5.1. On the use of surface-waves for multi-dimensional models Chapter 5. Towards a higher-dimensional approach below ≈ 4 km narrower bounds were set (2.8 km/s to 3.2 km/s) due to the minimal sensitivity observed at these depths. Homogeneous half-space models were used as starting models for inversion.

The inversion for the first model i) one block 35 km, shows the main features recovered i.e. lateral change from two to three layers starting at x=25 km. In the second model ii) one block 20 km, the three-layer block is well resolved laterally. The third and fourth models (one block 10km and one block 5 km) are the poorest resolved. It is expected since at the data processing part, the group velocity maps did not seem to recover the lateral changes (Figure 5.9). From this short exercise we can say that lateral changes >20 km are able to be recovered unlike to the smaller size anomalies i.e. models iii) and iv) with the receiver spacing designed here (i.e. 2 km) and a maximum offset separation of 40 km.

Pseudo-2D inversion approach

Synthetic model responses here have been computed using a three-dimensional elastic code. This approach is feasible for obtaining observed data for sensitivity analysis, nevertheless the implementation of this approach in an inversion framework, where forward 5.1. On the use of surface-waves for multi-dimensional models responses are computed at each iteration, becomes highly expensive in computation time.

Therefore, for the inversions presented previously we use 1D forward responses.

Up to now, shear velocity inversions have been done locally (at each grid position) assuming local 1D models. In order to compute 2D joint inversions with MT data we need first to find a way to build practical 2D velocity models i.e. 2D models that share a common mesh. The method we describe below is only applied for surface waves.

We use here a similar approach as Ogunbo [2019] for building a pseudo-2D model re- The Jacobian matrix J 2d is full of zeros in the entries where no data sites (nsite) are available.

Figure 5.12: Schema explaining the data and Jacobian matrix for both one-dimensional (left) and pseudo two-dimensional modeling for SWD data.

Pseudo-2D inversion application

In order to show if reliable velocity models can be obtained from pseudo-2D inversions we applied the approach to the four models previously introduced. We built a 2D mesh for each of the four models in Figure 5.7 with a total of 5588 parameter cells. Using the data recovered in Figure 5.9 through regionalization we computed pseudo-2D inversions using the standard Occam approach that has already been explained, with the only variation of using 2D smoothing matrices (i.e. to smooth the vertical and horizontal directions). More details in the construction of the 2D smoothing matrices can be found in deGroot-Hedlin and Constable [1990].

Multi-physics approach

The results obtained from pseudo-2D inversions are displayed in Figure 5.13.A 15 % standard deviation was used for all three data terms as for 1D inversions. Similar observations as from the 1D inversions can be made from the pseudo-2D inversion models.

The inversion models for the first two cases: i) one block 35 km and ii) one block 20 km, can identify the lateral changes from the true models. It is not the case for iii) one block 10 km and iv) one block 5 km. In all four cases, it is expected to obtain smoother models, laterally and vertically, due to the nature of the application of the smoothing matrix in 2D. Note the small traces (column wise) of the positions were data is present, this is an effect of the way the Jacobian matrix was built with only values at those locations.

The purpose of this exercise was to show whether inversion could perform stably using this pseudo-2D approach. The method has shown to provide reasonable pseudo-2D velocity models. Now that the approach has been validated, we will use it in a 2D joint inversion framework with MT data from 2D resistivity models. This subject is discussed in the next Section.

Multi-physics approach

From the tests in the previous Sections, we have observed limitations from the estimated group velocity maps and translated to the inverted models obtained in both 1D and pseudo-2D inversions, specially the shear velocity models with the smallest anomalies.

Here we intend to study the sensitivity of MT data for resistivity models with similar structure geometries as those defined for shear velocity. We use only the one block 10km model since this is one where limitations on lateral sensitivity were observed for seismic data. This test will assess whether the inclusion of MT data in the joint inversion brings improvements.

We have described previously in Section 4.2 the basics of the joint inversion approach using CM in 1D. The same methodology can be applied for 2D models with few modifications as pointed out for the pseudo-2D inversion. The 2D smoothing matrices are built to smooth vertically and horizontally rather than only vertically as for 1D. In the next Sections, an expanded version of the already covered joint inversion using CM in 1D is applied to 2D models.

Resistivity model case 1: degree 2 relationship

We built a resistivity model using the same mesh used previously for the pseudo-2D modeling of velocity, that meant 5588 parameter cells of resistivity values. A similar structure as that of the shear velocity model was used to build the resistivity model as can be seen in Figure 5. Apparent resistivity curves were computed using the PW2D code of Wannamaker et al.

[1987] for a range of frequencies from 0.01 to 100 Hz. In the MT response from this Chapter 5. Towards a higher-dimensional approach model (Figure 5.15), the conductive anomaly is clearly detected, both at the start of the lateral change (station 5) x=25 km and at the end of it (station 10), x=35 km. It is interesting to look at the sensitivity in both modes (TE and TM). For the TM mode the transition is quite sharp, since the electric field is oriented across the strike of the conductive structure. On the contrary, the lateral sensitivity of the TE mode is more smooth since the magnetic field is oriented across the conductive structure. This is a clear example of how the TM mode is more effective than the TE mode at locating lateral interfaces between bodies of varying resistivity. In this case, x direction is the direction of the profile and z direction along the depth.

Multi-physics approach

joint inversion results dot not greatly benefit from the CM coupling term.

Resistivity model case 2: degree 1 relationship

Once we understood the limitations of building a relationship of a parabola type with only three points of model pairs, we attempted to build a simpler relationship to test if the joint inversion algorithm would perform stably and the anomaly of the velocity model could be better resolved.

A new resistivity model was built and is shown in Figure 5.17. The parameter relationship type in this case is a straight line type i.e. degree 1. A linear relationship of this type is at first approach simpler to fit even on the presence of only three model pairs values, since the space of acceptable parameter relationships is better resolved compared to a degree 2. We performed a joint inversion using CM with a degree 1 for the ĝ(m 1 , m 2 , â)

function. For convergence, the initial model used for joint inversion, was the second iteration of the resistivity model from separate inversion. The results in Figure 5.18 show that the velocity anomaly is importantly upgraded compared to separate inversions tests.

The correspondence maps plot closer to a straight line as it fits with the true relationship (black squares in the cross-plots).

Figure 5.17: Section 5.2.2 resistivity model (right) built using geometry from shear velocity iii) One block 10 km (left) and the corresponding parameter relationship plot below

Conclusion

Conclusion

In this Chapter we have implemented a framework for computing surface waves dispersion curves from synthetic seismograms computed with a 3D elastic code. Synthetic data shows limitations in the lateral resolution of surface waves for small anomalies. But show that a pseudo-2D inversion is feasible for obtaining 2D shear velocity models if the velocity anomalies are sufficiently large. We expect that the higher sensitivity of MT data to lateral contrasts will mitigate this problem in a joint inversion.

Attempts to perform joint inversion with synthetic data in 2D have shown several limitations when the relationship type is of a too-high degree. However we are able to improve from separate inversion models when using a degree 1 polynomial, even in the presence of only three parameter pairs. This is because a degree 1 type of relationship has fewer acceptable coefficients sets.

Complex 3D synthetic models are encouraged to be built for examples and apply a similar workflow for estimating group dispersion curves for assessing the joint inversion approach.

In summary:

• Dispersion curves obtained from 2D models are different than those from standard approaches in 1D.

• Limitations on lateral resolution of surface-wave dispersion curves are observed where MT data can be beneficial in a joint inversion framework for higher dimensional models.

• Joint inversion using CM is feasible when the dynamics of the model pairs are sufficient to define a a limited number of parameter relationships. Chapter 4 and Chapter 5 described the implementation of joint inversion with CM using mostly synthetic models. In this Chapter we compare subsurface images obtained from separate and joint inversions of field data in North Alsace. Group velocity curves processed by Lehujeur et al. [2018] and MT data processed by [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF] are used for inversion, along with new MT data acquired during this Thesis.

Bibliography

General overview

North Alsace region in France has extensively been studied for its natural resources, historically for oil and gas exploitation between 1740-1970s [START_REF] Préveraud | L'alsace, pionnière du pétrole[END_REF] [START_REF] Ch | Deep geothermal energy in western europe: the soultz-sous-forêts project[END_REF]. The EGS consists on extracting heat from hot water produced into one or more wells by injecting cool water into another well. This concept is efficient even when permeability is naturally low but can be enhanced by means of hydraulic and/or chemical simulations [START_REF] Lu | A global review of enhanced geothermal system (egs)[END_REF].

The successful development of deep geothermal systems relies on knowledge of key parameters: i) the relationship between depth and temperature, ii) deep-seated geology, iii) fracture distribution and iv) stress regime, since the maximum principal stress direction is essential for fluid circulation [START_REF] Ch | Deep geothermal energy in western europe: the soultz-sous-forêts project[END_REF]. Within the definition of the geology, a key element is the definition of the top and base of lithological units, specifically the reservoir formation. In Soultz-sous-Forêts this is the granitic basement, with its top between 3 km to 5 km the deepest in the middle of the graben.

Multiple geophysical studies have been performed to map the lithological structures of the area: ASNC [Lehujeur, 2015], gravity [START_REF] Portier | New results on the gravity monitoring (2014-2017) of soultz-sous-forêts and rittershoffen geothermal sites (france)[END_REF], MT [Geiermann and Schill, 2010a], [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF], MT and gravity [START_REF] Abdelfettah | Using highly accurate land gravity and 3d geologic modeling to discriminate potential geothermal areas: Application to the upper rhine graben, france[END_REF], [START_REF] Geiermann | 2-d magnetotellurics at the geothermal site at soultz-sous-forêts: Resistivity distribution to about 3000m depth[END_REF]. Lehujeur et al. [2018], based on dispersion curves recovered from ASNC, produced a pseudo-3D shear velocity model of the region from 1D inversions, and stressed the lack of correlation between the deeper sections of their model (approximately >1.5km) and known lithology. [START_REF] Schill | 2-d magnetotellurics and gravity at the geothermal site at soultz-sous-forêts[END_REF] produced 2D density and resistivity profiles from gravity and MT data separately and distinguished the top of the granitic basement. Studies in the area have shown high salinity levels in the brine [START_REF] Pauwells | Water-rock interactions during experiments within the geothermal Hot Dry Rock borehole GPK1, Soultz-sous-Fore^ts[END_REF] resulting in high conductivities (10 to 100 S/m) at a reservoir temperature of 200 • C [Geiermann and Schill, 2010a]. The nature of the MT method at detecting conductive structures makes it an excellent imaging tool for geothermal reservoirs which we intend to exploit along surface waves dispersion curves.

To the best of our knowledge, there is no study that has attempted to integrate MT data and surface-wave dispersion curves from ASNC to produce subsurface images with a joint inversion approach in the area.

Here we perform one and two-dimensional joint inversions using the same objective function as described in Chapter 4. Two profiles were studied in the area (Figure 6.1), using data recovered from multiple surveys. We first give a brief geological background of the area, next we cover the details of the data sets used and then we describe the main results from separate and joint inversions for the two studied profiles. [Lehujeur, 2015], orange square are the MT stations from ECOGI survey [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF], yellow squares are MT stations from the DEEP-EM survey and purple points are the MT-Two extended MT stations. These two last datasets where acquired during the course of this PhD in 2020 and 2022 respectively. Background layer is the geological map obtained from BRGM Figure 6.3: Rayleigh-wave group velocity dispersion maps from EstOF at selected periods for the fundamental mode. The amplitudes correspond to velocity variations relative to a reference group velocity specified in each subplot. This data was processed by Lehujeur [2015] 6.2. The datasets

MT data

The In house: Razorback Table 6.1: Acquisition parameters for the three MT surveys discussed here. Blue fill shows data used in Profile 1 and green fill the data used in Profile 2.

MT data: ECOGI

The MT data used in Profile 1 (NW-SE, Section 6.3) is composed of 10 stations that were acquired in 2018 as part of the ECOGI project and were processed by [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF]. The data was acquired for three days at a sample rate of 512 Hz and processed using one magnetic remote reference station located 70 km to the south west of the survey and operated by IPGS (now ITES) at the Welschbruch observatory. The data recovered from [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF] are displayed as pseudosections in Figure 6.6.

MT data: DEEP-EM

The MT data used in Profile 2 (W-E, Section 6.4) is recovered from two surveys that were acquired during the course of this Thesis. As mentioned previously these are DEEP- Here we stress that considering the restricted frequency range used in MT data (up to 10 s) for the 1D inversions (Figure 6.9), it is expected that the base of the conducting sediments will be difficult to image, as the sensitivity of the MT data in that frequency range is concentrated the the shallower part of the model. Several tests including a broader range of frequencies confirmed this. In Figure 6.10 we show a comparison of two resistivity models to illustrate this. In the model obtained using a restricted frequency band (left) the base of the sediments is not evident, while in the model recovered using the full frequency band (right) the limit of the base of the sediments is noticeable due to a better sampling of the resistivity basement. 2D inversions are encouraged using the full frequency range for the profile and are performed in the next Section 6.3.2. For the two-dimensional models a common mesh was designed formed by irregular cells, with finer cells in the area of interest. A total of 810 parameter cells were set for each model. Ten MT stations and 16 dispersion curves locations were used. The full frequency band for MT was used i.e. 0.01 to 64 s while we kept the same for dispersion curves (0.6 to 4.16 s). The starting models for both joint and separate inversions were half space models. Bounds were set for velocity after 2.8 km to 3.2 to 3.6 km/s as suggested by [Lehujeur et al., 2018] from reference well-log velocities and due to the low sensitivity of the data at these depths.

Profile 1: North west -South east

A polynomial relationship of degree 1 was discarded from the tests performed in 1D which suggested a higher-order of relationship. A polynomial relationship of degree 2 limited to 3 coefficients âij was used for the ĝ function. Equal weights were assigned for each of the three terms i.e. 1/3 for SWD, 1/3 for MT and 1/3 for the CM.

Figure 6.11 shows the resistivity and shear velocity models obtained from inversions.

The resulting relationship is given by m 2 = 0.61m 2 1 -2.43m 1 + 0.5, where m 2 is the log10 of electrical resistivity and m 1 is log10 of shear velocity. The cross-plots obtained from the joint inversion models are more concentrated than those of separate inversions. This is the effect of the CM in the joint inversion, trying to satisfy a relationship between the model pairs while fitting the data. The joint inversion velocity map is smoother below 1km compared to the separate inverted velocity model. This is the depth at which kernels from Figure 6.7 suggested a lower sensitivity. [START_REF] Maurer | On-going seismic monitoring of the rittershoffen egs project (alsace, france)[END_REF] and [START_REF] Aichholzer | New stratigraphic interpretation of the Soultz-sous-Forêts 30-year-old geothermal wells calibrated on the recent one from Rittershoffen (Upper Rhine Graben, France)[END_REF] identify the transition to the Paleozoic granitic basement at 2250 m in well GTR-1 located approximately at x=14 km of profile. This is in good agreement with the resistivity model that defines the boundary between the sediments and a more resistive granitic basement (Figure 6.12). This confirms the sensitivity of the MT method to thick, high-resistivity layers as is the case here. In the shear velocity map the contrast between sediments and basement is not very strong and the basement definition might not be straightforward. The basement depth from the joint inversion however is more realistic to that obtained from separate inversion. These observations hint that approaching the joint inverse problem from a petrophysical point of view can bring structural similarities in the subsurface images. This question is expanded in more detail later in this Chapter. Figure 6.12: 2D joint inversion results for Profile 1 in Section 6.3. Top of Basement horizon (black pointed line) is my own interpretation in the resistivity model and overlayed in the velocity image.

Weight analysis in the joint inversion using CM

CM introduces an additional term in the joint inversion objective function that already contains MT and SWD data fit terms -Occam finds the optimal regularisation weight (λ in equation 4.1) regardless of the other term weights. How can one determine the proper weights to each term and how can one assess their impact on the inversion results?

As a test, we introduce individual weights to compensate for the diverse amount of data from each of the three datasets, or simply to increase or decrease the influence of one term with respect to the others. Here we compare two combinations of weights given to each of the three objective function terms in the joint inversion (A D in equation 4.2).

1. Case 1 is the one previously described for where all terms were given equal weight i.e. 1/3, 1/3, 1/3 each.

2. Case 2, we increase the CM term i.e. the coupling ĝ(m 1 , m 2 , â) function, by assigning a weight of 1/2 to it and 1/4 to the other two terms.

The effect of using a larger weight to the CM term results in a sharper parameter relationship in the cross-plot (green in Figure 6.13). Note that in both joint inversion tests, the data fit is not compromised. Case 1 shows a global RMS of 1.91 compared to 1.89 for case 2. Moreover, the higher the weight is given to the CM term, the closer the ĝ function to -1 (observed data for CM), as would be expected.

Another way to compare the models is by looking at the differences between them.

For this, we computed the difference for each of the two joint inversion attempts models minus the separate inversion models. Here, a positive delta relates to an increase in either velocity or resistivity in the joint inversion models compared to the separate inversion models. In general we can say that both joint inversion models are very similar since the differences obtained are almost the same for both velocity and resistivity models (Figure 6.14).

Profile 2: West -East

The profile 2 is located in the northern part of the EstOF network as can be seen in Figure 6.1. It runs E-W with a length of ≈20 km. Group velocity maps from Figure 6.3 highlight a high velocity structure towards the northwestern part of the region. The presence of a major geological discontinuity that separates the northern Vosges massif from the sedimentary plain by a northeast-southwest running fault (Rhenish fault) is well documented [START_REF] Duringer | The complete lithostratigraphic section of the geothermal wells in rittershoffen (upper rhine graben, eastern france): a key for future geothermal wells[END_REF], [START_REF] Bauer | Architecture, fracture system, mechanical properties and permeability structure of a fault zone in lower triassic sandstone, upper rhine graben[END_REF]. The intention of this second profile is to compare this lateral contrast by a multi-physics approach. We aim to study if the velocity contrast observed in the dispersion curves is also detected with MT data.

One-dimensional inversions

Similarly to the approach in the previous profile, we first computed inversions in 1D (Figure 6.16). A total of 12 pairs of stations were inverted using a halfspace as an initial Compared to the previous profile, correspondence maps in this region (Figure 6.15)

do not show a concentrated pattern over all sites: on the contrary, the cross-plot shows a sparse behavior. This is a first sign that the petrophysical relationship might not be common along the whole profile. To explore further the origin of the behavior in the cross-plots, we split the profile in three zones, i) west side, ii) center and iii) east side.

The west side includes the first 5 stations, the center the next 2, and the east side the remaining 5. The cross-plots with each zone colored distinctly are displayed in Figure 6.17. They clearly show how the relationship between the two physical properties varies along the profile, with western stations having the highest resistivity and velocity values and the eastern stations the lowest. Although only part of the MT frequency band was used due to our restriction to one-dimensionality, they suggest a petrophysical variability along the profile.

Two-dimensional inversions

As has been suggested by the 1D inversion results above, the petrophysical behavior along this profile is more complex than for Profile 1. From a structural point of view, Profile 2 crosses the Rhenan fault that results in a major geological discontinuity between the northern Vosges massif and the sedimentary plain. Before performing joint inversions with our own codes, we analyzed the MT data separately with the aid of inversion code MARE2DEM [Key, 2016].

Separate inversion of resistivity using MARE2EDEM

MARE2DEM is a freely available code for 2D modeling and inverting MT and CSEM data [Key, 2016]. It uses an unstructured triangular mesh. The non-linear inversion is carried out using the Occam approach described before.

The codes developed for 2D inversions in this Thesis have the limitation that they do not take into account the topography, all receivers are placed at z=0 and models start from a flat surface, which is usually not the case. Profile 2 runs along a zone where topography changes significantly compared to profile 1, ≈ 200 m between the lowest and highest sites (Figure 6.18).

Before performing inversions with MARE2DEM, we tested the effect of topography sampling on the modelling results. Our test model was a homogeneous subsurface (1 Ωm) and topography sampled at 100 m, 500 m and 1000 m. We then compared the computed forward responses for both TE and TM modes to their expected values. We found that a topographic sampling of 1000 m produced the least artifacts in the responses. The largest artefacts in the topography sampled at 100 m had a 5 to 10 percent error compared to the expected response i.e. 1 Ωm.

After these observations we designed a mesh for inversion using the 1000 m topo- A second inversion was made in order to compare the model obtained previously with one without topography. A similar mesh was generated and a 1 Ωm initial model was used. The inverted resistivity model without topography is displayed in Figure 6.20 (below). The results show that the main structures are preserved and the differences between with and without topography models are very small.

Joint inversion

After the understanding of the subsurface based on the resistivity image, we designed a mesh to jointly invert the MT data used previously with surface-wave dispersion curves.

A finer mesh in the shallow part was built compared to that used for Profile 1, due to a greater lateral variability on this profile. Each model had a total of 1350 parameter cells.

12 stations were used for MT data plus 16 locations of dispersion curves that fall along the profile. As done in the inversions in 2D for Profile 1, a half space was used as initial Figure 6.20: 2D resistivity model inverted using MARE2DEM with topography (above) and without topography (below). Overlayed in black line is my own interpretation of the Top of the Basement.

6.5. Discussion and conclusions model for both separate inversions. A degree 2 was chosen for the polynomial relationship of the ĝ function. Results in Figure 6.21 show the difficulty for the CM to find a common relationship along the whole profile. This time, cross-plots are color-coded according to their site position along the profile, rather than depth as in Figure 6.11. We can correlate them to the observations made from the 1D inversions of Figure 6.17. The west side of the profile shows a shallower transition to resistive layers due to the uplift of the basement on that side of the fault, which correlates with the high-velocity anomaly. The east side of the profile shows a thicker low-resistivity package (sediments). This is expected as we know from structural observations that the oldest geological units outcrop in the west side.

The joint inversion using CM has difficulties to find a common relationship along the whole profile. An option could be to compartmentalize the profile into distinct zones, each with its own CM. This would require re-organizing the bottom part of the Jacobian matrix (eq. 4.7).

Discussion and conclusions

In this Chapter we inverted data for two profiles in the North of Alsace located in two structurally different zones of the Upper Rhine Graben. We dedicate this Section to discuss some observations from the models obtained and correlate with some geological and petrophysical information.

Profile 1: NW-SE

The results of this profile, both 1D and 2D inversions, show a singular parameter relationship throughout the whole profile. This is expected as the profile is located in the basin and no major discontinuities are present. Resistivity contrast between the sediments and the granite expose the limits of the top of the basement that correlate with known lithology from wells drilled close to the profile (GRT-1). For well GRT-1 the top of the basement was reached at 2212 m (MD) [START_REF] Aichholzer | New stratigraphic interpretation of the Soultz-sous-Forêts 30-year-old geothermal wells calibrated on the recent one from Rittershoffen (Upper Rhine Graben, France)[END_REF], [START_REF] Duringer | The complete lithostratigraphic section of the geothermal wells in rittershoffen (upper rhine graben, eastern france): a key for future geothermal wells[END_REF].

The shear wave velocity obtained at greater depths from the separate inversions is not straightforward comparable with the known lithology. It has been suggested as a cause of this poor vertical resolution, the low sensitivity of surface waves to deep layers (> 1.5 km) as observed in the kernels (Figure 6.7). Another reason, could be the addition of higher frequencies to better resolve the shallow part and inherently improve the deeper parts of the image [Lehujeur et al., 2018]. However, the sensitivity of the MT data to thick resistivity layers helps alleviate this issue in the velocity image without compromising the data fit. The shear velocity model obtained from the joint inversion is more coherent with the known lithology, especially below 1.5 km.

We correlated the values of the models obtained here with proxies from laboratory measurements on core samples by [START_REF] Kushnir | Characterizing the physical properties of rocks from the Paleozoic to Permo-Triassic transition in the Upper Rhine Graben[END_REF], who measured thermal conductivity as a function of several properties, including compressional velocity. Keeping in mind all the caveats and limitations of this approach, we use thermal conductivity as a Figure 6.21: 2D inversion results for Profile 2 discussed in Section 6.4. Top: shear velocity models, middle: resistivity models, lower: correspondence maps. Left: Initial models, middle: separate inversions, right: joint inversions.

6.5. Discussion and conclusions proxy for electrical conductivity and compressional velocity as a proxy for shear velocity.

The former is somehow more speculative but the latter goes without saying, although their proportionality (i.e. V p /V s ) can be affected by the porosity and the nature of the pore fluid. Figure 6.22 shows a comparison between the results of the joint inversion models with observed physical properties of samples in the zone. In general, properties in the sediments (red shading) vary in the same direction: the higher the velocity, the higher the thermal conductivity. However, the properties from samples in the basement (blue shading) vary in opposite directions: the higher the velocity, the lower the thermal conductivity. The patterns of the pair of subsurface models obtained from joint inversion (Figure 6.22, left) are consistent with their proxies measured in the laboratory (Figure 6.22, right).

Comparison with other joint inversion operators

During this Thesis we focused on the application of correspondence maps only, i.e. a petrophysical type of joint inversion. However, we were interested to compare joint inversion results using structural and petrophysical couplings. In this Section, we attempt to draw a comparison between these two approaches. Implementing structural joint inversion is beyond the scope of this Thesis, but for sake of comparison, we computed three types of operators often included in the objective functions: i) the cross-gradients (CG), ii) dot product-image gradients (DPG) and iii) the variation information (VI). The equations for each of these operators have been defined in Chapter 3.

The values of these operators were computed using the results of the models obtained from separate and joint inversions and compared ( Figures 6.23 and 6.24). They evidence that the three other operators of joint inversion are also greatly decreased. The cross- In Chapter 1, I listed the main questions I wanted to address during the course of my Doctoral studies:

1. Can we implement a joint inversion approach using ambient seismic noise and magnetotelluric data in one and higher-dimension models?

2. Can we obtain meaningful velocity-resistivity relationships from geophysical data?

3. How similar or different models from structural and parameter-relationship joint inversion approaches can be?

At this stage, after three years of research, I report on contributions regarding these questions, and propose future directions of research.

Conclusions Items 1-2: One-dimensional joint inversion using CM

The work in Chapter 4 presented the framework of the joint inversion method applied throughout this Thesis i.e. the correspondence maps (CM). We applied successfully a first time correspondence maps joint inversion in 1D between surface-wave dispersion curves and magnetotelluric data.

The first tests using synthetic data from known parameter relationships showed that we were able to recover these relationships by including CM in the joint inversion. The model space of resistivity and shear velocity is hence reduced, which contributes to the mitigation of the non-uniqueness in the inverse problem.

Conclusions

The time consumption of the forward modeling for any 3D mesh and the ensuing processing, of velocity/density models, limited our ability to run a lot of tests. However, the construction of a higher degree relationship using a model with greater dynamics (i.e. more layers and blocks of varying velocity and resistivity) would be an important step to assess the feasibility of the joint inversion at higher dimensions.

This work was part of an oral presentation [Aquino et al., 2022a].

Items 1-2-3: Application: North Alsace

The work in Chapter 6 focuses on the application of the methodology to data sets from North Alsace. The area has been subject to multiple exploration campaigns for natural resources. Several studies performed separately e.g. passive seismic, MT, gravity among some, show the potential of an integrated multi-physics image. However up to this Thesis, no work has been published that performs joint inversion of multiple geophysical methods in the area. We show it is possible to perform successfully integrated images between resistivity and velocity.

For Profile 1 (NW-SE, all well inside the basin), we observe an overall 1D structure with small lateral variations. Joint inversion improves the velocity model in the deeper section compared to the separate inversion. This can be explained by the limited period range of SWD and hence the overall better sensitivity of MT at greater depths.

Furthermore, we show that, at least for this example, joint inversion models using a petrophysical coupling (i.e. CM), can produce images structurally more similar than those from separate inversions. This petrophysical vs structural coupling relationship has been a subject of research over the years of development of multi-physics inversion methods [Gallardo and Meju, 2004], [Moorkamp et al., 2007], [START_REF] Haber | Model Fusion and Joint Inversion[END_REF], [Sun and Li, 2016], [START_REF] Colombo | High-resolution velocity modeling by seismic-airborne tem joint inversion: A new perspective for near-surface characterization[END_REF].

The resulting electrical resistivity and shear velocity models are consistent with trends observed on proxies from laboratory measurements (thermal conductivity and compressional velocities) for both sediment infill and basement. In addition, the resistivities in the sedimentary units are comparable to those observed on well-logs.

For Profile 2 (W-E near the edge of the URG), our MT acquisition includes sites outside of the North Alsace basin. A major result is the evidence of a resistivity increase coincident with a high-velocity anomaly found earlier in group dispersion maps. We are able to correlate the known location of the basement on the resistivity model. We observe a stark contrast in petrophysical properties (and hence geological structures) between the west and east parts of the profile, likely separated by the Rhenan fault.

Extra item: Alternate approaches to 1D SWD inversion

After finalizing the results of my work on 1D joint inversion using CM, I wanted to push farther the work and apply my Occam-based inversion to the whole EstOF data set to produce volumes for comparison with those produced by Lehujeur et al. [2018], who kindly provided them in usable format. I also took an interest in model-space sampling methods and so implemented a fast MCMC inversion that I applied to the same data set.

The comparison is presented in Appendix B.

I also looked at implementing L-BFGS optimization for separate and joint 1D inversions of SWD and MT data. My results were mixed and are presented in Appendix C.

Perspectives

I identified throughout the Thesis future directions on research that are worth to recapitulate in this section.

3D modeling

Generally, all components of the MT impedance tensor are non-zero and require a threedimensional modelling approach to be reproduced. We have limited in this Thesis to 2D models obtained from PW2D [Wannamaker et al., 1987]. The resistivity models of northern Alsace would be substantially improved by the inclusion of other stations that were acquired that fall outside the two profiles presented in this work in Chapter 6. It is therefore highly encouraged to expand the approach to 3D modeling of MT data, say by using the code of [START_REF] Mackie | Three-dimensional magnetotelluric modelling using difference equations -theory and comparisons to integral equation solutions[END_REF] or more recent ones like ModEM by [START_REF] Kelbert | Modem: A modular system for inversion of electromagnetic geophysical data[END_REF].

For the case of surface-wave data, even though we explored the sensitivity using a 3D modeling code, all of the models tested actually had a 2D geometry. Therefore it would be interesting to extend the synthetic approach for sensitivity analysis of real 3D

geometries. From the observations we had on the limitations of the lateral resolution of surface waves using the 2D synthetic models in Chapter 5, we think it is important to evaluate the strategy that has been used as standard in surface wave imaging i.e. using 1D modeling for the construction of 3D velocity models. From literature, we know that this approach is valid only when long wavelength lateral changes are studied [Montagner, 1986].

Data processing

The field data examples in Chapter 6 used only apparent resistivity curves of both modes (TE and TM) for the MT data, however phases were not included mainly due to the poor quality of the those obtained after the Razorback processing, especially at low frequencies, a frequent occurrence in North Alsace. A deeper look at the time series is encouraged, say by using an extra remote reference station, perhaps Wittstock (Germany) operated by GFZ Potsdam [START_REF] Ritter | Permanent magnetotelluric reference station in Wittstock[END_REF], to enhance the SNR.

I identified in Chapter 5, for surface-wave data, the necessity to have a closer analysis on the wave content of the time series. This is especially important for the shortest offsets where dispersion curves recovered from the synthetic numerical model responses were substantially different to reference dispersion curves from 1D model responses using sdisp96.

Multi-physics approach

In North Alsace, vast petrophysical information from wells and other geophysical methods have been acquired. Among the other passive geophysical methods, gravity and magnetic data are available. Several recent studies in geophysical exploration Carrillo and Gallardo [2018], Moorkamp [2021b] and Moorkamp [2021a] have shown the improvement of density models from MT data. Ars et al. [2019] implemented a constrained inversion between surface-wave dispersion curves and gravity data using a 3D separate inversion resistivity model as reference.

In addition to the passive methods, active seismic and CSEM (controlled source electromagnetic) data have recently been acquired and, if available, could improve the models obtained here.

One of the main limitations of the correspondence maps approach that I observed through the different tests performed in this Thesis, is the behavior of the coupling term when no obvious parameter relationship was observed between the model pairs or when multiple parameter relationships describe the model pairs. This situation is difficult to handle, one way is to lower the contribution of the CM term, but it inherently means less coupling between the two datasets.

Another option could be the implementation of multiple parameter relationships, as already suggested in Section 6.4.2.2 by partitioning the CMs. An interesting approach proposed by Sun and Li [2016] inputs a predefined number of clusters (C-means) into the joint inversion. The CM term has the advantage that no prior knowledge about the number of clusters is required but perhaps providing some initial guidance might be a step in the right direction. Recently, as mentioned earlier in this Thesis, Moorkamp [2021b] produced stable inversions using a new coupling term called VI. It appears to be more flexible in the cases for which no obvious relationship is present, which could be a future direction of research to try, especially in areas of complex structures, e.g. 3D models

where multiple parameter relationships can coexist.

Bibliography

Aquino, M., Marquis, G., and Vergne, J. (2022a). Joint inversion of magnetotelluric and ambient seismic noise data using correspondence maps. In EGU General Assembly 2022.

EGU22-5056. Aquino, M., Marquis, G., and Vergne, J. (2022b). Joint one-dimensional inversion of magnetotelluric data and surface-wave dispersion curves using correspondence maps.

Geophysical Prospecting.

Ars, J.-M., Tarits, P., Hautot, S., Bellanger, M., Coutant, O., and Maia, M. (2019). Joint

A.2.1 Razorback output

Below are the six stations processed with Razorback using one reference station i.e.

Welschbruch. The stations used in Profile 2 from the MT-Two survey are all 6 acquired stations i.e. sites 401 to 406

B.5 Discussion

As mentioned at the beginning, we intend to make a comparison between the obtained models here and a reference one i.e. [Lehujeur et al., 2018]. We remind that since the number of layers used in the two inversion models presented here and the one from Lehujeur et al. [2018] (19 layers compared to 100 layers respectively), then an exact comparison between depths is not entirely possible. However, each depth slice can already give a good first approach to make observations.

I N T RO D U C T I O N

Joint inversion of geophysical data has been proposed as a means for understanding subsurface structures, since taking advantage of different data sensitivities can reduce the range of acceptable models (i.e. the 'non-uniqueness') (Moorkamp et al., 2007). One can implement joint inversion in many different ways: (1) joint inversion for a single physical property from multiple datasets is sensitive to it, for example apparent resistivity from DC resistivity and magnetotelluric data (Jupp

The paper was presented at the 82nd EAGE Conference & Exhibition, Amsterdam, the Netherlands. * E-mail: monicaquino92@gmail.com and Vozoff, 1975) or shear velocity from receiver functions with surface wave dispersion curves (Julià et al., 2000); and (2) joint inversion of data sensitive to different physical parameters, for example resistivity and seismic velocity (Gallardo and Meju, 2003). For the first type, the joint inversion procedure is straightforward as both methods share a common subsurface model. The second type, however, is more challenging since sensing different physical properties involves questions about their mutual coupling (or lack thereof). Here we focus on the second type of joint inversion to obtain shear-wave velocity and electrical resistivity models from surface-wave dispersion curves and magnetotelluric data, respectively. The core element of any joint inversion methodology is, therefore, the hypothesis of a link between the different physical property models (Carrillo and Gallardo, 2018). In general, these links or features can be classified as structural or petrophysical.

Structurally coupled joint inversion assumes that different geophysical methods sense the same underlying geology; therefore, it is the structure that controls the distribution of petrophysical properties (Gallardo and Meju, 2003). With this approach, the structural similarity of the different physical property models is considered rather than petrophysical relationships between these properties. They require the specification of mathematical conditions at structural boundaries in a numerically stable way. Today, the most common structural coupling method is the cross-gradient constraint, introduced by Gallardo andMeju (2003, 2004) and for which the structural similarity is enforced by minimizing the magnitude of the cross-gradient vector. Other structural constraints include either curvature-based methods, for example Zhang and Morgan (1997) and Haber and Oldenburg (1997), or Gramian constraints, for example Zhdanov et al. (2012) and [START_REF] Ogunbo | Gramian constraints in electromagnetic multi-physics joint inversion[END_REF].

Structural coupling does not make any assumption about the petrophysical relationships between the parameters. If there is a priori information about such relationships, joint inversion results can be improved significantly. Hence, structural coupling might overlook the full potential of the joint inversion (Moorkamp, 2017). Joint inversion methods using parameter relationships are based on the fact that for some specific geological environments geophysical parameters can be related to some mathematical relationships. Carcione et al. (2007) showed how in sedimentary environments porosity and fluids of the rock matrix can provide a physical link that causes a correlation between seismic velocities and electrical conductivities. The so-called petrophysical joint inversion has generally been approached in two ways: [START_REF] Introduction Ogunbo | Gramian constraints in electromagnetic multi-physics joint inversion[END_REF] when there is a relationship between the parameters which are used to constrain the dependent parameters, in a sequential inversion mode, for example Dufréchou et al. (2018) and [START_REF] Gautier | Detailed lithospheric structure of an arc-continent collision beneath Taiwan revealed by joint inversion of seismological and gravity data[END_REF] and ( 2) when the parameter relationship is included as part of the joint inversion workflow. From this second group, Tiberi et al. (2008) retrieved linear relationships with the depth between the density and the compressional velocity in a [START_REF] Birch | The velocity of compressional waves in rocks to 10 kilobars: 1[END_REF] type relationship. Carrillo and Gallardo (2018) proposed a way to retrieve higher order polynomial relationships for a number of heterogeneous zones.

We present here a joint inversion methodology that includes the correspondence map approach described by Carrillo and Gallardo (2018), implemented as an Occam-type inversion (Constable et al., 1987). Our approach, in addition to inverting the data using geometric constraints, seeks a functional relationship simultaneously for any pair of model parameters. The relationship is treated as a random variable as part of the joint inversion algorithm. Carrillo and Gallardo (2018) propose a general form of the relationship between two model parameters (m 1 and m 2 ): p and q are the maximum power allowed for each parameter, and a i j are the polynomial coefficients linking m 1 and m 2 . They introduce the function g:

g(m 1 , m 2 , a) = q i=0 p j=0 a i j m i 1 m j 2 = 0. ( 1
)
Our work focuses on exploring the use of correspondence maps to invert for the one-dimensional shear velocity from dispersion curve data (from ambient seismic noise) and the resistivity (from magnetotelluric data). We first illustrate the method of synthetic data using linear (polynomial of degree 1) and non-linear relationships (polynomial of degree higher than 1). We then apply our methodology to field data. Throughout the paper, we use the same notation as Carrillo and Gallardo (2018) and Constable et al. (1987) to help the reader to relate our work to these seminal articles.

DATA S E T S Magnetotelluric data

Magnetotellurics (MT) is one of the most widely used geophysical methods to obtain the Earth's subsurface resistivity structures. The principle of MT is to acquire simultaneous measurements of the horizontal components of the electric field E and the magnetic field H at the Earth's surface. In the Fourier domain, the transfer function between these components of E and H yields the complex impedance tensor Z from which the apparent resistivity ρ a and phase φ a can be computed:

E(ω) = Z(ω)H(ω), (2)
ρ a (ω) = |Z(ω)| 2 μω , φ a (ω) = arg Z(ω), (3)
where ω is the angular frequency (i.e. 2π f ) and μ is the magnetic permeability. The impedance tensor is frequency dependent, the impedances at high frequencies correspond to shallow structures, while impedances at low frequencies correspond to deeper structures. The reader is referred to Vozoff (1990) for more details about MT. We calculate synthetic impedances for a one-dimensional isotropic layered earth by computing the effective impedance at the surface of the Earth using the recursive approach described in, for example, Ward and Hohmann (2012). To ensure stability and avoid negative values of resistivity during the inversion process, we use the natural logarithm of ρ as inversion parameters in all our computations.

Rayleigh wave dispersion data from ambient seismic noise correlation

For many events in a seismic record, the surface waves are the largest amplitude feature present in the seismogram. Surface Rayleigh waves are the result of interfering P-Sv waves. Particle motion of the fundamental mode of Rayleigh waves moving from left to right is elliptical in a retrograde (counterclockwise) direction. The motion is constrained to the vertical plane that is consistent with the direction of the wave propagation [START_REF] Xia | Estimation of near-surface shear-wave velocity by inversion of rayleigh waves[END_REF]. Surface waves are guided and dispersive. Just as is the case for the electromagnetic waves in MT, lower frequencies penetrate deeper than higher frequencies for a given mode. Analysis of the kernels of surface wave dispersion curves shows that they are primarily sensitive to shear wave velocity variations and less sensitive to compressional velocity and density [START_REF] Dorman | Numerical inversion of seismic surface wave dispersion data and crust-mantle structure in the New York-Pennsylvania area[END_REF][START_REF] Xia | Estimation of near-surface shear-wave velocity by inversion of rayleigh waves[END_REF]Julià et al., 2000;[START_REF] Song | Sensitivity of near-surface shear-wave velocity determination from Rayleigh and Love waves[END_REF].

Commonly, subsurface reservoir imaging is based on active seismic sources, but more recently the cross-correlations of ambient seismic noise records have shown the capability of empirically obtaining Green's function between a pair of receivers (Shapiro et al., 2005), which has become a standard technique for passive seismic imaging. In this paper, we invert only for shear wave velocity and compute compressional velocity from a given V p /V s ratio. The forward computation of dispersion curves is done using the disba python library (Luu, 2019), which is based on the code sdisp96 from Computer Programs in Seismology (Herrmann, 2013). Similar to MT data and resistivity models, we use the natural logarithm of shear velocity in all computations for inversion.

J O I N T I N V E R S I O N M E T H O D

We have chosen Occam's inversion method (Constable et al., 1987) to solve the joint inversion objective function under correspondence map constraints. From equation (1), g is equal to zero in the general case, it is convenient to manipulate it algebraically to avoid the trivial solution, that is when all the coefficients describing the polynomiae are equal to zero. Carrillo and Gallardo (2018) circumvent this problem by normalizing equation ( 1) by the independent term coefficient a 00 . The normalized equation ( 1

) becomes ĝ(m 1 , m 2 , a) = q i=0 p j=0 âij m i 1 m j 2 = 0, (4)
where âij = a i j a 00 ; a 00 = 0; â00 = 1.

(5) Below is given an example for degree 1, where p and q are equal to 1:

ĝ(m 1 , m 2 , a) = â00 + â01 m 0 1 m 1 2 + â10 m 1 1 m 0 2 + â11 m 1 1 m 1 2 = 0. (6)
Rearranging and substituting â00 of equation ( 5) in equation (6) give

ĝ(m 1 , m 2 , a) = â01 m 0 1 m 1 2 + â10 m 1 1 m 0 2 + â11 m 1 1 m 1 2 = -1. (7)
The joint inversion problem including correspondence maps using Occam's approach is defined as finding the simplest (i.e. smoothest) shear velocity (m 1 ) and resistivity (m 2 ) models, which can be correlated by a polynomial relationship and reproduce the observed data. The problem is set to find m 1 , m 2 and the normalized coefficients âij . Therefore, we propose the objective function as an expanded version of the matrices described in the main equation of Occam, that is

m = [λ∂ T ∂ + (WJ) T WJ] -1 (WJ) T W(d -F[m i ] + Jm i ). ( 8
)
Here m is the model parameter to find and m i is the model at iteration i, λ∂ T ∂ is the smoothing matrix (here the difference in a physical property between consecutive layers); λ is the commonly known Lagrange multiplier, a large λ is a smooth model; J is the sensitivity matrix, and d and F[m] are the observed and calculated data, respectively. W is the dataweighing matrix, modified from the one in Constable et al. (1987) by including a factor F w , to account for each of the three data types. This factor ensures that a data type with more elements than the others will not inherently bias the inversion. W is then defined W = diag (F w /σ ), where σ is the uncertainty associated with each data value and F w :

F w = A D N N D , ( 9
)
where A D is the weight given to data type D with 3 D=1 A D = 1. N D is the number of data for each type, and N = 3 D=1 N D is the total number of data. An analysis of the effect of different A D combinations is presented in Appendix C. Equation (8) is computed over several values of λ in order to find the smoothest model satisfying the misfit:

X = ||Wd -WF[m]||. ( 10
)
The modified model matrix m is defined as a matrix that comprises shear velocity, m 1 , resistivity, m 2 , and the normalized polynomial coefficients, âij

m = ⎛ ⎜ ⎝ m 1 m 2 â ⎞ ⎟ ⎠ . ( 11
)
Similarly, the observed data d and predicted F[m] are extended to

d = ⎛ ⎜ ⎝ d 1 d 2 -1 ⎞ ⎟ ⎠ , ( 12
)

F[m] = ⎛ ⎜ ⎝ F SW D [m 1 ] F MT [m 2 ] ĝ[m 1 , m 2 , â] ⎞ ⎟ ⎠ , ( 13
)
where J is the general Jacobian matrix that integrates the sensitivity of the given model parameters. More details regarding the Jacobian matrices of the correspondence map function ĝ can be found in Carrillo and Gallardo (2018). The complete J matrix given as the model parameters in equation ( 11) and data in equation ( 12) is

J = ⎛ ⎜ ⎝ J[m 1 ] 0 0 0 J[m 2 ] 0 J[g m1 ] J[g m2 ] J[g a ] ⎞ ⎟ ⎠ . ( 14
)
From equation ( 8), the first term (μ∂ T ∂) is related to the smoothness of the model parameters. In our case of joint inversion using correspondence maps, there is no physical or mathematical reason to smooth the coefficients, and hence a smoothness of zero is imposed on the coefficients term. The smoothing term in equation ( 8) is then as follows:

∂ = ⎛ ⎜ ⎝ ∂m 1 ∂m 2 ∂a ⎞ ⎟ ⎠ = ⎛ ⎜ ⎝ ∂m 1 ∂m 2 0 ⎞ ⎟ ⎠ . ( 15
)
The experiments presented in the next section include the retrieval of all model parameters described in equation ( 11). The aim is to compare separate versus joint inversions for shear velocity and resistivity using the fundamental mode of Rayleigh wave dispersion curves and magnetotelluric data. The workflow followed in all experiments is described in Figure 1.

S Y N T H E T I C T E S T S

To verify the reliability and stability of the inversion algorithm, we test several scenarios of correspondence map relationships using synthetic data: (1) linear relationship, (2) second-order relationship and (3) second-order relationship from well log data. All experiments start with a homogeneous model that does not have to be close to the true model. For the case of the coefficients, the initial model is a vector of ones. The inversions are not constrained by any upper or lower bounds, neither on shear velocity nor resistivity, enabling the correspondence maps to search over a wide range of values. In all experiments, m 1 = ln V s and m 2 = ln ρ and 5% random noise is systematically added to the synthetic data. We calculate the synthetic responses for the magnetotelluric (MT) data over logarithmically spaced frequencies between 0.01 and 100 Hz and surface-wave dispersion curves from 0.1 to 10 Hz. The synthetic data are denoted as Target in Figures 2, 3 and7.

Linear relationship: Synthetic test

Although most seismic velocity-resistivity relationships are generally not linear as shown in Carcione et al. (2007), we start the joint inversion tests with the most simple scenario of correspondence map, that is a linear relationship. Linear relationships were observed previously in field data: for example, Harmon et al. (2021) recovered linear relationships in the mantle for shear velocity and log resistivity near the Mid-Atlantic Ridge. The cross-plots of these properties, although not true everywhere, show an overall linear trend.

The first experiment consists of a resistivity model (m 2 ) computed from a set of shear velocity values (m 1 ) using a linear relationship. The true linear relationship from which m 2 is computed for this experiment is

m 2 = 6m 1 + 1.3. ( 16
)
The normalized correspondence map function to search has the form of equation ( 16), namely,

ĝ(m 1 , m 2 , a) = â01 m 2 + â10 m 1 = -1. (17)
Figure 2 shows the resulting models and responses from the joint and separate inversions. The shear velocity model from the joint inversion is much closer to the original model, especially below 1 km depth. Sensitivity kernels show that this is the depth (corresponding to periods between 5 and 10 s) from which sensitivity to shear velocity decreases. The correspondence map plots show how the models from the joint inversion align better with the true linear relationship than those obtained from separate inversions. The retrieved parameter relationship after 15 iterations is

m 2 = 6.13m 1 + 1.29, (18)
that is very close to the true relationship in equation ( 16).

Non-linear relationship: Synthetic test

The second experiment consists of a resistivity model (m 2 ) computed from a set of shear velocity values (m 1 ) using a degree 2 polynomial relationship. This type of relationship has been observed in geothermal environments, where shallow formations dominated by sediments show a decreasing trend in resistivity and an increment towards deeper granitic formations with a persistent increase in the velocity throughout the whole section (Glass et al., 2018), resulting in nonlinear parameter relationships. Carcione et al. (2007) showed several examples of non-linear relationships between velocity and conductivity in oil-saturated shales and sandstones.

The true relationship from which m 2 is computed for this experiment is For this experiment, the normalized correspondence map function to search has the following form:

m 2 = 26m 2 1 -23m 1 + 6. ( 19
)
ĝ(m 1 , m 2 , a) = â20 m 2 1 + â10 m 1 + â01 m 2 = -1. (20)
After 15 iterations, the recovered relationship is as follows:

m 2 = 25.7m 2 1 -21.5m 1 + 5.55. ( 21
)
The correspondence map plot in Figure 3 shows that the joint inversion (magenta) from equation ( 21) is more consistent with the true relationship (black) from equation ( 19) than the results from separate inversions after 25 iterations (blue), which are more scattered, especially for the deeper section of the shear velocity model.

Discussion of linear and non-linear relationships from synthetic tests

To describe how accurate the joint and separate inversions results from the synthetic models are, we assume three factors of the inversion algorithm that are inherently correlated: smoothing, correspondence maps (resulting models) and the data (misfit). In the joint inversion using correspondence maps we are able to retrieve parameter relationships without compromising the smoothing and the data fit.

In Figure 4, we can see how joint inversion results in magenta represent a closer shape to the true relationship of the models than those of separate inversions for both linear and second degree synthetic tests. In order to quantify how close the correspondence maps obtained from separate and joint inversions was to the true relationships, we compute the ĝ function using the pairs of models retrieved from each inversion and compare it to the target ĝ = -1 within a certain error or tolerance. The green shading in the background represents the ĝ value calculated for all (m 1 , m 2 ) pairs in a mesh using the true coefficients. We obtain, for the linear relationship test, that 100% of the joint inversion model pairs (shear velocityresistivity) have ĝ values, computed using the true relationship, between -1.05 and -0.95 versus only 15% of the separate inversion pairs. In the non-linear relationship synthetic test, 80% of the joint inversion pairs but 45% of the separate inversion pairs are observed. These joint inversion results lead to a meaningful improvement in the models' relationships without compromising the data fitting. Hence we can say that joint inversion using correspondence maps leads to a reduction in the number of acceptable models that could explain the same observed data ('reduction of non-uniqueness').

Non-linear relationship: Well log data

The third test consists of a numerical experiment using subsurface parameters extracted from well log data. The first step refers to the problem of defining a proper degree of the polynomial for joint inversion and the second step to the actual results of the joint inversion. The data from well API 15023214760000 were downloaded from the publicly available Kansas Geological Survey database (http://www.kgs.ku. edu/Magellan/Logs/). Among other logs, it includes deep resistivity and compressional slowness logs. For the purpose of this exercise, shear velocity is computed using a constant V p /V s ratio of 1.7. Contrary to the previous two experiments, the ρ-V s relationship here is unknown and hence is closer to a field data case.

Since the correspondence map joint inversion requires the definition of the maximum degree (p and q) in equation ( 4), the first step is to seek what degree of the correspondence map function is adequate for the joint inversion based on the well log data. For this, we perform an inversion for only the coeffi- By constrained, we mean that not all the terms of equation (4) are used, so that m 2 can be expressed as a function of m 1 . Therefore, we limit or constrain which terms of equation (4) to be used in the joint inversion problem. Appendix A shows the constrained and non-constrained form of ĝ functions used for this study. On the other hand, non-constrained means that the ĝ function uses all the terms of equation ( 4). The results of these four separate inversions of coefficients are displayed in Figure 5. The background colour in each subplot is ĝ(m 1 , m 2 , a) evaluated for each m 1 and m 2 pair from a mesh, using the coefficients recovered for the four inversion tests described above.

Figure 5(a,c) shows that more than half of physical property pairs fall outside of the recovered relationships, suggesting that the relationship between shear velocity and resistivity cannot be explained using a linear relationship (degree 1). Therefore, these two scenarios were discarded for the joint inversion.

Although Figure 5(b) covers most of the model parameter pairs, when framed into a joint inversion problem, it can be-come unstable due to the high number of coefficients to search (i.e. eight for the case of p = 2 and q = 2). Appendix A, as shown in Figure A.1, exemplifies this subject, which has also been reported in Carrillo and Gallardo (2018). Finally, it is scenario (d) which is used to test the joint inversion strategy since it is the relationship that covers most of the well log model pairs after (b) and explains more closely than (a) and (c) the behaviour of the well log data.

After analysing the degree of polynomial to be used, the second step consists of the actual joint inversion using correspondence maps. For modelling purposes, the well log models are binned into 150 m thick layers to make the computations faster. In Figure 6, it is noted that root mean square (RMS) misfits for the three data types (MT, surface wave dispersion curves and g) marginally decrease after the fifth iteration for the joint inversion test. Figure 7 shows that both joint inversion and separate inversion fit the data similarly. From the parameter side, the shallow section represents a model where velocity increases and resistivity decreases, while in the deeper section both properties increase. The joint inversion approach tries to couple these two behaviours as a parabolic form. The recovered relationship for the joint inversion scenario (d) after 15 iterations is as follows: m 2 = 5.32m 2 1 -4.87m 1 + 2.14. ( 22)

F I E L D E X A M P L E

The previous experiments were necessary to assess the reliability of our methodology. As a first application to a real data case, we selected data from northern Alsace, France. This area has historically been explored for hydrocarbons, and more recently hosts the development of deep geothermal systems. A vast amount of geophysical data have been acquired in the region, making it ideal for testing new interpretation and imaging techniques such as the one presented in this paper. The datasets tested in this paper consist of one pair of magnetotellurics (MT) sounding from the ECOGI project [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF] and fundamental mode Rayleigh wave group dispersion curves from the EstOF network (Lehujeur et al., 2015(Lehujeur et al., , 2018)). The selected pair of stations (MT-SWD) have a surface separation of at most 150 m. The processed MT data [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF] have a frequency range from 0.01 to 100 Hz. The seismological data processed by Lehujeur et al. (2015) and Lehujeur (2015) consist of fundamental and first overtone Rayleigh-wave dispersion curves in the range of 0.2-1.28 Hz (0.7-5 s).

Generally, all components of the MT impedance tensor are non-zero and require a three-dimensional modelling approach to be reproduced. Since the scope of this work is one dimensional (1D), we choose an MT station for which the phase tensor ellipse (Caldwell et al., 2004) is closest to a circle over the broad frequency range of interest, that is suitable for our 1D modelling and inversion approach.

A normalized correspondence map function like equation (20) is used for this dataset, since the results from separate inversions (blue in Fig. 8) show a parabola in a cross-plot. Models obtained by [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF] and Lehujeur et al. (2018) are displayed as reference in grey. Note that the resistivity model from [START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF] was obtained from a two-dimensional inversion using MARE2DEM Key (2016), Figure 8 displays the extracted model at the station selected in this joint inversion test. The shear velocity model displayed as reference was obtained from Lehujeur et al. (2018), who used a Monte Carlo inversion approach from where the median of the best 2000 models was retained as the solution of the inversion.

The results of the joint inversion (Fig. 8) confirm the use of a non-linear parameter relationship. The introduction of the correspondence map in the joint inversion algorithm results in a shear velocity model below 2.5 km depth closer to the shear velocity model from Lehujeur et al. (2018), compared to that obtained by separate inversion. Note, however, that in Lehujeur et al. (2018) bounds for shear velocity based on well log data were used as input due to the nature of the Bayesian approach used in their work. In our case, we did not provide any prior information or bounds for shear velocity and resistivity and half-space initial models were used for both parameters. The recovered relationship after 15 iterations was as follows:

m 2 = 2.49m 2

Figure 9

Results of joint inversion of data in Section 4.2 using degree 1 and degree 2, correspondence maps (left) and RMS misfit evolution (right). Non-convergence is observed when performing joint inversion using degree 1.

E F F E C T O F T H E D E G R E E O F g(m 1 , m 2 , a) IN T H E J O I N T I N V E R S I O N

We have seen that the degree of the polynomial relationship used for the joint inversion is an important parameter. Therefore, it is worth investigating how to choose this parameter and how the joint inversion is affected by it. This study started with joint inversions from synthetic models and responses, where predefined relationships were used, that is linear and non-linear. In each of these two synthetic tests, the chosen relationship for the joint inversion was inherently used as the 'correct one', for instance, the non-linear relationship joint inversion used a degree 2 relationship. To study the effect of using different degrees, we took the same synthetic models from Section 4.2 and perform joint inversions using purposely a wrong degree for relationship. That is, for Section 4.2, where the known relationship is non-linear (degree 2), we perform a joint inversion using degree 1. Figure 9 shows the cross-plot ρ-V s and the root mean square (RMS) evolution at each iteration for this test. We observe that the joint inversion using degree 1 does not converge versus using a degree 2, confirming that the ρ-V s relationship cannot be explained by a linear function, and therefore a higher order is required. This suggests that when the relationship is unknown, a higher order relationship should be initially tested, followed by decreasing degrees to simplify the solution. If the simpler relationship can explain the data, the joint inversion will converge; otherwise, the higher order relationship must be kept.

As logical approach to decide the degree of joint inversion, one can think of looking at the correspondence maps from separate inversions and use this as a starting point to estimate the degree to use. If we look at the field example in Section 5, both the reference models and the models from separate inversions describe a non-linear relationship and this was the reason why degree 2 was chosen for joint inversion. Appendix B shows the result of the joint inversion for this dataset but using a linear relationship instead. Figure B.1 shows how the fit of the linear relationship joint inversion for the three datasets (dispersion curves, MT: apparent resistivity/phase and ĝ) is significantly lower than when using degree 2.

C O N C L U S I O N S

To the best of our knowledge, we have successfully applied a correspondence map joint inversion between surface-wave dispersion curves and magnetotelluric data for the first time.

It is an effective way to find meaningful physical parameter relationships while retrieving the shear velocity and resistivity models. The inclusion of correspondence maps in the joint inversion problem has shown to converge without the definition of prior bounds as input and by using homogeneous initial models.

The results of the first two synthetic data examples show how the recovered relationships were essentially the same as the true relationships. They also prove that the models retrieved from joint inversion are importantly more similar to the true models than those obtained from separate inversions. We observed that when evaluating the ĝ function using model pairs from joint and separate inversions, a significantly higher percentage of model pairs from joint inversion had a ĝ between -1.05 and -0.95, reducing the number of acceptable models without compromising the data fit ('reduction of non-uniqueness'). The first part of the well log example, the assessment of separate inversions of only the coefficients, demonstrates an effective way to evaluate and select an appropriate polynomial degree for the joint inversion.

The field data experiment illustrates how joint inversion with correspondence maps achieves a velocity model similar to that obtained by a Bayesian approach (Lehujeur et al., 2018), without forcing boundaries to the model, but rather by imposing the existence of a parameter relationship. It exemplifies how, if a cross-property mathematical correspondence exists that can explain observed apparent resistivity and dispersion curves, the resulting models have a higher coupling than when not using the correspondence map. The addition of a proper correspondence map term has been shown not to compromise the fit of the data in the joint inversion algorithm. It is essential, however, to remind the reader that the convergence of the joint inversion is inherently related to the appropriate selection of the polynomial degree, as described in Section 6. In order to evaluate the stability of the joint inversion algorithm for a degree 2 using the full expression of ĝ (p = 2 and q = 2), we perform a joint inversion test using the data and model responses from the well log data in Section 4.4. We produce two joint inversion tests: (a) joint inversion using only three coefficients of equation ( 4), which are described in equation (A.1), and (b) using the full eight coefficients (p = 2 and q = 2) of equation ( 4), developed in equation (A.2). The experiments mentioned in this article as constrained use the normalized correspondence map function using a limited number of polynomial coefficients, like test (a). On the other hand, the experiments called non-constrained use the full expression of the correspondence map function as in equation ( 4 The results in Figure A.1 show that the joint inversion can be stable using a restricted number of coefficients, that is the three coefficients in equation (A.1), while it can take more iterations for the non-constrained form to converge or not even converge. For this test, after 15 iterations the non-constrained inversion did not converge to a solution as the constrained inversion did. Carrillo and Gallardo (2018) described instability in the joint inversion of two-dimensional models (gravity and magnetotellurics) when using a large number of coefficients, that is eight coefficients (p = 2 and q = 2). We are also interested in evaluating the influence of the chosen degree of the polynomial relationship (p and q values in equation ( 4)) on the correspondence map joint inversion. We hence perform two more joint inversion tests on the same field data: (i) correspondence map joint inversion using a maximum degree of 1 (linear relationship) and (ii) degree 2 (parabolic relationship). The parabolic relationship is described in Section 5. It is expected that using a higher order relationship would show better results than including a linear relationship in the joint inversion, since the reference models do not show any hint of linearity.

The results are consistent with our expectations. The linear relationship correspondence map is never able to converge into a solution that can simultaneously fit the data while finding a relationship. This in- Joint inversion using correspondence maps 15 dicates that even if the models can be improved via joint inversion using correspondence maps the convergence of the inversion is restricted to a proper selection of the degree of the relationship. The model parameters of this field dataset require a higher polynomial degree in order to fit both observed datasets (d 1 and d 1 ) while explaining a common relationship as observed in

A P P E N D I X C : W E I G H T A NA LY S I S O F T H E DATA I N T H E J O I N T I N V E R S I O N F R A M E WO R K

We investigate the effect of the weight, A in equation ( 9), of each dataset in the joint inversion framework. Four scenarios of distribution of weights are tested and described in Table C.

  d'intérêt (<10 s) et les données MT porraient ici améliorer les choses. Des tests synthétiques d'inversion conjointe en 2D utilisant une relation au premier ordre montrent que le modèle de vitesse est sensiblement amélioré par l'inclusion de la magnétotellurique et la relation vitesse-resistivité est mieux recouvrée qu'avec des inversions séparées. L'application de la méthode d'inversion conjointe à des données de terrain du Nord de l'Alsace illustre ses avantages et ses inconvénients. Pour une zone avec une subsurface aux propriétés quasi-1D, l'inclusion des cartes de correspondance dans l'inversion résulte en des modèles pétrophysiquement et structurellement plus cohérents. Dans le cas d'une zone plus complexe, où plusieurs relations pourraient être observées, le terme CM a n'est pas en mesure d'amener une relation unique valable pour toutes les zones structurellement différentes. Mots clés: Inversion jointe, cartes de correspondance, magnetotellurique, bruit ambiant sismique, ondes de surface. v "It is vain to do with more what can be done with less"

Figure 1 . 1 :

 11 Figure 1.1: Distribution of published articles of Magnetotelluric (orange) and Ambient Seismic Noise inversion (blue). Data obtained from https://www-webofscience-com, from 1977 until 2022.

Figure 1 . 2 :

 12 Figure 1.2: Distribution of joint inversion articles published according to different geophysical methods. Data obtained from https://www-webofscience-com, from 1977 until 2022.

  Figure 1.3: Manuscript Organization

  method (MT) . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Electric and magnetic fields in a homogeneous medium . . . . . 2.1.2 Relationship between electric and magnetic fields . . . . . . . . . 2.1.3 Apparent resistivity . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.4 Impedance tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 MT data workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 MT data acquisition . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 MT data processing . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 MT forward modeling codes . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 One-dimensional modeling: in-house code . . . . . . . . . . . . . 2.3.2 Two-dimensional modeling: PW2D . . . . . . . . . . . . . . . . . 2.4 Ambient seismic noise correlation (ASNC) . . . . . . . . . . . . . . . 2.4.1 Sources of ambient seismic noise . . . . . . . . . . . . . . . . . . 2.4.2 Cross-correlation of noise . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Practical limitations . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 ASNC data workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Surface wave dispersion curves (SWD) from ASNC: data processing 2.6 SWD forward modeling codes . . . . . . . . . . . . . . . . . . . . . . 2.6.1 One-dimensional modeling: disba module . . . . . . . . . . . . . 2.6.2 Three-dimensional modeling: SW4 . . . . . . . . . . . . . . . . . 2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 . 1 :

 21 Figure 2.1: Amplitude spectrum of geomagnetic variations. Reproduced from Constable[2015]

1 and substituting Chapter 2 .

 2 Passive geophysical methods the constitutive relation for the magnetic flux density B yields ∇ × (∇ × E) = -µ(∇ × ∂H ∂t )

Figure 2 . 2 :

 22 Figure 2.2: Three-layer earth section, resistivity varying in the vertical direction only.From Cagniard [1953]

2. 3 .

 3 The data presented in this Thesis have been recorded using channels Ex, Ey, Hx and Hy from Metronix (ADU-07e and ADU-08). The equipment used is part of the EMMOB National equipment pool hosted by our laboratory. Magnetic probes type MFS-06 and MFS-06e are buried for measuring the horizontal and vertical magnetic components. The electric channels are deployed as a 50 m dipole using EFP-06 non polarized electrodes.Data is recorded using 512 Hz and 4196 Hz sampling rates. An example of time series recorded with this configuration can be seen in Figure2.4. It shows the four components acquired (H x ,H y ,E x and E y ) for a 3 s window zoom.

Figure 2 .

 2 Figure 2.3: a) MT Diagram of disposition of electrodes and magnetometers and b) A typical MT field disposition.

Figure 2 . 5 :

 25 Figure 2.5: Bounded influence results using Razorback for a site with multiple reference stations, no error bars displayed. All combinations of RR stations are computed, with color code corresponding to the associated combination of RR stations. Reproduced from Smaï and Wawrzyniak [2020]

Figure 2 . 7 :

 27 Figure 2.7: Sensitivity kernels for resistivity (left) and phase (right) for 1D layered model in 2.6. Color bar is dimensionless.

Figure 2 . 8 :

 28 Figure 2.8: Example of a two-dimensional resistivity model used in PW2D.

Figure 2 .

 2 Figure 2.7 shows the sensitivity over frequency and depth for both apparent resistivity and phase.

Figure 2 .

 2 Figure 2.10: Power spectral density estimated over noise data from the year 1995 for the three components of all continental GEOSCOPE stations modified from Stutzmann et al.[2000]. Dashed black lines show the extreme noise models from[START_REF] Peterson | Observations and modeling of seismic background noise[END_REF]. The shaded blocks in green, black, red and blue represent different sources process described in the text.

Figure 2 .

 2 Figure 2.11 a) illustrates the field produced by sources located around a circle. The signals are cross-correlated with a reference point A(+). b), c) and d) are the crosscorrelations at t=-30, 0 and +30 s respectively. At t=0 all the energy is focused in A, as if the source was in A. For positive d) and negative b) correlation times a wavefront can be seen to diverge and converge respectively. They represent the causal and acausal parts of the Green's function between A and any point in the medium. Both the causal and acausal parts are approximations of the Green's function and in theory if the distribution of sources (either primary sources or "secondary" diffracted sources) is homogeneous we expect this symmetry of the correlation functions.

Figure 2 .

 2 Figure 2.11: (a) Configuration of the numerical experiment. One thousand sources S(X) are surrounding the reference point A (+). Dots indicate the point scatterers. (b) Snapshot of the cross-correlation between the field in A with the field at location (x, y) after averaging over the sources S for correlation time -30 s. A converging wave front is well defined and constitutes the anti-causal part of the Green function. (c) Snapshot for correlation time t=0 s. The wave front is focused on A. (d) Snapshot for t=30 s. The diverging wave front corresponds to the causal part of the Green function . Reproduced from Paulet al. [2005]

Figure 2 .

 2 Figure 2.13: 1D layered model and Rayleigh-wave group velocity obtained with disba

  Fig 2.13 shows the model response, here the fundamental model of rayleigh wave group velocity, for a simple 1D layered model. Sensitivity kernels of group velocity relative to P-wave, S-wave velocities and density can also be obtained. They are estimated by perturbing each of the model parameters in each layer and measuring the relative changes induced on the dispersion curves at each period.

Figure 2 .

 2 Figure 2.14: Sensitivity kernels (shear velocity, compressional velocity and density) for 1D layered model in 2.13. Color bars are dimensionless.

Figure 2 .

 2 15 exemplifies this step i.e. the effect of moving the source at each receiver each time.

Figure 2 .

 2 Figure 2.15: Schematic representation of SW4 modeling using only three receivers (gray triangles) as example, a) shows the receivers disposition. b), c) and d) show how the sources (yellow symbols) are displaced at different receivers locations each time and are never active at the same time (b), c), and d)).

Figure 2 .

 2 Figure 2.16: Schematic representation of the SW4 workflow illustrating the three main steps i) Seismograms obtained from SW4 modeling, ii) dispersion curves computations using FTAN and iii) Regionalization.

Figure 3 . 1 :

 31 Figure 3.1: Results of inversion to solve for coefficients for a set of model parameters described by a circular function as equation 3.9. In a) the parameters used as input m 1 and m 2 , with (red) and without noise (blue). In b) the Model (coefficients) recovered from the two tests of inversions compared to the true ones (black). In c) the data computed using the recovered coefficients in b) and the parameters in a). The coefficients in b) are detailed in Table3.1

Figure 3 . 2 :

 32 Figure 3.2: Left, input well-log data used for inversion test (top) and the cross-plots(bottom). Right, the data from inversion for the three combinations of pairs (resistivitydensity, density-resistivity and velocity-density), in orange using p = 1 and q = 1 and in green with p = 2 and q = 2 in equation 3.5, black line is the observed data, i.e. ĝ0 = -1

Figure 3 . 3 :

 33 Figure 3.3: Input random log data used for inversion test and results, description same as Figure 3.2

  Figure 3.4 shows two examples of parameter relationships and their respective MI and VI values. Case 1

Figure 3 . 4 :

 34 Figure 3.4: Examples of computation of variation information. Case 1 shows a linear relationship case between properties m 1 and m 2 and Case 2 exemplifies a randomized relationship between the two variables. Upper panels are the cross-plots and lower panels the histograms for both cases. Color bar represents the number of m 1 -m 2 pairs that fall in each cell i.e. frequency. The VI computation code was provided within the 3D Earth Spring School that I attended during my PhD Thesis organized by ESA and 3D Earth group https://www.3dearth.uni-kiel.de/en/3d-earth-summerschool

5

 5 [START_REF] Haber | Model Fusion and Joint Inversion[END_REF]. Model Fusion and Joint Inversion. Surveys inGeophysics, 34:675-695. Mandolesi, E. and Jones, A. G. (2014). Magnetotelluric inversion based on mutual information. Geophysical JournalInternational, 199(1):242-252. Molodtsov, D., Kashtan, B., and Rostov, Y. (2011). Joint inversion of seismic and magnetotelluric data with structural constraint based on dot product of image gradients. In SEG Technical Program Expanded Abstracts 2011, pages 740-744. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 The method: joint inversion using CM . . . . . . . . . . . . . . . . . 4.3 Synthetic tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 First-order relationship: synthetic test . . . . . . . . . . . . . . . 4.3.2 Higher-order relationship: synthetic test . . . . . . . . . . . . . . 4.3.3 Discussion on first and higher-order relationships from synthetic tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.4 Unknown relationship: well-log data . . . . . . . . . . . . . . . . Effect of the degree of ĝ(m 1 , m 2 , â) in the joint inversion . . . . vs higher-order correspondence map joint inversion in field data 4.C Weight analysis of the data in the joint inversion framework . . . . 4.D What if there is no relationship? . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 D=1N

 3 D is the total number of data. An analysis of the effect of different A D combinations is presented in Appendix 4.C. Equation 4.1 is computed over several values of λ in order to find the smoothest model satisfying the misfit:

Figure 4 . 1 :

 41 Figure 4.1: Workflow of joint inversion using correspondence maps; m 1 and m 2 are the log of shear velocity and log of resistivity respectively.

Figure 4 . 2 :

 42 Figure 4.2: Correspondence map Joint Inversion results of Section 4.3.1: shear velocity and resistivity models; synthetic data and models. Color codes: synthetic (black), starting model (red), separate inversions (blue) and joint inversion (magenta).

Figure 4 . 3 :

 43 Figure 4.3: Correspondence map Joint Inversion results of Section 4.3.2: shear velocity and resistivity models; synthetic data and models. Color codes: synthetic (black), starting model (red), separate inversions (blue) and joint inversion (magenta).

Figure 4

 4 Figure 4.4: Correspondence map plot for results in Section 4.3.1 in a) and b) and for Section 4.3.2 in c) and d). Background mesh in green-shades, is the corresponding true relationship g(m 1 , m 2 , a), i.e. equation 4.9 for a) and b) and equation 4.12 for c) and d), evaluated over multiple values of m 1 and m 2 . The resulting models for separate (blue) and joint inversions (magenta) are displayed above. 100% and 80% of the model-pairs in b) and d) respectively have ĝ values between -1.05 and -0.95, while only 15% and 45% of the pairs in a) and c).

  5, the first step is to seek what degree of the correspondence map function is adequate for the joint inversion based on the well-log data. For this, we perform an inversion for only the coefficients âij using the well-log data (shear velocity and resistivity) and evaluate what type of relationship covers most of the well-log physical property pairs. Four scenarios of coefficient inversions are tested, varying the degree and the terms used in the ĝ(m 1 , m 2 , a) function, either the full expression (non constrained) or a constrained version of the function. The four scenarios of ĝ(m 1 , m 2 , a) functions used for testing are described below: (a) Degree 1 non constrained (p = 1 and q = 1 in equation 3.5) (b) Degree 2 non constrained (p = 2 and q = 2 in equation 3.5) (c) Degree 1 constrained (equation 4.10) (d) Degree 2 constrained (equation 4.13) By constrained we mean that not all the terms of equation 3.5 are used, so that m 2 can be expressed as a function of m 1 . Therefore we limit or constrain which terms of equation 3.5 to be used in the joint inverse problem. Appendix 4.A shows the constrained and non constrained form of ĝ functions used for this study. On the other hand, non constrained means that the ĝ function uses all the terms of equation 3.5. The results of these four separate inversions of coefficients are displayed in Figure 4.5. The background color in each subplot is ĝ(m 1 , m 2 , a) evaluated for each m 1 and m 2 pair, using the coefficients recovered for the four inversion tests described above.

Figure 4 .

 4 Figure 4.5 a) and c) show that more than half of physical property pairs fall outside of the recovered relationships, suggesting that the relationship between shear velocity and resistivity can not be explained using a linear relationship (degree 1). Therefore these two scenarios were discarded for the joint inversion. Although Figure 4.5 b) covers most of the model parameters pairs, when framed into a joint inverse problem, it can become unstable due to the high number of coefficients to

Figure 4 . 5 :

 45 Figure 4.5: Analysis of degree for correspondence map function on well-log data. Left hand side show the well-log data (gray). Green-shades background in a), b), c) and d) show the evaluated g(m 1 , m 2 , a) functions obtained from separate inversions of the coefficients â using degree 1 and 2 (non constrained) and degree 1, degree 2 (constrained) respectively. Displayed in all plots are the well-log data (gray) for reference. The g functions in the four panels, a) to d), were obtained from the separate inversions tests using the well-log data (gray).

Figure 4 .

 4 Figure 4.7 shows that both joint inversion and separate inversion fit the data similarly.

Figure 4 . 7 :

 47 Figure 4.7: Correspondence map Joint Inversion results of Section 4.3.4: shear velocity and resistivity models; synthetic data and model responses. Color codes: true model (black), starting model (red), separate inversions (blue) and joint inversion (magenta).

  Figure 4.8: Correspondence map Joint Inversion results of Section 4.4: shear velocity and resistivity models; field data. Color codes: reference models (gray), starting model (red), separate inversions (blue) and joint inversion (magenta).

  3.4. We produce two joint inversion tests: a) joint inversion using only three coefficients of equation 3.5 which are described in equation 4.17 and b) using the full eight coefficients (p = 2 and q = 2) of equation 3.5, developed in equation 4.18. The experiments mentioned in this article as constrained, use the normalized correspondence map function using a limited number of polynomial coefficients, like test a). On the other hand, the experiments called non constrained use the full expression of the correspondence map function as in equation 3.5, like test b).

Figure

  Figure 4.11: Correspondence map Joint Inversion results for data in Section 4.4 using Linear and Non-linear relationships.Color codes: reference models (gray), starting model (red), correspondence map joint inversion using linear relationship (green) and correspondence map joint inversion using higher-order relationship (magenta)

Figure 4 .

 4 Figure 4.12: RMS evolution for the four weight tests described in Table 4.1. The dashed black line represents the threshold when RMS=1. The title of each subplot indicates the value of the weights of each data set (3 values are given), the order is given as: index 1: weight of dispersion curves data set, index 2: weight of the magnetotelluric data set and index 3: weight of g data set

Figure 4 . 1

 41 Figure 4.14: Random models used for testing joint inversion using correspondence maps

Figure 5 .

 5 2 b), c) and d) are the synthetic seismograms obtained from three different source locations (beginning, middle and end of the profile). An example of a file used for computation for a source located in x=18 km is shown below. More details on the construction of the files can be found in geodynamics.org. # This is an SW4 input file for a scaled-down version of the LOH.1 problem from the # SCEC test suite, see # "Steven Day et al., Tests of 3D Elastodynamic Codes: # Final report for lifelines project 1A01, Pacific Eartquake 5.1. On the use of surface-waves for multi-dimensional models width of the super grid damping layer (30 is the default) supergrid gp=30# first define the material properties everywhere (half-space) all recording stations are relative to the (0,0,0) origin # Also note that the stations do NOT fall on a grid point when z=0 file=sta400 sacformat=1 usgsformat=1 rec y=3000 x=42000 z=0 file=sta420 sacformat=1 usgsformat=1 rec y=3000 x=44000 z=0 file=sta440 sacformat=1 usgsformat=1

Figure 5 . 3 :

 53 Figure 5.3: Raw seismograms (left) and their respective dispersion curves estimated with FTAN (right) for the three-layer models shown in Figure 5.1.Three station pairs are used with source at x s =18 km, receivers at x r =16 km for a) , x r =36 km in b) and x r =56 km for c). Curves in black in plots d), e) and f) are the estimated group velocities using FTAN methodology and in red from the 1D code ofHerrmann [2013]

Figure 5 .

 5 Figure 5.4: Source -receiver pairs contribution per period. Each panel shows the pairs contributing for an specific period 10 s,4.8 s, 2.19 s, 0.9 s, 0.44 s and 0.2 s. The x axis is the sources positions and the y axis the receivers. The station pairs with an offset < 2 wavelengths (λ) are discarded. Note that for T=10s no pair is selected.

Figure 5 . 5 :

 55 Figure 5.5: Group velocity map for the profile of model in Figure 5.1. Left: response using sdisp96 code from Herrmann [2013] and right: dispersion curves obtained from regionalization described in 5.1.2.1 of our synthetic data.

Figure 5 .

 5 Figure 5.6: 1D inversion of shear velocity using data of Figure 5.5 b). Cyan color crosses are the center of the grid cells used during regionalization for obtaining group velocity maps, and where 1D inversions were run.

Figure 5 . 7 :

 57 Figure 5.7: Shear velocity models used for synthetic sensitivity analysis. Black squares are the receiver stations i.e. the same as the sources. Cyan color crosses are the center of the grid cells used during regionalization for obtaining group velocity maps.

Figure 5 . 8 :

 58 Figure 5.8: Group dispersion maps for models in Figure5.7 using standard 1D model approach from sdisp96 codes[Herrmann, 2013]

Figure 5 .

 5 Figure 5.11: Shear velocity 1D inversion for four models of Figure 5.7 using data of Figure 5.9 b). Cyan color crosses are the center of the grid cells used during regionalization and obtaining group velocity maps, and where 1D inversions were done.

  sponse for the dispersion curves. A 2D model (m nz,nx ) is defined by n x and n z cells in the x and z directions. To build the pseudo-2D data vector, forward responses from 1D models positioned directly below the cells of the regionalization grid (nsite) are calculated. Figure 5.12 exemplifies the differences between 1D model compared to a pseudo-2D model.

Figure 5 .

 5 Figure 5.13: Shear velocity pseudo-2D inversion for four models of Figure 5.7 using data of Figure 5.9 b). Black dots are the positions were forward responses from 1D models are computed, equivalent to the cells defined during regionalization.

  14. The shear velocity model is composed of three values only i.e. 1.5, 2.0 and 3.0 km/s. The resistivity model has the following three values respectively: 250, 8 and 1500 Ωm. The position of the MT stations is set to be the same as the receivers deployed for the velocity models i.e. every 2 km (black squares in Figure 5.14).

Figure 5 .Figure 5 .

 55 Figure 5.14: Section 5.2.1 resistivity model (right) built using geometry from shear velocity iii) One block 10 km (left) and the corresponding parameter relationship plot below

Figure 5 .

 5 Figure 5.18: Separate and joint inversions results for models in Section 5.2.2. Black squares in the first two rows represent the receivers positions. Cyan crosses in velocity models are the center of the cells used for rationalizations i.e. the position of the dispersion curves used for inversions.

1 1

 11 Constable, S. C., Parker, R. L., and Constable, C. G. (1987). Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. GEO-PHYSICS, 52(3):289-300. deGroot-Hedlin, C. and Constable, S. (1990). Occam's inversion to generate smooth, two-dimensional models from magnetotelluric data. GEOPHYSICS, 55(12):1613-1624. Gaudot, I., Beucler, E., Mocquet, A., Drilleau, M., Haugmard, M., Bonnin, M., Aertgeerts, G., and Leparoux, D. (2020). 3-d crustal vs model of western france and the surrounding regions using monte carlo inversion of seismic noise cross-correlation dispersion diagrams. Geophysical Journal International, 224(3):2173-2188. Herrmann, R. B. (2013). Computer programs in seismology: An evolving tool for instruction and research. Seismological Research Letters, 84:1081-1088. General overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.1.1 Geological background . . . . . . . . . . . . . . . . . . . . . . . . 94 6.2 The datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.2.1 SWD data: EstOF project . . . . . . . . . . . . . . . . . . . . . . 94 6.2.2 MT data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . -dimensional inversions . . . . . . . . . . . . . . . . . . . . . 100 6.3.2 Two-dimensional inversions . . . . . . . . . . . . . . . . . . . -dimensional inversions . . . . . . . . . . . . . . . . . . . . . 106 6.4.2 Two-dimensional inversions . . . . . . . . . . . . . . . . . . . Profile 1: NW-SE . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.5.2 Profile 2: W-E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

  , and more recently for the development of deep geothermal systems. In 1987 the Soultz-sous-Forêts project started as a pilot test to exploit the potential energy of deep fractured rocks to produce electricity. It started as an HDR (Hot Dry Rock) concept and has over time migrated 91 Chapter 6. Application: North Alsace toward an EGS (Enhanced Geothermal System) concept [Dezayes

Figure 6 . 1 :

 61 Figure 6.1: Location of the two profiles studied in the northern Alsace area. White dots represent the location where local dispersion curves havebeen computed in the framework of the EstOF project[Lehujeur, 2015], orange square are the MT stations from ECOGI survey[START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF], yellow squares are MT stations from the DEEP-EM survey and purple points are the MT-Two extended MT stations. These two last datasets where acquired during the course of this PhD in 2020 and 2022 respectively. Background layer is the geological map obtained from BRGM

Figure 6

 6 Figure 6.4: Pseudosection of surface waves dispersion curves data used for Profile 1 (oriented NW-SE).

Figure 6 . 5 :

 65 Figure 6.5: Pseudosection of the SWD surface waves dispersion curves used for Profile 2 (oriented W-E).

  field MT data sets used in this Thesis are composed of three different surveys: i) ECOGI acquired in 2018, ii) DEEP-EM, acquired in October 2020 and iii) MT-Two, acquired in March 2022. Profile 1 uses only ECOGI data and Profile 2 uses data from DEEP-EM and MT-Two.

  EM and MT-Two. The DEEP-EM survey was acquired in October 2020 in the context of a broader electromagnetic campaign which included MT and CSEM (Controlled Source Electromagnetic) data. The survey benefited from the same remote reference as ECOGI plus a second magnetic remote reference station located about 40 km west of the survey area.The equipment used for the acquisition of these two surveys (DEEP-EM and MT-Two) was Metronix ADU-07e and ADU-08e and their sensors from the EMMOB National equipment pool hosted by our laboratory. We measured two components of electric field (E x and E y ) and two horizontal components of the magnetic field (H x and H y ). Two sequences were set up for recording, 512 Hz from 18:00 to 16:59 +1 day and 4096 Hz from 17:00 to 17:59. Schedule of the MT data acquisition, together with the apparent resistivity and phases obtained from Razorback can be found in Appendix A Section A.1.

Figure 6 . 6 :

 66 Figure 6.6: Pseudosection of the ECOGI MT data processed in Abdelfettah et al. [2019] and used for Profile 1 (oriented NW-SE).

Figure 6 . 7 :

 67 Figure6.7: Sensitivity kernels of Rayleigh-wave group velocity dispersion curves computed from a reference velocity model constructed from[START_REF] Maurer | On-going seismic monitoring of the rittershoffen egs project (alsace, france)[END_REF]

Figure 6

 6 Figure 6.8: Top-left: 1D correspondence maps, top-right: ĝ function and bottom: modelled vs observed data for separate (red) and joint (blue) inversions of Profile 1 discussed in Section 6.3

Figure 6 .Figure 6

 66 Figure 6.10: 1D inverted resistivity models for Profile 1 using only MT data. Left: inversion using a restricted frequency band (0.01-10 s) and right: using the full frequency band (0.01-64 s)

Figure 6 . 2 Figure 6

 626 Figure 6.13: 2D correspondence maps and data for Profile 1 in Section 6.3. Joint 1 label in blue refers to the joint inversion with weights equally distributed in the three data terms i.e. SWD:1/3, MT:1/3, CM:1/3. Joint 2 label in color green refers to the joint inversion with weights SWD:1/4, MT:1/4, CM:1/2

Figure 6 .

 6 Figure 6.15: Top-left: 1D correspondence maps, top-right: g function and bottom: modelled vs observed data for separate (red) and joint (blue) inversions of Profile 2 discussed in Section 6.4

Figure 6 .

 6 Figure6.17: 1D correspondence maps of Profile 2 in Section 6.4 same as in Figure6.15. Stations have been color coded according to their zone along the profile, west side (red), center (black) and east side (blue)

Figure 6 .

 6 Figure 6.19: Mesh designed for MARE2DEM inversions, a total of 8879 free parameters are used.

Figure 6 .

 6 Figure6.22: cross-plots of several physical properties. Left: correspondence maps of 2D joint inversion models of Profile 1 (shear velocity vs conductivity). Right: thermal conductivity as a function of compressional velocity for rock samples in the zone of study, reproduced from[START_REF] Kushnir | Characterizing the physical properties of rocks from the Paleozoic to Permo-Triassic transition in the Upper Rhine Graben[END_REF]

Figure 6 .

 6 Figure 6.24: Joint inversion operators computed for separate inversion models and the two tests of joint inversion. CG stands for cross gradients, DPG for the dot productimage gradients and VI for the variation information operator. Joint1 inversion refers to the case of weights equally distributed in the three data terms i.e. SWD:1/3, MT:1/3, CM:1/3. Joint2 inversion refers to the weights SWD:1/4, MT:1/4, CM:1/2

  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.

Figure A. 3 :

 3 Figure A.3: Apparent resistivity and phase curves obtained from Razorback for TE and TM modes for station 202 before and after corrected for time delay. The impedance tensor is importantly more coherent after the station is corrected by the time delay, and therefore apparent resistivity curves and phases are more similar to the reference stations, which were not affected by time delay.

Figure A. 4 :

 4 Figure A.4: Apparent resistivity and phase for site 204 computed with Razorback. R999 and R000 relate to the reference stations used i.e. Welschbruch and CIEL. The color code corresponds to an specific combination of RR used.

Figure A. 5 :

 5 Figure A.5: Apparent resistivity and phase for site 309 computed with Razorback. R999 and R000 relate to the reference stations used i.e. Welschbruch and CIEL. The color code corresponds to an specific combination of RR used.

Figure A. 6 :

 6 Figure A.6: Apparent resistivity and phase for site 318 computed with Razorback. R999 and R000 relate to the reference stations used i.e. Welschbruch and CIEL. The color code corresponds to an specific combination of RR used.

Figure A. 7 :

 7 Figure A.7: Apparent resistivity and phase for site 326 computed with Razorback. R999 and R000 relate to the reference stations used i.e. Welschbruch and CIEL. The color code corresponds to an specific combination of RR used.

Figure A. 8 :

 8 Figure A.8: Apparent resistivity and phase for site 331 computed with Razorback. R999 and R000 relate to the reference stations used i.e. Welschbruch and CIEL. The color code corresponds to an specific combination of RR used.

Figure A. 11 :

 11 Figure A.11: Apparent resistivity and phase for site 401 computed with Razorback. R999 relates to the reference station used i.e. Welschbruch. The color code corresponds to an specific combination of RR used.

Figure A. 12 :

 12 Figure A.12: Apparent resistivity and phase for site 402 computed with Razorback. R999 relates to the reference station used i.e. Welschbruch. The color code corresponds to an specific combination of RR used.

Figure B. 1 :

 1 Figure B.1: Example of 1D inversion using Occam algorithm described in Section B.3.Top left: shear velocity models (reference in black relates to that obtained from[Lehujeur et al., 2018]). Top right: RMS misfit evolution. Bottom left: data (group velocity).

Figure B. 8 :

 8 Figure B.8: Observed data, fundamental mode of Rayleigh-wave group velocity dispersion curves. The amplitudes correspond to velocity variations relative to a reference velocity indicated in the title of each subplot.

Figure B. 9 :

 9 Figure B.9: Forward response of models in Figure B.3. Same periods are displayed for reference as observed data in Figure B.8. The amplitudes correspond to velocity variations relative to a reference velocity indicated in the title of each subplot.

Figure 1

 1 Figure 1 Workflow of joint inversion using correspondence maps; m 1 and m 2 are the log of shear velocity and log of resistivity, respectively.

Figure 2

 2 Figure 2 Correspondence map joint inversion results of Section 4.1: shear velocity and resistivity models; synthetic data and models. Colour codes: synthetic (black), starting model (red), separate inversions (blue) and joint inversion (magenta).

Figure 3

 3 Figure 3 Correspondence map joint inversion results of Section 4.2: shear velocity and resistivity models; synthetic data and models. Colour codes: synthetic (black), starting model (red), separate inversions (blue) and joint inversion (magenta).

Figure 4

 4 Figure 4 Correspondence map plot for results in Section 4.1 in (a) and b) and for Section 4.2 in (c) and (d). The background mesh in green shading shows the corresponding true relationship g(m 1 , m 2 , a), that is equation (16) for (a) and (b) and equation (19) for (c) and (d), evaluating multiple values of m 1 and m 2 . The resulting models for separate (blue) and joint inversions (magenta) are displayed above. One hundred per cent and 80% of the model pairs in (b) and (d), respectively, have ĝ values between -1.05 and -0.95, while only 15% and 45% of the pairs in (a) and (c).

  cients âij using the well log data (shear velocity and resistivity) and evaluate what type of relationship cover most of the well log physical property pairs. Four scenarios of coefficient inversions are tested, varying the degree and the terms used in the ĝ(m 1 , m 2 , a) function, either the full expression (non-constrained) or a constrained version of the function. The four scenarios of ĝ(m 1 , m 2 , a) functions used for testing are described below: (a) Degree 1 non-constrained (p = 1 and q = 1 in equation (4)), (b) Degree 2 non-constrained (p = 2 and q = 2 in equation (4)), (c) Degree 1 constrained (equation (17)) and (d) Degree 2 constrained (equation (20)).

Figure 5

 5 Figure 5 Analysis of degree for correspondence map function on well log data. The left-hand side shows the well log data (grey). The greenshaded background in (a)-(d) shows the evaluated g(m 1 , m 2 , a) functions obtained from separate inversions of the coefficients â using degree 1 and 2 (non-constrained) and degree 1, degree 2 (constrained), respectively. Displayed in all plots are the well log data (grey) for reference. The g functions in the four panels, (a)-(d), were obtained from the separate inversions tests using the well log data (grey).

Figure 6

 6 Figure 6 As described in Section 4.4, RMS misfit evolution for joint inversion (left) and separate inversions (right).

Figure 7

 7 Figure 7 Correspondence map joint inversion results of Section 4.4: shear and resistivity models; synthetic data and model responses. Colour codes: true model (black), starting model (red), separate inversions (blue) and joint inversion (magenta).

Figure 8

 8 Figure 8 Correspondence map joint inversion results of Section 5: shear velocity and resistivity models; synthetic data and models. Colour codes: reference models (grey), starting model (red), separate inversions (blue) and joint inversion (magenta).

2 + â12 m 1

 21 ) as in test (b). The ĝ functions used in these two joint inversion tests are explicitly described below.(a) Constrained degree 2:ĝ(m 1 , m 2 , a) = â20 m 2 1 + â10 m 1 + â01 m 2 = -1 ( A . 1 ) (b) Non-constrained degree 2: ĝ(m 1 , m 2 , a) = â01 m 2 + â02 m 2 2 + â10 m 1 + â11 m 1 m

Figure A. 1

 1 Figure A.1 RMS misfit evolution for constrained (left) and non constrained (right) joint inversion of non-linear relationship test.

Figure B. 1

 1 Figure B.1 Correspondence map joint inversion results for data in Section 5 using linear and non-linear relationships. Colour codes: reference models (grey), starting model (red), correspondence map joint inversion using linear relationship (green) and correspondence map joint inversion using non-linear relationship (magenta).

Figure

1 .

 1 The observed datasets are based on the synthetic datasets of Section 4.2. The root mean square (RMS) plots in Figure C.1(a-c) show how similar the results of the first three tests are in terms of iterations needed to converge.

  Figure C.1(d) shows the effect of not giving any weight to the ĝ term, which results in a high RMS for this term. Figure C.2 shows the obtained models from these four tests. Overall, the four tests result in very similar models, with the last test showing a slight difference in the correspondence map plot compared to the rest, due to the non-inclusion of the correspondence map term in the joint inversion.

Figure C. 1

 1 Figure C.1 RMS evolution for the four weight tests described in Table C.1. The dashed black line represents the threshold when RMS = 1. The title of each subplot indicates the value of the weight of each dataset (three values are given), and the order is given as follows: index 1, the weight of the dispersion curve dataset; index 2, the weight of the magnetotelluric dataset; and index 3, the weight of the g dataset.

Figure C. 2

 2 Figure C.2 Resulting models and correspondence maps plots of joint inversion weight tests described in Table C.1. The label of each test indicates the value of the weights of each dataset (three values are given), the order is given as: index 1: weight of dispersion curves dataset, index 2: weight of the magnetotelluric dataset and index 3: weight of g dataset In black is shown the observed 'true' parameter relationship.

Chapter 3 Joint inversion methods Contents 3.1 Structural approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

	3.1.1 Cross-gradients . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
	3.1.2 Dot products of image gradients . . . . . . . . . . . . . . . . . . 36
	3.1.3 Gramian constraints . . . . . . . . . . . . . . . . . . . . . . . . . 37
	3
	Bibliography
	Bensen, G. D., Ritzwoller, M. H., Barmin, M. P., Levshin, A. L., Lin, F., Moschetti, M. P.,
	Shapiro, N. M., and Yang, Y. (2007). Processing seismic ambient noise data to ob-
	tain reliable broad-band surface wave dispersion measurements. Geophysical Journal
	International, 169(3):1239-1260.
	Cagniard, L. (1953). Basic theory of the magneto-telluric method of geophysical prospect-
	ing. Geophysics, 18(3):605-635.

.2 Petrophysical approach . . . . . . . . . . . . . . . . . . . . . . . . . . 37

	3.2.1 Correspondence Maps (CM) . . . . . . . . . . . . . . . . . . . . . 37
	3.2.2 Variation of Information (VI) . . . . . . . . . . . . . . . . . . . . 41
	3

.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

  Table6.1 gives a general overview of the acquisition parameters of the three MT surveys in the area and Figure6.1 the location of the stations for each of the three surveys. Note that only DEEP-EM and MT-Two were processed in-house, while ECOGI data were processed by[START_REF] Abdelfettah | Resistivity image under grt1-2 geothermal doublet of the rittershoffen egs project as revealed by magnetotelluric[END_REF].

	Survey	Year of	Stations	Reference	Sampling	Processing
		acquisi-	de-	stations	rate
		tion	ployed
	ECOGI	2018	10	1	512Hz	[Abdelfettah
						et al., 2019]
	DEEP-EM	2020	32	2	512Hz,	In house: Razor-
					4096Hz	back
	MT-Two	2022	6	1	512Hz,
					4096Hz

Table C . 1

 C1 Parameters of the tests for weight analysis of joint inversion using correspondence maps

		A 1 = Weight	A 2 = Weight	A 3 = Weight
	Test	(SWD)	(MT)	(Corr. Map)
	(a)	1/3	1/3	1/3
	(b)	1/2	1/4	1/4
	(c)	1/4	1/2	1/4
	(d)	1/2	1/2	0
	Abbreviations: Corr. Map: correspondence maps; MT, magnetotelluric; SWD,
	surface wave dispersion.

displaying the output of Razorback, apparent resistivity and phases, it is worth to mention an issue that occurred during this campaign that affected a random number of stations.A time difference was identified in the time series recorded during the DEEP-EM campaign, a problem related to the hardware of ADU units. This problem meant that the time stamp in the stations was not the same in all stations recording at the same time, which implied that when using reference stations for processing, common events would not be synchronized. Section A.1.1 describes the process implemented to work around this issue.

© 2022 The Authors. Geophysical Prospecting published by John Wiley & Sons Ltd on behalf of European Association of Geoscientists & Engineers., Geophysical Prospecting, 1-16

-1.51m 1 + 0.81. (23) ©

The Authors. Geophysical Prospecting published by John Wiley & Sons Ltd on behalf of European Association of Geoscientists & Engineers., Geophysical Prospecting, 1-16

Acknowledgments

I would like to start thanking my PhD supervisors Guy Marquis and Jérôme Vergne. To Guy for accepting me to start this PhD journey with you. For the patience you had while

AC K N OW L E D G E M E N T S

We are grateful to CONACyT for awarding a doctoral fellowship to MA. This work was supported by the TelluS Program of CNRS-INSU. Thanks to Jonathan Carrillo and Maximilien Lehujeur for their valuable comments on the implementation of correspondence maps and on the inversion of surface wave dispersion curves, respectively. We also thank two anonymous reviewers for their suggestions that resulted in a largely improved manuscript.

We are grateful to CONACyT for awarding a Doctoral Fellowship to MA. This work was supported by the TelluS Program of CNRS-INSU. We also thank the ECOGI/ES group (Exploitation de la Chaleur d'Origine Géothermale pour l'Industrie/Electricité de Strasbourg) for providing access to MT data and geological information. Thanks to EOST for the availability of passive seismic data (EstOF network). The EstOF dataset has been acquired in the framework of the LABEX ANR-11-LABX-0050-G-EAU-THERMIE-PROFONDE https://dx.doi. org/10.25577/2014-ESTOF. Finally, thanks to CCUS (Centre de Calcul de l'Université de Strasbourg) for their support on HPC (High Performance Computing).

DATA AVA I L A B I L I T Y S TAT E M E N T

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

Chapter 5. Towards a higher-dimensional approach rec y=3000 x=46000 z=0 file=sta460 sacformat=1 usgsformat=1 rec y=3000 x=48000 z=0 file=sta480 sacformat=1 usgsformat=1 rec y=3000 x=50000 z=0 file=sta500 sacformat=1 usgsformat=1 rec y=3000 x=52000 z=0 file=sta520 sacformat=1 usgsformat=1 rec y=3000 x=54000 z=0 file=sta540 sacformat=1 usgsformat=1 rec y=3000 x=56000 z=0 file=sta560 sacformat=1 usgsformat=1

Computation of surface-wave dispersion curves

In the previous step we computed model responses for all the source-receiver combina-Chapter 5. Towards a higher-dimensional approach Figure 5.9: Group dispersion maps for models in Figure 5.7 recovered from regionalization (Section 5. 1.3.1) of dispersion curves obtained from SW4 model responses.

Geological background

The NNE-SSW trending Upper Rhine Graben (URG), located in NE France -SW Germany, is a segment of the European Cenozoic Rift System that extends from the North Sea up to the Mediterranean [START_REF] Duringer | The complete lithostratigraphic section of the geothermal wells in rittershoffen (upper rhine graben, eastern france): a key for future geothermal wells[END_REF]. The URG is approximately 300 km long and 30 km wide, limited by the Rhenish Massif to the north and the Jura Mountains to the south. It is bounded by prominent fault systems (Figure 6.2): the Western and Eastern Rhine Graben faults) [START_REF] Bauer | Architecture, fracture system, mechanical properties and permeability structure of a fault zone in lower triassic sandstone, upper rhine graben[END_REF], [START_REF] Schumacher | Upper rhine graben: Role of preexisting structures during rift evolution[END_REF], [START_REF] Brun | Crustal-scale structure of the southern rhinegraben from ecors-dekorp seismic reflection data[END_REF], [START_REF] Duringer | The complete lithostratigraphic section of the geothermal wells in rittershoffen (upper rhine graben, eastern france): a key for future geothermal wells[END_REF]. The URG is an asymmetrical graben with thicker deposits in the north than in the south, the thickest Tertiary sediments can reach 3.5 km in the deepest parts (north of the URG).

The two profiles presented in this Chapter fall in two structurally different zones of the URG. Profile 1 is located entirely in the graben, as can be seen in Figure 6.1 while Profile 2 is located in a transition zone between the graben (east side of the line) and the Rhenan fault towards the Hochwald Horst (west side of the line). Figure 6.2: Schematic geological cross section through the Upper Rhine Graben at the latitude of Soultz-sous-Forêts and Rittershoffen geothermal plants. Reproduced from [START_REF] Duringer | The complete lithostratigraphic section of the geothermal wells in rittershoffen (upper rhine graben, eastern france): a key for future geothermal wells[END_REF] 6.2 The datasets

SWD data: EstOF project

The dispersion curves in this Thesis were acquired and processed as part of the EstOF project [Lehujeur, 2015] [ Lehujeur et al., 2018]. The EstOF survey was acquired in 2014 over 37 days. The main network was composed of 259 stations with an average spacing of 1.4 km, covering an area of about 490 sqkm. The stations consisted of 10 Hz one component (vertical) nodes with a sampling frequency of 250 Hz.

The datasets

The data was processed by Maximilien Lehujeur (details in Lehujeur [2015] and Lehujeur et al. [2018]) using standard processing steps up to regionalization [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF]. A total of 41,284 correlations between the station pairs were computed. A grid of 488 positions spaced about every 1 km was designed for regionalization. The white dots in Figure 6.1 represent these positions i.e. where group velocity curves were recovered from regionalization. Figure 6.3 shows the maps of the Rayleigh-wave group velocity as a function of period. Data consists of fundamental-mode group velocity curves for periods between 0.6 and 4.1 s. First overtone group velocity curves were provided for short periods (0.7 to 1.7 s), but they were not used here in the joint inversions here.

The first feature that comes out by looking at the data is a high-velocity anomaly toward the NW part of the area (Hochwald Horst), which can be seen too in the pseudosection of Profile 2, this will be discussed later when we analyze this profile. Figure 6.4

shows the pseudosection for Profile 1.

In this Chapter, all dispersion curves data come from the work in Lehujeur [2015] and Lehujeur et al. [2018], without additional processing. Details on the processing of the dispersion curves can be found in their work. Throughout this Chapter, the term EstOF is used to refer to this data set from Figure 6.3.

Chapter 6. Application: North Alsace

MT data: MT-Two

The MT-Two survey was acquired in March 2022 and was designed as an extension of the DEEP-EM survey in order to map the high-velocity anomaly observed in group velocity maps in the northwestern part (Figure 6.3). Here again, Welschbruch data were used as remote reference.

This area was not covered by the DEEP-EM survey and therefore it seems relevant to acquire MT data there to test whether this lateral variation is also present in the resistivity model. Six stations were deployed over four days with an acquisition sequence similar to the one used in DEEP-EM. The schedule of the acquisition along with the apparent resistivity and phase curves obtained from Razorback for all six stations can be found in Appendix A Section A.2.

Profile 1: North west -South east

This profile runs on NW-SE direction in the middle of the area covered by the EstOF survey as can be seen in Figure 6.1, it has a length of ≈15 km. For this profile we selected 10 MT stations from the ECOGI survey (Figure 6.6) with a frequency band from 1.56e-2 to 100 Hz (i.e. 0.01 to 64 s) and 16 locations of group velocity dispersion curves with frequency content from 0.24 to 1.62 Hz (i.e. 0.6 to 4.16 s). Sensitivity kernels for Rayleigh-wave group velocity dispersion curves were computed using the velocity model of well GRT-1, here reconstructed from the one presented in [START_REF] Maurer | On-going seismic monitoring of the rittershoffen egs project (alsace, france)[END_REF], which is located near the profile. The reference model consists of shear and compressional velocities built from VSP (Vertical Seismic Profile) data, sonic logs and litho-stratigraphic logs. Density was computed from a relationship built based on petrophysical data from [START_REF] Kushnir | Characterizing the physical properties of rocks from the Paleozoic to Permo-Triassic transition in the Upper Rhine Graben[END_REF], values vary from 2 to 2.65 g/cm3. Figure 6.7 shows that periods shorter than 2 s are sensitive to the first 500 m. Periods longer than 2 s are sensitive to a greater depth range down to the first kilometer. Note the lower sensitivity observed below approximately 1.5 km, even at the longest periods (4.16 s). This indicates that the surface waves might not have the sensitivity to the deeper formations. This observation has already been made by Lehujeur et al. [2018] as an explanation for the lack of correlation between their deeper velocity model and known lithology.

One-dimensional inversions

We attempt a first approach at joint inversion with 1D models to understand the relationship of the two properties, i.e. resistivity and shear velocity, before performing joint inversions in 2D. We take the locations of dispersion curves closest to the MT stations and performed separate and joint inversion for ten pairs of stations. The MT data in Figure 6.6 show that both XY and Y X apparent resistivity curves are similar and hence indicate a one-dimensional behavior up to about 10 s (0.1 Hz) and therefore the 1D inversions computed in this Section are restricted to this frequency range. Half-spaces were used as starting models for both methods and a 10 % error was used for the ĝ(m 1 , m 2 , â) func-gradients and dot product operators have a greater number of null or very low values in the joint inversion models than those from separate inversions (Figure 6.23). A similar situation occurs for the VI values: the VI term in a joint inversion objective function is to be minimized . Even though is not completely minimized, the joint inversion models have lower VIs (2.19 and 1.95) than those from a separate inversion (VI=2.62). This comparison of several joint inversion operators is a possible indicator of the synergy in the area between structural and petrophysical joint inversion approaches.

Profile 2: W-E

The second profile crosses a major discontinuity associated to the Rhenan fault between the sedimentary graben and the Vosgian massif. This discontinuity is very well identified in the resistivity model obtained from MARE2DEM. Contrary to that observed in profile 1, no singular parameter relationship (shear velocity-resistivity) for all the profile is found, suggesting different petrophysical behaviors on each side of the profile.

Public data in the area shows oil wells drilled close to the profile between [1955][1956][1957][1958][1959][1960] but none of them reach the Paleozoic granite basement ( BSS00MJSM, BSS000MJSP, BSS000MJSN). Stratigraphic descriptions show that these wells ended in the Jurassic marls at about 1100 m depth below surface. In the eastern part of the profile, the resistivity observed from the induction logs (BSS00MJSM) shows values that correspond to the inverted resistivity image below 1000 m i.e. <20 Ωm. Figure 6.25 shows an schema of the structural features in the joint inversion velocity and resistivity models. Note that even when a lateral change is observed in the velocity models, the position of the basement is not as well defined especially in the east side of the profile, as it is in the resistivity model. There is still an improvement however when comparing the model from the joint inversion to the one from the separate inversion: the high velocities in the east are reduced and pushed down by the joint inversion.

We try to correlate the major structural components from this profile with the schematic representation [START_REF] Duringer | The complete lithostratigraphic section of the geothermal wells in rittershoffen (upper rhine graben, eastern france): a key for future geothermal wells[END_REF]. Towards the east side, the thicker package of sed-Chapter 6. Application: North Alsace In summary:

• We implemented multi-physics imaging in the Upper Rhine Graben using MT and SWD data to produce resistivity and velocity models consistent with proxies measured in the laboratory.

• Profile 1 located well inside the URG shows an overall 1D structure with small lateral variations. The velocity model is improved in the deeper section when jointly inverted SWD and MT data, largely due to the higher sensitivity of MT data to thick resistors and better sensitivity at greater depths.

• Profile 2 located near the edge of the URG shows a stark contrast in petrophysical properties (and hence geological structures) between the west and east parts of the profile, likely separated by the Rhenan fault.

Chapter 7. Conclusion and Perspectives

We analyzed the importance of the choice of the polynomial degree to use in the ĝ(m 1 , m 2 , a) function for the convergence of the joint inversion. We modified the proposed correspondence maps objective function of Carrillo and Gallardo [2018] by adding weights to each data term, to account for the bias introduced when the datasets have very different sizes.

This work was published as Aquino et al. [2022b]. The full article is attached in Appendix D.

Item 1: Higher-dimensional joint inversion using CM

This part proved to be fairly challenging. Magnetotellurics practitioners have routinely used publicly-available 2D codes since the late 1980s and 3D codes from the mid 1990s, thanks to [START_REF] Weaver | A comparision of analytical and numerical results for a two-dimensional control model in electromagnetic induction -i. b-polarization calculations[END_REF], Wannamaker et al. [1987], [START_REF] Mackie | Three-dimensional magnetotelluric modelling using difference equations -theory and comparisons to integral equation solutions[END_REF] and others. It is however not the case for surface-wave dispersion curve practitioners who usually work with 1D codes (such as the ones I used here). Now there are plenty of fine reasons for this and I got no argument against the practice. I was however curious as to what would happen if we generated synthetic seismograms from a higher-dimensional velocity/density model and put them through a standard processing sequence to obtain the SWD curves and compare the responses to those obtained from 1D codes.

Sensitivity analysis of surface-wave dispersion curves of a 2D model

In Chapter 5 we built a workflow to estimate surface-wave dispersion curves using seismograms computed with a 3D modeling code. We first validated it on a one-dimensional model by comparing its results to those obtained by a 1D code. We show that it is a reliable approach to obtain group velocity maps from synthetic seismograms. We show that for 2D models, i.e. with lateral velocity variations, the group velocities obtained from our numerical approach are not equivalent to the ones we get from 1D models. It suggests revisiting the standard approach of using 1D inversions to produce 2D or 3D velocity models, at least for cases of complex subsurface.

2D joint inversion

We expose limitations on the lateral resolution of surface waves for small anomalies (compared to the wavelength) for which MT data is more sensitive. This is especially true for near-surface anomalies as those presented in this Thesis.

Our joint inversion attempts show limitations when the number of parameter relationships that can fit the data is very broad, especially if the models are too simple. This was the case for our attempts to fit a degree 2 relationship using a 3-block model, so with only 3 velocity-resistivity pairs. The joint inversion was however successful when we imposed a degree 1 relationship, which has inherently a smaller model space of parameter relationships that can fit the data with these three velocity-resistivity pairs.

Appendix A

MT processing

The MT data used in Profile 2 of Chapter 6 is composed of two acquisition campaigns: DEEP-EM and MT-Two. Below we present the apparent resistivity and phases with frequency output from Razorback and a workaround process encountered during the acquisition of DEEP-EM.

A.1 DEEP-EM processing

The calendar of the acquisition of DEEP-EM survey is described in Step 2 is a zoom in to one event, three pannels are shown, row one is H x for the station analyzed (i.e. site 113), row two is the cross correlation of station and CIEL reference station and row 3 is the cross correlation of CIEL and Welschbruch. Row two shows how the event detected in both reference stations (red star) at time t is not observed at the same time for site 113 and has a delay of 9 sec.

In order to identify the stations affected we developed a workflow which consisted in the use of the two reference stations i.e. CIEL and Welschbruch located 43 km and 72 km away from the area of the survey in order to find common events between these two and compare them with the rest of the stations. A priori the two reference stations are not affected by time shift. The process involved two steps: (i) Identify common magnetic events using the time series of the two reference stations.

(ii) Compare visually with the time series of each station in order to determine if there was or not a time delay and quantify it.

The first step (i) confirms that the reference stations have the correct time stamp and flags common events, i.e. if an event is registered in reference station 1 at time t and it is also registered in reference station 2 at time t, then both reference stations have the same time stamp and therefore common events could be flagged and then compared to the rest of the stations. This is based on the fact that the two reference stations were not affected by the time shift issue. This step is done firstly by identifying events with amplitudes greater than a certain threshold at each of the reference stations and then with the cross correlation to amplify the common events to be later compared with each of the stations. This step is illustrated in the upper graph 1) of Figure A.2.

The second step (ii) involves the visual comparison of the previously flagged common events in the reference stations to the station of study (s 1 ). This is done by computing the cross correlation of station (s 1 ) with one of the reference stations. In the case the station s 1 has no time delay, the events flagged from the reference stations in step (i) should overlay at the same time as events in station (s 1 ), otherwise there is a time shift and it has to be quantified. The second graph of Figure A.2 shows an example of how the delay is calculated between the event that was flagged in the reference stations and the station of analysis (s 1 ).

In Figure A.3 we show one example of site 202 processed with Razorback before and after time correction applied to the time series. It can be seen the more scattered the curves are when reference stations are used before the correction for time shift, compared to after time shift correction where the observations are more coherent to one another.

A.1.2 Razorback output

Below are the stations processed with Razorback using two reference stations i.e. Welschbruch and CIEL. The stations used in Profile 2 from the DEEP-EM survey are: site 309, 204, 318, 326, 331, 336. As can be seen in the schedule of the acquisition, more stations were processed, but are not displayed here since they do not fall along the profiles presented in this Thesis.

A.2 MT-Two processing

A.3 Conclusions

The work in this appendix shows the processing done for MT data using a very recent open source module in python i.e. Razorback [START_REF] Smaï | Razorback, an open source python library for robust processing of magnetotelluric data[END_REF]. We took advantage of multiple reference stations combinations that the code can handle. Together with the simplicity of using it, makes it a promising module for robust MT processing which we encourage to explore. Bibliography [START_REF] Smaï | Razorback, an open source python library for robust processing of magnetotelluric data[END_REF]. Razorback, an open source python library for robust processing of magnetotelluric data. Frontiers in Earth Science, 8.

Appendix B

Shear velocity maps in North Alsace: a comparison of different inversion approaches

In this appendix, we present two approaches of inversion of surface-wave dispersion curves that were developed as an aside work of the main scope of this Thesis and compare them with the results from Lehujeur et al. [2018] using the EstOF data set presented in Chapter 6.

B.1 Introduction

During this Thesis, we developed and tested several approaches of inversion of surfacewave dispersion curves that are worth to present as part of an aside result of the main scope of this Thesis i.e. joint inversion. Two inversion approaches were developed, i) the Occam approach (Section B.3) that has been covered in several Chapters of this Thesis and ii) a simple Monte Carlo approach (Section B.4). We compare the results from these two methods, to the shear velocity model obtained by Lehujeur et al. [2018] who followed a robust MCMC (Markov-chains-Monte-Carlo) inversion approach using the EstOF data.

It is important to note that while Lehujeur et al. [2018] used bimodal data (fundamental model and first overtone group velocity dispersion curves), the inversions performed here used only fundamental mode dispersion curves.

B.2 Reference velocity model context

As mentioned before, the models produced in this appendix are compared to a reference velocity model produced by Lehujeur et al. [2018]. Here we intend to give a brief description of the inversion approach they used for context purpose only. These models were computed using a Monte Carlo inversion approach. For each position, they inverted Rayleigh-wave group velocity dispersion curves (fundamental mode and first overtone) using 12 independent Markov chains in parallel, and finish each chain until 2000 best approaches models are kept. A total of ≈ 48000 models were produced. Next, the median of the best models was taken as the 1D shear velocity value for that pixel. Each 1D inversion was parameterized by 9 layers with upper and lower bounds for Vs, Vp/Vs and density (Table 1 in Lehujeur et al. [2018]). Shear velocity, density and Vp/Vs were inverted simultaneously in this study.

We acknowledge the robustness of this approach and here we only intend to compare the results with much simpler approaches used in this Thesis.

B.3 Occam inversion

As described in Chapter 3 an expanded version of the Occam algorithm was used to frame the joint inversion approach. Occam is a variant of non linear Gauss-Newton minimization problem [Constable et al., 1987] widely known in geophysical applications, specially in electromagnetic methods. It linearizes the functional about a starting model and accomplishes the minimization iteratively using the objective function defined in equation 4.1. It is an approach that rather than fitting the data as well as possible, the smoothest model which fits the data within a certain tolerance is sought [Constable et al., 1987] [deGroot-Hedlin and Constable, 1990]. We used an homogeneous model as initial model of 19 layers as has been done previously in mostly all inversions within this Thesis. Iterations ran until the condition of RMS<1 was reached or a maximum of 15 iterations allowed. An example for inversion at one data position can be seen in Figure B.1. The convergence is generally very fast, in this case at iteration 7 the condition was achieved.

B.4 Monte Carlo inversion: Metropolis sampling

A second inversion algorithm was developed to test the data set a side of the Occam approach. We designed a simple Metropolis sampling approach. The Metropolis algorithm is a Markov chain Monte Carlo method (MCMC), i.e. it is random and with no memory, each step depends only on the previous one. The basic idea is to perform a random walk that samples some initial probability distribution and then use a probabilistic rule to modify the walk in such a way that the modified walk samples the target distribution [Tarantola, 2005]. Here we use a probabilistic rule defined by the value of the likelihood Something interesting to note is the smoothing that occurs lateral and vertically in the models obtained using the Occam algorithm (Figure B.4). We remind that since these models are obtained using a 1D inversions at each pixel, only a vertical smooth is applied, however the models show a clear smooth laterally as well. Note that this is the core of the Occam razor, to produce the smoothest (simplest) solution that explains the data.

B.5.2 Deep part: >1.5 km

While the general structure of the shallow section is similar in the three models compared, this is not the case for the deeper section. Differences can be observed straightforward in the three models in depths below 1.5 km (Figure B.6). The models produced by Lehujeur et al. [2018] have overall higher velocities compared to those using Occam (2.0 km and 2.5 km panels in Figure B.6). One explanation for this is the influence of the strict priors set during the inversions for these depths within the approach followed by Lehujeur et al. Carlo approaches, the general lower velocity of these periods (4.2 and 2.0 s) match better to the observed data. These results mean that these two algorithms performed following the general tendency of the observed data, which was a lower velocity variation at 4.2 s, Based on the observations above, we can say that a big factor to achieve the high Appendix B. Shear velocity maps in North Alsace: a comparison of different inversion approaches velocities in the deeper sections in the models from [Lehujeur et al., 2018] is importantly guided by the priors defined. In Chapter 6 we pointed out the lower sensitivity at these depths (<1.5 km), (Figure 6.7). From the geological knowledge it is expected to have higher velocities at these depths, due to the transition to the basement It is however, not what is seen in the data at the longest periods.

B.6 Conclusions

We had access to a large database of processed surface-wave dispersion curves kindly provided by Lehujeur et al. [2018]. We took the initiative to test two in-house inversion approaches with this data and compare with the results obtained by Lehujeur et al. [2018].

The Occam algorithm is a very good approach to find a smooth solution fast, with simple initial models and without known bounds. We encourage the use this method that is an standard optimization approach in electromagnetic problems, however much less known in seismology literature.

It is straightforward the robustness of the algorithm implemented in Lehujeur et al.

[2018] that translates to more detailed models. However, it would be interesting to investigate why data wise, the longest periods (>3 s) show a consistent decrease in velocity while the expected lithology is to increase. It could certainly be due to the uncertainty of the data which has been acknowledged as difficult to estimate. Appendix C

Inversions using L-BFGS-B optimization attempts

C.1 Definition of the objective function

Multiple ways to approach the joint inversion problem exists. During this Thesis we mostly focused on the implementation of the Occam algorithm [Constable et al., 1987].

However, few attempts were made to develop an inversion using the L-BFGS-B optimization method. L-BFGS-B is a limited-memory algorithm for solving large non linear optimization problems subject to simple bounds on the variables [START_REF] Zhu | Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF]. It is an extension of the limited-memory algorithm (L-BFGS) for unconstrained optimization described in [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF], with the main improvement being the ability to deal with bounds on the variables. We do not intend to describe the algorithm, more details can be found in [START_REF] Nocedal | Large-Scale Unconstrained Optimization[END_REF].

Following the definition of the objective function as described by Moorkamp et al.

[2011] we have:

where ϕ d (m) is the data misfit term (equation C.2) and ϕ reg (m) the regularization term (equation C.3).

Here, g(m) is the synthetic data, d obs the vector of observed data, C -1 d the inverse of the data covariance matrix. α is a weight contribution, W is the smoothing matrix , m 0 is an a priori model, C M is a diagonal model covariance (similarly as the data covariance matrix)

C.2 Implementation of separate inversions

C.3 Towards a joint inversion objective function

The objective function defined in equation C.1 can be slightly modified to add a coupling term to link two or more model parameters. If we assume, two model parameters m 1 and m 2 , C.2 becomes:

The 2.8 obtained using the PW2D code of [Wannamaker et al., 1987]. Four components are output: apparent

5.15

Model response using PW2D [Wannamaker et al., 1987]