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 non-localized trial functions are 5 illustrated.


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., 2015Sabaka et al., , 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 re-cently 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 R. Telschow et al.: Approximation of spatial structures of tidal magnetic fields given by the time-harmonic Maxwell equations: where we assume to have knowledge of B main , the underlying conductivity σ , and the frequency ω.Furthermore, E denotes the electric field and µ 0 is the vacuum permeability.Instead of using a fixed velocity field u, we substitute it by a set of functions {u } =1,...,L (e.g., vectorial Slepian functions as in Plattner andSimons (2014, 2015) that are localized over the oceans) to obtain a set of corresponding trial functions {B } =1,...,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 {B } =1,...,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 {u } =1,...,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 {B } =1,...,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.

Method and function systems
Given a so-called dictionary D of trial functions, we use the Regularized (Orthogonal) Functional Matching Pursuit (cf.Fischer and Michel, 2013;Michel andTelschow, 2014, 2016 for details) for the approximation of B oc .Briefly speaking, this is a greedy-type algorithm that yields an approximation (2) denotes the residual between the data b ∈ R M and the approximation after i − 1 iterations.In this particular setup, F represents the linear operator that evaluates a function at the M locations where data are provided, and H is a suitable Hilbert space for the regularization of the problem 1 .The parameter λ controls the trade-off between the data misfit R i 2 R M and the regularizing term B i 2 H , which imposes a certain property to 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 D (or combinations of such) from which the approximant 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 D. In particular, Sect.2.4 describes in more detail the aforementioned trial functions {B } =1,...,L that contain temporal and spatial constraints tailored for ocean-tide-induced magnetic fields.

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 S r = {x ∈ R 3 : |x| = r}, we denote the sphere of radius r, while S = S 1 stands for the unit sphere.Every unit vector ξ = ξ(t, ϕ) ∈ 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 1 In our case, we use the norm ) , but other norms can be used as well depending on the property one wants to impose on f .By f(n,k) we denote the Fourier coefficients of f as indicated in Eq. ( 3).
The involved associated Legendre functions are, for t ∈ [−1, 1], defined as If f is a scalar-valued square-integrable function on S, then, for every degree n ∈ N 0 and order k = −n, . .., n, the values 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 with scalar-valued functions f 1 , f 2 , f 3 and the radial unit . By the surface gradient ∇ S , we denote the tangential component of the usual Euclidean gradient ∇, i.e., with unit vectors e t = −t cos(ϕ), −t sin(ϕ), T and e ϕ = (− sin(ϕ), cos(ϕ), 0) T .Moreover, the surface curl gradient L S is defined by , where × is the usual cross product in R 3 .In other words, Hence, we define three types of vector spherical harmonics: the radial for degrees n ≥ 0 and orders k = −n, . .., n, as well as the tangential for degrees n ≥ 1 and orders k = −n, . .., n.Note that, for convenience, we set y (2) 0,0 = y (3) 0,0 = 0.It should further be noted that the vector spherical harmonics in Eq. ( 4) are surface-curl-free while those in Eq. ( 5) are surfacedivergence-free.In analogy to Eq. (3), we can now define the Fourier coefficients 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 S.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: for r = |x| > a and ξ = x |x| ∈ S, where a is the radius of a reference sphere, e.g., Earth's mean radius.

