# Seismoelectric wave conversions at an interface: a quantitative comparison between laboratory data and full-waveform modelling

Victor Martins Gomes, D Brito, Stéphane Garambois, Michel Dietrich, Clarisse Bordes, Hélène Barucq

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

## To cite this version:

Victor Martins Gomes, D Brito, Stéphane Garambois, Michel Dietrich, Clarisse Bordes, et al. Seismoelectric wave conversions at an interface: a quantitative comparison between laboratory data and full-waveform modelling. Geophysical Journal International, 2023, 235 (3), pp.2992-3011. 10.1093/gji/ggad409. hal-04396286

HAL Id: hal-04396286

https://inria.hal.science/hal-04396286v1

Submitted on 18 Jan 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Seismoelectric wave conversions at an interface: a quantitative comparison between laboratory data and full-waveform modelling

V. Martins-Gomes,<fnref n="1" />,<fnref n="2" /> D. Brito<fnref n="1" />,<fnref n="1" /> S. Garambois<fnref n="1" />,<fnref n="3" /> M. Dietrich,<fnref n="3" /> C. Bordes<fnref n="1" /> and H. Barucq<fnref n="2" />

Accepted 2023 October 12. Received 2023 October 7; in original form 2023 May 22

## SUMMARY

Seismo-electromagnetic phenomena, electrokinetic in nature, take place whenever a seismic wave propagates in fluid-bearing media, its energy depending mainly on the electrical properties of the fluid and the hydraulic properties of the porous medium. They result from a conversion of mechanical into electromagnetic (EM) energy due to the transient ionic interactions occurring at the pore scale. Two of these phenomena are usually studied: the electric field accompanying seismic waves, and the EM field that travels independently, generated at discontinuities of physicochemical properties in the porous medium. Although the first event is sensitive to physical parameters of the surrounding medium, the second catches information about interfaces in the subsurface, with the resolution of seismic methods, making it very attractive to near surface exploration. In this context, we propose a new experimental setup where both phenomena can be simultaneously studied. At first, we use a porous medium composed of homogeneous water-saturated sand and study the characteristics of the coseismic electric field. Afterwards, a thin layer of Vosges sandstone is inserted into the sand, which allows the study of the EM waves generated at the two closely spaced sand-sandstone interfaces. We record the seismic displacement field at the upper surface of the sand volume using a laser vibrometer, and use stainless steel electrodes buried in the sand to acquire individual electric potentials rather than electric fields, seeking to favour the measurement of the EM interface-generated signals. With the help of direct numerical simulations, we compare experimental measurements and theoretical predictions, based on a well established set of seismoelectric governing equations. In both types of experiments, this comparison shows very good agreements between experimental and numerical waveforms, thus confirming the relevant theory. The electric potential data also show that the EM signals generated at interfaces are clearly recorded at distances of about 10 seismic *P* wavelengths away from the interface. By contrast, the same events are barely noticeable near the inserted layer when measured using classical electric dipolar arrays.

**Key words:** Electromagnetic theory; Numerical modelling; Hydrogeophysics; Wave propagation.

## 1 INTRODUCTION

Research on seismoelectric (SE) and its reciprocal electroseismic (ES) phenomena, electrokinetic in nature, dates back to the 1930s (Thompson 1936). However, dedicated laboratory studies only began around 1960 (see Parkhomenko 1971). Attempts to establish the governing equations of such phenomena were first made by Frenkel (1944), based on Ivanov's field observations (Ivanov 1939, 1940) and later by Neev & Yeatts (1989). Nonetheless, it was Pride (1994), complemented by Pride & Haartsen (1996) and Haartsen & Pride

(1997), that developed a sound theory by coupling Biot's theory for the propagation of seismic waves in porous media and Maxwell's electromagnetic (EM) formulation via frequency-dependent flux-force relations. More recently, alternative electrokinetic models were developed by Revil & Mahardika (2013), Jougnot et al. (2013), Jardani & Revil (2015) and Rosas-Carbajal et al. (2020), seeking to describe seismoelectromagnetic phenomena for a broader class of porous media.

EM signals, created by electrokinetic phenomena, originate in fluid-bearing porous media whenever a gradient of pressure forces

<footnote n="1">¹Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, TotalEnergies, LFCR, 64000 Pau, France. E-mail: daniel.brito@univ-pau.fr</footnote>

<footnote n="2">²Makutu, Inria, TotalEnergies, Université de Pau et des Pays de l'Adour, LMAP, CNRS, 64000 Pau, France</footnote>

<footnote n="3">³Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERre, 38000 Grenoble, France</footnote>
the pore-fluid to flow inside the pores or fractures. Microscopically (i.e. for grain size dimensions), it is the time-varying interaction between charges adsorbed on the surface of the grains and the ions moving with the pore fluid that governs all electrokinetic effects. This class of EM signals can be induced by global scale displacement, at a constant rate, of the fluids inside the rocks (streaming potential) or by the propagation of seismic waves (seismoelectromagnetic phenomena; Jouniaux & Bordes 2012). In the latter case, which is the focus of this work, a frequency-dependent coefficient is used to describe the coupling between seismic and EM wavefields. Two formulations of the frequency-dependent coupling coefficient are often mentioned in the literature: that of Packard (1953), conceptualized for a medium made of capillary tubes and that of Pride (1994), discussed hereinafter (see Reppert et al. 2001 for a comparison). More recently, Jougnot & Solazzi (2021) describes an upscaling procedure, adapted from Jougnot et al. (2012), used to formulate the frequency-dependent coupling coefficient as a function of salinity when the porous media is approximated by a bundle of capillaries. The proposed model is then successfully compared to that of Packard (1953) and that of Revil & Mahardika (2013). Shi et al. (2018) and Thanh et al. (2021) have also studied the frequency-dependency of the coupling coefficient for capillaries, with both proposing original formulations: the first one focusing on the normalized radius of the pore (the ratio between the pore radius and the Debye length) and the latter on the pore size distribution.

Following Pride's electrokinetic theory, research has been mainly directed towards two effects of seismic-to-EM wave conversion. For a propagating compressional seismic wave, the first effect is a co-seismic electric field whose support is strictly bound to that of the mechanical wave (Pride & Haartsen 1996). The second effect stems from the discontinuities of the medium's electrical and/or mechanical properties across interfaces, and results in an EM disturbance (EM-disturbance) that propagates at a much higher speed than the generating seismic wave (Haartsen & Pride 1997). Both effects have been studied experimentally at the laboratory scale (see Zhu et al. 1994, 2008; Bordes et al. 2008; Schakel et al. 2011a). Additionally, Pride & Garambois (2002) demonstrated that Biot slow waves play a significant role on the generation of EM signals at geological interfaces. They have quantitatively established that accurate modelling of this signal requires the utilization of the complete set of Biot's equations, including the proper consideration of Biot slow waves, in addition to *P* and *S* waves.

Numerical simulations of the seismoelectromagnetic response using Pride's formulation were first published by Haartsen & Pride (1994) for stratified porous media, with further details given in Haartsen (1995) and Haartsen & Pride (1997). These authors adapted the global matrix method (Chin et al. 1984) to solve the electrokinetic governing equations for point sources in layered porous media. They were followed by Garambois (1999) and Garambois & Dietrich (2002) who introduced a slightly different methodology by extending the generalized reflection and transmission matrix method (Kennett & Kerry 1979) to compute the SE and ES point-source responses. Ren et al. (2010) and Grobbe et al. (2016) also developed numerical techniques to simulate the 3-D propagation of SE wavefields in layered media. Complementarily, purely numerical methods for solving partial differential equations were used by Haines & Pride (2006) (finite-difference method), Zyserman et al. (2010) (finite-element method), Morency et al. (2021) (spectral-element method) and Meyer (2021) (discontinuous Galerkin method), to tackle the problem of simulating the SE phenomena in complex 2-D or 3-D media.

Experimental studies of the coseismic wavefield, in the kilohertz range, have mainly been focused on the influence of physical properties on the coseismic electric field amplitude. Thus, Zhu et al. (2000) investigated the effect of the electrical conductivity of the saturating fluid, Bordes et al. (2015) concentrated on the saturation rate and Holzhauer et al. (2017) analysed the results of varying both conductivity and saturation. Laboratory experiments devoted to the observation of interfacial SE responses were first reported by Zhu & Toksöz (1996). In their work, EM waves were produced by a compressional seismic wave crossing an interface separating a fluid-saturated lucite sandstone block and a man-made sandy porous medium. Until now, most of the experimental work performed either used a water-filled fracture borehole model, illuminated with a crosshole acquisition geometry (Zhu & Toksöz 2003a, b, 2005), a sand-based model whose layers differ mainly in terms of saturating fluids (Chen & Mu 2005; Block & Harris 2006; Liu et al. 2008) or a rock sample (real or synthetic) immersed in a water tank (Zhu et al. 2008, 2016; Charara et al. 2009; Schakel et al. 2011a, b, 2012; Smeulders et al. 2014; Peng et al. 2016, 2017, 2018, 2019) in which the seismoelectromagnetic phenomena are generated at a fluid-rock interface.

In this work, we perform a quantitative study of the SE phenomena based on the above mentioned electrokinetic theory (Pride 1994; Pride & Haartsen 1996), with a particular interest in the EM waves created at interfaces between geological materials. For that we have implemented a new laboratory set-up, designed to ease the comparisons between experimental measurements and numerical modelling results, enabling the simultaneous recording of the co-seismic and the interface responses. We use the computer code of Garambois & Dietrich (2002) for the numerical simulations, with its extensions to compute the electric potentials (Dietrich et al. 2018) and further adaptations to account for the source-receiver geometry and characteristics of our experiments. We first consider a homogeneous medium to investigate the limitations of the electrokinetic theory and the modelling code when the latter are used to predict the experimental data. Subsequently, we insert a thin layer in our sandbox to study the interface-generated EM signals, since thin layers can act as natural amplifiers for these signals (Pride & Garambois 2005; Grobbe & Slob 2016). At last, we compare the experimental data with numerical simulations in both cases: the one with a homogeneous medium and the one in which a thin layer is present.

## 2 THEORETICAL BASIS

The formulation proposed in Pride (1994) and Pride & Haartsen (1996) regarding seismoelectromagnetic effects taking place in fluid-saturated porous media will be used in this study. It assumes that pores are filled with an electrically conducting solution, which results in the electrokinetic coupling between seismic and EM fields. At the microscopic level (i.e. for grain size dimensions), SE phenomena are produced by the time-varying interactions between charges adsorbed on the surface of the grains and the ions moving in the pore fluid.

