3 Subsurface Properties and Propagation

In this section we expand upon Question 2 and 3: what is the measured signal, and how does one deduce the earth’s response.

1. What physical properties does one wish to measure?
2. How does one get a transfer function from the measured signal?
3. How does one get the desired properties from the transfer function?

At various times Vendor has claimed his method can distinguish both conductance and dielectric properties with extraordinary resolution in depth. At other times he has backed off a bit on those claims. We begin with a short digression on the electrical properties of selected rock and their influence upon subsurface propagation of the EM field, then continue with a brief discussion of the frequency domain reflection coefficient $R\left(\omega \right)$.

3.1 Stones in the Pathway

Maxwell’s Equations (see e.g. [6, ] or Appendix A), depend upon three material properties: permittivity (dielectric constant) $ϵ$, conductance $\sigma$, and magnetic permeability $\mu$. Outside of some truly fabulous ferromagnetic deposits, magnetic permeability is fairly constant and may be approximated by that of free space, ${\mu }_o=1.257\times 10^{-6}\text{ H/m}$ eq. (3). We shall see that, for most formations of hydrocarbon interest and practical signal frequencies, the effect of conductivity $\sigma$ dwarfs permittivity. However, we do not take that as an a priori assumption.

3.1.1 Permittivity

3.1.2 Conductance

Some sedimentarilly relevant approximate conductivities (from [6, Figs. 3.17, 3.39]):

As an (optimistic) example, take ${\sigma }_t=10^{-4}$ S/m, $f=10^4$ Hz, and ${ϵ}^{\mathrm{rel}}=100$. The loss factor is

which is still greater than unity. More realistically, when ${\sigma }_t\sim 10^{-3}$ S/m, $f\sim 10^4$ Hz, and ${ϵ}^{\mathrm{rel}}\sim 50$., then $\delta \sim 36$, twenty times larger.

3.1.3 Propagation in Lossy Material

The loss factor $\delta$ is a fundamental consequence of Maxwell’s Equations in lossy media (Appendix A), and its importance cannot be overstated. Large values of $\delta$ transform the subsurface EM propagation from purely wave-like in the atmosphere, to something diffusive, midway between the heat equation and Schrödinger’s. Large values of $\delta$ may decrease the phase velocity, nominally given as $V_p=\omega ∕\text{ Re}\left(k\right)$ ([6, eq. 12.38 ff]) from its speed-of-light value $3\times 10^8$ m/s in free space immensely. From the plane-wave representation ([6, eq. 12.30 ff]):

Equation for wave front:

(See [6, 12.42]) Since $\delta$ is a function of $\omega$, so is $V_p$: conductive media are dispersive. The penetration factor ${\delta }_e$ is defined by

from which

so that a signal from the surface and reflecting from a depth ${\delta }_e$ will return with a strength $1∕\left(e^2\right)=0.135$ of its surface value, a difference presumably one can distinguish. For large values of $\delta$, wavelength $\lambda \simeq 2\pi {\delta }_e$.

Consider an extreme example of MT signal with period 24 hours. $\omega =2\pi ∕\left(24\ast 3600\right)=7.27\times 10^{-5}$. Take ${ϵ}^{\mathrm{rel}}=50$ and $\sigma =0.01$ as volume averages. Then

Those last two numbers are why subsurface EM propagation may be termed ”diffusive”. A signal from the diurnal variation creeps along at 370 kph, taking nearly 4 hours to reach its penetration depth in the mid-mantle. As a rule-of-thumb

where ${\rho }_a$ is the bulk average resistivity ($\rho =1∕\sigma$).

Eq. (44) is something to keep in mind when planning an EM survey or evaluating its results. At the September survey site, for example, the geologists estimate apparent resistivity ${\rho }_a\sim 2\Omega$m, and Vendor states his pass-band as 700 Hz – 10 KHz. In which case ${\delta }_e\approx 500\sqrt{2∕700}\approx 27$m. But depth to shallowest target is $\sim 1$ km. And $e^{-1000∕27}\sim e^{-37}\sim 10^{-16}$, quite effectively zero. (That’s just the attentuation of surface source signal to target depth. Square that number for whatever you want back.) At such low resistivity, 1 km depth requires signal frequencies of $\sim 1$ Hz. Preferably less. Even if the data receiver was capable of recording sub-audio frequencies, the sixty second acquisition window is short of what is recommended for 1 Hz, by nearly a factor of two.

