TY - JOUR
AU - Antil, H
AB - SUMMARY A growing body of applied mathematics literature in recent years has focused on the application of fractional calculus to problems of anomalous transport. In these analyses, the anomalous transport (of charge, tracers, fluid, etc.) is presumed attributable to long-range correlations of material properties within an inherently complex, and in some cases self-similar, conducting medium. Rather than considering an exquisitely discretized (and computationally intractable) representation of the medium, the complex and spatially correlated heterogeneity is represented through reformulation of the governing equation for the relevant transport physics such that its coefficients are, instead, smooth but paired with fractional-order space derivatives. Here we apply these concepts to the scalar Helmholtz equation and its use in electromagnetic interrogation of Earth’s interior through the magnetotelluric method. We outline a practical algorithm for solving the Helmholtz equation using spectral methods coupled with finite element discretizations. Execution of this algorithm for the magnetotelluric problem reveals several interesting features observable in field data: long-range correlation of the predicted electromagnetic fields; a power-law relationship between the squared impedance amplitude and squared wavenumber whose slope is a function of the fractional exponent within the governing Helmholtz equation; and, a non-constant apparent resistivity spectrum whose variability arises solely from the fractional exponent. In geological settings characterized by self-similarity (e.g. fracture systems; thick and richly textured sedimentary sequences, etc.) we posit that these diagnostics are useful for geological characterization of features far below the typical resolution limit of electromagnetic methods in geophysics. Electrical properties, Electromagnetic theory, Magnetotellurics, Numerical modelling 1 INTRODUCTION Anomalous diffusion has been at the heart of considerable research directed at understanding non-standard, or ‘anomalous’, transport behaviour where the mean squared displacement of random walk particles no longer adheres to a linear relationship with time. As a result, such systems reveal power laws indicative of sub- and superdiffusive behaviour. Anomalous diffusion can be described through random walks endowed by heavy tail distributions and can be captured through non-integer exponents on time and space derivatives. Fractional derivatives in space model superdiffusion and are related to long power-law particle jumps, whereas fractional time derivatives model subdiffusion and reflect long waiting times between particle jumps. Such behaviour has been observed in many applications including reaction-diffusion, quantum kinetics, flow through porous media, plasma transport, magnetic fields, molecular collisions and geophysics applications. We refer the reader to Metzler & Klafter (2000) for a detailed description of anomalous diffusion, including a comprehensive list of applications. The underlying cause for anomalous behaviour in these applications is the presence of complex structures or mechanisms that either promote sub- or superdiffusion. For instance, in fluid flow through porous media, complex permeability fields cause subdiffusive transport through trapping mechanisms or solid/liquid surface-chemistry kinetics, ultimately inducing a memory-type effect (Caputo 2000; Deseri & Zingales 2015). Superdiffusive responses have been experimentally observed in diffusion-reaction system where the variance of a chemical wave exceeds a linear temporal relationship (von Kameke et al. 2010). This is interpreted to be a result of non-local interactions over distances beyond that to the nearest neighbours. In the diffusion-reaction case, vortices in a chaotic velocity field introduce flow paths that exceed standard diffusion rates. A range of natural phenomena can be described as processes in which a physical quantity is constantly undergoing small, random fluctuations. Such Brownian motion can be interpreted as a random walk that, according to the central limit theorem, approaches a normal distribution as the number of steps increases. A macroscopic manifestation of Brownian motion is defined as diffusion whereby a collection of microscopic quantities tends to spread steadily from regions of high concentration to regions of lower concentration, a connection between the micro- and macroscopic first described over a century ago (Einstein 1905). Through Fick’s first and second laws, macroscopic particle movement can be captured by the familiar diffusion equation, the solution of which is a normal distribution corresponding to the random walk probability density. These well known concepts provide the underpinning to investigate phenomena that violate the standard diffusive regime. An application space which has received relatively little attention but is poised for further exploration is low-frequency electromagnetic imaging and interrogation of Earth’s subsurface, a classic geophysical exploration technique premised on diffusive physics, through the mobility of solid-state defects in crystalline materials and free electrons in metals, and electrolytic conduction in fluids (Palacky 1987; Evans 2012; Karato & Wang 2013). Capturing anomalous diffusion in partial differential equations (PDEs) poses considerable mathematical and numerical challenges, particularly in the area of (1) imposing non-zero boundary conditions, (2) validating fractional behaviour for different physics and (3) achieving computational efficiency to realize scalable performance. In representing the fractional Laplacian operator, an integral or a spectral definition can be considered. The choice of method however remains an open question, especially in the presence of non-zero boundary conditions. In this paper, we consider the spectral definition and justify this choice based on the authors’ previous developments. Validating fractional PDEs against field observations, laboratory measurement or analytic solutions is difficult, in part, because fractional calculus development has been somewhat isolated from engineering and science applications community. In this paper, we offer results to help bridge that gap between the observational science and mathematics communities. In particular, fractional concepts are applied to geophysical electromagnetics to better characterize the subsurface and subsequently validated through field observations, as well as geophysical insight. A key challenge in simulating fractional PDEs is achieving computational efficiency. In our work we pursue an approach that leverages the Dunford-type integral representation, which in the case of the fractional Poisson equation, is computationally very attractive because the problem is easily transformed into an uncoupled set of classical (integer order) Poisson solves that can be computed in parallel and then aggregated afterward for the fractional-order solution. This parallelism however does not map to our Helmholtz implementation because of coupling cross-terms, a fact which impacts our eventual desire to solve the geophysical inverse problem by the methods of PDE-constrained optimization (e.g. Antil & Leykekhman 2018) whose adjoint equations require utilization of (the transpose of) the resulting linear system we derive here. Nonetheless, we begin with the solution strategy in Bonito & Pasciak (2015, 2017), with its resultant Helmholtz coupling, and augment it with the lifting (splitting) strategy of Antil et al. (2018b) for handling non-zero boundary conditions (eq. 7 and surrounding text for specific definition) to arrive at a unique approach to solving the fractional Helmholtz equation, regardless of its particular application here in geophysical electromagnetics. The aim of this paper is to explore anomalous diffusion in the context of geophysical electromagnetics and to derive mathematical and algorithmic strategies for practical simulation capabilities. Earth’s subsurface is known to be ripe with complex geological heterogeneity of hierarchical structures and self-similar geometries arising from its generative petrophysical processes and subsequent tectonic and physiochemical experiences. In many cases, low-frequency electromagnetic energy incident upon such media is well-known to produce a measured response incompatible with that of an ‘equivalent’ homogeneous and isotropic medium, or even practical approximations of piecewise-constant material assemblages. For example, non-local electromagnetic effects due to fractures and stratigraphic layering have been analysed, and in some cases observed, in the context of near-surface geotechnical engineering through analysis of spatial statistics (Everett & Weiss 2002; Weiss & Everett 2007; Ge et al. 2015) and transient decay rates (Vallianatos 2017, 2018; Vallianatos et al. 2018), as well as in regional tectonic fabric (Kumar et al. 2018). This perspective on electromagnetic data analysis was hugely influenced by speculative analogy with fluid flow in porous media and the fractional calculus framework developed therein (Benson et al. 2000; Caputo 2000; Deseri & Zingales 2015), and is bolstered by a growing body of observational evidence of electromagnetic phenomena elsewhere (e.g. the response of composite metamaterials in Elser et al. 2007). Ultimately, as geophysicists, we are interested in detecting small scale features in the geological subsurface (e.g. fine-scale stratigraphic laminations or regions permeated by fractures) that may aggregate into a hierarchical ‘meta-material’, but our interest is tempered by practical necessity, to avoid the detailed and computationally explosive discretization required to represent each of them in a given numerical simulation. To mitigate this problem, we propose to dispense with the exquisite discretization just described and replace it by a piecewise blocky medium endowed with a non-local response captured through a fractional space derivative in the material constitutive relation. This is similar to prior approaches for studying anomalous diffusion in hydrology (Caputo 2000) and geophysical electromagnetics (Everett & Weiss 2002; Weiss & Everett 2007; Ge et al. 2015), where, instead, a time-fractional derivative was considered.The challenge here is deriving a solution to the space-fractional Helmholtz equation, a task not obviously suited to the Fourier/Laplace spectral approach readily adaptable for its time-fractional counterpart (Everett 2009; Ge et al. 2015). Our main contributions consist of the following: (1) deriving the fractional Helmholtz equation by introducing a fractional space derivative into Ohm’s law to account for non-local conductivity and recognizing that the solution process for this fractional PDE requires a decomposition to separate boundary conditions on the fractional Laplacian operator; (2) implementing computationally efficient methods through a combination of spectral characterization of the Laplacian, finite element discretization, and a Dunford-type integral representation, followed by a reformulation allowing for sparse, scaled Jacobians with much improved conditioning and (3) validating EM behaviour through exemplar problems in magnetotellurics (MT). The remainder of the paper is organized by first deriving the fractional Helmholtz equation through a fractional gradient relationship between the magnetic field and the underlying electric conductivity. Next the fractional Helmholtz equation with non-zero boundary conditions is decomposed to separate non-zero right-hand side and non-zero boundary conditions. The separation provides a convenient solution strategy and leverages the Dunford-type integral approach to numerically solve one of the remaining equations with a fractional Laplacian. After the mathematical formulation, our finite element implementation is verified through the methods of manufactured solutions. Finally, our numerical capability is demonstrated on a relevant magnetotullerics application. Our numerical results are validated through field measurements and geophysical insight. To supplement the focused ‘Mathematical Formulation’ Section 2 that follows, an extended Appendix is also included so that intermediate mathematical steps, their connection to previous work by Everett (2009) on time-fractional diffusion, and any corresponding physical connections are laid bare. 2 MATHEMATICAL FORMULATION The formulation that follows consists of multiple steps. We start by motivating the fractional derivative operator in the context of geophysical electromagnetics and define Ohm’s constitutive law in terms of a fractional space derivative, later moving that derivative to the Laplacian term of the Helmholtz PDE. Given a fractional Helmholtz equation with non-zero boundary conditions, a two-equation decomposition provides a grouping of boundary conditions, source terms and fractional operators that allow for convenient solution strategies. One equation with non-homogeneous boundary conditions is transformed to non-fractional form by deriving the ‘very weak form’ so that standard solution techniques can be applied. Recall that the ‘weak’ or ‘variational’ form of a governing partial differential equation derives its solution with respect to an arbitrary test function, each of which reside in a function (typically, Sobolev) space with sufficient continuity of the function and its derivatives, provided Dirichlet boundary conditions (if present) also possess similar continuity. The ‘very weak’ form is an extension of this variational problem that allows for ‘rougher’ boundary conditions which reside outside this function space (Berggren 2004). For the remaining equation with a fractional Laplacian and homogeneous boundary conditions, we appeal to a spectral representation and resolvent formalism whereby the fractional Laplacian is transformed to a summation of standard Laplacians using a Dunford-type integral with appropriate quadrature. The solution to the final system of equations is detailed in Section 2, in which a finite element discretization of both equations results in a large and dense coefficient matrix that requires further manipulation to achieve efficient solution performance. An expanded derivation of the fractional Helmholtz equation we consider in the remainder of this paper is presented in the Appendix, along with reference to the relevant motivating prior work in electromagnetic geophysics. The details and intermediate steps of the derivation therein are only briefly mentioned below so that the primary, original mathematical development—a numerical solution to the fractional Helmholtz equations—may commence without too much distracting front-matter. To summarize the appendix, we start with Faraday’s law in the frequency domain with the Fourier convention of time derivatives ∂t mapping to the frequency domain ω as ∂t↦iω and assuming constant magnetic permeability μ0 = 4π × 10−7 H/m: $$\begin{eqnarray*} \nabla \times {\bf E} = -\text{i}\omega \mu _0{\bf H}, \end{eqnarray*}$$(1) relating the curl of electric field E to time variations in magnetic field |$\bf H$|⁠. Paired with this is Ampère’s Law, ∇ × H = J, where J is the total electric current density—the sum of Ohmic currents, Maxwell’s displacement current iωεE, and, any impressed external currents due to natural sources or engineered antennas. Typically, for simple linear, isotropic materials the Ohmic currents are described by the product of electrical conductivity σ and electric field and at sufficiently low frequencies σ > >ωε Maxwell’s displacement current can safely be ignored (MHz and below for geophysical applications). In a similar fashion as Caputo (2000), where he replaced the permeability in Darcy’s equation with a time-fractional derivative, we replace the simple conductivity/field product in Ohm’s law with a space-fractional derivative.Preserving the symmetry of both positive and negative power law jumps in the z direction for a two-sided stable diffusion process requires both positive and negative fractional derivatives (Meerschaert & Sikorskii 2012, p. 15) |$\mathcal {D}^\alpha _z = {1\over 2\cos ({1\over 2}\pi \alpha )}{\left( {\partial ^\alpha \over {\partial z^\alpha }} + {\partial ^\alpha \over \partial (-z)^\alpha } \right)}$| which we normalize by the factor |$\cos ({1\over 2}\pi \alpha )$| to preserve magnitude invariance under α. As a consequence, the space-Fourier transform for this operator maps to the wavenumber ν domain as |$D_z^\alpha \mapsto |\nu |^\alpha$|⁠. Note that had a one-sided derivative with Fourier mapping ( ± iν)α been used, the unit-magnitude prefactor (± i)α could be interpreted as rotating the electrical conductivity into the complex plane, effectively reintroducing the Maxwell displacement current and turning the Maxwell derivation into a mixed diffusion/wave propagation problem rather than a strictly diffusive one. Unlike Caputo’s attempt to emulate memory effects in permeability, our hypothesis is that certain non-local conductivity properties create superdiffusive behaviour and can be represented by a spatial fractional derivative. With this fractional Ohm’s law in the low-frequency limit (and no external sources) inserted into Ampère’s law, the curl of (1) is thereby: $$\begin{eqnarray*} \nabla \times \nabla \times {\bf E} = -\text{i}\omega \mu _0\mathcal {D}_z^\alpha \left[ \sigma _{\alpha ,z}{\bf E}\right], \end{eqnarray*}$$(2) where |$\mathcal {D}_z^\alpha$| is the α-order fractional derivative in the z direction and σα, z is the electrical conductivity in units of S m1−α. For an Earth model whose conductivity varies only as a function of vertical coordinate z, subject to a vertically incident electric field oriented in the horizontal x direction, the electric fields in the Earth are everywhere horizontal such that |${\bf E}=\hat{x} \, u(z)$| is the primary state variable that needs further consideration. Furthermore, for dimensional consistency in the fractional calculus methodology described in the following section, we non-dimensionalize with respect to the z coordinate such that z↦ζ = z/z* to arrive at $$\begin{eqnarray*} -{\text{d}^2u\over \text{d}\zeta ^2}\left( {1\over z^*}\right)^2 + \text{i}\omega \mu _0\mathcal {D}_\zeta ^\alpha \left[ \sigma _{\alpha ,\zeta }\, u(\zeta ) \right]= 0 \end{eqnarray*}$$(3) which after action by |$\mathcal {D}_\zeta ^{-\alpha }$| and generalization to 3-D, becomes $$\begin{eqnarray*} \left(-\Delta _\zeta \right)^s u + \text{i}\kappa ^2\, u(\zeta ) = 0, \end{eqnarray*}$$(4) where (− Δζ)s is the fractional-order Laplacian in dimensionless coordinate ζ, |$s = 1-{1\over 2}\alpha$| and κ2 the dimensionless squared wavenumber ωμ0σα, ζ(z*)2. Note that in eq. (3) the conductivity σα, z possesses fractional length dimensions to retain consistency with the fractional derivative operator |$\mathcal {D}_z^\alpha$|⁠. However, through the non-dimensionalization process transforming eq. (3) to eq. (4) we see that the conductivity σα, ζ reclaims its familiar, integer-ordered units of S m–1, thus avoiding awkward, fractional-dimensioned conductivities reported elsewhere (Everett 2009; Ge et al. 2015). We observe that the fractional exponent is on the Laplacian term and in combination with the Helmholtz term motivate the challenge of a solution strategy. An additional complication is the need to incorporate non-trivial boundary conditions, such as special radiation or self-absorbing boundary conditions. We address these issues through the use of linear decomposition, a Dunford type integral formulation and the very weak form for finite element discretizations. We write the generalized fractional-order Helmholtz equation as $$\begin{eqnarray*} \begin{aligned} \left(-\Delta \right)^{s}u -k^2 u & = f \quad \mbox{in}\, \, \Omega , \\ u & = g \quad \mbox{on}\, \, \Gamma , \end{aligned} \end{eqnarray*}$$(5) where Ω is the spatial domain with boundary Γ over which u is solved. The symbol k2 = −iκ2 is introduced to simplify notation for our electromagnetic problem, but in fact, transcends this particular choice of physics, and (−Δ)s is understood to be the spectral fractional Laplacian operator for non-zero Dirichlet boundary conditions $$\begin{eqnarray*} \left( -\Delta \right)^s u(\boldsymbol {x}):= \sum _{k=1}^{\infty } \left( \lambda _k^s\int _\Omega u\, \varphi _k \, {\text{d}}\Omega + \lambda _k^{s-1}\int _\Gamma u\, \partial _n\varphi _k\, {\text{d}}\Gamma \right) \varphi _k(\boldsymbol {x}), \end{eqnarray*}$$(6) where φk are the eigenfunctions of the Laplacian with corresponding eigenvalues λk (Antil et al. 2018b, def. 2.3) and |$\boldsymbol {x}\in \Omega$| is the coordinate of interest. Moreover, u is assumed to be sufficiently smooth. As point of mathematical rigor, we have stated this definition for smooth functions, as is customary in PDE theory, but point out that by standard density arguments it immediately extends to Sobolev spaces and we refer to Antil et al. (2018b) for details. Recall that In addition, we emphasize that when u = 0 on the boundary Γ, the definition above is simply the standard spectral fractional Laplacian (− Δ0)s with zero boundary conditions [i.e. no surface integral terms in (6)]. We will omit the subscript 0 when it is clear from the context. Following this earlier work on fractional Poisson equation, we extend the basic approach to fractional Helmholtz and split u into two parts (a.k.a. ‘lift’) thusly: Let v solve $$\begin{eqnarray*} \begin{aligned} \left(-\Delta \right)^{s}v -k^2 \left(v + w\right) & = f \quad \mbox{in}\, \, \Omega , \\ v & = 0 \quad \mbox{on}\, \, \Gamma , \end{aligned} \end{eqnarray*}$$(7) and let w solve $$\begin{eqnarray*} \begin{aligned} \left(-\Delta \right)^{s}w &= 0 \quad \mbox{in}\, \, \Omega , \\ w & = g \quad \mbox{on}\, \, \Gamma . \end{aligned} \end{eqnarray*}$$(8) Summing (7) and (8) it is evident that u = v + w. This splitting is motivated by the presence of the homogeneous boundary condition on (7) which enables use of the spectral representation of the fractional-power Laplacian operator (6). Furthermore, it has been shown Antil et al. (2018b) that solving (8) is equivalent to solving the standard, integer-power Laplacian equation in the very-weak sense (cf. Berggren 2004; Lions & Magenes 1972), which in the case of smooth g is simply the more-familiar weak sense. A simple algebraic manipulation of the spectral decomposition results in a Laplacian with integer exponents [see theorem 4.1 and the subsequent proof in Antil et al. (2018b) for additional details]. Hence, we may replace (8) with the following: $$\begin{eqnarray*} \begin{aligned} \left(-\Delta \right)w &= 0 \quad \mbox{in}\, \, \Omega , \\ w & = g \quad \mbox{on}\, \, \Gamma , \end{aligned} \end{eqnarray*}$$(9) which, when solved simultaneously with (7), yields the solution to the original eq. (5). To solve eq. (7), we follow others (e.g. Bonito & Pasciak 2015; Antil & Pfefferer 2017) in using spectral analysis of linear operators and resolvent formalism. Specifically, we start with the Kato (1960) definition of fractional powers for linear operators A, which states $$\begin{eqnarray*} \left(\lambda + A^s \right)^{-1} = {\sin s\pi \over \pi } \int _0^\infty {r^s \over \lambda ^2 + 2\lambda r^s \cos s \pi + r^{2s}}\left( r + A\right)^{-1}\, {\text{d}}r, \end{eqnarray*}$$ where for the present analysis, the linear operator A is the Laplacian −Δ. This general definition is simplified by setting Kato’s λ coefficient to zero. This definition due to Kato coincides with (6) when the function values are zero on the boundary, as in this case the surface integral over Γ vanishes which is indeed the case in (7). We remark that the Kato definition above derives from a more general starting point, the Dunford integral representation, but only emerges after significant derivation and deep reliance on theoretical functional analysis which we feel would distract from the present work. Nonetheless, we adopt the phrase ‘Dunford-type’ here because of its accuracy and its common use in the applied mathematics literature on fractional calculus. Curious readers are invited to consult Yosida (1995), chapter VIII–IX, for details on the connection between the Dunford integral, the necessary definition of its resolvant, and the ultimate connection to Kato (1960) by way of Balakrishnan (1960), among others. After a change of variable, the Kato definition results in a symmetric integral, which is approximated through quadrature as $$\begin{eqnarray*} \begin{aligned} \left(-\Delta \right)^{-s} &= {\sin s\pi \over \pi } \int _{-\infty }^\infty e^{(1-s)y} \left( e^y - \Delta \right)^{-1}\, {\text{d}}y, \\ & \approx {\sin s\pi \over \pi } m \sum _{\ell =-N^-}^{N^+} e^{(1-s)y_\ell } \left( e^{y_\ell } - \Delta \right)^{-1}, \end{aligned} \end{eqnarray*}$$(10) where the quadrature nodes are distributed uniformly as yl = mℓ. Accuracy of the quadrature representation of this continuous integral is a function of the constants m, N− and N+ and has been shown to be exponentially convergent (Bonito & Pasciak 2015). The constants are chosen such that the quadrature error is balanced with the error in the spatial discretization of (7) (cf. w = 0 case in Bonito & Pasciak 2015). In the case of a finite element solution with linear nodal basis functions on the unit interval and node spacing h, they are $$\begin{eqnarray*} m = {1\over \ln {1\over h}} , \quad N^+ = \Big \lceil {\pi ^2\over 4s\, m^2}\Big \rceil , \quad \text{and}\quad N^- = \Big \lceil {\pi ^2\over 4(1-s)\, m^2}\Big \rceil . \end{eqnarray*}$$(11) The use of the ‘ceiling’ operators ⌈ · ⌉ in (11) (the smallest integer greater than its argument) ensure N−N+ are integer valued, as required. We remark that in using (10), we avoid the costly (and in many cases, inaccurate) precalculation of the eigenspectrum for the Laplacian over an arbitrary spatial domain Ω with Dirichlet condition |$u\vert _\Gamma = g$|⁠. Even in cases where calculation of the eigenspectrum bears an acceptable computational cost, there still remains the outstanding question of just how much of the spectrum is required for computing the Laplacian by this method to acceptable accuracy. For these reasons, our eq. (10) is far more practical. Rewriting (7) as |$v - \left(-\Delta \right)^{-s}\, k^2 (v+w) = \left(-\Delta \right)^{-s}f$|⁠, we may write v as a Kato-style expansion $$\begin{eqnarray*} v = {\sin s\pi \over \pi } m \sum _{\ell =-N^-}^{N^+} e^{(1-s)y_\ell } v_\ell \end{eqnarray*}$$(12) and equate each of the ℓ terms to arrive at the coupled equation $$\begin{eqnarray*} v_\ell - \left( e^{y_\ell }-\Delta \right)^{-1}\, k^2\left(v + w\right) = \left( e^{y_\ell } - \Delta \right)^{-1}f, \end{eqnarray*}$$(13) or $$\begin{eqnarray*} \left( e^{y_\ell } - \Delta \right) v_\ell - k^2\left(v + w\right) = f. \end{eqnarray*}$$(14) Observe that there are N− + N+ + 1 of these equations and that the ℓth equation fully couples the function vℓ into w and all remaining functions |$v_{\ell ^\prime \ne \ell }$| of the expansion (12). Enforcement of the modified boundary condition vℓ = 0 guarantees enforcement of v = 0 via (12). Hence, with the inclusion of (9), we have the complete differential problem statement consisting of a coupled system of N− + N+ + 2 equations with unknown functions |$v_{N^-},\ldots ,v_{N^+}, w$| over the domain Ω. 3 NUMERICAL IMPLEMENTATION The method of solution for eqs (9) and (14) (including corresponding boundary conditions) is to first transform the differential problem statement into an equivalent variational problem statement for the appropriate infinite dimensional function spaces and then approximate its solution by the optimal one, in a Sobolev norm sense, taken from finite dimensional space of linear, nodal finite elements over some discretization. In doing so, we introduce the test function ξℓ and construct the weak form of (14) for all ℓ from −N− to N+: $$\begin{eqnarray*} \int _\Omega \left(-\xi _\ell \, \Delta v_\ell + e^{y_\ell }\, \xi _\ell v_\ell - k^2\, \xi _\ell \left( v + w\right) \right)\, {\text{d}}\Omega = \int _\Omega \xi _\ell f\, {\text{d}}\Omega , \end{eqnarray*}$$(15) recalling that v is given by the expansion (12). The test function ζ is used in the weak form of the Laplace eq. (9) as, $$\begin{eqnarray*} \int _\Omega \nabla \zeta \cdot \nabla w\, {\text{d}}\Omega = 0. \end{eqnarray*}$$(16) Combining the left-hand sides of (15) and (16), the bilinear form A( ·, ·) is therefore given as $$\begin{eqnarray*} A(\left\lbrace \xi _\ell \right\rbrace ,\zeta ;\left\lbrace v_\ell \right\rbrace ,w) = \sum _{\ell =-N^-}^{N^+} \int _\Omega \left(\nabla \xi _\ell \cdot \nabla v_\ell + e^{y_\ell }\, \xi _\ell v_\ell - k^2\, \xi _\ell \left( v + w \right) \right) \, {\text{d}}\Omega + \int _\Omega \nabla \zeta \cdot \nabla w\, {\text{d}}\Omega , \end{eqnarray*}$$(17) the combined right-hand sides are denoted as $$\begin{eqnarray*} F\left( \left\lbrace \xi _\ell \right\rbrace ,\zeta \right) = \sum _{\ell =-N^-}^{N^+} \int _\Omega \xi _\ell f\, {\text{d}}\Omega , \end{eqnarray*}$$(18) and v is understood to be expanded in terms of vℓ according to (12). The variational problem statement equivalent to the differential problem statement in (9) and (14) is therefore: Find {vℓ} ∈ V0; w ∈ V, such that $$\begin{eqnarray*} A(\left\lbrace \xi _\ell \right\rbrace ,\zeta ;\left\lbrace v_\ell \right\rbrace ,w) = F\left( \left\lbrace \xi _\ell \right\rbrace ,\zeta \right)\quad \forall \left\lbrace \xi _\ell \right\rbrace ,\zeta \in V_0, \end{eqnarray*}$$(19) where V is the space of L2 functions on Ω with first order weak derivatives also in L2 of Ω and inhomogeneous Dirichlet boundary conditions (8), and V0⊂V taking homogeneous boundary conditions as in (7). The next step is to choose a finite-dimensional space Vh⊂V from which the approximate solutions vh ≈ v and wh ≈ w will be drawn. To further simplify notation we will drop the h subscript and only re-introduce it as needed in relation to the true (weak) solutions v and w. Let |$\phi _1(\boldsymbol {x}), \phi _2(\boldsymbol {x}), \ldots , \phi _N(\boldsymbol {x})$| be the basis functions in Vh as a function of spatial coordinate |$\boldsymbol x$|⁠, which in our implementation will be linear, nodal finite elements . We write in bold the column vector |$\boldsymbol {\phi }$| of basis functions (ϕ1, ϕ2, …, ϕN)T, and |$\boldsymbol {v}_\ell =\left( v_\ell ^1, v_\ell ^2,\ldots ,v_\ell ^N \right)^T$|⁠. By construction, it follows that |$v_\ell (\boldsymbol {x}) = \boldsymbol {\phi }^T\boldsymbol {v}_\ell = \sum _{i=1}^N \phi _i (\boldsymbol {x}) v_\ell ^i$|⁠. Likewise, |$\boldsymbol {\xi }_\ell = \left(\xi _\ell ^1, \xi _\ell ^2,\ldots ,\xi _\ell ^N \right)^T$| and |$\boldsymbol {w}=\left( w_1, w_2,\ldots ,w_N \right)^T$| yield |$\xi _\ell = \boldsymbol {\xi }_\ell ^T\boldsymbol {\phi }$| and |$w = \boldsymbol {\phi }^T\boldsymbol {w}$|⁠, respectively. As such, construction of the linear system (19) requires the following block matrices built by volume integration of the basis functions and their spatial derivatives: $$\begin{eqnarray*} {{\bf K}}= \int _\Omega \left( \nabla \boldsymbol {\phi } \right)^T \left( \nabla \boldsymbol {\phi }\right) \, {\text{d}}\Omega , \quad {{\bf M}}_1 = \int _\Omega \boldsymbol {\phi }\boldsymbol {\phi }^T\, {\text{d}}\Omega , \quad \mbox{and}\quad {{\bf M}}_2 = -\int _\Omega k^2\, \boldsymbol {\phi }\boldsymbol {\phi }^T\, {\text{d}}\Omega , \end{eqnarray*}$$(20) where the integrands are understood to be outer products, each yielding a symmetric matrix of dimension N × N. To simplify notation, introduce the coefficients |$c_\ell = e^{y_{\ell -N^-}}$| and |$d_\ell = {1\over \pi } (\sin s\pi )\, m \, e^{(1-s)y_{\ell -N^-}}$|⁠, matrix Aℓ = K + cℓM1 + dℓM2, and sum L = N− + N+ so that we may compactly write A({ξℓ}, ζ; {vℓ}, w) as $$\begin{eqnarray*} \left( \boldsymbol {\xi }^T_{-N^-} \cdots \boldsymbol {\xi }^T_{N^+} \, \boldsymbol {\zeta }^T \right) \begin{pmatrix}{\bf A}_0 & \quad d_1\, {{\bf M}}_2 & \quad d_2\, {{\bf M}}_2 & \quad \cdots & \quad d_L{{\bf M}}_2 & \quad {{\bf M}}_2 \\ d_0\, {{\bf M}}_2 & \quad {{\bf A}}_1 & \quad d_2\, {{\bf M}}_2 & \quad \cdots & \quad d_L{{\bf M}}_2 & \quad {{\bf M}}_2 \\ d_0\, {{\bf M}}_2 & \quad d_1\, {{\bf M}}_2 & \quad {{\bf A}}_2 & \quad \cdots & \quad d_L{{\bf M}}_2 & \quad {{\bf M}}_2 \\ \vdots & \quad \vdots & \quad \vdots & \quad \ddots & \quad \vdots & \quad \vdots \\ d_0\, {{\bf M}}_2 & \quad d_1\, {{\bf M}}_2 & \quad d_2\, {{\bf M}}_2 & \quad \cdots & \quad {{\bf A}}_L& \quad {{\bf M}}_2 \\ {\bf 0} & \quad {\bf 0} & \quad {\bf 0} & \quad \cdots & \quad {\bf 0}& \quad {{\bf K}} \\ \end{pmatrix} \begin{pmatrix}\boldsymbol {v}_{-N^-}\\ \boldsymbol {v}_{1-N^-}\\ \boldsymbol {v}_{2-N^-}\\ \vdots \\ \boldsymbol {v}_{N^+}\\ \boldsymbol {w}\\ \end{pmatrix}. \end{eqnarray*}$$(21) Lastly, the right-hand side of (19) follows: $$\begin{eqnarray*} \left( \boldsymbol {\xi }^T_{-N^-} \cdots \boldsymbol {\xi }^T_{N^+} \, \boldsymbol {\zeta }^T \right) \begin{pmatrix}\boldsymbol {f}\\ \boldsymbol {f}\\ \boldsymbol {f}\\ \vdots \\ \boldsymbol {f}\\ \boldsymbol {0} \\ \end{pmatrix} \end{eqnarray*}$$(22) with column vector |$\boldsymbol {f} = \int _\Omega \boldsymbol {\phi } f\, {\text{d}}\Omega$|⁠. In equating (21) with (22) as required by the variational problem statement (19), we see that the coefficient vector |$\left( \boldsymbol {\xi }^T_{-N^-} \cdots \boldsymbol {\xi }^T_{N^+} \, \boldsymbol {\zeta }^T \right)$| is common to both the left- and right-hand sides, and may therefore be divided out, thus leaving a N(L + 2) × N(L + 2) system of linear equations for the unknown coefficients |$\left( \boldsymbol {v}^T_{-N^-} \cdots \boldsymbol {v}^T_{N^+} \, \boldsymbol {w}^T \right)$| which holds for all functions {ξℓ}, ζ ∈ Vh. Upon solution of the linear system, aggregation of the coefficient vectors |$\boldsymbol {v}_\ell$| according to (12) plus the vector |$\boldsymbol {w}$| completes the sum v + w, which we recognize as the discrete, approximate solution to the original differential eq. (5). Because the matrix in (21) is complex-valued, large and non-symmetric, the solution strategy for the linear system equating (21) and (22) must be carefully chosen for scalability and economy of compute resources. As such, we solve this linear system using stabilized bi-conjugate gradients (van der Vorst 1992): The algorithm is easily parallelizable; has a minimum number of working vectors; and requires only two matrix–vector products per iterative step. The latter is especially important for reducing computational resource burdens because these products can be computed cheaply and quickly ‘on the fly’ as needed and without explicit storage of entire system matrix. Note, however, that the matrix in (17) is block dense and a large number of floating point operations is required for a single matrix–vector multiply—operations which may significantly increase the time required to perform the multiplications. To remedy this, we modify the variational formulation (19) to include the function v in addition to the vectors vℓ and w, and augment A( · ) by weak enforcement of the compatibility expansion (12) between vectors |$\boldsymbol {v}_\ell$| and |$\boldsymbol {v}$|⁠.That is, in addition to {vℓ} and w, we introduce the additional unknown v and find {vℓ}, w and v such that $$\begin{eqnarray*} {\skew{4}\tilde{A}}(\left\lbrace \xi _\ell \right\rbrace ,\zeta ,\eta ;\left\lbrace v_\ell \right\rbrace ,w,v) = {\skew{4}\tilde{F}}\left( \left\lbrace \xi _\ell \right\rbrace ,\zeta ,\eta \right)\quad \forall \left\lbrace \xi _\ell \right\rbrace ,\zeta , \text{ and } \eta , \end{eqnarray*}$$(23) where $$\begin{eqnarray*} {\skew{4}\tilde{A}}(\left\lbrace \xi _\ell \right\rbrace ,\zeta ,\eta ; \left\lbrace v_\ell \right\rbrace ,w,v) &=& \sum _{\ell =N^-}^{N^+} \int _\Omega \left(\nabla \xi _\ell \cdot \nabla v_\ell + e^{y_\ell }\, \xi _\ell v_\ell - k^2\, \xi _\ell \left( v + w \right) + \nabla \zeta \cdot \nabla w \right) \, {\text{d}}\Omega \nonumber \\ && - {\sin s\pi \over \pi } m \sum _{\ell =-N^-}^{N^+} e^{(1-s)y_\ell } \int _\Omega \eta \, v_\ell \, {\text{d}}\Omega + \int _\Omega \eta \, v\, {\text{d}}\Omega \end{eqnarray*}$$(24) and $$\begin{eqnarray*} {\skew{4}\tilde{F}}\left( \left\lbrace \xi _\ell \right\rbrace ,\zeta ,\eta \right) = \sum _{\ell =N^-}^{N^+} \int _\Omega \xi _\ell f\, {\text{d}}\Omega . \end{eqnarray*}$$(25) The resulting sparse linear system is thus, $$\begin{eqnarray*} \begin{pmatrix}\tilde{{\bf A}}_0& \quad {\bf 0}& \quad \cdots & \quad \cdots & \quad {\bf 0} & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ {\bf 0} & \quad \tilde{{\bf A}}_1 & \quad {\bf 0} & \quad \cdots & \quad {\bf 0} & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ \vdots & \quad & \quad \ddots & \quad & \quad \vdots & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ {\bf 0} & \quad \cdots & \quad {\bf 0} & \quad \tilde{{\bf A}}_{L-1} & \quad {\bf 0} & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ {\bf 0} & \quad \cdots & \quad \cdots & \quad {\bf 0} & \quad \tilde{{\bf A}}_L & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ {\bf 0} & \quad \cdots & \quad \cdots & \quad \cdots & \quad {\bf 0} & \quad {{\bf K}} & \quad {\bf 0} \\ -d_0{{\bf M}}_1 & \quad \cdots & \quad \cdots & \quad \cdots & \quad -d_L{{\bf M}}_1 & \quad {\bf 0} & \quad {{\bf M}}_1 \\ \end{pmatrix} \begin{pmatrix}\boldsymbol {v}_{-N^-}\\ \boldsymbol {v}_{1-N^-}\\ \boldsymbol {v}_{2-N^-}\\ \vdots \\ \boldsymbol {v}_{N^+}\\ \boldsymbol {w}\\ \boldsymbol {v}\\ \end{pmatrix} = \begin{pmatrix}\boldsymbol {f}\\ \boldsymbol {f}\\ \boldsymbol {f}\\ \vdots \\ \boldsymbol {f}\\ \boldsymbol {0}\\ \boldsymbol {0}\\ \end{pmatrix} \\ \end{eqnarray*}$$(26) with |$\tilde{{\bf A}}_\ell = {{\bf K}} + c_\ell {{\bf M}}_1$|⁠. We refer to the sparsified system of linear equations in (26) as the v-formulation for the discretized, variational form of original differential eq. (5). Inspection of the pre-factors dℓ from the Kato expansion, however, suggests that their exponential decay with respect to ℓ may lead to ill-conditioning of the coefficient matrix in (26) by their presence in the last block-row. As such, we recast (26) as the scaled v-formulation: $$\begin{eqnarray*} \begin{pmatrix}\tilde{{\bf A}}^\prime _0& \quad {\bf 0}& \quad \cdots & \quad \cdots & \quad {\bf 0} & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ {\bf 0} & \quad \tilde{{\bf A}}^\prime _1 & \quad {\bf 0} & \quad \cdots & \quad {\bf 0} & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ \vdots & \quad & \quad \ddots & \quad & \quad \vdots & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ {\bf 0} & \quad \cdots & \quad {\bf 0} & \quad \tilde{{\bf A}}^\prime _{L-1} & \quad {\bf 0} & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ {\bf 0} & \quad \cdots & \quad \cdots & \quad {\bf 0} & \quad \tilde{{\bf A}}^\prime _L & \quad {{\bf M}}_2 & \quad {{\bf M}}_2 \\ {\bf 0} & \quad \cdots & \quad \cdots & \quad \cdots & \quad {\bf 0} & \quad {{\bf K}} & \quad {\bf 0} \\ -{{\bf M}}_1 & \quad \cdots & \quad \cdots & \quad \cdots & \quad -{{\bf M}}_1 & \quad {\bf 0} & \quad {{\bf M}}_1 \\ \end{pmatrix} \begin{pmatrix}\boldsymbol {v}^\prime _{-N^-}\\ \boldsymbol {v}^\prime _{1-N^-}\\ \boldsymbol {v}^\prime _{2-N^-}\\ \vdots \\ \boldsymbol {v}^\prime _{N^+}\\ \boldsymbol {w}\\ \boldsymbol {v}\\ \end{pmatrix} = \begin{pmatrix}\boldsymbol {f}\\ \boldsymbol {f}\\ \boldsymbol {f}\\ \vdots \\ \boldsymbol {f}\\ \boldsymbol {0}\\ \boldsymbol {0}\\ \end{pmatrix} \\ \end{eqnarray*}$$(27) with |$\tilde{{\bf A}}^\prime _\ell = \left( {{\bf K}} + c_\ell {{\bf M}}_1\right)/d_\ell$| and |$\boldsymbol {v}^\prime _\ell = \boldsymbol {v}_\ell d_{\ell + N^-}$|⁠. In the scaled system (27), the Kato scale factors dℓ are implicit in the unknown vectors |$\boldsymbol {v}^\prime _\ell$| and act within block diagonal matrices |$\skew{4}\tilde{\boldsymbol A}^\prime _\ell$| alone. As a closing remark on the theory and algorithms just described, observe that in (9), (26) and (27) the solution of the scalar Laplacian equation for w is fully decoupled from the solution for v and vℓ. Hence, one option for solving the full system of equations is a two-step procedure, where first the solution for w is obtained, and then used as a sourcing term for the remaining vℓ equations. We have, instead, chosen to solve the full system simultaneously. This has some advantages. First, in looking ahead to the implementation of Robin boundary conditions (e.g. a Sommerfeld radiation condition), we anticipate that the Laplacian equations will couple directly into the v (or, equivalently, vℓ) equations, which would consequently eliminate the convenience of solving for w a priori. We wish this coupling to modify our existing algorithm/code structure as little as possible and therefore retain the Laplacian equation for w in the full system matrix. Secondly, including the Laplacian comes at an increased cost of only N degrees of freedom on top of the existing cost of |$L\, N$| for the vℓ equations. Because L is typically on the order of 100 or more for adequately refined meshes (Fig. 1), this added cost is objectively minimal. Firstly, looking further ahead towards PDE constrained optimization where we might invert for s or σα, ζ, or for the design problem (e.g. optimal sensor sensor placement), it is more convenient to create and solve for the corresponding adjoint objects. Figure 1. Open in new tabDownload slide Total number of quadrature points L as a function of Laplacian exponent s and node spacing h. Sinc quadrature summation is over the range of indices ℓ = −N−, …, N+. 3.1 Non-locality of the fractional Laplacian operator We next provide insight as to why eqs (10), (6) are non-local operators. As pointed out in Song & Vondraček (2008), Caffarelli & Stinga (2016), Antil et al. (2017) and Antil et al.(2019), the spectral fractional Laplacian with zero boundary conditions can be equivalently written as: $$\begin{eqnarray*} (-\Delta _0)^s u(\boldsymbol {x}) = \int _\Omega ( u(\boldsymbol {x}) - u(\boldsymbol {x}^\prime ) ) \mathcal {K}_s(\boldsymbol {x},\boldsymbol {x}^\prime )\, \text{d}\Omega ^\prime + u(\boldsymbol {x}) B_s(\boldsymbol {x}) , \end{eqnarray*}$$(28) where |$\mathcal {K}_s$| and Bs are appropriate Kernel functions. One can see from the equivalent definition that, in order to evaluate ( − Δ0)su at a point |$\boldsymbol {x} \in \Omega$|⁠, we need information about u on the entire domain Ω, thus making ( − Δ0)s a nonlocal operator. Moreover, let |$\mathcal {O}$| be an open set contained in Ω such that u ≡ 0 on |$\mathcal {O}$|⁠. Then a classical Laplacian implies |$-\Delta u(\boldsymbol {x}) = 0$| for all |$\boldsymbol {x} \in \mathcal {O}$|⁠. However, this is not the case when we deal with ( − Δ0)s. Indeed, let |$\boldsymbol {x} \in \mathcal {O}$| and use the equivalent definition (28). Since u ≡ 0 on |$\mathcal {O}$|⁠, we obtain $$\begin{eqnarray*} (-\Delta _0)^s u(\boldsymbol {x}) &=& \int _\Omega \left( u(\boldsymbol {x}) - u(\boldsymbol {x}^\prime ) \right) \mathcal {K}_s(\boldsymbol {x},\boldsymbol {x}^\prime ) {\text{d}}\Omega ^\prime \nonumber \\ &=& \int _{\mathcal {O}} \left( u(\boldsymbol {x}) - u(\boldsymbol {x}^\prime ) \right) \mathcal {K}_s(\boldsymbol {x},\boldsymbol {x}^\prime )\, {\text{d}}\Omega ^\prime + \int _{\Omega \setminus \mathcal {O}} \left( - u(\boldsymbol {x}^\prime ) \right) \mathcal {K}_s(\boldsymbol {x},\boldsymbol {x}^\prime ) \, {\text{d}}\Omega ^\prime \nonumber \\ &=& \int _{\Omega \setminus \mathcal {O}} \left( - u(\boldsymbol {x}^\prime ) \right) \mathcal {K}_s(\boldsymbol {x},\boldsymbol {x}^\prime ) \, {\text{d}}\Omega ^\prime \end{eqnarray*}$$(29) which is not necessarily zero. This is unlike the local case. 4 NUMERICAL VERIFICATION To verify the implementation of the fractional Helmholtz equations with inhomogeneous Dirichlet boundary conditions we adopt the method of manufactured solutions (MMS, Roache 1998; Salar & Knupp 2000). In the MMS method, a proposed solution is substituted into the governing differential equation, after which the corresponding boundary conditions and sourcing functions are derived. Upon discretization, the recovered numerical solution is then compared to the known analytical solution. The MMS solution used here over the interval x ∈ [0, 1] takes the following form: v = sin (2πx), w = 1↦u = 1 + sin (2πx). Inspection of (9) results in setting g(0) = g(1) = 1, whereas recognizing that sin (2πx) is an eigenfunction for homogeneous ( − Δ)s results in f = ((2π)s − k2)sin (2πx) − k2 according to (7). Hence, we have constructed for arbitrary s the requisite source terms and boundary conditions for our posited solution u solving an inhomogeneous fractional Helmholtz with non-zero Dirichlet boundary conditions. Note that MMS solution u is independent of both s and k and is thus powerful test of fractional Helmholtz algorithm. Numerical evaluation of the MMS problem just described is done using linear, nodal finite elements with uniform node spacing for the case of s = 0.25 and k arbitrarily set to unity. Linear system (27) is solved to high accuracy using the stabilized bi-conjugate gradient van der Vorst (1992) iterative scheme with simple Jacobi scaling to a tolerance of 10−16 reduction in normalized residual. Over the range of node spacing 0.001 ≤ h ≤ 0.01 the MMS solution shows the expected h2 convergence in error between the recovered finite element and known analytic solutions (Fig. 2) in L2(Ω)-norm. For reference, the size Ntotal of the linear system (27) grows roughly as N1.4 over the corresponding range in h, resulting, for example, in L = 629 quadrature points for N = 1001 finite element nodes and a total of Ntotal = 631 631 unknowns in the linear system (27). Convergence of the bi-conjugate gradient residual error as a function of iteration count (Fig. 3) is generally well behaved, with only minor localized excursions from monotonicity. Furthermore, the error in simultaneously solving each of the three sets of coupled equations—fractional Helmholtz for vℓ; Laplacian for w; and, compatibility between v and vℓ—decreases synchronously with iteration count, with error for the compatibility equation approximately a factor 100 less than the error for the remaining two. Figure 2. Open in new tabDownload slide Convergence of MMS solution u = sin (2πx) + 1 for s = 0.25 as a function of mesh size N with corresponding node spacing (N − 1)−1. In symbols is the RMS residual; black lines, the curve h2; and in red, the total number of degrees of freedom in the discretized linear system (27). Linear system (27) is solved using BiCG-STAB to a tolerance of 10−16 in RHS-normalized residual. Figure 3. Open in new tabDownload slide Convergence of BiCG-STAB algorithm with Jacobi pre-conditioning for the N = 101 MMS solution in Fig. 1. In the solid line is the average residual for equations in blocks 0, …, L in (27) corresponding to fractional Helmholtz equations on vℓ; dashed, the residual for equations in block L + 1 corresponding to solution of the Laplace equation for w; and dotted, the residual for the final (L + 2) block of equations enforcing compatibility between v and vℓ. Finally, we confirm that the choice |$m=1/\log {1\over h}$| for quadrature spacing (and by extension, the number L of vℓ equations) is nearly optimal by examining the effect on RMS of varying m. As a representative example, the N = 101, s = 0.25 discretization of the MMS problem is solved for a range of m values around |$1/\log {1\over h}$|⁠. In this example the asymptotic limit for minimum RMS value is achieved at m roughly 90 per cent of its optimal value, where the asymptotic limit is driven by the error of the finite element discretization itself (Fig. 4). In contrast, choices of m larger than the optimal value result in a rapidly increasing RMS, consistent with the exponential convergence of quadrature error reported elsewhere (Bonito & Pasciak 2015). Figure 4. Open in new tabDownload slide Convergence of RMS error as a function of quadrature spacing m for the MMS problem with N = 101, h = 0.01 finite element nodes (see Fig. 2). Optimal quadrature interval m* is given by (11) yields RMS = 1.25 × 10−4, a value close to the asymptotic limit when m < m*. Note the rapidly increasing RMS error as m* < m. 5 NUMERICAL RESULTS Our fractional Helmholtz system is numerically demonstrated in the context of MT. This is a geophysical surveying method that measures naturally occurring, time-varying magnetic and electric fields. Resistivity estimates of the subsurface can be inferred from the very near surface to hundreds of kilometres that are applied to subsurface characterization for myriad, broad-reaching geoscience applications such as hydrocarbon extraction, geothermal energy harvesting and carbon sequestration, as well as studies into Earth’s deep tectonic history. The MT signal is caused by the interaction of the solar wind with the earth’s magnetic field and world-wide thunderstorms which, at frequencies greater than 1 Hz, augment those signals of magnetospheric origin. Fig. 5 provides a conceptual diagram of MT in which an idealized ionospheric ‘sheet current’ is the electromagnetic source (see Section 2) for inducing secondary electromagnet eddy currents in the subsurface. Because, in practice, the magnitude of this source current is unknown, the fundamental quantity for MT analysis is the impedance tensor mapping electric correlated electric and magnetic fields. Computed in the frequency domain, the impedance tensor is an estimate of the Earth ‘filter’ mapping magnetic to electric fields—in other words, it is an expression of Earth’s conductivity distribution. Common in preliminary MT analysis is the assumption of locally 1-D (depth dependent) electrical structure and excitation by a vertically incident plane wave, such as described in Section 2, an assumption founded in the earliest works on the subject, but later criticized for its physical consistency (see ch. 2 in Chave & Jones 2012). Nonetheless, we adopt these modest assumptions in our investigation of MT data—in particular, data collected by the decadal, transcontinental USArray/EarthScope project—and find examples where MT data are consistent with predicted impedances for a fractionally diffusing electromagnetic Earth. Figure 5. Open in new tabDownload slide Overview of magnetotelluric experiment and data reduction. Left-hand panel: collocated time-series of horizontal electric field, measured by pairs of grounded electrodes, and magnetic field, measured by induction coils or fluxgates, measure Earth’s inductive response to ionspheric source currents. Right-hand panel: time-series are windowed, filtered and transformed into the frequency domain, from which the impedance tensor is estimated, containing information on the distribution of electricical conductivity variations in Earth’s subsurface (Chave & Jones 2012). Because high-frequency fields decay more rapidly with depth than low frequency fields, frequency can loosely be interpreted as a proxy for depth, and hence an impedance spectrum is a coarse measure of the local, depth variations in electrical conductivity. Because of the novelty in applying fractional derivative concepts to electromagnetic geophysics, the first question that draws our attention is simply: How does a fractionally diffusing field, as described by (3) and (4), compare to a field derived from the classical Helmholtz equation? To address this question we solve (4) on the dimensionless unit interval ζ ∈ [0, 1] with unit amplitude Dirichlet conditions |$u(0)=\sqrt{2\text{i}}$| and u(1) = 0 on the horizontal electric field and choose the dimensioned scaling factor z* = 1000 m to represent the physical domain z ∈ [0, 1000] m. Choice of homogeneous Dirichlet condition at ζ = 1 is commonly known as ‘perfectly conducting’ boundary condition, representing the presence of an infinitely conductive region for ζ > 1, but is used here strictly out of computational convenience. Secondary fields from this interface back to Earth’s surface ζ = 0 will be negligible as long as the frequency ω in eq. (4) is sufficiently high that the electric field at depth is essentially zero. The unit interval is discretized with 501 evenly distributed nodes, on which the electric field is drawn from the finite dimensional vector space of linear nodal finite elements. Hence, node spacing is h = 0.002, which, when s = 0.7 for example, leads to N− = 318 and N+ = 137 according to (11) and a linear system (27) with 501(3 + N− + N+) = 229 458 equations. Comparable to the error tolerances on the BiCG-STAB solver specified previously for the MMS problem, the iterative sequence is terminated once the normalized residual is reduced by 10−12 over its starting value. The horizontal electric fields in Fig. 6 show depth-dependent behaviour that is clearly also s-dependent: increased curvature in the near-surface and decreased curvature at depth in comparison with their classical s = 1.0 counterpart. This suggests that the effect of the fractional Laplacian in (3) and (4) over a uniform σα, ζ earth model is, at first blush, in some ways similar to that of a classical Laplacian over a layered Earth which is conductive in the near surface and resistive at depth. However, closer inspection of the fractional response (see, for example the s = 0.60 curves) reveals that the damped oscillations, characteristic of classical Helmholtz, are simply not present as s decreases from unity, and instead then are replaced with a steady non-oscillatory decay with depth. Figure 6. Open in new tabDownload slide For a range of fractional exponent values s = 0.6, 0.7, 0.8, 0.9, 0.995, decay of unit-amplitude real (top panel) imaginary (bottom panel) components of horizontal electric field Ex = u as a function of depth z into a uniform σα, ζ = 0.01 S m–1 medium at frequency f = 1 kHz, corresponding to dimensionless wavenumber κ ≈ 8.89. In red is the analytic solution for the corresponding classical s = 1 Helmholtz. See text for additional details on boundary conditions and scaling to the physical domain from the dimensionless unit interval. There is a dramatic manifestation of this fractional Helmholtz response in observable MT data through calculation of the impedance spectrum (Fig. 7). Amplitude of the impedance spectrum, reported here as the familiar apparent resistivity $$\begin{eqnarray*} \rho _a = {1\over \omega \mu _0} \Big \vert {E_x \over H_y} \Big \vert ^2_{z=0} = \omega \mu _0 (z^*)^2 \Big \vert {u\over \partial _\zeta u}\Big \vert ^2_{\zeta = 0} \end{eqnarray*}$$(30) and complex phase angle θ of the ratio −u/∂ζu, show a clear s-dependence at frequencies above 1 Hz (Fig. 7). Decay of the apparent resistivity as frequency approaches zero can be understood as a consequence of the perfect electric conductor boundary condition at z = z*, where at these low frequencies the reciprocal wavenumber |${1\over \kappa }\gt \gt z^*$| and hence the apparent resistivity approaches that of the perfect conductor, zero, in the region z > z*. Furthermore, in the limit of zero frequency, the fractional Helmholtz equation asymptotes to the fractional Laplacian equation [analogous to (8)] with inhomogeneous Dirichlet boundary conditions, whose solution has already been established Antil et al. (2018b) as equivalent to the classical Laplacian equation, leaving the ratio −u/∂ζu = 1, or equivalently θ = 0. Figure 7. Open in new tabDownload slide Magnetotelluric sounding curves for a uniform σα, ζ earth model over the depth domain z ∈ [0, z*] for a range of fractional exponent values s = 0.5, …, 1.0 with a perfect conductor boundary condition at z = z*. Apparent resistivity (top panel); complex phase (bottom panel). See text for description. In red are the classical s = 1 Helmholtz solutions, computed analytically. The decrease in apparent resistivity at high frequencies when s ≠ 1 can further be understood by examination of the electric field gradient at z = 0 (Fig. 8). Although there is a slight decrease in the vertical gradient of the imaginary component of electric field when s ≠ 1, the magnitude of the real component increases dramatically in comparison to the s = 1 case. This overall rise in vertical gradient at the air–earth interface for a fractional earth model decreases the value of the quotient in (30), thereby leading to a decreased estimate of the apparent resistivity at large frequencies. Figure 8. Open in new tabDownload slide Real (top panel) and imaginary (bottom panel) components of horizontal electric field as a function of depth in a uniform 0.01 S m–1 Earth underlain by a perfect conductor for frequencies f = 316, 1000 and 3162 Hz, corresponding to the high-frequency region of the magnetotelluric apparent resistivity spectrum (Fig. 6) with approximate s dependent power law behaviour.Curves for classical s = 1 (heavy lines) and fractional s = 0.7 response (light lines) are shown.The decrease in apparent resistivity is evidently due to the strongly increased vertical gradient of Real component of electric field at the air–Earth interface z = 0 for fractional Helmholtz. Recall that from eq. (4.1) that the vertical gradient of electric fields resides in the denominator of the of the apparent resistivity estimator. 5.1 Validation through USArray data We have made progress towards validating our hypothesis of ‘fractional Helmholtz leading to new geophysical interpretation’ though geophysical insight of numerical experiments. These numerical results suggest that a uniform subsurface subject to superdiffusing electromagnetic plane waves exhibits a high apparent conductivity in the near surface while appearing more resistive at depth (Fig. 7). This is consistent with actual field data observed in the mid-continent of the USArray footprint. In these data, apparent resistivity and phase angle spectra from USArray MT station for KSP34 located NW of Kansas City, KS, USA are consistent with the non-local behaviour underlying our numerical experiments (Fig. 9). An in-depth study of the geology in the Kansas City region would surely drive additional geological (and source physics) realism in subsequent analyses but, clearly, is beyond the scope of this paper. Nonetheless, these field data correlations provide compelling motivation to support additional algorithmic development for fractional electromagnetics. Such development would naturally complement the concurrent investigation of fractional calculus processes within the MT source fields themselves, as evidenced recently by careful analysis of non-Gaussian, heavy-tailed MT data residuals (Chave 2014). Partitioning the observable data between signals due to fractional Earth and magnetospheric effects may be possible by exploiting the differences in lateral complexity of these two regions, respectively, where the former is likely (at least in continental regions) far rougher that the latter. Figure 9. Open in new tabDownload slide Apparent resistivity and phase angle data from USArray MT station for KSP34 located NW of Kansas City, KS, USA from the US Array: (a) Apparent resisistivity spectrum based on Zxy (blue) and Zyx (red) elements of the 2 × 2 impedance tensor Z. The similarity in the curves, especially at high frequencies, is indicative a locally 1-D conductivity profile beneath the observation point. (b) Location map for USArray MT station KSP34. (c) Complex phase angle of the ratios Ex/Hy (Zxy, blue) and Ey/Hx (Zyx, red), again generally similar and indicating a locally 1-D conductivity profile beneath the station. (d) Depth profile of electrical resistivity beneath station KSP34 estimated by 3-D inversion of all sites in subpanel (b) [Fig. 3, subpanel (a), Yang et al., EPSL 2015]. 5.2 Strategies for a spatially variable fractional exponent An initial assumption in problem statement (5) is the spatial invariance of the fractional exponent s over the spatial domain Ω. However, if s is interpreted to represent via non-locality some degree of long-range correlation of underlying material properties (e.g. electrical conductivity), then it is relevant to consider how spatial variability in this correlation is accommodated in the architecture of the fractional calculus paradigm. In addition, variability in s enables us to truly capture the non-smooth effects such as fractures by prescribing variable degree of smoothness across the scales. A detailed analysis for variable s, where the authors have created a time-cylinder based approach, has been recently carried out in Antil & Rautenberg (2019). For a precise definition of the fractional Laplacian with variable s we refer to Antil et al. (2018a). In the case of a piecewise constant s, a conceptually simple strategy is to decompose the domain Ω into subdomains on which s is constant and impose our Kato method over each of the subdomains. Note that the solution for w in (8) is independent of s and may be obtained without any need for domain decomposition. Although differences in s among domains means that the number of functions vℓ = 0, …, L also varies among domains, the boundary condition v = 0 on each of the subdomains ensures continuity of v, and therefore continuity of u = v + w throughout Ω. Observe that computation of {vℓ} in one subdomain is independent of its calculation in another, and hence, {vℓ} over each of the subdomains can computed in parallel with no message passing or inter-domain communication required once w is solved for and shared globally throughout Ω. That said, several issues need to be resolved before this idea can be defensibly implemented. First, the suggestion of zero (subdomain) boundary conditions on {vℓ} needs to be physically justified. If found to be unsound, the embarrassingly parallel structure just described will instead require interdomain communication and potentially interpolation. Secondly, because the fractional Laplacian is inherently non-local, its support extends over the global domain Ω. Ensuring global extent of non-locality in the context of subdomains requires further analysis. Lastly, function and flux continuity for a given vℓ within a subdomain is guaranteed; the conditions for such guarantees, in a general sense, at subdomain boundaries have yet to be determined. Because of these complexities, further analysis of this domain decomposition concept is deferred to future publication. 5.3 Fractional time derivatives Prior work in electromagnetic geophysics in contemplation of Ohm’s constitutive law being represented in terms of fractional calculus have focused on fractional time derivatives, rather than the fractional space derivatives described here (Weiss & Everett 2007; Everett 2009; Ge et al. 2015). Such analyses are comparatively simple in that the fractional space derivatives |$D_z^\alpha$| of (1.2) are replaced by time derivatives |$D_t^\beta$|⁠, thus modifying the complex wavenumber as k2 = −(iω)1 − βμ0σ. Solutions to (5) in layered media when s = 1 (equivalently, α = 0 since |$s=1-{\alpha \over 2}$|⁠) follow the usual method of posing characteristic solutions exp (± kz) in each of the layers, coefficients for which are determined through enforcement of boundary condition (5) along with continuity of u and ∂zu at layer boundaries. Solving this time-fractional Helmholtz equation on the domain Ω: z ∈ [0, z* = 1000] m with u(0) = 1 + i, u(z*) = 0 and σ = 0.01 S m–1(rad s–1)−β yields a characteristic MT response (Fig. 10) distinct from that obtained in the case of space-fractional derivatives s ≠ 1 (Fig. 7). As noted in Ge et al. (2015), imposing the time-fractional derivative in this way is equivalent to recasting real-valued electrical conductivity σ as a frequency-dependent, complex-valued conductivity σ(iω)β. The quasi-linear power-law behaviour in apparent resistivity and phase angle (Fig. 10) seen at high frequencies (f > 100 Hz) is objectively distinct from that computed for the space-fractional Helmholtz system (Fig. 7) and offers an unambiguous diagnostic for discriminating between the two. These differences have their origin in the how anomalous power-law diffusion is captured by each. In the case of fractional time derivatives of order 1 − β, as considered in this latest example, the system is considered subdiffusive and consistent with an anomalously high likelihood of long wait times between successive jumps of charge carriers in a continuous time random walk, analogous to that suggested for fluid transport in a porous medium (Metzler & Klafter 2000). Instead, the space-fractional derivatives which occupy the primary focus of this study capture long-range interactions (spatial non-locality) of charge carriers as a superdiffusive system, perhaps through inductive coupling (a phenomena absent in the physics of fluid flow in porous media). This contrast—super- versus subdiffusion—is the essence of the causative physics behind the different MT responses predicted by (Figs 7 and 10). Figure 10. Open in new tabDownload slide Apparent resistivity spectra for classical Helmholtz equation s = 1 with iκ2 = (iω)1 − βμ0σ, where the β terms arise in the 1-D magnetotelluric case from Ohm’s law with a fractional, non-local time dependence attributable to subdiffusion of electric charge following a continuous time random walk with a heavy tailed distribution of waiting times. Compare to Fig. 7 for the space-fractional case describing superdiffusion, where the heavy-tailed distribution of step length (a.k.a. Lévy flights) captures long-range interactions between charges. 6 CONCLUSIONS We have presented a novel, practical solution to the fractional Helmholtz equation based on the Kato formulation of the fractional Laplacian operator, a lifting (splitting) strategy to handle non-zero Dirichlet boundary conditions, and finite element discretization of the spatial domain. This specific finite element discretization derives from our statement of the variational problem, from which alternative discretizations (orthogonal polynomials, wavelets, spectral functions, etc.) offer an interesting direction for future research. Whereas the analogous Kato/lifting strategy for solving the fractional Poisson equation leads to decoupled system of integer Laplace solves which can be done in parallel with no intersolve communication, solution of the fractional Helmholtz admits no such decoupling. This leads to a large, block-dense system of linear equations upon discretization which significantly increases the resource requirements for obtaining a numerical solution. In response, we augment the variational problem by introducing an additional unknown which collapses the L − 1 block coupling matrices in a given block-row into a single block matrix, at the expense of only one additional (dense) block-row in the linear system. For typical problems where with L > >100, the added computational burden of this compatibility equation is inconsequential, yet the reduction in matrix storage is significant, going from L2 to simply 2L. Thus, a key feature of this augmented variational problem is the extreme block-sparsity of the resulting linear system of equations, a feature which is independent of the choice of discretization and important for efficient solution of large-scale systems. Validation of the algorithm for linear, nodal finite elements shows h2 reduction in RMS error for an MMS test problem—demonstrating that our formulation of the fractional Helmholtz problem does not corrupt the convergence behaviour expected from solution of integer-order Helmholtz. We apply this formulation for fractional Helmholtz to the growing body of observational evidence of anomalous diffusion in nature—here, asking the question, ‘Does the Earth, with its incalculable geological complexity, respond to electromagnetic stimulation in a way that is consistent with fractional diffusion and the non-locality that is central to the differential operators of the governing physics?’ Whereas temporal non-locality of Maxwell’s equations has previously been observed as subdiffusive propagation, the fractional Helmholtz equation studied here describes superdiffusion by attributing fractional derivatives directly to the spatial distribution of material properties in Ohm’s constitutive law. Earth electromagnetic response is computed in the context MT analysis—a classic geophysical exploration technique dating back to to middle 20th century—and comparison with the EarthScope USArray database. We find qualitative agreement between the predicted fractional Helmholtz response functions and those observed at a middle North American measurement site. This congruence in electromagnetic response thereby offers an alternative interpretation of the MT data at the site, one where the classical interpretation of a layered Earth geology with deep resistive rocks overlain by a conductive overburden is contrasted with new interpretation suggesting complex, geological texture consistent with the site’s proximity significant deep crustal tectonic structure. Outstanding issues for future research therefore lie in two fundamental areas: arriving at a clearer mapping between the value of a fractional exponent s and the material heterogeneity it’s intended to represent; and, extension of the computational tools to higher dimension with parallel implementation, including spatially variable and/or anisotropic s values. The former may be informed, as we have done here, by reinterpretation of existing observational data through fractional calculus concepts, but augmented by detailed material analysis. The latter naturally feeds into ongoing efforts in PDE-constrained optimization for material property estimation, now augmented with the desire to recover s, too, as a measure of material complexity or subgrid structure. Algorithmic advances in multilevel domain decomposition (decomposition over physical domain in addition to decomposition over the functional blocks of global system matrix) will also be required for full exploration of fractional Helmholtz concepts on large, 3-D domains. ACKNOWLEDGEMENTS This interdisciplinary, collaborative work is supported by Sandia National Laboratories LDRD program, NSF grants DMS-1521590 and DMS-1818772, and Air Force Office of Scientific Research under Award NO: FA9550-19-1-0036. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. REFERENCES Antil H. , Leykekhman D., 2018 . A brief introduction to PDE constrained optimization , in Frontiers in PDE-Constrained Optimization , Vol. 163 of IMA Volumes in Mathematics and Computation , p. 434 , eds Antil H., Kouri D.P., Lacasse M.D., Ridzal D., Springer-Verlag . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Antil H. , Pfefferer J., 2017 . A short MATLAB implementation of fractional Poisson equation with nonzero boundary conditions , Technical Report , Available at: http://math.gmu.edu/ hantil/Tech_Report/HAntil_JPfefferer_2017a.pdf . OpenURL Placeholder Text WorldCat Antil H. , Rautenberg C., 2019 . Sobolev spaces with non-Muckenhoupt weights, fractional elliptic operators, and applications , SIAM J. Math. Anal. , 51 ( 3 ), 2479 – 2503 . Google Scholar Crossref Search ADS WorldCat Antil H. , Verma D., Warma M., 2019 . Optimal control of fractional semilinear PDEs , ESAIM: Control, Optimisation and Calculus of Variations (ESAIM: COCV) , in press, doi:10.1051/cocv/2019003 . Google Scholar OpenURL Placeholder Text WorldCat Antil H. , Pfefferer J., Warma M., 2017 . A note on semilinear fractional elliptic equation: analysis and discretization , Math. Modelli. Numer. Anal. , 51 ( 6 ), 2049 – 2067 . Google Scholar Crossref Search ADS WorldCat Antil H. , Ceretani A., Rautenberg C.N., 2018a . The fractional Laplacian with spatially variable order , in preparation . OpenURL Placeholder Text WorldCat Antil H. , Pfefferer J., Rogovs S., 2018b . Fractional operators with inhomogeneous boundary conditions: Analysis, control, and discretization , Commun. Math. Sci. (CMS) , 16 ( 5 ), 1395 – 1426 . Google Scholar Crossref Search ADS WorldCat Baeumer M. , Meerschaert M.M., Nane E., 2009 . Space–time duality for fractional diffusion , J. Appl. Prob. , 46 , 1100 – 1115 . Google Scholar Crossref Search ADS WorldCat Balakrishnan V. , 1960 . Fractional powers of closed operators and the semi-groups generated by them , Pacific J. Math. , 10 ( 2 ), 419 – 437 . Google Scholar Crossref Search ADS WorldCat Benson D.A. , Wheatcraft S.W., Meerschaert M.M., 2000 . The fractional-order governing equations of Lévy motion , Wat. Res. Res. , 36 ( 6 ), 1413 – 1423 . Google Scholar Crossref Search ADS WorldCat Berggren M. , 2004 . Approximations of very weak solutions to boundary–value problems , SIAM J. Numer. Anal. , 42 , 860 – 877 . Google Scholar Crossref Search ADS WorldCat Bonito A. , Pasciak J., 2017 . Numerical approximation of the fractional powers of regularly accretive operators , IMA J. Numer. Anal. , 37 , 1245 – 1273 . Google Scholar Crossref Search ADS WorldCat Bonito A. , Pasciak J.E., 2015 . Numerical approximation of fractional powers of elliptic operators , Math. Comp. , 84 ( 295 ), 2083 – 2110 . Google Scholar Crossref Search ADS WorldCat Caffarelli L.A. , Stinga P.R., 2016 . Fractional elliptic equations, Caccioppoli estimates and regularity , Annales de l’Institut Henri Poincaré. Analyse Non Linéaire , 33 ( 3 ), 767 – 807 . Google Scholar Crossref Search ADS WorldCat Caputo M. , 2000 . Models of fluids in porous media with memory , Wat. Res. Res. , 36 ( 3 ), 693 – 705 . Google Scholar Crossref Search ADS WorldCat Chave A.D. , 2014 . Matnetotelluric data, stable distributions, and impropriety: and existential combination , Geophys. J. Int. , 198 , 622 – 636 . Google Scholar Crossref Search ADS WorldCat Chave A.D. , Jones A.G., 2012 . The Magnetotelluric Method: Theory and Practice , Cambridge Univ. Press . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Deseri L. , Zingales M., 2015 . A mechanical picutre of frational-order Darcy equation , Commun. Nonlin. Sci. Numer. Simulat. , 20 , 940 – 949 . Google Scholar Crossref Search ADS WorldCat Einstein A. , 1905 . Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teichlen , Annalen der Physik , 322 , 549 – 560 . Google Scholar Crossref Search ADS WorldCat Elser J. , Podolskiy V.A., Salakhutdinov I., Avrutsky I., 2007 . Nonlocal effects in effective–medium response of nanolayered metamaterials , Appl. Phys. Lett. , 90 , 191109 , doi: 10.1063/1.2737935 . Google Scholar Crossref Search ADS WorldCat Evans R.L. , 2012 . Conductivity of earth materials , in Magnetotelluric Method: Theory and Practice , chap. 3A , pp. 50 – 95 ., eds Chave A.D., Jones A.G., Cambridge Univ. Press . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Everett M.E. , 2009 . Transient electromagnetic response of a loop source over a rough geologic medium , J. geophys. Int. , 177 , 421 – 429 . Google Scholar Crossref Search ADS WorldCat Everett M.E. , Weiss C.J., 2002 . Geological noise in near-surface electromagnetic induction data , Geophys. Res. Lett. , 29 ( 1 ), 1010 , doi: 10.1029/2001GL014049 . Google Scholar Crossref Search ADS WorldCat Ge J. , Everett M.E., Weiss C.J., 2015 . Fractional diffusion analysis of the electromagnetic field in fractured media – Part 2: 3D approach , Geophysics , 80 , E175 – E185 . Google Scholar Crossref Search ADS WorldCat Karato S.-I. , Wang D., 2013 . Electrical conductivity in rocks and minerals , in Physics and Chemistry of the Deep Earth , chap. 5 , pp. 145 – 184 ., ed. Karato S.-I., John Wiley & Sons, Ltd ., doi: 10.1002/9781118529492.ch5 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Kato T. , 1960 . Note on fractional powers of linear operators , Proc. Jpn. Acad. , 36 ( 3) , 94–96 . Google Scholar OpenURL Placeholder Text WorldCat Kumar G.P. , Chaudhary I., Nagar M., Chouhan A.K., Prizomwala S.P., Mahesh P., Chopra S., 2018 . Transient electromagnetic investigations in a tectonic domain of the Kachchh intraplate region, western India: a morphotectonic study of the Kachchh Mainland Fault , Tectonics , 37 , 4239 – 4260 . Google Scholar Crossref Search ADS WorldCat Lions J.L. , Magenes E., 1972 . Non-Homogeneous Boundary Value Problems and Applications , Springer–Verlag . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Meerschaert M.M. , Sikorskii A., 2012 . Stochastic Models for Fractional Calculus , Walter de Gruyter GmbH & Co. KG . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Metzler R. , Klafter J., 2000 . The random walk’s guide to anomalous diffusion: a fractional dynamics approach , Phys. Rep. , 339 , 1 – 77 . Google Scholar Crossref Search ADS WorldCat Palacky G.J. , 1987 . Resistivity properties of geologic targets , in Electromagnetic Methods in Applied Geophysics , Vol. 3 of Investigations in Geophysics , p. 513 , ed. Nabighian M.N., Society of Exploration Geophysicists . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Roache P.J. , 1998 . Verification of codes and calculations , Am. Inst. Aeron. Astron. J. , 36 , 696 – 702 . Google Scholar Crossref Search ADS WorldCat Salar K. , Knupp P., 2000 . Code verification by the method of manufactured solutions , Tech. rep., Sandia National Laboratories, SAND2000-1444 . OpenURL Placeholder Text WorldCat Song R. , Vondraček Z., 2008 . On the relationship between subordinate killed and killed subordinate processes , Electr. Commun. Prob. , 13 , 325 – 336 . Google Scholar Crossref Search ADS WorldCat Vallianatos F. , 2017 . Transient electromagnetic method in the Keritis basin (Crete, Greece): evidence of hierarchy in a complex geological structure in view of Tsallis distribution , Ann. Geophys. , 60 , GM675 . Google Scholar Crossref Search ADS WorldCat Vallianatos F. , 2018 . A non extensive view of electrical resistivity spatial distribution estimated using inverted Transient Electromagnetic responses in a karstified formation (Keritis basin, Crete, Greece) , Phys. A , 505 , 171 – 178 . Google Scholar Crossref Search ADS WorldCat Vallianatos F. , Kouli M., Kalisperi D., 2018 . Evidence of hierarchy in the complex fracture system of Geropotamos (Crete, Greece), as extracted from transient electromagnetic responses , Pure appl. Geophys. , 175 , 1 – 4 . Google Scholar Crossref Search ADS WorldCat van der Vorst H.A. , 1992 . BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems , SIAM J. Scient. Stat. Comput. , 13 , 631 – 644 . Google Scholar Crossref Search ADS WorldCat von Kameke A. , Huhn F., Fenandez-Garcia G., Munuzuri A., Perez-Munuzuri V., 2010 . Propagation of a chemical wave front in a quad-two-dimensional superdiffusive flow , Phys. Rev. E , 81 ( 6 Pt 2 ), 066211. Epub 2010 Jun 18 . Google Scholar OpenURL Placeholder Text WorldCat Weiss C.J. , Everett M.E., 2007 . Anomalous diffusion of electromagnetic Eddy currents in geologic formations , J. geophys. Res. , 112 ( B8 ), doi: 10.1029/2006JB004475 . Google Scholar OpenURL Placeholder Text WorldCat Crossref Yosida K. , 1995 . Functional analysis , in Classics in Mathematics , Springer-Verlag , doi: 10.1007/978-3-642-61859-8 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crossref APPENDIX Faraday’s law of induction equates the circulation of electric field |$\bf E$| to the negative time rate of change of the magnetic induction |$\bf B$|⁠, $$\begin{eqnarray*} \nabla \times {\bf E} = -\partial _t{\bf B}. \end{eqnarray*}$$(A1) Assuming non-magnetic media where magnetic permeability μ is equal to that of free space μ = μ0 = 4π × 10−7 H m–1, the right-hand side of (A1) can be written as −μ0∂tH, where |$\bf H$| is the magnetic field. Similarly, Ampère’s Law equates the circulation of magnetic field to the total electric current density, Jtot, $$\begin{eqnarray*} \nabla \times {\bf H} = {\bf J}_\text{tot}, \end{eqnarray*}$$(A2) which in the absence of prescribed source currents—naturally occurring or engineered—is given by Ohm’s constitutive law, a convolution of the electrical conductivity σ and electric field over both space |$\boldsymbol r$| and time t: $$\begin{eqnarray*} {\bf J}_\text{tot} = \sigma * {\bf E}. \end{eqnarray*}$$(A3) The convolutional nature of the empirically derived constitutive laws of electromagnetics reflects not only the linearity between the vector fields such as |$\bf J$| and |$\bf E$|⁠, but also the fact that the effect of material properties may be non-local in space and/or time. The latter is representative of materials with ‘memory’, whose behaviour at a given time is dependent on the time-history of the material. This has been studied previously in the context of geophysical electromagnetics (Everett 2009) and motivated by observational evidence of anomalous power-law scaling (Weiss & Everett 2007) from an engineered loop antenna. The motivation for fractional derivatives, either in space or time, comes from the observation of anomalous diffusion in the moveout response of loop–loop CSEM data (Weiss & Everett 2007). A classically diffusing EM pulse results in an arrival time τ that is proportional to squared offset L2. Thus a log–log plot of these two quantities has slope equal to 1. Weiss & Everett (2007) present experimental evidence in which such a plot has slope not equal to 1, a result that is consistent with a governing PDE with fractional derivatives. Similar anomalous power laws have also been observed in the context of diffusive fluid transport in porous media (e.g. Benson et al. 2000). In attributing the anomalous diffusion to memory effects alone, Everett (2009) defines the spatially uniform but temporally non-local Ohm’s law as the convolution $$\begin{eqnarray*} \left( \sigma * {\bf E}\right) (t) \equiv {1\over \Gamma (\beta )}\int _0^t {\sigma _\beta {\bf E}(t^\prime )\over (t-t^\prime )^{1-\beta }}\, \text{d}t^\prime , \end{eqnarray*}$$(A4) where Γ( · ) is the Gamma function, σβ is electrical conductivity (with fractional units for dimensional consistency) and β describes the anomalous power-law relationship. The exact mapping between (piecewise?) constant σβ and the presumed, causative ‘rough geological medium’ σ remains an open research question. Nonetheless, substituting (A4) into (A3) and (A2) eliminates the magnetic field variable and introduces via (A3) and (A4) the term $$\begin{eqnarray*} {\partial \over \partial t}\left( \sigma * {\bf E}\right) (t) = {\partial \over \partial t}{1\over \Gamma (\beta )}\int _0^t {\sigma _\beta {\bf E}(t^\prime )\over (t-t^\prime )^{1-\beta }}\, \text{d}t^\prime . \end{eqnarray*}$$(A5) The term on the right is recognized as the Riemann–Liouville fractional derivative, in time, of order 1 − β, thusly motivated by a particular representation of Ohm’s convolutional constitutive law for consistency with observations of anomalous, fractional power law in geophysical electromagnetics. In such cases the governing diffusion equation possesses integer order derivatives in space coupled with fractional-order derivatives in time. Alternatively, an anomalous power law can also be captured through a governing diffusion equation whose space derivatives are of fractional order and the time derivatives are integer order. Indeed, mixed fractional space–time derivatives are now an accepted element of the mathematical architecture for quantifying anomalous hydrologic diffusion (e.g. Benson et al. 2000; Meerschaert & Sikorskii 2012). To motivate the proposed fractional space derivatives, we return to the fundamental definition of Ohm’s constitutive law and write it as local in time, but non-local in space. For the sake of notational simplicity, consider non-locality in the z direction alone and write analogously to (A4) $$\begin{eqnarray*} \left( \sigma * {\bf E}\right) ({z}) \propto {\partial \over \partial z}{1\over \Gamma (1-\alpha )}\int _{-\infty }^{+\infty } {\sigma _{\alpha ,z} (z^\prime ) {\bf E}(z^\prime )\text{sgn}(z-z^\prime )\over \vert z-z^\prime \vert ^\alpha }\, \text{d}z^\prime , \end{eqnarray*}$$(A6) where the signum function, sgn( · ), is equal to ±1, consistent with the sign of its real argument. A couple of differences between our space-convolution (A6) and the (Everett 2009) time-convolution (A4) merit some comment. First, the limits of integration in (A6) extend from ±∞ and in practice would be truncated at the bounds of some model domain Ω under consideration, whereas the limits on the time-convolution are truncated at time t to enforce causality. Because of this we may further decompose the convolution into an integration over two domains, one with z′ < z and the other with z′ > z: $$\begin{eqnarray*} \left( \sigma * {\bf E}\right) ({z}) \propto {\partial \over \partial z}{1\over \Gamma (1-\alpha )} \left[ \int _{-\infty }^z {\sigma _{\alpha ,z} (z^\prime ) {\bf E}(z^\prime )\over (z-z^\prime )^\alpha }\, \text{d}z^\prime - \int _z^{+\infty } {\sigma _{\alpha ,z} (z^\prime ) {\bf E}(z^\prime )\over (z^\prime -z)^\alpha }\, \text{d}z^\prime \right]. \end{eqnarray*}$$(A7) Thus, unlike the causality constraint for time derivatives, the spatial convolution is free to consider regions to both ‘forward’ and ‘backward’ from z. The second major difference between the proposed spatial convolution and the (Everett 2009) time-convolution is the leading derivative with respect to z. However, we point out that it’s trivial by invoking the Fourier transform to show that the derivative of a convolution remains a convolution, nonetheless. Thus, our expression still retains the familiar linearity properties of Ohm’s law while admitting the non-locality we desire. The first term in (A7) is recognizable (from above) as the fractional derivative of order α with respect to z. The second term, we find, is the fractional derivative of order α with respect to −z (Baeumer et al. 2009). Like derivatives of integer order, the Fourier transform of a α-order fractional derivative with respect to z of a function f(z) is simply |$(-i\nu )^\alpha \bar{f}(\nu )$|⁠, whereas the derivative with respect to −z yields Fourier transform |$(i\nu )^\alpha \bar{f}(\nu )$|⁠. Consequently, Fourier transform of the sum, as we have in (A7), has Fourier transform |$2\cos ({1\over 2}\alpha \pi )\vert -\nu \vert ^\alpha \bar{f}(z)$|⁠.We thereby define a fractional differential operator |$\mathcal {D}_z^\alpha$| with Fourier transform | − ν|α for the left/right derivatives in (A7) as $$\begin{eqnarray*} \mathcal {D}_z^\alpha = {1\over 2\cos ({1\over 2}\alpha \pi )} \left( {\partial ^\alpha \over \partial z^\alpha } + {\partial ^\alpha \over \partial (-z)^\alpha } \right). \end{eqnarray*}$$(A8) Note that as was done previously with σβ in (Everett 2009), the connection between the conductivity function σα,z the underlying presumably rough geology σ left undefined. This still remains an outstanding research question, but before too much effort is spent in pursuit of its answer, we argue the such answers are fruitless if the governing spatial space–fractional diffusion equation is too burdensome to solve or yields results inconsistent with observation. It is the latter which defines the focus of the present research. To do so, we consider the 1-D magnetotelluric problem where there is, again, a presumably rough and depth-dependent conductivity function σ(z), and the electric field is fully represented by its horizontal component, written here as u(z). Under these simplifying assumptions, the combination of Ampère’s and Faraday’s law, as described above through elimination of the magnetic field, yields a simple ordinary differential equation in the ω frequency domain: $$\begin{eqnarray*} -{\text{d}^2u(z)\over \text{d}z^2} = -\text{i}\omega \mu _0\, \mathcal {D}_z^\alpha \left[\sigma _{\alpha ,z} u(z)\right]. \end{eqnarray*}$$(A9) Recall and observe above that for dimensional consistency, the function σα, z necessarily has fractional physical units of S m(1 − α)–1. This awkwardness can be avoided by non-diminensionalizing the space coordinate z↦ζ = z/z* where z* is taken here to be the range of z in our model domain Ω, analogous to the τ parameter in Everett (2009) for scaling of the time variable t. In doing so, the fractional differential operator is written as |$\mathcal {D}_\zeta ^\alpha$| to represent its action with respect to the dimensionless variable ζ and the classically dimensioned conductivity σα, ζ = (z*)ασα, z is introduced such that $$\begin{eqnarray*} -{\text{d}^2u(\zeta )\over \text{d}\zeta ^2}\left({1\over z^*}\right)^2 = -\text{i}\omega \mu _0\, \mathcal {D}_z^\alpha \left[\sigma _{\alpha ,z} u(z)\right]. \end{eqnarray*}$$(A10) As pointed out in the discussion preceding (A8), the Fourier transform of our differential operator |$\mathcal {D}_z^\alpha$| is | − ν|α and the same holds true for |$\mathcal {D}_\zeta ^\alpha$| under ζ↦ν. As a consequence, the Fourier transform |$\mathscr {F}$| of (A10) can be written as $$\begin{eqnarray*} \vert -\nu \vert ^2 \mathscr {F}\left[ u \right] + \text{i}\omega \mu _0 \left(z^*\right)^2\vert -\nu \vert ^\alpha \mathscr {F}\left[ \sigma _{\alpha ,\zeta } u\right] = 0, \end{eqnarray*}$$(A11) which invites combining the ν terms to write $$\begin{eqnarray*} \vert -\nu \vert ^{2-\alpha } \mathscr {F}\left[ u \right] + \text{i}\omega \mu _0 \left( z^*\right)^2 \mathscr {F}\left[ \sigma _{\alpha ,\zeta } u\right] = 0, \end{eqnarray*}$$(A12) an expression in the ζ domain that reads $$\begin{eqnarray*} \left( -{\text{d}^2 u(\zeta )\over \text{d}z^2 }\right)^s + \text{i}\omega \mu _0 \left( z^*\right)^2 \sigma _{\alpha ,\zeta } u(\zeta ) = 0 \end{eqnarray*}$$(A13) with |$s=1-{1\over 2}\alpha$|⁠.The 3-D analogue of (A13) follows immediately in the form of a generalized fractional-order Helmholtz equation ( − Δ)su − k2u = 0, whose fractional Laplacian operator and methods for working with it constitutes the bulk of the algorithmic development in the main body of this work. Published by Oxford University Press on behalf of The Royal Astronomical Society 2019. This work is written by (a) US Government employee(s) and is in the public domain in the US. Published by Oxford University Press on behalf of The Royal Astronomical Society 2019.
TI - Fractional operators applied to geophysical electromagnetics
JF - Geophysical Journal International
DO - 10.1093/gji/ggz516
DA - 2020-02-01
UR - https://www.deepdyve.com/lp/oxford-university-press/fractional-operators-applied-to-geophysical-electromagnetics-UGtNXs6bXh
SP - 1242
EP - 1259
VL - 220
IS - 2
DP - DeepDyve
ER -