Background Theory¶

The background theory presented here is copied entirely from the Master’s thesis “Forward modeling and inversion of viscous remanent magnetization responses in the time domain” (Cowan, 2016). See this thesis for a more thorough treatment of the background theory.

General¶

To first order, and during the off-time, the Earth’s inductive and viscous remanent magnetization (VRM) responses are separable so long as the Earth’s magnetic susceptibility is sufficiently small (Druyts et al., 2009); i.e. \(\chi < 0.1\). This was first postulated by Buselli (1982); who used various SiroTEM configurations to perform field experiments. The approximation was later validated numerically by Pasion (2007), among others. At early times, the response is equivalent to that of a purely conductive Earth (red); at late times, the response is equivalent to that of a purely viscous Earth (purple). By summing the responses from these two models, the response from a conductive and magnetically viscous Earth can be approximated. Thus:

\[\vec{b}(t) = \vec{b}_{em}(t) + \vec{b}_{vrm}(t)\]

Vertical step-off response from a horizontal circular current loop over a conductive and magnetically viscous halfspace. (a) \(b_z\). (b) \(\partial b_z/\partial t\).¶

Here, we develop a methodology for modeling \(\mathbf{b_{vrm}}\); as methodologies currently exist for modeling \(\mathbf{b_{em}}\). According to Blakely (1996), the anomalous magnetic response \(\vec{b} (P ,t)\) at location \(P\) outside of a magnetized region is given by the following integral equation:

(1)¶\[\vec{b}(P,t) = \frac{\mu_0}{4\pi} \int_{V} \nabla \nabla \frac{1}{r} \cdot \vec{M} (Q , t) dV\]

where \(\vec{M} (Q,t)\) is the magnetization at locations \(Q\) within the source region \(V\) and \(r\) is given by:

\[r = \sqrt{\big ( x_P - x_Q \big )^2 + \big ( y_P - y_Q \big )^2 + \big ( z_P - z_Q \big )^2 }\]

SimPEG methods for modeling VRM responses in the time-domain are numerical solutions to the above equation.

Dispersive Magnetic Susceptibility¶

Magnetic susceptibility defines the degree of induced magnetization experienced by a material, proportional to the strength of an applied magnetic field. In magnetically viscous rocks/soils, the induced magnetization depends on the frequency of the applied field (Lee, 1984; Dabas et al., 1992; Fannin and Charles, 1995; Meglich et al., 2008), thus:

\[\vec{M} (\omega )= \chi (\omega) \vec{H}_0 (\omega)\]

where \(\vec M\) is the magnetization, \(\vec{H}_0\) is the applied magnetic field and \(\chi\) is a frequency-dependent magnetic susceptibility. Commonly used models for the frequency-dependent magnetic susceptibility are derived by applying a weighting function \(f(\tau )\) and integrating over all Debye models:

\[\chi (\omega) = \chi_\infty + \Delta \chi \int_0^\infty \frac{f (\tau )}{1 +i \omega\tau} d\tau\]

where \(\tau\) is the time-relaxation constant, \(\chi_\infty\) is the instantaneous magnetic susceptibility and \(\Delta \chi\) is the DC contribution due to magnetic viscosity.

Frequency-dependent magnetic susceptibility assuming various distributions of time-relaxation constants.¶

As we can see, the magnetic viscosity is defined by the choice in weighting function. The majority of soil samples can be adequately fit by assuming a log-uniform distribution of time relaxation constants between a set of finite limits \(\tau_1\) and \(\tau_2\) (Dabas et al., 1992; Worm, 1998; Igel et al., 2012). This leads to a magnetic susceptibility of the form:

\[\chi (\omega) = \chi_\infty + \Delta \chi \Bigg [ 1 - \frac{1}{\textrm{ln} (\tau_2 / \tau_1) } \textrm{ln} \Bigg ( \frac{1 + i\omega\tau_2}{1 + i\omega\tau_1} \Bigg ) \Bigg ]\]

In the time-domain, the relationship between the applied magnetic field and the resulting magnetization becomes a convolution:

(2)¶\[\vec{M}(t) = \chi (t) \otimes \vec{h}_0 (t) = \int_{-\infty}^\infty \chi (\xi) \vec{h}_0 (t - \xi) d\xi\]

