Investigation of regional variation in core flow models using spherical Slepian functions

Abstract

By assuming that changes in the magnetic field in the Earth’s outer core are advection-dominated on short timescales, models of the core surface flow can be deduced from secular variation. Such models are known to be under-determined and thus require other assumptions to produce feasible flows. There are regions where poor knowledge of the core flow dynamics gives rise to further uncertainty, such as within the tangent cylinder, and assumptions about the nature of the flow may lead to ambiguous patches, such as if it is assumed to be strongly tangentially geostrophic. We use spherical Slepian functions to spatially and spectrally separate core flow models, confining the flow to either inside or outside these regions of interest. In each region we examine the properties of the flow and analyze its contribution to the overall model. We use three forms of flow model: (a) synthetic models from randomly generated coefficients with blue, red and white energy spectra, (b) a snapshot of a numerical geodynamo simulation and (c) a model inverted from satellite magnetic field measurements. We find that the Slepian decomposition generates unwanted spatial leakage which partially obscures flow in the region of interest, particularly along the boundaries. Possible reasons for this include the use of spherical Slepian functions to decompose a scalar quantity that is then differentiated to give the vector function of interest, and the spectral frequency content of the models. These results will guide subsequent investigation of flow within localized regions, including applying vector Slepian decomposition methods.

Introduction

Over 95% of the geomagnetic field observed at Earth’s surface arises from the Earth’s fluid outer core (e.g., Jacobs 1987). The main field is generated by the geodynamo which converts fluid motion of the liquid iron into electric and magnetic energy. The geodynamo is a complex and self-sustaining magnetohydrodynamic system which is still relatively poorly understood in detail (e.g., Jones 2015). Several geomagnetic issues contribute to this limited knowledge, including the loss of resolution due to upward continuation of the magnetic field from the core–mantle boundary to the Earth’s surface, the masking effects of the intervening weakly conductive mantle and the imposition of the signals from the shallow magnetized crust (e.g., Hide 1969; Holme 1998). Due to these restrictions we are unable to image the small spatial scales and the very rapid temporal changes of the core field. Though main field models are continuously improving, thanks to better datasets from both ground-based observatories and satellite missions, there is an inherent limit to the achievable spatial accuracy of the core field around spherical harmonic degree and order 20, as observed at the surface (Hulot et al. 2009).

Large-scale changes in the magnetic field at the Earth’s surface over the period of months to years (termed secular variation or SV) can be used as a ‘tracer’ for the flow of the liquid at the core–mantle boundary (CMB), by assuming the conductivity of the core fluid is sufficiently high such that diffusion is negligible on periods of years to decades; this is known as the ‘frozen-flux’ hypothesis (Roberts and Scott 1965). By assuming frozen-flux, the SV can be inverted for the flow along the CMB. This allows us to probe some of the dynamical properties of the core, though at the cost of requiring additional assumptions about the structure of the flow itself. This is in order to reduce the ambiguity involved in inverting for two parameters (the eastward and northward flow components) from a single set of observations of the radial SV component (e.g., Backus 1968).

Fig. 1

Tangentially geostrophic ambiguous patches used in decomposition. Contours of $B_r/\cos \theta$ for values less than $10^6\hbox { nT}$. The shaded regions are the area considered to be ambiguous when the tangentially geostrophic assumption is applied. These regions were defined by considering the locations of closed contours which do not connect to the equator. The $B_r$ values are from the 2015 IGRF field model (Thébault et al. 2015). This and all other figures in this projection are centered on longitude $90^{\circ }$

Fig. 2

Location of the tangent cylinder and tangentially geostrophic ambiguous patches. The tangent cylinder can be described by two caps subtending an angle of $21^\circ$ at the Earth’s center (red) projecting the inner core to the surface of the core–mantle boundary (green) along the rotation axis. A schematic illustration of the tangentially geostrophic ambiguous patches is also shown in blue

Fig. 3

Decomposition of three synthetic spectra. Three synthetic flow models, with different energy spectra, decomposed using the Shannon number ($K=1504$) for the ambiguous patches. The top row (left: increasing; middle: flat; right: decreasing) shows maps of the azimuthal ($u \phi$) component of the toroidal flow. The second and third rows illustrate the spatial decomposition. The fourth row shows the summation of the decomposed parts of the flow model, which matches the first row. The bottom row shows the spectral decomposition. Note the different magnitudes of aliasing along the boundaries of each separated region depending on shape of spectra

Fig. 4

Decomposition of the numerical dynamo flow simulation in eastward direction ($u_{\phi}$ component). Decomposition of a snapshot flow from a numerical model with coefficients up to degree 60 (3721 coefficients) using the Shannon number ($K= 1504$) to represent the flow within the ambiguous patches

Fig. 5

