Detecting induced polarisation effects in time-domain data - a modelling study using stretched exponentials

1Discretization¶

In this section, we discuss important elements about discretizing Maxwell’s equations in the time-domain with the convolution term shown in eq. (9), to simulate IP effects in time-domain EM data. Appendix A.1 illustrates how convolutionary time-domain Maxwell’s equations can be discretized. Appendix A.2 describes how the singularity of SE conductivity function at $t=0$ is handled. Most of key challenges about this discretization are tackled in Marchant (2015) (see page 21), and we have extended his work, applied for Cole-Cole conductivity, to SE conductivity.

1.1Maxwell’s equations¶

The stretched exponential (SE) conductivity provided in eq. (8) in the time-domain can be rewritten as

where $\delta(t)$ is a Dirac-Delta function and $\dsig(t)$ is

Considering a time-dependent conductivity, Ohm’s Law can be written as

and substituting eq. (1) yields

Using the Backward Euler method, we discretize Maxwell’s equations in eqs. (4) and (5) in time:

where $\triangle t^{(n)} = t^{(n)}- t^{(n-1)}$. To discretize the integral in eq. (4), we use the trapezoidal rule:

Fig. 1 shows a conceptual diagram for this discrete convolution procedure. Hence eq. (4) can be discretized as

This can be rewritten as

where the polarization current, $\j_{pol}^{(n-1)}$ is

For the simplest case when ($c_{se}=1$), then $\dsig(t=0)$ is well defined and $\gamma (\triangle t^{(n)})$ and $\kappa (\triangle t^{(n)})$ are respectively:

However, when $c_{se}\neq1$, $\dsig(t=0)$ is singular and hence it requires special numerical treatment; this is described in Appendix A.2.

For the discretization, we use a staggered mimetic finite volume approach Hyman et al. (2002). Here, boldface with uppercase and lowercase indicate matrices and column vectors, respectively. Further details about the discretization can be found in Haber (2014) (see page 31). Discretizing eqs. (5), (6), and (9) yields

where

Here, $\mathbf{C}$ is the discrete edge-curl operator; $\mathbf{M}^e$ and $\mathbf{M}^f$ are the edge and face inner-product matrices, respectively. For an inner-product matrix, the subscript indicates corresponding physical property (e.g. $M^{f}_{\mu^{-1}}$: the face inner-product matrix for $\mu^{-1}$).

Rearranging the above equations to solve for $\de$ yields:

By solving the above equation at each time step, we obtain $\de$. The measured data for AEM are often $-db/dt$ , which can be computed as

The measured data at a receiver loop can be expressed as

where $\mathbf{P}$ is an interpolation matrix, which projects $\mathbf{db/dt}$ fields, defined in a 3D domain, to a receiver location, and samples those fields at the measured time channels. For discretization of eqs. ( 17) to ( 19) we use, SimPEG’s mesh toolbox. The developed code is open-source as a SimPEG-EMIP package (https://github.com/sgkang/simpegEMIP)

1.2Handling the singularity at $\sigma(t=0)$¶

The SE conductivity, $\sigma_{se}(t)$ at $t$=0, is singular, whereas its integral is well-defined, as shown in eq. (6). When discretizing eq. (3), this singularity will be problematic. In particular, the issue occurs at the last time segment ($k=n$) of the convolution term in eq. (8), which can be written in continuous form:

This problem also occurs when the Cole-Cole function is used. Marchant (2015) (see page 31) tackled this issue by approximating $\e$ at this time segment as a linear function:

Then by substituting this in to eq. (20), and evaluating the integration, the discrete form of eq. (20) is obtained:

To obtain $\gamma(\triangle t^{(n)})$ and $\kappa(\triangle t^{(n)})$, we use the same trick. Integration of $\dsig(t)$ is not possible, so by Taylor expanding, we obtain an approximate form of $\dsig(t)$ which is valid for small $t$:

By substituting eq. (23) into eq. (20) and evaluating the integral, we finally obtain

2Analytic test¶

To test the developed SimPEG-EMIP code, we compare our numerical solution with an analytic solution. A halfspace earth is assumed. The conductivity of the halfspace is 0.05 S/m and its SE parameters are: $\eta_{se}$=0.7, $\tau_{se}$=4ms, $c_{se}$=0.6. Corresponding Cole-Cole parameters are: $\eta_{cc}$=0.8, $\tau_{cc}$=0.005s, $c_{cc}$=0.6. For the spatial discretization, a 2D cylindrically symmetric mesh is used; the smallest cell size is 6.5m × 5m. A horizontal source loop is located 30m above the surface. A step-off waveform is used for the input current and a horizontal receiver loop measuring the voltage (equivalent to -$db_z/dt$) is coincident with the source loop. Data are measured in the off-time over the time-range: 10^-2-10 ms. Fig. 2 shows comparison between analytic and numerical solutions; they match well except for a small shift in the time of the zero-crossing, the two solutions are in good agreement.