# Simultaneous joint inversion of synthetic seismic and ground penetrating radar data with petrophysical variable change

Quentin Didier, Slimane Arhab, Gaëlle Lefeuve-Mesgouez

► To cite this version:

Quentin Didier, Slimane Arhab, Gaëlle Lefeuve-Mesgouez. Simultaneous joint inversion of synthetic seismic and ground penetrating radar data with petrophysical variable change. Geophysical Journal International, 2026, 244 (1), 10.1093/gji/ggaf423. hal-05385830

HAL Id: hal-05385830

https://hal.inrae.fr/hal-05385830v1

Submitted on 27 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Simultaneous joint inversion of synthetic seismic and ground penetrating radar data with petrophysical variable change

Quentin Didier<sup>1, ∗</sup>, Slimane Arhab and Gaëlle Lefeuve-Mesgouez

Avignon Université, INRAE, UMR EMMAH, F-84000 Avignon, France. E-mail: quentin didier@univ-avignon.fr

Accepted 2025 October 13. Received 2025 September 15; in original form 2025 January 13

## SUMMARY

In this work, we address the characterization of porosity and water saturation in a synthetic model of a shallow alluvial subsurface using frequency electromagnetic and seismic data. The inversion method employs a Gauss–Newton scheme, where the Jacobian of the merged seismic and electromagnetic data is formulated with respect to the spatially heterogeneous petrophysical parameters. This is made possible by introducing realistic petrophysical relationships, which significantly reduce the number of unknowns in the inverse problem and incorporate a strong prior correlation between the information contained in both data types regarding the subsurface composition. The results obtained show that this simultaneous joint petrophysical inversion produces reconstructions with clear improvements compared to independent petrophysical inversion. Indeed, it greatly enhances the spatial resolution of subsurface mapping, as well as the quantitative estimation of porosity and saturation.

**Key words:** Ground penetrating radar; Inverse theory; Joint inversion; Numerical modelling; Seismic tomography; Waveform inversion.

## 1 INTRODUCTION

Water is among the resources most affected by climate change, as rising temperatures contribute to its scarcity, particularly in its potable form. The majority of drinking water sources on Earth are either rivers/streams or underground aquifers. In this context, it becomes crucial to develop methods for quantifying the water present in the subsoil and its exchange with surface reservoirs, to enable informed management of this vital resource. Both destructive and non-destructive methods can be used for quantifying subsoil water. Among non-destructive techniques, seismic reflection (SR) and ground penetrating radar (GPR) are promising measurement methods, along with electrical resistivity tomography (ERT) and controlled source electromagnetism (CSEM). Full waveform inversion (FWI) of seismic reflection data was first proposed by Tarantola (1984) and Mora (1987). This inversion method presents a great potential to build high resolution map of multiparameter subsurface model both in the electromagnetic (Feng et al. 2019) and seismic (Brossier et al. 2009) cases and it is therefore chosen in this work, but as part of a strategy that combines these two physics. Indeed joint inversion is a promising tool for subsoil geophysical prospecting, widely applied in other research fields, such as biomedical imaging (Song et al. 2019; Qin et al. 2021), due to its potential to enhance spatial resolution and to its versatility.

At the reservoir scale, previous studies have shown that seismic data are more sensitive to porosity, while electromagnetic data of the CSEM type are more sensitive to fluid saturation (Gao et al. 2012). Thus each physics is complementary to the other making

joint inversion a strategy of interest. As previously stated, for subsoil prospection, different kinds of geophysical data can be jointly inverted to characterize this complex natural medium through its geophysical mechanical parameters (density, speeds of compressional and shear waves), geophysical electromagnetic parameters (relative dielectric permittivity, electric conductivity) and petrophysical parameters (porosity, fluid saturation). GPR is most often chosen over ERT due to its higher reconstruction accuracy, as evidenced in the literature (Loiseau et al. 2023), and its suitability for shallow subsurface imaging. Additionally, incorporating petrophysical relations in the inversion process allowed for the correlation of geophysical parameters with each other, particularly through a joint inversion of seismic (tomography borehole or seismic refraction) and electric data (DC resistivity) (Haber & Holtzman 2013), with an application to permafrost imaging (Wagner et al. 2019). When the petrophysical relations are reliable, a better joint inversion than using a cross-gradient constraint is obtained (Abubakar et al. 2012). Cross-gradient constraint have been used in different configurations, for instance, by Gallardo & Meju (2003). Joint inversions of the FWI type have a validity domain beyond the Born approximation, which suffers from limitations, especially when multiple scattering of waves in heterogeneous media occurs (Lan et al. 2018). These approaches combine a rigorous forward model with optimization algorithms of the Gauss–Newton or quasi-Newton types (Habashy & Abubakar 2004; Nocedal & Wright 2006). For instance, in microwave tomography, an inversion using a Gauss–Newton algorithm stabilized with standard Tikhonov regularization (Tikhonov & Arsenin 1977), and constrained by a variable change,
allowed for high spatial resolution reconstructions (Didier et al. 2024).

The same approaches were deployed on a smaller scale for the characterization of shallow subsurface from seismic and GPR electromagnetic data. They were applied to the prospection of latent cavities and tunnels using a truncated-Newton algorithm and cross-gradient constraint (Feng et al. 2017), with seismic data inverted under the acoustic approximation. However, surface Love and Rayleigh waves are predominant in the near-surface seismic signal. Therefore, taking them into account determines the quality of inversions. To integrate this information, a joint inversion using petrophysical relations (Archie 1942; Biot 1951a, b; Gassmann 1951; Birchak et al. 1974; Nur et al. 1998; Garambois et al. 2002) performed on GPR data and seismic data interpreted as elastic waves efficiently reconstructed the associated geophysical and petrophysical parameters. This inversion was carried out with a conjugate gradient algorithm, exploiting the high sensitivity of seismic data to porosity and electromagnetic data to water saturation (Qin et al. 2022, 2024). It was performed sequentially at each step of the iterative process, using a separate inversion sequence for each type of physical data, which is why no mixing parameter between GPR and seismic data was introduced. In each case, petrophysical parameters were estimated based on the reconstruction of the corresponding geophysical ones.

In this work, we address the characterization of porosity and water saturation of a shallow alluvial subsurface using synthetic seismic and GPR data. The inversion method uses a Gauss–Newton scheme, where the Jacobian of the merged seismic and electromagnetic data is calculated with respect to the two sought petrophysical parameters. This approach drastically reduces the number of unknowns of the inverse problem, and simultaneously estimates geophysical and petrophysical parameters at each step of the iterative process. Additionally, we take into account the presence of an outer subsoil layer with low consolidation, which is not the case in Qin et al. (2024), using the relationship from Gassmann (1951), Biot (1951a, b), Nur et al. (1998) and Pride (2005), and integrating the conductivity of the rock matrix in addition to the conductivity of water (Glover 2010). This choice is made because of a non-negligible fraction of clay in the zone under study and some other conductive rocks. This reduces the overestimation of porosity from which other models suffer and helps to discriminate between the contributions of water and rock matrix in the evaluation of conductivity. Furthermore, the approach proposed here differs from the conventional method used by Qin et al. (2022, 2024), in that it updates at each step of the iterative process the petrophysical parameters by simultaneously processing GPR and seismic data. This simultaneous joint petrophysical inversion (SiJPI) is compared to an independent petrophysical inversion (IPI) to demonstrate the benefits of joint data processing.

Following this introduction, the paper is organized as follows: Section 2 covers the study configuration, the seismic and electromagnetic forward models, as well as an overview of the petrophysical relationships used in this work. Section 3 presents the two inversion processes, IPI and SiJPI, which are compared, particularly in Section 4 dedicated to the numerical results. Finally, the paper concludes with a summary and some prospects for future work.

# 2 STUDY CONFIGURATION AND FORWARD MODEL

## 2.1 Study configuration

This section outlines the configuration used in this study, focusing on the measurement setup and the subsurface composition. Our

objective is to achieve a highly accurate reconstruction of the subsurface geophysical and petrophysical parameters. To this end, both GPR and seismic data are used.