Decomposition of numerical dynamo flow simulation in southward direction ($u_{\theta}$). Decomposition of a snapshot flow from a numerical model using the coefficients up to degree 60 (3721 coefficients) using the Shannon number ($K= 1504$) to represent the flow within the ambiguous patches

Here we localize such flow models, represented by spherical harmonic coefficients, to a particular region of the CMB surface by decomposing them with scalar spherical Slepian functions (Simons and Dahlen 2006). Spherical Slepian functions have been applied to a number of geophysical and astronomical problems where only partial or noisy datasets covering a fraction of the surface of a sphere are available. The functions provide a best estimate for incomplete spatio-spectral signals, offering an optimal trade-off between spatial and spectral leakage (Simons et al. 2006). They have been used for determining local gravitational changes arising from large earthquakes, the spectral structure of the cosmic microwave background radiation, the study of geodetic datasets containing temporal gaps, and to investigate the global crustal magnetic field (Beggan et al. 2013; Dahlen and Simons 2008; Harig and Simons 2012; Kim and von Frese 2017; Simons et al. 2009).

These studies have all successfully retrieved useful information about the geophysical system under consideration often through operation on a global spherical harmonic model. For example, Harig and Simons (2015, 2016) examine relatively small polar regions within the GRACE and GOCE global gravity models to deduce ice-loss over time. We note there are alternative methods for isolating regions of a spherical harmonic model on a spherical surface, such as spherical cap harmonics or localized basis functions which have been used for modeling the main magnetic field on Earth and other planets (Lesur 2006; Thébault et al. 2006, 2018). Spherical Slepian functions are the only functions that achieve regional separation in a fully analytical, and easily computable framework, depending only on the geometry of the region (Beggan et al. 2013).The double orthogonality of the Slepian functions, over the region of interest and the sphere, is a property that is convenient and very welcome on statistical grounds, for example, when inversions for the source or estimations of the power spectral density of the field components or the overall potential are being made on the basis of actual satellite data (Dahlen and Simons 2008; Plattner and Simons 2013; Simons et al. 2006).

Fig. 6

Varying K number of functions representing ‘inside’ the patches for the numerical dynamo flow simulation decomposition. Decomposition of the numerical dynamo model snapshot into ‘inside’ and ‘outside’ patch. The K value is the number of functions included inside the patches and $3721-(K+1)$ is the number of functions outside the patches. The $\phi$-component of the toroidal flow is shown

Fig. 7

Decomposition of inverted flow from satellite data in eastward direction ($u_{\phi}$). Decomposition of steady flow model obtained from satellite data (up to degree 20) using an ‘optimum’ value of K = 199 (coefficients out of 441) to represent flow within the patch. The Shannon number ($K=178$) produced a similar decomposition, but the unwanted boundary signals were slightly larger compared to the optimum K decomposition, which is shown in this image

Fig. 8

Decomposition of inverted flow from satellite data in southward direction ($u_{\theta}$). Decomposition of steady flow model obtained from satellite data (up to degree 20) using an ‘optimum’ value of $K = 199$ (coefficients out of 441) to represent flow within the patches. The Shannon number ($K=178$) produced a similar decomposition, but the unwanted aliasing was slightly larger compared to our chosen ‘optimum’ K decomposition, which is shown in this image

Fig. 9

Example of coefficient splitting from tangentially geostrophic ambiguous patch for inverted satellite data flows. An example of coefficient values for the input, inside tangentially geostrophic ambiguous patches, the complementary outside region and the summed coefficients of the decomposition for the Shannon number of the inverted satellite data. The input and summed values are exactly the same for each coefficient, but there are some coefficients which separate into much larger absolute values than their corresponding input

Fig. 10

Decomposition of numerical dynamo simulation only using coefficients that are larger than the input value. Flow maps of the decomposition shown in Fig. 7 when only the absolute value of each coefficients that are bigger than the input coefficients is plotted. There is no longer a distinction between inside and outside the region of interest, and the flow is concentrated in strong bands along the boundaries of the regions. The $\phi$-component of the toroidal flow is shown

Fig. 11

Decomposition of inverted flow from satellite data plotted as velocity vectors, stream function and poloidal potential. A decomposition of the toroidal and poloidal components of the inverted flow from satellite data, seen in Figs. 7 and 8, as velocity (column one), stream function (column two) and poloidal potential (third column). The input, top row, is split into inside and outside the tangentially geostrophic ambiguous patches, rows two and three, before being summed together in the final row. As with previous decompositions, the input and the summed decomposition are identical and flow is mostly restricted to the region of interest with the decomposed maps

