The inverse problem for electroseismic conversion: Stable recovery of the conductivity and the electrokinetic mobility parameter

Pride (1994, Phys. Rev. B 50 15678-96) derived the governing model of 
electroseismic conversion, in which Maxwell's equations are coupled 
with Biot's equations through an electrokinetic mobility 
parameter. The inverse problem of electroseismic conversion was first 
studied by Chen and Yang (2013, Inverse Problem 29 115006). By 
following the construction of Complex Geometrical Optics (CGO) 
solutions to a matrix Schrodinger equation introduced by Ola and 
Somersalo (1996, SIAM J. Appl. Math. 56 No. 4 1129-1145), we analyze 
the recovering of conductivity, permittivity and the 
electrokinetic mobility parameter in Maxwell's equations with internal 
measurements, while allowing the magnetic permeability $\mu$ to be a 
variable function. We show that knowledge of two internal data sets 
associated with well-chosen boundary electrical sources uniquely 
determines these parameters. Moreover, a Lipschitz-type stability is 
obtained based on the same set.


1.
Introduction. In fluid-saturated porous media, an electrical double layer (EDL) is formed at the interface (the pore boundaries) of the fluid and solid rock. The fluid side of the EDL is charged with positive ions and the solid side with negative electrons. When electric or magnetic fields impinge on the EDL, the electrokinetic phenomenon causes movement of the fluid relative to rock frame and thus emits seismic waves, which can be remotely detected. This effect is named electroseismic conversion. Conversely, a seismic wave can cause the separation of charges and therefore generate electromagnetic fields. At the intersection of two formation layers, the singularities in electrical, hydraulic and mechanical properties will lead to a discontinuity in the induced electromagnetic fields, and emit electro-magnetic waves, which can be remotely observed. This effect is named seismoelectric conversion. In fact, these two ways of conversion always happen simultaneously.
The wave propagation in fluid-saturated porous media was studied by Biot [3,4]. In Biot's theory, in addition to the conventional compressional and shear waves, a compressional slow wave appears. The first experimental observation of this slow wave was obtained by Plona [14].
Some basic properties of the coupling effect were studied by Pride and Haartsen in [16].
We note that the term L(−∇p + ω 2 ρ f u) in (2) and the term LE in (4) model the seismoelectric effect and electroseismic effect, respectively. Under the assumption that the coupling is weak, depending on the source type, we can focus on one way coupling and ignore the second coupling which is weaker as it is caused by induced waves. For instance, in the case when a seismic source is employed, seismoelectric conversion will dominate and the governing equations are given as above but with the term LE removed from (4). Starting with electrical sources, electroseismic conversion is dominant and seismoelectric conversion is negligible, in which case the governing model is as above but with the term L(−∇p + ω 2 ρ f u) removed in (2).
In 2005, White [21] established a forward model for the electroseismic method by solving Maxwell's equations and the elastic wave equation sequentially, while the initial amplitude of elastic waves is calculated according to Pride's equations with the high-frequency asymptotic theory.
In 2007, Thompson et al. at ExxonMobil [19] presented three field experiments of electroseismic conversion in Texas and Canada. Two tests were applied to shallow gas sands. By using specially designed electrical source waveforms at dominant frequencies of 8Hz, 18Hz and 25Hz, these tests showed that electroseismic conversion at gas sands up to 1000 m deep can be detected with geophones placed on the surface of the earth. The tests also succeeded in identifying several gas sands and delineating hydrocarbon accumulations. The third test detected the second-order conversion from a carbonate oil reservoir greater than 1500 m deep. Their well tests indicated clear electroseismic conversion in the formation and were proved useful in diagnosing well properties, including cement integrity and near-well formation. Schakel's [17] designed laboratory experiments to demonstrated electroseismic conversion in the frequency range 50-150 kHz and observed that electroseismic amplitudes as a function of frequency maximize at 70 kHz.
The inverse problem of electroseismic conversion consists of two steps. The first step is the inverse source problem of Biot's equations, which recovers the coupling term Σ := LE in (4) from boundary measurements of seismic waves. The second step is the inversion of Maxwell's equations, which recovers the conductivity σ and the coupling coefficient L from the internal data Σ. The critical advantage of exploiting electroseismic conversion is that the first step provides certain internal data Σ for the inversion of the (time-harmonic) Maxwell's equations and, as we will show, lead to well-posedness of this problem.
By substitution, we obtain the curl-curl form of Maxwell's equations The boundary source is expressed in terms of the boundary tangential components of the electric field, where tE is the tangent component of E. We also assume that the internal data of the form Σ = LE ∈ H s−1 (Ω) is given.
The inverse problem we study in this paper is that, with a number of properly chosen boundary sources of the electrical fields G, can we recover the coupling coefficient L and the complex parameter γ from the internal data Σ?
The inverse problem of Maxwell's equations with boundary Cauchy data or impedance map are studied in [5,12,13], and that with different types of internal data are studied in [2,7,9]. One reconstruction methodology of the inverse problem with internal data, inspired by Bal and Uhlmann's study of photo-acoustic tomography (PAT) [1], primarily consists of converting the governing equation to a transport equation and constructing the Complex Geometrical Optics (CGO) solutions to the governing equations. The explicitly constructed CGO solutions then lead to the unique and stable solvability of the transport equation. The same methodology was followed by Chen and Yang [7] when they proved the uniqueness and stability of the reconstruction of the coupling coefficient L and conductivity σ in the case of a constant magnetic permeability µ. In this paper, we extend their results to general cases to allow variable µ which approaches the physics while providing an insight in the source design for experiments. We follow a different method to construct CGO solutions.
The construction of CGO solutions for a variable µ employs the idea of converting Maxwell's equations to a matrix Schrödinger equation, to which CGO solutions can be constructed. The matrix Schrödinger equation formulation was first developed by Ola and Somersalo [13] to study the inverse boundary value problem in electromagnetics; the relevant analysis was simplified in [12].
We remark that the solvability of Maxwell's equations by CGO solutions is validated by the explicit construction. We also remark that, for any preselected frequency ω, the construction argument requires µ > 0 and γ = + iσ/ω = 0. Thus, the method of CGO solutions works fine in the diffusive regime when ω is close to zero and hence γ can be approximated by its second term even though is not necessarily small. This paper is organized as follows: in Section 2, we introduce the matrix Schrödinger equation, of which the CGO solutions are constructed in Section 3. Our main theorems are stated and proven in Section 4. We derive the transport equation and prove the uniqueness result in Section 4.1. We prove the stability result in Sections 4.2. In Section 5, we analyze the temporal behavior of our CGO solutions used in the stable recovery of L and γ.
The CGO solutions to (11) were constructed by Colton [8] and extended to high order Sobolev spaces by Chen and Yang [7]. However, their construction fails for non-constant µ.
Here, we study the case when µ is variable. We follow the argument in [12,5] to convert the first-order system of Maxwell's equations to a matrix Schrödinger equation and construct the corresponding CGO solutions.
We write (13) as a 8 × 8 matrix system, We define which satisfy We also define Here A, B are 4 × 4 diagonal complex-valued matrices. We now introduce scalar fields Φ and Ψ, and X = (Φ, H, Ψ, E) t , so that the matrix system in (14) attains the form Solutions to this system with Φ = Ψ = 0 correspond to solutions of the original Maxwell system. By changing variables so that A straightforward calculation leads to the following ). One has Here, the matrix potentials,Q µ,γ andQ µ,γ , are given bỹ . We can extend µ and γ to R 3 so that for some nonzero constants µ 0 and γ 0 , µ − µ 0 and γ − γ 0 are compactly supported. We also assume that Z is the solution to (17) (P ∓ (D) − W t µ,γ )Z = Y. Lemma 2.1 then implies that Z solves the matrix Schrödinger equation where k = ω(µ 0 γ 0 ) 1/2 and Q µ,γ = k 2 I 8 +Q µ,γ . It follows that Q µ,γ is compactly supported.
In the following section, we will construct solutions to (18) and thus solutions to Maxwell's equations according to (15) and (17). [20] constructed CGO solutions to the scalar Schrödinger equation. Ola and Somersalo [12] followed the Sylvester-Uhlmann method to construct CGO solutions to the matrix Schrödinger equation given in (18). The CGO solutions by Ola and Somersalo are in a weighted L 2 space. In the present work, we apply the Sylvester-Uhlmann method to construct CGO solutions in higher order Sobolev spaces.

