On the approximation of spatial structures of global tidal magnetic field models

The extraction of the magnetic signal induced by the oceanic M2 tide is typically based solely on the temporal periodicity of the signal. Here, we propose a system of tailored trial functions that additionally takes the spatial constraint into account that the sources of the signal are localized within the oceans. This construction requires knowledge of the underlying conductivity model but not of the inducing tidal current velocity. Approximations of existing tidal magnetic field models with these trial functions and comparisons with approximations based on other localized and nonlocalized trial functions are illustrated.

1 Introduction

Conductive seawater moving through the ambient Earth's main magnetic field B_main induces a secondary magnetic field B_oc. Due to their periodic nature, magnetic signals generated by ocean tides are particularly easy to detect and have been studied in observatory data as early as, e.g., Malin (1970). However, the extraction of global models for magnetic fields induced by the dominating M2 tide from satellite data has only recently become possible (e.g., Sabaka et al., 2015, 2016; Tyler et al., 2003). Although the extraction procedures used are solely based on the temporal periodicity of the tidal signal (and not on further information on the spatial localization of the sources), they seem, by visual inspection, to coincide very well with results obtained by forward models such as in Kuvshinov and Olsen (2005). A more extensive comparison of forward models of electromagnetic ocean tidal signals based on different ocean tide models has recently been published in Saynisch et al. (2018). In that work, it was shown that the residuals between the different models can exceed the nominal noise level of the Swarm satellite mission. The ability to extract M2 tidal magnetic field signals in satellite data more precisely can therefore help in constraining ocean tide models. In Grayver et al. (2016) and Schnepf et al. (2015) it has been shown that an M2 tidal magnetic field model can also be used to constrain 1-D models of the Earth's mantle conductivity, and forward studies in Irrgang et al. (2016) have shown that lateral variations in the conductivity of the ocean water itself should have a detectable influence on the measured magnetic field (although the latter study was performed for general ocean circulation and not for tidal current systems).

In this short paper, we want to illustrate the effects that different (spatially localized) sets of trial functions can have on the approximation of the magnetic field induced by the M2 tide in the first place.