Based on the success of the aforementioned studies, we apply similar techniques to a series of global core surface flow models. This paper focuses on the use of scalar spherical Slepian functions to decompose the flow models in order to investigate how localizing flow within and outside the tangentially geostrophic ambiguous regions (Backus and Le Mouël 1986) affects the overall spatial and spectral distribution of the flow. In Additional file 1, we also show the flow localized to the surface of the cylinder tangent to the inner core and in its complement. It is believed that the cylinder tangent to the inner core has a different flow regime (compared to the lower latitude parts of the outer core) where velocities are significantly enhanced locally (Jones 2015; Livermore et al. 2016).

The motivation for this study was to act as a preliminary investigation into the application of scalar spherical Slepian functions to outer core surface flow. Localizing the scalar potentials of the flow, rather than the vector flow itself, simplifies the matrix multiplications to perform the separation. These investigations should inform subsequent inversions of SV data for flows over localized patches of the CMB. We had hoped that forecasts of the geomagnetic field based on flow advection (Beggan and Whaler 2010) could be improved by using flows confined to regions where the dynamics are ‘known’ or the flow is less affected by non-uniqueness. Unfortunately, this study indicates that this aim will be difficult to achieve in this manner because strong leakage is observed which is too severe to make meaningful forecasts. Instead, we summarize our attempts to mitigate the leakage while noting that care must be taken when attempting further work in the application of spherical Slepian functions to core surface flows.

In the next section we explain the background and methodology for decomposing flow models using spherical scalar Slepian functions. In “Flow models: synthetic, numerical and inverted” section we describe the three flow model types decomposed—synthetic, numerical and inverted from satellite data—and the results are shown in “Decomposition of flow models” section. Synthetic flow models (maximum spherical harmonic degree and order, $L=60$) are initially decomposed to appraise the technique. We proceed to test the decomposition on a high degree ($L=60$) model from a geodynamo simulation and then a low degree ($L=20$) model from inverting satellite SV data. We discuss our findings in “Analysis of the flow separation” section with consideration of the issues of strong aliasing that we have discovered. We focus on the tangentially ambiguous regions in the paper, but the analysis of the tangent cylinder decomposition is detailed in Additional file 1.

Localizing flow models with Slepian functions

Representing flow models

A compact and convenient way to represent core surface flow is through the use of spherical harmonics, with model coefficients for the toroidal and poloidal scalars describing the flow (e.g., Roberts and Scott 1965). Because the velocity is non-divergent and its radial component across the boundary vanishes, the horizontal velocity vector $\mathbf{u }_H$ can be expressed in terms of the poloidal and toroidal scalars, S and T, which can be expanded in spherical harmonics, in a spherical polar coordinate system $(r, \theta , \phi )$:

$$\begin{aligned} \mathbf{u }_H = \nabla \times (T\mathbf{r }) + \nabla _H (rS) \end{aligned}$$

(1)

where

$$\begin{aligned} T(\theta , \phi ) = \sum _{l,m} t^m_l Y^m_l(\theta , \phi ) \quad \text {and} \quad S(\theta , \phi ) = \sum _{l,m} s^m_l Y^m_l(\theta , \phi ) \end{aligned}$$

(2)

The coefficients $\{t^m_l, s^m_l\}$ are the flow model coefficients, $Y^m_l(\theta , \phi )$ are the Schmidt quasi-normalized real spherical harmonics, and l and m are the degree and order, respectively.

The north–south, $u_\theta$, and west–east, $u_\phi$, components of the horizontal velocity are

$$\begin{aligned} u_{\theta }= & {} \sum _{l=1}^{L} \sum _{m=0}^{l} \left[ - t ^{m(c)}_l \sin (m\phi ) + t ^{m(s)}_l \cos (m\phi )\right] \frac{m}{\sin \theta } P^m_l(\cos \theta ) \nonumber \\&+ \sum _{l=1}^{L} \sum _{m=0}^{l} \left[ s ^{m(c)}_l \cos (m\phi ) + s ^{m(s)}_l \sin (m\phi )\right] \frac{dP^m_l(\cos \theta )}{d\theta }\nonumber \\ u_{\phi }= & {} \sum _{l=1}^{L} \sum _{m=0}^{l} \left[ - t ^{m(c)}_l \cos (m\phi ) - t ^{m(s)}_l \sin (m\phi )\right] \frac{dP^m_l(\cos \theta )}{d\theta } \nonumber \\&+ \sum _{l=1}^{L} \sum _{m=0}^{l} \left[ - s ^{m(c)}_l \sin (m\phi ) + s ^{m(s)}_l \cos (m\phi )\right] \frac{m}{\sin \theta } P^m_l(\cos \theta ) \end{aligned}$$

(3)

where $P^m_l$ are the associated Legendre polynomials for degree (l) and order (m) and superscripts (c) and (s) denote the coefficients of T and S multiplying $\cos (m\phi )$ and $\sin (m\phi )$, respectively.

