# Experimental Validation of the Electrokinetic Theory and Development of Seismoelectric Interferometry by Cross-Correlation

Authors: F. C. Schoemaker; N. Grobbe; M. D. Schakel; S. A. L. de Ridder; E. C. Slob; D. M. J. Smeulders
Journal: International Journal of Geophysics
 (2012)

---

## Abstract

We experimentally validate a relatively recent electrokinetic formulation of the streaming potential (SP) coefficient as developed by Pride (1994). The start of our investigation focuses on the streaming potential coefficient, which gives rise to the coupling of mechanical and electromagnetic fields. It is found that the theoretical amplitude values of this dynamic SP coefficient are in good agreement with the normalized experimental results over a wide frequency range, assuming no frequency dependence of the bulk conductivity. By adopting the full set of electrokinetic equations, a full-waveform wave propagation model is formulated. We compare the model predictions, neglecting the interface response and modeling only the coseismic fields, with laboratory measurements of a seismic wave of frequency 500 kHz that generates electromagnetic signals. Agreement is observed between measurement and electrokinetic theory regarding the coseismic electric field. The governing equations are subsequently adopted to study the applicability of seismoelectric interferometry. It is shown that seismic sources at a single boundary location are sufficient to retrieve the 1D seismoelectric responses, both for the coseismic and interface components, in a layered model.


---

## Body

## 1. Introduction

The first observation of coupling between electromagnetic and mechanical effects (also known as electroosmosis, which is one of the electrokinetic effects) dates back to the beginning of the 19th century. In 1809, Reuss [1] was the first to report on an experiment where a direct current was applied to a clay-sand-water mixture. The experiment was performed with a U-tube, filled with quartz at the bottom. Application of an electric current caused the water to rise in the leg containing the negative electrode [2].The electrokinetic effect works as follows. In a fully fluid-saturated porous medium, a charged nanolayer at the solid-liquid interface is present (see Figure1). The origin of this charged nanolayer lies in the presence of an aqueous solution, typically a negatively charged silane grain surface. The resulting interface potential is called the zeta-potential, which is typically negative and on the order of a few tens of millivolts [9]. The counterions in the fluid reorganize in a layer that is bound to the surface by electrostatic forces (Stern layer) and a diffuse layer that is free to flow. In the diffuse layer two types of physical phenomena are competing, the electrostatic forces between the ions and the Brownian motion of the particles. This effectively results in an exponentially decreasing potential away from the solid-liquid interface towards the bulk of the pore (see Figure 1). The characteristic length over which the EDL exponentially decays, known as the Debye length, is of the order of a few tens of nanometers for typical reservoir rocks. The Stern layer and the diffuse layer together are usually called the electric double-layer (EDL), see Figure 1. The Debye length is considerably thinner than any viscous boundary layer that normally develops in pore fluid transport phenomena [3]. Quincke [10] performed electroosmosis experiments on glass capillaries. The simple geometry used, allowed for controlled experimental conditions. Linearity between the electroosmotic volume and the applied electric field was observed. Another electrokinetic effect, the physical phenomenon of electrophoretic mobility, where particles are mobilized by electrical fields, was described by Quincke [10] together with Reuss [1]. A mathematical description of both phenomena (electroosmotic and electrophoretic mobility) was later derived by Helmhlotz [11]. However he did not consider the electric permittivity. Von Smoluchowski [12] derived the well-known Helmholtz-Smoluchowski equation, in which the electric permittivity is incorporated. Smoluchowski also recognized reciprocity between electroosmotic flow and streaming potential phenomena (mechanical to electromagnetic effect), later described by Onsager [13, 14].Figure 1
Electric double-layer according to the Stern model. The inner and outer Helmholtz planes are indicated as IH and OH, respectively. The slipping plane is denoted by S and its charge is characterized by theζ-potential (modified from [2–8]).Gouy [16] and Chapman [17] improved the theoretical model by including the diffuse layer of counterions in the model, thereby relating the thickness of the diffuse layer to the ionic strength of the solution [3]. To overcome limitations with highly charged electric double layers, Stern [18] added another layer to the model, see Figure 1. This theorem was some years later perfected by Derjaguin and Landau [19], and also by Verwey and Overbeek [4] in the “DLVO” theory, which describes in even more detail the forces between charged surfaces interacting through an electrolyte.In 1936, Thompson [20] suggested that the electrokinetic effect could be used for geophysical prospecting. The Russian physicist Yacov Il’ich Frenkel [21] developed a theory for wave propagation of electrokinetic phenomena in fluid-saturated porous media, in which he predicted the slow compressional wave and the seismoelectric effect (thereby he made a marginal error in the development of the Biot-Gassmann constants, he also only considered the electric effect and not the full Maxwell equations [22]). In 1959, Martner and Sparks [23] were the first to report that an electric potential difference generated in the subsurface by the passage of seismic waves could be detected by electrode pairs. Somewhat later, an experimental programme was undertaken to evaluate the electroseismic effect as a possible means for detecting underground nuclear tests. The goal was to develop long range systems for detection of nuclear blasts, see for example, Broding et al. [24] and Long and Rivers [25]. Due to insensitive technical equipment, lack of computing power, and the success of conventional seismic and electromagnetic methods, electrokinetics never gained much attention in geophysical exploration, at least until the 1990s. Moreover, the majority of field tests up to that time were concerned with the seismoelectric effect while the reciprocal electroseismic effect was underexposed. Extended field tests were only performed recently [26].Regarding wave modeling, Neev and Yeatts [28] were the first in recent history (since Frenkel) to postulate a set of equations, which attempted to model the interaction between mechanical waves and electric fields due to electrokinetics. Their model did not include the Maxwell equations and frequency-dependence of the transport laws. A possible way to include all effects is by volume averaging the continuum equations for solid grains and electrolyte fluids. Using this approach, Pride [29] obtained the governing equations for coupled electromagnetics and elastodynamics of porous media.The governing equations of Pride describe coupled seismic and electromagnetic wave propagation effects. A schematic description of the coseismic and interface response effects is given in Figure2. Figures 2(a) and 2(c) show a cross-section of the subsurface, with corresponding seismogram and seismoelectrograms in Figures 2(b) and 2(d), respectively. The subsurface consists of two layers. Geophones and electrodes are positioned at the surface. In Figure 2(a), a local pressure disturbance is initiated at t=0. Due to the mechanical pressure source, a longitudinal wave is created (labelled 1 in Figures 2(a)–2(d)). The seismic wave creates a fluid pressure gradient within the pulse that induces pore fluid flow. Excess electrical charge in the double layer is transported by this flow. The net flow of charge relative to the grains is known as the streaming electric current. The induced conduction current leads to the electric field known as the “coseismic field” [5, 6, 30]. The coseismic field travels along with the seismic wave, giving it the same velocity as the compressional wave (compare the slope of event 1 in the left and right part of Figure 2(b)). When the pressure wave encounters the interface (with changing medium parameters) between porous layers 1 and 2, this results in a local asymmetry in the charge distribution. This will induce an oscillating electric dipole (Figure 2(c)). The associated independent electromagnetic field will travel almost immediately to the receiver electrodes (Figure 2(d), right part). This seismoelectric effect is known as the “interface response field”. The coseismic and interface response fields were measured in the laboratory (e.g., [7, 31–34]) and in the field (e.g., [25, 35–40]). Zhu and Toksöz [41] and Bordes et al. [34, 42] reported on coseismic magnetic field measurements associated with a Stoneley wave and a shear wave, respectively. The dynamic SP coefficient, that links the mechanical and electromagnetic fields in Pride’s set of equations, was measured by Reppert et al. [8] and another validation is presented in this paper. Also, full-waveform seismoelectric models that adopt Pride’s theory were compared with measurements. Mikhailov et al. [36] and Haines et al. [39] compare seismoelectric synthetic sections with field measurements and find qualitative agreement. Zhu et al. [32] found kinematic agreement between full-waveform seismoelectric predictions and laboratory measurements. Block and Harris [7] compared amplitudes of coseismic wave fields in sands with numerically predicted amplitudes and fitted their measurements to Pride’s theory by incorporating an additional surface conductivity term. Charara et al. [43] found agreement between modeled and measured seismoelectric waveforms and amplitudes at a fluid/porous-medium interface in a laboratory setup. Schakel et al. [44, 45] found agreement between laboratory measurements of the coseismic and interface response fields and full-waveform and spatial seismoelectric predictions in terms of traveltime, waveform, and spatial amplitude pattern. Seismic waves can image to great depths but at the cost of resolution. Electromagnetic waves are sensitive to additional material properties and can therefore provide us with information about the pore fluid content. Seismoelectric conversion methods in field studies can combine seismic resolution and electromagnetic hydrocarbon sensitivity [39].Schematic of a “standard” geometry for a seismoelectric survey (modified from Haines [5], Kroeger [6, 15]). One fluid-saturated porous layer (layer 1) overlies a porous half-space (layer 2), with contrasting medium parameters (panels (a) and (c)). The seismic source generates longitudinal wave (labelled 1) that is recorded by surface geophones (left part of panel (b)). This wave also creates a coseismic electric field that is registered by the surface electrodes (right part of panel (b)). When the longitudinal wave arrives at the interface between the porous layers, the interface field is generated (labelled 2), that is also recorded by the surface electrodes (right part of panel (d)).
(a)(b)(c)(d)However, in seismoelectric surveys, the interface response is known to be very weak, that is, the response suffers from a very low signal-to-noise ratio. Therefore, the sources in classical seismoelectric surveys need to be strong. This is not always possible and therefore it is beneficial to be able to replace those strong sources by receivers: the principle of interferometry. In addition, by doing interferometry, stacking inherently takes place with a possible improvement of the signal-to-noise ratio as a result [46]. From an imaging point of view, the principle of interferometry has already been proven useful for a wide class of phenomena, for example in seismic and electromagnetic systems (e.g., Wapenaar et al. [47], Slob et al. [48]). Seismic interferometry is a seismological technique which makes use of the cross-correlation of responses at different receivers in order to obtain the Green’s function between these receivers [49]. It can include both passive and active sources. Due to the fact that the cross-correlation generates new data from measured data, it may allow for improved imaging compared to the situations where imaging algorithms are applied to the measured data only.The foundations of the principle of interferometry were lain in 1968 by Claerbout who showed that by using the autocorrelation of the 1D transmission response of a horizontally layered medium (bounded by a free surface), the reflection response of this medium can be obtained [50, 51]. Later, Claerbout conjectured that this relation could also be extended for 3D inhomogeneous situations, which was proven by Wapenaar [52]. By cross-correlating the recorded noise at two locations on the surface, it is possible to construct the wavefield that would be recorded at one of the locations as if there was a source at the other [53]. For a detailed overview of the theory of interferometry (e.g., stationary phase arguments, controlled-sources, interferometric imaging), the reader is referred to Wapenaar et al. [49, 54, 55] and Schuster [56]. Wapenaar et al. [57] showed the link between the principle of reciprocity and seismic interferometry. Using the reciprocity theorem of the correlation type, they generalized Claerbout’s relation between transmission and reflection responses to 3D inhomogeneous acoustic and elastic media. This theory was confirmed with numerically modeled seismic data in laterally varying media [58]. Wapenaar et al. [47] have shown that using cross-correlation to retrieve the Green’s function response between two stations is in principle not limited to seismic systems but holds for a wide class of phenomena, including seismoelectromagnetic wave propagation. We take the principle to the next level by numerically simulating seismoelectric interferometry by cross-correlation. de Ridder et al. [46] have already shown, with three numerical examples, that it is indeed possible (under certain conditions) to obtain accurate Green’s functions from boundary sources only. Here, we will increase the complexity of the numerical configuration by adding an extra layer to the system, to investigate the Green’s function retrieval for a 1D, three-layered system bounded by a free-surface.Although the individual constituents of Pride’s model (i.e., Biot’s theory and Maxwell’s theory) have been experimentally validated, the dynamicSP coefficient that links these theories has been rarely studied (for a review see Jouniaux and Bordes, this issue). Also, despite the experimental verification of the coseismic and interface response fields, direct comparisons between electrokinetic wave propagation theory and measurements are scarce. In this paper we (1) validate electrokinetic theory by measurements and (2) investigate the applicability of correlation imaging with coupled seismic and electromagnetic wave propagation. First we present Pride’s electrokinetic governing equations. Second, the theoretical dynamic SP coefficient is compared against normalized measurements. Third, a seismoelectric wave model is formulated and model predictions are compared against seismoelectric wave propagation measurements.It is shown that measurements of both the dynamicSP coefficient and the coseismic wave field are adequately described by the electrokinetic theory. This theory is subsequently adopted, when we numerically investigate the applicability of correlation imaging with seismoelectromagnetic waves.


## 2. Governing Equations