In particular, we describe a possible setup for the inclusion of spatial localization constraints (in addition to the constraint of temporal periodicity) for the approximation of ocean-tide-induced magnetic fields. Clearly, the velocity field u that is responsible for the generation of the corresponding secondary magnetic field B_oc vanishes over the continents. The precise connection between u and B_oc is given by the time-harmonic Maxwell equations:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{\nabla }\times {\mathit{B}}_{\mathrm{oc}}={\mathit{\mu }}_{\mathrm{0}}\mathit{\sigma }(\mathit{E}+\mathit{u}\times {\mathit{B}}_{\mathrm{main}}),\\ \text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\times \mathit{E}=i\mathit{\omega }{\mathit{B}}_{\mathrm{oc}},{\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot {\mathit{B}}_{\mathrm{oc}}=\mathrm{0},\end{array}$$

where we assume to have knowledge of B_main, the underlying conductivity σ, and the frequency ω. Furthermore, E denotes the electric field and μ₀ is the vacuum permeability. Instead of using a fixed velocity field u, we substitute it by a set of functions $\mathit{\{}{\mathit{u}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$ (e.g., vectorial Slepian functions as in Plattner and Simons (2014, 2015) that are localized over the oceans) to obtain a set of corresponding trial functions $\mathit{\{}{\mathit{B}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$ that each solve Eq. (1). The latter is suitable for the approximation of B_oc and reflects the spatial localization of the sources of the induced magnetic signal in the oceans. Thus, a magnetic field model that is based on an expansion of the signal in terms of the function system $\mathit{\{}{\mathit{B}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$ automatically reflects the spatial origin of the signal as well as its temporal periodicity (described by the frequency ω). Additionally, due to the linear connection between u and B_oc, an approximation of B_oc directly yields an approximation of the underlying tidal current velocity u in terms of the functions $\mathit{\{}{\mathit{u}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$. However, a model of the underlying conductivity σ has to be assumed for the construction of the B_ℓ. Throughout this paper, we fix the underlying conductivity, meaning that we do not test the influence of a variation of the conductivity model on the approximation of B_oc. The goal of the paper is rather the illustration of the effect of the general constraint that the (unknown) underlying u is restricted to the oceans. In a forthcoming study, the simultaneous reconstruction of u and approximation of B_oc, and a comparison with existing models, shall be investigated more thoroughly. Since the connection between σ and B_oc is nonlinear, a simultaneous determination of σ and B_oc (assuming a fixed velocity field model for u) is not as straightforward. A detailed description of the trial functions is provided in Sect. 2.

In Sect. 3, we illustrate our approach with input data derived from the (satellite- and observatory-data-based) CM5 model of Sabaka et al. (2015) and from data derived from a forward model based on the X3DG solver from Kuvshinov (2008). We approximate these input data sets separately in terms of time-periodic vector spherical harmonics, a system of spatially localized trial functions that contains no particular information on the underlying sources (in our case, Abel–Poisson kernels), and the new set of trial functions indicated in the previous paragraph, respectively. We also include an example with artificial continental noise. The residuals with respect to the input data show that the use of the function system $\mathit{\{}{\mathit{B}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$ can filter out undesired contributions to the M2 tidal magnetic field over the continents, without neglecting data over the continents. These residuals can reach up to 15% of the maximal signal strength and have a magnitude that should be detectable at satellite altitude.

Figure 1 The kernel $K(r\cdot ,a{\mathit{\eta }}_{\mathrm{1}})$ for $\frac{a}{r}=\mathrm{0.91}$ (a) and $\frac{a}{r}=\mathrm{0.67}$ (b). The fixed nodal point η₁∈𝕊 is marked by a white cross.

2 Method and function systems

Given a so-called dictionary 𝒟 of trial functions, we use the Regularized (Orthogonal) Functional Matching Pursuit (cf. Fischer and Michel, 2013; Michel and Telschow, 2014, 2016 for details) for the approximation of B_oc. Briefly speaking, this is a greedy-type algorithm that yields an approximation

$${\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{B}}}_N=\sum _{i=\mathrm{1}}^N{\mathit{\alpha }}_i{\mathit{d}}_i$$

of B_oc by iteratively choosing coefficients α_i∈ℝ and dictionary elements d_i∈𝒟 via

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathrm{argmin}}_{\mathit{\alpha },\mathit{d}}\left(\Vert {\mathit{R}}_{i-\mathrm{1}}-\mathit{\alpha }\mathcal{F}\mathit{d}{\Vert }_{{\mathbb{R}}^M}^{\mathrm{2}}+\mathit{\lambda }\Vert {\stackrel{\mathrm{‾}}{\mathit{B}}}_{i-\mathrm{1}}+\mathit{\alpha }\mathit{d}{\Vert }_{\mathcal{H}}^{\mathrm{2}}\right).\end{array}$$

${\mathit{R}}_{i-\mathrm{1}}=\mathit{b}-\mathcal{F}{\stackrel{\mathrm{‾}}{\mathit{B}}}_{i-\mathrm{1}}$ denotes the residual between the data b∈ℝ^M and the approximation after i−1 iterations. In this particular setup, ℱ represents the linear operator that evaluates a function at the M locations where data are provided, and ℋ is a suitable Hilbert space for the regularization of the problem¹. The parameter λ controls the trade-off between the data misfit $\Vert {\mathit{R}}_i{\Vert }_{{\mathbb{R}}^M}^{\mathrm{2}}$ and the regularizing term $\Vert {\stackrel{\mathrm{‾}}{\mathit{B}}}_i{\Vert }_{\mathcal{H}}^{\mathrm{2}}$, which imposes a certain property to ${\stackrel{\mathrm{‾}}{\mathit{B}}}_i$ such as smoothness (as in our case). The Regularized (Orthogonal) Functional Matching Pursuit, in general, has the advantage that it can easily deal with different dictionaries 𝒟 (or combinations of such) from which the approximant ${\stackrel{\mathrm{‾}}{\mathit{B}}}_N$ is built. However, any other approximation method could be used as well with the proposed function systems. In this paper, we use the term “dictionary” simply to describe a set of arbitrary functions that we consider suitable for our purposes. These functions do not necessarily have to satisfy particular mathematical properties such as orthogonality or completeness. Therefore, we call such functions “trial functions” rather than, e.g., “basis functions”.

In the following, we briefly introduce some function systems that can be used for the constitution of 𝒟. In particular, Sect. 2.4 describes in more detail the aforementioned trial functions $\mathit{\{}{\mathit{B}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$ that contain temporal and spatial constraints tailored for ocean-tide-induced magnetic fields.

Figure 2 Absolute value of the vectorial Slepian function ${\mathit{g}}_{\mathrm{50}}^{\left(\mathrm{3}\right)}$ with 50th best localization over the oceans (a) and ${\mathit{g}}_{\mathrm{1630}}^{\left(\mathrm{3}\right)}$ with the 50th worst localization over the oceans (b), for bandlimit N=40.

2.1 Vector spherical harmonics

We briefly recapitulate the notion of classical vector spherical harmonics in a form that we need at a few occasions later on. By ${\mathbb{S}}_r=\mathit{\{}x\in {\mathbb{R}}^{\mathrm{3}}:|x|=r\mathit{\}}$, we denote the sphere of radius r, while 𝕊=𝕊₁ stands for the unit sphere. Every unit vector $\mathit{\xi }=\mathit{\xi }(t,\mathit{\phi })\in \mathbb{S}$ can be expressed in spherical coordinates with longitude φ and polar distance t=cos(ϑ), where ϑ is the corresponding co-latitude. By Y_n,k, we denote fully normalized spherical harmonics of degree n and order k: for every n∈ℕ₀ and $k=-n,\mathrm{\dots },n$,

$${\displaystyle }{Y}_{n,k}\left(\mathit{\xi }\right)=\sqrt{{\displaystyle \frac{\mathrm{2}n+\mathrm{1}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{(n-|k\left|\right)\mathrm{!}}{(n+|k\left|\right)\mathrm{!}}}}{P}_n^{\left|k\right|}\left(t\right)\left\{\begin{array}{ll}\sqrt{\mathrm{2}}\cos \left(k\mathit{\phi }\right),& k<\mathrm{0},\\ \mathrm{1},& k=\mathrm{0},\\ \sqrt{\mathrm{2}}\sin \left(k\mathit{\phi }\right),& k>\mathrm{0}.\end{array}\right.$$

The involved associated Legendre functions are, for $t\in [-\mathrm{1},\mathrm{1}]$, defined as

$${\displaystyle }{P}_n^k\left(t\right)={\displaystyle \frac{(-\mathrm{1}{)}^k}{{\mathrm{2}}^{n}n\mathrm{!}}}{\left(\mathrm{1}-t^{\mathrm{2}}\right)}^{k/\mathrm{2}}{\left({\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\right)}^{n+k}(t^{\mathrm{2}}-\mathrm{1}{)}^n.$$

If f is a scalar-valued square-integrable function on 𝕊, then, for every degree n∈ℕ₀ and order $k=-n,\mathrm{\dots },n$, the values

$$\begin{array}{}\text{(3)}& {\displaystyle }{\widehat{f}}_{(n,k)}=\underset{\mathbb{S}}{\int }f\left(\mathit{\eta }\right)Y_{n,k}\left(\mathit{\eta }\right)\mathrm{d}S\left(\mathit{\eta }\right)\end{array}$$

are called the Fourier coefficients of the function f.

Going over to the vectorial setting, it is well known that every square-integrable vector field f on the unit sphere can be uniquely decomposed into its radial and two tangential components such that

$${\displaystyle }\mathit{f}={\mathit{e}}_r{f}_{\mathrm{1}}+{\mathrm{\nabla }}_{\mathbb{S}}{f}_{\mathrm{2}}+{\mathit{L}}_{\mathbb{S}}{f}_{\mathrm{3}},$$

with scalar-valued functions $f_{\mathrm{1}},f_{\mathrm{2}},f_{\mathrm{3}}$ and the radial unit vector ${\mathit{e}}_r={\left(\sqrt{\mathrm{1}-t^{\mathrm{2}}}\cos \left(\mathit{\phi }\right),\sqrt{\mathrm{1}-t^{\mathrm{2}}}\sin \left(\mathit{\phi }\right),t\right)}^T$. By the surface gradient ∇_𝕊, we denote the tangential component of the usual Euclidean gradient ∇, i.e.,

$${\displaystyle }{\mathrm{\nabla }}_{\mathbb{S}}={\mathit{e}}_{\mathit{\phi }}{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}-t^{\mathrm{2}}}}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}+{\mathit{e}}_t\sqrt{\mathrm{1}-t^{\mathrm{2}}}{\displaystyle \frac{\partial }{\partial t}},$$

with unit vectors ${\mathit{e}}_t={\left(-t\cos \left(\mathit{\phi }\right),-t\sin \left(\mathit{\phi }\right),\sqrt{\mathrm{1}-t^{\mathrm{2}}}\right)}^T$ and ${\mathit{e}}_{\mathit{\phi }}={\left(-\sin \left(\mathit{\phi }\right),\cos \left(\mathit{\phi }\right),\mathrm{0}\right)}^T$. Moreover, the surface curl gradient L_𝕊 is defined by ${\mathit{L}}_{\mathbb{S}}f\left(\mathit{\xi }\right)=\mathit{\xi }\times {\mathrm{\nabla }}_{\mathbb{S}}f\left(\mathit{\xi }\right)$, where × is the usual cross product in ℝ³. In other words,

$${\displaystyle }{\mathit{L}}_{\mathbb{S}}=-{\mathit{e}}_{\mathit{\phi }}\sqrt{\mathrm{1}-t^{\mathrm{2}}}{\displaystyle \frac{\partial }{\partial t}}+{\mathit{e}}_t{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}-t^{\mathrm{2}}}}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}.$$

Hence, we define three types of vector spherical harmonics: the radial

$${\displaystyle }{\mathit{y}}_{n,k}^{\left(\mathrm{1}\right)}{\displaystyle }={\mathit{e}}_r{Y}_{n,k},$$

for degrees n≥0 and orders $k=-n,\mathrm{\dots },n$, as well as the tangential

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{\mathit{y}}_{n,k}^{\left(\mathrm{2}\right)}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n(n+\mathrm{1})}}}{\mathrm{\nabla }}_{\mathbb{S}}{Y}_{n,k},\text{(5)}& {\displaystyle }& {\displaystyle }{\mathit{y}}_{n,k}^{\left(\mathrm{3}\right)}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n(n+\mathrm{1})}}}{\mathit{L}}_{\mathbb{S}}{Y}_{n,k},\end{array}$$

for degrees n≥1 and orders $k=-n,\mathrm{\dots },n$. Note that, for convenience, we set ${\mathit{y}}_{\mathrm{0},\mathrm{0}}^{\left(\mathrm{2}\right)}={\mathit{y}}_{\mathrm{0},\mathrm{0}}^{\left(\mathrm{3}\right)}=\mathrm{0}$. It should further be noted that the vector spherical harmonics in Eq. (4) are surface-curl-free while those in Eq. (5) are surface-divergence-free. In analogy to Eq. (3), we can now define the Fourier coefficients

$${\displaystyle }{\widehat{f}}_{(n,k)}^{\left(i\right)}=\underset{\mathbb{S}}{\int }\mathit{f}\left(\mathit{\eta }\right)\cdot {\mathit{y}}_{n,k}^{\left(i\right)}\left(\mathit{\eta }\right)\mathrm{d}S\left(\mathit{\eta }\right)$$

of square-integrable vector fields f on the unit sphere.

The vector spherical harmonics from above are defined solely on the unit sphere and can, therefore, only be used for the expansion of vector-valued functions on 𝕊. However, for the approximation of satellite potential field data it is necessary to have related functions that are also defined in the exterior of a sphere. For that purpose, we define the following gradients of harmonic extensions of (scalar) spherical harmonics:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{h}}_{n,k}\left(x\right)={\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\left({\displaystyle \frac{a}{r}}\right)}^n\left({\mathrm{\nabla }}_{\mathbb{S}}{Y}_{n,k}\left(\mathit{\xi }\right)-\mathit{\xi }(n+\mathrm{1})Y_{n,k}\left(\mathit{\xi }\right)\right),\end{array}$$

for $r=\left|x\right|>a$ and $\mathit{\xi }=\frac{x}{\left|x\right|}\in \mathbb{S}$, where a is the radius of a reference sphere, e.g., Earth's mean radius.

2.2 Vectorial Abel–Poisson kernel

While the set of functions $\mathit{\{}{\mathit{h}}_{n,k}{\mathit{\}}}_{n\in {\mathbb{N}}_{\mathrm{0}},k=-n,\mathrm{\dots },n}$ from Eq. (6) is suitable for the global approximation of potential field data, we are also interested in localized functions. One possible choice is the Abel–Poisson kernel (see, e.g., Freeden and Gerhards, 2012; Freeden et al., 1998). For $x,y\in {\mathbb{R}}^{\mathrm{3}}$, $\left|x\right|>\left|y\right|$, it is defined by

$${\displaystyle }K(x,y)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{|x{|}^{\mathrm{2}}-|y{|}^{\mathrm{2}}}{|x-y{|}^{\mathrm{3}}}}.$$