Complex Geometrical Optics (CGO) solutions. Sylvester and Uhlmann
We first introduce some notations. Let the space L 2 δ for δ ∈ R be the completion of C ∞ 0 (R 3 ) with respect to the norm · L 2 δ defined by We also define the space H s δ for s > 0 as the completion of C ∞ 0 (R 3 ) with respect to the norm · H s δ defined by Here (I − ∆) s 2 u is defined as the inverse Fourier transform of ξ sû (ξ), whereû(ξ) is the Fourier transform of u(x). When δ = 0, this is the standard Sobolev space It is proven in [20] that for |ζ| ≥ c > 0 and v ∈ L 2 δ+1 with −1 < δ < 0, equation (19) ( for some constant C = C(δ, c). We note that (−∆+2ζ ·D) and (I −∆) s are constant coefficient operators and hence commute. We deduce that for v ∈ H s δ+1 , for s ≥ 0, One defines the integral operator G ζ : where F −1 is the inverse Fourier transform. Clearly, for |ζ| ≥ c > 0, G ζ is the inverse operator of (−∆ + 2ζ · D) and G ζ is bounded by for some constant C = C(δ, c).
As in [5], we choose ζ ∈ C 3 such that ζ · ζ = k 2 and |ζ| large. Compared to [12,5], while requiring more regularity in µ, γ, the following proposition extends the construction of CGO solutions to H s δ (Ω), which is the restriction of H s δ (R 3 ) to Ω. Proposition 1. Let s > 3 2 and −1 < δ < 0. Assume that µ, γ ∈ H s+2 (Ω) and µ, γ are constant outside a compact set. There exists a CGO solution in H s δ (Ω) of the form to equation (18), such that Z 0 ∈ C 3 is a constant vector and is an algebra. By the assumption, Q µ,γ ∈ H s (Ω) is compactly supported and thus multiplication by Q µ,γ is a bounded operator mapping H s δ (R 3 ) to H s δ+1 (R 3 ). In particular, Q µ,γ Z 0 ∈ H s δ+1 (R 3 ) for any constant Z 0 . The estimate in (22) implies that I + G ζ (Q µ,γ ·) is well defined and is an invertible mapping for |ζ| sufficiently large. We can then define (25) The operator Z r also satisfies (23) and (24). Thus e iζ·x (Z 0 + Z r ) ∈ H s δ (R 3 ) is a solution to (18). Using (25) we conclude that With CGO solutions to the matrix Schrödinger equation, we can now construct solutions to the original Maxwell equations. We define Y by (17) and thus Y satisfies (16). Note that the solution Y can be written as Here, Y 0 is a constant vector depending on the choice of ζ and Z 0 . We denote the components of Y by We recall that (16) is equivalent to the original Maxwell equations only if Y Φ = Y Ψ = 0. This condition can be satisfied with properly chosen Z 0 according to the following lemma, which is cited from [5] without proof.

Proposition 2. The Maxwell equations have a solution of the form
The proof of the above proposition follows directly from (15) and Lemma 3.1. An explicit choice of Z 0 satisfying Lemma 3.1 is given by where a, b are any vectors in C 3 . It follows that We choose b ∈ C 3 such that |b| = O(|ζ|), and define a CGO solution as . Compared to the CGO solutions in [7], our new construction applies to variable µ, while the required regularity of parameters µ, γ is one order higher. Here and below, we use a check sign, for exampleĚ, to indicate corresponding field variables computed from the CGO solutions. 4. Inversion of Maxwell's equations with internal data. In this section, we state and prove the uniqueness and stability results of the inverse problem of Maxwell's equations with internal data. In addition to Sobolev spaces H s (Ω), we also work on continuous function spaces C d (Ω). Here s and d are related by s = 3 2 + d + 2 + ι for some small ι > 0. We assume that µ ∈ H Sources generating the data: Our proof makes use of CGO solutions, constructed in the previous section, of the form (31). For a fixed frequency, the electric boundary sources are controlled to be close to the traces of the CGO solutions. We cannot consider arbitrary sources. When h is small, (32) indicates that the dominant term, e iζ·xĚ 0 , of a CGO solution is characterized by parameters ζ, a, b and the coefficient γ. We assume γ is known on the domain boundary. Therefore, for a particular choice of ζ, a and b and the given γ on the boundary of the domain, the restriction of e iζ·xĚ 0 to the boundary gives an explicit distribution of the controlled source in an experiment. In our inversion method, we construct two CGO solutions, according to which we define two tangential electric field boundary values, denoted by (G 1 , G 2 ).
Here and in the following, we shall abuse the notation and use C d (Ω) to denote either a set of complex-valued functions or a set of vector-valued functions the elements of which have up to d-th order continuous derivatives. The function space (H d+3+ι (∂Ω)) 2 is an abbreviation of the product space H d+3+ι (∂Ω)×H d+3+ι (∂Ω)).
The proof of the uniqueness theorem establishes an explicit reconstruction. Following this reconstruction we also prove a Lipschitz stability result. To consider the stability of the reconstruction, we restrict the proof to a subset of Ω. Let ζ 0 be a constant unit vector close to ζ 1 /|ζ 1 |, where ζ 1 occurs in the CGO solution. An explicit choice of ζ 0 and ζ 1 is given in (40) and (41). Define Ω 1 to be the subset of Ω obtained by removing a neighborhood of every point x 0 ∈ ∂Ω such that n(x 0 ) · ζ 0 = 0, where n(x 0 ) is the outward normal of ∂Ω at x 0 .
There is a non-empty open set of G in (H d+3+ι (∂Ω)) 2 , defined as a neighborhood of the trace of CGO solutions given by (31) and (40), for small ι > 0 such that for some constant C.
Similar to the proof in [7], if we have at least 6 sets of complex measurements, we can prove the same stability results with the subset requirement removed and obtain (35) onΩ. The proof relies on the unique solvability of a transport equation, the vector field of which can be controlled by properly constructed CGO solutions of Maxwell's equations.
Remark. In (40) we construct CGO solutions with desired asymptotic properties when |ζ| → ∞, while |ζ| can be chosen large and independent of the preselected frequency ω, see (40). Thereby, for any fixed frequency ω, we can still choose ζ so that Ě r H s−1 δ (Ω) is small. This also validates the application of our results to the low frequency diffusive EM fields.
We recall that E = Σ/L. By some calculations, we get Here χ(x) is any smooth nonzero complex-valued function.
By the method of characteristics, the transport equation (37) has a unique solution only if the integral curves of φ connect every point in Ω to a point on ∂Ω. We will prove, that with properly chosen CGO solutions, the integral curves of φ will be close to straight lines, and thus connect every internal point to two boundary points.
Let ω be fixed. Recall that ζ · ζ = k 2 = ω 2 γ 0 µ 0 and large |ζ| leads to small residual term Z r in CGO solutions due to Proposition 1. We make specific choices for ζ, a, b in (31): with h > 0 a free parameter. The choice of h is independent of ω. Moreover, |ζ| ≈ 1/h for small h independent of ω; precisely, By choosing any a j ∈ C 3 and b 1 = b 2 = (0, 0, 1/h), we can verify that CGO solutionš E 1 andĚ 2 defined in (31) satisfy We recall that the check sign, for exampleĚ j , indicates fields corresponding with the CGO solutions. We denoteΣ j = LĚ j = ϑe iζ·x (η j + R j ) with ϑ = Lk γ 1/2 , η j = γ 1/2Ě j 0 /k and R = γ 1/2Ě R /k. Then η 1 · η 2 = 1 + O(h). We now evaluateφ. We find that Similarly, on a bounded domain, Moreover, We also find that We let s = 3 2 + d + 2 + ι for d > 0 and small ι > 0. The Sobolev embedding theorem implies that If we choose b 1 = b 2 = (0, 0, h ε−1 ) and χ(x) = − h 1−2ε 4 e −i(ζ1+ζ2)·x with 0 < ε < 1/2, the same calculations as above give So far we proved that if the electric fields are exactly given by the CGO solutions as we constructed, the vector fieldφ in (37) can be approximated by (43), which implies that every integral curve ofφ is almost a straight line and thus connects every internal point to two boundary points. Next, we prove that we can perturb the particular CGO solutions so that an estimate similar to (43) still holds. To follow the dependencies of vector fields on boundary conditions, we introduce a regularity theorem for Maxwell's equations. Let tE be the boundary condition for the tangential component of E. The following function spaces were introduced in It is clear that t(H l Div (Ω)) = T H l−1/2 Div (∂Ω). We also observe that (46) E H l (Ω) ≤ E H l Div (Ω) and G T H l Div (∂Ω) ≤ G H l+1 (∂Ω) . Proposition 3 ( [11]). Let , µ ∈ C l , l > 2, be positive functions. There is a discrete subset Ξ ⊂ C, such that if ω / ∈ Ξ, then one has a solution E ∈ H l Div to (7) given any tangential boundary condition G ∈ T H l−1/2 Div (∂Ω). The solution satisfies with C independent of G.
We now estimate vector field φ in (37) when the boundary value of the electric field is given by a perturbation of the trace of a CGO solution.
Let s = 3 2 + d + 2 + ι, for small ι > 0. Then this proposition follows from Proposition 3.6 in [7] directly. To be self-contained in this section, we summarize the proof here.
Proof. According to the Sobolev embedding theorem, Proposition 3 and (46), we have that where various constants are all named "C". Hence, We now define boundary conditions G j ∈ H d+3+ι (∂Ω), j = 1, 2, such that for some ε > 0 sufficiently small. Let E j be the solution to the Maxwell equations (9) with tE j = G j . By (49), we thus have for some positive constant C. We introduce the complex-valued internal data, Σ j = LE j and conclude that We obtain the estimate (cf. (38)) where χ(x) is defined above (43). Therefore, (48) follows from (44) and (53).
We recall that M is the parameter space of (L, γ) defined in (34). We now prove Theorem 4.1.
Proof of Theorem 4.1. Let d ≥ 3. By Proposition 4, we choose the set of boundary conditions for electrical fields to be a neighborhood of (Ǧ j ) = (tĚ j ) in H d+3+ι (∂Ω)) 2 . By assuming that the measurements coincide, that is, Σ =Σ, we have that φ =φ and ψ =ψ by (38)-(39). Thus, L andL solve the same transport equation (37) while L =L = Σ/G on ∂Ω. As φ satisfies (48), we deduce that L =L since integral curves of φ map any x ∈ Ω to two boundary points. More precisely, consider the flow θ x (t) associated with the imaginary part of φ, that is, θ x (t) is the solution tȯ By the Picard-Lindelöf theorem, (54) admits a unique solution since φ is of class C 1 (Ω). Also by (48), for all x ∈ Ω, there exist x ± (x) ∈ ∂Ω and t ± (x) > 0 such that By the method of characteristics, the solution L to the transport equation (37) is given by where L 0 := L| ∂Ω is the restriction of L on the boundary. The solutionL is given by the same formula since θ x (t) =θ x (t). This implies that L =L and thus E j =Ẽ j = Σ j /L =Σ j /L, j = 1, 2. By the choice of illuminations, we also have |E j | = 0 due to (51) and |Ě j | = 0. Finally, by substituting E j andẼ j into (9), we can solve for γ andγ by We then conclude that γ =γ in C d−1 (Ω) for d ≥ 1.