The governing equations for seismoelectric and electroseismic wave propagation in a fluid-saturated porous medium are derived from the compilation of Biot’s theory [60, 61] together with Maxwell’s theory. The Biot equations describe the acoustic side of electrokinetic phenomena. They are a combination of momentum equations and the stress strain relationships for an isotropic material, together with the continuity equations [62–64].Expressing the expanded Biot equations, for the solid as well as the fluid and adopting aneiωt time dependence, yields the following linearized set of governing equations(1)G∇2ûs+(A+G)∇(∇⋅ûs)+Q∇(∇⋅ûf)=-ω2[ρ11(ω)ûs+ρ12(ω)ûf]+ηϕL(ω)k(ω)Ê,(2)Q∇(∇⋅ûs)+R∇(∇⋅ûf)=-ω2[ρ12(ω)ûs+ρ22(ω)ûf]-ηϕL(ω)k(ω)Ê,
where A, Q, R are the Biot Gassmann constants [62], G the shear modulus, ρf is the fluid density, ρs is the solid density, α∞ is the tortuosity, η is the fluid viscosity, k0 is the (static) permeability, Ê is the electric field, and L(ω) is the dynamic electrokinetic coupling [29](3)L(ω)L0=[1+iωωc2m(1-2dΛ)2(1+diωρfη)2]-1/2,
where Λ is a characteristic pore size parameter and m is the shape factor. Please note that (3) is written in a slightly different form than in [29], because we used Johnson’s definition of the shape factor [65]: m=8α∞k0/(ϕΛ2). The characteristic (or rollover) frequency ωc is defined as ωc=ηϕ/(α∞k0ρf) [59]. The Debye length is denoted by d (see, e.g., [29]) and L0 represents the static electrokinetic coupling for a porous medium(4)L0=-ϕα∞ɛ0ϵrfζη(1-2dΛ),
where ɛ0 is the vacuum permittivity, ϵrf is the pore fluid relative permittivity, and ζ is the zeta-potential. We note that Pride [29] uses an additional relaxation mechanism when the complex viscous skin depth η/(ωρf) becomes smaller than the Debye length. However, due to the fact that the Debye length d for most salinity cases [22] is much smaller than Λ, Pride’s relaxation mechanism can often be neglected. The dynamic permeability is closely related to the viscous correction factor(5)k(ω)k0=(F(ω)+iωωc)-1,
the viscous correction factor is defined by Johnson et al. [59] as(6)F(ω)=1+miω2ωc.
The coefficients ρ11(ω), ρ12(ω), and ρ22(ω) are the so-called generalized effective density functions [65](7)ρ11(ω)=(1-ϕ)ρs-ρ12(ω),ρ12(ω)=ϕρf[1+iηϕωρfk(ω)],ρ22(ω)=ϕρf-ρ12(ω).
Considering the definitions for ρ12(ω) and ρ22(ω), (2) can be written as(8)iωŵ=k(ω)η(-∇p̂+ω2ρfûs)+L(ω)Ê,
where ŵ=ϕ(ûf-ûs) is the relative displacement. Pride [29] developed the following equation coupling the streaming and the conduction currents(9)Ĵ=L(ω)(-∇p̂+ω2ρfûs)+σ(ω)Ê,
where Ĵ is the electric current density and σ(ω) the dynamic bulk conductivity. We recognize that the electrokinetic coupling is present in the mechanical and the electromagnetic equations (8) and (9) (see [29, 66]). The dynamic bulk conductivity σ(ω) for a porous medium of arbitrary pore structure is assumed to be independent of the frequency [22, 29] so that(10)σ(ω)≈ϕα∞σf,
or(11)σ0=ϕα∞σf,
where σ0 represents the bulk electric conductivity and σf the pore-fluid conductivity. Closely related to the dynamic electrokinetic coupling (3) is the dynamic SP coefficient, defined as L(ω)=C(ω)σ(ω) [8]. Using this mutual relationship together with the hypothesis of frequency independence of the dynamic bulk conductivity (11), the measured dynamic SP coefficient and dynamic coupling are mutually related in their normalized form by(12)C(ω)C0=L(ω)L0.
Eliminating (-∇p̂+ω2ρfûs) from (8) and (9), we obtain(13)iωŵ=L(ω)Ê+k(ω)ηL(ω)(J-σ(ω)Ê).
The Maxwell relation for the magnetic field is given by Ampère’s Circuit Law(14)Ĵ=∇×Ĥ-iωɛÊ,
with Ĥ the magnetic field and ɛ the electric permittivity for a fluid-saturated porous medium(15)ɛ=[ϕα∞(ϵrf-ϵrs)+ϵrs]ɛ0,
where ϵrs is the solid relative permittivity. Faraday’s induction law states that(16)iωμĤ=-∇×Ê,
with μ the magnetic permeability. Substitution of (14) in (13) results in(17)iωŵ=L(ω)Ê+k(ω)ηL(ω)(∇×H-[iωɛ+σ(ω)]Ê).
Substituting the cross-product of Faraday’s law (16) [29, 67] into (17) yields(18)μɛ¯(ω)ω2Ê+ω2ημk(ω)L(ω)ŵ=∇(∇⋅Ê)-∇2Ê,
where ɛ¯(ω) is the effective electric permittivity [67] of the porous continuum(19)ɛ¯(ω)=ɛ-iσ(ω)ω+iηL2(ω)ωk(ω).
Here -iσ(ω)/ω is a term accounting for the energy losses. The electrokinetic effect manifests itself in ɛ¯(ω) as an energy gain that is quadratic in L(ω) (third term in the right-hand side of (19)). Equations (1), (2), and (18) form a closed set of equations necessary to describe electrokinetic phenomena, for the displacements ûs, ûf (mechanical part of the equations), and electric fields Ê (the electromagnetic part).


## 3. Experimental Validation of the Dynamic Coupling Coefficients

We experimentally validatek(ω)/k0 and C(ω)/C0. The experiments are performed with the dynamic Darcy cell (DDC) as shown in Figure 3, within a steel cylinder (see [27]). At the bottom of the DDC an oscillating pressure is applied (generated by HP Agilent 33120A waveform Generator). A power amplified (Gearing and Watson) vibrating exciter (GW V20) drives a rubber membrane which induces an oscillating pressure. Vibrations are induced in a frequency band ranging from 5 Hz up to 150 Hz. Two identical piezoelectric transducers (PCB 116 Druck) are used to measure the pressure drop across the sample, one at the bottom inside the cylinder and the other mounted just above the sample in the center of the cylinder. On the top and bottom of the porous medium, electrodes are placed from which the streaming potential gradient is measured. These electrodes are sintered plates of Monel (an alloy primarily composed of Nickel and Copper). The signals from the two piezoelectric transducers are modified by means of amplifiers (Kistler 5011), and the signals of the electrodes amplified (Tektronix AM 502). The porous sample (parameters given in Table 1) consists of tubes of glass (borosilicate), which are glued together with an epoxy resin (Figure 4) and oriented in the flow direction. The combination of sintered plates together with a large surface area of the glass capillary tubes makes it possible to measure a relative strong signal. The sample is carefully saturated with degassed, demineralised water with a small amount of sodium chloride (with a density of ρf=1·103 kg/m3, a viscosity of η=0.9·10-3 Pa s obtained from [68], and a measured pore fluid conductivity of σf=1.3·10-3 S/m), whereafter the setup is left until equilibrium of the salt solution is reached.Table 1
Sample properties.
PropertySymbolValueUnitPermeabilityak02.1·10-10[m2]Shape factorbm1.75[—]Porositycϕ0.093[—]Tortuositydα∞1.8[—]Debye lengthed2.9·10-8[m]Weighted pore volume-to-surface ratiofΛ1.3·10-4[m]aThe permeability is measured directly. b, dThe shape factor and the tortuosity are derived from an independent dynamic head experiment [27, 69], by means of curve fitting. cThe porosity is computed from [70, 71]. e, fThe Debye length and the characteristic pore size are computed from theory (see [29, 59], respectively).Figure 3
Schematic of the dynamic Darcy cell with borosilicate sample and Monel disks (modified from [27]).Figure 4
Exploded view of the capillary core.The 50 Hz electromagnetic frequency radiating from the equipment is suppressed by shielding the setup and its wires (therefore use has been made of shielded twisted cable pairs). To reduce uncorrelated noise the data are averaged multiple times.In Figures5 and 6, normalized amplitude and phase values of the dynamic permeability are plotted for the theoretical solution (5) together with the laboratory measurements. At low pulsation frequencies (viscosity dominated), the normalized dynamic permeability necessarily tends to its steady-state value, whereas above the characteristic pulsation frequency (the area where viscous dominated flow switches to inertia dominated flow [65]) a strong decline can be observed. The theory correlates well with the measurement. The offsets in the lower frequency range are caused by limitations of the equipment, while in the high frequency area this difference is mainly caused by resonance of the setup.Figure 5
Amplitude of the normalized dynamic permeability. Theory of Johnson et al. [59] versus measurements, using the parameters as shown in Table 1. Experiments 1–4 denote repeated experiments at different times.Figure 6
Phase value of the dynamic permeability. Theory of Johnson et al. [59] versus measurements, using the parameters as shown in Table 1. Experiments 1–4 denote repeated experiments at different times.The measurements of the normalized dynamic SP coefficient (normalized to the measured value at 11 Hz, whereC0=1.7·10-5 V/Pa) shown in Figures 7 and 8 (using the parameters shown in Table 1), are performed by measuring the potential difference and the pressure difference across the sample between the Monel disks (see Figure 4). The rigid glass capillary tubes make it possible to assume no solid displacement us=0. Using (9) for a conservative (irrotational) electric field E=-∇U (with U the streaming potential difference) in a setup where the electric current density is equal to J=0, we obtain(20)C(ω)=-∇U(ω)∇p(ω),
with C(ω) being the dynamicSP coefficient. The dynamicSP coefficient theory agrees well with the measurement regarding the normalized amplitudes. The phase values show a large offset for the low as well as the high-frequency range. The offsets in the lower frequency range are caused by limitations of the equipment, while in the high frequency area this difference mainly is caused by resonance of the setup. This could be counteracted by applying notch filters at these higher frequencies. Due to the layering of the sample, the theory agrees well with the measurements. It is seen in measurements from [69, 72–74], that with a single capillary [8] it is possible to obtain remarkably consistent results between theory and measurement. However, the experimental setup (a set of capillaries combined with Monel disks) gives a more accurate representation of capillary networks in natural environments than a single capillary tube.Figure 7
Amplitude of the normalized dynamic SP coefficient. Theory of Pride [29] versus measurements, using the parameters as shown in Table 1  and assuming σ(ω)=σ0. Experiments 1–4 denote repeated experiments at different times.Figure 8
Phase values of the dynamic SP coefficient. Theory of Pride [29] versus measurements, using the parameters as shown in Table 1  and assuming σ(ω)=σ0. Experiments 1–4 denote repeated experiments at different times.The difference between measurement and theory in the high-frequency range can be caused by the possibility of the system to function as a capacitor [8]. To prevent the capacitor effect, using insulating plates and electrodes perforated in them may be a solution. The impedance of the system can be determined using a two or four electrode method. The amplitude and phase of the impedance of the system can be determined and be used for data correction [8]. This can uplift particularly the phase values in the higher frequency range [8]. Also some offsets can be caused by the relative low permeability of the applied sample structure, especially the two Monel plates disturb the flow for high frequencies (which can also be observed in Figure 6). This limits the possibility of measuring samples with even lower permeability, in the current setup.


## 4. Seismoelectric Wave Propagation

### 4.1. Seismoelectric Wave Propagation Theory