where \(\vec M (t)\) is the induced magnetization, \(\vec{h}_0 (t)\) is the applied magnetic field, and \(\chi (t)\) represents the magnetization’s impulse response. For a log-uniform distribution of time-relaxation constants, this is given by:

(3)¶\[\begin{split}\begin{align} \chi (t) & = \chi_\infty \delta (t) + \Delta \chi \, u(t) \int_0^\infty f(\tau) \frac{e^{-t/\!\tau}}{\tau} d\tau \\ & = \chi_\infty \delta (t) + \frac{\Delta \chi}{\textrm{ln} (\tau_2 / \tau_2)} \Bigg [ \frac{e^{-t/\!\tau_2} - e^{-t/\!\tau_2}}{t} \Bigg ] \, u(t) \end{align}\end{split}\]

where \(u(t)\) is the unit step function and \(\delta (t)\) is the Dirac delta function. In the previous equation, the first term characterizes instantaneous magnetization and the second term characterizes magnetic viscosity.

General shape of the impulse response.¶

Numerical Solution¶

The anomalous magnetic response outside of a magnetized region is given by Eq. (1). Assuming the magnetic response is the result of a set of uniformly magnetized rectangular cells, Eq. (1) is approximated by the following (Bhattacharyya, 1964; Varga, 1996):

(4)¶\[\vec{b} (P , t) = \sum_{k=1}^K G^k \cdot \vec{M}^k (t)\]

where \(\vec{b} (P , t)\) is the time-dependent magnetization observed at \(P\), \(\vec{M} (t)\) contains the vector components of induced magnetization within each cell \(k\) at time \(t\) and \(G\) are \(3 \times 3\) tensors representing a linear operator between the magnetization of the cell \(k\) and the anomalous magnetic response at location \(P\). Thus \(\vec{b} (P , t)\), \(G\) and \(\vec{M} (t)\) are given by:

\[\begin{split}\vec{b} (P , t) = \begin{bmatrix} b_x (P,t) \\ b_y (P,t) \\ b_z (P,t) \end{bmatrix} \;\; , \;\; \vec{M} (t) = \begin{bmatrix} M_x (t) \\ M_y (t) \\ M_z (t) \end{bmatrix} \;\; \textrm{and} \;\; G = \begin{bmatrix} g_{xx} & g_{xy} & g_{xz} \\ g_{yx} & g_{yy} & g_{yz} \\ g_{zx} & g_{zy} & g_{zz} \end{bmatrix}\end{split}\]

Eq. (4) can be re-expressed as a linear system. And in the case that there are multiple observation locations, this system may be augmented and expressed as:

(5)¶\[\mathbf{b}(t) = \mathbf{G} \, \vec{\mathbf{M}}(t)\]

where \(\mathbf{b}(t)\), \(\vec{\mathbf{M}}(t)\) and \(\mathbf{G}\) take the form:

\[\begin{split}\mathbf{b}(t) = \begin{bmatrix} \mathbf{b_x}(t) \\ \mathbf{b_y}(t) \\ \mathbf{b_z}(t) \end{bmatrix} \;\; , \;\; \vec{\mathbf{M}} (t) = \begin{bmatrix} \mathbf{M_x}(t) \\ \mathbf{M_y}(t) \\ \mathbf{M_z}(t) \end{bmatrix} \;\; \textrm{and} \;\; \mathbf{G} = \begin{bmatrix} G_{xx} & G_{xy} & G_{xz} \\ G_{yx} & G_{yy} & G_{yz} \\ G_{zx} & G_{zy} & G_{zz} \end{bmatrix}\end{split}\]

Thus \(\mathbf{b}(t)\) contains the cartesian components of the magnetic response at all observation locations, \(\vec{\mathbf{M}}(t)\) contains the cartesian components of induced magnetization within all of the cells at time \(t\) and \(\mathbf{G}\) is a linear operator. The time-derivative of the response is given by:

\[\frac{\partial \mathbf{b}}{\partial t} = \mathbf{G} \, \frac{\partial \vec{\mathbf{M}} }{\partial t}\]

The magnetization experienced by each cell \(k\) is given by the following convolution:

(6)¶\[\vec{M}^k (t) \approx \chi^k (t) \otimes \vec{h}_0^k (t)\]

where \(\vec{h}_0^k (t)\) is the free-space field at the center cell \(k\) due to an inductive source and \(\chi^k (t)\) is the impulse response for cell \(k\). Here we assume 2 things:

That the inducing field is approximately homogeneous throughout the cell

That because we are modeling a purely viscous Earth, attenuation of the primary and secondary signal are negligible

Linear Problem¶

The linear formulation allows for rapid computation of VRM responses during the off-time; see Cowan, (2016). Assuming the distribution of time-relaxation constants is sufficiently broad for a log-uniform distribution (i.e. \(\tau_1 \ll t \ll \tau_2\)), then the off-time VRM experienced by each cell is given by:

(7)¶\[\vec{M}^k (t) \approx m^k \, \vec{h}_{max}^k \, \eta (t)\]

where \(\vec{h}_{max}^k\) is the inducing field at its maximum amplitude, \(\eta (t)\) is a purely time-dependent function which defines the decay for any cell and \(m\) is an indirect measure of magnetic viscosity and is related to \(\Delta \chi\), \(\tau_1\) and \(\tau_2\) by:

\[m = \frac{\Delta \chi}{ \textrm{ln} (\tau_2 / \tau_1)}\]

Using Eq. (7), \(\vec{\mathbf{M}}(t)\) in Eq. (5) becomes:

(8)¶\[\vec{\mathbf{M}}(t) = \eta (t) \, \mathbf{H_0 \, m}\]

where \(\mathbf{H_0}\) is a sparse matrix containing the cartesian components of \(\vec{h}_{max}^k\) for every cell and \(\mathbf{m}\) contains the parameter \(m^k\) for every cell. Substituting Eq. (8) into Eq. (5):

\[\begin{split}\begin{align} \mathbf{b}(t) &= \eta (t) \mathbf{G \, H_0 \, m} \\ & = \eta (t) \mathbf{A \, m} \end{align}\end{split}\]

where \(\mathbf{A}\) is a dense linear operator containing the geometric sensitivity and

\[\frac{\partial \mathbf{b}}{\partial t} = \frac{d \eta}{dt} \mathbf{A \, m}\]

In the case that the predicted data vector \(\mathbb{F} [\mathbf{m}]\) contains field observations at multiple times and locations, \(\eta\) evaluated at different times can be used to construct a sparse matrix \(\mathbf{T}\) such that:

(9)¶\[\mathbb{F} [\mathbf{m}] = \mathbf{T \, A \, m}\]

Although this problem is linear, storing the product of \(\mathbf{T}\) and \(\mathbf{A}\) can become prohibitively large as the number of observation times increases. As a result, the two operators are stored separately and matrix-vector products are taken as needed.

The Jacobian is given by:

\[\mathbf{J} = \mathbf{T \, A}\]

however it is never stored explicitly. Instead, optimization approach need only compute the output of \(\mathbf{Ju}\) and \(\mathbf{J^T v}\) where \(\mathbf{u}\) and \(\mathbf{v}\) are vectors.

Log-Uniform Problem¶

The general forward problem is non-linear, as the VRM does not decay in the same fashion within all cells. In this case however, we are not limited to computing the magnetic response during the off-time. The magnetization experienced by each cells is given by Eq. (6). However since the geometric and time-dependent terms are separable, this can be re-expressed as:

(10)¶\[\vec{M}^k (t) = \vec{h}_{max}^k \, \zeta^k\]

where \(\vec{h}_{max}^k\) is the inducing field at its maximum amplitude and \(\zeta^k\) is the impulse response convolved with the normalized source waveform; \(\zeta\) depends on \(t\) and the values of \(\chi_\infty, \Delta \chi, \tau_1\) and \(\tau_2\) for each cell. Using Eq. (10), \(\vec{\mathbf{M}}(t)\) in Eq. (5) becomes:

(11)¶\[\vec{\mathbf{M}}(t) = \mathbf{H_0 \,} \boldsymbol{ \zeta}(t)\]

where \(\mathbf{H_0}\) is a sparse matrix containing the cartesian components of \(\vec{h}_{max}^k\) for every cell and \(\boldsymbol{\zeta}(t)\) contains the results of a convolution for every cell at time \(t\). Substituting Eq. (11) into Eq. (5):

\[\begin{split}\begin{align} \mathbf{b}(t) &= \mathbf{G \, H_0 \,} \boldsymbol{\zeta} (t) \\ & = \mathbf{A \,} \boldsymbol{\zeta} (t) \end{align}\end{split}\]