Forward and inverse modeling - The Seismoelectric Method - page 105

respectively, where d denotes the distance between S and
the interface and

simulation are reported in Table 4.1. Figure 4.1 shows
the two components (horizontal and vertical) of the
displacement of the solid skeleton and the two
components of the relative displacement w . The results
are shown at 15 Hz. To benchmark the numerical code,
we compare the numerical solution with the analytical
solution given by Dai et al. (1995). As shown in
Figure 4.2, both solutions are in excellent agreement.

λ S corresponds to the wavelength of the
seismic wave. The approximation in Equations (4.47)
and (4.48) are obtained by assuming that d >> λ S , yielding

r SE ≈ 2 r S . In Equation (4.48), it is further assumed that

the diffusion of the electromagnetic disturbances are
much faster than the propagation of the seismic energy,
a very good approximation as discussed in Chapters 2 and
3. Therefore, the lateral resolution of surface seismoelec-
tric data is poorer when compared with the lateral reso-
lution of surface seismic data. Fourie (2003) also showed
that for a horizontal interface, the first seismoelectric
Fresnel zone is nearly circular and centered beneath
the shot point S. For the seismic case, the first Fresnel
zone is an ellipse centered halfway between the shot
point and the observation point P.

4.2 Synthetic case study

Next, we show an application of the previously presented
model to the detection of a NAPL(oil)/water encroach-
ment front during the remediation of a NAPL-
contaminated aquifer by flooding the aquifer with water.
Such a remediation process can also be enhanced with
the use of surfactants (e.g., Mercier & Cohen, 1990; Pope
& Wade, 1995; Londergan et al., 2001). We will use the
subscript

4.1.6 Benchmark test of the code
We check the reliability of the finite-element formulation
we use in this topic through the use of a 2D benchmark
test. In this test, we simulate the fast and slow P-waves
associated with an explosive source in a homogeneous
porous material filled with a Newtonian fluid (water).
The dimensions of the 2D domain are 800 m × 800 m,
where the reference of the Cartesian coordinate system,
O(0, 0), is at the bottom left corner. To put seismic energy
into the system,weusea time-dependent source function,
F ( t ), that is a Ricker wavelet with a dominant frequency of
10 Hz, located at the source point S( x , y ) = (400 m, 400m).
The position of the receiver (observation point) is P( x , y )=
(200 m, 300 m). The four edges are absorbing boundaries
for which we use C-PML boundary conditions described
in Section 4.1.2. The material properties used for this

“ o ” to characterize the properties of oil.

4.2.1 Simulation of waterflooding
of a NAPL-contaminated aquifer
We use the following two steps to simulate the
remediation of a NAPL-contaminated aquifer.
Step 1. We use the approach developed in Karaoulis
et al . (2012) to generate a 2D heterogeneous aquifer in
terms of porosity and permeability (Figure 4.3). A random
field for the clay content was generated with the SGeMS
library (Stanford University). We used an isotropic
semivariogram to compute the clay content distribution
(Karaoulis et al . , 2012). The porosity and permeability
were then computed according to the petrophysical
model defined byRevil andCathles (1999). This heteroge-
neous aquifer is assumed to be initially saturated with
75%of lightmotor oil resulting fromanoil spill. The initial
water saturation of the aquifer is therefore s w =025,
which corresponds to the irreducible water saturation, s r .
The aquifer is located between two wells, Wells A and B.
Well B is located 250 m away fromWell A. The reference
position, O(

Table 4.1 Material properties used for the numerical benchmark.

Parameter Value Units

kg m − 3 ρ s 2650

kg m − 3 ρ f 1000

80,30), for the coordinate system is located
at the upper left corner of this domain.
Step 2. Waterflooding of this aquifer is simulated in
2.5D by injecting water inWell A (constant injection rate)
while removing the NAPL withWell B (constant pressure
condition). The computations are done in two-phase flow
conditions following the same equations as in Karaoulis
et al . (2012). The properties of the NAPL and water are

− K s 35 GPa

K fr 5 GPa

G 11 GPa

K f 2.25 GPa

1×10 − 3

η f Pa s

1×10 − 11 m 2

k 0

ϕ 0.30

—

1×10 − 2 Sm − 1 σ