Even if this work is purely numerical, it is based on a real experimental test site, known as the ‘Lysimeter,’ located in Avignon on the INRAE campus, near the Durance River. This site is of great interest because of its extensive instrumentation, including piezometric measurements and numerous geotechnical studies. The lysimeter experiment is specifically designed to deepen our understanding of water infiltration processes in the soil while also serving as a benchmark for advancing non-invasive geophysical techniques. The subsurface of the plot is primarily composed of alluvial deposits, with a shallow water table situated at a depth of approximately 4 to 5 m. We created a digital twin of the experimental test site to generate realistic synthetic data. Using this model, we generated the reference data to be inverted, to which we added noise with an SNR of 25 dB. This numerical model benefits from geotechnical surveys providing detailed information on the mineral composition, porosity and water saturation of each subsurface layer. On this basis, we treated the chosen petrophysical relationships and their tuning parameters as part of the a priori information for the inverse problem. Fig. 1 illustrates the geometry and parameters of each layer, along with the multisource, multi-offset measurement configuration. The configuration under study is restricted to a 2-D one with an invariance axis along the z direction. The measurement line spans 18 m, with 1 m of spacing between two adjacent receivers, which can be either GPR antennas or geophones. For GPR data, the transmitter can be any antenna operating in emission mode, while for seismic data, a hammer strike is performed vertically on the soil surface between two successive geophones. The depth of the medium is approximately 11 m and consists of four distinct layers. The first layer is predominantly composed of clay and sand, the second consists of sandy clay with some pebbles, the third layer, which contains the water table, is made up of coarse sand and large pebbles. Finally, the fourth layer, located beneath the water table, is composed of marl.

As presented in Fig. 1, porosity decreases with depth, while water saturation increases. It reaches its maximum value of one below the water table and is lower above it. In this study, we based our approach on a Van Genuchten law (van Genuchten 1980) to establish a piecewise constant saturation profile, with a minimum value of 0.4. Although this lower limit may seem high, it is actually due to the residual water trapped in the soil’s closed pores.

## 2.2 Forward model

Two forward models, based on the finite element method, are used to describe the propagation of electromagnetic and seismic waves through the subsurface. The latter, which represents a heterogeneous medium within the investigation domain enclosed by a dashed line (see Fig. 1), is characterized by its conductivity $\sigma(\mathbf{r})$, relative dielectric permittivity $\epsilon_r(\mathbf{r})$, density $\rho(\mathbf{r})$ and Lamé parameters $\lambda(\mathbf{r})$ and $\mu(\mathbf{r})$, all of which vary with the spatial position $\mathbf{r}$. These parameters, referred to as geophysical parameters, are assumed to be frequency-independent, meaning that we do not account for any dispersion effects. This assumption is justified by the frequency ranges and their relative narrowness used for the electromagnetic and seismic waves in this study (Lavoué 2014). From now on, we choose to address each forward problem in the frequency domain, considering convention $e^{-j\omega t}$ and where j is the imaginary unit ($j^2 = -1$). This is in line with our intention to solve the inverse problem in
<figure><img src="image_2.png" /><figcaption>Figure 1. Study configuration with geophysical and petrophysical parameters of each stratum composing the medium under investigation.</figcaption></figure>

the frequency domain too. Of course, this assumption is not restric-
tive as we can always invert a time-domain signal using the Fourier
transform.

For both physics, due to the assumption of a 2-D configuration,
with an axis of invariance along the **z** direction, the geophysical
properties of the medium and the excitation source only depend
on the position vector spanning the plane of incidence **r** = **x** **x** +
**y** **y**. More specifically, we address the case of E// polarization in
electromagnetism, where the non-zero component of the electric
field aligns with the axis of invariance. Seismic waves, on the other
hand, are modelled within the framework of viscoelasticity theory
and the PS-V mode, where the displacement field lies in the plane
of incidence. Thus, two partial differential equations are obtained,
one scalar and the other vectorial:

$$
\Delta E_s^k + \frac{\omega_k^2}{c_0^2} \left[ \varepsilon_r - j \frac{\sigma}{\varepsilon_0 \omega_k} \right] E_s^k = -j \omega_k J_s \delta_{rs} \quad (1)
$$

$$
\left[ (\lambda + 2\mu) \nabla \nabla \cdot \mathbf{u}_s^l - \mu \nabla \times \nabla \times \mathbf{u}_s^l + \omega_l^2 \rho \mathbf{u}_s^l \right] + \left[ (\nabla \cdot \mathbf{u}_s^l) \nabla \lambda + 2 (\nabla \mu) \cdot (\nabla^S \mathbf{u}_s^l) \right] = \mathbf{f}_s \delta_{rs}, \quad (2)
$$

where $\delta_{rs}$ represents the Dirac delta distribution localized at point $\mathbf{r}_s$. It is used to define a source term by associating it with an electric current amplitude $J_s$ (amperes) for the electromagnetic case and a linear force density $\mathbf{f}_s$ (Newton per metre) for the mechanical case. The operators $\nabla$, $\nabla \cdot$, $\nabla \times$, $\Delta$ and $\nabla^S$ correspond, respectively, to the gradient, divergence, curl, Laplacian and symmetric gradient. Constants $c_0$ and $\varepsilon_0$ denote, respectively, the speed of light and the dielectric permittivity in vacuum. The non-zero component of the electric field, which aligns with the axis of invariance, is expressed as $E_s^k = E(\mathbf{r}, \mathbf{r}_s, \omega_k)$. This indicates that this field is evaluated at point $\mathbf{r}$ for a source located at point $\mathbf{r}_s$ and emitting at frequency $\omega_k$. The displacement field lying within the plane of incidence is formulated analogously as $\mathbf{u}_s^l = \mathbf{u}(\mathbf{r}, \mathbf{r}_s, \omega_l)$. The relationships between

the compressional and shear wave velocities and the Lamé parame-
ters ($v_p = \sqrt{\frac{\lambda+2\mu}{\rho}}$ and $v_s = \sqrt{\frac{\mu}{\rho}}$) are used to perform the inversion
on the velocities, while the Lamé parameters are employed during
the forward modelling step. Under the assumption of a viscoelastic
medium, attenuation is introduced through a non-dispersive hys-
teretic parameter $\eta$, which adds an imaginary component to the first
and second Lamé parameters as $\lambda = \lambda(1 + j\eta)$ and $\mu = \mu(1 + j\eta)$.
It is assumed to be homogeneous and known during the resolution
step of the inverse problem.

Eqs (1) and (2) are solved using the finite element method with a triangular mesh using second-order Lagrange polynomial basis functions. Since the medium is assumed to be infinite, perfectly matched layers (PML) are applied to the lower, left and right boundaries of the simulation domain to eliminate artificial reflections in both physical cases. At the upper boundary, a free-surface condition is used to model the air/ground interface for the mechanical case, while a layer of air with a PML on top is implemented for the electromagnetic case. These PMLs are based on the complex stretching coordinate method (Chew et al. 1997). The forward models are run using the Comsol Multiphysics® software.

2.3 Petrophysial relationships

The previously mentioned parameters are referred to as geophys-
ical parameters because they appear explicitly in PDEs (1,2) and
thus directly govern wave propagation in the subsurface. The petro-
physical parameters of interest in this work, water saturation $S_w$
and porosity $\phi$, also influence wave propagation by modifying the
values of the geophysical parameters. These two sets of parameters
are connected through petrophysical relationships. The following
subsections present the chosen ones to constrain the solutions of the
inverse problem.
### 2.3.1 Electromagnetic petrophysical relationships

The medium under study contains a non-negligible percentage of clay in its composition. Consequently the Glover law has been selected to relate the conductivity of the medium to the petrophysical parameters (Glover 2010), with a biphasic medium containing water and rock:

$$ \sigma = (1 - \phi)^p \sigma_s + \sigma_w \phi^m S_w^k \quad (3) $$

where the adjustment coefficients *p*, *m* and *k* are used to calibrate the relationship to the specific type of soil under study. Here, *m* and *k* refer, respectively, to the cementation and saturation exponents. This relationship takes into account the conductivity of the rock including clay *σ*<sub>s</sub> as well as the conductivity of water *σ*<sub>w</sub>.

The relative dielectric permittivity is linked to the petrophysical parameters thanks to the complex refractive index method (CRIM), (Garambois et al. 2002), as follows:

$$ \varepsilon_r = \left[ (1 - \phi) \sqrt{\varepsilon_r^s} + \phi \left( S_w \sqrt{\varepsilon_r^w} + (1 - S_w) \sqrt{\varepsilon_r^{air}} \right) \right]^2 \quad (4) $$

where $\varepsilon_r^s$ and $\varepsilon_r^w$ denote the relative dielectric permittivity of the rock and water, respectively. $\varepsilon_r^{air}$ is assumed to be equal to 1, as in a vacuum. Petrophysical relationships (3) and (4) are established empirically and therefore cannot be applied to all types of soil. However, they are applicable here, as they faithfully reproduce the behaviour of the digital twin of the experimental test site studied in this work. Moreover these relations are selected over the Waxman law because of their simplicity and good ability to take into account clay in soil, without having too much adjustment coefficients as mentioned by Glover (2010). In addition, Archie's law (Archie 1942) could have been used to link conductivity with porosity and saturation, as shown by Abubakar et al. (2012), as well as Topp's law for permittivity. However, since they are determined on pure sand, they are not well suited for the present configuration, which contains a significant amount of clay in the soil.

### 2.3.2 Mechanical petrophysical relationships