Decomposing flow models

Spherical harmonics are global functions, but they can be converted by linear transformation into spherical Slepian basis functions that are localized onto two or more regions of the sphere. We summarize the methodology here for clarity but refer the reader to Simons (2010) for a detailed derivation.

Spherical harmonics up to degree and order L are typically expressed as a vector of $(L+1)^2$ elements, each of which is a function of position $(\theta ,\phi )$ on the unit sphere:

$$\begin{aligned} \mathbf{y }({\theta ,\phi })=\left[ \begin{array}{ccccc} Y_{0}^{0}({\theta ,\phi })&\cdots&Y_{l}^{m}({\theta ,\phi })&\cdots&Y_{L}^{L}({\theta ,\phi }) \end{array}\right] ^{{\mathsf {T}}}. \end{aligned}$$

(4)

where the ordering of the spherical harmonics $Y_{l}^{m}$ is arbitrary. The coefficient multiplying the monopole harmonic ($Y_{0}^{0}$) is usually ignored (or set to zero) in geomagnetic and core flow studies, but is included here to preserve generality.

On a unit sphere, a potential $V({\theta , \phi })$ up to degree L is represented in a spherical harmonic basis by a single $(L+1)^{2}$-dimensional vector of coefficients, $\mathbf{v }$. The potential on the surface is obtained from these coefficients as:

$$\begin{aligned} V(\theta ,\phi )=\mathbf{v }\cdot \mathbf{y }(\theta ,\phi ). \end{aligned}$$

(5)

Spherical Slepian functions provide a different set of orthonormal basis functions written as:

$$\begin{aligned} \mathbf{g }_{\alpha }(\theta ,\phi )=\left[ \begin{array}{ccccc} g_{1}(\theta ,\phi )&\cdots&g_{\alpha }(\theta ,\phi )&\cdots&g_{(L+1)^{2}}(\theta ,\phi ) \end{array}\right] ^{{\mathsf {T}}}. \end{aligned}$$

(6)

Each of these basis functions is linearly related to spherical harmonics by the expansion

$$\begin{aligned} g_{\alpha }(\theta ,\phi ) = \mathbf{g }_{\alpha }\cdot \mathbf{y }(\theta ,\phi ). \end{aligned}$$

(7)

$\mathbf{g }$ is produced from the spherical surface harmonic basis by multiplying $\mathbf{y }({\theta , \phi })$ by a unitary matrix

$$\begin{aligned} \mathbf{G }^{{\mathsf {T}}}=\left[ \begin{array}{c} \mathbf{g }_{1}^{{\mathsf {T}}}\\ \vdots \\ \mathbf{g }_{(L+1)^{2}}^{{\mathsf {T}}} \end{array}\right] . \end{aligned}$$

(8)

The matrix $\mathbf{G }$ is constructed by optimization to localize the solution over specified regions (and their complements) for a given maximum spherical harmonic degree L. The procedure determines a complete set of basis functions, which are ordered in terms of contribution to the region considered and, finally, split into the ‘in’ and ‘out’ of region sections as two distinct domains:

$$\begin{aligned} \mathbf{G }^{{\mathsf {T}}}\mathbf{y }(\theta ,\phi )=\left[ \begin{array}{c} \mathbf{G }_{{in}}^{{\mathsf {T}}}\mathbf{y }(\theta ,\phi )\\ \mathbf{G }_{{out}}^{{\mathsf {T}}}\mathbf{y }(\theta ,\phi ) \end{array}\right] =\left[ \begin{array}{c} g_{1}(\theta ,\phi )\\ \vdots \\ g_{K}(\theta ,\phi )\\ g_{K+1}(\theta ,\phi )\\ \vdots \\ g_{(L+1)^{2}}(\theta ,\phi ) \end{array}\right] =\left[ \begin{array}{c} \mathbf{g }_{\text {in}}(\theta ,\phi )\\ \mathbf{g }_{\text {out}}(\theta ,\phi ) \end{array}\right] , \end{aligned}$$

(9)

where K indicates the last element of the functions primarily concentrated in the first domain (in this case ‘in’), and $K+1$ labels the beginning function for the second (‘out’) domain. If an optimal decomposition has occurred, the result of summing the basis functions for the ‘in’ and ‘out’ regions will be identical to the input function and the ‘in’ region will be fully recreated by using all functions up to the K value with no signal outside the region of interest. An optimal decomposition is less likely when a band-limited signal is decomposed because the bandwidth, L, controls spatial resolution. Therefore, the resulting decomposition of a model with restricted spherical harmonic degree is more likely to exhibit spatial aliasing and unwanted leakage outside the region of interest.