That is, with unit vectors $\mathit{\xi },\mathit{\eta }\in \mathbb{S}$ and radii $r>a>\mathrm{0}$, we have

$${\displaystyle }K(r\mathit{\xi },a\mathit{\eta })={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{r^{\mathrm{2}}-a^{\mathrm{2}}}{{\left(a^{\mathrm{2}}+r^{\mathrm{2}}-\mathrm{2}ar(\mathit{\xi }\cdot \mathit{\eta })\right)}^{\mathrm{3}/\mathrm{2}}}},$$

which shows that K only depends on the spherical distance between ξ and η, since $|\mathit{\xi }-\mathit{\eta }{|}^{\mathrm{2}}=\mathrm{2}(\mathrm{1}-\mathit{\xi }\cdot \mathit{\eta })$. Therefore, the kernel is radially symmetric if we keep one of the variables fixed (we strictly keep the second argument fixed, here aη). The degree of localization is determined by the ratio $\frac{a}{r}$. The closer it is to 1, i.e., the smaller the difference between the radii a and r, the better the spatial localization of $K(r\cdot ,a\mathit{\eta })$ around η. In our case, we choose a to be the Earth's mean radius, and r is the radius of the sphere at which we evaluate the kernel. An illustration of the kernel is provided in Fig. 1.

