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.
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 Figure 1). 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].
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 Figure 2. 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 . 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].
(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 dynamic SP 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 dynamic SP 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 an time dependence, yields the following linearized set of governing equations where , , are the Biot Gassmann constants [62], the shear modulus, is the fluid density, is the solid density, is the tortuosity, is the fluid viscosity, is the (static) permeability, is the electric field, and is the dynamic electrokinetic coupling [29] where is a characteristic pore size parameter and 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]: . The characteristic (or rollover) frequency is defined as [59]. The Debye length is denoted by (see, e.g., [29]) and represents the static electrokinetic coupling for a porous medium where is the vacuum permittivity, 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 becomes smaller than the Debye length. However, due to the fact that the Debye length 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 the viscous correction factor is defined by Johnson et al. [59] as The coefficients , , and are the so-called generalized effective density functions [65] Considering the definitions for and , (2) can be written as where is the relative displacement. Pride [29] developed the following equation coupling the streaming and the conduction currents where 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 or where represents the bulk electric conductivity and the pore-fluid conductivity. Closely related to the dynamic electrokinetic coupling (3) is the dynamic SP coefficient, defined as [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 Eliminating from (8) and (9), we obtain The Maxwell relation for the magnetic field is given by Ampère’s Circuit Law with the magnetic field and the electric permittivity for a fluid-saturated porous medium where is the solid relative permittivity. Faraday’s induction law states that with the magnetic permeability. Substitution of (14) in (13) results in Substituting the cross-product of Faraday’s law (16) [29, 67] into (17) yields where is the effective electric permittivity [67] of the porous continuum Here is a term accounting for the energy losses. The electrokinetic effect manifests itself in as an energy gain that is quadratic in (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 , (mechanical part of the equations), and electric fields (the electromagnetic part).
3. Experimental Validation of the Dynamic Coupling Coefficients
We experimentally validate and . 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 kg/m3, a viscosity of Pa s obtained from [68], and a measured pore fluid conductivity of S/m), whereafter the setup is left until equilibrium of the salt solution is reached.
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 Figures 5 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.
The measurements of the normalized dynamic SP coefficient (normalized to the measured value at 11 Hz, where 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 . Using (9) for a conservative (irrotational) electric field (with the streaming potential difference) in a setup where the electric current density is equal to , we obtain with being the dynamic SP coefficient. The dynamic SP 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.
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 slow -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 from (1) and (18), we obtain two modified momentum equations for the fields and where complex effective densities , , and , containing the electrokinetic coupling factor , are defined as follows Employing Helmholtz decomposition for the fields and leads to Substituting expressions (23) into (21) yields where . For the longitudinal waves, associated with potentials and , the first terms in square brackets of (24) are set equal to zero from which we obtain where we used that , and . Applying a spatial Fourier transformation and recasting (25) into an eigenvalue problem lead towhere 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 () and slow () longitudinal waves , , where The slowness yields the wave mode speed and intrinsic attenuation (see, e.g., [67]). For the transversal waves, associated with potentials and , the second term in square brackets of (24) are set equal to the zero vector which giveswhere a spatial Fourier transformation is applied. Nontrivial solutions for 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 () and seismic shear () waves. The dispersion relations are given in (27) for where and where we used that and . Note that 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 for . By writing the vector potentials as = and , for , in (29) we obtain for the transversal fluid-solid ratios 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 where we note that the Helmholtz decomposition of the electric field is , and . 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
We now model coseismic electric potentials generated within a porous medium due to a fast -wave, using its electric-solid ratio , 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 -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]) where is the distance to the source, is the angle of incidence, is the amplitude spectrum, and is the acoustic fluid wavenumber, where the fluid -wave slowness is given by . The directivity function , which characterizes the radiation pattern of the source, is given by Here, is the Bessel function of the first kind and first order and 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 at receiver position for , , , where and are the radial and vertical components of , respectively, and is the vertical component of the fast -wave wavenumber. The fast -wave wavenumber is , where the fast -wave slowness is given by the dispersion relations (27)-(28). Note that the factor 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 relates the incident acoustic wavefield to the transmitted fast -wave signal. The transmitted signal generates a coseismic potential at . It also reflects [] at , and travels back to the receiver position, where it generates a second coseismic potential. The transmission coefficient is derived from substituting plane wave expressions into the following (open-pore) boundary conditions [76] with subscript denoting the -component of the vectors and where 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) for . Hence we consider an incident () acoustic wave that reflects () and transmits as -, -, and vertical shear ()-waves. Displacement fields are obtained from these expressions as follows Fluid pressure is related to fluid displacement by , with . For the poroelastic medium, solid displacement and pore-fluid displacement are obtained as follows 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 so that we arrive at the following linear system of equations where the elements of matrix are given in the appendix. By solving (44) and (A.3) we obtain and , respectively (see appendix).
For the geometry of Figure 9, where a source is located at cm, and where the receiver is located at cm, we numerically evaluate the integral of (38). An experimentally recorded 500 kHz single sine pressure waveform is used for the amplitude spectrum . The incident pressure is related to the mechanical displacement potential in the denominator of (see (34)) by the factor , which arises from the relation . 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 -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 -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.
(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 cm recorded several pulses. The first (IRF1) is the interface response generated at the front () 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 -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 -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 -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 -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 Figure 10(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 -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 Figure 10(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 (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].
(a)
(b)
(c)
(d)
It is possible to model all interface responses and coseismic effects of Figure 10(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 and . Therefore, in the above, we only adopted Biot’s poroelastic theory to solve for and used the electric-solid ratio 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.
Figure 12 shows a possible seismoelectric interferometry setting. The cross-correlation of electric () and acoustic signals () 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 () 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.
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 where matrix contains the space-dependent material parameters, represents the field vector (in the space-frequency domain), matrix contains the spatial differential operator , denotes the source vector, and where 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 time s and at location , in (45) is replaced by . As a consequence, the field vector can be replaced by a Green’s matrix . In this way, (45) is changed to where the Green’s matrix is given byThe first superscript () in denotes the type of response measured at location , resulting from the type of impulsive source located at , which is denoted by the second superscript ().
Starting from the general interferometric Green’s function representation (48) as derived by Wapenaar et al. [47] where it is assumed that the two reciprocity states have the same medium parameters and where represents the normal vector matrix containing the components of a normal vector , arranged in the same way as the partial spatial derivative in the matrix 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 To arrive at this form, they have chosen the element of the 1D SH-TE seismoelectric Green’s matrix (47) and expanded (48) using this element. Here, denotes the power spectrum of the emitted source signal and .
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 Hence, the left-hand side of (50) can be rewritten as . This signal will be antisymmetric around 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 generated by an elastic force source located at , is obtained by cross-correlating the registered electric fields at with the registered particle velocities at , 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.
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 Table 3). 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 and are located at m and m, respectively. The upper boundary is called and the lower boundary , located at m and m, respectively.
Figure 14 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 , the electric field response registered at due to an impulsive seismic source located at . 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 , the particle velocity response measured at due to an impulsive electrical current source at . These strong volume source contributions correct the polarity of the single boundary term contribution.
Several events can be recognized in Figure 14. The purely diffusive electromagnetic field is the first event to arrive, with its maximum at approximately ms ( ms for the time-reversed causal signal). The second arrival, at approximately ms corresponds to a direct coseismic shear wave event (labelled ). Its time-reversed causal equivalent arrives around ms. The source-side ghost of this direct coseismic shear wave event (labelled ) arrives at approximately s. The overlapping causal and time-reversed causal electromagnetic events at 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), we can omit the volume source contributions 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 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 are electric fields which, without any wavetype conversion, only contribute significantly at (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 ( and ), ending up with the following reduced interferometric seismoelectric integral representation 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 and the Green’s function retrieved by using (54), that is, by considering only seismic boundary source contributions ( and ) at . 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 ( and ). As is visible, the electromagnetic boundary sources have a negligible contribution to the retrieved Green’s function in the positive time window . 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, . The two red-dashed peaks at roughly s and s correspond to the spurious events and , 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 and 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.
Looking at Figure 14, several other events are present as well. For seismoelectric exploration purposes, the events arriving at roughly s and s are of major interest. These represent the interface response (labelled 1) of the most shallow interface and its source-side ghost (labelled ). 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 ). 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 and . As can be seen in Figure 14, the spurious event generated by the boundary sources is equal but opposite in sign to the spurious event from the volume sources. The same holds for spurious event . 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, and 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 Figure 15. 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 and , 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), -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 fast -wave. The model was simplified by neglecting the interface response. The predictions of fast -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 -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.
Appendix
Substituting plane wave expressions into the poroelastic boundary conditions (39) for an incident acoustic wave from the fluid which impinges on a fluid/poroelastic-medium boundary leads to the following linear system of equations where the elements of matrix are and , and where , for .
Substituting plane wave expressions into the poroelastic boundary conditions (39) for an incident fast -wave from a poroelastic-medium which impinges on a fluid/poroelastic-medium boundary leads to the following linear system of equationswhere the elements of matrix are and . The reflection and transmission coefficients are defined as The potentials occur in the following plane wave expressions for is , or , . The subscript and refer to the incident and reflected wave, respectively. These potentials are related to the displacement fields as follows Fluid pressure is related to fluid displacement by . Pore-fluid pressure and intergranular stresses are obtained following the basic equations described in [60–64].
Acknowledgments
The research was performed at Delft University of Technology. The research was partly funded as Shell-FOM (Fundamental Research on Matter) projects within the research programs “The physics of fluids and sound propagation” and “Innovative physics for oil and gas.” The authors would like to thank Kees Wapenaar for useful comments, as well as the editor and two anonymous reviewers for constructive remarks and suggestions to improve this journal paper.