The Slepian functions span a linear subspace of $\mathbf{y } (\theta , \phi )$ in which the sum-squared function value–in this case, the energy—over the chosen region, R, is maximized. This ‘localization’ matrix is symmetric, and the subspace of maximum energy is readily obtained by eigenvalue decomposition. We compute the Gram matrix of energy in R as:

$$\begin{aligned} \mathbf{D } = \int _R \mathbf{y }(\theta , \phi ) \mathbf{y }^T (\theta ,\phi ) {\text {d}}\Omega = \int _R \left[ \begin{array}{ccc} Y_0^0Y_0^0 &{} \cdots &{} Y_0^0Y_L^L\\ \vdots &{} \ddots &{} \vdots \\ Y_0^0Y_L^L &{} \cdots &{} Y_L^LY_L^L \end{array}\right] {\text {d}}\Omega \end{aligned}$$

(10)

where the eigenvalues and eigenvectors of $\mathbf{D }$ are defined as:

$$\begin{aligned} \mathbf{DG } = \mathbf{G }{{\varvec{\Lambda }}} \end{aligned}$$

(11)

Each column of $\mathbf{G }$ contains one eigenvector, and ${{\varvec{\Lambda }}}$ is a diagonal matrix with the corresponding eigenvalues:

$$\begin{aligned} {{\varvec{\Lambda }}} = {\text {diag}}(\lambda _1, \ldots , \lambda _\alpha , \ldots , \lambda _{(L+1)^2}). \end{aligned}$$

(12)

Due to the symmetry of $\mathbf{D }$, all of its eigenvalues are positive (or zero) and real, and its eigenvectors are orthogonal, which makes $\mathbf{G }$ unitary. K, called the Shannon number, is typically chosen from inspection of the eigenvalues which contribute less than 50% to the ‘in’ patch (Simons et al. 2006). Note that in this study we are localizing the potentials of the flow, T and S of Eq. (1), not the flows inverted from them. The implications of this will be discussed in “Analysis of the flow separation” section.

Localization regions

To apply spherical Slepian decomposition, a region of interest R must be defined in order to separate the signal into two complementary parts. In this study, the tangentially geostrophic unambiguous regions are chosen as an area of particular interest as well as the two polar caps subtended by the intersection of the cylinder tangent to the inner core.

Tangentially geostrophic ambiguous region

The geostrophic flow assumption is one of the most frequently used within flow modeling to reduce the inherent ambiguity. Le Mouël (1984) and Hills (1979) originally proposed the constraint, which assumes that the pressure gradient, the Coriolis and the buoyancy forces balance in the Navier–Stokes (momentum) equation. Hence, by assuming that gravity is purely radial, the radial component of the thermal wind equation vanishes, giving

$$\begin{aligned} \nabla _{H} \cdot \ (\mathbf{u }_H \cos \theta ) = 0. \end{aligned}$$

(13)

Substituting this into the radial component of the frozen-flux induction equation gives

$$\begin{aligned} \dot{B}_r + \cos \theta \mathbf{u }_H \cdot \nabla _{H} (B_r / \cos \theta ) = 0 \end{aligned}$$

(14)

where $\dot{B}_r$ is the first time derivative of the radial magnetic field (e.g., Holme 2015). The flow is unique at all points along a contour of $B_r / \cos \theta$ that intersects the equator, but the flow elsewhere in the ambiguous patches is only determined in the direction perpendicular to the $B_r/\cos \theta$ contours (Backus and Le Mouël 1986).

Assuming tangential geostrophy, the flow ambiguity disappears over the majority of the outer core surface except in ‘ambiguous patches,’ where only one component is determined. Figure 1 shows the three ambiguous patches, from contours of $B_r/ \cos \theta$ that do not connect to the equator either directly or by saddle points. The shaded regions are the ambiguous patches and correspond to $40\%$ of the surface area of the CMB, similar to values quoted elsewhere in the literature (e.g., Amit and Pais 2013).

Cylinder tangent to the inner core

It is believed that the flow within the polar regions of the outer core inside a cylinder tangential to the inner core (and parallel to the rotation axis, as shown in Fig. 2) is substantially different to the flow outside this zone (Aubert et al. 2013; Hollerbach and Gubbins 2007; Jones 2015; Pais and Jault 2008). One theory to explain the difference in convection within and outside the tangent cylinder is that inside the tangent cylinder the gravity and the rotation vectors are largely parallel, whereas outside they are largely perpendicular (Hollerbach and Gubbins 2007).