In order to link the petrophysical parameters to the mechanical ones, we use the Biot (1951a, b) and Gassmann (1951) theories with the assumption that only water and air are present in the pores of the subsurface. This leads to the following relationships:

$$ v_P = \sqrt{\frac{K_{sat} + \frac{4}{3}\mu_{sat}}{\rho_{sat}}} \qquad v_S = \sqrt{\frac{\mu_{sat}}{\rho_{sat}}} \quad (5) $$

$$ K_{sat} = (1 - \beta)K_{frame} + \beta^2 M \qquad \mu_{sat} = (1 - \beta)\mu_{frame} \quad (6) $$

$$ \frac{1}{M} = \frac{\beta - \phi}{K_{frame}} + \phi \left( \frac{S_w}{K_w} + \frac{1 - S_w}{K_{air}} \right) \quad (7) $$

$$ \rho_{sat} = (1 - \phi)\rho_s + \phi(S_w\rho_w + (1 - S_w)\rho_{air}), \quad (8) $$

where $K_{frame}$, $\mu_{frame}$ are respectively the bulk and shear modulus of the dry frame, and $\rho_s$ the density of grains. Subscript *w* refers to water and subscript air refers to air. In this approach we suppose that the medium is fully saturated by a fluid which is an equivalent of the mix between water and air. Consequently the bulk and shear modulus as well the density of the fluid saturated rock matrix are noted $K_{sat}$, $\mu_{sat}$ and $\rho_{sat}$. This formulation of geophysical parameters is based on the critical porosity theory (Nur et al. 1998), which introduces a coefficient called $\beta$ that depends on a critical porosity beyond which the soil no longer behaves as a solid but rather as a

<table><caption>**Table 1.** The parameters of the petrophysical relationships assumed to be known during the inversion process.</caption><tr><td>K<sub>s</sub></td><td>2.2 × 10<sup>10</sup> GPa</td></tr><tr><td>μ<sub>s</sub></td><td>1.6 × 10<sup>10</sup> GPa</td></tr><tr><td>ρ<sub>s</sub></td><td>2695 kg m<sup>-3</sup></td></tr><tr><td>σ<sub>s</sub></td><td>1.7 × 10<sup>-4</sup> S m<sup>-1</sup></td></tr><tr><td>ε<sub>r</sub><sup>s</sup></td><td>5</td></tr><tr><td>k</td><td>2</td></tr><tr><td>m</td><td>0.9</td></tr><tr><td>p</td><td>1.1</td></tr><tr><td>c</td><td>10</td></tr><tr><td>φ<sub>c</sub></td><td>0.4</td></tr><tr><td>ε<sub>r</sub><sup>w</sup></td><td>80</td></tr><tr><td>ε<sub>r</sub><sup>air</sup></td><td>1</td></tr><tr><td>σ<sub>w</sub></td><td>0.01 S m<sup>-1</sup></td></tr><tr><td>ρ<sub>w</sub></td><td>1000 kg m<sup>-3</sup></td></tr><tr><td>K<sub>w</sub></td><td>2.3 × 10<sup>9</sup> GPa</td></tr><tr><td>K<sub>air</sub></td><td>1 × 10<sup>5</sup> GPa</td></tr></table>

suspension:

$$ \beta = \begin{cases} \frac{\phi}{\phi_c} & \text{si } 0 \le \phi \le \phi_c \\ 1 & \text{si } \phi > \phi_c \end{cases} \quad (9) $$

Here, the critical porosity $\phi_c$ is typically set to 0.4. Furthermore, the first few centimetres of a crop field are not fully consolidated. This may be taken into account by introducing a consolidation factor *c* as proposed by Pride (2005):

$$ K_{frame} = K_s \frac{1 - \phi}{1 + c\phi} \quad (10) $$

$$ \mu_{frame} = \mu_s \frac{1 - \phi}{1 + \frac{3}{2}c\phi}, \quad (11) $$

where $K_s$ and $\mu_s$ are the bulk and shear moduli of the grains. The combination of all these relationships establishes a link between geophysical and petrophysical parameters. This link depends on the knowledge of several adjustment parameters and on the mechanical and electromagnetic properties of the two liquid phases and the solid phase that make up the porous medium. Their values, provided in Table 1, are chosen based on geotechnical studies conducted at the experimental test site, and are assumed to remain constant throughout the inversion process.

## 3 INVERSION SCHEME

To reconstruct both petrophysical and geophysical parameters, FWI is applied to the observed data. This method minimizes a cost function that measures the misfit between observed and simulated data at each iteration. Two approaches are compared: the *independent petrophysical inversion (IPI)* scheme, which involves inverting mechanical and electromagnetic data separately to reconstruct geophysical parameters. Afterward, a nearest neighbour algorithm maps geophysical to petrophysical parameters; and the *SiJPI* scheme, which simultaneously inverts both mechanical and electromagnetic data. During this process, the petrophysical parameters are first estimated iteratively and then used at each step to update the geophysical ones.

### 3.1 Independent petrophysical inversion (IPI)

We begin by separately inverting the seismic and electromagnetic data using a Gauss–Newton algorithm. Being a second-order optimization method, Gauss–Newton offers a significant advantage
over first-order approaches such as the conjugate gradient method: it updates the medium parameters with a richer descent direction that incorporates additional information from the Jacobian, leading to more accurate reconstructions (Nocedal & Wright 2006). To do so, let us start with the definition of the two independent cost functionals:

$$ \mathcal{F}_E(\mathbf{m}_E) = \frac{\|E^{\text{obs}} - E^{\text{sim}}(\mathbf{m}_E)\|_2^2}{\|E^{\text{obs}}\|_2^2}, \quad \mathbf{m}_E = [\varepsilon_r, \sigma]^T \quad (12) $$

$$ \mathcal{F}_M(\mathbf{m}_M) = \frac{\|\mathbf{u}^{\text{obs}} - \mathbf{u}^{\text{sim}}(\mathbf{m}_M)\|_2^2}{\|\mathbf{u}^{\text{obs}}\|_2^2}, \quad \mathbf{m}_M = [\rho, v_P, v_S]^T \quad (13) $$

where $\mathcal{F}_E$ and $\mathcal{F}_M$ are the cost functional for the electromagnetic part and the seismic part, respectively. $E^{\text{sim}}$ represents the simulated electric field and $E^{\text{obs}}$ the observed one, while $\mathbf{u}^{\text{sim}}$ is the simulated displacement field and $\mathbf{u}^{\text{obs}}$ the observed one. The two model vectors, $\mathbf{m}_E$ and $\mathbf{m}_M$, contain the geophysical electromagnetic and seismic parameters, respectively. The iterative inversion procedure produces the calculation at each inversion step $n$ of a small correction of the parameter such as $\mathbf{m}^{n+1} = \mathbf{m}^n + \delta\mathbf{m}^n$ and $\mathcal{F}(\mathbf{m}^{n+1}) < \mathcal{F}(\mathbf{m}^n)$. In order to calculate these corrections, we solve the following linear system:

$$ \left\{ \begin{aligned} & \left[ E^{\text{obs}}(\mathbf{r}_d, \mathbf{r}_s, \omega_k) - E_n^{\text{sim}}(\mathbf{r}_d, \mathbf{r}_s, \omega_k) \right] J_s = \\ & \quad j \omega_k \varepsilon_0 \int_{\Omega} \mathcal{J}_{\varepsilon_r}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k) \delta\varepsilon_r^n(\mathbf{r}) \, d\mathbf{r} \\ & + \int_{\Omega} \mathcal{J}_{\sigma}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k) \delta\sigma^n(\mathbf{r}) \, d\mathbf{r} \\ & \quad \text{written in the compact form as:} \\ & \quad \delta E^n(\mathbf{r}_d, \mathbf{r}_s, \omega_k) J_s = (\mathbf{J}_{\varepsilon_r}^n \mathbf{J}_{\sigma}^n) \left( \frac{\delta\varepsilon_r^n}{\delta\sigma^n} \right), \end{aligned} \right. \quad (14) $$

with $\delta E^n(\mathbf{r}_d, \mathbf{r}_s, \omega_k) = E^{\text{obs}}(\mathbf{r}_d, \mathbf{r}_s, \omega_k) - E_n^{\text{sim}}(\mathbf{r}_d, \mathbf{r}_s, \omega_k)$, the misfit between measured and simulated data at the current iteration for a source located at $\mathbf{r}_s$ and a receiver at $\mathbf{r}_d$. Here $\mathbf{J}_{\varepsilon_r}^n$ and $\mathbf{J}_{\sigma}^n$ correspond to the Fréchet derivatives of the electric field with respect to $\varepsilon_r$ and $\sigma$. The expressions for $\mathcal{J}_{\varepsilon_r}^n$ and $\mathcal{J}_{\sigma}^n$ are detailed in Appendix A. Eq. (14) is derived using a formulation based on the generalized reciprocity theorem (Roger 1983; Arhab et al. 2018), which is equivalent to the well-known adjoint method (Norton 1999). The same process can be applied to seismic data, which leads to:

$$ \left\{ \begin{aligned} & \left[ \mathbf{u}^{\text{obs}}(\mathbf{r}_d, \mathbf{r}_s, \omega_k) - \mathbf{u}_n^{\text{sim}}(\mathbf{r}_d, \mathbf{r}_s, \omega_k) \right] \cdot \mathbf{f}_s = \\ & \quad \int_{\Omega} \mathcal{J}_{\rho}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k) \delta\rho^n(\mathbf{r}) \, d\mathbf{r} + \int_{\Omega} \mathcal{J}_{v_P}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k) \delta v_P^n(\mathbf{r}) \, d\mathbf{r} \\ & + \int_{\Omega} \mathcal{J}_{v_S}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k) \delta v_S^n(\mathbf{r}) \, d\mathbf{r} \\ & \quad \text{written in the compact form as:} \\ & \quad \delta \mathbf{u}^n(\mathbf{r}_d, \mathbf{r}_s, \omega_k) \cdot \mathbf{f}_s = (\mathbf{J}_{\rho}^n \mathbf{J}_{v_P}^n \mathbf{J}_{v_S}^n) \left( \frac{\delta \rho^n}{\delta v_P^n} \right), \end{aligned} \right. \quad (15) $$

with $\delta \mathbf{u}^n(\mathbf{r}_d, \mathbf{r}_s, \omega_k) = \mathbf{u}^{\text{obs}}(\mathbf{r}_d, \mathbf{r}_s, \omega_k) - \mathbf{u}_n^{\text{sim}}(\mathbf{r}_d, \mathbf{r}_s, \omega_k)$, the misfit between measured and simulated displacement field at the current iteration. Here $\mathbf{J}_{\rho}^n, \mathbf{J}_{v_P}^n$ and $\mathbf{J}_{v_S}^n$ correspond, respectively, to the Fréchet derivatives of the displacement field with respect to density, compressional and shear wave speeds. Expressions for $\mathcal{J}_{\rho}^n, \mathcal{J}_{v_P}^n$ and $\mathcal{J}_{v_S}^n$ are detailed in Appendix A. Similarly to eq. (14), eq. (15) is derived using the same formulation (Tromp et al. 2005; Arhab et al. 2017).

The matrix systems in eqs (14) and (15) are ill-conditioned and cannot be solved as they stand. This problem is amplified in (14) by the variations in permittivity and conductivity, which are not of the

same scale. This difficulty is addressed by using a normalization factor called $\sigma_0$ as proposed by Lavoué et al. (2015). To tackle the ill-posed nature of the problem, a standard Tikhonov regularization (Tikhonov & Arsenin 1977) is used. The two linear systems can then be rewritten in the following forms:

$$ \mathbf{J}_E^{n\dagger} \delta E^n(\mathbf{r}_d, \mathbf{r}_s, \omega_k) J_s = \left[ \mathbf{J}_E^{n\dagger} \mathbf{J}_E^n + \begin{pmatrix} \beta_{\varepsilon_r} & 0 \\ 0 & \beta_{\sigma} \end{pmatrix} \right] \left( \frac{\delta \varepsilon_r^n}{\delta \sigma_0^n} \right) \quad (16) $$

$$ \mathbf{J}_M^n [\delta \mathbf{u}^n(\mathbf{r}_d, \mathbf{r}_s, \omega_k) \cdot \mathbf{f}_s] = \left[ \mathbf{J}_M^n \mathbf{J}_M^n + \begin{pmatrix} \beta_\rho & 0 & 0 \\ 0 & \beta_{v_P} & 0 \\ 0 & 0 & \beta_{v_S} \end{pmatrix} \right] \left( \frac{\delta \rho^n}{\delta v_P^n} \right), \quad (17) $$

where $\mathbf{J}_E^n = (\mathbf{J}_{\varepsilon_r}^n \sigma_0 \delta \rho^n)$ and $\mathbf{J}_M^n = (\mathbf{J}_{\rho}^n \mathbf{J}_{v_P}^n \mathbf{J}_{v_S}^n)$ are the Jacobian matrices of the cost functionals $\mathcal{F}_E$ and $\mathcal{F}_M$, respectively. The normalization factor for conductivity is constructed as $\sigma_0 = \varepsilon_0 \omega_k$ when only one frequency is involved, and $\sigma_0 = \frac{\varepsilon_0}{N_\omega} \sum_{k=1}^{N_\omega} \omega_k$ when $N_\omega$ frequencies are considered. Here, the parameters $\beta_{\varepsilon_r}, \beta_{\sigma}, \beta_\rho, \beta_{v_P}$ and $\beta_{v_S}$ are the Tikhonov regularization parameters for each class of geophysical parameters. Consequently, the entire linear system to be solved is formed by concatenating for each physics the different Fréchet derivatives computed for each pair of excitation and measurement points and for each frequency. The same concatenation is applied to the differences between simulated and measured fields. Frequencies are used in successive packets of two during the inversion process, progressing from low to high values, which has been shown to be an effective strategy in Sirgue & Pratt (2004) and Didier et al. (2022).

Consequently, due to all these assumptions, the two independent linear systems to solve (16) and (17) can be extremely large in terms of RAM usage. Therefore, constructing these normal equations in their entirety is impractical. To address this issue, a conjugate gradient squared (CGS) algorithm (Sonneveld 1989) is used to compute small corrections without fully building these systems. After performing the two separate inversions, we apply a nearest neighbours algorithm (Taunk et al. 2019) to retrieve the petrophysical parameters from the geophysical ones, by exploiting electromagnetic and mechanical petrophysical relationships. This entire process is summarized in Fig. 2, where $\varepsilon_r^0, \sigma^0, \rho^0, v_P^0$ and $v_S^0$ are the initial guess of the geophysical parameters, with $\mathcal{F}_E^0$ and $\mathcal{F}_M^0$ the corresponding values of the cost functionals.

## 3.2 Simultaneous joint petrophysical inversion (SiJPI)

The second approach involves coupling the two inversions by summing the two cost functions and then applying a variable change to directly retrieve the petrophysical parameters (Wagner et al. 2019). In this way, a strong correlation is established between mechanical and electromagnetic data while reducing the number of unknowns, as the medium is reconstructed solely based on the two petrophysical parameters. We refer to this approach as SiJPI, and the cost function to minimize is written as follows:

$$ \mathcal{F}_{\text{SiJPI}}(\mathbf{m}_P) = \mathcal{F}_E(\mathbf{m}_P) + \nu \mathcal{F}_M(\mathbf{m}_P), \quad \mathbf{m}_P = [\phi, S_w]^T. \quad (18) $$

It is necessary to normalize each cost function so that they have the same influence on the global one. To do so, the weighting parameter $\nu$ is used to balance the contributions of the two physics, and its value is determined at the first iteration, after which it remains constant throughout the entire inversion process.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** IPI workflow for inversion with: $\frac{\partial \mathcal{F}}{\partial m}$ the Fréchet derivative of the cost function $\mathcal{F}$ with respect to parameter $m$ and $\gamma$ the stopping criterion.</figcaption></figure>

The minimization of cost functional (18) requires the computation of the Fréchet derivatives of the electric and displacement fields with respect to the petrophysical parameters. These new derivatives are obtained by utilizing the Fréchet derivatives provided in eqs (14) and (15) with respect to the geophysical parameters and by applying the chain rule, which involves differentiating the geophysical parameters with respect to porosity and saturation through the petrophysical relationships. Consequently, they are expressed in the following forms:

$$ \mathbf{J}_{E, \phi}^{n} = \mathbf{J}_{\varepsilon_{r}}^{n} \frac{\partial \varepsilon_{r}}{\partial \phi} + \mathbf{J}_{\sigma}^{n} \frac{\partial \sigma}{\partial \phi} \quad (19) $$

$$ \mathbf{J}_{E, S_w}^{n} = \mathbf{J}_{\varepsilon_r}^{n} \frac{\partial \varepsilon_r}{\partial S_w} + \mathbf{J}_{\sigma}^{n} \frac{\partial \sigma}{\partial S_w} \quad (20) $$

$$ \mathbf{J}_{M, \phi}^{n} = \mathbf{J}_{\rho}^{n} \frac{\partial \rho}{\partial \phi} + \mathbf{J}_{v_P}^{n} \frac{\partial v_P}{\partial \phi} + \mathbf{J}_{v_S}^{n} \frac{\partial v_S}{\partial \phi} \quad (21) $$

