Finite element analysis of shallow-watermarine controlled-source electromagnetic signals for hydrocarbon exploration
Mark E. Everett
Department of Geology and Geophysics
TexasA&MUniversity, College StationTX77843USA
Introduction
Offshore geophysical exploration of hydrocarbon reservoirs using the marine controlled-source electromagnetic(mCSEM) method (Edwards, 2005) has its share oftechnical issues which are currently being addressed by various groups within the petroleum industry and academia. This paper describesthe development of a 3-D finiteelement analysis for addressing some of the significant problemsthat are associated with forward modeling of shallow-water mCSEM signals. This work is important because the availability of accurate and detailed 3-D forward modeling tools isa pre-requisite to reliable geological interpretation of mCSEM survey data.
A typical source-receiver survey configuration is shown in Figure 1. The cartoon shows the three major pathways for diffusion of the electromagnetic signal: through the seabed, through the seawater, through the air. It is useful to review the basic mCSEM physics directly in the time domain.
Figure 1. mCSEM: marine controlled-source electromagnetic exploration
of a hydrocarbon reservoir; image from Barker (2004).
In a transientmCSEM experiment, the source transmits a controlled signal. The diffusion time τ of the signal to a given receiver scales with the electrical conductivity σ of the medium and the source-receiver separation L by τ~μσL2, where μ is magnetic permeability. Since σ~0 in air, typically in shallow water mCSEM situations the first arrivalat the receiver is the “up-over-down” air-wave signal which can overwhelmslower signals diffusing through the seabed (σ1~0.1—1.0 S/m). The latter carrypotentially valuable information about the resistive hydrocarbon reservoir (σ2~10—2—10—4 S/m). The direct seawater signal is typically not problematic, since it arrives considerably later due to the high electrical conductivity (σ0~3.2 S/m) of seawater. The air-wave and seabed signals cannot be unmixed by data processing but must be modeled together by numerical simulation of the underlying Maxwell equations. Note that the air-wave signal is always the first arrival if σ0d2σ1D2where d is the water depth and D is the depth to the reservoir beneath the seafloor. For a reservoir burial depth ofD~500 m within conductive sediments, this corresponds to water depthd~300 m. We use this condition to defineshallow water.
Of particular concern to mCSEM practitioners are survey parameters that include shallow water depths, thin resistive reservoirs of finite extent, and large source-receiver separations. Computing high-frequency responses at large L in a shallow-water environment in the presence of thin resistors is challenging for several reasons: the receiver is effectively located at a significant number of skin depths from the source; the first-arrival is always the air wave; the signal amplitude from the reservoir is small.
The outline of the paper is as follows. First, the 3-D finite element methodology for computing frequency-domain (mCSEM) signals is briefly described, followed by a discussion of execution speed, memory requirements, and convergence rates. Then, a special modification of the modeling algorithm that is required to handlea shallow-water environment is discussed. After this, the capabilities of the finite element analysis are demonstrated for a few simple 3-D reservoir models. Lastly, an outlook is provided of the mCSEM forward modeling challenges that still lie ahead.
3-D finite element analysis
ThemCSEMgeophysical prospecting method is based on electromagnetic induction in the diffusive limit. The governing physics is typically castas a second-order differential equation in terms of the electric Eor magnetic B field. An equivalent formulation is made in terms of a vector and scalar potential pair (A,ψ) which is related to the fields by B=curl A and E=iω(A+grad ψ). The Coulomb gauge condition divA=0 is applied to uniquely specify the potentials. The equations satisfied by (A,ψ) are
(1) div2A + iωμ0σ(r) (A+grad ψ ) = —μ0JS ;
(2) div[ iωμ0σ(r) (A+grad ψ ) ] = 0 ,
where JS describes the source current density and σ(r) denotes an arbitrary three-dimensional conductivity structure. Our finite element analysis solves the weak form of (1) and (2) after separating the (A,ψ) potentials and the conductivity σ(r) into primary and secondary components. The unknowns are the Cartesian componentsof a piecewise linear approximation to (A,ψ) on the nodes of a cylindrical-shaped mesh composed of tetrahedral elements. The governinglinear system of equations is solved using the efficient QMR method (Freund et al., 1992). Readersinterested in further details and the evolution of our(A,ψ) finite element formulation are advised to pursue the chain of publicationsBiro and Preis (1989),Everett and Schultz (1996), Badea et al. (2001), and Stalnaker et al. (2006).
Figure 2. (left) SUB2, SUB4, and SUB8 tetrahedron refinement. The SUB2 and
SUB4 refinements are used to properly connect a refined region to the unrefined
region;(center) a 2-D horizontal projection of the 3-D cylindrical-shaped mesh;
(right)a vertical mesh projection.
A critical componentof mCSEMfinite element computations is the art of mesh design. There is a huge tradeoff between accuracy of the solution and cost and storage requirements, since these all increase with the number of nodes.The placement of the outer boundary of the mesh, upon which homogeneous Dirichlet boundary conditions are applied, is also very important. In an effort to increase the solution efficiency, we have implemented the local refinement algorithm ofLiu and Joe (1996) which readily permits multiple nested refinements of the mesh tetrahedra, as shown in Figure 2. The mesh is typically refined in the vicinity of a buried target at whose edges sharp gradients occur in the secondary (A,ψ) potentials. Somewhat surprisingly, we find it important to refine the mesh also in the vicinity of the transmitter (Stalnaker et al., 2006) where sharp gradients occur in the primary (A,ψ) potentials that drive the problem. Local mesh refinement to include arbitrary bathymetric variations in the forward model is readilyaccommodated.
Figure 3. (left) QMR convergence as a function of frequency; (center) mCSEM response
at 30 Hz for different numbers of QMR iterations; (right) effect of initial solution guess on
QMR convergence of the mCSEM response at 100 Hz.
The performance of the finite element code has been evaluated using a 3-D test structure relevant to mCSEM exploration consisting of a resistive disk(σ2=0.01 S/m) of radius ρ=8 km, thickness T=100 m and burial depth D=1.6 km. This simple reservoir target is excited by a horizontal electric dipole (HED) source located directly above the disk on the seafloor in shallow water, d=50 m. The host sediment is conductive, σ1=1.0 S/m. The outer radius of the cylindrical-shaped mesh is placed at ρMAX=16 km. The mesh extends 475 m above and 1.9 km below the seafloor. There are 325,000 nodes and 1.9 million tetrahedra in the mesh, which has not been locally refined. The CPU time to compute the mCSEM signals (EXEX-SEC, or the in-line electric dipole-dipole secondary response) associated with the resistive disk target is 10—30 minutes per frequency on a 2.4 GHz computer and the storage requirement is 16 GB. The convergence of the QMR linear solver, which occupies the bulk of program execution time, is displayed in Figure 3. It can be concluded from an inspection of the figure that the algorithm convergence slows down with increasing frequency due to degradation of matrix conditioning; the algorithm requires about 1200 iterations to converge at a nominal 30 Hz frequency; the convergence is faster if a good, non-trivial initial guess is chosen.
Shallow-water modification
The 3-D mCSEM finite element analysis, described above,is formulated in terms of Coulomb-gauged secondary potentials (A,ψ).The use of the Coulomb gauge is recommended since it leads to a symmetric finite-element matrix,unlike other possible gauges such as the Lorentz gauge, defined by div A=—iωμ0σΨ. The Coulomb-gauged primary potentials have been defined, thus far, to be the analytic (A,ψ) potentials appropriate for a uniform seabed halfspace overlain by a uniform seawater halfspace.
To incorporate the shallow water requirement in an efficient manner, we have re-defined the primary potentials to be those associated with a uniform seabed halfspace overlain by a seawater layer of finite thickness. Standard analytic solution techniques for the layered Earth problem can not be used since they would yield Lorentz-gauged primary potentials. To find Coulomb-gauge solutions, we have modified the Sommerfeld(1964, p.236ff) technique to determine the Hertz vectorΠand thenfind the (A,ψ) potentials by means of
(3) A = μ0σ Π +grad M ;
(4) —iω (M + Ψ) = div Π.
The novel feature is the introduction of an auxiliary scalar potentialMsatisfying
(5) div2M = — μ0σ divΠ ,
which fixes the Coulomb gauge.Without introducing the auxiliary scalar function M,
the shallow-water potentials (A,ψ) areimproperly gauged. A complete derivation of the shallow-water potentials shall appear elsewhere (Everett, 2006).
Simple 3-D reservoir models
To demonstrate the capabilities of the code as a powerful tool for offshore hydrocarbon geophysical exploration,we have modeled the 3-D shallow-watermCSEM response of various resistive targets which could represent an idealized petroleum or gas hydrate reservoir. In these examples, we show components of the vector potential A from which the measured (E,B) fields are readily found by differentiation. The source used is an infinitesimal horizontal electric dipole of unit strength lying on the seafloor, oriented in the x-direction. The reader should note that we have not made any attempt to model the details of any specific instrumentation nor have we tried to predict signal-to-noise ratios that would be measured during an actual mCSEM survey. These are important practical tasks that could be accomplished using the finite element code. The main point here is to show the seafloor patterns generated by the sub-seafloor targets and to show differences in these patterns as various target and host parameters are varied.
Figure 4. The secondary Ayshallow-water mCSEM response at 1 Hz
for three resistive targets. The HED source location is indicated
by the red arrow in the plan view image of the targets.
The effect of the shape of a resistive disk target on shallow-water (d=50 m)mCSEM responses is shown in Figure 4. The quantity plotted is the seafloor distribution on a 6 km × km grid of the logarithm of the secondary y-component of the vector potential A.The excitation frequency is 1.0 Hz. The disk, of radius ρ=8 km, thickness T=100 m and conductivity σ2=0.01 S/m, is buried at D=1.6 km depth in sediments of conductivity σ1=1.0 S/m. Notice that the full disk (right panel) generates a symmetric four-lobed quadrupole response, the half-disk (center) generates an asymmetric dipole response while the quarter-disk (left) generates an off-centered monopole response. These patterns are consistent with a scenario in which the induced current associated with each quarter of the disk acts as a secondary source of EM field.
The x-component of the vector potential A for the same set of targets is shown in Figure 5. These patterns can be interpreted in a similar fashion as for the y-component. Notice in particular that the quarter-disk response (left) is similar in the upper-right seafloor quadrant to the full-disk response (right) but attenuates rapidly in the other three seafloor quadrants due to the lack of a secondary source directly beneath those locations.
The effect of the size of the resistive quarter-disk target on shallow-water secondary Ax mCSEM responses is shown in Figure 6. As the quarter-disk radius decreases from ρ=8 km to ρ=2 km, going from left to right in Figure 6, the mCSEM response pattern decreases in size and becomes symmetric about the center of the target, which is offset from the source location. This trend appears to confirm the plausible scenario that fewer details about the target may be resolved from the seafloor mCSEM-response distribution as the target decreases in size.
Figure 5. The secondary Ax shallow-water mCSEM response at 1 Hz
for three resistive targets.
Figure 6. The secondary Ax shallow-water mCSEM response at 1 Hz
for the quarter-disk target: effect of target radius. The left panel contains the
same information as the left panel in the previous figure.
The effect of excitation frequency on shallow-water secondary Ax mCSEM responses of the resistive quarter-disk target is shown in Figure 7. As the frequency decreases from f=0.01 Hz to f=10 Hz, the mCSEM response pattern decreases in size accompanied by minor changes in shape. The response pattern does not change too much over the range of lower frequencies (top row), suggesting that a signal saturation has occurred since the mCSEM system is probing mainly beneath the target. Over the range of higher frequencies (bottom row) the response pattern diminishes rapidly since the mCSEM system is probing mainly above the target and hence cannot detect it. These conclusions are borne out by skin depth arguments: at 0.01 Hz, the disk is in the electromagnetic near-surface (D~7δ) while at 10 Hz it is electromagnetically deep (D~δ/4); where δ is the skin depth.
Figure 7.The secondary Ax shallow-water mCSEM response
for the quarter-disk target: effect of excitation frequency. The bottom left
panel contains the same information as the left panel of the previous figure.
The shallow-water secondary Ax mCSEM responses of a composite disk target are shown in Figure 8. This target, which is supposed to represent an inhomogeneous hydrocarbon reservoir, consists of a basal full-disk overlain by a half-disk and capped by a quarter-disk, as shown in plan view along the bottom of the figure. Each disk is T=100 m thick.The most substantial part of the reservoir, in this case the upper right quadrant of the target, is presumably the best drilling site. A glance at Figure 8 indicates that the symmetry of the mCSEM response pattern (left) changes as the half-disk is added to the familiar basal full-disk (center). The response pattern does not change appreciablywhen the quarter-disk is then added (right). The latter result clearly illustrates that minor alterations to the hydrocarbon reservoir model generate correspondingly minor changes to the mCSEM response pattern. Whether an actual mCSEM experiment has sufficient resolving capability to detect such subtle variations in the thickness of a hydrocarbon reservoir is beyond the scope of this article.
Figure 8. The secondary Ax shallow-water mCSEM response
for composite disk targets. The left panel contains the same information
as the right panel of Figure 5.
Figure 9. Comparison of the secondary Ax shallow-water mCSEM response
for two targets: quarter-disk (left) and prism (right).
The shallow-water secondary Ax mCSEM responses of two fundamentally different targets is compared in Figure 9. The response at left is due to the familiar quarter-disk. The response at right is due to a buried prism (D=800 m) of dimensions 6 km × 6m × 800 m thick. The size, shape and the location of the maximum of the mCSEM response pattern, in both cases, clearly reflects the size, depth of burial, and location of the target relative to the source.
Figure 10. The secondary Ax shallow-water mCSEM response
for a composite target consisting of a separated prism and half-disk.
A final example is provided to further explore the capabilities of the finite element analysis. In Figure 10, shallow-water secondary Ax mCSEM responses are displayed for a target that consists of two disconnected components: a prism and a half-disk. The mCSEM response pattern reflects the symmetry of the underlying target. The half-disk generates a dipole-like response in the lower two seafloor quadrants while the off-centered prism generates a single lobe in the upper left quadrant. This example illustrates that geometrically complex targets generate complex mCSEM responses. It is clear that 3-D forward modeling is required to interpret the mCSEMseafloor response.
Conclusions
This paper has described recent developments in 3-D finite element analysis for modelingmCSEM responsesof geometrically complex resistive targets embedded in the conductive seafloor beneath shallow seawater. The finite element algorithm is formulated in terms of Coulomb-gauged (A,ψ) electromagnetic potentials, from which the field vectors (E,B) can be recovered as required by appropriate numerical differentiation.
The range of resistive targets examined in this paper provides a sense of the capabilities of the 3-D finite element analysis for interpreting data from actual mCSEM experiments on offshore reservoirs. The code permits users to model complex geologicalscenarios such as stratigraphic or structuraltraps, salt dome traps, sand/gravel paleochannelsoverlain by non-reservoir rocks and bathymetry.For example, complex, mature multi-phase oil/gas/water reservoirs can be simulated as finely-discretized composite resistor/conductor targets using suitable local meshrefinement strategies
It is important to assess the performance of the 3-D finite element algorithm against other numerical methods. In general, the CPU time for differential-equation-based modeling (FE or FD) is strongly determined by the performanceof the linear solver. The iterative QMR method used here iscompetitive with the linear solvers found in the most advanced finite difference codes. Other numerical methods such as integral equations and spectral methods have been used for mCSEM forward modeling with varying degrees of success.
The CPU time requirement to run a 3-D forward modelon a 2.4 GHz platformis ~10-30 minutes per frequency. It is evident that parallel grid/cluster computing shall be needed to
interpret experimental data sets, regardless of whether trial-and-error forward modeling or a full inversion is utilized. Computer memory is a limitation to simulate mCSEMresponses.Modeling a detailed oil reservoir with realistic bathymetry requires ~16 GB.
Our forward modeling experience has highlighted the importance of careful mesh design. Local mesh refinement can be demonstrated to be a very effective tool to increase solution accuracy for moderate TX-RX ranges extending to several skin depths. For larger TX-RX ranges (>10 skin depths) however, local mesh refinement becomes ineffective and new techniques,such as multigrid Dirichlet conditions, are presently under investigation for improving solution accuracy.