The tangent cylinder is closely linked to rapid changes of the magnetic field, within and along the boundaries of the spherical caps, with recent research suggesting fast-moving features can be observed with high resolution field modeling from satellite data (Finlay et al. 2016; Livermore et al. 2016). The Taylor–Proudman theorem states that rotating fluids perturbed by a solid body tend to form non-axisymmetric, z-invariant fluid columns parallel to the axis of rotation, called Taylor columns. These are thought to operate in the Earth’s outer core (Glatzmaier and Roberts 1996). Aurnou et al. (2015) summarized how advanced asymptotically reduced theoretical models, efficient Cartesian direct numerical simulations and laboratory experiments show good agreement that axially coherent, helical convection columns are present and break into three-dimensional geostrophic turbulence. These dynamics are not present within the tangent cylinder. Instead, observations of the Earth’s magnetic field suggest that there are anticyclonic, axisymmetric, z-variant polar vortices inside the tangent cylinder (Hulot et al. 2002; Olson and Aurnou 1999; Sreenivasan and Jones 2005). Cao et al. (2018) recently suggested three alternative mechanisms for the variations in the cylinder tangent to the inner core based on inertia-free, axisymmetric numerical simulations but concluded further work was required to conduct quantitative assessment under Earth’s core conditions.

The area on the outer core surface which lies within the tangent cylinder is described by two spherical caps, each subtending a half-angle of $21^{\circ }$ with respect to the Earth’s rotation axis at the north and south poles. This investigation provides a contrast with other polar caps studies using spherical Slepian functions (Plattner and Simons 2015; Simons and Dahlen 2006). For brevity, the tangent cylinder results are included in Additional file 1.

Flow models: synthetic, numerical and inverted

The flow models on which the investigation of Slepian decomposition are focused come from three sources: synthetic flow models up to degree and order 60, a flow model extracted from a numerical dynamo simulation of Aubert et al. (2013), again to degree and order 60, and a steady flow model inverted from three years of satellite magnetic data to degree and order 20. With these, we seek to illustrate the capability and limitations of the scalar Slepian decomposition of core flow models, which will inform subsequent inversions of SV data for CMB flow in localized regions.

Synthetic flow models

Three sets of synthetic flow models with different kinetic energy spectral ‘colors’ (approximately red, white and blue) were generated to test the robustness of the decomposition methodology. Flow coefficients were randomly sampled from a normal distribution, ${\mathcal {N}}(0,1)$. In order to produce a ‘flat’ (white) spectrum, each coefficient was divided by its degree, l. In the decreasing (blue) power spectrum model, the coefficients were divided by $l \sqrt{l+2}$. The coefficients of the model with increasing (red) spectral energy required no further modification.

The kinetic energy, $E_{\mathrm{Kinetic}}$, for each degree is:

$$\begin{aligned} E_{\mathrm{Kinetic}}(l) = \frac{1}{d^2} \frac{l(l+1)}{2l+1} \sum _{m=0}^{l} \left\{ (t_l^m)^2 + (s_l^m)^2 \right\} \end{aligned}$$

(15)

where d is the radius of the surface considered (here, the CMB) and $t_l^m$ and $s_l^m$ are spherical harmonic coefficients representing the stream function and/or the velocity potential, i.e., corresponding to the toroidal and/or poloidal component (Le Mouël et al. 1985). Therefore, the kinetic energy for each model can be represented on a ‘power spectrum’ or ‘degree variance’ given by Eq. 15, which removes phase information and sums over all orders to give a value at each individual degree (Holme 2015).

Flow from a numerical dynamo simulation

To investigate the decomposition of a more realistic flow model ($L=60$), a snapshot (single time step) from the Coupled Earth numerical geodynamo model was used (J. Aubert, pers. comm., 2017). The model’s spatial and temporal variations are comparable with observations of real fields and contain some Earth-like aspects of the SV (Aubert et al. 2013; Christensen et al. 2010). Further details of the geodynamo simulation and its behavior can be found in Aubert et al. (2013).

Flow inverted from satellite data

Large satellite vector data sets of the Earth’s magnetic field have become available in recent years which give better spatial coverage compared to ground observatory datasets. By combining the data into a set of ‘virtual observatories’ (VO) in space, they can mimic ground-based observatories (Barrois et al. 2018; Beggan et al. 2009; Mandea and Olsen 2006). Improved VO calculations used sums and differences of CHAMP along-track measurements to calculate time series at 500 equally spaced VO, based on quiet data selected from all local times, corrected for external, induced and crustal fields (Hammer 2018). Data from within 200 km of the VO were reduced to the central point using a cubic expansion of the potential. VO SV data and associated $3\times 3$ data covariance matrix for each VO were inverted for a weakly regularized flow to $L=20$ with minimal month-to-month time variation using the algorithm described by Whaler et al. (2016). The average flow from thirty-six months between 2007.0 and 2010.0 has energy roughly equally partitioned between the toroidal and poloidal components, and between the geostrophic and ageostrophic components, in contrast to the more typical predominantly toroidal and geostrophic flows obtained with stronger regularization.