Figure 3 Absolute value of the trial function ${\mathit{B}}_{\mathrm{50}}^{\mathrm{re}}$ corresponding to $\mathit{u}={\mathit{g}}_{\mathrm{50}}^{\left(\mathrm{3}\right)}$ at time t=0 (a) and the trial function ${\mathit{B}}_{\mathrm{1630}}^{\mathrm{re}}$ corresponding to $\mathit{u}={\mathit{g}}_{\mathrm{1630}}^{\left(\mathrm{3}\right)}$ at time t=0 (b). Models of surface shell conductance on 𝕊_a and absolute value of B_main on 𝕊_a used for the generation of the trial functions (c, d).

The corresponding vectorial Abel–Poisson kernel is simply defined by

$${\displaystyle }\mathit{k}(x,y){\displaystyle }={\mathrm{\nabla }}_{x}K(x,y)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\sum _{n=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=-n}^n{Y}_{n,k}\left({\displaystyle \frac{y}{\left|y\right|}}\right)h_{n,k}\left(x\right).$$

Further calculations show

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathit{k}(x,y)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\left({\displaystyle \frac{\mathrm{2}}{|x-y{|}^{\mathrm{3}}}}x-\mathrm{3}{\displaystyle \frac{|x{|}^{\mathrm{2}}-|y{|}^{\mathrm{2}}}{|x-y{|}^{\mathrm{5}}}}(x-y)\right).\end{array}$$

A spatially localized alternative to $\mathit{\{}{\mathit{h}}_{n,k}{\mathit{\}}}_{n\in {\mathbb{N}}_{\mathrm{0}},k=-n,\mathrm{\dots },n}$ could then be defined by the set of functions $\mathit{\left\{}\mathit{k}\right(\cdot ,a{\mathit{\eta }}_i){\mathit{\}}}_{i=\mathrm{1},\mathrm{\dots },M}$, where ${\mathit{\eta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\eta }}_M\in \mathbb{S}$ is a fixed set of adequately distributed nodal points.

2.3 Spherical Slepian functions

While the localization of Abel–Poisson kernels is of radially symmetric nature, one is often interested in regions of more complex geometry, e.g., continents or oceans. Spherical Slepian functions, for instance, provide an orthonormal system of functions that can reflect localization in such general predefined regions Γ⊂𝕊 (see, e.g., Plattner and Simons, 2015, 2017a; Simons and Plattner, 2015; Simons et al., 2006 for details).

Specifically, the function f showing the best localization in Γ, is the one that maximizes the energy ratio

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{\Gamma }}\left(\mathit{f}\right)={\displaystyle \frac{{\int }_{\mathrm{\Gamma }}\left|\mathit{f}\right(\mathit{\eta }\left){|}^{\mathrm{2}}\mathrm{d}S\right(\mathit{\eta })}{{\int }_{\mathbb{S}}\left|\mathit{f}\right(\mathit{\eta }\left){|}^{\mathrm{2}}\mathrm{d}S\right(\mathit{\eta })}},\end{array}$$