$$ \mathbf{J}_{M, S_w}^{n} = \mathbf{J}_{\rho}^{n} \frac{\partial \rho}{\partial S_w} + \mathbf{J}_{v_P}^{n} \frac{\partial v_P}{\partial S_w} + \mathbf{J}_{v_S}^{n} \frac{\partial v_S}{\partial S_w} \quad (22) $$

These four Fréchet derivatives are involved in the construction of the gradient $G_{\text{petro}}^{E}$ and Hessian $H_{\text{petro}}^{E}$ of the electrical data with respect to the petrophysical parameters. Similarly an equivalent methodology is conducted for the displacement field data, with gradient $G_{\text{petro}}^{M}$ and Hessian $H_{\text{petro}}^{M}$. All these expressions are detailed in Appendix B. In a compact form, the linearized system to be solved at each step of the inversion process is expressed as follows:

$$ \begin{aligned} & \left[ H_{\text{petro}}^{E} + \nu H_{\text{petro}}^{M} + \begin{pmatrix} \beta_{\phi} & 0 \\ 0 & \beta_{S_w} \end{pmatrix} \right] \begin{pmatrix} \delta\phi \\ \delta S_w \end{pmatrix} \\ & = G_{\text{petro}}^{E} + \nu G_{\text{petro}}^{M}. \end{aligned} \quad (23) $$

As can be observed, this type of inversion introduces strong coupling between the two physics during the inversion process because the problem is formulated for both physics in terms of the same petrophysical parameters. The Tikhonov parameters $\beta_{\phi}$ and $\beta_{S_w}$ regularize the inversions of the petrophysical parameters by constraining their respective norms, while $\nu$ balances the contributions of the two physics. As shown in the SiJPI workflow in Fig. 3, after updating the petrophysical parameters at each step of the iterative process, the geophysical parameters are calculated using the petrophysical relationships.

# 4 NUMERICAL RESULTS

In this section, we present two different series of results. The first one concerns a comparison between the IPI and SiJPI strategies for the reconstruction of the physical parameter mapping applied to the synthetic medium presented in Fig. 1. We demonstrate an improvement in reconstruction in both qualitative and quantitative terms. The second series is a time-lapse result, showing the ability of our inverse approach to effectively reconstruct the water table depth and porosity over time. Since a synthetic medium is studied, quality indicators of the reconstruction can be computed, such as normalized residual error (NRE) and correlation coefficient (CC) which were selected in this work. The NRE is linked to the model error at iteration $n$ but normalized with the initial error. To compute the misfit between the two models, we use a L2 norm. The CC measures the likelihood between the reconstruction and the true model.

## 4.1 Inversion in a fixed configuration

We choose a frequency domain FWI instead of a temporal FWI for both its efficiency and memory management. The choice of frequencies to be used during the inversion process has to be done carefully. The selection is based on the criterion presented by Sirgue & Pratt (2004) and the characteristics of the experimental tools. Since geophones have a common frequency bandwidth ranging from 4.5 to 250 Hz, the set of frequencies for seismic inversion process is chosen as [2, 3, 10, 15, 25, 50] Hz. Similarly, GPR antennas are designed for a central frequency. For geophysical prospection, a GPR antenna with a central frequency of 50 MHz can be used. Therefore, the set of frequencies for electromagnetic inversion process is chosen as [4, 8, 11, 17, 22, 26] MHz. Higher frequencies are not necessary in this configuration because the medium does not contain very fine details and they are also affected by a damping effect. The seismic data inversion process begins with a low frequency (2 Hz), and the initial model is perfectly homogeneous, with a porosity value of 0.35 and a water saturation value of 0.4. These petrophysical initial parameters lead to an initial homogeneous medium with a permittivity of 13, a conductivity of 1 mS m⁻¹, a density of 2200 kg m⁻³, a $v_P$ of 800 m s⁻¹ and a $v_S$ of 300 m s⁻¹. For field experimental data, it would be more efficient to start with a gradient model or preliminary results obtained through traveltime inversion or other approximate inverse models to accelerate the convergence of the FWI process and prevent the solution from becoming trapped in local minima. Starting inversion with low frequency is crucial as well as having a good starting model in order to prevent FWI from converging into local minima. In this case, using a starting frequency of 10 Hz in mechanics and 8 MHz in electromagnetism would lead the FWI to be trapped in a local minimum and thus to a reconstructed medium which would not resemble to the real one. Moreover as shown in
<figure><img src="image_4.png" /><figcaption>**Figure 3.** SiJPI workflow for inversion with: $\frac{\partial \mathcal{F}}{\partial m}$ the Fréchet derivative of the cost function $\mathcal{F}$ with respect to parameter $m$ and $\gamma$ the stopping criterion.</figcaption></figure>

the previous section, we use regularization to prevent our algorithm from becoming unstable and to deal with the ill-conditionness of our problem. These prerequisites are essential to obtain a final solution well representing the targeted medium.

Concerning the present work, Fig. 4 shows the inversion results for the reconstruction of the geophysical parameters obtained using the IPI and SiJPI inversion schemes. The first column presents the true/target medium. The data used for the inversion process were generated with this medium. Nevertheless, to avoid the inverse crime as much as possible (Wirgin 2004), noise was added, with a signal-to-noise ratio (SNR) of 25 dB, and the mesh generating the data was different from the mesh used for the inversion. The second column of Fig. 4 shows the results obtained separately for each physics using the IPI strategy. The last column presents the reconstructed geophysical parameters obtained from the simultaneous inversion of electromagnetic and seismic data using the SiJPI strategy.

From a qualitative point of view, the different layers of the subsoil can be distinguished in the GPR reconstructions making the first conclusion quite positive for both strategies. Note that the SiJPI

approach produces better results, as the boundaries of each layer are more accurately retrieved, and the values of the parameters are also well estimated. However, this conclusion no longer applies to the seismic parameters estimation. Indeed the IPI strategy does not yield satisfactory results, and it is not possible to retrieve any of the subsoil layers, except for an approximate water table. Moreover, results are really bad for the density parameter. This can be explained because seismic data are not very sensitive to this parameter (Tarantola 1986). The use of the SiJPI strategy produces significantly better results for the $v_S$ reconstruction, which leads to the visualization of all the layers of the subsoil as well as a reasonable estimation of the velocity values. To a lesser extent, the SiJPI strategy also improves the reconstruction of $v_P$ and $\rho$. In general, the SiJPI strategy enhances the results for both physics, allowing for a better definition of each layer of the subsoil and a more accurate estimation of the spatial distribution of the parameter values.

Table 2 yields a quantitative analysis. A significant improvement is obtained for the reconstruction of the density values when employing the SiJPI strategy. The correlation coefficient (CC) is
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Maps of the reconstructed geophysical parameters, with the synthetic medium (True Model) maps used to generate the observed data shown in the left column, the maps reconstructed using the iterative IPI scheme in the middle column, and those reconstructed using the SiJPI scheme in the right column.</figcaption></figure>

<table><caption>**Table 2.** Quantitative results for both strategies for the geophysical parameters reconstruction.</caption><thead><tr><th rowspan="2"></th><th colspan="2">IPI</th><th colspan="2">SiJPI</th></tr><tr><th>CC</th><th>NRE</th><th>CC</th><th>NRE</th></tr></thead><tbody><tr><td>ρ</td><td>56.7 per cent</td><td>0.75</td><td>87.9 per cent</td><td>0.40</td></tr><tr><td>v<sub>P</sub></td><td>89.39 per cent</td><td>0.36</td><td>86.1 per cent</td><td>0.54</td></tr><tr><td>v<sub>S</sub></td><td>88.3 per cent</td><td>0.27</td><td>87.9 per cent</td><td>0.45</td></tr><tr><td>ε<sub>r</sub></td><td>90.0 per cent</td><td>0.66</td><td>83.2 per cent</td><td>0.35</td></tr><tr><td>σ</td><td>79.8 per cent</td><td>0.59</td><td>82.1 per cent</td><td>0.38</td></tr></tbody></table>

improved of about 30 per cent, which is substantial. The normalized residual errors (NRE) are lower for the SiJPI reconstructions of the GPR parameters and density, while they are higher for the seismic velocities. In fact, the gains achieved in the overall reconstruction lead to a reduction of the reconstruction precision of some geophysical parameters. The main objective of this approach is not to reconstruct the geophysical parameters, but rather to reconstruct the petrophysical ones. As a result, porosity and water saturation are the key characteristics of the SiJPI strategy.

The reconstruction of petrophysical parameters is presented in Fig. 5. The same format as Fig. 4 is used, with the first column corresponding to the true/target model and the other two columns representing models reconstructed with the IPI and SiJPI strategies. The IPI strategy does not yield accurately reconstructions of the porosity of the subsoil. In contrast, the SiJPI strategy shows a quite good mapping, as every layer of the subsoil can be distinguished easily. An accurate estimation of porosity is of real importance, as

it is a crucial parameter for estimating the water content in the soil and for assessing its compaction.