Vectorial Abel-Poisson kernel
While the set of functions {h n,k } n∈N 0 ,k=−n,...,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) That is, with unit vectors ξ, η ∈ S and radii r > a > 0, we have which shows that K only depends on the spherical distance between ξ and η, since |ξ − η|2 = 2(1 − ξ • η).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 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•, aη) 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.
The corresponding vectorial Abel-Poisson kernel is simply defined by Further calculations show A spatially localized alternative to {h n,k } n∈N 0 ,k=−n,...,n could then be defined by the set of functions {k(•, aη i )} i=1,...,M , where η 1 , . .., η M ∈ S is a fixed set of adequately distributed nodal points.

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 ⊂ S (see, e.g., Plattner andSimons, 2015, 2017a;Simons et al., 2006;Simons and Plattner, 2015 for details).
Specifically, the function f showing the best localization in , is the one that maximizes the energy ratio i.e., the one with an energy ratio closest to 1. Let us now assume that g (i) is a bandlimited vectorial function of type i with bandlimit N; i.e., it can be expanded as Further, the matrix P = (P (n,k),(m,j ) ) ∈ R (N +1) 2 ×(N +1) 2 contains (properly sorted 2 ) all of the appearing inner products ) 2 , with n = 0, . .., N and k = −n, . .., n.If we now restrict ourselves to normalized functions g (i) (i.e., S |g (i) (η)| 2 dS(η) = ĝT ĝ = 1), one obtains the simple expression λ (g (i) ) = ĝT P ĝ.Eventually, the maximization of the energy ratio Eq. ( 8) leads to the eigenvalue problem The eigenvalues λ are the possible energy ratios and the corresponding eigenvectors ĝ contain the Fourier coefficients of bandlimited functions g (i) attaining the energy ratio i) .The set of functions {g In typical scenarios, it turns out that the eigenvalues are clustered close to 1 and close to zero.Those eigenvalues λ 1 , . .., λ L which are closer to 1 determine the subset {g (i) } =1,...,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.

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 B a , and at the surface S 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 S 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., y n,k or g (3) , respectively.
For the generation of the tailored trial functions {B } =1,...,L , we therefore substitute u by a set of surfacedivergence-free functions {u } =1,...,L that reflect spatial lo-calization within the oceans.More precisely, we choose u = g ( 3) , where g 3) 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 B .It should be noted that in order to obtain Maxwell's equation in the timeharmonic 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 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 Thus, each choice of u yields two functions B re and B 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 ω = 2π 12.42h can be found in Fig. 3.For the computation of the B as solutions of Eq. ( 1), we have used the X3DG solver from Kuvshinov (2008).
Figure 4 shows the accumulated energy L =1 |u (ξ )| 2 , for ξ ∈ S, of the underlying functions u that describe the velocity field and the accumulated energy L =1 |B re (x, t)| 2 + |B im (x, t)| 2 , for x ∈ S r with r = a + 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.

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 B CM5 oc 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 B X3DG oc and evaluate it on the same point grid as before.These samples are used as input data b ∈ R 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 D) on the approximation of B CM5 oc and B X3DG oc .For that purpose, we choose three different dictionaries: the spherical-harmonic-based  with ω = 2π 12.42h and the functions h n,k from Eq. ( 6); the Abel-Poisson-kernel-based D 2 = {cos(ωt)k(x, aη i ), sin(ωt)k(x, aη i )} i=1,...,M p , where a = 6371.2km, {η i } i=1,...,M p is a Reuter grid on S with M p = 6201 nearly equally distributed points (see, e.g., Michel (2013), p. 137), and k given as in Eq. ( 7); and with B re and B im , the physics-based trial functions from Eqs. ( 9) and ( 10).
The actual signals that we want to approximate are indicated in Fig. 5 In the case of B CM5 oc as the underlying signal, it can be seen in Fig. 6 that the dictionary D 1 yields the overall best approximation, which, however, is not surprising since we try to fit a spherical-harmonic-based model with a spherical-harmonicbased dictionary.The result for dictionary D 2 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 B CM5 oc .The situation for dictionary D 3 of physics-based trial functions is different.The agreement of the approximation with B CM5 oc is good over the oceans but significant deviations exist over the continents.The latter could be an indication that the original model B CM5 oc 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 B via the X3DG solver) might not be accurate in these areas.Nonetheless, the resid-  uals 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 B X3DG oc 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 B X3DG oc 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 where the Fourier coefficients ê(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 B X3DG oc .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 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 D 1 , 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 D 2 , the (undesired) with the three different dictionaries.Left-hand columns show errors for the approximation of B X3DG oc (compare Fig. 7), while center columns show errors for the approximation of B e oc (compare Fig. 9).The right-hand columns display the RMSEs of the difference between the respective approximations (as shown in Fig. 10).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 B N and B e N , respectively, can be found in Table 1.In both cases, we compared the approximations to the undisturbed data B X3DG oc in order to emphasize the impact of continental noise on the overall approximation.The errors over continental and oceanic regions are provided separately.

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.
Data availability.The authors would like to thank the following providers of software and models: Alexey Kuvshinov for the code of the X3DG solver and a processed version of the depth-integrated tidal ocean velocities from TPXO8-ATLAS, Nils Olsen for the coefficients of the M2 tidal contribution in the CM5 model, Alain Plattner for the code for the generation of vectorial Slepian functions.
Competing interests.The authors declare that they have no conflict of interest.Special issue statement.This article is part of the special issue "Dynamics and interaction of processes in the Earth and its space environment: the perspective from low Earth orbiting satellites and beyond".It is not associated with a conference.

Figure 1 .
Figure 1.The kernel K(r•, aη 1 ) for a r = 0.91 (a) and a r = 0.67 (b).The fixed nodal point η 1 ∈ S is marked by a white cross.

Figure 2 .
Figure 2. Absolute value of the vectorial Slepian function g (3) 50 with 50th best localization over the oceans (a) and g (3) 1630 with the 50th worst localization over the oceans (b), for bandlimit N = 40.

Figure 3 .
Figure 3. Absolute value of the trial function B re 50 corresponding to u = g (3) 50 at time t = 0 (a) and the trial function B re 1630 corresponding to u = g (3) 1630 at time t = 0 (b).Models of surface shell conductance on S a and absolute value of B main on S a used for the generation of the trial functions (c, d).

Figure 5 .
Figure 5. Absolute value of the radial part of the tidal model B CM5 oc as well as the forward model B X3DG oc at an altitude of 300 km above the Earth's surface.

Figure 6 .
Figure 6.Absolute value of the radial part of approximations of B CM5 oc based on dictionary D 1 (a), dictionary D 2 (b), and dictionary D 3 (c), as well as the corresponding residuals with respect to B CM5 oc (d, e, f).Note the different scales in the bottom row which are chosen in order to emphasize the spatial distribution of the residuals.

Figure 7 .
Figure 7. Absolute value of the radial part of approximations of B X3DG oc based on dictionary D 1 (a), dictionary D 2 (b), and dictionary D 3 (c), as well as the corresponding residuals with respect to B X3DG oc

Figure 8 .
Figure 8. Absolute value of the radial part of B X3DG oc (a) and model of continental "noise" e (b), as well as superposition B e oc of both of these (c).
. The approximations B N of B CM5 oc together with the residuals |B N − B CM5 oc | for each of the three dictionaries above are shown in Fig. 6, whereas the respective approximations of B X3DG oc and corresponding residuals |B N − B X3DG oc | are displayed in Fig. 7.

Figure 9 .
Figure 9. Absolute value of the radial part of approximations of superposition B e oc based on dictionary D 1 (a), dictionary D 2 (b), and dictionary D 3 (c), as well as the corresponding residuals with respect to original B X3DG oc

Figure 10 .
Figure 10.The difference between the approximations B N of the undisturbed B X3DG oc (as shown in Fig. 7a, b, and c) and the approximations B e N of the noisy B e oc (as shown in Fig. 9a, b, and c).

N
instead of B N .The latter still represents the approximation of B X3DG oc without extra continental noise e.

Table 1 .
Root mean square errors (RMSEs) corresponding to the approximations of the undisturbed forward model B X3DG With the proposed physics-based functions in dictionary D 3 , 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 B e oc in coastal areas.A closer look at the differences |B N −B e N | between the approximations of noisy and undisturbed data is given in Fig.