Working the other way, suppose one were to propose an EM sounding (of any kind) to 2 Km depth in an area where one suspects ${\rho }_a\sim 2$m. Then $T_{\text{ max}}\simeq {\left(d∕500\right)}^2∕{\rho }_a={\left(2000∕500\right)}^{∕}2=8$ seconds. (Of course, one would want to acquire higher frequencies, and possibly lower, as well.) A rule-of-thumb minimum acquisition time is 130 times this (for the fourier transform): 1040 seconds, or 17 minutes. The 30 minute time window for the field variations exemplified in Fig. (2)  becomes somehat more relevant.

3.2 Reflection Coefficient

Vendor has sometimes described his processing flow as

At other times he has mentioned signal smoothing with a ”Ten (or twenty or forty) foot filter”. Although we do not pretend to deduce – and certainly do not reproduce – Vendor’s actual algorithm, both of these statements are consistent with performing a Fourier Transform of the recorded data time series, assuming some relation between frequency and depth, and applying a simple sliding window smoothing filter. We have not yet reviewed the patents, and do not know that is what is being done, but Fourier transformation to frequency domain is an extremely common (almost universal) step in signal processing. In particular, if one knows the time waveform of the source signal, or equivalently its Fourier transform, then Fourier transformation of the recorded signal is the first step (after analog filtering) toward deconvolving the recorded data from the source to obtain the Earth’s transfer function or response. For measurements made at or above the Earth’s surface, this response is expressed in the frequency domain by the Reflection Coefficient $R\left(\omega \right)$.

3.2.1 Reflection at a planar interface

Computing the reflection coefficient for a plane wave normally incident upon an homogeneous layered earth model is straightforward. Reflection at a planar interface is determined by the condition that the tangential components of $\overrightarrow{E}$ and $\overrightarrow{H}$ be continuous across the boundary. There are two immediate results: the reflection coefficient, and Snell’s Law that gives the angle of refraction $\psi$ in terms of the angle of incidence $\theta$ upon boundary between material 1 and material 2 ([2, eqs 7.3]):

where ${ϵ}_m^{\star }={ϵ}_m-i{\sigma }_m∕{\omega }_m={ϵ}_m\left(1-i{\delta }_m\right)$ for dissipative (conductive) media in layer $m$. Snell’s Law of refraction becomes a bit complex, but for dielectric layer (e.g. air) above strong dissipative media (typical earth) where ${\delta }_2\gg 1$, one obtains ([2, pg. 305]):

since ${\mu }_1\simeq {\mu }_2\simeq {\mu }_0$, with the result that $\psi \to 0$ as $\sigma \gg 2\omega {ϵ}_1$ (but the square root slows the effect.) This is a somewhat good thing, as it means that plane waves of even moderately normal incidence (e.g. 45 or 50 deg.) may be treated as essentially normal in the subsurface of the homogeneous layered earth model. This condition holds for MT, and probably most AMT events, but about power grid noise, particularly near transmissions line, I don’t yet know. I hope Vendor can supply some guidance. (See Question 1.) We assume normal incidence in the present analysis. The limiting condition $R\to -1$ as $\omega \to 0$ (eq. 48) means that any conductor completely reflects DC (or near DC) fields, out-of-phase, so that $\overrightarrow{E_{\parallel }}\left(\omega =0\right)=0$ . This relates to the DC assumption that there is no steady-state homogeneous subsurface current. (There may be near-DC tellurics, but these cannot be homogeneous and fall outside the 1-D model.) The reflection coefficient for a plane wave normally incident upon a uniform half-space is plotted in Figure 3.3.

3.2.2 Reflection from a slab

Reflection from a uniform slab embedded between two homogeneous half-spaces is also straightforward. We take our results from [2, Sec. 6.9], to which we refer for details. The result is the electromagnetic field components on one boundary (e.g. the the earth surface) may be obtained from their values at the next via a simple propagator:

where the propagator $K_m$ is a $2\times 2$ matrix:

At normal incidence ($\theta =0$)

Introduce the electromagnetic impedance $Z\left(\omega \right)=E\left(\omega \right)∕H\left(\omega \right)$. The downward-going radiation condition in the bottom layer $M$ requires $Z_M={\eta }_M$, from which

where for $M>2$ the propagator is the product matrix

For a single slab between two half-spaces $M=2$ (two boundaries), and $K$ takes the simple form (51). The impedance $Z$ and reflection coefficient $R$ are mathematically equivalent:

(See [2, eq 9.9].) If region 1 is air, then

Reflection coefficients for resistive slab are plotted in Figures 3.3 and 3.3, as that is what is (at best) measured in Vendor’s experiment. In this context ”mathematical equivalence” can be misleading. The impedance $Z$ describes the ratio of the electric and magnetic fields, whereas the reflection coefficient $R$ may be applied to each field, electric or magnetic, individually. An astute observer might note that, if one were to measure the total electric field and the total magnetic fields separately but simultaneously at the earth’s surface, then nearly all the noise and spectral density problems posed in Question 1 disappear when one takes their ratio to obtain the impedance $Z$.

A.N. Tikhonov and L. Cagniard were ”astute” observers.

3.3 Illustrating the Point

Our numerical simulations are quite simple and assume permittivity and conductance are independent of frequency. This is not strictly true. These effects are discussed in [6, Chapter 3] and, to whatever degree present in the audio band, may be expected to modulate our modeled reflection coefficients to an as yet unknown degree.

In the following figures we parametrize the layered earth model by $Frac:0\le Frac\le 1$, the fraction of the model’s $\sigma$ values (the conductances of each layer) used in the sub-model for a particular curve. This helps illustrate the effects of increasing conductance upon the reflection coefficients of an earth model. The important result is that variations in $R\left(\omega \right)$ decrease markedly with increasing conductance, such that at $\sigma =0.001$ the reflection coefficient becomes an essentially monotonic function of $\omega$. This suggests that unless Vendor can propose a (geologically) reasonably conductive earth model that gives a significantly different result, his claim the recorded signal contains detailed information within its second derivative should be scrutinized for the possible influence of noise.

3.4 Free Modes of Decay

An earth model’s internal modes or free modes are defined by radiative boundary conditions: (1) the mode may propagate only downward at great depth, and (2) may propagate only upward at the earth’s surface. The ”downward at depth” condition is a requirement of all realizable Maxwell solutions under the no deep source current” assumption, but the ”upward only” condition at the surface means the free modes have no above-surface excitation either, i.e. $A\left(\omega \right)\equiv 0$ in eq. (1). Therefore these modes may be excited only by subsurface sources (but not really deep ones). This condition is readily realized in controlled-source Transient Electromagnetic Methods (TEM), wherein e.g. a high-current DC magnetic dipole is energized at the surface for a sufficient time for its magnetic field to penetrate to depth. The source current is then abruptly turned off (not as easy as it sounds). The magnetic field the source has set up in the earth then decays to zero, and induces subsurface current transients as it does so. These currents depend on the conductance structure of the earth, and temporarily excite the transient modes, or free modes, of decay.

Free modes of decay are nothing new, being discussed 30 years ago by Schmucker and Weidelt [4, Section 2.2]. Those authors proved that the free mode eigenfrequencies are all purely imaginary for very general earth models, not just for uniform layers. Such modes are purely damped: they have no oscillatory component.

A few items to note:

• First, Schmucker and Weidelt’s ”purely damped” proof was obtained for Maxwell’s Equations in their inductive limit where permittivity is ignored in deference to the relatively large conductivity and loss factor $\delta$. One might suggest this approximation is not necessarily valid, and for relatively dry earth models $\left(\sigma <0.0001\right)$ one might be right.
• Second, the controlled step-source described above is very broad-band, the step function at current cut-off ideally having infinite frequency content. Nothing is ideal, of course, and it is not possible to instantaneously stop the current flow through an inductor. Any other external source with proper frequency content may also excite free modes as it propagates from the surface downward through the earth. Such sources include MT, AMT, and presumably ”power-grid noise” as well. However, it seems unlikely the magnitude of the excited free modes may ever exceed that of the excitation source as they reflect back at the surface: in the frequency domain their effect is completely described by the reflection coefficient $R\left(\omega \right)$. Referring to eqs (1) and (2), we see that $A\left({\omega }_{\text{ f}}\right)=0\Rightarrow R\left({\omega }_{\text{ f}}\right)=\infty$ i.e., the free mode eigenfrequencies ${\omega }_{\text{ f}}$ occur at the complex poles of the reflection coefficient. Importantly, $R$ is an analytic function of $\omega$.
• Third, in his presentations Vendor has frequently stressed an analogy between his method and Singularity Expansion Methods (SEM) used in Radar and Sonar target characterization. The long dimensions of an aircraft’s wings and fuselage may support damped creeping mode resonances at radar frequencies. A ship’s hull may exhibit damped ringing at audio frequencies characterized by its material and dimensions, much like a bell. Such modes may be excited by either transient or continuous source. If the received signal can be deconvolved from the source, the singular eigenfrequencies and their relative excitations can yield information about the character of the target.

In sufficient absence of noise, Prony-type fits of the deconvolved time-series have been used to obtain the singular frequencies and their excitation magnitudes from radar and sonar targets. Pending further research, we cannot say whether Prony techniques are similarly applicable to purely-damped TEM data: there is a branch cut on the $\text{ Im}\left(\omega \right)$ axis that might be expected to contribute as well, complicating the interpretation.

Singular mode expansions may be illustrated by a simple astrophysical analog: it turns out the Regge-Wheeler potential surrounding a black hole supports free gravitational modes that may be excited by in-falling matter. Such mass may have been orbiting close to the hole for some time, and in this sense truly does constitute a source that is ”internal” to the potential region. The Regge-Wheeler potential also supports electromagnetic free modes, and while the distribution of a black hole’s eigenfrequencies in the complex frequency plane differs somewhat from those of free modes in the earth or singular modes of Radar/Sonar targets, they all support under-damped oscillatory modes (for a sufficiently dry earth).

The free-mode contribution to the Green’s Function propagator may be obtained by similar methods for each case, and in this regard we cite reference [3, ]. That paper’s Figure 2 (reproduced here as Fig. 3.4) illustrates a numerical time-domain reflection experiment that is in some ways analogous to radar reflectometry. Here a very simple impulse of finite duration (Gaussian peak in time) is directed at a Regge-Wheeler potential and the reflected waveform observed. The reflection is of lower amplitude and energy content than the source impulse. The complex eigenfrequencies of the quasi-normal (free-mode) ringing – if they can be deduced – may be used to characterize the details of the potential and/or the source. In radar and sonar the source waveform is known with high precision, and any ringing in the returned signal may be deconvolved from the source and used to characterize the target. However, in the particular case of Vendor’s experiment, the source is not of finite duration and its waveform is (thus far) not known. Deconvolution is therefore (probably) not possible, and even if the reflected signal contained free-mode ringing (and for a reasonably conductive earth the free modes are damped exponentials and not oscillatory), it will be (thoroughly) mixed with subsequent incoming (unknown) signal at the time it is reflected back to the receiving antenna at the earth’s surface. Deconvolution of a signal from a single antenna is not possible unless one knows something about the source.

Figure 16: Reproduction of [3, Fig. 2], showing results of a numerical time-domain reflection experiment. Initial value problems are properly addressed via Laplace Transform. Here $A_{in}\left(s\right),A_{out}\left(s\right)$ correspond to our present $A\left(\omega \right),B\left(\omega \right)$, where $\omega =-is$.

3.4.1 Free Mode Eigenfrequencies by Newton Root Search