In their implementation of the theory, Pride & Haartsen (1996) present solutions for point-sources and plane waves in homogeneous isotropic whole spaces, inferring the existence of four wave types in the so-called PSVTM wave system: two longitudinal seismic waves (slow- and fast-*P*), a transverse seismic wave (shear or *S*) and a transverse EM wave. They also derive the boundary conditions that hold at the interface between two distinct porous media. Using their
<table><caption>**Table 1.** Physical parameters adopted for the sand and sandstone materials used in the SE experiments and numerical simulations.</caption><thead><tr><th rowspan="2">Name, symbol, units</th><th colspan="2">Sand</th><th colspan="2">References</th></tr><tr><th>Sand</th><th>Sandstone</th><th>Sand</th><th>Sandstone</th></tr></thead><tbody><tr><td>Porosity (φ) [ ]<sup>a</sup></td><td>0.4</td><td>0.24</td><td>Barrière et al. (2012)</td><td>Measured</td></tr><tr><td>Saturation (S<sub>w</sub>) [ ]<sup>a</sup></td><td></td><td>1</td><td></td><td></td></tr><tr><td>DC permeability (k<sub>0</sub>) [m<sup>2</sup>]</td><td>1 × 10<sup>-11</sup></td><td>7 × 10<sup>-15</sup></td><td>Barrière et al. (2012)</td><td>Lorne et al. (1999)</td></tr><tr><td>Bulk (solid) modulus (K<sub>s</sub>) [Pa]</td><td>3.6 × 10<sup>10</sup></td><td>3.6 × 10<sup>10</sup></td><td></td><td>b</td></tr><tr><td>Shear (solid) modulus (G<sub>s</sub>) [Pa]</td><td>4.4 × 10<sup>10</sup></td><td>4.4 × 10<sup>10</sup></td><td></td><td>b</td></tr><tr><td>Bulk (fluid) modulus (K<sub>f</sub>) [Pa]</td><td>2.5 × 10<sup>9</sup></td><td>2.5 × 10<sup>9</sup></td><td></td><td>b</td></tr><tr><td>Bulk (frame) modulus (K<sub>fr</sub>) [Pa]</td><td>2.38 × 10<sup>7</sup></td><td>3.41 × 10<sup>8</sup></td><td>Pride (2005)</td><td>Seismic velocities</td></tr><tr><td>Shear (frame) modulus (G<sub>fr</sub>) [Pa]</td><td>1.43 × 10<sup>7</sup></td><td>2.10 × 10<sup>8</sup></td><td></td><td>Pride (2005)</td></tr><tr><td>P-wave velocity (v<sub>p</sub>) [m s<sup>-1</sup>]</td><td>1780</td><td>2200</td><td></td><td>From the first breaks</td></tr><tr><td>EM-wave velocity (v<sub>EM</sub>) [m s<sup>-1</sup>]</td><td>5.4 × 10<sup>7</sup></td><td>8.5 × 10<sup>7</sup></td><td></td><td>Pride & Haartsen (1996)</td></tr><tr><td>S-wave velocity (v<sub>S</sub>) [m s<sup>-1</sup>]</td><td>90.4</td><td>336.8</td><td></td><td>Pride & Haartsen (1996)</td></tr><tr><td>Fluid viscosity (η) [Pa·s]</td><td></td><td>0.001</td><td></td><td>b</td></tr><tr><td>Solid density (ρ<sub>s</sub>) [kg m<sup>-3</sup>]</td><td>2650</td><td>2120</td><td>b</td><td>Measured</td></tr><tr><td>Fluid density (ρ<sub>f</sub>) [kg m<sup>-3</sup>]</td><td></td><td>1000</td><td></td><td>b</td></tr><tr><td>Fluid conductivity (σ<sub>f</sub>) [mS m<sup>-1</sup>]</td><td></td><td>1.2–1.8</td><td></td><td>Measured</td></tr><tr><td>Molarity (C<sub>0</sub>) [mmol l<sup>-1</sup>]</td><td></td><td>0.128–0.195</td><td></td><td>Pride (1994): C<sub>0</sub> = σ<sub>f</sub>/(e<sup>2</sup>(b<sub>Na</sub> + b<sub>Cl</sub>)N<sub>A</sub>)</td></tr><tr><td>Surface conductivity (σ<sub>n</sub>) [mS m<sup>-1</sup>]</td><td>0.083</td><td>2</td><td>Revil & Glover (1998)</td><td>Lorne et al. (1999)</td></tr><tr><td>Temperature (T) [K]</td><td></td><td>288.35</td><td></td><td>Measured</td></tr><tr><td>Fluid permittivity (κ<sub>f</sub>) [ ]<sup>a</sup></td><td></td><td>80</td><td></td><td>b</td></tr><tr><td>Solid permittivity (κ<sub>s</sub>) [ ]<sup>a</sup></td><td></td><td>4</td><td></td><td>b</td></tr><tr><td>Tortuosity (α<sub>∞</sub>) [ ]<sup>a</sup></td><td>1.65</td><td>3</td><td>Holzhauer (2015)</td><td>Lorne et al. (1999)</td></tr><tr><td>Pore geometry number (m) [ ]<sup>a</sup></td><td>6.79–12.91</td><td>6</td><td>Pride (1994)</td><td>Haartsen & Pride (1997)</td></tr><tr><td>Zeta potential (ζ) [mV]</td><td></td><td></td><td>Guichet et al. (2003): ζ = -0.0146ln(σ<sub>f</sub>) - 0.0854</td><td>Lorne et al. (1999)</td></tr><tr><td>Electrokinetic coefficient (C<sub>ek</sub>) [μV Pa<sup>-1</sup>]</td><td>2.18–5.78</td><td>0.113–0.218</td><td></td><td>Guichet et al. (2003)</td></tr><tr><td colspan="5"><strong>Constants</strong></td></tr><tr><td>Avogadro's number (N<sub>A</sub>) [mol<sup>-1</sup>]</td><td></td><td>6.02 × 10<sup>23</sup></td><td></td><td></td></tr><tr><td>Elementary charge (e) [C]</td><td></td><td>1.60 × 10<sup>-19</sup></td><td></td><td></td></tr><tr><td>Ionic mobility (b<sub>l=Na=Cl</sub>) [m (sN)<sup>-1</sup>]</td><td></td><td>3 × 10<sup>11</sup></td><td></td><td>Approximation, after Pride (1994), for typical inorganic ion</td></tr><tr><td>Vacuum permittivity (ε<sub>0</sub>) [F m<sup>-1</sup>]</td><td></td><td>8.85 × 10<sup>-12</sup></td><td></td><td></td></tr><tr><td>Vacuum permeability (μ<sub>0</sub>) [H m<sup>-1</sup>]</td><td></td><td>4π × 10<sup>-7</sup></td><td></td><td></td></tr><tr><td>Boltzman constant (k<sub>b</sub>) [J K<sup>-1</sup>]</td><td></td><td>1.38 × 10<sup>-23</sup></td><td></td><td></td></tr></tbody></table>

Relevant constants are also included in the list.

<sup>a</sup>Refers to dimensionless properties.

<sup>b</sup>Refers to standard values.

definition of frequency-dependent phase slowness (s<sub>ξ</sub>) for each ξ wave type, we can calculate the real phase velocity v<sub>ξ</sub> and the intrinsic attenuation α<sub>ξ</sub> (Biot 1956a, b) with:

$$v_{\xi=P,S,EM} = 1/\text{Re}\{s_{\xi=P,S,EM}\}, \quad (1a)$$

$$\alpha_{\xi=P,S,EM} = \omega\text{Im}\{s_{\xi=P,S,EM}\}; \quad (1b)$$

where ω is the pulsation, related to the frequency *f* through ω = 2π*f*. In what follows, after computing v<sub>ξ</sub> and α<sub>ξ</sub> with these equations, the values found will be compared with those obtained from the experimental data. ξ is equal to *P* for fast-compressional waves, *S* for transverse (shear) seismic waves and EM for transverse EM waves.

The parameters which are relevant to the quantitative analysis of SE phenomena in this study are listed in Table 1.

### 3 EXPERIMENTAL SETUP

The experiment takes place in a parallel epipedal box filled with sand as shown in Fig. 1. A piezoelectric transducer (PZT), which acts as the seismic source, is placed against the 24.7 (length) × 11.5 (height) × 11 (width) cm plastic box filled with water-saturated sand. The walls of this box are 5.2 mm thick, except the one against

which the PZT is placed, which has a thickness of 2.5 mm. A comb of electrodes is inserted in the sand to measure the electric potentials (noted *V* hereafter), with the help of a reference electrode positioned in an auxiliary plastic container also filled with sand, but free of all seismic disturbances created by the active source. A laser vibrometer (Bodet et al. 2010; Shen et al. 2022) measures the vertical seismic displacement *u*<sub>y</sub> with the laser beam pointing towards the top of the sand (see Devi et al. 2018, for a similar configuration). In the experiments designed to measure the EM interface response, a water saturated sandstone layer is inserted vertically in the sand, approximately in the middle of the box, as shown in Fig. 1.

In our experiments, when choosing the orientation of the Cartesian coordinate system, we sought to closely reproduce field situations while also accounting for the operational constraints of the laboratory equipment. Therefore, since the sandstone layer is vertically inserted in the sandbox, the *z*-axis (depth in field surveys) was chosen to be perpendicular to the slab, that is, along the horizontal axis parallel to the largest dimension of the sandbox, as shown in Fig. 1. The *x*-axis is parallel to the second horizontal direction, and the *y*-axis follows the vertical direction, pointing downwards. Finally, the origin *O* of the Cartesian system corresponds to a point on the top of the sand as indicated in Fig. 1.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Picture of the experimental set-up, highlighting the main electronic components of the acquisition system connected to the sandbox, and detailing the coordinate system used.</figcaption></figure>

## 3.1 The porous media

Both porous media, sand and sandstone, are kept saturated with distilled water whose fluid conductivity ranges between 1.2 and 1.8 mS m⁻¹. The fluid conductivity was measured before and after each experiment to confirm that the values remained within this range. To ensure full saturation, a vacuum chamber is used to pump the air out of the pores during 2.5 hr. In our protocol, the experiment is set in place and left to stabilize while, at the same time, the amplitudes of seismic and electric signals are monitored during test measurements. The actual experiment begins when two successive measurements show that amplitudes are stable (after a period of about 2 hr).

### 3.1.1 Sand

The main material used in the experiments is a homogeneous and well-sorted sand from Southwest France, mainly composed of quartz (>98 per cent) and whose mean grain diameter is 250 µm (e.g. Barrière *et al.* 2012; Bordes *et al.* 2015; Holzhauer *et al.* 2017).

Concerning the mechanical properties of the sand, the drained bulk $K_{fr}$ and shear $G_{fr}$ moduli (see Table 1) are not measured. They are estimated using the modified Walton's model (Walton 1987), presented in Pride (2005), and adjusted so that the velocities predicted with the moduli correspond to the direct P-wave velocities measured from first break (FB) picking in the experimental data (see Sections 5 and 6).

Regarding the electric properties, Guichet *et al.* (2003) studied a sand whose properties are similar to the one used here, therefore the zeta-potential $\zeta$ was calculated following their empirical law as, indicated in Table 1. Additionally, to account for a non-negligible surface conductivity (see values in Table 1), Guichet's equation (Guichet *et al.* 2003, 2006) for the electrokinetic coupling coefficient $C_{ek}$ is adopted (see eq. A4 in Appendix A).

### 3.1.2 Sandstone layer