Decomposition of flow models

We present results from the decomposition of each of the flow models for the region representing the tangentially geostrophic unambiguous area of Fig. 1, based on synthetic data, a numerical dynamo and satellite data. The synthetic models are used as a test of the method and give insight into how the shape (i.e., color) of the energy spectrum affects the decomposition. The numerical dynamo flow model allows us to investigate the decomposition of a more realistic flow to a relatively high degree. Although the Shannon number provides the optimal value for the spectral-spatial separation, we also consider varying the number of functions, K, which we use to represent the flow inside the patches and in its complement. We show both the $\phi$ (or west-east) and $\theta$ (or north–south) components.

Synthetic flow models

Figure 3 shows the decomposition of the three synthetic models. The Shannon number of $K = 1504$, out of 3721 coefficients, was used. When the spatial values of the two regions (second and third rows) are added together, the sum (fourth row) creates a map almost identical to the input (first row). The absolute maximum difference between the input and the summed result of the decomposition is $<5 \times 10^{-13}$ km/year, which can be ascribed to rounding error. The bottom row of Fig. 3 shows the power spectra of the input (red), inside (green) and outside (blue) regions. As they match exactly, the summed spectra (black dashed) overlap the input spectra (red).

Although the Slepian decomposition offers an optimal trade-off between spatial and spectral fidelity, there are a number of unavoidable side effects. The plots of the ‘inside’ and ‘outside’ ambiguous patches show two such features. Firstly, there is strong aliasing of the signal along the boundary between the regions. Secondly, the strength of the aliasing depends on the slope of the input spectrum. As can be seen from Fig. 3, more leakage occurs when the spectral energy decreases with degree (right hand column), the typical spectral slope of actual flows.

Numerical dynamo flow model

The numerical dynamo flow decomposition is shown in Figs. 4 and 5. The comparison between the input and sum of the decomposed regions shows very few differences in the spectral and spatial domains. However, when the coefficients inside and outside the ambiguous patches are plotted, strong flows appear along the boundary, particularly in the toroidal $\phi$ (eastward) and the poloidal $\theta$ (southward) components.

It has been noted in previous studies that the Shannon number, dependent on the proportion of the surface area in the patches and the spherical harmonic degree of the model, is not always the best parameter for separating the model into the ‘inside’ and ‘outside’ patches (Beggan et al. 2013; Plattner and Simons 2013). The effect of altering the K value on the decomposition of the dynamo simulation flow is shown in Fig. 6. Reducing K to less than the Shannon number (in this case, $K=1504$), the magnitude of the flow within the region is reduced along with some of the aliasing. It is only by using small values of K that the aliased signals are strongly reduced but, at these low K values, flow within the patches is a poor recreation of the input. This result suggests that although leakage within the inside and outside regions can be ameliorated, artificially large flow values along the region boundaries remain.

Core flow inverted from satellite data

The maximum degree of the flow model inverted from satellite SV data is lower than the synthetic or geodynamo models at $L=20$. As a result, the flow features are larger scale compared to the numerical dynamo or synthetic data. For consistency only the tangentially geostrophic component of the flow was used in the decomposition.

The optimum solution to minimize the boundary leakage was found to be when the K value was 199 out of 441 coefficients, compared to a Shannon number of 178 (Figs. 7, 8). The summed coefficients again produce an exact recreation of the input coefficients, regardless of K value. This is consistent with the experiments using the numerical dynamo model, even though the bandwidth and spectral content of the two are rather different. The separation of the flow into the two parts is not ideal, as energy from the low degrees leaks across the patch boundaries, causing strong flows along the edges.

Analysis of the flow separation

The results from synthetic flow models indicate that all models, irrespective of the slope of their energy spectrum, can be readily divided into multiple regions such that the difference between the addition of the flows from the separated regions compared to the input signal is negligible. However, while the method enables near-perfect reconstruction of the input model from the decomposed regions, there are obvious issues with aliasing/spectral leakage within the individual regions as seen in Figs. 3 to 8. The strong flows seen at the boundaries of the ambiguous patches are particularly apparent in the azimuthal component of the toroidal part of the flow reconstructed both inside and outside it.

To probe the source of the leakage, we examined how individual coefficients are split to effect the flow separation. Some input coefficients are split into a larger magnitude value for inside the patches and an oppositely signed coefficient (of similar magnitude) for outside the patches (as shown in Fig. 9). This is especially noticeable for small coefficients. The majority of the coefficients affected by this splitting are in the higher degrees (though they occur throughout the model). We tested to see whether the artificially enlarged coefficients were the cause of the leakage by plotting flow maps using only those coefficients which had greater absolute values compared to their input coefficient (Fig. 10). The maps have a very distinctive concentration of rapid flow at the region edges and are much weaker away from these boundaries. Thus, the large values produce unrealistic but oppositely directed flows along the patch edges which cancel when the flows from inside and outside patches are summed.