It should be noted that the final radical in eq. (54) introduces a branch cut between $\omega =0$ and $\omega =i{\left(\sigma ∕ϵ\right)}_{\text{ max}}$, where ”max” indicates the maximum of the ${\sigma }_m∕{ϵ}_m$ taken over all layers. This branch cut might introduce complexities in attempts to interpret purely damped free-modes whose eigenfrequencies fall upon (or near) the branch. One obtains an implicit equation for the free mode eigenfrequencies by setting $Z_0=-Z_{\text{ free}}$ in eq. 56 or by setting $R_0=\infty$ in eq. 60. These two equations suggest two approaches for finding the eigenfrequencies: either (1) find the zeros of $f\left(\omega \right)=Z_0\left(\omega \right)+Z_{\text{ free}}$, or (2) find the poles of $R_0\left(\omega \right)$ directly. The first involves non-linear root search requiring numerical differentiation. The second involves numerical contour integration in the complex $\omega$ plane. Both have their advantages and disadvantages: we have implemented both.

The free mode eigenfrequencies for the slab-between-two-half-spaces may be obtained by setting $Z_0=-Z_{\text{ free}}$ in eq. (56) with $K$ given by the single-layer expression (51) with $m=1$ for the slab. After some algebra one obtains

where

For a slab embedded in air ${\sigma }_0={\sigma }_2=0$ and ${ϵ}_2={ϵ}_0$. There are two limiting cases: ${\sigma }_1=0$ (dielectric slab), and ${\sigma }_1\gg {ϵ}_1$ (conductive slab). In the first ${ϵ}_m^{\star }={ϵ}_m$ has no $\omega$ dependence and (63) gives the eigenfrequencies directly:

Regrettably (or not), we have not yet been able to find a similarly simple expression in the $\text{ Re}\left(\omega \right)\to 0$ limit. If ${\delta }_1\gg 1$ then one might hope to approximate ${ϵ}_1^{\star }\simeq -i{\sigma }_1∕\omega$ and ${\xi }_1\simeq d_1\sqrt{-i\omega {\sigma }_1{\mu }_1}$, and rearrange eq. (63) to give

suggesting the eigenfrequencies ${\omega }_n$ in the ${\sigma }_1\gg {ϵ}_1$ limit (71) indeed become purely imaginary, in agreement with Schmucker and Weidelt’s general proof. But there are problems in the numerical verification. These arise because as long as $\text{ Re}\left(\omega \right)>0$, numerical root search always finds $\left|\omega \right|$ large enough that ${\delta }_1$ is on the order of unity or less, so that the requisite ${\delta }_1\gg 1$ limit for (71) to hold never obtains. In fact, the numerical root search result for a simple conducting slab embedded in air suggests that, for this arrangement, the magnitudes of the eigenfrequencies remain essentially constant, as illustrated in Figure (3.4.1).

In Figure (3.4.1) we plot the trajectories the Free Mode Eigenfrequencies $f_f$ make in the complex frequency plane for the 4-layer resistive-slab-overlaying-thin-insulator model of Fig. 3.3. The trajectories are again parametrized by $Frac:\left(0\le Frac\le 1\right)$, with the conductances of all layers being multiplied $Frac$. The permittivities of the deeper layers are scaled such that they all match that of the thin insulator at $Frac=0$, and progress to their final desired values as $Frac$ increases to 1. In other words, the conductivity of the top 1km thick layer increases from 0 to 0.001, while that of the damp basement increases from 0 to 0.0001. The permittivity of the basement increases from 9 (that of the thin insulator) to 25. This allows the trajectories to start from the simple analytic result of eq. (69). The method is relatively rapid and applicable to a fairly large frequency range, with the caveat that beginning frequency searches at the simple three-layer model’s known values may miss solutions for many-layer models that have no simple-model counterpart. (See Sec. (3.4.2)).

To give a flavor for how the eigenfrequencies change with $Frac$, cross-hatches are made at values of $Frac:\left(0,0.1,0.2,...1.0\right)$, with the $Frac=0$ (non-conducting, $\sigma =0$) values lying closest to the real $f$ axis in agreement with eq. (69). As $\sigma$ increases, $\text{ Im}\left(f_{in}\right)$ increases, while $\text{ Re}\left(f_{in}\right)$ decreases. The imaginary parts of each harmonic are essentially the same.