Concerning water saturation both reconstructions are satisfactory. The distinction between the saturated and the unsaturated zones is clear. The water table can be computed accurately with each strategy. However, the SiJPI methodology leads to a better reconstruction, with a more precise identification of the boundary between the water table and the unsaturated zone. The boundary between the first two layers is also more accurately determined. Nevertheless the bottom of the model is not well reconstructed in either strategy. This issue arises from the limited aperture of the acquisition setup. Moreover, electromagnetic waves have limited depth penetration due to the skin effect. Therefore, these poor reconstructions were expected regardless of the selected strategy.

Table 3 presents indicators from which a quantitative analysis can be proposed. The results highlight a significant improvement in the reconstruction of the porosity parameter with the SiJPI strategy. The associated correlation coefficient has increased by 32 per cent, and the normalized residual error has decreased by 0.4 compared to the IPI strategy. This represents a substantial improvement. Regarding the water saturation, the IPI methodology shows a slightly better indicator. However, the decrease in precision for the estimation of the water saturation with the SiJPI strategy is really low compared to the improvement obtained for the porosity. The SiJPI methodology enables better identification of the water content, as it is determined by the product of porosity and water saturation. The trade-off between these two parameters in relation to the data can explain that the quality of reconstruction for porosity improves in SiJPI while there is a slight decrease for water saturation.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Maps of the reconstructed petrophysical parameters, with the synthetic medium (True Model) maps used to generate the observed data shown in the left column, the maps reconstructed using the iterative IPI scheme in the middle column, and those reconstructed using the SiJPI scheme in the right column.</figcaption></figure>

<table><caption>**Table 3.** Quantitative results for both strategies for the petrophysical parameters reconstruction.</caption><thead><tr><th rowspan="2"></th><th colspan="2">IPI</th><th colspan="2">SiJPI</th></tr><tr><th>CC</th><th>NRE</th><th>CC</th><th>NRE</th></tr></thead><tbody><tr><td>φ</td><td>57.3 per cent</td><td>0.75</td><td>89.3 per cent</td><td>0.34</td></tr><tr><td>S<sub>w</sub></td><td>89.3 per cent</td><td>0.36</td><td>84.2 per cent</td><td>0.34</td></tr></tbody></table>

<table><caption>**Table 4.** Result on the convergency speed and final misfit of the two strategies.</caption><thead><tr><th>Strategy</th><th>IPI</th><th>SiJPI</th></tr></thead><tbody><tr><td>Number of iteration</td><td>110</td><td>45</td></tr><tr><td>Final misfit</td><td>-8 dB</td><td>-12 dB</td></tr></tbody></table>

Specifically, the SiJPI strategy leverages the strengths of each physical measurement, compensates for their weaknesses and thus produces improved reconstructions. On one hand, while the GPR tool is limited by its depth penetration, the seismic measurements yield depth information on the medium. On the other hand, seismic data are not sensitive to the near-surface content, whereas the GPR data are specifically suited for getting information on this part. Therefore, coupling the two inversion processes allows us to take into account more information about the whole medium than using either physical method separately.

Additionally, the sensitivity of geophysical parameters to petrophysical parameters through petrophysical relationships is of importance. While GPR parameters are sensitive to both petrophysical parameters but show higher sensitivity to water saturation, seismic parameters are not sensitive to water saturation but are highly sensitive to porosity, as presented in Appendix C. Thus, coupling both physical inversions yields better results by using the strengths of each physical measurement to mitigate their individual weaknesses. The SiJPI strategy also reduces the number of unknowns in the inverse problem to the estimation of only the petrophysical parameters, while both physical data sets are simultaneously inverted.

Moreover this strategy leads to a notable convergence rate and also a better fit of the data as shown in Table 4. Indeed, we can see that the SiJPI inversion scheme tends toward a final discrepancy of -12 dB

after 45 iterations, whereas the IPI scheme only converges after 110 iterations, reaching a higher value of -8 dB. Even though in both cases a Gauss-Newton algorithm with Tikhonov regularization was implemented, the fusion of mechanical and electromagnetic data in the SiJPI gives the cost function a favourable structure, with stronger descent directions that enable faster convergence toward the minimum of the cost function. Finally, note that the matrix sizes used during inversion are about 3 GB for SiJPI inversion and 6 GB for IPI. It is also important to note that all the results presented in this work were obtained on a computer with 128 GB of RAM and an AMD Ryzen Threadripper PRO 3955WX 16-Core CPU. The total inversion time for the SiJPI strategy is about 6 hr instead of 20 hr for IPI.

To conclude, the use of the SiJPI strategy produces better results than the use of IPI strategy, more specifically for the reconstruction of petrophysical parameters. The disadvantage of IPI stems from the independent process, with no information exchange. Information is shared only when the nearest neighbours algorithm is applied after the inversion process to retrieve petrophysical parameters from the geophysical ones. This step is not the source of the reduced resolution and precision of the reconstruction, as the nearest neighbours algorithm is a very powerful process ensuring the retrieved petrophysical parameters best fit to the geophysical ones. Next step is to ensure that the SiJPI strategy is not medium-dependent. To do so a type of time-lapse inversion using an evolving medium is conducted in the next section.

## 4.2 Time lapse inversion with variable water table depth

In this part, we perform a time-lapse-like inversion with a varying water table level. The porosity is supposed not to evolve over time, which is a reasonable assumption if there is no anthropogenic mechanical action. In contrast, the water table level varies with respect to time due to rainfall or drought. Therefore, two configurations, called ‘high water table’ and ‘low water table’, are considered. The water table is situated at approximately -2 m for the higher one and -4.5 m for the lower one.

Figs 6 and 7 show the reconstruction of petrophysical parameters for the two configurations using the SiJPI strategy. The water table level is retrieved precisely in both cases. The low resolution at the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Maps of the reconstructed petrophysical parameters in the case of a lower water table, with the synthetic medium (True Model) maps used to generate the observed data shown in the left column and the maps reconstructed using the iterative SiJPI scheme shown in the right column.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Maps of the reconstructed petrophysical parameters in the case of a higher water table, with the synthetic medium (True Model) maps used to generate the observed data shown in the left column and the maps reconstructed using the iterative SiJPI scheme shown in the right column.</figcaption></figure>

bottom of the medium, due to the limited aperture of the acquisition setup, can be observed. Although the porosity reconstructions are satisfactory, they tend to follow the water saturation level. However, the wavering boundaries are not very well reconstructed. Overall, this strategy enables to obtain a very good estimation of the water saturation in the subsoil. The water table level as well as the boundaries between unsaturated layers are accurately estimated. As a consequence we can track an infiltration front within the subsoil. These time-lapse inversions are independent: no information exchange occurs between the higher and lower configurations. One way to further improve the inversion is by coupling time-lapse inversions with the assumption of consistent porosity across each configuration. This approach would reduce the number of unknowns in the inverse problem, resulting in a more accurate estimation of water saturation.

Table 5 present quantitative results that confirm the ability of the SiJPI method to reconstruct the near subsurface with high precision for both configurations. The porosity reconstruction shows a very high correlation coefficient (CC) with a low normalized residual error (NRE), indicating good reconstruction quality. Concerning the water saturation, the quality of reconstruction is slightly decreased with the depth of the water table. Since only GPR data

<table><caption>**Table 5.** Quantitative results for both configurations for the petrophysical parameters reconstructions.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Higher</th><th colspan="2">Lower</th></tr><tr><th>CC</th><th>NRE</th><th>CC</th><th>NRE</th></tr></thead><tbody><tr><td>φ</td><td>92.9 per cent</td><td>0.35</td><td>92.6 per cent</td><td>0.34</td></tr><tr><td>S<sub>w</sub></td><td>85.3 per cent</td><td>0.35</td><td>81.0 per cent</td><td>0.35</td></tr></tbody></table>

are sensitive to water saturation, this result is attributed to the skin effect of the soil which prevents the electromagnetic waves from reaching deeper positions in the medium. However, the results for the low water table level are after all quite good. Thus the water content can be computed with a high degree of confidence in both cases.

To summarize, the SiJPI strategy is a good approach to reconstruct petrophysical parameters in various situations, especially when the water table level changes. The results remain satisfactory, even though they are slightly worse when the water table is very low. The skin effect of the soil highlights the need for low-frequency electromagnetic waves to explore the deeper parts of the model. Seismic data are essential for retrieving the soil porosity, as they are primarily sensitive to this parameter.
## 5 CONCLUSION