Electrokinetic theory in isotropic, homogeneous, and fluid-saturated poroelastic media predicts the existence of a fast and a slowP-wave, a shear wave, and an electromagnetic wave. In this section, we derive wave speeds and attenuations (the dispersion relations) from the momentum equations (1), (2), and (18), for each of these waves. This derivation also yields the so-called fluid-to-solid and electric-to-solid field ratios. The fluid-to-solid ratio describes the fluid-to-solid displacement amplitude ratio, while the electric-to-solid field ratio describes the strength of the electric field with respect to the solid displacement field. These ratios and the dispersion relations are subsequently used to solve a boundary value problem and to formulate a full-waveform seismoelectric model.Using (2) to eliminate the electric field Ê from (1) and (18), we obtain two modified momentum equations for the fields ûs and ûf(21)G∇2ûs+(A+G+Q)∇(∇⋅ûs)+(Q+R)∇(∇⋅ûf)=-ω2[(ρ11(ω)+ρ12(ω))ûs+(ρ12(ω)+ρ22(ω))ûf],Q∇∇⋅ûs+R∇∇⋅ûf=-ω2(ρ¯12(ω)ûs+ρ¯22(ω)ûf)+ρ12(ω)μɛ¯(ω)(∇(∇⋅ûs)-∇2ûs)+ρ22(ω)μɛ¯(ω)(∇(∇⋅ûf)-∇2ûf),
where complex effective densities ρ¯11(ω), ρ¯12(ω), and ρ¯22(ω), containing the electrokinetic coupling factor EK(ω), are defined as follows
(22)ρ¯11(ω)=ρ11(ω)-EK(ω),ρ¯12(ω)=ρ12(ω)+EK(ω),ρ¯22(ω)=ρ22(ω)-EK(ω),EK(ω)=η2ϕ2L2(ω)ω2k2(ω)ɛ¯(ω).
Employing Helmholtz decomposition for the fields ûs and ûf leads to
(23)ûs=∇φ̂s+∇×Ψ̂s,ûf=∇φ̂f+∇×Ψ̂f.
Substituting expressions (23) into (21) yields
(24)∇[((P+Q)∇2+ω2(1-ϕ)ρs)φ̂s+((Q+R)∇2+ω2ϕρf)φ̂f]+∇×[(G∇2+ω2(1-ϕ)ρs)Ψ̂s+ω2ϕρfΨ̂f]=0,∇[(Q∇2+ω2ρ¯12(ω))φ̂s+(R∇2+ω2ρ¯22(ω))φ̂f]+∇×[(ω2ρ¯12(ω)+ρ12(ω)μɛ¯(ω)∇2)Ψ̂s+(ω2ρ¯22(ω)+ρ22(ω)μɛ¯(ω)∇2)Ψ̂f]=0,
where P=A+2G. For the longitudinal waves, associated with potentials φ̂s and φ̂f, the first terms in square brackets of (24) are set equal to zero from which we obtain
(25)(PQQR)∇2(φ̂sφ̂f)=-ω2(ρ¯11(ω)ρ¯12(ω)ρ¯12(ω)ρ¯22(ω))(φ̂sφ̂f),
where we used that (1-ϕ)ρs-ρ¯12(ω)=ρ¯11(ω), and ϕρf-ρ¯12(ω)=ρ¯22(ω). Applying a spatial Fourier transformation and recasting (25) into an eigenvalue problem lead to(26)1PR-Q2(ρ¯11(ω)R-ρ¯12(ω)Qρ¯12(ω)R-ρ¯22(ω)Qρ¯12(ω)P-ρ¯11(ω)Qρ¯22(ω)P-ρ¯12(ω)Q)(φ̃sφ̃f)=k⋅kω2(φ̃sφ̃f),where k is the wavenumber vector and tildes over a potential indicate frequency-wavenumber domain quantities. The complex eigenvalues correspond with the slownesses squared of the fast (Pf) and slow (Ps) longitudinal waves sk2(ω), k=Pf,Ps,
(27)sk2(ω)=-d1(ω)2d2∓d1(ω)2d21-4d0(ω)d2d12(ω),
where
(28)d0(ω)=ρ¯11(ω)ρ¯22(ω)-[ρ¯12(ω)]2,d1(ω)=-[ρ¯22(ω)P+ρ¯11(ω)R-2ρ¯12(ω)Q],d2=PR-Q2.
The slowness yields the wave mode speed and intrinsic attenuation (see, e.g., [67]). For the transversal waves, associated with potentials Ψ̃s and Ψ̃f, the second term in square brackets of (24) are set equal to the zero vector which gives(29)([-Gk⋅kω2+(1-ϕ)ρs]I=ϕρfI=[-ρ12(ω)μɛ¯(ω)k⋅kω2+ρ¯12(ω)]I=[-ρ22(ω)μɛ¯(ω)k⋅kω2+ρ¯22(ω)]I=)(Ψ̃sΨ̃f)=(00),where a spatial Fourier transformation is applied. Nontrivial solutions for k·k/ω2 are obtained by requiring the determinant of the matrix in (29) to be equal to zero. The solutions correspond with squared complex slownesses of the electromagnetic (EM) and seismic shear (S) waves. The dispersion relations are given in (27) for k=EM,S where
(30)d0(ω)=μɛ¯(ω)ρ¯11(ω)ρ¯22(ω)-[ρ¯12(ω)]2G,d1(ω)=-μɛ¯(ω)ρ¯22(ω)-ρ11(ω)ρ22(ω)-[ρ12(ω)]2G,d2(ω)=ρ22(ω),
and where we used that (1-ϕ)ρsρ22(ω)-ϕρfρ12(ω)=ρ11(ω)ρ22(ω)-[ρ12(ω)]2 and (1-ϕ)ρsρ¯22(ω)-ϕρfρ¯12(ω)=ρ¯11(ω)ρ¯22(ω)-[ρ¯12(ω)]2. Note that d2 in (27) is now frequency-dependent. Dispersion relations given by (27), (28), and (30) are equal to the expressions given by Pride and Haartsen [67].The longitudinal fluid-solid ratio, which describes the fluid-to-solid displacement amplitude ratio, is derived from the first row in (25). By applying a spatial Fourier transformation we obtain for the longitudinal fluid-solid ratios
(31)βm(ω)=φ̃mfφ̃ms=ρ¯11(ω)-Psm2(ω)Qsm2(ω)-ρ¯12(ω),
for m=Pf,Ps. By writing the vector potentials as Ψ̃ns = (0,ψ̃ns,0)T and Ψ̃nf=(0,ψ̃nf,0)T, for n=EM,S, in (29) we obtain for the transversal fluid-solid ratios
(32)βn(ω)=ψ̃nfψ̃ns=Gsn2(ω)-(1-ϕ)ρsϕρf.
The electric solid ratios, which describe the strength of the electric field with respect to the solid displacement field, are derived by applying Helmholtz decomposition (see (23)) to the fields in (18). This yields
(33)∇[ω2μɛ¯(ω)φ̂E+ω2ηϕμL(ω)k(ω)(φ̂f-φ̂s)]+∇×[(ω2μɛ¯(ω)+∇2)Ψ̂E+ω2ηϕμL(ω)k(ω)(Ψ̂f-Ψ̂s)]=0,
where we note that the Helmholtz decomposition of the electric field is Ê=∇φ̂E+∇×Ψ̂E, and Ψ̂nE=(0,ψ̃nE,0)T. Again, the scalar potentials are associated with longitudinal wave behavior and the vector potentials with transversal wave behavior. By applying a spatial Fourier transformation to (33) we obtain
(34)αm(ω)=φ̃mEφ̃ms=ηϕL(ω)k(ω)ɛ¯(ω)[1-βm(ω)],(35)αn(ω)=ψ̃nEψ̃ns=ηϕμL(ω)k(ω)[μɛ¯(ω)-sn2(ω)][1-βn(ω)].We now model coseismic electric potentials generated within a porous medium due to a fastP-wave, using its electric-solid ratio αPf(ω), for the geometry of Figure 9. The interface field responses are not modelled, which simplifies the expressions. In the forthcoming, it will be shown that this simplified model describes the measured coseismic electric potentials adequately. An acoustic wavefield from a source in a compressible fluid impinges on an interface between the fluid and an isotropic, homogeneous, and fluid-saturated poroelastic medium. It transmits as a fast P-wave in the poroelastic medium, where it generates coseismic electric potentials. We model a transducer (piezoelectric) source, as it is used in the experiment described in the following section. The acoustic pressure due to the transducer is modeled as (see [44, 75])
(36)p̂fl(ω,Rs,θ)=A(ω)D(θ)Rse-ikflRs,
where Rs=(r-rs)2+(z-zs)2 is the distance to the source, θ is the angle of incidence, A(ω) is the amplitude spectrum, and kfl=ωsP is the acoustic fluid wavenumber, where the fluid P-wave slowness is given by sP=ρf/Kf. The directivity function D(θ), which characterizes the radiation pattern of the source, is given by
(37)D(θ)=J1(kflasin⁡θ)kflasin⁡θ.
Here, J1 is the Bessel function of the first kind and first order and a is the radius of the transducer. Schakel et al. [44] show that seismoelectric effects can be modelled by expanding the source pressure wavefield into conical waves, which leads to the so-called Sommerfeld integral, and by relating acoustic potentials to electric signals with reflection/transmission coefficients as well as electric-solid ratios. While Schakel et al. [44] model both coseismic and interface field responses, we only model the coseismic fields. We arrive at the following Sommerfeld integral for the coseismic electric potential φ̂(ω,rr,zr) at receiver position (rr,zr) for rs=0, zs<0, zr>0,
(38)φ̂(ω,rr,zr)=-iA(ω)ω2ρf∫0∞krkzflD(kr)J0(krrr)eikzflzsαPf(ω)TPf(kr)(e-ikzPfzr+RPf(kr)e-ikzPf(2ws-zr))dkr,
where kr=kfl sinθ and kzfl=kfl cosθ are the radial and vertical components of kfl, respectively, and kzPf is the vertical component of the fast P-wave wavenumber. The fast P-wave wavenumber is kPf=ωsPf(ω), where the fast P-wave slowness is given by the dispersion relations (27)-(28). Note that the factor ω2ρf in the denominator of (38) is absent in Schakel et al. [44], because their reflection and transmission (conversion) coefficients are pressure normalized, whereas here they are displacement potential normalized. The transmission coefficient TPf(kr) relates the incident acoustic wavefield to the transmitted fast P-wave signal. The transmitted signal generates a coseismic potential at (r,z)=(rr,zr). It also reflects [RPf(kr)] at z=ws, and travels back to the receiver position, where it generates a second coseismic potential. The transmission coefficient TPf(kr) is derived from substituting plane wave expressions into the following (open-pore) boundary conditions [76]
(39)ûs,z+ŵz=ûzfl,p̂=p̂fl,σ̂13=σ̂33=0,
with subscript z denoting the z-component of the vectors and where ûfl denotes the fluid displacement. By only solving the mechanical (Biot) boundary value problem (no electrokinetic coupling is present in (39)), the interface field responses are neglected. Pride and Garambois [77] discussed the influence of the Biot slow wave in the generation of interface response seismoelectric amplitudes and numerically showed that when the Biot slow wave is neglected, the amplitudes can easily be off by as much as an order of magnitude. In our approach, that aims to model coseismic fields rather than interface responses, the Biot slow wave is taken into account in the solution of the boundary value problem (39). Its coseismic field is not modeled. For the parameters of Table 2, the slow wave skin depth is approximately 5 mm at 500 kHz and is unlikely to cause any appreciable coseismic signal for larger distances. We substitute the following plane wave expressions into (39)
(40)φ̂lfl=φ̃lfle-i(krrr±kzflz),φ̂ms=φ̃mse-i(krrr+kzmz),Ψ̂SVs=(0,ψ̃SVse-i(krrr+kzSVz),0)T,
for l=I,R. Hence we consider an incident (I) acoustic wave that reflects (R) and transmits as Pf-, Ps-, and vertical shear (SV)-waves. Displacement fields are obtained from these expressions as follows
(41)ûlfl=∇φ̂lfl,ûm=∇φ̂ms,ûSV=∇×Ψ̂SVs.
Fluid pressure is related to fluid displacement by p̂fl=-Kf∇·ûfl, with ûfl=ûIfl+ûRfl. For the poroelastic medium, solid displacement and pore-fluid displacement are obtained as follows
(42)ûs=ûs,Pf+ûs,Ps+ûs,SV,ûf=βPfûs,Pf+βPsûs,Ps+βSûs,SV.
Following the basic equations described in [60–64], the pore-fluid pressure and intergranular stresses are obtained. We define the reflection and transmission coefficients as
(43)RP=φ̃Rflφ̃Ifl,TPf=φ̃Pfsφ̃Ifl,TPs=φ̃Pssφ̃Ifl,TSV=ψ̃SVsφ̃Ifl,
so that we arrive at the following linear system of equations
(44)A⋅(RP,TPf,TPs,TSV)T=(kzfl,ϕρf,0,0)T,
where the elements of matrix A are given in the appendix. By solving (44) and (A.3) we obtain TPf and RPf, respectively (see appendix).Table 2
(Modified from [44].) Parameters of the poroelastic medium and fluid layer. The temperature is 293.15 K.
PropertySymbolValueUnitBulk modulus skeleton grainsaKfr50·109[Pa]Bulk modulus (pore) fluidbKf2.2·109[Pa]Bulk modulus framework of grainscKs0.93·109[Pa]Shear modulus framework of grainscG0.88·109[Pa]Pore fluid viscositybη1·10-3[kg/(m s)]Pore fluid densitybρf1000[kg/m3]Solid densitycρs2570[kg/m3]Weighted pore volume-to-surface ratiodΛ9.4·10-6[m]Porosity of the porous mediumcϕ0.52[—]Permeabilityck03.4·10-12[m2]Tortuositycα∞1.7[—]Sample widthews3.21·10-2[m]Relative permittivity of the (pore) fluidbϵrf80.1[—]Relative permittivity of the solidbϵrs4[—]Fluid magnetic permeability (=μ0)μf4π·10-7[H/m](Pore) fluid conductivityeσf4.8·10-2[S/m]Zeta-potentialfζ-4.0·10-2[V]a[78], b[68]. We take the value of Pyrex 7070 glass for the solid permittivity. csee N5B in [79], d[29, 59, 80, 81], emeasured values, and fsee [22]. We assume that conductivity is due to a NaCl salt solution and pH = 6.Figure 9
(Modified from [44].) Geometry of seismoelectric model and experiment. An acoustic wavefield emitted from (rs, zs) transmits into the poroelastic medium at z=0, where it generates coseismic electric potentials. Fluid/poroelastic-medium interfaces are at z=0 and z=3.21 cm. The electric receiver is at 1 cm from the front interface, while the dominant wavelength of the fast P-wave is roughly 4 mm. The electric potential recordings of the electric receiver are with respect to ground level.For the geometry of Figure9, where a source is located at (rs,zs)=(0,-15) cm, and where the receiver is located at (rr,zr)=(0,1) cm, we numerically evaluate the integral of (38). An experimentally recorded 500 kHz single sine pressure waveform is used for the amplitude spectrum A(ω). The incident pressure is related to the mechanical displacement potential ψ̃Pfs in the denominator of αPf(ω) (see (34)) by the factor ω2ρf, which arises from the relation p̂fl=-Kf∇·ûfl. The parameters of Table 2 are used and a 144–896 kHz numerical band-pass filter is applied. Figure 10(a) shows the resulting coseismic electric potentials caused by the fast P-wave. The first (CSP1) arrives at around 0.106 ms. This is the travel time of the acoustic wave from the source to the interface (approximately 0.101 ms) plus the travel time of the fast P-wave from the interface to the receiver location (approximately 0.005 ms). The predicted amplitude of the coseismic electric potential is approximately 0.5 mV, for an incident pressure amplitude of approximately 50 kPa. The second coseismic potential CSP2 arrives at around 0.130 ms and has an amplitude of approximately 0.15 mV. We conclude that coseismic electric potentials can be straightforwardly modelled in layered geometries by electric-solid ratios and solutions to mechanical boundary value problems.(Modified from [44].) Model of coseismic responses (a) and seismoelectric wave propagation measurements (b). Pulse abbreviations are: interface response (from the) front 1 (IRF1), coseismic response 1 (CS1), interface response (from the) back (IRB), coseismic response 2 (CS2), and interface response (from the) front 2 (IRF2).
(a)(b)


### 4.2. Seismoelectric Wave Propagation Experiment