In the audio range the low-harmonics have all collapsed to the $\text{ Im}\left(f_{in}\right)$ axis. In the trans-audio and RF range the real part of the higher-overtone $f_{in}$ persists, but the damping caused by the imaginary part is extreme. In any event, the penetration depth of a 50 KHz signal in the model’s top-most (${ϵ}^{\mathrm{rel}}=100,,\rho =1000,d=1$ km layer is only ${\delta }_e\simeq 500\sqrt{\rho ∕f}\sim 70$m.

3.4.2 Free Mode Eigenfrequencies by Contour Integration

To illustrate the pitfalls inherent in the ‘Start from a known solution and boldly set forth into the Great Unknown” approach described above, and from the form of the propagator product eq (57) rather suspecting bifurcation of the eigenfrequency trajectories of Fig. (3.4.1) as additional layers are introduced, out of an abundance of caution we have additionally sought reflection coefficient poles by dividing the audio region of the complex $R\left(\omega \right)$ plane into small squares, initially 120 Hz per side, and integrating around each in search of a residue within. Squares with non-zero residue were further subdivided to obtain $\pm 6$ Hz accuracy. We present the result in Fig. (3.4.2), which shows what transpires when one adds an additional thin insulating slab beneath the conducting slab embedded between air and a slightly conducting basement of Fig. (3.4.2).

With the introduction of the thin insulator we see the appearance of a new highly-damped trajectory branch, starting at $\omega =\left(4414,50942\right)$ with $Frac=0.18$ , and continuing rightward and upward with increasing Frac until the branch disappears off scale at $\left(20157,62327\right)$ with $Frac=0.28$. A few caveats:

1. The magnitudes of residues of the poles on this retrograde branch start at $\sim 0.1\%$ the magnitudes of the lesser-damped poles, and decrease with increasing $Frac$, which does not bode well for strong contribution from these modes to the Green’s Function.
2. We have not followed this branch to higher values of $\sigma$ and $Frac$. The Schmucker-Weidelt theorem only says the free mode eigenfrequencies have no non-zero real part in the $\delta \gg 1$ limit, but gives no details upon how they get there. This retrograde branch might circle back to the imaginary $\omega$ axis with increasing $\sigma$, or the magnitudes of the pole eigenfrequencies might continue to increase such that the $\delta \gg 1$ limit is never obtained. Or the branch may die out if the residues at its poles fall to zero. Our particular contour integration routine becomes very inefficient when the imaginary part of $\omega$ becomes large, which is why we’ve limited this plot to $\left(32000,64000\right)$. The non-linear Newton search used in Fig. (3.4.1) could be employed to extend the frequency range, but there is probably little point: we are talking here about very highly damped modes with very small residues (less than unity magnitude) that appear to be headed well outside the audio band, and all this at still geologically tiny conductance. But the question remains.

3.4.3 Summary

The behavior of the Reflection Coefficient poles in the complex $\omega$ plane can be deceptively fascinating. It must be kept firmly in mind that data is recorded in the time domain. While it is usually very useful to decompose the time-domain signal into its Fourier frequency components, the measured Fourier frequencies are real. Alternatively, in transient-source cases such as sonar and radar when one knows a priori that the deconvolved response signal has significant damped sinusoid components, or as in TEM methods when one knows a priori the deconvolved response has significant sums of exponential decay (purely damped inductive free modes), a Prony-type fit of the deconvolved time series may be of benefit. But those cases do not apply here, where Vendor’s method – as do MT and VLF techniques – relies upon a continuous source signal. Now Fourier analysis of the recoded signal is most appropriate, and if the recorded signal may be deconvolved from the source, the reflection coefficient obtained for real $\omega$. Such reflection coefficients are plotted in Figures (3.3 – 3.3). Complex pole analysis is useful only when $R$ shows local extrema for real $\omega$, as exemplified by the $\rho \ge 10,000\Omega$-m portions of the Figures. However, resistivities in sedimentary basins are typically less than 100 ohm-meter; the reflection coefficent for the very simple earth models we have devised here do not show such variation even for resististance ten times this value. We hope Vendor can suggest models that do.