The thin layer inserted in the sand is a Vosges sandstone (Eastern France) with dimensions: 0.74 (thickness) × 5.7 (height) × 9.3 (width) cm. In order to estimate the density of this rock, samples from the same material but with different sizes were used. Their weight was determined using a precision balance, their volume was measured with the fluid displacement method, and from that a density of 2120 kg m⁻³ was estimated for the rock, with a standard deviation of 161.7 kg m⁻³. Using the same rock samples, a value of 0.24 was estimated for the porosity via the imbibition method, considering distilled water, with a standard deviation of 2.4 × 10⁻². The fluid conductivity of the sample used in the experiments was maintained constant and equal to that of the fluid saturating the sand through periodic monitoring and adjustment with distilled water. Complementarily, the frame moduli were estimated from the model discussed in Pride (2005) for consolidated media, and adjusted so that the direct P-wave velocities predicted for the sand (estimated as mentioned in the previous section) and for the sandstone layer lead to traveltimes that correspond to the FBs picked in the experimental data.

The sample used here has a fair amount of feldspar, mica, clay, and oxides. It is similar to the Adamswiller sandstone sample studied in David *et al.* (1994) and Lorne *et al.* (1999). Permeability, tortuosity (derived from the definition of formation factor, $F = \alpha_{\infty}/\phi$), surface conductivity, and zeta potential are finally taken from Lorne *et al.* (1999) (see Table 1).

## 3.2 Seismic source

A PZT model Olympus V1012 whose effective diameter is 3.8 cm and nominal frequency is 250 kHz, is used as a source of seismic waves in the experimental setup. It is positioned outside the plastic box containing the sand such that its upper end coincides with the sand surface (see Fig. 1) and with the origin of the Cartesian system (x, y, z) used in this work.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Source signature of the PZT measured in the centre of its surface using a laser vibrometer. (a) Source wavelet of the seismic displacement perpendicular to the surface of the PZT. (b) Amplitude spectrum of the source wavelet.</figcaption></figure>

The source signal is generated as a 250 kHz half-sine burst, with an amplitude of 8 V, using a RIGOL DG4102 waveform generator, sent to an amplifier (Falcon Systems WMA 300) and then to the PZT (see Fig. 2). Using a laser vibrometer, we have verified the spatial homogeneity of the PZT signature on its entire surface. It is clear from the time-series and the amplitude spectrum in Fig. 2 that the early part of the source wavelet has a 180 kHz peak frequency, a half-sine shape, and is followed by a series of lower amplitude oscillations that last until about 65 µs.

## 3.3 Seismic measurements description

A laser vibrometry system was used to measure the particle displacements of the sand. It consists of a OFV-505 head sensor that detects particle movement and a OFV-5000 controller which analyses and transmits the recorded displacement using a DD-300 decoder. The system bandwidth is [0.3 kHz–24 MHz], and the upper bound for the absolute displacement is 75 nm.

In the experiment, measurements are carried out on the surface of the sand in such a way that the displacement component $u_y$ is orthogonal to the axis of the seismic source (y-direction in Fig. 1). To increase the signal-to-noise ratio, a 4 mm² piece of retro-reflective film was placed at each measuring point on top of the sand. The water-level was kept constant and equal or slightly lower than the sand-level in the box to make sure the displacement detected is as close as possible to the sand grain displacement, assuming that the reflective paper has no influence on the measurements at our working frequencies (Hasanian & Lissenden 2017; Shen 2020).

During experiments, the laser vibrometry measurements are performed following a central line along the box's length (z-axis in Fig. 1), covering the ranges $z = [1–10.75]$ cm and $z = [13–22.75]$ cm, with a 0.125 cm spacing between points, resulting in 159 source–receiver offsets. Finally, each seismic trace is stacked 3000 times and recorded during 10 ms.

## 3.4 SE measurements description

In SE studies, the measurements of the electric field are usually obtained from potential differences $\Delta V$ using dipolar arrays of electrodes (Garambois & Dietrich 2001; Bordes 2005; Dupuis 2008;

Holzhauer *et al.* 2017). However, Dietrich *et al.* (2018) and Devi *et al.* (2018) questioned this choice of arrays when the goal is to record EM waves generated at interfaces, because electric dipoles strongly attenuate signals whose arrival times are practically synchronous at all electrodes. Alternatively, they proposed configurations using an odd number of electrodes because such configurations make it possible to estimate the electric potential $V$, which they consider as the ideal choice to measure the interface response. Measuring electric potentials implies that a reference electrode is used, which is often impracticable in field surveys due to noise issues (Devi *et al.* 2018). Conversely, SE laboratory studies that intend to detect the EM interfacial signal generally adopt a reference electrode which is either taken as the common-ground (Zhu *et al.* 1999; Smeulders *et al.* 2014; Zhu *et al.* 2016; Peng *et al.* 2017), as a point inside the experimental setup (Chen & Mu 2005; Zhu & Toksöz 2005; Devi *et al.* 2018), or as an adapted/modified common-ground (Block & Harris 2006).

The measurement of SE potentials $V_{exp}$ are performed with an electrode comb made of 32 stainless steel electrodes with an inter-electrode distance of 0.25 cm. Each electrode has a diameter of 1.1 ± 0.2 mm and a length of 10.8 cm, 5.7 ± 0.2 cm buried in the sand. A total of 80 measurements are carried out, 40 in each side of the box, since the sandstone layer is placed around the middle of the box (see Fig. 1). The electric measuring points cover exactly the same ranges in the z-direction than the seismic measurements performed with the laser. Because the comb has only 32 electrodes, it was moved to get eight additional measurements on each side of the box.

The precise measurement of the potential $V_{exp}$ at each electrode uses the reference electrode inserted in the second container filled with the same sand as the one in the main plastic box (pink container in the bottom-right part of Fig. 1). A copper wire connected to a 10 MΩ resistance connects the two plastic containers, to ensure that the potential in the reference container stays close to zero at all times. Each electrode in the comb is connected to a preamplifier (model SR560) in differential mode with the reference electrode, and a [0.3 kHz–1 MHz] band-pass filter is applied. Finally, a study carried out by Martins-Gomes (2022) confirmed that the 100 MΩ input impedance of the preamplifier is high enough to ensure our
measurements are consistent, in terms of amplitude, with the absolute potential voltages mainly due to the electrokinetic phenomena taking place inside the sandbox.

# 4 NUMERICAL MODELLING TOOL

Numerical solutions of Pride's equations (eqs B1 in Appendix B) have been implemented by various research groups using different numerical techniques. Here, our numerical simulations of the laboratory measurements are carried out with the numerical code developed by Garambois & Dietrich (2002), which is capable of computing the mechanical and EM fields, as well as their coupling, due to a point source exciting a homogeneous poro-elastic layered medium. It also allows a straightforward modification of the reflection and transmission coefficients, which is a convenient property to enhance small amplitude, hard to detect events, such as the EM interface response.

## 4.1 Extension of the existing numerical code

To be consistent with the geometry of the experiment, the $E_z$ component of the electric field is needed for our comparisons. Because this component was not considered in the original code of Garambois & Dietrich (2002), the calculation of $E_z$ is detailed in Appendix B and was implemented in the numerical simulation code.

Additionally, since the electric potentials $V_{exp}$ are measured in the experiments, this quantity was also calculated with the procedure developed in Dietrich et al. (2018). In this approach, the electric field **E** has no transverse component for sources that are axially symmetric (explosion, or point forces oriented perpendicularly to the layering), implying that the solution for **V** can be derived from the cylindrical radial component of the electric field using a double integral. Thus, in cylindrical coordinates:

$$V(\rho, z, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-i\omega t} d\omega \int_{0}^{\infty} \hat{E}_{\rho}(\omega, k, z) J_{0}(k\rho) dk, \quad (2)$$

where $J_0$ is the zeroth order Bessel function, $\hat{E}_{\rho}$ is the radial component of the electric field in the Fourier-Hankel domain, k is the horizontal wavenumber, and $\rho$ is the radial distance from the z-axis.

## 4.2 Particularities of the numerical modelling of the experiments

### 4.2.1 Numerical modelling of the PZT source

The signal measured in the centre of the PZT surface, using the laser vibrometer, is used as the input seismic source in the numerical simulations of the experiments. Since the surface area of the PZT source is relatively large compared to the dimensions of the experimental setup, the seismic source is discretized and modelled as a group of $N_{sources}$ synchronous explosive point sources, distributed along the line defined by the coordinates $x = 0, z = 0$ and $y \in [0, 3.8]$ cm. It corresponds to a y-oriented linear distribution of point sources equally spaced along the PZT diameter. Following a comparative study between the numerical results output by a point source at $x = 0, z = 0$ and $y = 1.9$ cm, a y-oriented line of sources, two perpendicular source lines (x and y-oriented) and a 2-D distribution of sources along the whole PZT surface, we have found that the y-oriented linear distribution provided the best compromise between an accurate solution and computational cost for the range

of offsets we explore when comparing numerical simulations and experimental data.

It follows that, for a point on top of the sand, the numerically computed vertical displacement $u_{y-num}$, accounting for contributions from all sources, will be obtained from a summation of the displacements generated by each individual source $u_{y-n}$, weighted by the number sources:

$$u_{y-num}(x = 0, y = 0, z, t) = \sum_{n=1}^{N_{sources}} u_{y-n}(x = 0, y = 0, z, t) \delta y, \quad (3)$$

with $\delta y = l_{PZT}/N_{sources}$ the spatial discretization step, and $l_{PZT} = 3.8$ cm the diameter of the PZT. After a numerical parametric study on the influence of the value of $N_{sources}$ on the total displacement $u_{y-num}$ (Martins-Gomes 2022), $N_{sources} = 80$ was deemed appropriate, thus the value for $\delta y$ used in eq. (3) was $0.475 \times 10^{-3}$ m, which corresponds approximately to $\lambda/20$, where $\lambda$ is the P wavelength calculated for peak frequencies around 180 kHz in the sand.

### 4.2.2 Numerical modelling of the experimental data recorded by elongated electrodes buried in the sand

Another correction is required to correctly evaluate the electric potential using the computer code because the experimental measurements involve relatively long electrodes that do not satisfy the point receiver assumption, which is intrinsic to the simulation program. Therefore, recalling that the length of the buried part of the electrode is $l_{elec} = 5.7$ cm (see Section 3.4), the effective values of the potential $V_{num}$, at a specified z-position, must be estimated by integrating the electric potential along the buried part of the electrode at that location, that is

$$V_{num}(x = 0, z, t) = \int_{0}^{l_{elec}} V_{num}(x = 0, y, z, t) dy. \quad (4)$$

To formulate a discrete version of eq. (4), each electrode is discretized using $N_{elec} = 120$ points. This criterion was chosen to ensure that $\delta y$, the spatial sampling step of the electrodes, is equal to the discretization used to represent the PZT source (see eq. 3). Finally, the discrete form of eq. (4), accounting for the $N_{sources} = 80$ seismic sources is:

$$V_{num}(x = 0, z, t) = \sum_{p=1}^{N_{elec}} \sum_{n=1}^{N_{sources}} V_{n,p}(x = 0, z, t)(\delta y)^2, \quad (5)$$

where $V_{n,p}$ is the electric potential signal computed at a given point p of an electrode, and generated by an elementary source n in the PZT discretization.

# 5 SEISMOELECTRICS IN A HOMOGENOUS MEDIUM

## 5.1 Experimental data