We next checked whether removing these coefficients produced a better decomposition of the flow models. However, the clear distinction between the inside and outside patches disappeared, and the flow within the region of interest was no longer similar to the input. Thus, all coefficients are required for an accurate decomposition of the flow models, and we cannot remove the spectral leakage by ignoring the ‘poorly behaved’ subset of coefficients.

We also investigated the cross-spectral leakage plots by examining the $\mathbf{G } \mathbf{G }^T$ matrix, which should ideally be the identity matrix. The elements of $\mathbf{G } \mathbf{G }^T$ are related to size of the region of interest, its shape, the degree resolution of the model and the truncation level of the basis. The complex shapes of the tangentially geostrophic ambiguous patches produce departures from the identity matrix, but they are no worse than those found by Beggan et al. (2013) for their crustal magnetic decomposition.

Finally, we reran the decomposition of the flow models using the ocean-continent regions for the crustal magnetic field taken from Beggan et al. (2013) and found the same aliasing along the boundaries. This suggests the cause of the leakage is a combination of the energy distribution of the flow model (predominantly red spectra for realistic models) and the use of scalar Slepian decomposition on a vector field given by spatial derivatives of potential functions expanded in spherical harmonics.

Our method was also applied to the regions inside and outside the tangent cylinder (shown in Additional file 1), represented by two caps of $21^\circ$ half-angle (Fig. 2). As with the tangentially geostrophic ambiguous patches decomposition, large aliasing can be observed along the boundaries. As the tangent cylinder region is much smaller than the tangentially geostrophic ambiguous patches, the leakage obscures most of the flow structure within the region, making investigation of these separated flows difficult.

Discussion

The results in “Decomposition of flow models” section allow some conclusions to be drawn about the influence of the tangentially geostrophic unambiguous patch on global flow models. Our experiments suggest that the spatial leakage seen in the flow maps is exaggerated by taking the spatial derivative $\frac{{\text {d}}P^m_l(\cos \theta )}{{\text {d}}\theta }$ of the associated Legendre polynomials when considering the west–east component of the toroidal flow and the north–south component of the poloidal flow (see Eq. (2)). When plotted as vector flow maps, the leakage can be seen but is less obvious, as shown in Fig. 11.

The energy contribution of the tangentially geostrophic unambiguous patch in the satellite flow model is 41% of the toroidal component and 67% of the poloidal component input energy, which is different from the percentage surface area of the patch on the whole CMB (60%). However, the leakage at the boundaries is likely to be affecting the true distribution of kinetic energy. Hence, without further work to minimize the leakage, there remains a large uncertainty. We note that the decomposition tends to be better behaved in the Pacific Ocean, due to the lower flow complexity and magnitude in this region.

The aliasing strongly affects the analysis of the tangent cylinder flow, which Livermore et al. (2016) indicated was of particular interest due to significantly enhanced flow velocities. Spectra indicate that the tangent cylinder contains $\sim 16\%$ of the energy of the summed decomposition flow, which is greater than the $\sim 6\%$ of the energy in the original flow. We also note that without spherical harmonic models to higher degree it is not easy to isolate the detailed core surface flow inside or along the boundary of the tangent cylinder.

This study has highlighted the nature of the aliasing we have found when using Slepian decomposition with vectors represented by scalar potentials. These factors will likely affect direct inversions for flow from SV data into a select region, as undertaken for the crustal magnetic field by Plattner and Simons (2015). This will inform investigations of the different regional dynamics at the CMB, such as in the tangent cylinder or the large low-shear-velocity provinces (LLSVP). LLSVPs are features of low seismic velocity at the base of the mantle with distinct edges (Garnero et al. 2016). Possibilities for this seismic signal include a thermochemical pile or super-plume feature, which could be anomalously hot and/or dense. Should these structures exist, the regional differences on the top of the outer core are likely to influence the flow structures within it, making it an interesting region to study.

Conclusion

Using the tangentially geostrophic ambiguous patches and the tangent cylinder as example regions, we have decomposed global flow models into these spatial patches and their complements using scalar Slepian functions. We analyzed three different types of flow models including a model inverted from satellite magnetic data. However, although the flows in the region of interest and its complement successfully sum to the input, the method produces strong aliasing along the region boundaries, to the extent that flows are partially obscured.

We examined the reasons for this and conclude that further work is required to reduce spectral leakage between the regions. We intend to apply spherical Slepian decomposition to the flow treated as a vector quantity in areas of interest, rather than decomposing the globally representative scalar potential, and to undertake direct inversion of SV data for localized flow.