Schakel et al. [44] report on a seismoelectric wave propagation experiment in which coseismic electric and interface field responses are measured. The results are reproduced in Figure 10(b). The geometry of the experiment is that of Figure 9. A 500 kHz single sine pulse generated by a waveform generator (Agilent Technologies 33220A) was used as input to the source. The second interface corresponds with the back of a porous sample. The receiver located at (rr,zr)=(0,1) cm recorded several pulses. The first (IRF1) is the interface response generated at the front (z=0) of the sample (see also Figure 2(c)). It arrives at around 0.100 ms, which corresponds with the acoustic wave travel time from the source to the interface. The travel time of the fast P-wave from the interface to the receiver location is approximately 0.005 ms. Therefore, the next pulse, labelled CS1, is the coseismic (electrical) response caused by the fast P-wave (see also Figure 2(a)). This wave also generates an interface response when it arrives at the back of the sample (IRB). It also reflects as a fast P-wave. When the reflected fast P-wave passes the receiver location for the second time, it generates another coseismic response (CS2). The last significant pulse, labelled (IRF2), is the interface response caused by the reflected fast P-wave when it arrives at the front of the sample. These experimental data were obtained using a 3.21 cm thick sample and a 500 kHz single sine pulse. It takes about 20 μs for the fast wave to arrive at the second interface, while the (measured) pulse period does not exceed 5 μs. Thus the pulses are clearly separated in time and do not cause amplitude and waveform changes.By comparing the model for the first coseismic response (CS1 in Figure10(a)) with the measurement (CS1 in Figure 10(b)) we observe agreement in travel time, waveform and amplitude. Small differences in waveform, such as the onset of the modelled waveform which is absent in the recording, are probably related to geometric misalignment and/or inaccuracies in the model/parameters (Table 2). The scale of Figure 10(a) is different from that of Figure 10(b). This amplitude difference is probably also related to geometric misalignment and/or inaccuracies in the model/parameters. For example, the model predictions are sensitive to the zeta-potential. This parameter was not directly measured but is obtained from an empirical relationship (see Table 2). For general field geometries the seismoelectric amplitudes of radiation generated at interfaces is significantly smaller than the coseismic amplitudes. For field geometries, electric receivers are typically positioned at several seismic wavelengths from the interfaces that generate seismoelectric conversion. In our experiment, the electric receiver is at 1 cm from the front interface, while the fast P-wave wavelength is roughly 4 mm. For this configuration, the measurements are as shown in Figure 10(b).The model for the second coseismic response (CS2 in Figure10(a)) shows less agreement with its corresponding measurement (CS2 in Figure 10(b)). We investigate the reason for this observation by matching the theory to the measurement for CS1. A frequency filter is constructed from the selected theoretical and measured CS1 pulses. This filter is subsequently applied to the selected theoretical CS2. The results are shown in Figure 11. The filtered theoretical CS1 fits the measurements exactly because it is forced to coincide with the measured CS1. The filtered theoretical CS2 now shows better agreement in terms of waveform and amplitude (Figures 11(a) and 11(b)). However there also remains to be mismatch, particularly the measured CS2 has its energy distributed over smaller frequencies than the filtered theoretical CS2. The latter fact is illustrated in Figures 11(c) and 11(d). The filtered theoretical CS2 differs from the filtered theoretical CS1 by the term RPf(kr)e-2ikzPf(ws-zr) (see (38)). Thus this observation indicates that the theory underpredicts the amount of seismic attenuation. It is well known that Biot’s theory can underestimate seismic attenuation [82]. However, the observation of Figure 11 could also be related to geometric misalignment in the experimental setup. We note that the possibility of underestimation by the electrokinetic coupling ratio α(ω) is excluded as it is effectively removed by the filter. In this paper we focus only on comparing theoretical and measured coseismic amplitudes rather than the seismoelectric responses at interfaces. The receiver is located at a constant distance from the interface, so that we do not compare the amplitudes to those generated by a (vertical) dipole located at the interface. A thorough comparison of seismoelectric amplitudes radiated from interfaces as a function of distance towards the interface with the pattern due to a dipole is given by [45].Model of filtered coseismic responses in time (a) and frequency domain (c) and coseismic measurements in time (b) and frequency domain (d). Pulse abbreviations are: coseismic response 1 (CS1) and coseismic response 2 (CS2).
(a)(b)(c)(d)It is possible to model all interface responses and coseismic effects of Figure10(b) by adopting full electrokinetic theory for the poroelastic medium in the boundary value problem [44]. This results in complicated expressions for the so-called seismoelectric reflection and conversion coefficients, which describe the interface responses, and also for the transmission coefficient TPf and RPf. Therefore, in the above, we only adopted Biot’s poroelastic theory to solve for TPf and used the electric-solid ratio αPf(ω) to describe the coseismic electric potential of Figure 10(a). The disadvantage of the approach is that interface response effects cannot be modelled. On the other hand, it results in simpler expressions for the coseismic fields.


## 4.1. Seismoelectric Wave Propagation Theory

Electrokinetic theory in isotropic, homogeneous, and fluid-saturated poroelastic media predicts the existence of a fast and a slowP-wave, a shear wave, and an electromagnetic wave. In this section, we derive wave speeds and attenuations (the dispersion relations) from the momentum equations (1), (2), and (18), for each of these waves. This derivation also yields the so-called fluid-to-solid and electric-to-solid field ratios. The fluid-to-solid ratio describes the fluid-to-solid displacement amplitude ratio, while the electric-to-solid field ratio describes the strength of the electric field with respect to the solid displacement field. These ratios and the dispersion relations are subsequently used to solve a boundary value problem and to formulate a full-waveform seismoelectric model.Using (2) to eliminate the electric field Ê from (1) and (18), we obtain two modified momentum equations for the fields ûs and ûf(21)G∇2ûs+(A+G+Q)∇(∇⋅ûs)+(Q+R)∇(∇⋅ûf)=-ω2[(ρ11(ω)+ρ12(ω))ûs+(ρ12(ω)+ρ22(ω))ûf],Q∇∇⋅ûs+R∇∇⋅ûf=-ω2(ρ¯12(ω)ûs+ρ¯22(ω)ûf)+ρ12(ω)μɛ¯(ω)(∇(∇⋅ûs)-∇2ûs)+ρ22(ω)μɛ¯(ω)(∇(∇⋅ûf)-∇2ûf),
where complex effective densities ρ¯11(ω), ρ¯12(ω), and ρ¯22(ω), containing the electrokinetic coupling factor EK(ω), are defined as follows
(22)ρ¯11(ω)=ρ11(ω)-EK(ω),ρ¯12(ω)=ρ12(ω)+EK(ω),ρ¯22(ω)=ρ22(ω)-EK(ω),EK(ω)=η2ϕ2L2(ω)ω2k2(ω)ɛ¯(ω).
Employing Helmholtz decomposition for the fields ûs and ûf leads to
(23)ûs=∇φ̂s+∇×Ψ̂s,ûf=∇φ̂f+∇×Ψ̂f.
Substituting expressions (23) into (21) yields
(24)∇[((P+Q)∇2+ω2(1-ϕ)ρs)φ̂s+((Q+R)∇2+ω2ϕρf)φ̂f]+∇×[(G∇2+ω2(1-ϕ)ρs)Ψ̂s+ω2ϕρfΨ̂f]=0,∇[(Q∇2+ω2ρ¯12(ω))φ̂s+(R∇2+ω2ρ¯22(ω))φ̂f]+∇×[(ω2ρ¯12(ω)+ρ12(ω)μɛ¯(ω)∇2)Ψ̂s+(ω2ρ¯22(ω)+ρ22(ω)μɛ¯(ω)∇2)Ψ̂f]=0,
where P=A+2G. For the longitudinal waves, associated with potentials φ̂s and φ̂f, the first terms in square brackets of (24) are set equal to zero from which we obtain
(25)(PQQR)∇2(φ̂sφ̂f)=-ω2(ρ¯11(ω)ρ¯12(ω)ρ¯12(ω)ρ¯22(ω))(φ̂sφ̂f),
where we used that (1-ϕ)ρs-ρ¯12(ω)=ρ¯11(ω), and ϕρf-ρ¯12(ω)=ρ¯22(ω). Applying a spatial Fourier transformation and recasting (25) into an eigenvalue problem lead to(26)1PR-Q2(ρ¯11(ω)R-ρ¯12(ω)Qρ¯12(ω)R-ρ¯22(ω)Qρ¯12(ω)P-ρ¯11(ω)Qρ¯22(ω)P-ρ¯12(ω)Q)(φ̃sφ̃f)=k⋅kω2(φ̃sφ̃f),where k is the wavenumber vector and tildes over a potential indicate frequency-wavenumber domain quantities. The complex eigenvalues correspond with the slownesses squared of the fast (Pf) and slow (Ps) longitudinal waves sk2(ω), k=Pf,Ps,
(27)sk2(ω)=-d1(ω)2d2∓d1(ω)2d21-4d0(ω)d2d12(ω),
where
(28)d0(ω)=ρ¯11(ω)ρ¯22(ω)-[ρ¯12(ω)]2,d1(ω)=-[ρ¯22(ω)P+ρ¯11(ω)R-2ρ¯12(ω)Q],d2=PR-Q2.
The slowness yields the wave mode speed and intrinsic attenuation (see, e.g., [67]). For the transversal waves, associated with potentials Ψ̃s and Ψ̃f, the second term in square brackets of (24) are set equal to the zero vector which gives(29)([-Gk⋅kω2+(1-ϕ)ρs]I=ϕρfI=[-ρ12(ω)μɛ¯(ω)k⋅kω2+ρ¯12(ω)]I=[-ρ22(ω)μɛ¯(ω)k⋅kω2+ρ¯22(ω)]I=)(Ψ̃sΨ̃f)=(00),where a spatial Fourier transformation is applied. Nontrivial solutions for k·k/ω2 are obtained by requiring the determinant of the matrix in (29) to be equal to zero. The solutions correspond with squared complex slownesses of the electromagnetic (EM) and seismic shear (S) waves. The dispersion relations are given in (27) for k=EM,S where
(30)d0(ω)=μɛ¯(ω)ρ¯11(ω)ρ¯22(ω)-[ρ¯12(ω)]2G,d1(ω)=-μɛ¯(ω)ρ¯22(ω)-ρ11(ω)ρ22(ω)-[ρ12(ω)]2G,d2(ω)=ρ22(ω),
and where we used that (1-ϕ)ρsρ22(ω)-ϕρfρ12(ω)=ρ11(ω)ρ22(ω)-[ρ12(ω)]2 and (1-ϕ)ρsρ¯22(ω)-ϕρfρ¯12(ω)=ρ¯11(ω)ρ¯22(ω)-[ρ¯12(ω)]2. Note that d2 in (27) is now frequency-dependent. Dispersion relations given by (27), (28), and (30) are equal to the expressions given by Pride and Haartsen [67].The longitudinal fluid-solid ratio, which describes the fluid-to-solid displacement amplitude ratio, is derived from the first row in (25). By applying a spatial Fourier transformation we obtain for the longitudinal fluid-solid ratios
(31)βm(ω)=φ̃mfφ̃ms=ρ¯11(ω)-Psm2(ω)Qsm2(ω)-ρ¯12(ω),
for m=Pf,Ps. By writing the vector potentials as Ψ̃ns = (0,ψ̃ns,0)T and Ψ̃nf=(0,ψ̃nf,0)T, for n=EM,S, in (29) we obtain for the transversal fluid-solid ratios
(32)βn(ω)=ψ̃nfψ̃ns=Gsn2(ω)-(1-ϕ)ρsϕρf.
The electric solid ratios, which describe the strength of the electric field with respect to the solid displacement field, are derived by applying Helmholtz decomposition (see (23)) to the fields in (18). This yields
(33)∇[ω2μɛ¯(ω)φ̂E+ω2ηϕμL(ω)k(ω)(φ̂f-φ̂s)]+∇×[(ω2μɛ¯(ω)+∇2)Ψ̂E+ω2ηϕμL(ω)k(ω)(Ψ̂f-Ψ̂s)]=0,
where we note that the Helmholtz decomposition of the electric field is Ê=∇φ̂E+∇×Ψ̂E, and Ψ̂nE=(0,ψ̃nE,0)T. Again, the scalar potentials are associated with longitudinal wave behavior and the vector potentials with transversal wave behavior. By applying a spatial Fourier transformation to (33) we obtain
(34)αm(ω)=φ̃mEφ̃ms=ηϕL(ω)k(ω)ɛ¯(ω)[1-βm(ω)],(35)αn(ω)=ψ̃nEψ̃ns=ηϕμL(ω)k(ω)[μɛ¯(ω)-sn2(ω)][1-βn(ω)].We now model coseismic electric potentials generated within a porous medium due to a fastP-wave, using its electric-solid ratio αPf(ω), for the geometry of Figure 9. The interface field responses are not modelled, which simplifies the expressions. In the forthcoming, it will be shown that this simplified model describes the measured coseismic electric potentials adequately. An acoustic wavefield from a source in a compressible fluid impinges on an interface between the fluid and an isotropic, homogeneous, and fluid-saturated poroelastic medium. It transmits as a fast P-wave in the poroelastic medium, where it generates coseismic electric potentials. We model a transducer (piezoelectric) source, as it is used in the experiment described in the following section. The acoustic pressure due to the transducer is modeled as (see [44, 75])
(36)p̂fl(ω,Rs,θ)=A(ω)D(θ)Rse-ikflRs,
where Rs=(r-rs)2+(z-zs)2 is the distance to the source, θ is the angle of incidence, A(ω) is the amplitude spectrum, and kfl=ωsP is the acoustic fluid wavenumber, where the fluid P-wave slowness is given by sP=ρf/Kf. The directivity function D(θ), which characterizes the radiation pattern of the source, is given by
(37)D(θ)=J1(kflasin⁡θ)kflasin⁡θ.
Here, J1 is the Bessel function of the first kind and first order and a is the radius of the transducer. Schakel et al. [44] show that seismoelectric effects can be modelled by expanding the source pressure wavefield into conical waves, which leads to the so-called Sommerfeld integral, and by relating acoustic potentials to electric signals with reflection/transmission coefficients as well as electric-solid ratios. While Schakel et al. [44] model both coseismic and interface field responses, we only model the coseismic fields. We arrive at the following Sommerfeld integral for the coseismic electric potential φ̂(ω,rr,zr) at receiver position (rr,zr) for rs=0, zs<0, zr>0,
(38)φ̂(ω,rr,zr)=-iA(ω)ω2ρf∫0∞krkzflD(kr)J0(krrr)eikzflzsαPf(ω)TPf(kr)(e-ikzPfzr+RPf(kr)e-ikzPf(2ws-zr))dkr,
where kr=kfl sinθ and kzfl=kfl cosθ are the radial and vertical components of kfl, respectively, and kzPf is the vertical component of the fast P-wave wavenumber. The fast P-wave wavenumber is kPf=ωsPf(ω), where the fast P-wave slowness is given by the dispersion relations (27)-(28). Note that the factor ω2ρf in the denominator of (38) is absent in Schakel et al. [44], because their reflection and transmission (conversion) coefficients are pressure normalized, whereas here they are displacement potential normalized. The transmission coefficient TPf(kr) relates the incident acoustic wavefield to the transmitted fast P-wave signal. The transmitted signal generates a coseismic potential at (r,z)=(rr,zr). It also reflects [RPf(kr)] at z=ws, and travels back to the receiver position, where it generates a second coseismic potential. The transmission coefficient TPf(kr) is derived from substituting plane wave expressions into the following (open-pore) boundary conditions [76]
(39)ûs,z+ŵz=ûzfl,p̂=p̂fl,σ̂13=σ̂33=0,
with subscript z denoting the z-component of the vectors and where ûfl denotes the fluid displacement. By only solving the mechanical (Biot) boundary value problem (no electrokinetic coupling is present in (39)), the interface field responses are neglected. Pride and Garambois [77] discussed the influence of the Biot slow wave in the generation of interface response seismoelectric amplitudes and numerically showed that when the Biot slow wave is neglected, the amplitudes can easily be off by as much as an order of magnitude. In our approach, that aims to model coseismic fields rather than interface responses, the Biot slow wave is taken into account in the solution of the boundary value problem (39). Its coseismic field is not modeled. For the parameters of Table 2, the slow wave skin depth is approximately 5 mm at 500 kHz and is unlikely to cause any appreciable coseismic signal for larger distances. We substitute the following plane wave expressions into (39)
(40)φ̂lfl=φ̃lfle-i(krrr±kzflz),φ̂ms=φ̃mse-i(krrr+kzmz),Ψ̂SVs=(0,ψ̃SVse-i(krrr+kzSVz),0)T,
for l=I,R. Hence we consider an incident (I) acoustic wave that reflects (R) and transmits as Pf-, Ps-, and vertical shear (SV)-waves. Displacement fields are obtained from these expressions as follows
(41)ûlfl=∇φ̂lfl,ûm=∇φ̂ms,ûSV=∇×Ψ̂SVs.
Fluid pressure is related to fluid displacement by p̂fl=-Kf∇·ûfl, with ûfl=ûIfl+ûRfl. For the poroelastic medium, solid displacement and pore-fluid displacement are obtained as follows
(42)ûs=ûs,Pf+ûs,Ps+ûs,SV,ûf=βPfûs,Pf+βPsûs,Ps+βSûs,SV.
Following the basic equations described in [60–64], the pore-fluid pressure and intergranular stresses are obtained. We define the reflection and transmission coefficients as
(43)RP=φ̃Rflφ̃Ifl,TPf=φ̃Pfsφ̃Ifl,TPs=φ̃Pssφ̃Ifl,TSV=ψ̃SVsφ̃Ifl,
so that we arrive at the following linear system of equations
(44)A⋅(RP,TPf,TPs,TSV)T=(kzfl,ϕρf,0,0)T,
where the elements of matrix A are given in the appendix. By solving (44) and (A.3) we obtain TPf and RPf, respectively (see appendix).Table 2
(Modified from [44].) Parameters of the poroelastic medium and fluid layer. The temperature is 293.15 K.
PropertySymbolValueUnitBulk modulus skeleton grainsaKfr50·109[Pa]Bulk modulus (pore) fluidbKf2.2·109[Pa]Bulk modulus framework of grainscKs0.93·109[Pa]Shear modulus framework of grainscG0.88·109[Pa]Pore fluid viscositybη1·10-3[kg/(m s)]Pore fluid densitybρf1000[kg/m3]Solid densitycρs2570[kg/m3]Weighted pore volume-to-surface ratiodΛ9.4·10-6[m]Porosity of the porous mediumcϕ0.52[—]Permeabilityck03.4·10-12[m2]Tortuositycα∞1.7[—]Sample widthews3.21·10-2[m]Relative permittivity of the (pore) fluidbϵrf80.1[—]Relative permittivity of the solidbϵrs4[—]Fluid magnetic permeability (=μ0)μf4π·10-7[H/m](Pore) fluid conductivityeσf4.8·10-2[S/m]Zeta-potentialfζ-4.0·10-2[V]a[78], b[68]. We take the value of Pyrex 7070 glass for the solid permittivity. csee N5B in [79], d[29, 59, 80, 81], emeasured values, and fsee [22]. We assume that conductivity is due to a NaCl salt solution and pH = 6.Figure 9
(Modified from [44].) Geometry of seismoelectric model and experiment. An acoustic wavefield emitted from (rs, zs) transmits into the poroelastic medium at z=0, where it generates coseismic electric potentials. Fluid/poroelastic-medium interfaces are at z=0 and z=3.21 cm. The electric receiver is at 1 cm from the front interface, while the dominant wavelength of the fast P-wave is roughly 4 mm. The electric potential recordings of the electric receiver are with respect to ground level.For the geometry of Figure9, where a source is located at (rs,zs)=(0,-15) cm, and where the receiver is located at (rr,zr)=(0,1) cm, we numerically evaluate the integral of (38). An experimentally recorded 500 kHz single sine pressure waveform is used for the amplitude spectrum A(ω). The incident pressure is related to the mechanical displacement potential ψ̃Pfs in the denominator of αPf(ω) (see (34)) by the factor ω2ρf, which arises from the relation p̂fl=-Kf∇·ûfl. The parameters of Table 2 are used and a 144–896 kHz numerical band-pass filter is applied. Figure 10(a) shows the resulting coseismic electric potentials caused by the fast P-wave. The first (CSP1) arrives at around 0.106 ms. This is the travel time of the acoustic wave from the source to the interface (approximately 0.101 ms) plus the travel time of the fast P-wave from the interface to the receiver location (approximately 0.005 ms). The predicted amplitude of the coseismic electric potential is approximately 0.5 mV, for an incident pressure amplitude of approximately 50 kPa. The second coseismic potential CSP2 arrives at around 0.130 ms and has an amplitude of approximately 0.15 mV. We conclude that coseismic electric potentials can be straightforwardly modelled in layered geometries by electric-solid ratios and solutions to mechanical boundary value problems.(Modified from [44].) Model of coseismic responses (a) and seismoelectric wave propagation measurements (b). Pulse abbreviations are: interface response (from the) front 1 (IRF1), coseismic response 1 (CS1), interface response (from the) back (IRB), coseismic response 2 (CS2), and interface response (from the) front 2 (IRF2).
(a)(b)