Fig. 3 displays the y-oriented displacement traces $u_{y-exp}$ and the electrograms of potential $V_{exp}$, measured during the experiment performed with the homogeneous saturated sand model, after low-pass filtering using a 1 MHz cut-off frequency. The gap between offsets 10.75 and 13 cm was intentionally planned to anticipate the insertion of a sandstone layer (Section 6). Concerning the displacement traces (Fig. 3a), the dominant direct compressional P wave is clearly observed with its linear trend of arrival times (dashed red line); it is followed by reflected P waves, coming from all sides of the sandbox
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Experimental seismoelectric response of the homogeneous sand model, as a function of time and source-receiver offset. (a) Measured particle displacement traces $u_{y-exp}$. (b) Measured electric potential $V_{exp}$. The front and back walls of the sandbox are shown for convenience. The direct wave emitted by the seismic source and its coseismic signature are indicated by the red dashed lines. Some of their reflections in the sandbox are indicated by dashed blue lines.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Amplitude versus source-receiver offset distributions of the experimental seismoelectric response obtained with the homogeneous sand model. (a) Representation of the maxima and minima of the measured particle displacement $u_{y-exp}$. (b) Same for the recorded electric potential $V_{exp}$. The plotted values are associated with the first arrivals of the two wavefields. The dashed curves superimposed on the experimental data (stars) are mathematical representations of the amplitude decay as a function of offset (see text). (c) Peak frequency $f_{peak}$ of the first arrival of the seismic displacement $u_{y-exp}$ as a function of offset, calculated from the time delay separating the minimum and maximum amplitudes (orange circles), and from the amplitude spectrum (black circles). (d) Same for the recorded electric potential $V_{exp}$.</figcaption></figure>

and from the sand-air interface, which are highlighted by the blue dashed lines. Except for the *P* wave reflected on the back wall (opposite to the PZT source), all other reflected *P* waves are difficult to interpret as to their exact origin.

In the electric measurements (Fig. 3b), the event near $t = 0$ is EM noise inherent to PZT sources, and the synchronous events around 0.05 and 0.12 ms are artefacts due to the presence of the electrode comb. One can easily notice the coseismic waves travelling in the sand with the same velocity $v_p$ as the seismic *P* wave. The first-break (FB) of the direct coseismic event and its reflections on the back of the box are indicated by the blue dashed lines in Fig. 3(b).

The *P*-wave velocity of the medium $v_p$ was estimated to be about $1782 \text{ m s}^{-1}$ from the linear regression of the FBs picked on the electric potential data. The FBs of both $u_{y-exp}$ and $V_{exp}$ data sets in Fig. 3 were then consistently adjusted so that their values match exactly the $v_p$ velocity estimated for the sand. This correction will be very pertinent when comparing experimental and numerical data, to ensure that the simulated waveforms match the arrival times of their experimental counterparts.

All discrepancies observed between the seismic and electric responses, notably the absence of reflected waves in Fig. 3(b), are explained by the polarization of the various arrivals and the difference in acquisition geometry and receiver size. The coseismic signature of the direct *P* wave is essentially propagating along the *z*-axis, therefore the inline electrode array will favour its detection. By contrast, the reflections coming from the laterals and bottom of the sandbox reach the electrode array under oblique incidence,

but more importantly, the reflections from the side walls tend to cancel each other due to the symmetry of the electrode array, positioned in the middle of the sandbox ($x \simeq 0$) and oriented along the *z*-axis. Additionally, while the seismic measurements record the *y*-oriented mechanical disturbances in a considerably small region, with respect to the dimensions of the experiment, the electrode, due to its size, measures the electric potential averaged throughout its length. This averaging will attenuate small amplitude and short period events that exist as localized disturbances inside the sandbox and, therefore, only reach a small region of the electrode.

In order to compare the amplitude decay of the recorded SE signals as a function of source-receiver offset with the predictions of Pride's theory, the first minima and maxima of $u_{y-exp}$ and $V_{exp}$ were selected and displayed versus offsets in Figs 4(a) and (b). Subsequently, to extract information about the geometrical spreading and the intrinsic attenuation of the direct *P* waves, and its corresponding coseismic arrival, the Gauss-Newton algorithm was used to fit functions of the form $u_0 e^{-\alpha r} \cdot r^{-g}$ and $V_0 e^{-\alpha r} \cdot r^{-g}$ to the data in Figs 4(a) and (b), respectively. In these expressions, *r* is the radial source-receiver distance at the top of the sand (distance *r* between *O* and the *z*-offset of each measurement point) and $\alpha_p$ is the intrinsic attenuation of the *P* wave, as predicted from Biot's theory for the range of frequencies studied and the parameters displayed in Table 1. The parameters estimated by the curve fitting are $u_0$ and $V_0$, which represent the amplitudes, and *g*, corresponding to the geometrical spreading coefficient. To be consistent with the frequency range [150–300] kHz of the experimental seismic data, the curve
fit considered three theoretical predictions for $\alpha_p$ (calculated with eq. 1b), corresponding to the frequencies 150, 225 and 300 kHz, respectively: 0.060, 0.075 and 0.089 m$^{-1}$. In Fig. 4(a) the difference between the three curves is so small that only that corresponding to 225 kHz is shown.

For the geometrical spreading, assuming a pure 3-D seismic wave propagation in the far-field, $g = 1$, whereas for a 2-D propagation, $g = 0.5$, and for a 1-D propagation, $g = 0$ (Aki & Richards 2002, chap. 4). In Fig. 4(a), the curve fitting of the experimental maxima and minima estimated $g \simeq 1$, which suggests a 3-D wave propagation. Conversely, the curve fitting of $V_{\text{exp}}$ in Fig. 4(b) estimated values close to zero for the geometrical spreading coefficient. This means that the recorded coseismic signal has a near 1-D propagation, likely due to the use of a long electrode positioned at the centre of the box combined with the large surface of the PZT source. Note that such conditions favour the measurement of coseismic waves associated with z-oriented seismic wavefields, since they have more energy than their x and y counterparts.

Focusing now on the frequency content of the recorded seismic signals, we present in Figs 4(c) and (d) the evolution of the peak frequency of the first arrivals ($f_{\text{peak}}$) with increasing offset. In black, $f_{\text{peak}}$ was estimated from the magnitude spectra obtained with a fast Fourier transform (FFT) computed in a 0.011 ms window starting with the FB. In orange, $f_{\text{peak}}$ was computed supposing that the first arrival has a sinusoidal shape passing through the first through and peak, and therefore, its frequency is the inverse of the maximum-to-minimum time lag multiplied by two (i.e. the lag corresponds to the half-period).

No matter how $f_{\text{peak}}$ is calculated, for $u_{y-\text{exp}}$, its value increases with offset, which is an uncommon behaviour that is attributed to interferences of the P wavefields emerging from all the surface elements of the PZT source, as discussed in Fig. S2. Indeed, the summation of the various source contributions is more likely to destroy the high frequency components of the signal at short distances from the source because the seismic path difference is more pronounced at near offsets than at far offsets.

At last, Fig. 4(d) shows the peak frequency of the electric potential data. One can note that the $f_{\text{peak}}$ estimated with the minima-maxima methodology, orange dots in Fig. 4(d), displays an almost constant trend, with values varying mainly in the range [200, 230] kHz. Moreover, the $f_{\text{peak}}$ values estimated with a FFT (black dots) increase at offsets lower than about 8 cm and afterwards they reach a plateau, with values around 180 kHz. This increase in frequency is probably due to wave interferences at offsets smaller than 8 cm, and to the influence of the EM radiation intrinsic to the PZT source. Finally, the use of a long electrode as a receiver can possibly explain that the $f_{\text{peak}}$ values do not vary much for offsets greater than 8 cm.

## 5.2 Comparison between experiments and theory

The results of the numerical simulation of the experiments will now be presented for the homogeneous sand medium. Both amplitudes and waveforms of the displacement and electric potential wavefields obtained experimentally are compared with synthetic traces based on Pride's theory for the input parameters listed in Table 1, and fluid conductivity $\sigma_f = 1.31 \text{ mS m}^{-1}$, measured right after the measurements. All numerical calculations account for the actual size of the seismic source and the length of the electrodes buried in the sand.

Figs 5(a) and (b) display comparisons of the experimental and simulated seismograms ($u_y$) and electrograms ($V$), respectively. Four

specific traces were chosen to emphasize the overall good agreement between real and synthetic waveforms.

The seismic displacement section in Fig. 5(a) shows that the simulated traces correctly reproduce the first arrivals of the experimental records in terms of waveform and frequency content. However, because none of the sandbox's walls are considered in the modelling, the reflections noticed in Fig. 3(b) are not seen in the simulations. Regarding the $V$ electrograms in Fig. 5(b), they also show a good match between the numerical and experimental first arrivals, provided a 0.44 correction factor is applied to the amplitudes.

Fig. 6 presents the minimum and maximum amplitudes of the first arrivals, obtained from experiments and from simulations, for the displacement field $u_y$ and for the electric potential $V$. Fig. 6(a), for $u_y$, shows a fair correspondance between numerical and experimental data, and, as in Fig. 4(a), only offsets larger than 8.5 cm are considered. Fig. 6(b), for $V$, shows a good agreement between experimental and synthetic traces at offsets larger than 8 cm, when near-field effects and/or interferences due to the extended source are less significant. It is important to note, however, that the amplitudes of the simulated wavefields were scaled to obtain a good match with the experimental seismic data, which led to slightly over-amplified $V_{\text{num}}$ fields. For that reason, the amplitudes of $V_{\text{num}}$ were multiplied by a factor of 0.44 in Fig. 6(b) to adjust them with the amplitudes of $V_{\text{exp}}$.

Regarding the need for a 0.44 factor to correct the amplitudes of the synthetic data, Block & Harris (2006) and Schakel et al. (2011a) have noticed that for the same range of fluid conductivities used here, there is a discrepancy between Pride's theory and experimental measurements which they attribute to surface conductivity effects not accounted for in the SE coupling. In here, even though we have introduced the surface conductivity through the $C_{\text{ek}}$ formulation of Guichet et al. (2003), the choice of the surface conductivity values in Table 1 may not correspond to the actual values of the media, or this formulation of $C_{\text{ek}}$ may not be adapted to our case, probably due to the range of fluid conductivities we consider.

In summary, the comparison of experimental data with their numerical simulations is generally a difficult task given all the constraints related to each of the two approaches. Nevertheless, and except for the limitations mentioned above, the overall good agreement observed for (the early part of) the time-series obtained experimentally and numerically in a homogeneous sand model reinforces Pride's electrokinetic theory. Our quantitative analysis goes beyond the qualitative results of Devi et al. (2018) and may pave the way for yet-to-be developed inversion algorithms.

# 6 SEISMOELECTRICS IN A LAYERED MEDIUM

To study the SE signals generated at an interface between two porous media, a layered medium was designed by inserting a 0.74-cm-thick sandstone layer in the middle of the sandbox, around $z = 11.88$ cm. Table 1 lists the physical properties of the sandstone; since the P-wave velocity was experimentally estimated as $v_P = 2200 \text{ m s}^{-1}$, the layer is thinner than the P-wavelength ($\lambda$) in the sand medium, which ranges between 0.6 and 1.2 cm for the range of frequencies considered (150–300 kHz). Both seismic and electric data were acquired with the geometry displayed in Fig. 1. The first-breaks were again time-adjusted in the seismic and electric sections to match the theoretical arrival times calculated from the velocities experimentally estimated for sand and sandstone. Moreover, data were low-pass filtered using a cut-off frequency of 1 MHz.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison between experimental (black) and numerical (orange) seismograms and electrograms for the homogeneous sand model. (a) Seismic displacement traces $u_{y-exp}$ and $u_{y-num}$. (b) Electric potential traces $V_{exp}$ and $V_{num}$. The $V_{num}$ time-series were scaled by a factor of 0.44 to adjust them to the $V_{exp}$ traces obtained experimentally. The 4 panels displayed in (a) and (b) are close-up views of the first arrivals at the offsets specified in each panel.</figcaption></figure>