In conclusion, the research work presented in this article has underscored the effectiveness of the SiJPI strategy in reconstructing geophysical and petrophysical parameters of the subsurface, in the context of geophysical applications. By leveraging a coupled inversion approach, the SiJPI methodology has provided significant improvements over traditional methods, notably in the accurate reconstruction of petrophysical parameters such as porosity and water saturation. Moreover, our proposed inverse approach successfully retrieved the spatial distribution of the water table, even in a time-lapse configuration. The ability to effectively integrate data from different geophysical methods enhances our understanding of subsurface characteristics, and can be of interest to facilitate more informed decision-making in various applications, including environmental monitoring and resource management.

Future investigations will focus on refining the inversion process within the framework of time-lapse analyses to track changes in subsurface conditions, further enhancing the robustness of the developed models. A natural next step in this work will be the application of this method to experimental field data, and to do so, the SiJPI inverse model will need to be adapted to account for the 3-D geometric nature of both the excitation source and the subsurface. It should also take into account noise, both intrinsic to the measurement devices and extrinsic, such as errors in the positioning of the detectors. Among our current work, we are developing a method for characterizing the seismic source, which accurate knowledge is essential for the proper functioning of the joint inversion approach.

## DATA AVAILABILITY

The software and data that support the findings of this study are available from SA (e-mail: slimane.arhab@univ-avignon.fr).

## REFERENCES

Abubakar, A., Gao, G., Habashy, T. & Liu, L., 2012. Joint inversion approaches for geophysical electromagnetic and elastic full-waveform data, *Inverse Probl.*, **28**, 055016.

Arhab, S., Lefeuve-Mesgouez, G. & Mesgouez, A., 2017. Principe de réciprocité appliqué au calcul de la dérivée de Fréchet de l'opérateur non linéaire de propagation d'ondes dans les milieux élastiques et viscoélastiques, *23ème Congrès Français de Mécanique*.

Arhab, S., Anagnostou, D. & Maminirina, J., 2018. High-order functional derivatives of the scattered field according to the permittivity-contrast function, *Wave Motion*, **83**, 67–69.

Archie, G., 1992. The electrical resistivity log as an aid in determining some reservoir characteristics, *Trans. AIME*, **146**, 54–62.

Biot, M., 1956. Theory of elastic waves in a fluid-saturated porous solid. I. Low frequency range, *J. acoust. Soc. Am.*, **28**, 168–178.

Biot, M., 1956. Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range, *J. acoust. Soc. Am.*, **28**, 179–191.

Birchak, J., Gardner, C., Hipp, J. & Victor, J., 1974. High dielectric constant microwave probes for sensing soil moisture, *Proc. IEEE*, **62**, 93–98.

Brossier, R., Operto, S. & Virieux, J., 2009. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion, *Geophysics*, **74**, 105–118.

Chew, W., Jin, J. & Michielssen, E., 1997. Complex coordinate stretching as a generalized absorbing boundary condition, *Microwave Opt. Technol. Lett.*, **15**, 363–369.

Didier, Q., Arhab, S. & Lefeuve-Mesgouez, G., 2022. Comparative study between different frequency strategies for relative dielectric permittivity and electrical conductivity reconstruction. Application to near subsurface imaging, *EEE MTT-S International Conference on Numerical Electromagnetic and Multiphysics Modeling and Optimization (NEMO)*, pp. 1–4.

Didier, Q., Arhab, S. & Lefeuve-Mesgouez, G., 2024. Introducing a priori information with variable changes for a two-parameter reconstruction from experimental Fresnel Institute electromagnetic data, *IEEE Antennas Wireless Propag. Lett.*, **23**, 1774–1778.

Feng, X., Ren, Q. & Liu, C., 2017. Quantitative imaging for civil engineering by joint full waveform inversion of surface-based GPR and shallow seismic reflection data, *Construct. Build. Mater.*, **154**, 1173–1182.

Feng, D., Cao, C. & Wang, X., 2019. Multiscale full-waveform dual-parameter inversion based on total variation regularization to on-ground GPR data, *IEEE Trans. Geosci. Remote Sens.*, **57**, 9450–9465.

Gallardo, L. & Meju, M., 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of dc resistivity and seismic data, *Geophys. Res. Lett.*, **30**, 1658.

Gao, G., Abubakar, A. & Habashy, T., 2012. Joint petrophysical inversion of electromagnetic and full-waveform seismic data, *Geophysics*, **77**, 3–18.

Garambois, S., Sénéchal, P. & Perroud, H., 2002. On the use of combined geophysical methods to assess water content and water conductivity of near-surface formations, *J. Hydrol.*, **259**, 32–48.

Gassmann, F., 1951. Über die elaztizität poroser medien, *Vierteljahrsschrift der Naturfor-schenden Gesellschaft in Zurich*, pp. 1–23.

van Genuchten, M., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, *Soil Sci. Soc. Am. J.*, **44**, 892–898.

Glover, P., 2010. A generalized Archie's law for n phases, *Geophysics*, **75**, 247–265.

Habashy, T. & Abubakar, A., 2004. A general framework for constraint minimization for the inversion of electromagnetic measurements, *Progr. Electromagn. Res.*, **46**, 265–312.

Haber, E. & Holtzman, M., 2013. Model fusion and joint inversion, *Surv. Geophys.*, **34**, 675–695.

Lan, T., Liu, H., Liu, N., Li, J., Han, F. & Liu, Q., 2018. Joint inversion of electromagnetic and seismic data based on structural constraints using variational Born iteration method, *IEEE Trans. Geosci. Remote Sens.*, **56**, 436–445.

Lavoué, F., 2014. 2D Full Waveform Inversion of ground penetrating radar data: towards multiparameter imaging from surface data, *Earth Science, Université de Grenobles, NNT: 2014GRENU050*.

Lavoué, F., Frossier, R., Métivier, L. & Garambois, S., 2015. Frequency-domain modelling and inversion of electromagnetic data for 2D permittivity and conductivity imaging: An application to the Institut Fresnel experimental dataset, *Near Surf. Geophys.*, **13**, 227–241.

Loiseau, B., Carrière, S., Jougnot, D., Singha, K., Mary, B., Delpierre, N., Guérin, R. & Martin-St-Paul, N., 2023. The geophysical toolbox applied to forest ecosystems? A review, *Sci. Total Environ.*, **899**, 165503.

Mora, P., 1987. Nonlinear two-dimensional elastic inversion of multiooffset seismic data, *Geophysics*, **52**, 1211–1228.

Nocedal, J. & Wright, S., 2006. *Numerical Optimization*, Springer science & Business media.

Norton, S., 1999. Iterative inverse scattering algorithms: Methods of computing Fréchet derivatives, *J. acoust. Soc. Am.*, **106**, 2653–2660.

Nur, A., Mavko, G., Dvorkin, J. & Galmudi, D., 1998. Critical porosity: a key to relating physical properties to porosity in rocks, *Leading Edge*, **17**, 357–362.

Pride, S., 2005. Relationships between seismic and hydrological properties, *Hydrogeophysics*, **253**, 290–290.

Qin, Y., Rodet, T., Lambert, M. & Lesselier, D., 2021. Joint inversion of electromagnetic and acoustic data with edge-preserving regularization for breast imaging, *IEEE Trans. Comput. Imaging*, **7**, 349–360.

Qin, T., Bohlen, T. & Pan, Y., 2022. Indirect joint petrophysical inversion of synthetic shallow-seismic and multi-offset ground-penetrating radar data, *J. geophys. Int.*, **229**, 1770–1784.
Qin, T., Bohlen, T. & Pan, Y., 2024. Indirect joint petrophysical inversion of shallow-seismic and multi-offset ground-penetrating radar field data, *J. geophys. Int.*, **237**, 974–988.

Roger, A., 1983. Generalized reciprocity for gratings of finite conductivity, *Optica acta: Int. J. Optics*, **30**, 575–585.

Sirgue, L. & Pratt, R., 2004. Efficient waveform inversion and imaging: a strategy for selecting temporal frequencies, *Geophysics*, **69**, 231–248.

Song, X., Li, M., Yang, F., Xu, S. & Abubakar, A., 2019. Study on joint inversion algorithm of acoustic and electromagnetic data in biomedical imaging, *IEEE J. Multiscale Multiphys. Comput. Techn.*, **4**, 2–11.

Sonneveld, P., 1989. A fast Laczos-type solver for non symmetric linear systems, *J. Sci. Stat. Comput.*, **10**, 36–52.

Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation, *Geophysics*, **49**, 1259–1266.

Tarantola, A., 1986. A strategy for nonlinear elastic inversion of seismic reflection data, *Geophysics*, **51**, 1893–1903.

Taunk, K., De, S., Verma, S. & Swetapadma, A., 2019. A brief review of nearest neighbor algorithm for learning and classification, in *International Conference on Intelligent Computing and Control Systems (ICCS)*, IEEE, Madurai, India, pp. 1255–1260.

Tikhonov, A. & Arsenin, V., 1977. *Solutions of Ill-Posed Problems*, Winston.

Tromp, J., Tape, C. & Liu, Q., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, *J. geophys. Int.*, **160**, 195–216.