## 4.2. Seismoelectric Wave Propagation Experiment

Schakel et al. [44] report on a seismoelectric wave propagation experiment in which coseismic electric and interface field responses are measured. The results are reproduced in Figure 10(b). The geometry of the experiment is that of Figure 9. A 500 kHz single sine pulse generated by a waveform generator (Agilent Technologies 33220A) was used as input to the source. The second interface corresponds with the back of a porous sample. The receiver located at (rr,zr)=(0,1) cm recorded several pulses. The first (IRF1) is the interface response generated at the front (z=0) of the sample (see also Figure 2(c)). It arrives at around 0.100 ms, which corresponds with the acoustic wave travel time from the source to the interface. The travel time of the fast P-wave from the interface to the receiver location is approximately 0.005 ms. Therefore, the next pulse, labelled CS1, is the coseismic (electrical) response caused by the fast P-wave (see also Figure 2(a)). This wave also generates an interface response when it arrives at the back of the sample (IRB). It also reflects as a fast P-wave. When the reflected fast P-wave passes the receiver location for the second time, it generates another coseismic response (CS2). The last significant pulse, labelled (IRF2), is the interface response caused by the reflected fast P-wave when it arrives at the front of the sample. These experimental data were obtained using a 3.21 cm thick sample and a 500 kHz single sine pulse. It takes about 20 μs for the fast wave to arrive at the second interface, while the (measured) pulse period does not exceed 5 μs. Thus the pulses are clearly separated in time and do not cause amplitude and waveform changes.By comparing the model for the first coseismic response (CS1 in Figure10(a)) with the measurement (CS1 in Figure 10(b)) we observe agreement in travel time, waveform and amplitude. Small differences in waveform, such as the onset of the modelled waveform which is absent in the recording, are probably related to geometric misalignment and/or inaccuracies in the model/parameters (Table 2). The scale of Figure 10(a) is different from that of Figure 10(b). This amplitude difference is probably also related to geometric misalignment and/or inaccuracies in the model/parameters. For example, the model predictions are sensitive to the zeta-potential. This parameter was not directly measured but is obtained from an empirical relationship (see Table 2). For general field geometries the seismoelectric amplitudes of radiation generated at interfaces is significantly smaller than the coseismic amplitudes. For field geometries, electric receivers are typically positioned at several seismic wavelengths from the interfaces that generate seismoelectric conversion. In our experiment, the electric receiver is at 1 cm from the front interface, while the fast P-wave wavelength is roughly 4 mm. For this configuration, the measurements are as shown in Figure 10(b).The model for the second coseismic response (CS2 in Figure10(a)) shows less agreement with its corresponding measurement (CS2 in Figure 10(b)). We investigate the reason for this observation by matching the theory to the measurement for CS1. A frequency filter is constructed from the selected theoretical and measured CS1 pulses. This filter is subsequently applied to the selected theoretical CS2. The results are shown in Figure 11. The filtered theoretical CS1 fits the measurements exactly because it is forced to coincide with the measured CS1. The filtered theoretical CS2 now shows better agreement in terms of waveform and amplitude (Figures 11(a) and 11(b)). However there also remains to be mismatch, particularly the measured CS2 has its energy distributed over smaller frequencies than the filtered theoretical CS2. The latter fact is illustrated in Figures 11(c) and 11(d). The filtered theoretical CS2 differs from the filtered theoretical CS1 by the term RPf(kr)e-2ikzPf(ws-zr) (see (38)). Thus this observation indicates that the theory underpredicts the amount of seismic attenuation. It is well known that Biot’s theory can underestimate seismic attenuation [82]. However, the observation of Figure 11 could also be related to geometric misalignment in the experimental setup. We note that the possibility of underestimation by the electrokinetic coupling ratio α(ω) is excluded as it is effectively removed by the filter. In this paper we focus only on comparing theoretical and measured coseismic amplitudes rather than the seismoelectric responses at interfaces. The receiver is located at a constant distance from the interface, so that we do not compare the amplitudes to those generated by a (vertical) dipole located at the interface. A thorough comparison of seismoelectric amplitudes radiated from interfaces as a function of distance towards the interface with the pattern due to a dipole is given by [45].Model of filtered coseismic responses in time (a) and frequency domain (c) and coseismic measurements in time (b) and frequency domain (d). Pulse abbreviations are: coseismic response 1 (CS1) and coseismic response 2 (CS2).
(a)(b)(c)(d)It is possible to model all interface responses and coseismic effects of Figure10(b) by adopting full electrokinetic theory for the poroelastic medium in the boundary value problem [44]. This results in complicated expressions for the so-called seismoelectric reflection and conversion coefficients, which describe the interface responses, and also for the transmission coefficient TPf and RPf. Therefore, in the above, we only adopted Biot’s poroelastic theory to solve for TPf and used the electric-solid ratio αPf(ω) to describe the coseismic electric potential of Figure 10(a). The disadvantage of the approach is that interface response effects cannot be modelled. On the other hand, it results in simpler expressions for the coseismic fields.


## 5. Seismoelectric Interferometry

Considering the combined character of seismo-electromagnetic waves it can be very beneficial to use them for a wide range of applications. (The application for oil-field exploration has already been shown by Thompson et al. [83].) From an imaging point of view, the principle of interferometry has already been proven useful for a wide class of phenomena, for example in seismic systems or electromagnetic systems (e.g., [47, 48]). Hence, we are taking this principle to the next level: correlation imaging with seismo-electromagnetic waves. Before showing some examples, the principle of interferometry will be explained first.


### 5.1. Theory