## 6.1 Experimental data

Experimental measurements of the displacement $u_{y-exp}$ and of the electric potential $V_{exp}$ are displayed in Figs 7(a) and (b), respectively. As expected, in both cases the main events are similar to the ones

observed in the homogeneous case (Fig. 3), with the notable addition of the P-wave reflection on the sandstone layer (highlighted with the magenta line in Fig. 7a). Furthermore, on the electrogram section, we can observe two new arrivals generated by the presence of a
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Amplitude versus offset distributions of the experimentally measured (black asterisks) and numerically computed (red circles) seismoelectric time-series obtained with the homogeneous sand model. (a) Representation of the maxima and minima of the seismic particle displacements $u_{y-exp}$ and $u_{y-num}$. (b) Same for the electric potentials $V_{exp}$ and $V_{num}$. The $V_{num}$ time-series were scaled by a factor of 0.44.</figcaption></figure>

sandstone layer: the coseismic event created by the P-wave reflection on the layer (magenta) and the EM-wave (green line) emerging from the interface with near-infinite apparent velocity along the z-axis. Fig. S2(a) displays an enlargement of the EM interfacial wave between offsets 13 and 17.75 cm, to highlight the regularity of its decay. In addition, Fig. S2(b) shows the similarity between the peak-frequencies of the interface response and of the direct coseismic wave (Fig. 4d).

In Fig. 8(a), the maximum and minimum amplitudes of the first coseismic arrival are compared with those of the homogeneous sand experiment, shown in Fig. 4(b). There is a good agreement for the traces located at near offsets, before the sandstone layer, but it appears that the amplitudes of the traces located beyond the sandstone layer are somewhat smaller, most likely because of a loss of energy due to the reflections and attenuation in the sandstone layer.

Fig. 8(b) displays the absolute maximum amplitude versus offsets of the EM sandstone interface response $V_{exp-EM}$. The absolute maxima were deemed more appropriate in this representation because of the amplitude jitter of the EM interface response caused by strong interference with coseismic events close to the sandstone layer (the most affected traces are represented by cyan-coloured diamond markers in Fig. 8b). The most striking feature of the $V_{exp-EM}$ amplitudes in Fig. 8(b) is the polarity change of the EM wave on either side of the sandstone layer. This is a known feature of EM interface responses, that has been reported in all numerical simulations (e.g. Garambois & Dietrich 2002; Haines & Pride 2006) and in most experiments (Dupuis et al. 2009; Schakel et al. 2011b; Peng et al. 2017; Liu et al. 2018).

Another important characteristic of the data in Fig. 8(b) is the decay rate of the $V_{exp-EM}$ amplitudes vs offset: it is clearly different on each side of the sandstone layer. As in the homogeneous case, we superimposed the Gauss-Newton curve fits based on the intrinsic attenuation of the EM wave $\alpha_{EM} = -0.00029 \text{ m}^{-1}$ (for a 225 kHz frequency), as predicted by Pride's theory (eq. 1b) for the sand parameters listed in Table 1. The curves fitted to the data in Fig. 8(b) show a slower decay, with a $r^{-0.6}$ dependence, for the EM wave decay on the source side of the sandstone layer, and a $r^{-1.1}$ dependence beyond the sandstone layer. This difference in decay rate can also be noticed in Fig. 7(b) and it may, again, be explained by constructive and destructive wave interferences. Indeed, for offsets between the

source and the sandstone layer, the backscattered $V_{exp-EM}$ traces are contaminated by the tail of the coseismic event, whereas beyond the sandstone layer, the forward-scattered EM wave is the first arrival, therefore free from any perturbation.

## 6.2 Numerical modelling

As in the homogeneous case, we calculated the theoretical SE response of the medium considered, this time with a sandstone layer. The physical parameters used are, as before, those given in Table 1 and the fluid conductivity considered is the same one used in the homogeneous case, namely $\sigma_f = 1.31 \text{ mS m}^{-1}$. Because the difference between both experiments is mainly due to the EM interface response, we will focus our attention on the analysis of electric potentials and electric fields.

The electric potential electrograms obtained from the experiment $V_{exp}$ and from the simulation $V_{num}$ are superimposed in Fig. 9(a), where $V_{num}$ was multiplied by the same 0.44 amplitude correction factor used in the homogeneous case. As for the coseismic wavefields, the waveform of the numerically modelled EM interfacial response agrees well with the measurements. In general, it can be observed in Fig. 9(a) that the polarity, the decay rate, and the amplitudes of the $V_{exp-EM}$ event are fully consistent with the corresponding numerical simulations beyond the location of the sandstone slab, where the EM wave is not strongly affected by coseismic waves, except for traces measured near the interface.

An alternative option to study the interface response, experimentally and numerically, is to analyse the vector electric field $E$. In our laboratory experiments, $E_z$ can be obtained from the difference of the electric potential measured at two neighbouring electrodes. Accordingly, Fig. 9(b) displays the complete electrograms of the experimental and simulated $E_z$. From the trace superposition, it becomes evident that, again, there is a good agreement between experiments and simulations, for both the coseismic field and the EM interface response even if the latter is barely visible on the $E_z$ traces, in particular those measured in the laboratory which are strongly corrupted by noise.

Indeed, the most important information contained in the $E_z$ electrogram of Fig. 9(b) is that the amplitudes of the interface response, relative to the coseismic signals, are significantly smaller when compared with the electric potential records. This was an expected
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Experimental seismoelectric response of the model with a sandstone layer in the middle of the sandbox, as a function of time and source-receiver offset. (a) Measured particle displacement traces $u_y-\exp$. (b) Measured electric potential $V_{\exp}$. The front and back walls of the sandbox and the sandstone layer are shown for convenience. The direct wave emitted by the seismic source and its coseismic signature are indicated by red dashed lines. Some of their reflections in the sandbox are indicated by red dash-dotted lines. The main reflections generated by the sandstone layer are indicated with the dashed magenta lines. In (b), the EM event created by the presence of the sandstone layer is highlighted by the dashed green line.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Amplitude versus offset distributions of the seismoelectric potential. (a) Representation of the maxima and minima of the coseismic signature of the direct P-wave obtained experimentally in the homogeneous sand model (empty circles) and in the sand-sandstone layered model (filled circles). (b) Representation of the absolute maxima of the EM interfacial wave generated by the sandstone layer, when recorded experimentally (black stars) or when computed numerically (empty green circles). The dashed red lines are curve fits based on the theoretical EM amplitude decay predicted by Pride's theory: for offsets smaller than 11 cm, the curve equation is $0.28e^{-0.00029 \cdot r} \cdot r^{-0.6}$ while for offsets larger than 12 cm, the curve equation is $-0.39e^{-0.00029 \cdot r} \cdot r^{-1.1}$. The points nearest to the layer, highlighted with cyan diamonds, were not taken into account in the curve fit.</figcaption></figure>

result, since the electric field is calculated as a difference of the electric potentials at two adjacent traces, a procedure that strongly impacts any signal arriving simultaneously at all electrodes with approximately the same amplitudes, as already noticed by Devi et al. (2018). The electric potential traces of Fig. 9(a) show that the EM response can be followed visually over a distance of about 10 cm from the sandstone layer on the far side of the source. On the other hand, the EM response can only be followed up to 1–2 cm from the interface when looking at the electric field records in Fig. 9(b). Eventually, these observations highlight the fact that electric potential measurements are far superior to recordings when it comes to study the interface response, as formally demonstrated by Dietrich et al. (2018).

Furthermore, varying the salinity (molarity) of the fluids saturating the sand and sandstone in the numerical simulations, within the range given in Table 1, did not lead to significant changes, even though it has been shown that the interface response is very sensitive to salinity contrasts (Garambois & Dietrich 2002).

## 6.3 Interpreting the interface response

Seeking to understand both the physical phenomena involved in the generation of EM-waves at the sand-sandstone interface, and the physical factors controlling their amplitude decay as a function of distance from the interface, we ran additional numerical tests, with the same layered medium as in the previous simulation, but using, for simplicity, a single seismic point source placed at the origin of the coordinate system. Point receivers measuring the displacement $u_{\text{num}}$ and the electric potential $V_{\text{num}}$ are distributed in the x–z plane, in the z range [1 to 23] cm and the x range [−8 to 8] cm. In this particular simulation, all P-to-EM conversions were amplified 1000 times (i.e. the respective reflection and transmission coefficients were artificially increased) in order to favour the interface response over the coseismic signals, so that they become the most prominent events in the electrograms.

At the dominant frequency of the seismic source, the wavelength of the EM interface response is about three orders of magnitude larger than the dimensions of the sandbox model considered. Therefore, the PZT source and the electrodes may be considered as point

source and point receivers with no loss of generality. This remark is important for the validity of the following simulations.

In Fig. 10, we display four snapshots of these simulations. The four columns displayed for each selected time show, from left to right: the displacement $u_{z-num}$, the vertical electric field $E_{z-num}$, the electric potential $V_{num}$, and the electric potential evaluated along the yellow $x=0$ line in the middle of the model. We will focus on what happens inside the sandstone layer and at its boundaries with the sand, which are represented with the two horizontal green lines at z-offsets 11.88 and 12.62 cm.

The four time values associated with the snapshots in Fig. 10 were carefully chosen to image the passage of the P wave front through the sandstone slab. The first one (67.97 µs) corresponds to the arrival of the first trough of the P wavelet at the sand-sandstone interface, at offset 11.88 cm. The second one (70.01 µs) is when the main peak reaches that interface. Likewise, the third (71.42 µs) and fourth (73.4 µs) times are respectively associated with the arrivals of the main trough and peak of the P wavelet at the sandstone-sand boundary, at offset 12.62 cm.

Regarding $u_{z-num}$ in the first column of Fig. 10, we observe that the perfectly spherical P wave front originating from the point source undergoes only minor changes when crossing the sand-sandstone interface, since the P-wave velocities in the saturated sand and sandstone are relatively close (1780 m s⁻¹, 2200 m s⁻¹). As in the previous simulations, the source time function mimics the waveform of the PZT source depicted in Fig. 2, with its rather long duration and secondary wiggles. The full complexity of the source's signature can be observed in the $u_{z-num}$ and $E_{z-num}$ snapshots of Fig. 10.

A remarkable feature seen in Fig. 10, at time 67.97 µs, is that an electric charge separation occurs across the interface as soon as the P wave reaches the sandstone (Pride & Haartsen 1996). This leads to a discontinuity in the z-profile of the electric potential. This charge imbalance can be readily noticed in the second, third, and fourth columns of Fig. 10, and manifests itself as a dipole-like source. Moreover, the existence of a positive lobe and a negative lobe in $V_{num}$, as well as two negative lobes in $E_{z-num}$, indicates that the charge separation takes place preferentially along the z-axis. In the second snapshot (70.01 µs), the second lobe and main peak
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Comparisons of electrograms obtained experimentally (black) and numerically (orange) for the sandstone layer model. (a) Electric potential section. (b) z-component of the electric field. As for the homogeneous sand model, the time-series obtained by numerical simulations were scaled by a factor of 0.44 to obtain a good match with the experimental traces. The 4 panels displayed in (a) and (b) are close-up views of the interface response at the offsets specified in each panel.</figcaption></figure>

