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

Research Paper

# PyLGRIM: Modelling 3D-ERI with infinite elements in complex topography context

Antoine Tonnoir<sup>a</sup>, Cyrille Fauchard<sup>b,*</sup>, Yannick Fargier<sup>c</sup>, Vincent Guilbert<sup>b</sup>, Raphael Antoine<sup>b</sup>

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

Check for
updates

## ARTICLE INFO

Dataset link: https://github.com/atonnoir/PyLGRIM.git

Keywords:

Electrical resistivity imaging

Topography

Finite element

Infinite element

## ABSTRACT

Electrical Resistivity Imaging (ERI) is one of the most used techniques in geophysics. As for many imaging methods, Digital Elevation Models (DEMs) are required to consider complex topography conditions. In this paper, we present some developments implemented into a new 3D-ERI software optimized in this context. The article focuses on the forward problem and discusses (i) the meshing methodology that directly consider DEMs in the processing and several profiles where electrodes are not necessarily aligned and (ii) new aspects for taking into account the unbounded domain. Indeed, defining boundary conditions of a numerical modelling problem arises as one of the most important issues into solving Partial Differential Equations (PDE). In order to solve the 3D-ERI forward problem, we propose an original implementation of the infinite elements, together with conventional finite elements. This methodology is first validated on synthetic case reproducing cliffs and, then, on a real case study presenting Badlands-like cliffs. Our results show that both the meshing procedure as well as the use of infinite elements enhance the efficiency of the forward problem as well as the accuracy of the inverse problem. In particular, this allows to reproduce more closely the local geology in complex environments than with a conventional 2D approach.

## 1. Introduction

Electric Resistivity Imaging has become increasingly important these last fifteen years (Loke et al., 2019) to provide useful subsurface information in civil-engineering (Dezert et al., 2019), archaeology (Fauchard et al., 2018) or Earth Sciences (Portal et al., 2019), due to its complementarity with other methods. ERI requires efficient modelling tools able to invert thousands of parameters in order to assess the distribution of electrical resistivity. To avoid artifacts and therefore misinterpretations, 3D informations such as surface topography have to be properly taken into account within the inversion process (Penz, 2012; Marescot, 2004; Rücker et al., 2017; Fargier et al., 2021). Nevertheless, such procedure may be particularly fastidious (1) due to the difficulties of meshing for complex structures and (2) due to the important cost of time and memory during the numerical inversion itself (Günther et al., 2006). Therefore, this challenge motivates the development of open source softwares for interdisciplinary and flexible modelling and inversion approaches.

Modelling ERI requires an accurate description of ground topography, especially around electrodes, where current injection and potential measurements are carried out (Rücker et al., 2006; Marescot et al.,

2006). In particular, taking into account several electrodes profiles together, where electrodes are not necessarily located on a straight line, is not classical and requires to consider the full 3D geometry of the domain of interest, as explained in the work of Udpuay et al. (2011). Also, modelling methods need consistent boundary conditions (BC) to formulate a well-posed problem (Hadamard, 1902). While the surrounding medium of a survey can be considered infinite, artificial boundaries for the numerical model, by definition, have to be set (Wood, 1976). To limit artifacts due to boundary conditions, a classical technique in ERI consists in putting homogeneous Dirichlet conditions far enough from the source injections (since the electric field decreases to infinity (Marescot, 2004)). This technique implies a large distance extension of the domain (and very stretched tetrahedron) increasing the computational cost of the forward problem. Other techniques also exist in the literature. Let us mention the so-called transparent boundary conditions (Bayliss et al., 1982) or absorbing boundary conditions (ABC), see Dey and Morrison (1979) for classical mixed boundary conditions or Meeker (2014) and Sugahara (2017) for more recent ABC. These methods are usually based on an analytical resolution of the problem in the exterior domain using polar (or ellipsoidal) coordinates. Yet, the analytical solution cannot be obtained

<footnote><sup>a</sup> INSA Rouen Normandie, LMI, France</footnote>

<footnote><sup>b</sup> Cerema, Research Team ENDSUM, France</footnote>

<footnote><sup>c</sup> Université Gustave Eiffel/GERS, France</footnote>

<footnote>* Corresponding author.
E-mail address: cyrille.fauchard@cerema.fr (C. Fauchard).</footnote>
when considering complex geometry. Finally, let us mention the infinite element method, see Wood (1976), Bettess (1977, 1980), Gerdes (1998, 2000) and Blome et al. (2009) for a more recent implementation of the method. This method consists in enriching the classical finite element basis with functions with infinite support (that decay at infinity). We expect these additional functions to catch the behaviour of the solution at infinity. Finally, this method is the most versatile and does not require strong assumptions on the DEM.

In this paper, we present a new Python based code called PyLGRIM, which stands for Python Language-based for Geoelectrical Resistivity Imaging and the methodology implemented for solving ERI problems. This work originates in the many measurements we have made in complex topography environment, such as vertical cliffs (Unique et al., 2024; Letortu et al., 2022), Badlands formations (Fauchard et al., 2023) or underground quarries (Fargier et al., 2017). It is also motivated by the desire to master and to freely share the entire process of modelling the direct problem as well as solving the inverse problem. The main features of the code are the following:

* We propose a quite simple procedure based on the use of GMSH (Geuzaine and Remacle, 2009) to integrate any digital elevation modelling (DEM) into the model calculated from photogrammetric or Lidar methods.

* PyLGRIM can handle several electrodes profiles to generate a mesh that exactly take into account the 3D position of the profiles. Also, we can invert together the data of several electrodes profiles at once.

* Finally, we propose a new infinite elements approach to bound the computational domain. In particular, the proposed method is cheap/efficient and robust to complex topography.

In a first section, we recall (i) the basic mathematical modelling of ERI methods, (ii) the procedure to take into account steep topography and several electrodes profiles, (iii) the discretization of the problem with infinite elements approach as well as (iv) the inverse problem. The second section focuses on (i) numerical examples and results (direct model and inversion) on flat case and a synthetic cliff, followed by (ii) a real case study located in France, on Badlands-like cliffs. In conclusion, we review the main results of this work and indicate the challenges that still need to be met in order to improve PyLGRIM. Finally, we specify the tools available for reproducing the work presented and using the software freely.

## 2. Mathematical modelling of ERI experiments/geoelectrical surveys

### 2.1. Mathematical formulation of the direct problem

The classical mathematical model for the simulation of ERI experiments/geoelectrical surveys is here described and shows how the equations are obtained (Zhdanov and Keller, 1994). ERI is based on the study of the circulation of an electric current inside the subsoil. The fundamental equation is then the Ohm's law:

$$ \vec{J} = \sigma \vec{E}, \qquad (1) $$

where $\vec{J}$ represents the density of the current (with unit in A m⁻²), $\sigma$ is the electrical conductivity (S/m) and $\vec{E}$ is the electric field (V m⁻¹). It is common to introduce also $\rho = \sigma^{-1}$ which corresponds to the resistivity (Ω m).

The electric field satisfies the static Maxwell-Faraday equation so that $\vec{\nabla} \times \vec{E} = \vec{0}$, where $\vec{\nabla}$ denotes the gradient and × is the cross product. Since $\vec{\nabla} \times \vec{E} = \vec{0}$, $\vec{E}$ derives from a scalar potential V (by Helmholtz decomposition theorem) and $\vec{E} = \vec{\nabla}V$. Moreover, without source term in any volume D, the conservation of charge is written as follow:

$$ - \int_{\partial D} \vec{J} \cdot \vec{n} \, d\Gamma = 0, \qquad (2) $$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Representation of the subsoil V and the ground topography S. $V_b$ shows the volume of interest that contains the profile of electrodes at positions $\varepsilon$.</figcaption></figure>

where $\vec{n}$ is the outward unit normal, meaning that the intensity flux over the boundary of D is null, and $\Gamma$ is the boundary of the volume D. By flux-divergence theorem and since the above result is true for any volume D, we deduce that:

$$ - \vec{\nabla} \cdot \vec{J} = 0. \qquad (3) $$

Finally, with (1) and $\vec{E} = \vec{\nabla}V$, the equation satisfied by the potential V is given by:

$$ - \vec{\nabla} \cdot (\sigma \vec{\nabla} V) = 0. \qquad (4) $$

This equation is set in an open domain V corresponding to the subsoil (see Fig. 1). To formulate a well-posed problem, it must be equipped with boundary conditions. First, at infinity, we impose the potential to vanish. Second, on the boundary corresponding to the floor, we impose the intensity flux to be zero, that is to say $\sigma \vec{\nabla} V \cdot \vec{n} = 0$, which amounts to consider that the air is a perfect insulator. It remains to model the source term which corresponds to the current injection with an electrode. Classically, this is represented by a Dirac point source located on the surface. To sum up, the problem to solve is the following: Find the potential V satisfying the following PDE system:

$$ \left\{ \begin{array}{ll} -\vec{\nabla} \cdot (\sigma \vec{\nabla} V) = c_{\varepsilon} \delta_{x_{\varepsilon}, y_{\varepsilon}, z_{\varepsilon}} & \text{in } V \\ \sigma \vec{\nabla} V \cdot \vec{n} = 0 & \text{on } S, \\ \lim_{\varepsilon \to \infty} V = 0 & \end{array} \right. \qquad (5) $$

where $\delta_{x_{\varepsilon}, y_{\varepsilon}, z_{\varepsilon}}$ is the Dirac distribution located at the electrode position $(x_{\varepsilon}, y_{\varepsilon}, z_{\varepsilon})$ and $c_{\varepsilon}$ is a constant that represents the current density injected on the electrode $\varepsilon$. It is known that this problem has a unique solution V (Hadamard, 1902).

Since the computational domain must be bounded, we denote by $V_b$ the domain of interest (see Fig. 1). More precisely, we have:

$$ V_b = \{ (x, y) \in [-L_x - x_{min}, x_{max} + L_x] \times [-L_y - y_{min}, y_{max} + L_y] \text{ and } z \in [z_{min} - L_z, S(x, y)] \} \qquad (6) $$

where $S(x, y)$ is the function that describes the topography of the ground, $L_x > 0$, $L_y > 0$, $L_z > 0$ are three parameters that give the position of the artificial boundary of $V_b$ (in particular $L_z$ describes the depth to model). The parameters $x_{min} < x_{max}$, $y_{min} < y_{max}$ and $z_{min} < z_{max}$ represent the bounds of the domain of interest and will be clarify in the next section. In the following, we denote by $V_{\infty} = V \setminus V_b$ the (infinite) exterior domain.

### 2.2. Mesh topography

Complex environment DEMs can be obtained by Lidar (Brodu and Lague, 2012; Hobbs et al., 2019; Milan et al., 2007; Roulland et al., 2021) technique or by photogrammetry (Stocker et al., 2015; Medjkane et al., 2018; Fauchard et al., 2023). Typically, this leads to get a *.xyz file listing the xyz-coordinates of a point cloud representing the 3D
<figure><img src="image_4.png" /><figcaption>Fig. 2. Example of mesh generated by script generateMsh.py of PyLGRIM / GMSH in the case of 3 electrodes profiles (P2, P4, P5, see Table 1). In Figures (a) and (b), the scale in x, y and z are not the same. The line separating the two colours of the mesh in Figure (c) corresponds to the recursively computed path linking all the electrodes.</figcaption></figure>

geometry of the ground. The next step consists in building from this DEM a volume mesh representing the subsoil in depth. This operation is not classical in standard photogrammetry softwares and we will explain here the procedure implemented in PyLGRIM. This procedure, based on mesh deformation as in Loke (2000), Rücker et al. (2006), Udhpuay et al. (2011) and Blanchy et al. (2020) is detailed here. Given an *.xyz file representing the topography, i.e. a list of 3D coordinates (xᵢ, yᵢ, zᵢ = S(xᵢ, yᵢ)ᵢ=1,...,N) (N being the number of points), we proceed in two steps:

1. Denoting by (x_min, x_max), (y_min, y_max) and (z_min, z_max) the minima and maxima of the coordinates (xᵢ, yᵢ, zᵢ), we build a rectangular cube of dimension [x_min, x_max] × [y_min, y_max] × [z_min - L_z, z_min], where the parameter L_z > 0 describes the depth of the basement one wishes to consider. Using GMSH (Geuzaine and Remacle, 2009), it is straightforward to create and mesh this cube. The integration of the electrode profiles positions is yet less obvious and it is discussed hereafter.

2. This mesh is then deformed by changing the z_p coordinate of each point into $\bar{z}_p$ thanks to the following transformation:

$$ \bar{z}_p = z_p + \frac{z_p - (z_{min} - L_z)}{L_z} \times (z_M - z_{min}), \quad (7) $$

where $z_M$ is the z coordinate of the nearest point in the list of points of the DEM. With this transformation, all points of the cube with z-coordinate $z_p = z_{min}$ get as new coordinate $z_M$, so that the top of the cube is deformed to become the DEM.

To take into account the exact electrodes positions ($x'_j, y'_j, z'_j$)ᵢ=1,...,Nᵢ (where Nᵢ is the number of electrodes on electrode profile l), the coordinates ($x'_j, y'_j, z'_j$)ᵢ=1,...,Nᵢ are added to the list of the DEM coordinates. Also, we must add the positions ($x'_j, y'_j, z_{min}$)ᵢ=1,...,Nᵢ to the mesh of the cube. In GMSH, this requires to construct a path linking the electrodes. In particular, if two electrodes profiles intersect, we must split each profiles of electrodes into two profiles so that no new electrodes profiles intersect. The computation of this path is automatically done in PyLGRIM and it is implemented through a recursive algorithm, which works efficiently up to 5 possibly intersecting electrodes profiles (for more profiles of 64 electrodes the computation times can become very long). This also enables to build a mesh which is refined around the electrodes and under the electrodes profiles and therefore in better adequacy with the sensitivity of a 2D profile. This enables also a better comparison between 3D inversion results and independent 2D inversion results since the slice of the 3D mesh under the electrode profile corresponds directly to a 2D mesh (if the electrodes line is straight). For more than 5 electrodes profiles, a simpler (but faster) algorithm is implemented which leads to a mesh which is not refined under electrodes profiles.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Representation of the support of an infinite element test function lying on the boundary {x = x_max + L_x}. In the interior domain $V_b$, the tetrahedron is represented.</figcaption></figure>

The Fig. 2 presents an example of application of this procedure with PyLGRIM on a case of study: the “Vaches Noires” cliffs (VNCs) (Villers-sur-Mer, Normandy, France). On this example, the required data to build the mesh are the DEM points cloud describing the ground topography, and three *.xyz files corresponding to electrodes positions for each profiles (see the documentation in https://github.com/atonnoir/PyLGRIM.git). PyLGRIM then generates three files: meshDEMFlat.geo describing the geometrical positions of the electrodes on the cube, meshDEMFlat.msh corresponding to the 3D mesh of the cube (without the ground topography), and meshDEM.msh corresponding to the 3D deformed mesh (taking into account the 3D ground topography).

## 2.3. Discretization aspects

Let us now discuss the discretization aspects. Formally, the variational formulation of problem (5) is given by: For any test function $Ψ$ we seek for V solution to:

$$ a(V, \Psi) = c_{\varepsilon} \Psi(x_{\varepsilon}, y_{\varepsilon}, z_{\varepsilon}), \quad \text{where} \quad a(V, \Psi) = \int_V \sigma \nabla V \cdot \nabla \Psi \, dV. \quad (8) $$

To numerically solve this problem, a finite element space $V$ is considered to discretize the unknown $V$. In the bounded domain $V_b$, based on classical mesh, usual Lagrange finite elements discretization are used.

In the unbounded domain $V_\infty$, the test functions are defined as the tensor product of classical Lagrange finite elements on triangle, representing the faces of the tetrahedron on artificial boundaries of $V_b$, and a decay function with infinite support (see Fig. 3).

More precisely, we consider $Ψ$ a Lagrange test function in $V_b$ such that $Ψ$ is not null on the boundary {x = x_max + L_x} for instance, and let
<figure><img src="image_6.png" /><figcaption>Fig. 4. Representation of the support of an infinite element test function on the edge {$x = x_{max} + L_x$} × {$y = y_{min} - L_y$} (on the left) and on the corner {$x = x_{max} + L_x$} × {$y = y_{min} - L_y$} × {$z = z_{min} - L_z$} (on the right).</figcaption></figure>

P(y, z) = Ψ(x_{max} + L_x, y, z). Then, we define the associated function as follows:

$$N(x, y, z) = P(y, z)D(|x - x_{max} - L_x|), \quad \text{with} \quad \left\{ \begin{array}{l} D(0) = 1 \\ D(x) \underset{x \to +\infty}{\longrightarrow} 0 \cdot \end{array} \right. \qquad (9)$$

Note that we impose $D(0) = 1$ to ensure the continuity of the test function. Similarly, we can defined the infinite elements on the boundaries {$x = x_{min} - L_x$}, {$y = y_{min} - L_y$}, {$y = y_{max} + L_y$} and {$z = z_{min} - L_z$} of $V_b$. We must also handle carefully the element on the edges and on the corners of the artificial boundary of $V_b$. Indeed, to define continuous functions, we need to construct element on the edge and on the corner as a tensor product. For instance, on the line {$x = x_{max} + L_x$} × {$y = y_{min} - L_y$}, we have:

$$N^e(x, y, z) = P(z) \left( D(|x - x_{max} - L_x|)D(|y - y_{min} + L_y|) + D^e(|x - x_{max} - L_x|, |y - y_{min} + L_y|) \right), \qquad (10)$$

where P is the (trace of the) Lagrange function on a vertical segment, see Fig. 4, and $D^e$ is a chosen function s.t. $D^e(0, \cdot) = D^e(\cdot, 0) = 0$. This last condition ensures the continuity of the test function $N^e$ (for instance, $N^e(x, y_{min} - L_y, z) = P(z)D(|x - x_{max} - L_x|)$ which corresponds to the test function (9) on the same boundary $x = x_{max} + L_x$). In the same way, we have, for instance, on the corner {$x = x_{max} + L_x$} × {$y = y_{min} - L_y$} × {$z = z_{min} - L_z$}:

$$N^c(x, y, z) = D(|x - x_{max} - L_x|)D(|y - y_{min} + L_y|)D(|z - z_{min} + L_z|) + D^c(|x - x_{max} - L_x|, |y - y_{min} + L_y|, |z - z_{min} + L_z|), \qquad (11)$$

where $D^c$ is a chosen function ensuring once again the continuity of the test function, i.e. $D^c(0, \cdot, \cdot) = D^c(\cdot, 0, \cdot) = D^c(\cdot, \cdot, 0) = 0$. As we will see below, the precise choice of the functions $D$, $D^e$ and $D^c$ is not essential (the main point will be their integrals) as soon as they are integrable.

Let us now detail the computation of the integral terms in the variational formulation for infinite element test functions. For simplicity, we focus on a particular triangle $\mathcal{T}$ on boundary {$x = x_{max} + L_x$}. We denote by $P_i(y, z)$, $i \in \{1, ..., m\}$, the set of Lagrange polynomial functions in this triangle, such that:

$$\left\{ \begin{array}{l} P_i(y, z) = a_i + b_i y + c_i z + \cdots + d_i y^p + e_i z^p; \quad i = 1...m, \\ P_i(y_j, z_j) = \delta_{ij} \end{array} \right. ,$$

with $m$ the degree of freedom (DOF) in the triangle. We also denote by $N_i = P_i D(|x - x_{max} - L_x|)$ the associated infinite test functions defined as in (9). To compute the integral term of the variational formulation

associated to this function, we need to compute:

$$\begin{align} & \int_{x_{max}+L_x}^{+\infty} \int_{\mathcal{T}} \vec{\nabla} N_i(x, y, z) \cdot (\sigma \vec{\nabla} N_j(x, y, z)) \, dxdydz \nonumber \\ &= \int_{0}^{+\infty} D^2(x)dx \int_{\mathcal{T}} \vec{\nabla} P_i(y, z) \cdot \vec{\nabla} P_j(y, z)dydz \nonumber \\ & \quad + \int_{0}^{+\infty} (D'(x))^2 dx \int_{\mathcal{T}} P_i(y, z)P_j(y, z)dydz \nonumber \\ &= \alpha \mathbb{K}_{i,j}^{2D} + \alpha'_{e} \mathbb{M}_{i,j}^{2D} \tag{12} \end{align}$$

where $\alpha > 0$ and $\alpha' > 0$ are the integrals of $D^2$ and $(D')^2$ on the semi-line $[0, +\infty]$ respectively, and $\mathbb{K}^{2D}$ and $\mathbb{M}^{2D}$ are the classical rigidity and mass Lagrange finite element matrices in 2D. In the same way, we get for the functions (10):

$$\begin{align} & \int_{x_{max}+L_x}^{+\infty} \int_{-\infty}^{y_{min}-L_y} \int_{\mathcal{L}} \vec{\nabla} N_i^e(x, y, z) \cdot (\sigma \vec{\nabla} N_j^e(x, y, z)) \, dxdydz \nonumber \\ &= \alpha_e \mathbb{K}_{i,j}^{1D} + \alpha'_{e} \mathbb{M}_{i,j}^{1D} \tag{13} \end{align}$$

where $\alpha_e$ and $\alpha'_{e}$ are the integrals of $(D(\cdot)D(\cdot)+D^e(\cdot, \cdot))^2$ and $(D'(\cdot)D'(\cdot)+D^e(\cdot)'(\cdot))^2$, and for the functions (11):

$$\begin{align} & \int_{x_{max}+L_x}^{+\infty} \int_{-\infty}^{y_{min}-L_y} \int_{-\infty}^{z_{min}-L_z} \vec{\nabla} N^e(x, y, z) \cdot (\sigma \vec{\nabla} N^e(x, y, z)) \, dxdydz \nonumber \\ &= \alpha'_{e} \tag{14} \end{align}$$

In other words, to compute the integral terms with these infinite elements, we simply need to compute classical finite element matrices, and to add it to the global rigidity matrix multiplied by the coefficients $\alpha$, $\alpha_{e,c}$ or $\alpha'_{e,c}$. These coefficients are the parameters we shall choose. In fact, we can also notice that it corresponds to consider ABC with compatibility relations at the edges and the corners, in the same spirit as in Modave et al. (2020) for wave equation.

## 2.4. Recall on the inverse problem

The discretization of the direct problem (8) with infinite elements leads to solve a linear system of the form:

$$K(\vec{\sigma})\vec{V} = \vec{F}, \qquad (13)$$

where $\vec{V} \in \mathbb{R}^n$ represents the potential solution on each nodes of the mesh, $\vec{F}$ the source term and $\vec{\sigma} \in \mathbb{R}^p$ the conductivity in each cell of the mesh, supposed to be constant ($p$ being the number of cells).

Let us define now the inverse problem we are interested in. The purpose is to reconstruct the conductivity $\vec{\sigma}$ in each cell given apparent resistivity measurements on the surface, that we will denote by a vector $\vec{d}_{exp} \in \mathbb{R}^m$, where $m$ is the number of data. Let us also introduce the observation operator $\odot$ represented by an $m \times n$ matrix that maps the
solution $\vec{V}$ of the (discretized) direct problem (13) to the observations. Then, the data $\vec{d}_{exp}$ can be expressed as follows:

$$ 
\begin{aligned}
\vec{d}_{exp} &= \mathcal{O}(\vec{V} + \vec{\eta}), \\
&= \mathcal{O}([K(\vec{\sigma}_{ex})]^{-1} \vec{F} + \vec{\eta}),
\end{aligned} 
\quad (14) $$

where $\vec{\eta}$ represents the noise and $\vec{\sigma}_{ex}$ is the exact parameter we are looking for. Classically, the idea is to formulate the inverse problem as a minimization problem: Find the conductivity $\vec{\sigma}$, that minimizes

$$ J(\sigma) = \frac{1}{2} \|\vec{d}_{exp} - \mathcal{O}([K(\vec{\sigma})]^{-1} \vec{F})\|_{2}^{2} + \frac{\epsilon}{2} R(\vec{\sigma}), \quad (15) $$

where the first term above corresponds to the discrepancy to the data, and the second term $R(\vec{\sigma})$ to the regularization term (which accounts for our a priori knowledge on the medium). The parameter $\epsilon > 0$ corresponds to the weight given to the regularization part. In our case, we consider as regularization term a smoothing matrix $L$ so that

$$ R(\vec{\sigma}) = \|L\vec{\sigma}\|_{2}^{2}, $$

where

$$ L_{i,j} = \begin{cases} -1 & \text{if } j \in N^i, \\ \#N^i & \text{if } i = j, \\ 0 & \text{elsewhere} \end{cases} $$

and $N^i$ is the set of neighbouring cells to cell $i$. A second smoothing matrix is also implemented in the code to enforced stratified media:

$$ L_{i,j}^{\text{strat}} = \begin{cases} -\theta_{i,j} & \text{if } j \in N^i, \\ \sum_{j \in N^i} \theta_{i,j} & \text{if } i = j, \\ 0 & \text{elsewhere} \end{cases} $$

where $\theta_{i,j}$ is a parameter that gradually varies between 0, if the barycenters of cells $i$ and $j$ have the same positions $(x, y)$, and 1, if the barycenters of cells $i$ and $j$ have the same altitudes $z$. Thought this is only empirical, this smoothing matrix penalizes variations in direction $x$ and $y$, and not in direction $z$.

To solve the minimization problem (15), we have implemented a Gauss–Newton algorithm where the regularization parameter $\epsilon$ varies during the iterations. The main difficulty then is to compute efficiently the sensitivity matrix (or Jacobian matrix) of $\mathcal{O}([K(\vec{\sigma})]^{-1} \vec{F}$. This is done through the adjoint state method which amounts to solve direct problems. For the seek of completeness, let us recall the method. We start from the fact that $K(\vec{\sigma} + \vec{\delta}\sigma) = K(\vec{\sigma}) + K(\vec{\delta}\sigma)$, so that using first order Taylor expansion we get

$$ 
\begin{aligned}
\vec{V}(\vec{\sigma} + \vec{\delta}\sigma) &= [K(\vec{\sigma} + \vec{\delta}\sigma)]^{-1} \vec{F}, \\
&= [\mathbb{I} + [K(\vec{\sigma})]^{-1} K(\vec{\delta}\sigma)]^{-1} \underbrace{[K(\vec{\sigma})]^{-1} \vec{F}}_{=\vec{V}}, \\
&= \vec{V}(\vec{\sigma}) - [K(\vec{\sigma})]^{-1} K(\vec{\delta}\sigma) \vec{V}(\vec{\sigma}) + \vec{O}(\|\vec{\delta}\sigma\|^2).
\end{aligned} 
$$

In other words, the variation of the potential $\vec{\delta}\vec{V}$ considering a variation of the conductivity $\vec{\delta}\sigma$ is given by

$$ \vec{\delta}\vec{V} = - [K(\vec{\sigma})]^{-1} K(\vec{\delta}\sigma) \vec{V}(\vec{\sigma}) + \vec{O}(\|\vec{\delta}\sigma\|^2). $$

Taking the scalar product with $\mathcal{O}^t \vec{e}_i$, where $\vec{e}_i$ is the $i$th canonical vector of $\mathbb{R}^m$, and using the fact that $[K(\vec{\sigma})]^{-1}$ is symmetric, we get

$$ 
\begin{aligned}
(\mathcal{O}\vec{\delta}\vec{V}, \vec{e}_i) &= -(K(\vec{\delta}\sigma)\vec{V}, [K(\vec{\sigma})]^{-1}\mathcal{O}^t\vec{e}_i) + \vec{O}(\|\vec{\delta}\sigma\|_2^2) \\
&= -(K(\vec{\delta}\sigma)\vec{V}, \vec{P}_i) + \vec{O}(\|\vec{\delta}\sigma\|_2^2)
\end{aligned} 
$$

where $\vec{P}_i$ is the adjoint variable that we need to compute by solving the direct problem $[K(\vec{\sigma})]^{-1} \mathcal{O}^t \vec{e}_i$. To conclude, we simply need to recall that the conductivity $\sigma$ is constant on each cell, so that the sensitivity with respect to cell $k$ of the $i$th observation is given by:

$$ \frac{\partial(\mathcal{O}\vec{V}, \vec{e}_i)}{\partial\sigma_k} = -(K(\vec{1}_k)\vec{V}, \vec{P}_i) \quad (16) $$

<figure><img src="image_7.png" /><figcaption>Fig. 5. Mesh used for the flat case study (containing 106 453 nodes and 82 182 tetrahedra). The electrodes are located on the straight line on the top surface in the middle (at the interface between the two coloured meshes in green and orange). The dimension of the cube are $V_b = [-0.5, 0.5] \times [-0.5, 0.5] \times [-1, 0]$.</figcaption></figure>

where $\vec{1}_k$ is the vector of $\mathbb{R}^p$ equal to 0 except at position $k$ where it is equal to 1.

In the implementation in the code, to reduce the computational cost, a storage in a compact form (i.e. storing only non zero values) of all the local matrix $K(\vec{1}_k)$ is used. Doing so, at each step of the Gauss–Newton algorithm, only the adjoint states and the scalar product (16) must be computed. Finally, the Gauss–Newton iterative algorithm of minimization is stopped when we reach an admissible error, i.e. when the functional $J$ is lesser than a given tolerance.

### 3. Numerical examples & results

In this last section, numerical results obtained with PyLGRIM are presented as well as some examples of use of the codes. The whole codes are written in Python and are decomposed into three main scripts (see the documentation and the codes in https://github.com/atonnoir/PyLGRIM.git):

1. DEM: This script aims at constructing the 3D mesh of the domain of study. The methodology has been described in Section 2.2.

2. Simu ERI: This script aims at solving the direct problem using infinite elements described previously. Also, other boundary conditions can be used: simple homogeneous Dirichlet conditions or mixed boundary conditions.

3. Inv ERI: This script aims at solving the inverse problem of resistivity reconstruction based on electrode measurements. At the moment, the data must be formatted as apparent resistivity.

The following numerical examples can be reproduced thanks to the tutorial given with the software. In the next examples, the first purpose will be to demonstrate the interest of the proposed BC in the case of a cliff topography. To do so, we will first focus on the direct problem. Then, we will show some inversion results taking into account several profiles of electrodes.

#### 3.1. Synthetic flat case (Direct problem)

We consider first the situation of a half space $V = \mathbb{R}^2 \times \{z \le 0\}$ perfectly homogeneous and isotropic with resistivity $\rho = 1$ $\Omega$ m
<figure><img src="image_8.png" /><figcaption>Fig. 6. Comparison of the matrices E defined by (18) corresponding to the error percentage between the approximate solution and the analytical one using IE (left) and MBC (right). The colours represent the value of the entries of the matrix E$_{i,j}$ at row i and column j.</figcaption></figure>

(we consider for simplicity adimensional quantity for the synthetic examples) and a profile of 64 electrodes equally spaced on the segment {(-0.4, 0, 0), (0.4, 0, 0)} on the floor, with an inter electrode spacing d = 0.0125 (see Fig. 1). In that case, considering a current injection on one electrode i, we know that the analytical solution V$_{ex}$ is given by:

$$V_{ex}^{l}(\vec{x}) = \frac{1}{2\pi \|\vec{x} - \vec{x}_{j}\|}$$

where $\vec{x}_j$ is the location of electrode i. Then, we can construct the matrix $P_a$ which entry $[P_a]_{ij}$ corresponds to the measurement of electrical potential $V_{ex}^l$ at electrode i ≠ j due to a current injection at electrode j ≠ i (by convention, we will set $[P_a]_{ii} = 0$), or in other words:

$$[P_a]_{ij} = \frac{1}{2\pi d_{ij}} \qquad (17)$$

where $d_{ij} = \|\vec{x}_i - \vec{x}_j\|$ is the distance between the two electrodes i and j. This analytical solution will allow us to evaluate in this simple configuration the effect of the boundary conditions.

For the numerical computations, we have considered the cube $V_b = [-0.4 - L_x, 0.4 + L_x] \times [-L_y, L_y] \times [-L_z, 0]$. The mesh, see Fig. 5, has been generated with the first script (DEM) and the generated mesh is refined around electrodes so that we have approximately 5 tetrahedrons between two electrodes. The numerical results were obtained using second order tetrahedral finite elements. Let us note that in the code first order elements are also available. Also, the conversion from linear to quadratic finite element is implemented in the first script (see the Tutorial with the code for more details). The electrical current injected on each electrode $\mathcal{E}$ is $c_{\mathcal{E}} = 1$ A.

Taking $L_x = 0.1$, $L_y = 0.5$ and $L_z = 1.0$, we have compared the so-called mixed boundary conditions (MBC) with the Infinite Element (IE). To do so, we computed the error matrix:

$$E = \frac{|P_a - P|}{P_a} \times 100 \qquad (18)$$

where P is the matrix of electrical potential [V] measurement numerically obtained. On Fig. 6, we have represented the error matrix in the two cases. The problem being commutative the matrix P and $P_a$ are symmetric, and so is the matrix E. Also, as we can expect, the error is larger for electrodes located close to the artificial boundaries (electrode 1 and 64). For the two methods, the results are quite similar. More precisely, the maximum error $\max_i E_{i,j}$ is smaller with IE (10.8% versus 15.5%) whereas it is greater in $L^2$-norm $\sum_{i \neq j} E_{i,j}^2 / 64^2$ (5.7% versus 3.5%). As a first conclusion, we can say that, for a similar computational cost (the two methods involve similar computations since they both correspond to a Robin type boundary condition), we get similar performances.

## 3.2. Synthetic cliff case (Direct problem)

Let us now consider a simple cliff topography and the domain $V = \mathbb{R}^2 \times \{z \le f(x, y)\}$ where

$$f(x, y) = \begin{cases} 0 & \text{if } x \le 0.3 \\ \frac{x-0.3}{0.3} & \text{if } x \in [0.3, 0.6] \\ 0.5 & \text{if } x \ge 0.6 \end{cases} \qquad (19)$$

is a function that described the ground topography, see Fig. 7. This time, we supposed the medium to be heterogeneous. More precisely, we consider that:

$$\rho = \begin{cases} 100 & \text{if } z \le -0.25 \ \& \ (x, y, z) \notin B_1 \cup B_2 \\ 40 & \text{if } z > -0.25 \ \& \ (x, y, z) \notin B_1 \cup B_2 \end{cases}$$

where $B_1$ and $B_2$ are two sphere of centre (0.65, -0.2, -0.0) and (0.2, 0.35, -0.25) and radii 0.15 and 0.1 respectively. In the two spheres, we have set $\rho = 5$ Ω m. The choice of the sphere is a classic example in geophysics to study the response of methods to an anomaly for example, see the response of a sphere in the electrical method presented by Telford et al. (1990), p. 530. Let us mention that the inhomogeneities are taken into account in PyLGRIM as a simple text file given in each cell the value of the resistivity, supposed to be constant inside the element.

In this configuration, we have computed the matrix P of potential measurement (described previously in flat case configuration) using a profile of 64 electrodes equally spaces between (0.1, 0, 0) and (0.9, 0, 0.5). The computational domain $V_b$ is given by $[0.1 - L_x, 0.9 + L_x] \times [-L_y, L_y] \times [-L_z, 0]$, $f(x, y)$. Since no analytical solution is available in this case, to measure the quality of the results, we have computed the matrix P twice: first with a “small” domain, taking $L_x = 0.1$, $L_y = 0.5$ and $L_z = 1$, and then with a “large” domain, taking $L_x = 0.6$, $L_y = 1$ and $L_z = 1.5$. The corresponding mesh are presented on Fig. 8.

Denoting by $P_s$ and $P_l$ the obtained matrices in each case (’s’ for small and ’l’ form large), we have computed the relative error matrix:

$$E_{diff} = \left| \frac{P_s - P_l}{P_l} \right| \qquad (20)$$

where we assume that the matrix $P_l$ is closer to the true matrix of potential measurement. One expect of course that if the boundary condition on the artificial boundary exactly takes into account the infinite exterior domain, then we would have $E_{diff} = 0$ (in other words, the solution does not depend on the size of the computational domain if the boundary condition is “transparent”). This matrix $E_{diff}$ has been computed using IE and using MBC. As we can see on Fig. 9, the error is smaller when using IE (the maximum error is 12.74% with IE whereas it is 19.44% with MBC). Let us mention that we have also compared the computed matrix $P_l$ (in the large domain) using IE or MBC, and the maximum relative difference is 3.5%, which means that in both case, the computation in the large domain gives similar results (and the impact of the artificial boundary is small at this distance).
<figure><img src="image_9.png" /><figcaption>Fig. 7. Representation of the synthetic cliff topography and the subsoil $V$ containing the two spheres and the two layers. On the right, we have represented with colours the value of the resistivity of the heterogeneities.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Meshes used for the cliff case study (on the left, for the small domain, with 106453 nodes and 84182 elements, and on the right, for the large domain, with 185135 nodes and 142692 elements). The dimension of the small domain is $[0, 1] \times [-0.5, 0.5] \times [-1, f(x, y)]$ and the dimension of the large domain is $[-0.5, 1.5] \times [-1, 1] \times [-1.5, f(x, y)]$, where the topography is given by the function $f$ defined in (19).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Comparison of the matrices $E_{diff}$ defined by (20) corresponding to the error percentage between the solution computed in the small domain and the solution computed in the large domain using IE (on the left) and MBC (on the right). The colours represent the value of the entries of the matrix $E_{i,j}$ at row $i$ and column $j$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 10. Representation of the mesh for data inversion in the synthetic cliff topography (259 937 Nodes and 202 249 elements). The dimension of the domain is [−0.35, 1.25]× [−0.6, 0.7]× [−1, f(x, y)], where the topography is given by the function f defined in (19).</figcaption></figure>

### 3.3. Synthetic cliff case (Inverse problem)

In this section, we still consider the topography and configuration presented in the previous section, see Fig. 7, but this time for the inverse problem. Using PyLGRIM, we have generated synthetic data of resistivity on two profiles as illustrated on Fig. 10. The two profiles are composed of 32 electrodes equally spaced between points (−0.15, −0.2, 0) and (0.85, −0.2, 0.5) for the first profile, and between (0.1, 0.4, 0) and (0.75, −0.3, 0.5) for the second profile. The data have been generated on a very large domain (using second order finite elements) so that IE BC and MBC give very close apparent resistivity of data resistivity. To get apparent resistivity data, we solve the direct problem as mentioned in the very large domain and knowing the heterogeneities. Thus, we get the potential V. We compute also the solution V<sup>hom</sup> when considering a homogeneous medium with resistivity 1. Then, the dimensionless apparent resistivity is given for a quadripole ABMN by

$$ \rho_a = \frac{V_{AM} - V_{BM} - V_{AN} + V_{BN}}{V_{AM}^{hom} - V_{BM}^{hom} - V_{AN}^{hom} + V_{BN}^{hom}} $$

where V<sub>ij</sub> corresponds to the electrical potential at electrode i injecting at electrode j.

The inversion process yet is carried out on the smaller domain (to avoid “inverse crime” and to get a reasonable computational time) represented on Fig. 10 using also second order finite elements. The dataset used for the inversion is composed of two Wenner acquisition array about 1110 quadripoles, one for each profile, without cross profiles quadripoles. Each inversion starts with a prior homogeneous model whose resistivity equal the mean of measured apparent resistivities. We have compared the inversion process using either IE BC and MBC. On Fig. 11, we have represented the convergence curves in each case. As we can see, both BC give similar RMS error.

On Fig. 12, we have represented the reconstructed solution at different iterations of the inversion process. As we can see, the choice of the BC has a non negligible impact on the reconstructed resistivity distribution. In particular, we can observe on Fig. 13 (representing slices of the solution under each profile of electrodes) that the two anomalies are better identified when using IE for the BC, although the interface between the two layers is more straight when using the MBC. Let us recall that the inverse problem of conductivity reconstruction is very ill-posed and usually has not a unique solution (even with a regularization term). This is why the two reconstructed solutions can be different although both have a similar RMS error.

<figure><img src="image_13.png" /><figcaption>Fig. 11. Convergence curves of the RMS error using MBC or IE versus the iteration of Gauss–Newton algorithm.</figcaption></figure>

### 3.4. Case study: the “Vaches Noires Cliffs”, a Badlands formation

#### 3.4.1. Context and measurements

In order to illustrate the PyLGRIM capabilities in 3D environment, an electrical survey was carried out on coastal context, a Badlands formation presenting steep crests and deep valleys made of successive layers of clays and marls, with a thick layer of weathered limestone. The Vaches Noires Cliffs (VNCs, in Normandy, France), face the sea and are 4.5 km long from Villers-sur-Mer to Houlegate. This site has already been the subject of many publications (Roulland et al., 2021; Maquaire et al., 2013; Hassen et al., 2021). The locale geology is shown in Fig. 14. As can be seen, the local geology is tabular. The expected resistivity images should show alternating highly conductive materials (marls, clays) and more resistant materials (limestones). Locally we should also find areas of debris due to the strong erosion caused by the plateau hydrology. The details of the geophysical results and the diachronic study of eroded and accumulated volumes of materials over 4 years, carried out by drone and photogrammetry and with help of ERI methods have been recently presented in Fauchard et al. (2023). One of the prospects of this previous work was to improve the geophysical inversion results by taking into account the very rugged topography of the studied area. Indeed, the inversions made in two dimensions with commercial software, even if they included the 2D topography of the profiles with the required parameters (adapted distorted grid and flatness options), hardly restore the local multilayered geology. The goal of this example is to show the ability of PyLGRIM to produce more realistic results than other classical approaches.

We focus on a section of the VNCs, 200 m long, 75 m height and 300 m large. Six ERI profiles were performed using an ABEM Terrameter LS2 in Wenner and Dipole-dipole configurations for P1, P2, P5 and P6 profiles and only in Wenner configuration for P3 and P4 profiles: their technical characteristics are presented in the Table 1. A survey was performed by drone flying at about 100 m height, taking dozens of photos in order to generate a point cloud, precisely georeferenced with twenty Differential GPS targets on the ground (Guilbert et al., 2020). The obtained DEM has a point density around 555 points/m². The ERI profiles location and the drone flight planning are represented with an aerial view of the studied site in the Fig. 15, as well as the DEM in Fig. 16, (a).

#### 3.4.2. Results and interpretation

The inverted 3D data block obtained with PyLGRIM is represented with the Paraview software (Ayacht, 2015) in Fig. 16, (b). The block is composed of 441 426 tetrahedrons and the inversion is performed (i) taking into account of all data (Wenner and Dipole-dipole) and (ii) carried out with a computer composed of a Intel core i9 9900KF at 3.6 GHz and DDR4 RAM 64 Go. The inversion process stopped after
<figure><img src="image_14.png" /><figcaption>Fig. 12. Resistivity reconstruction at iterations {1,3,5,7} using IE BC (left) or MBC (right). The grey sphere and plane represent the position of the two anomalies and the interface between the two layers in the domain we wish to recover.</figcaption></figure>

<caption>**Table 1**</caption>

Characteristics of the ERI profiles. Profiles P1, P2 and P5, P6 have been performed in both Wenner and Dipole-dipole configuration while P3 and P4 in Wenner configuration only.

<table><thead><tr><th>Name</th><th>Init. nb.<br>of quad.<br>W/DD</th><th>Final nb.<br>of quad.<br>W/DD</th><th>Stacking<br>Err. at 1%</th><th>Nb. elec.</th><th>Length (m)<br>(m)</th><th>Spacing (m)<br>(m)</th></tr></thead><tbody><tr><td>P1</td><td>555/888</td><td>555/888</td><td>2</td><td>64</td><td>126</td><td>2</td></tr><tr><td>P2</td><td>555/888</td><td>555/888</td><td>2</td><td>64</td><td>126</td><td>2</td></tr><tr><td>P3</td><td>555/Ø</td><td>555/Ø</td><td>2</td><td>64</td><td>126</td><td>2</td></tr><tr><td>P4</td><td>555/Ø</td><td>555/Ø</td><td>2</td><td>64</td><td>126</td><td>2</td></tr><tr><td>P5</td><td>1110/1776</td><td>1110/1758</td><td>2</td><td>80</td><td>316</td><td>4</td></tr><tr><td>P6</td><td>555/888</td><td>480/798</td><td>2</td><td>59</td><td>232</td><td>4</td></tr></tbody></table>

more than 24 h at the fourth iteration (L2-norm) where the root mean square error $\epsilon_{RMS}$ is equal to 0.16 for IE and 0.2 for MBC. This is of course something we need to improve in the future, using intensive computing methods. From this 3D-block, six vertical slices are extracted according to the location of the corresponding six ERI profiles. They are

represented into the DEM in Fig. 16, (c). The Fig. 17 shows a special view of the ERI profiles obtained with PyLGRIM, juxtaposing the DEM with the inverted profiles. The advantage of the 3D inversion is the perfect continuity between the profiles parallel to the coastline and the transverse profiles running from the plateau to the beach, despite
<figure><img src="image_15.png" /><figcaption>Fig. 13. Slices of the resistivity reconstruction at iterations {7} using IE BC (left) or MBC (right) under the two profiles of electrodes.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Locale geology of the VNCs (Normandy, France), reproduced from Fauchard et al. (2023). Expected ERI results should show alternating layers of conductive (marls, clay) and resistive (limestone) materials with local deposition of materials (debris) due to the strong erosion caused by the plateau hydrology.</figcaption></figure>

the difference in the electrode spacing (2 m for P1 to P4, 4 m for P5 and P6). The high resistivity at the bottom of the inverted profile represents a non significative value (over 300 Ω m) where the sensitivity is high and is located far from the depth of investigation of ERI profiles, classically ranging to 1/6 to 1/8 of the profile length (Roy and Apparao,

1971). DOI limit and sensitivity matrix representation are necessary and form parts of the future priorities to be achieved.

To better compare the approaches presented in Section 2 using either MBC or IE, the P5 and P6 profiles are shown in Figs. 18 and 19, respectively. Both inversions are carried out with a smoothing matrix
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Location of the six ERI profiles (P1-6) transverse and longitudinal to the shoreline (red dots = position of electrodes) with the flight planning (white line) and the takeoff point (blue point) of the drone (Normandy, France).</figcaption></figure>

*Source:* Adapted from Fauchard et al. (2023).

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** (a) DEM of the VNCs performed with a drone survey, adapted from Fauchard et al. (2023). (b) DEM in low opacity and 3D-block of inverted data with PyLGRIM obtained with the six ERI profiles and superimposed. (c) DEM in low opacity and the six inverted ERI profiles, obtained by slicing vertically the 3D-block resistivity at the electrodes position with the Paraview software (Ayachit, 2015).</figcaption></figure>

$L^{strat}$ depicted in Section 2.4. The PyLGRIM results are compared with a famous commercial software which is a reference in the domain. The objective is not to evaluate both approaches since this well-known

solution offers very complete interpretation and inversion scenarios in ERI measurements. Our aims is to underline how a complete 3D approach that takes full account of the topography as well as the entire
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Six ERI profiles extracted by slicing vertically the 3d-block shown in Fig. 16. The global overview shows a perfect concordance between transverse and parallel profiles to the shoreline.</figcaption></figure>

ERI measurements of all the profiles (Dipole-dipole and Wenner) can provide a more accurate representation of the subsurface. P5 (Fig. 18) and P6 (Fig. 19) profiles were carried out from the plateau and passed in the valley between crests. Since no DOI or sensitivity matrix is proposed until now, the contour of the Wenner configuration results of the commercial software is extracted and superimposed on PyLGRIM results to discussed comparisons. The P5 profile is 316 m long and stops at approximately 40 m on the beach from the basal scarp while the P6 profile Fig. 18 is 232 m long and stops at the limit between the basal scarp and the beach. In the P5 profile, the limestone formation (see Fig. 14 for the locale geology) of the top cliffs are described by two resistive layers separated by a less resistive formation (the brown clay). This description is not underline with the commercial software, where this upper part of the cliff is only described by the upper grey marl formation followed by a more resistive layer corresponding to limestone. Conversely, PyLGRIM offers a less precise description of surface anomalies after the first escarpment, where are present debris composed of both marls and limestone materials. In the slope after the first escarpment, only PyLGRIM describes a layered structure corresponding to the marl of Villers where the commercial software seems to diverge in function of depth. In the P6 profile (Fig. 19), the same conclusion can be addressed. The upper formation of limestone and grey marls are ill-described with commercial software while very near surface anomalies in the first scarp seems more precisely delineated. In depth, both Dipole-dipole and Wenner measurements with the commercial software describe perpendicular formations while PyLGRIM tends to a more realistic description of the geology.

## 4. Conclusion

In this paper, we present a new software called PyLGRIM that allows to process 3D ERI modelling and inversion. PyLGRIM is based on scripts in Python language that contains all elements shown in the present paper and integrating the GMSH software for mesh issue. The first outcome is the capability of PyLGRIM to input any DEM performed on the field to take into account complex topographies encountered in geophysical surveys. This approach based on a nearest neighbour algorithm between a 2D projection of the points cloud of the DEM and the points of the mesh cube is certainly well known, but it is clearly explained here. The second main outcome is an algorithm that has been

developed to make the mesh of profiles which may intersect with each other consistent with the mesh of the DEM. The third main results deals with a new formulation of infinite elements which offers an original absorbing boundary condition, whose performances are demonstrated to be more efficient than the common boundary conditions in particular for complex topography. In order to present some of the PyLGRIM performances, an analytical validation illustrates the aforementioned meshing and boundary discussion. Finally, a real case study in Badlands context shows how taking into account the complex topography in geophysical processing and both all the profiles and configurations carried out on the field provides a more relevant description of materials distribution in the subsurface than simply process measurements with classical 2D approaches.

In terms of perspectives, many improvements remain to be made. The use of complementary inversion methods is currently being implemented. Similarly, the need to calculate the DOI index will enable better interpretation of the inverted data as well as the sensitivity matrix. Work is in progress on the inversion of 3D data collected on complex structures such as quarry pillars, vertical cliffs or parallelepiped laboratory samples. Ensuring the utilization of high-performance calculation methods stands as a paramount objective. Notably, PyLGRIM, a freely available software, extends an open invitation to users to actively participate in enhancing and optimizing the software.

### CRediT authorship contribution statement

**Antoine Tonnoir:** Writing – original draft, Software, Formal analysis, Conceptualization. **Cyrille Fauchard:** Writing – original draft, Visualization, Supervision, Project administration, Methodology, Funding acquisition. **Yannick Fargier:** Writing – review & editing, Methodology, Conceptualization. **Vincent Guilbert:** Visualization, Validation. **Raphael Antoine:** Writing – review & editing.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Cyrille Fauchard reports administrative support, article publishing charges, equipment, drugs, or supplies, and travel were provided by Normandy Region. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
## Data availability

The PyLGRIM code is written in Python and is available in the following GitHub repository: https://github.com/atonnoir/PyLGRIM.git It contains all the directory and the documentation “Tuto basic.pdf” and “PyLGRIM Link To Github Software and Data.pdf” required to run PyLGRIM as well as to reproduce the present work. Furthermore, all geophysical data (DEM with ERI measurements in VTK files) are available and can be displayed with the needed Paraview (version 5.9.1 and later) project.

## Acknowledgements

This work has been conducted in the context of the French Regional Project TéléDETaC, (“Télédéction par drone du trait de côte”) dealing with the Unmanned Aerial Vehicle (or drone) teledetection of the coastal line and its erosion. Its development is continuing within the framework of the DEFHY3GEO project (Detection and Study of Fracturing by HYdrology, GEOmorphodynamics, Geology and Geophysics, 2022–2025). The authors would like to thank the Normandy Region, France for funding these two research projects. The

authors also thank the council of Villers-sur-Mer who gave access to the protected area of the VNCs, in order to perform ERI measurements and drone flights. All measurements and ERI data processing were carried out with helps of the both research teams ENDSUM (Cerema), the LMI (UR 3226, Insa Rouen Normandy), and IDEES (UMR CNRS 6266, Caen University).

## Code and data availability section

The PyLGRIM code is written in Python and is available in the following GitHub repository:

https://github.com/atonnoir/PyLGRIM.git

It contains all the directory and the documentation “Tuto basic.pdf” and “PyLGRIM Link To Github Software and Data.pdf” required to run PyLGRIM as well as to reproduce the present work. Furthermore, all geophysical data (DEM with ERI measurements in VTK files) are available and can be displayed with the needed Paraview (version 5.9.1 and later) project.

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** P5 profile (316 m long, 64 electrodes, spacing 4 m, see the Table 1 for further details) with (a) the Wenner (7 iteration (it.), ε<sub>RMS</sub> = 2.5%) and (b) Dipole–dipole (it = 5, ε<sub>RMS</sub> = 7.3%) configuration inverted with commercial software compared with (c) IE and (d) MBC obtained with PyLGRIM (both it = 4, ε<sub>RMS</sub> < 0.2%. The described materials 1 to 5 correspond to the local geology shown in Fig. 14.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. 19.** P6 profile (236 m long, 64 electrodes, spacing 4 m, see the Table 1 for further details) with (a) the Wenner (6 iteration (it.), ε<sub>RMS</sub> = 3.5%) and (b) Dipole–dipole (it = 7, ε<sub>RMS</sub> = 5%) configuration inverted with commercial software compared with (c) IE and (d) MBC obtained with PyLGRIM (both it = 4, ε<sub>RMS</sub> < 0.2%. The described materials 1 to 5 correspond to the local geology shown in Fig. 14.</figcaption></figure>

References

Ayachit, U., 2015. The ParaView Guide: A Parallel Visualization Application. Kitware, Inc., Clifton Park, NY, USA.

Bayliss, A., Gunzburger, M., Turkel, E., 1982. Boundary conditions for the numerical solution of elliptic equations in exterior regions. SIAM J. Appl. Math. 42 (2), 430–451, URL: http://www.jstor.org/stable/2101222.

Bettess, P., 1977. Infinite elements. Int. J. Numer. Methods Eng. 11 (1), 53–64. http://dx.doi.org/10.1002/nme.1620110107.

Bettess, P., 1980. More on infinite elements. Internat. J. Numer. Methods Engrg. 15 (11), 1613–1626. http://dx.doi.org/10.1002/nme.1620151105.

Blanchy, G., Saneiyan, S., Boyd, J., McLachlan, P., Binley, A., 2020. ResIPy, an intu-
itive open source software for complex geoelectrical inversion/modeling. Comput.
Geosci. 137, 104423. http://dx.doi.org/10.1016/j.cageo.2020.104423, URL: https:
//www.sciencedirect.com/science/article/pii/S0098300419308192.

Blome, M., Maurer, H., Schmidt, K., 2009. Advances in three-dimensional geoelectric forward solver techniques. Geophys. J. Int. 176 (3), 740–752. http://dx.doi.org/10.1111/j.1365-246X.2008.04006.x.

Brodu, N., Lague, D., 2012. 3D terrestrial lidar data classification of complex natural scenes using a multi-scale dimensionality criterion: Applications in geomorphology. ISPRS J. Photogramm. Remote Sens. 68, 121–134. http://dx.doi.org/10.1016/j.isprsjprs.2012.01.006, URL: https://www.sciencedirect.com/science/article/pii/S0924271612000330.

Dey, A., Morrison, H., 1979. Resistivity modelling for arbitrarily shaped two-
dimensional structures. Geophys. Prospect. 27 (1), 106–136. http://dx.doi.org/10.
1111/j.1365-2478.1979.tb00961.x.

Dezert, T., Lopes, S.P., Fargier, Y., Côte, P., 2019. Combination of geophysical and geotechnical data using belief functions: Assessment with numerical and laboratory data. J. Appl. Geophys. 170, 103824. http://dx.doi.org/10.1016/
j.jappgeo.2019.103824, URL: http://www.sciencedirect.com/science/article/pii/
S0926985118310413.

Fargier, Y., Antoine, R., Dore, L., Lopes, S.P., Fauchard, C., 2017. 3D assessment of an
underground mine pillar by combination of photogrammetric and geoelectric meth-
odscombination of photogrammetry and 3D-ERI. Geophysics 82 (4), E143–E153.
http://dx.doi.org/10.1190/geo2016-0274.1.

Fargier, Y., Dezert, T., Antoine, R., Tonnoir, A., Fauchard, C., 2021. Impact of the dem's resolution on 2d and 3d-eri common practice. In: NSG2021 2nd Conference on Geophysics for Infrastructure Planning, Monitoring and BIM. European Association of Geoscientists & Engineers, pp. 1–5.

Fauchard, C., Guilbert, V., Antoine, R., Ledun, C., Beaucamp, B., Maquaire, O., Costa, S., Medjkane, M., Roulland, T., 2023. Diachronic UAV study of coastal badlands supported by geophysical imaging in the context of accelerated erosion processes. Landslides http://dx.doi.org/10.1007/s10346-022-02006-2.

Fauchard, C., Saley, A.D., Camerlynck, C., Fargier, Y., Antoine, R., Thérain, P.-F., 2018. Discovery of the Romanesque church of the Abbey of our lady of Bec (Le Bec-Hellouin, Normandy, France) by means of geophysical methods. Archaeol. Prospect. 25 (4), 315–328. http://dx.doi.org/10.1002/arp.1711, arXiv:https://onlinelibrary.
wiley.com/doi/pdf/10.1002/arp.1711, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/arp.1711.

Gerdes, K., 1998. Solution of the 3D-Helmholtz equation in exterior domains of arbitrary shape using hp-finite-infinite elements. Finite Elem. Anal. Des. 29 (1), 1-20. http://dx.doi.org/10.1016/j.camwa.2012.06.010.

Gerdes, K., 2000. A review of infinite element methods for exterior Helmholtz problems. J. Comput. Acoust. - J COMPUT ACOUST 8, 43-62. http://dx.doi.org/10.1016/S0218-396X(00)00004-2.

Geuzaine, C., Remacle, J.-F., 2009. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. Internat. J. Numer. Methods Engrg. 79 (11), 1309-1331. http://dx.doi.org/10.1002/nme.2579, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.2579, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2579.

Guilbert, V., Antoine, R., Heinkelé, C., Maquaire, O., Costa, S., Gout, C., Davidson, R., Sorin, J.-L., Beaucamp, B., Fauchard, C., 2020. Fusion of thermal and visible point clouds: Application to the vaches noires landslide normandy, France. ISPRS - Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci. XLIII-B2-2020, 227-232. http://dx.doi.org/10.5194/isprs-archives-XLIII-B2-2020-227-2020, URL: https://www.int-arch-photogramm-remote-sens-spatial-inf-sci.net/XLIII-B2-2020/227/2020/.

Günther, T., Rücker, C., Spitzer, K., 2006. Three-dimensional modelling and inversion of dc resistivity data incorporating topography — II. Inversion. Geophys. J. Int. 166 (2), 506-517. http://dx.doi.org/10.1111/j.1365-246X.2006.03011.x, arXiv:https://academic.oup.com/gji/article-pdf/166/2/506/1498513/166-2-506.pdf.

Hadamard, J., 1902. Sur Les Problèmes Aux Dérivées Partielles Et Leur Signification Physique. Princeton University Bulletin, pp. 49-52.

Hassen, I., Fauchard, C., Antoine, R., Roulland, T., Maquaire, O., Costa, S., Dugué, O., 2021. 3D geological modelling of a coastal area: case study of the Vaches Noires Cliffs, Normandy, France. Bull. Eng. Geol. Environ. 80 (2), 1375-1388. http://dx.doi.org/10.1007/s10064-020-01955-z.

Hobbs, P.R.N., Jones, L.D., Kirkham, M.P., Pennington, C.V.L., Morgan, D.J.R., Dashwood, C., 2019. Coastal landslide monitoring at Aldbrough, East Riding of Yorkshire, UK. Q. J. Eng. Geol. Hydrogeol. 53 (1), 101-116. http://dx.doi.org/10.1144/jegjeh2018-210, arXiv:https://pubs.geoscienceworld.org/jegjeh/article-pdf/53/1/101/4941954/jegjeh2018-210.pdf.

Junique, T., Antoine, R., Guilbert, V., Beaucamp, B., Ledun, C., Costa, S., Maquaire, O., Davidson, R., Fauchard, C., 2024. Investigation of the geological and hydrogeological structure of chalk cliffs with visible, thermal infrared and electrical resistivity imaging. J. Hydrol. (ISSN: 0022-1694) 630, 130642. http://dx.doi.org/10.1016/j.jhydrol.2024.130642, https://www.sciencedirect.com/science/article/pii/S0022169424000362.

Letortu, P., Le Dantec, N., Augereau, E., Costa, S., Maquaire, O., Davidson, R., Fauchard, C., Antoine, R., Flahaut, R., Guirriecc, Y., Longuevergne, L., de la Bernardie, J., David, L., 2022. Experimental field study on the fatigue and failure mechanisms of coastal chalk cliffs: Implementation of a multi-parameter monitoring system (Sainte-Marguerite-sur-Mer, France). Geomorphology 408, 108211. http://dx.doi.org/10.1016/j geomorph.2022.108211, URL: https://www.sciencedirect.com/science/article/pii/S0169555X22001040.

Loke, M., 2000. Topographic modelling in resistivity imaging inversion. In: 62nd EAGE Conference & Technical Exhibition Extended Abstracts, D-2.

Loke, M.H., Rucker, D., Dahlin, T., Chambers, J.E., 2019. Recent advances in the geoelectrical method and new challenges: a software perspective. FastTimes 24 (4), 56-62.

Maquaire, O., Affchain, P., Costa, S., Lissak, C., Fressard, M., Letortu, P., Davidson, R., 2013. Evolution à long terme des falaises des "Vaches Noires" et occurrence des glissements (Calvados, Basse-Normandie, France). In: Jounées Aléas Gravitaires. URL: https://isterre.fr/IMG/pdf/maquaire.pdf.

Marescot, L., 2004. Modélisation Directe Et Inverse En Prospection Électrique Sur Les Structures 3D Complexes Par La Méthode Des Éléments Finis (Ph.D. thesis). Universités de Nantes et de Lausanne.

Marescot, L., Rigobert, S., Palma Lopes, S., Lagabrielle, R., Chapellier, D., 2006. A general approach for DC apparent resistivity evaluation on arbitrarily shaped 3D structures. J. Appl. Geophys. 60 (1), 55-67. http://dx.doi.org/10.1016/j.jappgeo.2005.12.003, URL: http://www.sciencedirect.com/science/article/pii/S0926985106000024.

Medjkane, M., Maquaire, O., Costa, S., Roulland, T., Letortu, P., Fauchard, C., Antoine, R., Davidson, R., 2018. High-resolution monitoring of complex coastal morphology changes: cross-efficiency of SfM and TLS-based survey (Vaches-Noires cliffs, Normandy, France). Landslides http://dx.doi.org/10.1007/s10346-017-0942-4.

Meeker, D.C., 2014. Improvised asymptotic boundary conditions for electrostatic finite elements. IEEE Trans. Magn. 50 (6), 1-9. http://dx.doi.org/10.1109/TMAG.2014.2300196.

Milan, D.J., Heritage, G.L., Hetherington, D., 2007. Application of a 3D laser scanner in the assessment of erosion and deposition volumes and channel change in a proglacial river. Earth Surf. Process. Landf. 32 (11), 1657-1674. http://dx.doi.org/10.1002/esp.1592.

Modave, A., Geuzaine, C., Antoine, X., 2020. Corner treatments for high-order local absorbing boundary conditions in high-frequency acoustic scattering. J. Comput. Phys. 401, 109029. http://dx.doi.org/10.1016/j.jcp.2019.109029, URL: https://www.sciencedirect.com/science/article/pii/S0021999119307351.

Penz, S., 2012. Modélisation Et Inversion De Données Électriques En Courant Continu: Vers Une Prise En Compte Efficace De La Topographie (Ph.D. thesis). Ecole Nationale Supérieure des Mines de Paris.

Portal, A., Fargier, Y., Labazuy, P., Lénat, J.-F., Boivin, P., Miallier, D., 2019. 3D electrical imaging of the inner structure of a complex lava dome, Puy de Dôme volcano (French Massif Central, France). J. Volcanol. Geotherm. Res. 373, 97-107. http://dx.doi.org/10.1016/j.jvolgeores.2019.01.019, URL: http://www.sciencedirect.com/science/article/pii/S0377027318303561.

Roulland, T., Maquaire, O., Costa, S., Medjkane, M., Davidson, R., Fauchard, C., Raphaël, A., 2021. Seasonal activity quantification of coast badlands by TLS monitoring over five years at the "Vaches Noires" cliffs (Normandy, France). Geomorphology 10803. http://dx.doi.org/10.1016/j geomorph.2021.108083, URL: https://www.sciencedirect.com/science/article/pii/S0169555X21004918.

Roy, A., Apparao, A., 1971. Depth of investigation in direct current methods. Geophysics 36 (5), 943-959. http://dx.doi.org/10.1190/1.1440226, arXiv:https://doi.org/10.1190/1.1440226.

Rücker, C., Günther, T., Spitzer, K., 2006. Three-dimensional modelling and inversion of dc resistivity data incorporating topography – I. Modelling. Geophys. J. Int. 166 (2), 495-505. http://dx.doi.org/10.1111/j.1365-246X.2006.03010.x.

Rücker, C., Günther, T., Wagner, F.M., 2017. pyGIMLi: An open-source library for modelling and inversion in geophysics. Comput. Geosci. 109, 106-123. http://dx.doi.org/10.1016/j.cageo.2017.07.011, URL: https://www.sciencedirect.com/science/article/pii/S0098300417300584.

Stocker, C., Eltner, A., Karrasch, P., 2015. Measuring gullies by synergetic application of UAV and close range photogrammetry — A case study from Andalusia, Spain. Catena 132, 1-11. http://dx.doi.org/10.1016/j.catena.2015.04.004, URL: http://www.sciencedirect.com/science/article/pii/S0341816215001228.

Sugahara, K., 2017. Improvised asymptotic boundary conditions for magnetostatic field problems in ellipsoidal and elliptic cylindrical domains. IEEE Trans. Magn. 53 (6), 1-4. http://dx.doi.org/10.1109/TMAG.2017.2665645.

Telford, W.M., Geldart, L.P., Sheriff, R.E., 1990. Applied Geophysics, second ed. Cambridge University Press, Cambridge, http://dx.doi.org/10.1017/CBO9781139167932, URL: https://www.cambridge.org/core/product/FA576F8F33FCA4D56975BC58CE07851E.

Udphuay, S., Günther, T., Everett, M.E., Warden, R.R., Briaud, J.-L., 2011. Three-dimensional resistivity tomography in extreme coastal terrain amidst dense cultural signals: application to cliff stability assessment at the historic D-Day site. Geophys. J. Int. 185 (1), 201-220. http://dx.doi.org/10.1111/j.1365-246X.2010.04915.x, arXiv: http://oup.prod.sis.lan/gji/article-pdf/185/1/201/1641817/185-1-201.pdf.

Wood, W.L., 1976. On the finite element solution of an exterior boundary value problem. Internat. J. Numer. Methods Engrg. 10 (4), 885-891. http://dx.doi.org/10.1002/nme.1620100415, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620100415, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620100415.

Zhdanov, M.S., Keller, G.V., 1994. The geoelectrical methods in geophysical exploration, Elsevier ed. Methods in Geochemistry and Geophysics, Elsevier.