Interferometry makes use of the cross-correlation of responses at different receivers in order to obtain the Green’s function of the field response between these stations. In other words, it is the deterministic response from one station to the other.Figure12 shows a possible seismoelectric interferometry setting. The cross-correlation of electric (Ex) and acoustic signals (vz) from sources located at the surface (Figure 12(a)) or in the bulk (Figure 12(c)) results in the direct electric response of an acoustic source (fz) generating a seismoelectric wave (Figure 12(b)). The known challenging problems in using seismo-electromagnetics as a geophysical exploration tool can potentially be addressed by applying interferometric Green’s function retrieval techniques to seismo-electromagnetic phenomena [46]. First of all, sources in “classical” seismoelectric surveys need to be strong. This is not always possible and therefore it is beneficial to be able to replace those strong sources by receivers.Figure 12
Schematic seismoelectric interferometry setting. Cross-correlation of electric (Ex) and acoustic signals (vz) from sources located at the surface (a) or in the bulk (c) results in the direct electric response of an acoustic source (fz) generating a seismoelectric wave (b).A second well-known problem in these conventional seismoelectric surveys is the very low signal-to-noise ratio. By doing interferometry, stacking inherently takes place with a possible improvement of the signal-to-noise ratio as a result. After deriving the system of equations for coupled seismic and electromagnetic waves in saturated porous media [29], the convolution-type reciprocity theorem and a power balance for seismoelectric waves was derived by Pride and Haartsen [67]. In 2003, this result was extended to a reciprocity theorem of the correlation-type for seismoelectric waves [84].Following Wapenaar and Fokkema [85], de Ridder et al. [46] showed that the 1D seismoelectric system for the SH-TE propagation mode can be captured in the following matrix-vector equation(45)iωÂû+Dxû=ŝ,
where matrix Â contains the space-dependent material parameters, û represents the field vector (in the space-frequency domain), matrix Dx contains the spatial differential operator ∂/∂x, ŝ denotes the source vector, and where iω arises due to Fourier transformation of the temporal derivative of a field. It is important to capture the 1D seismoelectric system in such a general diffusion, flow and wave equation, in order to employ the derived expressions for unified Green’s function retrieval by cross-correlation [47] and to finally end up with interferometric seismoelectric Green’s function representations.Next, considering the Fourier transform of an impulsive source acting at timet=0 s and at location x=xs, ŝ in (45) is replaced by Iδ(x-xs). As a consequence, the field vector û can be replaced by a Green’s matrix Ĝ(x,xs,ω). In this way, (45) is changed to(46)iωÂĜ(x,xs,ω)+DxĜ(x,xs,ω)=Iδ(x-xs),
where the Green’s matrix Ĝ(x,xs,ω) is given by(47)Ĝ(x,xs,ω)=(ĜE,Je(x,xs,ω)ĜE,Jm(x,xs,ω)ĜE,f(x,xs,ω)ĜE,hb(x,xs,ω)ĜH,Je(x,xs,ω)ĜH,Jm(x,xs,ω)ĜH,f(x,xs,ω)ĜH,hb(x,xs,ω)Ĝvs,Je(x,xs,ω)Ĝvs,Jm(x,xs,ω)Ĝvs,f(x,xs,ω)Ĝvs,hb(x,xs,ω)Ĝτb,Je(x,xs,ω)Ĝτb,Jm(x,xs,ω)Ĝτb,f(x,xs,ω)Ĝτb,hb(x,xs,ω)).The first superscript (vs) in Ĝvs,Je(x,xs,ω) denotes the type of response measured at location x, resulting from the type of impulsive source located at xs, which is denoted by the second superscript (Je).Starting from the general interferometric Green’s function representation (48) as derived by Wapenaar et al. [47](48)Ĝ(xB,xA,ω)+Ĝ†(xA,xB,ω)=-{Ĝ(xB,x,ω)NxĜ†(xA,x,ω)}|x1x2+∫D{Ĝ(xB,x,ω)[iω(Â-Â†)]Ĝ†(xA,x,ω)}dx,
where it is assumed that the two reciprocity states have the same medium parameters and where Nx represents the normal vector matrix containing the components of a normal vector n, arranged in the same way as the partial spatial derivative ∂/∂x in the matrix Dx(49)Nx=(0-100-100000010010),
de Ridder et al. [46] derived the following interferometric integral representation for one element of the seismoelectric SH-TE Green’s matrix in 1D (50), using seismoelectric reciprocity theorems(50){ĜE,f(xB,xA,ω)+Ĝvs,Je*(xA,xB,ω)}Ŝ={ĜE,Je(xB,x,ω)Ĝvs,Jm*(xA,x,ω)+ĜE,Jm(xB,x,ω)Ĝvs,Je*(xA,x,ω)-ĜE,f(xB,x,ω)Ĝvs,h*(xA,x,ω)-ĜE,h(xB,x,ω)Ĝvs,f*(xA,x,ω)}|x1x2Ŝ+2iω∫D{ĜE,Je(xB,x,ω)iI{ϵ̂}Ĝvs,Je*(xA,x,ω)+ĜE,f(xB,x,ω)R{ρfL̂0}Ĝvs,Je*(xA,x,ω)-ĜE,Je(xB,x,ω)R{ρfL̂0}Ĝvs,f*(xA,x,ω)+ĜE,f(xB,x,ω)iI{ρ̂c}Ĝvs,f*(xA,x,ω)}dxŜ.
To arrive at this form, they have chosen the {1,3} element of the 1D SH-TE seismoelectric Green’s matrix (47) and expanded (48) using this element. Here, Ŝ denotes the power spectrum of the emitted source signal and ρ̂c=ρb-iωk0(ρf)2/η.We can distinguish two terms in this integral representation. The first term on the right-hand side represents correlations of recorded responses of sources on the boundary of the domain of reciprocity, whereas the second term on the righthand-side represents correlations of recorded responses of sources throughout the reciprocity domain.As shown by de Ridder et al. [46], the following source-receiver reciprocity holds(51)Ĝvs,Je(xA,xB,ω)=-ĜE,f(xB,xA,ω).
Hence, the left-hand side of (50) can be rewritten as 2iℑ{ĜE,f(xB,xA,ω)}Ŝ. This signal will be antisymmetric around t=0 s in the time-domain.Looking at expression (50) in more detail it can be seen that the left-hand side, the electric field response registered at xB generated by an elastic force source located at xA, is obtained by cross-correlating the registered electric fields at xB with the registered particle velocities at xA, which are the result of four different types of boundary sources and two types of domain sources. The two types of domain sources, an electric current source and an elastic force source, are both weighted with two different medium parameters.Due to the fact that wave energy is dissipated during wave propagation, the domain sources are necessary to account for these losses. However, these sources are not likely to exist in reality or cannot be rewritten for practical applications and therefore we would like to be able to ignore their contributions.As is already shown in three examples by de Ridder et al. [46], it is indeed possible (under certain conditions) to obtain accurate Green’s functions from boundary sources only. The most complex situation considered by de Ridder et al. [46] was a medium consisting of two layers bounded by a vacuum. For this situation it was shown that the domain integral contribution could be neglected as long as the domain of reciprocity was chosen in such a way, that it included the heterogeneities (i.e., the interface between the two layers). Then, spurious events would only occur on one side of the symmetrized, retrieved Green’s function.In the following section, we will increase the complexity of the numerical configuration by adding an extra layer to the system, to investigate the Green’s function retrieval for a 1D, three-layered system bounded by a free-surface. In other words, we will look at the applicability of the interferometric seismoelectric Green’s function representation (50) when there are two interfaces located in the subsurface. We will consider a configuration where a medium B is sandwiched between two identical layers (medium A) with different medium parameters, as given in Table 3. We have chosen the medium parameters in such a way, that there is a very small seismic contrast between the layers; the porosity is the only contrasting seismic parameter between the layers. In this way, we are minimizing the dominant coseismic field response from the two subsurface interfaces and are able to focus mainly on the retrieval of the interface response field. The free-surface on the other hand, acts as a reflector for the seismic waves and therefore the coseismic fields related to this interface are still preserved. In this way, we are able to investigate separately the retrieval of both the coseismic field responses and the interface response fields. For field geometries, the amplitudes of the coseismic field responses related to the subsurface interfaces are often much higher than the interface response field amplitudes.Table 3
Overview of the relevant medium parameters for the 1D seismoelectric interferometry model.
PropertyUnitValue medium AValue medium BDimensionPorosityϕ0.40.2[—]Pore fluid densityρf1.0·1031.0·103[kg/m3]Solid densityρs2.7·1032.7·103[kg/m3]Shear modulus framework of grainsG9.0·1099.0·109[Pa]Pore fluid viscosityη1.0·10-31.0·10-3[kg/(m s)]Static permeabilityk01.3·10-121.6·10-12[m2]Static electrokinetic couplingℒ01.0·10-81.0·10-9[m2/(s V)]Tortuosityα∞3.03.0[—]Relative perm. of the (pore) fluidϵrf8080[—]Relative perm. of the solidϵrs44[—]Bulk electric conductivityσ01.0·10-11.0·10-1[S/m]


### 5.2. Results

We consider a three-layered 1D medium bounded by a vacuum half-space. The top and bottom layer consist of medium parameters belonging to medium A and the sandwiched layer has the properties of medium B (see Table3). The bottom layer is in fact a half-space. The whole three-layered system is bounded by a vacuum half-space in which only electromagnetic waves can propagate. The interface separating the subsurface from the vacuum is called the free-surface. According to Wapenaar and Fokkema [86], the free-surface acts as a mirror to both shear waves and electromagnetic waves (the latter due to the fact that a 1D geometry is considered here). Therefore, the sources on the domain boundary at the free surface can be neglected and also the contributions of the vacuum above the free-surface can be disregarded. The range of the domain integral contribution is from 0 to 650 m (see Figure 13 for an overview of the geometry). The receivers at xA and xB are located at x=100m and x=275m, respectively. The upper boundary is called x1 and the lower boundary x2, located at x=-150m and x=650m, respectively.Figure 13
The geometry of the 1D numerical experiment. PositionsxA and xB represent the receiver positions, denoted by the top-down triangles. The upper boundary is called x1 and the lower boundary x2, located at x=-150m and x=650m, respectively. The positions of the interfaces are visualized by the red bars. FS denotes the free-surface, whereas IF1 and IF2 correspond to the first and second subsurface interfaces, respectively. Furthermore, schematic ray paths of events in the seismoelectric Green’s function GE,f(xB,xA,t) are shown. The black arrows represent shear wave ray paths, whereas the green arrows depict electromagnetic ray paths. The labels correspond to the labelled events in Figures 14 and 15, where d denotes the direct event and dg its source-side ghost, which are both coseismic field responses, 1 represents the interface response of the most shallow interface and 1g its source-side ghost, 2 corresponds to the interface response of the second, deeper interface and event 2g represents again its source-side ghost.Figure14 shows the time-domain equivalent of the separated contributions of the domain integral and the boundary points to the retrieved Green’s functions. In other words, it shows the relative contributions of the two right-hand side terms in (50) to the retrieved Green’s functions. The positive time corresponds to the Green’s function GE,f(xB,xA,t), the electric field response registered at xB due to an impulsive seismic source located at xA. As is visible, the dominant contribution in the positive time window comes from just the boundary term. Therefore, it is shown that this Green’s function can be mainly reconstructed by using the boundary contribution only. In contrast, the negative time window contains strong domain integral contributions as well. The negative times correspond to the Green’s function Gvs,Je(xA,xB,-t), the particle velocity response measured at xA due to an impulsive electrical current source at xB. These strong volume source contributions correct the polarity of the single boundary term contribution.Figure 14
Separated contributions of the domain integral and the boundary points to the retrieved Green’s functions. In other words, it shows the relative contributions of the two right-hand side terms in (50) to the retrieved Green’s functions, for a three-layered medium bounded by a vacuum.Several events can be recognized in Figure14. The purely diffusive electromagnetic field is the first event to arrive, with its maximum at approximately t=0.58 ms (t=-0.58 ms for the time-reversed causal signal). The second arrival, at approximately t=83 ms corresponds to a direct coseismic shear wave event (labelled d). Its time-reversed causal equivalent arrives around t=-83 ms. The source-side ghost of this direct coseismic shear wave event (labelled dg) arrives at approximately t=0.18 s. The overlapping causal and time-reversed causal electromagnetic events at t=0 s are constructed completely by the sources in the domain integral. In contrast, the shear wave event is retrieved by mainly boundary source contributions. This makes sense, considering the sources of wave energy loss. As mentioned already, the electromagnetic event is, in the considered seismic frequency range, primarily a diffusive field. Therefore, volume sources are required to compensate for the wave energy loss. For the shear wave event, the amount of wave energy loss is relatively small. Hence, the need for volume source energy compensation is negligible. Considering this in terms of the interferometric seismoelectric integral representation (50),(52){ĜE,f(xB,xA,ω)+Ĝvs,Je*(xA,xB,ω)}Ŝ={ĜE,Je(xB,x,ω)Ĝvs,Jm*(xA,x,ω)+ĜE,Jm(xB,x,ω)Ĝvs,Je*(xA,x,ω)-ĜE,f(xB,x,ω)Ĝvs,h*(xA,x,ω)-ĜE,h(xB,x,ω)Ĝvs,f*(xA,x,ω)}|x1x2Ŝ+2iω∫D{ĜE,Je(xB,x,ω)iI{ϵ̂}Ĝvs,Je*(xA,x,ω)+ĜE,f(xB,x,ω)R{ρfL̂0}Ĝvs,Je*(xA,x,ω)-ĜE,Je(xB,x,ω)R{ρfL̂0}Ĝvs,f*(xA,x,ω)+ĜE,f(xB,x,ω)iI{ρ̂c}Ĝvs,f*(xA,x,ω)}dxŜ,
we can omit the volume source contributions(53)2iω∫D{ĜE,Je(xB,x,ω)iI{ϵ̂}Ĝvs,Je*(xA,x,ω)+ĜE,f(xB,x,ω)R{ρfL̂0}Ĝvs,Je*(xA,x,ω)-ĜE,Je(xB,x,ω)R{ρfL̂0}Ĝvs,f*(xA,x,ω)+ĜE,f(xB,x,ω)iI{ρ̂c}Ĝvs,f*(xA,x,ω)}dxŜ≈0.
The free-surface acts as a mirror to both shear waves and electromagnetic waves and therefore the sources on the domain boundary at the free-surface can be neglected and also the contributions of the vacuum medium above the free-surface can be disregarded. For (52), this effectively means that the contributions of the boundary sources at x1 can be omitted. The first two terms on the right-hand side of (52) represent the correlation products of the two fields generated by either electrical or magnetic current sources. The contribution of these two correlation products to the positive time window is very small, due to the fact that the causal fields registered at xB are electric fields which, without any wavetype conversion, only contribute significantly at t=0 (they arrive instantaneously). The contributions of the electric signals with positive seismic traveltimes are negligible, because these signals have encountered at least two wavetype conversions; this implies significant energy losses. Therefore, we can additionally neglect the contributions of the electromagnetic boundary sources (Je and Jm), ending up with the following reduced interferometric seismoelectric integral representation(54){ĜE,f(xB,xA,ω)+Ĝvs,Je*(xA,xB,ω)}Ŝ={-ĜE,f(xB,x2,ω)Ĝvs,h*(xA,x2,ω)-ĜE,h(xB,x2,ω)Ĝvs,f*(xA,x2,ω)}Ŝ.
Figure 16 shows the result of using this reduced interferometric seismoelectric integral representation. The figure displays a comparison between the exact Green’s function in the positive time window GE,f(xB,xA,t) and the Green’s function retrieved by using (54), that is, by considering only seismic boundary source contributions (f and h) at x2. As is visible, the amplitude errors between the exact and retrieved Green’s functions are still very small in this situation (about 10% or less), showing that neglecting these types of sources is allowed. When comparing these losses with Figure 14, it is visible that these amplitude losses are probably related to the fact that the volume source contributions (blue solid line in Figure 14) are neglected in Figure 16. In addition, Figure 17 shows the differences between the exact Green’s function and the Green’s function retrieved by considering only electromagnetic boundary source contributions (Je and Jm). As is visible, the electromagnetic boundary sources have a negligible contribution to the retrieved Green’s function in the positive time window GE,f(xB,xA,t). Hence, the amplitude losses visible in Figure 16 are indeed caused by ignoring the volume source contributions. However, the electromagnetic boundary sources do contribute to the Green’s function retrieval in the negative time window, that is, Gvs,Je(xA,xB,-t). The two red-dashed peaks at roughly t=-0.28 s and t=-0.38 s correspond to the spurious events B1 and B2, respectively. These spurious events result from the boundary and volume sources that are related to the edges of the modeling domain. They will remain present when considering only boundary sources or domain sources and will vanish when considering both. It is visible that the spurious events are not present in the exact case. Figure 17 clearly illustrates the contribution of the electromagnetic boundary sources in cancelling out the spurious events in the negative time window. As visible in Figure 15, the spurious events B1 and B2 are never stationary. Due to the bounded modeling domain, a contribution exists of sources at the edge of the modeling domain and that contribution needs to be compensated for by a source at that boundary surface.Figure 15
The obtained correlation gather of the domain integral for a three-layered medium bounded by a vacuum. The scale is taken as the logarithm of the absolute value of the amplitude. Summing this correlation gather panel yields the total contribution of the domain integral as shown in Figure14 by the blue volume line. Several events can be recognized.Figure 16
Comparison between the exact Green’s functionGE,f(xB,xA,t) and the Green’s function retrieved by considering only seismic boundary source contributions in the right-hand side of (52). In other words, the contributions of the domain sources and electromagnetic boundary sources are neglected. As can be seen, the amplitude errors between the exact and retrieved Green’s functions are acceptable (about 10% or less). Because only the Green’s function corresponding to the positive time window is considered here, no spurious events (which reside in the negative time window) are visible.Figure 17
Comparison between the exact Green’s functionGE,f(xB,xA,t)+Gvs,Je(xA,xB,-t) and the Green’s function retrieved by considering only electromagnetic boundary source contributions (Je and Jm) in the right-hand side of (52). In other words, the contributions of the domain sources and seismic boundary sources are neglected. As is visible, the electromagnetic boundary sources have a negligible contribution to the retrieved Green’s function in the positive time window GE,f(xB,xA,t). However, the electromagnetic boundary sources do contribute to the Green’s function retrieval in the negative time window, that is, Gvs,Je(xA,xB,-t). The two red-dashed peaks at roughly t=-0.28 s and t=-0.38 s correspond to the spurious events B1 and B2, respectively. It is visible that the spurious events are not present in the exact case. The contribution of the electromagnetic boundary sources in cancelling out the spurious events in the negative time window is clearly illustrated.Looking at Figure14, several other events are present as well. For seismoelectric exploration purposes, the events arriving at roughly t=0.14 s and t=0.23 s are of major interest. These represent the interface response (labelled 1) of the most shallow interface and its source-side ghost (labelled 1g). Similarly, the other two strong arrivals in the positive time window correspond to the interface response of the second, deeper interface (labelled 2) and its source-side ghost (labelled 2g). The schematic ray paths of these events are displayed in Figure 13. Especially worth noticing are two additional nonphysical events that reside in the negative time-window, labelled B1 and B2. As can be seen in Figure 14, the spurious event B1 generated by the boundary sources is equal but opposite in sign to the spurious event B1 from the volume sources. The same holds for spurious event B2. So, when retrieving the Green’s functions by using the complete right hand side of (52) these spurious events will vanish. However, when considering either boundary sources or domain sources, B1 and B2 will remain. The spurious events exist due to a correlation between a seismic and an electromagnetic wave event. Because the correlation implicitly subtracts the traveltime of the seismic event, which is relatively long compared with the traveltime of the electromagnetic wave (which arrives almost immediately), from the traveltime of the electromagnetic wave, the resulting spurious event resides in the negative time window.This is visible in Figure15. This figure represents the obtained correlation gather of the domain integral for a three-layered medium bounded by a vacuum. In other words, it represents the cross-correlation results for different source positions in the domain integral. The scale is taken as the logarithm of the amplitude. This, in order to be able to present the different events despite their large amplitude differences. Summing this correlation gather panel yields the total contribution of the domain integral as shown in Figure 14 by the blue volume line.As is visible, the correlation gather of this relatively simple 1D example already shows a great complexity of events. It contains lots of multiple arrivals and free-surface ghosts. Therefore, distinguishing all the different events is quite a task. Looking at the different events, some contributions are so-called non-stationary. That means that this contribution of a certain source position to a certain event shifts in time as a function of the source position [46]. For example, looking at the area in between the receiver positions xA and xB, all the non-horizontal events are non-stationary. However, outside the range enclosed by the two receivers, the contributions of the sources in the domain integral of the interferometric Green’s function representation are stationary. This combined with the slight amplitude losses visible in Figure 16 (about 10% or less), partly confirms both the analyses of Snieder [87] and Slob et al. [48]. They show that, for respectively the seismic interferometry and the electromagnetic interferometry, no spurious events will be created by neglecting the contribution of the domain integral in weakly dissipative media. Only the amplitudes of the retrieved events will be affected. Furthermore, the spurious events that are created in our modeling indeed only reside in the negative time window, as should be the case. Because Figure 16 only considers the Green’s function corresponding to the positive time window, no spurious events are visible.The numerical 1D SH-TE example presented here has shown that the presence of seismic sources only is sufficient to retrieve an accurate seismoelectric response. This means effectively that both seismic and electromagnetic signals are registered at different receivers (without the need of explicit electromagnetic sources) and that by cross-correlating these registered signals, the accurate seismoelectric Green’s function (less than 10% amplitude difference) is retrieved. In addition, it has been shown that the electromagnetic boundary source contribution to the Green’s function retrieval in the positive time window is negligible. However, the numerical example presented here is of course far from resembling a real Earth setting. Nevertheless, recent seismic interferometry studies performed on real data have shown that, for example, by using seismic noise sources (e.g. from microseisms),P-wave reflection responses can be correctly retrieved [88]. Here it is shown that for the seismoelectric case, the use of seismic sources only is sufficient to correctly retrieve the seismoelectric Green’s function response (for the coseismic field responses as well as the interface response fields). This seems promising for real applications of seismoelectric interferometry. We are currently investigating seismoelectric interferometry for both propagation modes (SH-TE and P-SV-TM) for 3D configurations.