Wagner, F., Mollaret, C., Günther, T., Kemna, A. & Hauck, C., 2019. Quantitative imaging of water, ice and air in permafrost systems through petrophysical joint inversion of seismic refraction and electrical resistivity data, *J. geophys. Int.*, **219**, 1866–1875.

Wigin, A., 2004. *The inverse crime*, preprint, hal-id: hal-00001084.

Saltelli, A.,..., &, 2010. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index, *Comput. Phys. Commun.*, **181**, 259–270.

## APPENDIX A: FRÉCHET DERIVATIVE

In Section 3 we derived the Fréchet derivatives using some functions to integrate. This appendix details these functions. First we start with $\mathcal{J}_{\varepsilon_r}^n$ and $\mathcal{J}_{\sigma}^n$ which are the kernels of the Fréchet derivatives with

respect to $\varepsilon_r$ and $\sigma$ at iteration $n$.

$$ \begin{aligned} \mathcal{J}_{\varepsilon_r}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_i) &= \mathcal{J}_{\sigma}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_i) \\ &= E_n^{\text{sim}}(\mathbf{r}, \mathbf{r}_s, \omega_k) E_n^{\text{sim}}(\mathbf{r}, \mathbf{r}_d, \omega_k) \end{aligned} \quad (A1) $$

We observe that the Fréchet derivatives use the same kernel for these two parameters. The Fréchet derivative kernels for the seismic data are derived in two steps. The first step derives the kernels for the Lamé parameters and density. Then, a variable transformation is applied to calculate the kernels for the two wave velocities and a new one for the density. We begin with the first step:

$$ \mathcal{J}_{\rho}^{n}(\mathbf{r}, \mathbf{r}_{s}, \mathbf{r}_{d}, \omega_{k}) = \omega^{2} \mathbf{u}_{n}^{\text{sim}}(\mathbf{r}, \mathbf{r}_{s}, \omega_{k}) \cdot \mathbf{u}_{n}^{\text{sim}}(\mathbf{r}, \mathbf{r}_{d}, \omega_{k}) \quad (A2) $$

$$ \mathcal{J}_{\lambda}^{n}(\mathbf{r}, \mathbf{r}_{s}, \mathbf{r}_{d}, \omega_{k}) = -\nabla \cdot \mathbf{u}_{n}^{\text{sim}}(\mathbf{r}, \mathbf{r}_{s}, \omega_{k}) \nabla \cdot \mathbf{u}_{n}^{\text{sim}}(\mathbf{r}, \mathbf{r}_{d}, \omega_{k}) \quad (A3) $$

$$ \mathcal{J}_{\mu}^{n}(\mathbf{r}, \mathbf{r}_{s}, \mathbf{r}_{d}, \omega_{k}) = -\nabla \mathbf{u}_{n}^{\text{sim}}(\mathbf{r}, \mathbf{r}_{s}, \omega_{k}) : \nabla^{S} \mathbf{u}_{n}^{\text{sim}}(\mathbf{r}, \mathbf{r}_{d}, \omega_{k}) \quad (A4) $$

Using the relations $v_P = \sqrt{\frac{\lambda+2\mu}{\rho}}$ and $v_S = \sqrt{\frac{\mu}{\rho}}$, we can relate the above kernels to those of the velocities by differentiating these expressions. This also yields a new kernel for density:

$$ \mathcal{J}_{v_P}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k) = 2v_P(\mathbf{r})\rho(\mathbf{r})\mathcal{J}_{\lambda}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k) \quad (A5) $$

$$ \mathcal{J}_{v_S}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k) = 2v_S(\mathbf{r})\rho(\mathbf{r}) \times [\mathcal{J}_{\mu}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k) - 2\mathcal{J}_{\lambda}^n(\mathbf{r}, \mathbf{r}_s, \mathbf{r}_d, \omega_k)] \quad (A6) $$

$$ \begin{aligned} \mathcal{J}_{\rho}^{n}(\mathbf{r}, \mathbf{r}_{s}, \mathbf{r}_{d}, \omega_{k}) &= \mathcal{J}_{\rho}^{n}(\mathbf{r}, \mathbf{r}_{s}, \mathbf{r}_{d}, \omega_{k}) \\ &+ (v_{P}(\mathbf{r})^{2} - 2v_{S}(\mathbf{r})^{2})\mathcal{J}_{\lambda}^{n}(\mathbf{r}, \mathbf{r}_{s}, \mathbf{r}_{d}, \omega_{k}) \\ &+ v_{S}(\mathbf{r})^{2}\mathcal{J}_{\mu}^{n}(\mathbf{r}, \mathbf{r}_{s}, \mathbf{r}_{d}, \omega_{k}) \end{aligned} \quad (A7) $$

These kernels are computed using the generalized reciprocity theorem method. This implies the computation of two forward models for pair of source and receiver. This step could be time consuming but thanks to the measurement configuration, the number of times the forward model needs to be called is reduced to the number of sources.
<figure><img src="image_9.png" /><figcaption>**Figure A1.** Sobol indices (V: First-order Sobol index; TV: Total Sobol index) for electrical conductivity (a) and relative dielectric permittivity (b) with respect to input parameters.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure A2.** Sobol indices (V: First-order Sobol index; TV: Total Sobol index) for density (a), compressional wave speed (b) and shear wave speed (c) with respect to input parameters.</figcaption></figure>

## APPENDIX B: NORMAL EQUATION FOR SIJPI STRATEGY

In this article, we use a joint simultaneous petrophysical inversion that utilizes a combined normal equation derived from the two physical models. In this appendix, we will describe the gradients and Hessians involved. Let's begin with the expression for the Hessians:

$$H_{\text{petro}}^{E} = \begin{pmatrix} (\mathbf{J}_{E,\phi}^{n})^{\dagger} \mathbf{J}_{E,\phi}^{n} & (\mathbf{J}_{E,\phi}^{n})^{\dagger} \mathbf{J}_{E,S_{w}}^{n} \\ (\mathbf{J}_{E,S_{w}}^{n})^{\dagger} \mathbf{J}_{E,\phi}^{n} & (\mathbf{J}_{E,S_{w}}^{n})^{\dagger} \mathbf{J}_{E,S_{w}}^{n} \end{pmatrix} \quad (B1)$$

$$H_{\text{petro}}^{M} = \begin{pmatrix} (\mathbf{J}_{M,\phi}^{n})^{\dagger} \mathbf{J}_{M,\phi}^{n} & (\mathbf{J}_{M,\phi}^{n})^{\dagger} \mathbf{J}_{M,S_{w}}^{n} \\ (\mathbf{J}_{M,S_{w}}^{n})^{\dagger} \mathbf{J}_{M,\phi}^{n} & (\mathbf{J}_{M,S_{w}}^{n})^{\dagger} \mathbf{J}_{M,S_{w}}^{n} \end{pmatrix} \quad (B2)$$

Then, the gradients are given by the following forms:

$$G_{\text{petro}}^{E} = \begin{pmatrix} (\mathbf{J}_{E,\phi}^{n})^{\dagger} \delta E^{n}(\mathbf{r}_{d}, \mathbf{r}_{s}, \omega_{k}) \\ (\mathbf{J}_{E,S_{w}}^{n})^{\dagger} \delta E^{n}(\mathbf{r}_{d}, \mathbf{r}_{s}, \omega_{k}) \end{pmatrix} \quad (B3)$$

$$G_{\text{petro}}^{M} = \begin{pmatrix} (\mathbf{J}_{M,\phi}^{n})^{\dagger} \delta \mathbf{u}^{n}(\mathbf{r}_{d}, \mathbf{r}_{s}, \omega_{k}) \\ (\mathbf{J}_{M,S_{w}}^{n})^{\dagger} \delta \mathbf{u}^{n}(\mathbf{r}_{d}, \mathbf{r}_{s}, \omega_{k}) \end{pmatrix} \quad (B4)$$

In fact, the entire linear system is too large to be computed in its entirety, so we only evaluate the matrix-vector product during the inversion process of the linear system.

## APPENDIX C: SENSITIVITY ANALYSIS OF PETROPHYSICAL RELATIONSHIP

Figs A1 and A2 present the results of a sensitivity analysis on the petrophysical relationships used in this study. This analysis has been performed using Sobol indices which measure the influence of each input on the outputs of a function. In this case, the outputs are the geophysical parameters, while the inputs include the petrophysical parameters as well as the parameters defining the petrophysical relationships. For performing such an analysis we used the algorithm of Saltelli (2010) using a variance based algorithm to estimate Sobol indices. Moreover we used realistic ranges for each petrophysical parameters. Thus, if *p* and *k* would have a strong impact in case of large ranges, it is not the case here since their ranges are very narrow and have low values. This explains why some analyses deviate from theoretical expectation.