of the *P* wave touches the sand-sandstone interface, with opposite sign and larger amplitude, and generates a new dipole-like charge separation, sign-reversed with respect to the previous one. These two successively generated dipoles are respectively responsible for

the first and second lobes of the EM wave, as seen in the *V*<sub>num</sub> time-series of Fig. 11, calculated with virtual point-receivers at positions given by the magenta triangles in Fig. 10. The other two time snapshots of Fig. 10 correspond to the times when the two main lobes
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Snapshots of the seismic particle displacement $u_{z-num}$, electric field $E_{z-num}$, and electric potential $V_{num}$ (from left to right) computed for a seismic point source at four successive propagation times. The curves on the rightmost column are profiles of the electric potential calculated along the dotted yellow line in the $V_{num}$ panels. The location of the sandstone layer is given the green lines. The white lines in the $E_{z-num}$ and $V_{num}$ panels are iso-value curves. The two magenta triangles in the $V_{num}$ snapshots indicate the position of virtual receivers whose records are shown in Fig. 11.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Electric potential time-series $V_{\text{num}}$ at source-receiver offsets 10.5 and 14 cm indicated by the magenta triangles in Fig. 10. (a) Trace computed on the side near the source. (b) Trace computed on the side far from the source. The symbols superimposed on the two curves mark the four times selected for the snapshots in Fig. 10.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Snapshots of the seismic particle displacement $u_{z-num}$, electric field $E_{z-num}$ and electric potential $V_{num}$ (from left to right) at time $t = 67.97$ µs. The simulation results presented account for the actual size of the PZT seismic source exciting the layered sand-sandstone model. The location of the sandstone layer is indicated by the green lines. The white lines in the $E_{z-num}$ and $V_{num}$ panels are iso-value curves.</figcaption></figure>

of the direct *P* wave hit the sandstone-sand (second) interface, also creating a charge separation at that discontinuity, and consequently, EM waves. These waves interact with the former EM-signals generated almost simultaneously at the first interface, because of the relatively long time duration of the source wavelet. These EM wave interferences dictate the shape of the electric potential signals displayed in Fig. 11 which end with the excitation of the mechanical waves, once the waves have completely crossed the sandstone layer. Eventually, because the thickness of the layer (0.74 cm) is comparable or smaller to the wavelengths ($\lambda \simeq 0.6$ to 1.2 cm), in both the experimental and numerical cases, it is difficult to distinguish the individual EM waves associated with each interface in the traces shown in Fig. 11. The first two lobes can be assigned to the first boundary, but the late arrivals are interfering EM waves generated at both interfaces.

In order to understand the effect of an extended seismic source (as opposed to a point source) on the global pattern of the EM interface response, we repeated the previous simulations (Fig. 10) by using five seismic point sources, equally spaced along the *y*-axis along the PZT diameter $l_{PZT}$, for a better analogy with the experiments. Fig. 12 shows a single snapshot of this simulation at 67.97 µs, that can be compared to the first row in Fig. 10. We observe that the extended seismic source generates a *P* wave front characterized by slightly lower frequencies than those of the point source because of

the implied summation of contributions. We also observe a spatially wider dipolar pattern in $V_{num}$ due to the enlarged first Fresnel Zone radius (Thompson 1936).

## 7 CONCLUSION

An experimental setup was proposed as a laboratory analogy to SE field experiments, allowing in particular to measure both the coseismic and the interface responses simultaneously. The experimental measurements were confronted with numerical simulations based on the electrokinetic theory of Pride (1994) by using the material properties, the data acquisition geometry, and the source time function of the experimental setup. The comparisons showed that the measured and simulated time-series have very consistent amplitudes, only differing by a factor of about 0.44 when comparing the electric time-series, while in terms of waveforms they agree satisfactorily. The same conclusions can be drawn when a thin layer is present, with an almost perfect match between waveforms. In the future, to further improve the results of the simulations, and consequently our understanding of Pride's theory, the electrokinetic coupling coefficient should be estimated for the whole range of frequencies we have considered here, using a methodology similar to Zhu & Nafi Toksöz (2013). Additionally, because the thin
sandstone layer contains a small percentage of clay, adapting the current numerical code to account for the surface conductivity of clay minerals or working with an electrokinetic theory adapted to the sandstone we use, may be a way to obtain numerical simulations that are more compatible with the experimental data, in terms of the physical phenomena happening in the interface between solid and liquid phases.

Both the experimental and the numerical data confirm that measuring electric potentials is more suited than recording electric fields (i.e. using dipolar arrays) for the detection and interpretation of EM waves generated at interfaces, as suggested by the analytical developments and numerical simulations of Dietrich *et al.* (2018) and Devi *et al.* (2018). As an alternative to the multi-electrode configuration proposed in Dietrich *et al.* (2018), the transposition to field surveys of our sandbox electrode assembly would need very long shielded cables to position the reference electrode very far from the measuring electrodes. Failing that, the EM noise picked up by the electrodes would probably not completely cancel out. Another option would be to properly shield the reference electrode to prevent the interfacial EM wave from reaching it.

Finally, we noticed that in our electric potential records, the amplitude of the interface response has decreased by almost 90 per cent at a distance of 10 *P* wavelengths from the interface. Therefore, in a hypothetical scaled up field survey, involving similar sand and sandstone formations with both source and receivers located at the surface, and using a *f*<sub>peak</sub> = 50 Hz seismic source (λ ≃ 36 m) the EM wave generated at a depth of 360 m would appear as a very weak signal. However, since contrasts in fluid conductivity have a significant impact on the amplitudes of the EM interface response (Garambois & Dietrich 2002), and were not tested in our experiments, we expect that in a real situation, at least the first tens of metres could be reasonably investigated by the SE method. This must be confirmed with future field surveys that use the electric potential rather than the electric field to properly capture the EM interface responses generated at depth.

## DATA AND SOFTWARE AVAILABILITY

Both the data and the software underlying this work can be shared upon reasonable request.

## ACKNOWLEDGMENTS

This work was supported by the Excellence Initiative of Université de Pau et des Pays de l'Adour – I-Site E2S UPPA (Project CHICKEA). We would like to express our sincere appreciation to the two reviewers, one that chose to remain anonymous, and Damien Jougnot, whose insightful comments and suggestions helped improve this paper.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

### suppl.data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## REFERENCES

* Aki, K. & Richards, P., 2002. *Quantitative Seismology, Geology Seismology, 2nd edn*, University Science Books.

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

* Barrière, J., Bordes, C., Brito, D. & Sénéchal, P., 2012. Laboratory monitoring of P waves in partially saturated sand, *Geophys. J. Int.*, **191**, 1152–1170.

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

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

* Block, G.I. & Harris, J.G., 2006. Conductivity dependence of seismoelectric wave phenomena in fluid-saturated sediments, *J. geophys. Res.*, **111**(B1), doi:10.1029/2005JB003798.

* Bodet, L., Jacob, X., Tournat, V., Mourgues, R. & Gusev, V., 2010. Elasticity profile of an unconsolidated granular medium inferred from guided waves: Toward acoustic monitoring of analogue models, *Tectonophysics*, **496**(1), 99–104.

* Bordes, C., 2005. Etude expérimentale des phénomènes transitoires sismoélectromagnétiques: Mise en oeuvre au Laboratoire Souterrain de Rustrel, Pays d’Apt, *PhD thesis*, Université Joseph-Fourier - Grenoble I.

* Bordes, C., Jouniaux, L., Garambois, S., Dietrich, M., Pozzi, J.-p. & Gaffet, S., 2008. Evidence of the theoretically predicted seismo-magnetic conversion, *Geophys. J. Int.*, **174**(2), 489–504.

* Bordes, C., Sénéchal, P., Barrière, J., Brito, D., Normandin, E. & Jougnot, D., 2015. Impact of water saturation on seismoelectric transfer functions: a laboratory study of coseismic phenomenon, *Geophys. J. Int.*, **200**(3), 1317–1335.

* Charara, M., Zaretsky, I. & Zamora, M., 2009. Seismoelectric modeling of laboratory experiments, in *Proceedings of the SPE Europec/71st EAGE Conference and Exhibition 2009*, Amsterdam, Netherlands, 8–11 June 2009.

* Chen, B. & Mu, Y., 2005. Experimental studies of seismoelectric effects in fluid-saturated porous media, *J. geophys. Eng.*, **2**(3), 222–230.

* Chin, R. C.Y., Hedstrom, G.W. & Thigpen, L., 1984. Matrix methods in synthetic seismograms, *Geophys. J. Int.*, **77**(2), 483–502.

* David, C., Wong, T.-F., Zhu, W. & Zhang, J., 1994. Laboratory measurement of compaction-induced permeability change in porous rocks: implications for the generation and maintenance of pore pressure excess in the crust, *Pure appl. Geophys.*, **143**(1–3), 425–456.

* Devi, M.S., Garambois, S., Brito, D., Dietrich, M., Poydenot, V. & Bordes, C., 2018. A novel approach for seismoelectric measurements using multielectrode arrangements: II—laboratory measurements, *Geophys. J. Int.*, **214**(3), 1783–1799.

* Dietrich, M., Devi, M.S., Garambois, S., Brito, D. & Bordes, C., 2018. A novel approach for seismoelectric measurements using multielectrode arrangements - I: theory and numerical experiments, *Geophys. J. Int.*, **215**(1), 61–80.

* Dupuis, C., 2008. Field measurements and analyses of electrokinetic seismoelectric signals generated in sedimentary environments, *PhD thesis*, University of New Brunswick.

* Dupuis, J.C., Butler, K.E., Kepic, A.W. & Harris, B.D., 2009. Anatomy of a seismoelectric conversion: measurements and conceptual modeling in boreholes penetrating a sandy aquifer, *J. geophys. Res.*, **114**(B10), doi:10.1029/2008JB005939.

* Frenkel, J., 1944. On the theory of seismic and seismoelectric phenomena in a moist soil, *Izvestiya Akademii Nauk SSSR*, **8**(4), 134–157.

* Garambois, S., 1999. Études expérimentales et théoriques des conversions d'ondes sismo-électriques dans les milieux poreux superficiels, *PhD thesis*, L'Université Joseph Fourier - Grenoble I.

* Garambois, S. & Dietrich, M., 2001. Seismoelectric wave conversions in porous media: field measurements and transfer function analysis, *Geophysics*, **66**(5), 1417–1430.
Garambois, S. & Dietrich, M., 2002. Full waveform numerical simulations of seismoelectromagnetic wave conversions in fluid-saturated stratified porous media, *J. geophys. Res.*, **107**(B7), ESE5–1-ESE5-18.

Grobbe, N. & Slob, E.C., 2016. Seismo-electromagnetic thin-bed responses: natural signal enhancements?, *J. geophys. Res.*, **121**(4), 2460–2479.