## 5.1. Theory

Interferometry makes use of the cross-correlation of responses at different receivers in order to obtain the Green’s function of the field response between these stations. In other words, it is the deterministic response from one station to the other.Figure12 shows a possible seismoelectric interferometry setting. The cross-correlation of electric (Ex) and acoustic signals (vz) from sources located at the surface (Figure 12(a)) or in the bulk (Figure 12(c)) results in the direct electric response of an acoustic source (fz) generating a seismoelectric wave (Figure 12(b)). The known challenging problems in using seismo-electromagnetics as a geophysical exploration tool can potentially be addressed by applying interferometric Green’s function retrieval techniques to seismo-electromagnetic phenomena [46]. First of all, sources in “classical” seismoelectric surveys need to be strong. This is not always possible and therefore it is beneficial to be able to replace those strong sources by receivers.Figure 12
Schematic seismoelectric interferometry setting. Cross-correlation of electric (Ex) and acoustic signals (vz) from sources located at the surface (a) or in the bulk (c) results in the direct electric response of an acoustic source (fz) generating a seismoelectric wave (b).A second well-known problem in these conventional seismoelectric surveys is the very low signal-to-noise ratio. By doing interferometry, stacking inherently takes place with a possible improvement of the signal-to-noise ratio as a result. After deriving the system of equations for coupled seismic and electromagnetic waves in saturated porous media [29], the convolution-type reciprocity theorem and a power balance for seismoelectric waves was derived by Pride and Haartsen [67]. In 2003, this result was extended to a reciprocity theorem of the correlation-type for seismoelectric waves [84].Following Wapenaar and Fokkema [85], de Ridder et al. [46] showed that the 1D seismoelectric system for the SH-TE propagation mode can be captured in the following matrix-vector equation(45)iωÂû+Dxû=ŝ,
where matrix Â contains the space-dependent material parameters, û represents the field vector (in the space-frequency domain), matrix Dx contains the spatial differential operator ∂/∂x, ŝ denotes the source vector, and where iω arises due to Fourier transformation of the temporal derivative of a field. It is important to capture the 1D seismoelectric system in such a general diffusion, flow and wave equation, in order to employ the derived expressions for unified Green’s function retrieval by cross-correlation [47] and to finally end up with interferometric seismoelectric Green’s function representations.Next, considering the Fourier transform of an impulsive source acting at timet=0 s and at location x=xs, ŝ in (45) is replaced by Iδ(x-xs). As a consequence, the field vector û can be replaced by a Green’s matrix Ĝ(x,xs,ω). In this way, (45) is changed to(46)iωÂĜ(x,xs,ω)+DxĜ(x,xs,ω)=Iδ(x-xs),
where the Green’s matrix Ĝ(x,xs,ω) is given by(47)Ĝ(x,xs,ω)=(ĜE,Je(x,xs,ω)ĜE,Jm(x,xs,ω)ĜE,f(x,xs,ω)ĜE,hb(x,xs,ω)ĜH,Je(x,xs,ω)ĜH,Jm(x,xs,ω)ĜH,f(x,xs,ω)ĜH,hb(x,xs,ω)Ĝvs,Je(x,xs,ω)Ĝvs,Jm(x,xs,ω)Ĝvs,f(x,xs,ω)Ĝvs,hb(x,xs,ω)Ĝτb,Je(x,xs,ω)Ĝτb,Jm(x,xs,ω)Ĝτb,f(x,xs,ω)Ĝτb,hb(x,xs,ω)).The first superscript (vs) in Ĝvs,Je(x,xs,ω) denotes the type of response measured at location x, resulting from the type of impulsive source located at xs, which is denoted by the second superscript (Je).Starting from the general interferometric Green’s function representation (48) as derived by Wapenaar et al. [47](48)Ĝ(xB,xA,ω)+Ĝ†(xA,xB,ω)=-{Ĝ(xB,x,ω)NxĜ†(xA,x,ω)}|x1x2+∫D{Ĝ(xB,x,ω)[iω(Â-Â†)]Ĝ†(xA,x,ω)}dx,
where it is assumed that the two reciprocity states have the same medium parameters and where Nx represents the normal vector matrix containing the components of a normal vector n, arranged in the same way as the partial spatial derivative ∂/∂x in the matrix Dx(49)Nx=(0-100-100000010010),
de Ridder et al. [46] derived the following interferometric integral representation for one element of the seismoelectric SH-TE Green’s matrix in 1D (50), using seismoelectric reciprocity theorems(50){ĜE,f(xB,xA,ω)+Ĝvs,Je*(xA,xB,ω)}Ŝ={ĜE,Je(xB,x,ω)Ĝvs,Jm*(xA,x,ω)+ĜE,Jm(xB,x,ω)Ĝvs,Je*(xA,x,ω)-ĜE,f(xB,x,ω)Ĝvs,h*(xA,x,ω)-ĜE,h(xB,x,ω)Ĝvs,f*(xA,x,ω)}|x1x2Ŝ+2iω∫D{ĜE,Je(xB,x,ω)iI{ϵ̂}Ĝvs,Je*(xA,x,ω)+ĜE,f(xB,x,ω)R{ρfL̂0}Ĝvs,Je*(xA,x,ω)-ĜE,Je(xB,x,ω)R{ρfL̂0}Ĝvs,f*(xA,x,ω)+ĜE,f(xB,x,ω)iI{ρ̂c}Ĝvs,f*(xA,x,ω)}dxŜ.
To arrive at this form, they have chosen the {1,3} element of the 1D SH-TE seismoelectric Green’s matrix (47) and expanded (48) using this element. Here, Ŝ denotes the power spectrum of the emitted source signal and ρ̂c=ρb-iωk0(ρf)2/η.We can distinguish two terms in this integral representation. The first term on the right-hand side represents correlations of recorded responses of sources on the boundary of the domain of reciprocity, whereas the second term on the righthand-side represents correlations of recorded responses of sources throughout the reciprocity domain.As shown by de Ridder et al. [46], the following source-receiver reciprocity holds(51)Ĝvs,Je(xA,xB,ω)=-ĜE,f(xB,xA,ω).
Hence, the left-hand side of (50) can be rewritten as 2iℑ{ĜE,f(xB,xA,ω)}Ŝ. This signal will be antisymmetric around t=0 s in the time-domain.Looking at expression (50) in more detail it can be seen that the left-hand side, the electric field response registered at xB generated by an elastic force source located at xA, is obtained by cross-correlating the registered electric fields at xB with the registered particle velocities at xA, which are the result of four different types of boundary sources and two types of domain sources. The two types of domain sources, an electric current source and an elastic force source, are both weighted with two different medium parameters.Due to the fact that wave energy is dissipated during wave propagation, the domain sources are necessary to account for these losses. However, these sources are not likely to exist in reality or cannot be rewritten for practical applications and therefore we would like to be able to ignore their contributions.As is already shown in three examples by de Ridder et al. [46], it is indeed possible (under certain conditions) to obtain accurate Green’s functions from boundary sources only. The most complex situation considered by de Ridder et al. [46] was a medium consisting of two layers bounded by a vacuum. For this situation it was shown that the domain integral contribution could be neglected as long as the domain of reciprocity was chosen in such a way, that it included the heterogeneities (i.e., the interface between the two layers). Then, spurious events would only occur on one side of the symmetrized, retrieved Green’s function.In the following section, we will increase the complexity of the numerical configuration by adding an extra layer to the system, to investigate the Green’s function retrieval for a 1D, three-layered system bounded by a free-surface. In other words, we will look at the applicability of the interferometric seismoelectric Green’s function representation (50) when there are two interfaces located in the subsurface. We will consider a configuration where a medium B is sandwiched between two identical layers (medium A) with different medium parameters, as given in Table 3. We have chosen the medium parameters in such a way, that there is a very small seismic contrast between the layers; the porosity is the only contrasting seismic parameter between the layers. In this way, we are minimizing the dominant coseismic field response from the two subsurface interfaces and are able to focus mainly on the retrieval of the interface response field. The free-surface on the other hand, acts as a reflector for the seismic waves and therefore the coseismic fields related to this interface are still preserved. In this way, we are able to investigate separately the retrieval of both the coseismic field responses and the interface response fields. For field geometries, the amplitudes of the coseismic field responses related to the subsurface interfaces are often much higher than the interface response field amplitudes.Table 3
Overview of the relevant medium parameters for the 1D seismoelectric interferometry model.
PropertyUnitValue medium AValue medium BDimensionPorosityϕ0.40.2[—]Pore fluid densityρf1.0·1031.0·103[kg/m3]Solid densityρs2.7·1032.7·103[kg/m3]Shear modulus framework of grainsG9.0·1099.0·109[Pa]Pore fluid viscosityη1.0·10-31.0·10-3[kg/(m s)]Static permeabilityk01.3·10-121.6·10-12[m2]Static electrokinetic couplingℒ01.0·10-81.0·10-9[m2/(s V)]Tortuosityα∞3.03.0[—]Relative perm. of the (pore) fluidϵrf8080[—]Relative perm. of the solidϵrs44[—]Bulk electric conductivityσ01.0·10-11.0·10-1[S/m]


## 5.2. Results