i.e., the one with an energy ratio closest to 1. Let us now assume that g⁽ⁱ⁾ is a bandlimited vectorial function of type i with bandlimit N; i.e., it can be expanded as

$${\displaystyle }{\mathit{g}}^{\left(i\right)}=\sum _{n=\mathrm{0}}^N\sum _{k=-n}^n{\widehat{g}}_{(n,k)}^{\left(i\right)}{\mathit{y}}_{n,k}^{\left(i\right)}.$$

Further, the matrix $\mathbf{P}=\left({\mathbf{P}}_{(n,k),(m,j)}\right)\in {\mathbb{R}}^{(N+\mathrm{1}{)}^{\mathrm{2}}\times (N+\mathrm{1}{)}^{\mathrm{2}}}$ contains (properly sorted²) all of the appearing inner products

$${\displaystyle }{\mathbf{P}}_{(n,k),(m,j)}=\underset{\mathrm{\Gamma }}{\int }{\mathit{y}}_{n,k}^{\left(i\right)}\left(\mathit{\eta }\right)\cdot {\mathit{y}}_{m,j}^{\left(i\right)}\left(\mathit{\eta }\right)\mathrm{d}S\left(\mathit{\eta }\right)$$

and $\widehat{\mathit{g}}=({\widehat{g}}_{(n,k)}^{\left(i\right)}{)}^{\mathrm{T}}\in {\mathbb{R}}^{(N+\mathrm{1}{)}^{\mathrm{2}}}$, with $n=\mathrm{0},\mathrm{\dots },N$ and $k=-n,\mathrm{\dots },n$. If we now restrict ourselves to normalized functions g⁽ⁱ⁾ (i.e., ${\int }_{\mathbb{S}}\left|{\mathit{g}}^{\left(i\right)}\right(\mathit{\eta }\left){|}^{\mathrm{2}}\mathrm{d}S\right(\mathit{\eta })={\widehat{\mathit{g}}}^{\mathrm{T}}\widehat{\mathit{g}}=\mathrm{1}$), one obtains the simple expression ${\mathit{\lambda }}_{\mathrm{\Gamma }}\left({\mathit{g}}^{\left(i\right)}\right)={\widehat{\mathit{g}}}^{\mathrm{T}}\mathbf{P}\widehat{\mathit{g}}$. Eventually, the maximization of the energy ratio Eq. (8) leads to the eigenvalue problem

$$\mathbf{P}\widehat{\mathit{g}}=\mathit{\lambda }\widehat{\mathit{g}}.$$

The eigenvalues λ_ℓ are the possible energy ratios and the corresponding eigenvectors ${\widehat{\mathit{g}}}_{\mathrm{\ell }}$ contain the Fourier coefficients of bandlimited functions ${\mathit{g}}_{\mathrm{\ell }}^{\left(i\right)}$ attaining the energy ratio ${\mathit{\lambda }}_{\mathrm{\ell }}={\mathit{\lambda }}_{\mathrm{\Gamma }}\left({\mathit{g}}_{\mathrm{\ell }}^{\left(i\right)}\right)$. The set of functions $\mathit{\{}{\mathit{g}}_{\mathrm{\ell }}^{\left(i\right)}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },(N+\mathrm{1}{)}^{\mathrm{2}}}$ is ordered such that $\mathrm{1}\ge {\mathit{\lambda }}_{\mathrm{1}}\ge {\mathit{\lambda }}_{\mathrm{2}}\ge \mathrm{\dots }\ge {\mathit{\lambda }}_{(N+\mathrm{1}{)}^{\mathrm{2}}}\ge \mathrm{0}$.