Stability result.
The proof of the stability theorem follows the argument of [7] with the estimate of the vector field replaced by (48).
Proposition 5. Let d ≥ 1. Let L andL be solutions to (37) corresponding to coefficients (φ, ψ) and (φ,ψ), respectively, where (48) holds for both φ andφ. Let L 0 = L| ∂Ω ,L 0 =L| ∂Ω , and L 0 ,L 0 ∈ C d−1 (∂Ω). We also assume that h is sufficiently small and that Ω is convex. Then there is a constant C such that upon restricting to Ω 1 , The choice of Ω 1 depends on the proof of the above proposition. Readers are referred to [7] for the details. Now we can prove the main stability theorem.
Proof of Theorem 4.2. From (38) and (39), it is straightforward to check that where C > 0 is a positive constant. The first part of (35) then follows directly from Proposition 5. To estimate the difference between γ andγ, we notice that Because L andL are non-vanishing, by the stability result for L we obtain By choosing the boundary values close to the boundary conditions of CGO solutions, (50) and (51) imply that E j is non-vanishing since the CGO solutions are nonvanishing. We recall that γ andγ are computed by (55). By taking the difference and using (57) we find that This completes the proof.
The required regularity of the different quantities is summarized in Table 1, where s = 3 2 + d + 2 + ι.
The measured internal data are then given by Σ j 1 , Σ j 2 . Proposition 3 shows that the vector field defined by (38) satisfies Following the similar analysis in the proof of Theorem 4.2, we can prove (35) on Ω when more internal measurements are available.
Remark. Although the choice of b j such that ζ j ∧ b j = 0 could simplify the above argument, it will destroy the proof in Section 4.1. Also, in the second choice of (a j , b j ), we would prefer large ε to get better approximations in (60). However, the proof in Section 4.1 requires larger 1 − 2ε in (45), that is, smaller ε, to obtain a better stability estimate. We need to choose ε to balance these two conditions.
Suppose α and ω are bounded and h is a fixed small number. Upon taking an inverse Fourier transform of (59) and (60), the time behaviors of our CGO solutions are of the form of a multiplier, f (t), which is concentrated at zero time. To produce the field experiments, f (t) is convolved with a waveform, the amplitude spectrum of which is centered at a particular frequency.
Our forward model is the same as that of White [21], that is to excite electroseismic conversion with electrical sources and detect the seismic responses. While White focused on the forward problem with the high frequency asymptotic analysis, we consider the inverse problem with temporal δ-type sources.