We consider a three-layered 1D medium bounded by a vacuum half-space. The top and bottom layer consist of medium parameters belonging to medium A and the sandwiched layer has the properties of medium B (see Table3). The bottom layer is in fact a half-space. The whole three-layered system is bounded by a vacuum half-space in which only electromagnetic waves can propagate. The interface separating the subsurface from the vacuum is called the free-surface. According to Wapenaar and Fokkema [86], the free-surface acts as a mirror to both shear waves and electromagnetic waves (the latter due to the fact that a 1D geometry is considered here). Therefore, the sources on the domain boundary at the free surface can be neglected and also the contributions of the vacuum above the free-surface can be disregarded. The range of the domain integral contribution is from 0 to 650 m (see Figure 13 for an overview of the geometry). The receivers at xA and xB are located at x=100m and x=275m, respectively. The upper boundary is called x1 and the lower boundary x2, located at x=-150m and x=650m, respectively.Figure 13
The geometry of the 1D numerical experiment. PositionsxA and xB represent the receiver positions, denoted by the top-down triangles. The upper boundary is called x1 and the lower boundary x2, located at x=-150m and x=650m, respectively. The positions of the interfaces are visualized by the red bars. FS denotes the free-surface, whereas IF1 and IF2 correspond to the first and second subsurface interfaces, respectively. Furthermore, schematic ray paths of events in the seismoelectric Green’s function GE,f(xB,xA,t) are shown. The black arrows represent shear wave ray paths, whereas the green arrows depict electromagnetic ray paths. The labels correspond to the labelled events in Figures 14 and 15, where d denotes the direct event and dg its source-side ghost, which are both coseismic field responses, 1 represents the interface response of the most shallow interface and 1g its source-side ghost, 2 corresponds to the interface response of the second, deeper interface and event 2g represents again its source-side ghost.Figure14 shows the time-domain equivalent of the separated contributions of the domain integral and the boundary points to the retrieved Green’s functions. In other words, it shows the relative contributions of the two right-hand side terms in (50) to the retrieved Green’s functions. The positive time corresponds to the Green’s function GE,f(xB,xA,t), the electric field response registered at xB due to an impulsive seismic source located at xA. As is visible, the dominant contribution in the positive time window comes from just the boundary term. Therefore, it is shown that this Green’s function can be mainly reconstructed by using the boundary contribution only. In contrast, the negative time window contains strong domain integral contributions as well. The negative times correspond to the Green’s function Gvs,Je(xA,xB,-t), the particle velocity response measured at xA due to an impulsive electrical current source at xB. These strong volume source contributions correct the polarity of the single boundary term contribution.Figure 14
Separated contributions of the domain integral and the boundary points to the retrieved Green’s functions. In other words, it shows the relative contributions of the two right-hand side terms in (50) to the retrieved Green’s functions, for a three-layered medium bounded by a vacuum.Several events can be recognized in Figure14. The purely diffusive electromagnetic field is the first event to arrive, with its maximum at approximately t=0.58 ms (t=-0.58 ms for the time-reversed causal signal). The second arrival, at approximately t=83 ms corresponds to a direct coseismic shear wave event (labelled d). Its time-reversed causal equivalent arrives around t=-83 ms. The source-side ghost of this direct coseismic shear wave event (labelled dg) arrives at approximately t=0.18 s. The overlapping causal and time-reversed causal electromagnetic events at t=0 s are constructed completely by the sources in the domain integral. In contrast, the shear wave event is retrieved by mainly boundary source contributions. This makes sense, considering the sources of wave energy loss. As mentioned already, the electromagnetic event is, in the considered seismic frequency range, primarily a diffusive field. Therefore, volume sources are required to compensate for the wave energy loss. For the shear wave event, the amount of wave energy loss is relatively small. Hence, the need for volume source energy compensation is negligible. Considering this in terms of the interferometric seismoelectric integral representation (50),(52){ĜE,f(xB,xA,ω)+Ĝvs,Je*(xA,xB,ω)}Ŝ={ĜE,Je(xB,x,ω)Ĝvs,Jm*(xA,x,ω)+ĜE,Jm(xB,x,ω)Ĝvs,Je*(xA,x,ω)-ĜE,f(xB,x,ω)Ĝvs,h*(xA,x,ω)-ĜE,h(xB,x,ω)Ĝvs,f*(xA,x,ω)}|x1x2Ŝ+2iω∫D{ĜE,Je(xB,x,ω)iI{ϵ̂}Ĝvs,Je*(xA,x,ω)+ĜE,f(xB,x,ω)R{ρfL̂0}Ĝvs,Je*(xA,x,ω)-ĜE,Je(xB,x,ω)R{ρfL̂0}Ĝvs,f*(xA,x,ω)+ĜE,f(xB,x,ω)iI{ρ̂c}Ĝvs,f*(xA,x,ω)}dxŜ,
we can omit the volume source contributions(53)2iω∫D{ĜE,Je(xB,x,ω)iI{ϵ̂}Ĝvs,Je*(xA,x,ω)+ĜE,f(xB,x,ω)R{ρfL̂0}Ĝvs,Je*(xA,x,ω)-ĜE,Je(xB,x,ω)R{ρfL̂0}Ĝvs,f*(xA,x,ω)+ĜE,f(xB,x,ω)iI{ρ̂c}Ĝvs,f*(xA,x,ω)}dxŜ≈0.
The free-surface acts as a mirror to both shear waves and electromagnetic waves and therefore the sources on the domain boundary at the free-surface can be neglected and also the contributions of the vacuum medium above the free-surface can be disregarded. For (52), this effectively means that the contributions of the boundary sources at x1 can be omitted. The first two terms on the right-hand side of (52) represent the correlation products of the two fields generated by either electrical or magnetic current sources. The contribution of these two correlation products to the positive time window is very small, due to the fact that the causal fields registered at xB are electric fields which, without any wavetype conversion, only contribute significantly at t=0 (they arrive instantaneously). The contributions of the electric signals with positive seismic traveltimes are negligible, because these signals have encountered at least two wavetype conversions; this implies significant energy losses. Therefore, we can additionally neglect the contributions of the electromagnetic boundary sources (Je and Jm), ending up with the following reduced interferometric seismoelectric integral representation(54){ĜE,f(xB,xA,ω)+Ĝvs,Je*(xA,xB,ω)}Ŝ={-ĜE,f(xB,x2,ω)Ĝvs,h*(xA,x2,ω)-ĜE,h(xB,x2,ω)Ĝvs,f*(xA,x2,ω)}Ŝ.
Figure 16 shows the result of using this reduced interferometric seismoelectric integral representation. The figure displays a comparison between the exact Green’s function in the positive time window GE,f(xB,xA,t) and the Green’s function retrieved by using (54), that is, by considering only seismic boundary source contributions (f and h) at x2. As is visible, the amplitude errors between the exact and retrieved Green’s functions are still very small in this situation (about 10% or less), showing that neglecting these types of sources is allowed. When comparing these losses with Figure 14, it is visible that these amplitude losses are probably related to the fact that the volume source contributions (blue solid line in Figure 14) are neglected in Figure 16. In addition, Figure 17 shows the differences between the exact Green’s function and the Green’s function retrieved by considering only electromagnetic boundary source contributions (Je and Jm). As is visible, the electromagnetic boundary sources have a negligible contribution to the retrieved Green’s function in the positive time window GE,f(xB,xA,t). Hence, the amplitude losses visible in Figure 16 are indeed caused by ignoring the volume source contributions. However, the electromagnetic boundary sources do contribute to the Green’s function retrieval in the negative time window, that is, Gvs,Je(xA,xB,-t). The two red-dashed peaks at roughly t=-0.28 s and t=-0.38 s correspond to the spurious events B1 and B2, respectively. These spurious events result from the boundary and volume sources that are related to the edges of the modeling domain. They will remain present when considering only boundary sources or domain sources and will vanish when considering both. It is visible that the spurious events are not present in the exact case. Figure 17 clearly illustrates the contribution of the electromagnetic boundary sources in cancelling out the spurious events in the negative time window. As visible in Figure 15, the spurious events B1 and B2 are never stationary. Due to the bounded modeling domain, a contribution exists of sources at the edge of the modeling domain and that contribution needs to be compensated for by a source at that boundary surface.Figure 15
The obtained correlation gather of the domain integral for a three-layered medium bounded by a vacuum. The scale is taken as the logarithm of the absolute value of the amplitude. Summing this correlation gather panel yields the total contribution of the domain integral as shown in Figure14 by the blue volume line. Several events can be recognized.Figure 16
Comparison between the exact Green’s functionGE,f(xB,xA,t) and the Green’s function retrieved by considering only seismic boundary source contributions in the right-hand side of (52). In other words, the contributions of the domain sources and electromagnetic boundary sources are neglected. As can be seen, the amplitude errors between the exact and retrieved Green’s functions are acceptable (about 10% or less). Because only the Green’s function corresponding to the positive time window is considered here, no spurious events (which reside in the negative time window) are visible.Figure 17
Comparison between the exact Green’s functionGE,f(xB,xA,t)+Gvs,Je(xA,xB,-t) and the Green’s function retrieved by considering only electromagnetic boundary source contributions (Je and Jm) in the right-hand side of (52). In other words, the contributions of the domain sources and seismic boundary sources are neglected. As is visible, the electromagnetic boundary sources have a negligible contribution to the retrieved Green’s function in the positive time window GE,f(xB,xA,t). However, the electromagnetic boundary sources do contribute to the Green’s function retrieval in the negative time window, that is, Gvs,Je(xA,xB,-t). The two red-dashed peaks at roughly t=-0.28 s and t=-0.38 s correspond to the spurious events B1 and B2, respectively. It is visible that the spurious events are not present in the exact case. The contribution of the electromagnetic boundary sources in cancelling out the spurious events in the negative time window is clearly illustrated.Looking at Figure14, several other events are present as well. For seismoelectric exploration purposes, the events arriving at roughly t=0.14 s and t=0.23 s are of major interest. These represent the interface response (labelled 1) of the most shallow interface and its source-side ghost (labelled 1g). Similarly, the other two strong arrivals in the positive time window correspond to the interface response of the second, deeper interface (labelled 2) and its source-side ghost (labelled 2g). The schematic ray paths of these events are displayed in Figure 13. Especially worth noticing are two additional nonphysical events that reside in the negative time-window, labelled B1 and B2. As can be seen in Figure 14, the spurious event B1 generated by the boundary sources is equal but opposite in sign to the spurious event B1 from the volume sources. The same holds for spurious event B2. So, when retrieving the Green’s functions by using the complete right hand side of (52) these spurious events will vanish. However, when considering either boundary sources or domain sources, B1 and B2 will remain. The spurious events exist due to a correlation between a seismic and an electromagnetic wave event. Because the correlation implicitly subtracts the traveltime of the seismic event, which is relatively long compared with the traveltime of the electromagnetic wave (which arrives almost immediately), from the traveltime of the electromagnetic wave, the resulting spurious event resides in the negative time window.This is visible in Figure15. This figure represents the obtained correlation gather of the domain integral for a three-layered medium bounded by a vacuum. In other words, it represents the cross-correlation results for different source positions in the domain integral. The scale is taken as the logarithm of the amplitude. This, in order to be able to present the different events despite their large amplitude differences. Summing this correlation gather panel yields the total contribution of the domain integral as shown in Figure 14 by the blue volume line.As is visible, the correlation gather of this relatively simple 1D example already shows a great complexity of events. It contains lots of multiple arrivals and free-surface ghosts. Therefore, distinguishing all the different events is quite a task. Looking at the different events, some contributions are so-called non-stationary. That means that this contribution of a certain source position to a certain event shifts in time as a function of the source position [46]. For example, looking at the area in between the receiver positions xA and xB, all the non-horizontal events are non-stationary. However, outside the range enclosed by the two receivers, the contributions of the sources in the domain integral of the interferometric Green’s function representation are stationary. This combined with the slight amplitude losses visible in Figure 16 (about 10% or less), partly confirms both the analyses of Snieder [87] and Slob et al. [48]. They show that, for respectively the seismic interferometry and the electromagnetic interferometry, no spurious events will be created by neglecting the contribution of the domain integral in weakly dissipative media. Only the amplitudes of the retrieved events will be affected. Furthermore, the spurious events that are created in our modeling indeed only reside in the negative time window, as should be the case. Because Figure 16 only considers the Green’s function corresponding to the positive time window, no spurious events are visible.The numerical 1D SH-TE example presented here has shown that the presence of seismic sources only is sufficient to retrieve an accurate seismoelectric response. This means effectively that both seismic and electromagnetic signals are registered at different receivers (without the need of explicit electromagnetic sources) and that by cross-correlating these registered signals, the accurate seismoelectric Green’s function (less than 10% amplitude difference) is retrieved. In addition, it has been shown that the electromagnetic boundary source contribution to the Green’s function retrieval in the positive time window is negligible. However, the numerical example presented here is of course far from resembling a real Earth setting. Nevertheless, recent seismic interferometry studies performed on real data have shown that, for example, by using seismic noise sources (e.g. from microseisms),P-wave reflection responses can be correctly retrieved [88]. Here it is shown that for the seismoelectric case, the use of seismic sources only is sufficient to correctly retrieve the seismoelectric Green’s function response (for the coseismic field responses as well as the interface response fields). This seems promising for real applications of seismoelectric interferometry. We are currently investigating seismoelectric interferometry for both propagation modes (SH-TE and P-SV-TM) for 3D configurations.


## 6. Conclusion and Prospects

It was shown that the computed amplitude and phase for the dynamic permeability correlate well with the normalized measurements, whereas for the dynamic SP coefficient, only the normalized amplitude correlates well with the predictions of the theory. This difference could be due to a capacitor effect of the set-up. To prevent the capacitor effect, using insulating plates and electrodes perforated in them may be a solution. In addition, this difference could be related to a slight frequency-dependence of the bulk conductivity. Using independent impedance measurements of the sample could also improve the results. A full-waveform seismoelectric model in a layered geometry was obtained from the solution of a mechanical boundary value problem and the electric-solid ratio of the fastP-wave. The model was simplified by neglecting the interface response. The predictions of fast P-wave coseismic fields were compared against coseismic field measurements. Agreement was found in terms of travel time and waveform, while predicted amplitudes fell within the range of the measured amplitudes. Further modeling indicated that the (Biot) theory underestimates the measured seismic attenuation. The experimental results confirm the existing electrokinetic theory for the seismoelectric wave effect. Moreover, it was shown that coseismic fields can be modeled in a relatively simple way. The electrokinetic theory was subsequently adopted to study the applicability of seismoelectric interferometry. It was shown that the 1D interferometric seismoelectric SH-TE Green’s function representation retrieves accurate results for a three-layered 1D medium bounded by a vacuum. From the numerical results it can be concluded that seismic sources at a single boundary location are sufficient to extract the 1D electric field response generated by an impulsive seismic source in a layered model, both for the coseismic field responses and the interface response fields. In addition, it has been shown that the electromagnetic boundary source contribution to the Green’s function retrieval in the positive time window is negligible. However, the numerical example presented here is of course far from resembling a real Earth setting. Nevertheless, recent seismic interferometry studies performed on real data have shown that by using seismic noise sources, for example P-wave reflection responses can be correctly retrieved. Here it is shown that for the seismoelectric case, the use of seismic sources only is sufficient to correctly retrieve the seismoelectric Green’s function response (for the coseismic field responses as well as the interface response fields). This seems promising for real applications of seismoelectric interferometry.



---

*Source: 514242-2012-07-26.xml*