Figure 4 Accumulated energy for $\mathit{\{}{\mathit{u}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$ (a) and $\mathit{\{}{\mathit{B}}_{\mathrm{\ell }}^{\mathrm{re}},{\mathit{B}}_{\mathrm{\ell }}^{\mathrm{im}}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$ (b), with L=1200.

In typical scenarios, it turns out that the eigenvalues are clustered close to 1 and close to zero. Those eigenvalues ${\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_L$ which are closer to 1 determine the subset $\mathit{\{}{\mathit{g}}_{\mathrm{\ell }}^{\left(i\right)}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$ of well-localized Slepian functions that should be used for approximation in Γ. The code for the generation of vectorial Slepian functions has been kindly supplied in Plattner and Simons (2017b). For our situation, where Γ denotes the region of (a spherical) Earth which is covered by oceans, an illustration is provided in Fig. 2.

Figure 5 Absolute value of the radial part of the tidal model ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$ as well as the forward model ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ at an altitude of 300 km above the Earth's surface.

2.4 Physics-based trial functions

We start with the time-harmonic Maxwell equations as already indicated in Eq. (1). For simplicity, we assume a 1-D (only radially varying) conductivity model for σ within the ball 𝔹_a, and at the surface 𝕊_a we allow a laterally varying conductivity (cf. the bottom left image in Fig. 3 for an illustration). Further, the magnetic field B_main is taken from the CHAOS-5 model (see Finlay et al., 2015) and u is supposed to denote a depth-integrated velocity field that is restricted to 𝕊_a (in fact, within the numerical framework of the X3DG solver, we assume a constant ocean depth of 1km, with u being tangential to the sphere and independent of the depth). Since we are mainly interested in tidal velocity fields, it is a reasonable assumption that u is surface-divergence-free for most parts of the oceans. The latter means that u can be expanded in terms of vector spherical harmonics or vectorial Slepian functions of type 3, i.e., ${\mathit{y}}_{n,k}^{\left(\mathrm{3}\right)}$ or ${\mathit{g}}_{\mathrm{\ell }}^{\left(\mathrm{3}\right)}$, respectively.

Figure 6 Absolute value of the radial part of approximations of ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$ based on dictionary 𝒟₁ (a), dictionary 𝒟₂ (b), and dictionary 𝒟₃ (c), as well as the corresponding residuals with respect to ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$ (d, e, f). Note the different scales in the bottom row which are chosen in order to emphasize the spatial distribution of the residuals.

For the generation of the tailored trial functions $\mathit{\{}{\mathit{B}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$, we therefore substitute u by a set of surface-divergence-free functions $\mathit{\{}{\mathit{u}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L}$ that reflect spatial localization within the oceans. More precisely, we choose

$${\mathit{u}}_{\mathrm{\ell }}={\mathit{g}}_{\mathrm{\ell }}^{\left(\mathrm{3}\right)},$$

where ${\mathit{g}}_{\mathrm{\ell }}^{\left(\mathrm{3}\right)}$ is the ℓth best localized vectorial Slepian function of type 3. The corresponding solution B_oc of Eq. (1) within this setup then provides an auxiliary function ${\stackrel{\mathrm{̃}}{\mathit{B}}}_{\mathrm{\ell }}$. It should be noted that in order to obtain Maxwell's equation in the time-harmonic form Eq. (1), one has to apply a Fourier transform in time. Therefore, for the actual trial function B_ℓ, we have to invert the Fourier transform and get

$${\displaystyle }{\mathit{B}}_{\mathrm{\ell }}(x,t)=e^{-i\mathit{\omega }t}{\stackrel{\mathrm{̃}}{\mathit{B}}}_{\mathrm{\ell }}\left(x\right),x\in {\mathbb{R}}^{\mathrm{3}},t\in \mathbb{R}.$$

For technical reasons, we choose to work in a real-valued framework, so that the real and imaginary part of B_ℓ each yield a trial function

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{B}}_{\mathrm{\ell }}^{\mathrm{re}}(x,t)=\cos \left(\mathit{\omega }t\right){\stackrel{\mathrm{̃}}{\mathit{B}}}_{\mathrm{\ell }}^{\mathrm{re}}\left(x\right)+\sin \left(\mathit{\omega }t\right){\stackrel{\mathrm{̃}}{\mathit{B}}}_{\mathrm{\ell }}^{\mathrm{im}}\left(x\right),\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{B}}_{\mathrm{\ell }}^{\mathrm{im}}(x,t)=\sin \left(\mathit{\omega }t\right){\stackrel{\mathrm{̃}}{\mathit{B}}}_{\mathrm{\ell }}^{\mathrm{re}}\left(x\right)-\cos \left(\mathit{\omega }t\right){\stackrel{\mathrm{̃}}{\mathit{B}}}_{\mathrm{\ell }}^{\mathrm{im}}\left(x\right).\end{array}$$

Thus, each choice of u_ℓ yields two functions ${\mathit{B}}_{\mathrm{\ell }}^{\mathrm{re}}$ and ${\mathit{B}}_{\mathrm{\ell }}^{\mathrm{im}}$ that reflect the temporal periodicity of the tidal magnetic field as well as the spatial localization of the sources within the oceans. An illustration for the M2 tide with $\mathit{\omega }=\frac{\mathrm{2}\mathit{\pi }}{\mathrm{12.42}h}$ can be found in Fig. 3. For the computation of the ${\stackrel{\mathrm{̃}}{\mathit{B}}}_{\mathrm{\ell }}$ as solutions of Eq. (1), we have used the X3DG solver from Kuvshinov (2008).

Figure 7 Absolute value of the radial part of approximations of ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ based on dictionary 𝒟₁ (a), dictionary 𝒟₂ (b), and dictionary 𝒟₃ (c), as well as the corresponding residuals with respect to ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ (d, e, f).

Figure 4 shows the accumulated energy ${\sum }_{\mathrm{\ell }=\mathrm{1}}^L\left|{\mathit{u}}_{\mathrm{\ell }}\right(\mathit{\xi }){|}^{\mathrm{2}}$, for ξ∈𝕊, of the underlying functions u_ℓ that describe the velocity field and the accumulated energy ${\sum }_{\mathrm{\ell }=\mathrm{1}}^L\left|{\mathit{B}}_{\mathrm{\ell }}^{\mathrm{re}}\right(x,t){|}^{\mathrm{2}}+|{\mathit{B}}_{\mathrm{\ell }}^{\mathrm{im}}(x,t){|}^{\mathrm{2}}$, for x∈𝕊_r with $r=a+\mathrm{300}$ km and time t=0, of the corresponding trial functions. In both cases, one can clearly see the spatial localization over the oceans. However, the accumulated energy of the trial functions additionally reflects the influence of the conductivity σ and the main/core magnetic field B_main indicated in Fig. 3.

Figure 8 Absolute value of the radial part of ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ (a) and model of continental “noise” e (b), as well as superposition ${\mathit{B}}_{\mathrm{oc}}^{\mathit{e}}$ of both of these (c).

3 Examples

For our experiments we rely on the CM5 geomagnetic field model (cf. Sabaka et al., 2015) and a forward model based on the M2 depth-integrated tidal velocity field from TPXO8-ATLAS (cf. Egbert and Erofeeva, 2002) that has also been used in Kuvshinov (2008). The contribution of CM5 that is due to the oceanic M2 tide is given as an expansion in terms of spherical harmonics up to degree 18; we denote it as ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$ for the remainder of this section and sample it at M=250 000 points which are taken from actual Swarm satellite tracks. The forward model has been computed via the X3DG solver based on the surface conductance and the main/core magnetic field model indicated in the bottom row of Fig. 3 and a depth-integrated M2 tidal velocity field from TPXO8-ATLAS. We denote it by ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ and evaluate it on the same point grid as before. These samples are used as input data b∈ℝ^M for the Regularized (Orthogonal) Functional Matching Pursuit, which works iteratively as indicated in Eq. (2). In the following, we want to illustrate the influence of the choice of different function systems (i.e., the choice of different dictionaries 𝒟) on the approximation of ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$ and ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$. For that purpose, we choose three different dictionaries: the spherical-harmonic-based

$${\displaystyle }{\mathcal{D}}_{\mathrm{1}}=\mathit{\left\{}\cos \right(\mathit{\omega }t\left){\mathit{h}}_{n,k}\right(x),\sin (\mathit{\omega }t\left){\mathit{h}}_{n,k}\right(x){\mathit{\}}}_{n=\mathrm{0},\mathrm{\dots },\mathrm{20},k=-n,\mathrm{\dots },n},$$

with $\mathit{\omega }=\frac{\mathrm{2}\mathit{\pi }}{\mathrm{12.42}h}$ and the functions h_n,k from Eq. (6); the Abel–Poisson-kernel-based

$${\displaystyle }{\mathcal{D}}_{\mathrm{2}}=\mathit{\left\{}\cos \right(\mathit{\omega }t\left)\mathit{k}\right(x,a{\mathit{\eta }}_i),\sin (\mathit{\omega }t\left)\mathit{k}\right(x,a{\mathit{\eta }}_i){\mathit{\}}}_{i=\mathrm{1},\mathrm{\dots },M_p},$$

where a=6371.2 km, $\mathit{\{}{\mathit{\eta }}_i{\mathit{\}}}_{i=\mathrm{1},\mathrm{\dots },M_p}$ is a Reuter grid on 𝕊 with M_p=6201 nearly equally distributed points (see, e.g., Michel (2013), p. 137), and k given as in Eq. (7); and

$${\displaystyle }{\mathcal{D}}_{\mathrm{3}}=\mathit{\left\{}{\mathit{B}}_{\mathrm{\ell }}^{\mathrm{re}}\right(x,t),{\mathit{B}}_{\mathrm{\ell }}^{\mathrm{im}}(x,t){\mathit{\}}}_{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },\mathrm{1200}},$$

with ${\mathit{B}}_{\mathrm{\ell }}^{\mathrm{re}}$ and ${\mathit{B}}_{\mathrm{\ell }}^{\mathrm{im}}$, the physics-based trial functions from Eqs. (9) and (10).

Figure 9 Absolute value of the radial part of approximations of superposition ${\mathit{B}}_{\mathrm{oc}}^{\mathit{e}}$ based on dictionary 𝒟₁ (a), dictionary 𝒟₂ (b), and dictionary 𝒟₃ (c), as well as the corresponding residuals with respect to original ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ without the continental “noise” e (d, e, f).

The actual signals that we want to approximate are indicated in Fig. 5. The approximations ${\stackrel{\mathrm{‾}}{\mathit{B}}}_N$ of ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$ together with the residuals $|{\stackrel{\mathrm{‾}}{\mathit{B}}}_N-{\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}|$ for each of the three dictionaries above are shown in Fig. 6, whereas the respective approximations of ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ and corresponding residuals $|{\stackrel{\mathrm{‾}}{\mathit{B}}}_N-{\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}|$ are displayed in Fig. 7.

Figure 10 The difference between the approximations ${\stackrel{\mathrm{‾}}{\mathit{B}}}_N$ of the undisturbed ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ (as shown in Fig. 7a, b, and c) and the approximations ${\stackrel{\mathrm{‾}}{\mathit{B}}}_N^{\mathit{e}}$ of the noisy ${\mathit{B}}_{\mathrm{oc}}^{\mathit{e}}$ (as shown in Fig. 9a, b, and c).

Table 1 Root mean square errors (RMSEs) corresponding to the approximations of the undisturbed forward model ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ as well as the noisy model ${\mathit{B}}_{\mathrm{oc}}^{\mathit{e}}={\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}+\mathit{e}$ compared to the “ground truth” ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ with the three different dictionaries. Left-hand columns show errors for the approximation of ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ (compare Fig. 7), while center columns show errors for the approximation of ${\mathit{B}}_{\mathrm{oc}}^{\mathit{e}}$ (compare Fig. 9). The right-hand columns display the RMSEs of the difference between the respective approximations (as shown in Fig. 10).

In the case of ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$ as the underlying signal, it can be seen in Fig. 6 that the dictionary 𝒟₁ yields the overall best approximation, which, however, is not surprising since we try to fit a spherical-harmonic-based model with a spherical-harmonic-based dictionary. The result for dictionary 𝒟₂ shows a more localized pattern in the residual, as is expected for the use of Abel–Poisson kernels. However, the maxima in the residual are not correlated to specific continental or oceanic structures but they mainly coincide with the maxima of the original signal ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$. The situation for dictionary 𝒟₃ of physics-based trial functions is different. The agreement of the approximation with ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$ is good over the oceans but significant deviations exist over the continents. The latter could be an indication that the original model ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{CM}\mathrm{5}}$ contains contributions over the continents whose physical origin is not due to induction by oceanic tides. Some smaller deviations over oceanic areas exist around southern Africa and the east of Australia. Since we are dealing with the approximation of a low-degree (up to degree 18) spherical-harmonic-based M2 tidal magnetic field model by localized trial functions, one cannot reliably say if the latter deviations are artifacts from the approximation procedure or if they have a physical origin. However, those are areas with a shallower ocean topography, so the assumption of surface-divergence-free depth-integrated tidal velocities (which we made for our choice of the underlying u_ℓ) and the assumption of a constant ocean depth (we chose a depth of 1 km for the generation of the ${\stackrel{\mathrm{̃}}{\mathit{B}}}_{\mathrm{\ell }}$ via the X3DG solver) might not be accurate in these areas. Nonetheless, the residuals over the continents show that the use of the adapted trial functions might eventually deliver improved tidal magnetic field models that correct unrealistic continental contributions without disregarding continental areas entirely.

The residuals of the approximations of the forward model ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ in Fig. 7, on the other hand, indicate that the quality of the approximations does not vary too much (at least on scales that are relevant for satellite data approximation) among the three tested function systems. This is mainly due to the fact that the input model ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ already reflects certain spatial localization properties over the oceans. In such a scenario (if additionally solely interested in the approximation of the signal and not the underlying velocity fields) it would, therefore, not be necessary to use the adapted trial functions that we have introduced. The crucial point, however, is that satellite data typically contain undesired contributions over the continents that are not due to ocean-tide-generated magnetic fields. In order to illustrate this behavior, we use the following additional example. We take randomly distributed Fourier coefficients to construct a “noise” function e with a bandlimit of degree 40, i.e.,

$$\mathit{e}=\sum _{n=\mathrm{1}}^{\mathrm{40}}\sum _{k=-n}^n{\widehat{e}}_{(n,k)}{\mathit{y}}_{n,k}^{\left(\mathrm{3}\right)},$$

where the Fourier coefficients ${\widehat{e}}_{(n,k)}$ are normally distributed with zero mean and a variance such that the amplitude of e is in the range of the oceanic signal ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$. This function e is then restricted to the continents and eventually superposed with the forward model (cf. Fig. 8). For the sake of clarity, we denote the approximations of the noisy data

$${\mathit{B}}_{\mathrm{oc}}^{\mathit{e}}={\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}+\mathit{e}$$

by ${\stackrel{\mathrm{‾}}{\mathit{B}}}_N^{\mathit{e}}$ instead of ${\stackrel{\mathrm{‾}}{\mathit{B}}}_N$. The latter still represents the approximation of ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ without extra continental noise e.

In Fig. 9 one can directly see the influence which the continental noise has on the approximation depending on the various dictionaries (Fig. 7 shows the same quantities for the approximations in the undisturbed setup). In the case of dictionary 𝒟₁, the spherical harmonics also approximate a part of the continental data which in turn also has some impact on the approximation in oceanic areas. Due to the localization of the kernels contained in dictionary 𝒟₂, the (undesired) reconstruction of the continental noise is even more accurate, while the reconstruction over the oceans only changes very slightly. With the proposed physics-based functions in dictionary 𝒟₃, however, the influence of the continental noise is much less apparent. The maxima occur very close to the coastline, which is most likely due to numerical issues stemming from discontinuities of the data ${\mathit{B}}_{\mathrm{oc}}^{\mathit{e}}$ in coastal areas. A closer look at the differences $|{\stackrel{\mathrm{‾}}{\mathit{B}}}_N-{\stackrel{\mathrm{‾}}{\mathit{B}}}_N^{\mathit{e}}|$ between the approximations of noisy and undisturbed data is given in Fig. 10. This shows again that the inclusion of continental noise has a smaller effect on the approximation via physics-based trial functions than on the approximations via the other tested trial functions. Moreover, the corresponding root mean square errors of the approximations ${\stackrel{\mathrm{‾}}{\mathit{B}}}_N$ and ${\stackrel{\mathrm{‾}}{\mathit{B}}}_N^{\mathit{e}}$, respectively, can be found in Table 1. In both cases, we compared the approximations to the undisturbed data ${\mathit{B}}_{\mathrm{oc}}^{\mathrm{X}\mathrm{3}\mathrm{DG}}$ in order to emphasize the impact of continental noise on the overall approximation. The errors over continental and oceanic regions are provided separately.

4 Conclusions

The main goal of this paper is to study the errors that are made by the approximation of tidal magnetic fields by use of different sets of trial functions. While, e.g., Saynisch et al. (2018) compared forward models for the M2 tidal magnetic field based on different tidal models, we aim at illustrating the effect of the involved trial functions on the possible extraction of the tidal magnetic field from satellite data in the first place. The indicated residuals for the synthetic examples show that the use of the presented adapted physics-based trial functions could have a detectable effect for the extraction of such signals in satellite data. These trial functions reflect the underlying physics in the sense that they satisfy the time-harmonic Maxwell equations and that they include knowledge of the ambient core magnetic field and the Earth's conductivity, but they are not designed to rely on detailed oceanographic models. The latter can be advantageous since the extracted magnetic field induced by ocean tides might eventually be used to make inferences on such models.