Grobbe, N., Slob, E.C. & Thorbecke, J.W., 2016. Comparison of eigenvectors for coupled seismo-electromagnetic layered-Earth modelling, *Geophys. J. Int.*, **206**(1), 152–190.

Guichet, X., Jouniaux, L. & Pozzi, J.-P., 2003. Streaming potential of a sand column in partial saturation conditions, *J. geophys. Res.*, **108**(B3), doi:10.1029/2001JB001517.

Guichet, X., Jouniaux, L. & Catel, N., 2006. Modification of streaming potential by precipitation of calcite in a sandwater system: laboratory measurements in the pH range from 4 to 12, *Geophys. J. Int.*, **166**(1), 445–460.

Haartsen, M.W., 1995. Coupled electromagnetic and acoustic wavefield modeling in poro-elastic media and its applications in geophysical exploration, *PhD thesis*, Massachusetts Institute of Technology.

Haartsen, M.W. & Pride, R.S., 1994. Modeling of coupled electroseismic wave propagation from point sources in layered media, in *SEG Technical Program Expanded Abstracts 1994*, pp. 1155–1158, Society of Exploration Geophysicists.

Haartsen, M.W. & Pride, S.R., 1997. Electroseismic waves from point sources in layered media, *J. geophys. Res.*, **102**(B11), 24745–24769.

Haines, S.S. & Pride, S.R., 2006. Seismoelectric numerical modeling on a grid, *Geophysics*, **71**(6), doi:10.1190/1.2357789.

Hasanian, M. & Lissenden, C.J., 2017. Assessment of coating layers on the accuracy of displacement measurement in laser Doppler vibrometry, *AIP Conf. Proc.*, **1806**(1), doi:10.1063/1.4974600.

Holzhauer, J., 2015. Seismic and seismoelectric monitoring of an unconsolidated porous medium, *PhD thesis*, Université de Pau et des Pays de l'Adour.

Holzhauer, J., Brito, D., Bordes, C., Brun, Y. & Guatarbes, B., 2017. Experimental quantification of the seismoelectric transfer function and its dependence on conductivity and saturation in loose sand, *Geophys. Prospect.*, **65**(4), 1097–1120.

Hudson, J.A., 1969. A quantitative evaluation of seismic signals at teleseismic distances—I. Radiation from point sources, *Geophys. J. R. astr. Soc.*, **18**(3), 233–249.

Ivanov, A., 1939. Effect of electrization of earth layers by elastic waves passing through them, *Proceedings of the USSR Academy of Sciences (Dokl. Akad. Nauk SSSR)*, **24**(1), 42–45.

Ivanov, A., 1940. Seismoelectric effect of second kind, *Izvd. Akad. Nauk SSSR, Ser. Geogr. i Geofiz.*, **5**, 699–726.

Jardani, A. & Revil, A., 2015. Seismoelectric couplings in a poroelastic material containing two immiscible fluid phases, *Geophys. J. Int.*, **202**(2), 850–870.

Jougnot, D. & Solazzi, S.G., 2021. Predicting the frequency-dependent effective excess charge density: a new upscaling approach for seismoelectric modeling, *Geophysics*, **86**(5), WB19–WB28.

Jougnot, D., Linde, N., Revil, A. & Doussan, C., 2012. Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils, *Vadose Zone J.*, **11**(1), doi:10.2136/vzj2011.0086.

Jougnot, D., Rubino, J., Rosas-Carbajal, M., Linde, N. & Holliger, K., 2013. Seismoelectric effects due to mesoscopic heterogeneities, *Geophys. Res. Lett.*, **40**, 2033–2037.

Jouniaux, L. & Bordes, C., 2012. Frequency-dependent streaming potentials: a review, *Int. J. Geophys.*, **2012**, 1–11.

Kennett, B., 2009. *Seismic Wave Propagation in Stratified Media, DOAB Directory of Open Access Books*, ANU E Press.

Kennett, B. L.N. & Kerry, N.J., 1979. Seismic waves in a stratified half space, *Geophys. J. Int.*, **57**(3), 557–583.

Liu, Y., Smeulders, D., Su, Y. & Tang, X., 2018. Seismoelectric interface electromagnetic wave characteristics for the finite offset vertical seismoelectric profiling configuration: theoretical modeling and experiment verification, *J. acoust. Soc. Am.*, **143**(1), EL13–EL18.

Liu, Z., Yuan, L., Zhang, X., Liu, Z. & Wu, H., 2008. A laboratory seismo-electric measurement for the permafrost model with a frozen-unfrozen interface, *Geophys. Res. Lett.*, **35**(21), 1–4.

Lorne, B., Perrier, F. & Avouac, J.-P., 1999. Streaming potential measurements: 2. Relationship between electrical and hydraulic flow patterns from rock samples during deformation, *J. geophys. Res.*, **104**(B8), 17879–17896.

Martins-Gomes, V., 2022. Experimental characterization and modelling of electromagnetic waves generated by seismoelectric conversion at porous media interfaces, *PhD thesis*, Université de Pau et des Pays de l'Adour.

Meyer, R.-C., 2021. Mathematical modeling and simulation of waves in conducting poroelastic media using HDG method, *PhD thesis*, Université de Pau et des Pays de l'Adour.

Morency, C. et al., 2021. Seismoelectric effects for geothermal resources assessment and monitoring, in *IMAGE 2021 - Society of Exploration Geophysicists Annual Meeting and Exposition*, Denver, United States.

Neev, J. & Yeatts, F.R., 1989. Electrokinetic effects in fluid-saturated poroelastic media, *Phys. Rev. B*, **40**(13), 9135–9141.

Packard, R.G., 1953. Streaming potentials across glass capillaries for sinusoidal pressure, *J. Chem. Phys.*, **21**(2), 303–307.

Parkhomenko, E.I., 1971. *Electrification Phenomena in Rocks*, Springer US.

Peng, R., Di, B., Wei, J., Hu, T. & Liu, Z., 2016. Seismoelectric effects: seismoelectric response of single interface, in *SPG/SEG 2016 International Geophysical Conference*, Beijing, China, 20-22 April 2016, pp. 658–661, Society of Exploration Geophysicists and Society of Petroleum Geophysicists.

Peng, R., Di, B., Wei, J., Ding, P., Zhao, J., Pan, X. & Liu, Z., 2017. Experimental study of the seismoelectric interface response in wedge and cavity models, *Geophys. J. Int.*, **210**(3), 1703–1720.

Peng, R., Di, B., Wei, J., Glover, P., Lorinczi, P., Ding, P. & Liu, Z., 2018. The seismoelectric coupling in shale, in *Proceedings of the 80th EAGE Conference and Exhibition 2018*, June 2018, Vol. 2018, pp. 1–5, European Association of Geoscientists & Engineers.

Peng, R., Di, B., Glover, P., Wei, J., Lorinczi, P., Ding, P., Liu, Z., Zhang, Y. & Wu, M., 2019. The effect of rock permeability and porosity on seismoelectric conversion: experiment and analytical modelling, *Geophys. J. Int.*, **219**(1), 328–345.

Pride, S., 1994. Governing equations for the coupled electromagnetics and acoustics of porous media, *Phys. Rev. B*, **50**(21), 15678–15696.

Pride, S.R., 2005. Relationships between seismic and hydrological properties, in *Hydrogeophysics*, pp. 253–290, eds Rubin, Y. & Hubbard, S.S., Springer Netherlands.

Pride, S.R. & Garambois, S., 2002. The role of Biot slow waves in electroseismic wave phenomena, *J. acoust. Soc. Am.*, **111**(2), 697–706.

Pride, S.R. & Garambois, S., 2005. Electroseismic wave theory of Frenkel and more recent developments, *J. Eng. Mech.*, **131**(9), 898–907.

Pride, S.R. & Haartsen, M.W., 1996. Electroseismic wave properties, *J. acoust. Soc. Am.*, **100**(3), 1301–1315.

Ren, H., Huang, Q. & Chen, X., 2010. A new numerical technique for simulating the coupled seismic and electromagnetic waves in layered porous media, *Earthq. Sci.*, **23**(2), 167–176.

Reppert, P.M., Morgan, F.D., Lesmes, D.P. & Jouniaux, L., 2001. Frequency-dependent streaming potentials, *J. Colloid Interf. Sci.*, **234**(1), 194–203.

Revil, A. & Glover, P. W.J., 1998. Nature of surface electrical conductivity in natural sands, sandstones, and clays, *Geophys. Res. Lett.*, **25**(5), 691–694.

Revil, A. & Mahardika, H., 2013. Coupled hydromechanical and electromagnetic disturbances in unsaturated porous materials, *Water Resour. Res.*, **49**(2), 744–766.

Rosas-Carbajal, M., Jougnot, D., Rubino, J.G., Monachesi, L., Linde, N. & Holliger, K., 2020. *Seismoelectric Signals Produced by Mesoscopic Heterogeneities*, Chapter 19, pp. 269–287, American Geophysical Union.

Schakel, M.D., Smeulders, D. M.J., Slob, E.C. & Heller, H. K.J., 2011a. Seismoelectric interface response: experimental results and forward model, *Geophysics*, **76**(4), N29–N36.

Schakel, M.D., Smeulders, D. M.J., Slob, E.C. & Heller, H.K.J., 2011b. Laboratory measurements and theoretical modeling of seismoelectric interface response and coseismic wave fields, *J. appl. Phys.*, **109**(7), doi:10.1063/1.3567945.
Schakel, M.D., Smeulders, D. M.J., Slob, E.C. & Heller, H. K.J., 2012. Seismoelectric fluid/porous-medium interface response model and measurements, Transport Porous Media, 93(2), 271–282.

Shen, C., 2020. Experimental and numerical studies of seismic wave propagation in carbonate rocks at the laboratory scale, PhD thesis, Université de Pau et des Pays de l'Adour.

Shen, C., Brito, D., Diaz, J., Sanjuan, F., Bordes, C. & Garambois, S., 2022. Pulsed-laser source characterization in laboratory seismic experiments, Geomech. Geophys. Geo-Ener. Geo-Resour., 8(1), 1–26.

Shi, P., Guan, W. & Hu, H., 2018. Dependence of dynamic electrokinetic-coupling-coefficient on the electric double layer thickness of fluid-filled porous formations, Ann. Geophys., 61(3), 1–12.

Smeulders, D., Grobbe, N., Heller, H. & Schakel, M., 2014. Seismoelectric conversion for the detection of porous medium interfaces between wetting and nonwetting fluids, seismoelectric conversion for the detection of porous medium interfaces between wetting and nonwetting fluids, Vadose Zone J., 13(5), 1–7.

Thanh, L.D., Jougnot, D., Solazzi, S.G., Van Nghia, N. & Van Do, P., 2021. Dynamic streaming potential coupling coefficient in porous media with different pore size distributions, Geophys. J. Int., 229(1), 720–735.

Thompson, R.R., 1936. The seismic electric effect, Geophysics, 1(3), 327–335.

Walton, K., 1987. The effective elastic moduli of a random packing of spheres, J. Mech. Phys. Solids, 35(2), 213–226.

Zhu, Z. & Nafi Toksöz, M., 2013. Experimental measurements of the stream- ing potential and seismoelectric conversion in Berea sandstone, Geophys. Prospect., 61(3), 688–700.

Zhu, Z. & Toksöz, M.N., 1996. Experimental studies of electrokinetic con- versions in fluid-saturated porous medium, in SEG Technical Program Expanded Abstracts 1996, pp. 1699–1702, Society of Exploration Geo- physicists.

Zhu, Z. & Toksöz, M.N., 2003a. Crosshole seismoelectric measurements in borehole models with fractures, Geophysics, 68(5), 1519–1524.

Zhu, Z. & Toksöz, M.N., 2003b. Effects of Saturant Conductivity on Seis- moelectric Conversion, Tech. rep., Massachusetts Institute of Technology, Earth Resources Laboratory.

Zhu, Z. & Toksöz, M.N., 2005. Seismoelectric and seismomag- netic measurements in fractured borehole models, Geophysics, 70(4), F45–F51.

Zhu, Z., Cheng, C.H. & Toksöz, M.N., 1994. Electrokinetic conversion in a fluid-saturated porous rock sample, in SEG Technical Program Expanded Abstracts 1994, pp. 1057–1060, Society of Exploration Geophysicists.

Zhu, Z., Haartsen, M.W. & Toksöz, M.N., 1999. Experimental stud- ies of electrokinetic conversions in fluid-saturated borehole models, Geophysics, 64(5), 1349–1356.

Zhu, Z., Haartsen, M.W. & Toksöz, M.N., 2000. Experimental studies of seismoelectric conversions in fluid-saturated porous media, J. geophys. Res., 105(B12), 28 055–28 064.

Zhu, Z., Toksöz, M.N. & Burns, D.R., 2008. Electroseismic and seismoelec- tric measurements of rock samples in a water tank, Geophysics, 73(5), E153–E164.

Zhu, Z., Toksöz, M. & Zhan, X., 2016. Seismoelectric measurements in a porous quartz-sand sample with anisotropic permeability, Geophys. Prospect., 64(3), 700–713.

Zyserman, F.I., Gauzellino, P.M. & Santos, J.E., 2010. Finite element mod- eling of SHTE and PSVTM electroseismics, J. appl. Geophys., 72(2), 79–91.

APPENDIX A: COEFFICIENTS CONTROLLING THE COUPLING BETWEEN MECHANICAL AND ELECTROMAGNETIC DISTURBANCES

In Pride's theory (Pride 1994), all frequency-dependant coupling between Biot's and Maxwell's equations is controlled by the

frequency-dependant coupling coefficient $L(\omega)$:

$$L(\omega) = L_0 \left[ 1 - i \frac{\omega}{\omega_t} \frac{m}{4} \left( 1 - 2 \frac{\tilde{d}}{\Lambda} \right)^2 \left( 1 - i^{3/2} \tilde{d} \sqrt{\frac{\omega \rho_f}{\eta}} \right)^2 \right]^{-1/2}, \quad (A1)$$

with $\tilde{d}$ taken here as the Debye length, according to eqs (29) and (210) of Pride (1994), which is a measure of the diffuse double layer thickness. $\Lambda$ is usually defined as the pore radius, and $\omega_t$ is the transitional angular frequency and refers to the transition from viscous (low-frequency) to inertial (high-frequency) flow:

$$\omega_t = \frac{\phi}{\alpha_{\infty} k_0} \frac{\eta}{\rho_f}. \quad (A2)$$

At last, $L_0$ is the low-frequency coupling coefficient defined as:

$$L_0 = - \frac{\phi \epsilon_0 k_f \zeta}{\alpha_{\infty} \eta} \left( 1 - 2 \frac{\tilde{d}}{\Lambda} \right) \Leftrightarrow - \frac{\phi}{\alpha_{\infty}} C_{ek} \sigma_f \left( 1 - 2 \frac{\tilde{d}}{\Lambda} \right), \quad (A3)$$

where $C_{ek}$ refers to the electrokinetic coupling obtained from streaming potential measurements.

In order to be consistent with the medium used in the experiments (Section 3.1.1) the following relation for $C_{ek}$, that considers surface conductivity ($\sigma_n$), will be used in eq. (A3) (Guichet et al. 2003, 2006):

$$C_{ek} = \frac{\epsilon_0 k_f \zeta}{\eta \sigma_f} \frac{1}{1 + m_c \left( \frac{\alpha_{\infty}}{\phi} - 1 \right)^{\frac{\sigma_n}{\sigma_f}}}. \quad (A4)$$

Parameter $m_c$ is the cementation exponent (see Archie 1942) which here is defined as $m_c = (\log(\phi) - \log(\alpha_{\infty})) / \log(\phi)$, using the defi- nition of formation factor.

To illustrate how $L(\omega)$ behaves for both sand and sandstone, Fig. A1 displays the variation of $|L(\omega)|$ as a function of frequency. To calculate both curves, we have used the values in Table 1.

APPENDIX B: DERIVING $E_z$

In what follows, we present the details on how to obtain $E_z$, the z-component of the electric field. The derivations presented below are necessary because the original computer code of Garambois & Dietrich (2002) was restricted to the computation of $u_x$, $u_y$, $u_z$, $E_x$, $E_y$, and $H_y$. We start with the governing equations of Pride (1994), classified according to the physical phenomena they represent.

Mechanical equations (Biot) read:

$$\nabla \cdot \tau = -\omega^2 (\rho u + \rho_f w) + \mathbf{F}, \quad (B1a)$$

$$\tau = (K_G \nabla \cdot \mathbf{u} + C \nabla \cdot \mathbf{w})\mathbf{I} + G_{fr} (\nabla \mathbf{u} + \nabla \mathbf{u}^T - \frac{2}{3} \nabla \cdot \mathbf{u}\mathbf{I}), \quad (B1b)$$

$$-P = C \nabla \cdot \mathbf{u} + M \nabla \cdot \mathbf{w}; \quad (B1c)$$

Maxwell's equations write:

$$\nabla \times \mathbf{E} = i\omega\mathbf{B} - \mathbf{M}, \quad (B1d)$$

$$\nabla \times \mathbf{H} = -i\omega\mathbf{D} + \mathbf{J} + \mathbf{C}, \quad (B1e)$$

$$\mathbf{D} = \epsilon \mathbf{E}, \quad (B1f)$$

$$\mathbf{B} = \mu \mathbf{H}; \quad (B1g)$$

Transport equations are:

$$\mathbf{J} = \sigma(\omega)\mathbf{E} + L(\omega)(-\nabla P + \omega^2 \rho_f \mathbf{u} + \mathbf{f}), \quad (B1h)$$
<figure><img src="image_16.png" /><figcaption>**Figure A1.** Variation of the coupling coefficient $L(\omega)$ as a function of frequency computed for sand and sandstone, the dashed blue vertical line indicates the transitional frequency in Hertz.</figcaption></figure>

$$-i\omega\mathbf{w} = L(\omega)\mathbf{E} + \frac{k(\omega)}{\eta}(-\nabla P + \omega^2\rho_f\mathbf{u} + \mathbf{f}). \quad (B1i)$$

In the equations above, $\tau$ is the stress tensor (acting on both fluid and solid phases), $\mathbf{u}$ the solid particle displacement, $\mathbf{w}$ the relative fluid-to-solid displacement ($\mathbf{w} = \phi(\mathbf{u}_{\text{fluid}} - \mathbf{u})$), $\mathbf{I}$ the identity matrix, $\mathbf{E}$ the electric field, $\mathbf{J}$ the current density, $\mathbf{D}$ the dielectric displacement field, $\mathbf{B}$ the magnetic-flux density, $\mathbf{H}$ the magnetic field, $\mathbf{P}$ the pressure, and $\mathbf{M}, \mathbf{C}, \mathbf{f}$ and $\mathbf{F}$ are source terms corresponding, respectively, to: magnetic-current, current-density, body force density acting on fluid and body force density acting on both phases. $\sigma(\omega)$, $k(\omega)$ and $L(\omega)$ are the frequency-dependent coefficients of the transport equations derived in Pride (1994) and correspond, respectively, to: the total conductivity of the porous medium, the dynamic permeability, and the coupling coefficient which controls all seismic-to-electromagnetic (and vice versa) energy conversion. $\epsilon$ is the electric permittivity of the porous medium, as defined in Pride (1994), and $\mu$ the magnetic permeability, which for all purposes is considered equal to that of the vacuum $\mu_0$.

To deduce $E_z$ from the governing equations above, the relations proposed in Hudson (1969) are employed to simplify the cylindrical representation $(\rho, \theta, z)$ of the wavefields' horizontal components. For the present case, the relations used are the following:

$$E_V = \rho^{-1}(\partial_\rho(\rho E_\rho) + \partial_\theta E_\theta), \quad (B2a)$$

$$H_H = \rho^{-1}(\partial_\rho(\rho H_\theta) - \partial_\theta H_\rho), \quad (B2b)$$

$$M_H = \rho^{-1}(\partial_\rho(\rho M_\theta) - \partial_\theta M_\rho), \quad (B2c)$$

$$\nabla_1^2 \varphi = \rho^{-1}(\partial_\rho(\rho \varphi_\rho) + \partial_\theta \varphi_\theta); \quad (B2d)$$

where $\varphi$ can represent for any vector field.

The first step to find $E_z$ consists in isolating the z-components of the cylindrical representation of eq. (B1e):

$$\begin{gathered} (\rho^{-1}\partial_{\theta}H_z - \partial_z H_{\theta})\rho + (\partial_z H_{\rho} - \partial_{\rho}H_{\rho})\theta + \rho^{-1}(\partial_{\rho}\rho H_{\theta} - \partial_{\theta}\rho H_{\rho}) \\ z = -i\omega\epsilon\mathbf{E} + \mathbf{J} + \mathbf{C} \\ \therefore \rho^{-1}(\partial_{\rho}\rho H_{\theta} - \partial_{\theta}H_{\rho}) = -i\omega\epsilon E_z + J_z + C_z. \end{gathered} \quad (B3)$$

Furthermore, the algebraic manipulation of the transport eqs (B1h and B1i) leads to:

$$J_z = E_z(\sigma - L^2 \frac{\eta}{k(\omega)}) - i\omega w_z \frac{\eta}{k(\omega)} L, \quad (B4)$$

which can be substituted in eqs (B3) together with (B2b) to find $E_z$:

$$\begin{aligned} H_H &= -i\omega\epsilon E_z(\epsilon + i\frac{\sigma}{\omega} - iL^2\frac{\eta}{\omega k(\omega)}) - i\omega L \frac{\eta}{k(\omega)} w_z + C_z \quad (B5) \\ &= -i\omega\tilde{\epsilon} E_z - \omega^2 \tilde{\rho} L w_z + C_z. \end{aligned}$$

Now, to be consistent with what is computed in the code, a discrete Fourier-Hankel (HF) transform is applied (see Haartsen & Pride 1997; Kennett 2009), and the definition $\hat{H}_2 = \hat{H}_H/ik$, where k represents the horizontal wavenumber, is used to ensure all wavefields have the same dimensionality. At last, eq. (B5) becomes:

$$\hat{E}_z = -\frac{k}{\omega\tilde{\epsilon}}\hat{H}_2 + \frac{i\omega\tilde{\rho}L}{\tilde{\epsilon}}\hat{w}_z - \frac{i}{\omega\tilde{\epsilon}}\hat{C}_z. \quad (B6)$$

Finally, because in this work only seismoelectric phenomena are considered, both C and M can be disregarded.