TY - JOUR
AU - Finlay, Christopher C
AB - SUMMARY We present a new model of time-dependent flow at low latitudes in the Earth’s core between 2000 and 2018 derived from magnetic field measurements made on board the Swarm and CHAMP satellites and at ground magnetic observatories. The model, called CoreFlo-LL.1, consists of a steady background flow without imposed symmetry plus a time-dependent flow that is dominated by geostrophic and quasi-geostrophic components but also allows weak departures from equatorial symmetry. Core flow mode amplitudes are determined at 4-month intervals by a robust least-squares fit to ground and satellite secular variation data. The l1 norm of the square root of geostrophic and inertial mode enstrophies, and the l2 norm of the flow acceleration, are minimized during the inversion procedure. We find that the equatorial region beneath the core–mantle boundary is a place of vigorous, localized, fluid motions; time-dependent flow focused at low latitudes close to the core surface is able to reproduce rapid field variations observed at non-polar latitudes at and above Earth’s surface. Magnetic field acceleration pulses are produced by alternating bursts of non-zonal azimuthal flow acceleration in this region. Such bursts are prominent in the longitudinal sectors from 80–130°E and 60–100°W throughout the period studied, but are also evident under the equatorial Pacific from 130°E to 150°W after 2012. We find a distinctive interannual alternation in the sign of the non-zonal azimuthal flow acceleration at some locations involving a rapid crossover between flow acceleration convergence and divergence. Such acceleration sign changes can occur within a year or less and, when the structures involved are of large spatial extent, they can give rise to geomagnetic jerks at the Earth’s surface. For example, in 2014, we find a change in the sign of the non-zonal azimuthal flow acceleration under the equatorial Pacific as a region of flow acceleration divergence near 130°E changes to a region of flow acceleration convergence. This occurs at a maximum in the amplitude of the time-varying azimuthal flow under the equatorial Pacific and corresponds to a geomagnetic jerk at the Earth’s surface. Satellite magnetics, Inverse Theory, Rapid time variations, Core 1 INTRODUCTION Time variations in the motion of liquid metal in Earth’s outer core produce fluctuations in the geomagnetic field. A striking example of this process is the occurrence of a sharp change in the field acceleration known as a geomagnetic jerk (Courtillot & Le Mouël 1984; Mandea et al.2010; Brown et al.2013). Improved satellite-data-based models in the late 2000s made it possible for the first time to reliably estimate the global secular acceleration (SA) and its time variations since 2000 (Olsen & Mandea 2008; Chulliat et al.2010; Lesur et al.2010b). Global maps of the field acceleration at the core surface revealed the existence of distinctive sub-decadal SA pulses, their prominence at low latitudes, their alternating sign, and their relationship with jerks observed at the Earth’s surface which can be seen as by-products of successive pulses (Chulliat et al.2010; Chulliat & Maus 2014; Torta et al.2015). An important step towards understanding the underlying physical mechanism producing the phenomenon of SA pulses is to determine how the core flow and its patterns of acceleration have altered during these events. Here, in an effort to clarify this process, we construct a new model of the time-varying core flow spanning the past 18 yr, specifically designed to focus on core flow time variations at low latitudes. The core flow is linked to geomagnetic field changes through the magnetic induction equation (Roberts & Scott 1965; Bloxham & Jackson 1991; Holme 2015). Under the approximation that magnetic diffusion plays a secondary role, core flow advects and stretches the magnetic field producing its observed changes. This enables field variations to be inverted for the responsible core flow (Kahle et al.1967) but unfortunately this procedure is highly non-unique since flow around contours of Br produces no observable field change (Backus 1968). The difficulties are exacerbated by the fact that observations constrain only the largest lengthscales of the field (Rau et al.2000; Roberts & Glatzmaier 2000; Amit et al.2007). In order to make progress, it is therefore necessary to make assumptions concerning the underlying nature of the flow (Gubbins 1982; Le Mouël 1984; Pais & Jault 2008; Aubert 2012) and to account for unobserved small-scale processes (Hulot et al.1992; Eymin & Hulot 2005; Gillet et al.2009, 2015a). Most studies of core flow time variations have been based on spherical harmonic field models derived primarily from ground observatory data. In a landmark study, Jackson (1997) constructed time-dependent flows spanning 1840–1990 and demonstrated that time variations were needed to fit observatory data to an acceptable level. Significant oscillatory motions were evident in these flows, especially at low latitudes (see his fig. 3 and related discussion). These equatorial flow oscillations have attracted surprisingly little attention in the intervening 20  years. Le Huy et al. (2000) investigated changes in flow accelerations at the time of geomagnetic jerks and found similar acceleration patterns, but with alternating signs, could explain the 1969, 1979 and 1992 jerks. Bloxham et al. (2002) argued that time variations in only the equatorially symmetric, zonal, toroidal flow components were sufficient to account for geomagnetic jerks. More recently, using observatory monthly mean data processed to isolate the core field signal, Whaler et al. (2016) showed that such simple time-dependent flows are unable to adequately fit the most recent observations, finding that temporal variability of the non-zonal azimuthal flow was required at low latitudes. A major observational advance in the past two decades has been the ability to monitor geomagnetic secular variation (SV) and SA from space. This has enabled studies of core flow time variations based on global data constraints. Building on initial studies of the high-resolution flow at a single epoch (e.g. Holme & Olsen 2006), Olsen & Mandea (2008) were the first to find evidence for sub-decadal bursts of flow acceleration. In 2003 they inferred a strong equatorward flow acceleration under Asia accompanied by westward acceleration at low latitudes under the Indian ocean. Also considering flow changes around 2003, Wardinski et al. (2008) found that time variations of both the equatorially symmetric, zonal, toroidal flow coefficients, and the non-zonal toroidal flows with azimuthal wavenumber m = 1 and 4, were required. Similarly, Silva & Hulot (2012) concluded that while changes in flow acceleration around 2003 could be purely equatorially symmetric, they could not be only zonal and toroidal, inferring a change in the sign of non-zonal flow accelerations under India and Indonesia. Lesur and colleagues (Lesur et al.2010a; Wardinski & Lesur 2012; Lesur et al.2015) have pursued an approach that involves the co-estimation of field and flow models from satellite data while Baerenzung et al. (2014, 2016, 2017) have developed probabilistic flow inversion methods. Taken together, these studies suggest that in the absence of strong physical constraints, primarily toroidal flows, with only a small amount of poloidal flow and a comparatively weak time-dependence, are adequate to explain the observations; this point was also made early on by Whaler (1986). The use of appropriate prior information is essential in mitigating flow non-uniqueness and in making the most of the information provided by geomagnetic observations. The most important factor controlling the form of fluid motions in Earth’s core is its rapid rotation; this results in a leading order geostrophic force balance between Coriolis forces and pressure gradients (e.g. Aubert et al.2017). Hide (1966) pointed out that slow flow perturbations about a state of geostrophic balance are expected to take a columnar form varying only slowly parallel to the rotation axis. Such motions are commonly referred to as quasi-geostrophic (QG). Busse (1970) showed that convection in a rotating spherical shell takes this form. Jault (2008) showed that transient flow disturbances are expected to take a QG form even in the presence of strong magnetic fields provided the inertial wave speed (associated with rotation) is much faster than the Alfvén wave speed. Further tests using 3-D simulations, with strong imposed magnetic fields, were performed by Gillet et al. (2012) who found that QG motions indeed dominate on large spatial scales and short timescales. Schaeffer & Pais (2011) further noted that although equatorially symmetric QG flows were dominant with equatorially symmetric forcing, in the presence of equatorially antisymmetric forcing, cross-equatorial flow could occur even in the rapidly rotating regime. In the past few years, it has become possible to carry out 3-D dynamo simulations that are beginning to approach the correct rapidly rotating regime; Schaeffer et al. (2017) found large-scale fluid motions outside the tangent cylinder remained primarily columnar, while Aubert (2018) showed that not only the flow, but also the flow acceleration remained primarily columnar in the rapidly rotating regime, provided the Alfvén wave timescale governing magnetic disturbances was much longer than the timescale of rotation. Bardsley & Davidson (2016) have also clarified how columnar motions can be built by propagating energy rapidly along the rotation axis at the group velocity of inertial Alfvén waves in the case that the disturbance wavevector is perpendicular to the rotation axis. Aubert (2018) has highlighted that once formed, such QG disturbances can result in non-zonal Alfvén waves due to the heterogeneous nature of a dynamo-generated magnetic field. Taken together, these studies point to a crucial role for QG motions in time-varying core flow. Many studies have attempted to incorporate constraints based on rapid rotation into the determination of core flow. Le Mouël (1984) suggested constraining motions to be tangentially geostrophic at the core surface such that horizontal pressure gradients and the Coriolis force are in balance. This approach has since been exploited by many workers (e.g. Jackson 1997; Holme & Olsen 2006). Jault et al. (1988) and Jackson et al. (1993) went further and extended the zonal, equatorially symmetric, part of such flows throughout the depth of the spherical shell showing that the resulting motions were in good agreement with observed changes in the length-of-day. Pais & Jault (2008) pointed out that QG flows should in addition be equatorially symmetric. Enforcing this condition enabled them to highlight the presence of an eccentric and planetary-scale anticyclonic gyre consisting of westward flow at mid-to-low latitudes in the Atlantic hemisphere and at high latitudes under the Pacific hemisphere that has now been broadly supported by many subsequent studies. Gillet et al. (2009) derived time-dependent QG flows for the periods 1960–2002 and 1997–2008 using an ensemble method to take into account the impact of unresolved scales. They found the time-varying part of their flow, which involved prominent interannual variations, was dominated by non-zonal flow structures. Gillet et al. (2012) constructed QG and tangentially geostrophic flows spanning 1840–2008 and found that their tangentially geostrophic flows become increasingly equatorially symmetric (i.e. increasingly QG) as the quality of the observational data improved, concluding that QG flows were able to account well for the best quality recent data. Pais et al. (2014) have studied the mean flow and major empirical modes of flow temporal variability on decade to century timescales within a QG framework. On top of the planetary-scale eccentric gyre, they found that a small number of empirical modes involving large-scale vortices and zonal motions coupled between low and high latitudes can account for 90 per cent of the variance in the observed SV and for decadal changes in the length-of-day. Gillet et al. (2015a) constructed QG flows from the COV-OBS field model, accounting for the first time for time-correlated errors due to unresolved scales, finding non-zonal, time-dependent structures in the azimuthal flow within 10° of the equator. Describing these features as non-zonal equatorial jets, and finding that they possessed significant power in periods between 4 and 9.5 yr, it was argued that these may be responsible for the low-latitude SA pulses. Finlay et al. (2016) obtained similar results using the same method to invert the higher resolution CHAOS-6 geomagnetic field model spanning 1999–2016. Schaeffer & Pais (2011) studied the effects of allowing more power in zonal flows (i.e. anisotropy) and departures from equatorial symmetry within the framework of predominantly QG flows. They found a preference for small-scale flow to cross the equator under Indonesia but concluded that the time-dependent large-scale core flow was probably dominated by QG motions on decadal timescales. Amit & Pais (2013) have investigated differences between tangentially geostrophic and QG (columnar) flows and in all cases found intense patterns of upwelling and downwelling at low latitudes, particularly under the Indian ocean. Barrois et al. (2017, 2018a) recently used a stochastic ensemble Kalman filter algorithm to derive core flows consistent with prior information from a numerical dynamo simulation; the prior flow in this case was predominantly equatorially symmetric but not strictly QG. The resulting flows were again dominated by the planetary-scale eccentric gyre and involved prominent upwelling and downwelling structures in the equatorial region. In a follow up studies Barrois et al. (2018a,b), applying the same technique to ground and satellite SV data from the past 20 yr, found low-latitude flow accelerations on interannual timescales, notably under the Pacific. Despite much progress, an observation-based study of rapid core flow dynamics at low latitudes, taking into account the effects of rapid rotation and using the latest ground and satellite data, with a focus on clarifying the mechanism of SA pulses and geomagnetic jerks, has been lacking. This is the aim of the study. We adopt an alternative approach to describing flows expected in a rapidly rotating spherical geometry. Solutions to the Navier–Stokes equation in an incompressible fluid sphere, subject to a non-penetration boundary condition, can in general be decomposed into a geostrophic mode plus inertial modes (Zhang & Liao 2017). Analytical solutions for these are available in a full sphere geometry (Zhang et al.2001; Zhang & Liao 2017). A special class of slow, equatorially symmetric, inertial modes, which we refer to as QG modes (Zhang et al.2001; Busse et al.2005; Maffei et al.2017), have been shown to efficiently describe rotating flow in a sphere at the onset of convection (Zhang & Liao 2004; Zhang et al.2007), and when combined with the geostrophic mode can also describe weakly nonlinear convection (Zhang & Liao 2004; Liao et al.2012). More generally, geostrophic and inertial modes can be used to describe the transient response of rotating spherical systems to a forcing (Liao & Zhang 2010) and in principle they provide a complete basis for representing flows in a rotating sphere (Cui et al.2014; Ivers et al.2015; Backus & Rieutord 2017). Of particular interest here is that these inertial modes are well-suited to describing motions at low latitudes. They do not suffer from a formal requirement for the boundary slope to be small as is the case in classical QG theory (e.g. Busse 1970; Pais & Jault 2008; Schaeffer & Cardin 2006), and they are able to efficiently describe equatorially trapped inertial wave motions (Zhang 1993). We use a combination of inertial and geostrophic modes to parametrize the core flow. In particular, our time-dependent flow consists primarily of QG (i.e. almost axially invariant, equatorially symmetric) inertial modes and a geostrophic mode. Our hypothesis is that such modes provide a concise description of the spatial structure of core flow. We therefore minimize the l1 norm of the square root of the mode enstrophies; this favours a sparse distribution of amplitude amongst the modes, while also penalizing small scales to ensure spectral convergence. In addition, we seek solutions for which the flow acceleration is minimized. Further details of our parametrization of the flow and the model estimation scheme are given in Section 3. We use the frozen-flux induction equation, together with a potential field formalism, to link the flow to SV time-series collected at ground observatories and at satellite altitude, via virtual observatories (Mandea & Olsen 2006; Olsen & Mandea 2007). Details of the data selection and processing are given in Section 2. Note that this enables us to avoid using spherical harmonic field model based descriptions of the field SV and SA that depend heavily on the chosen temporal regularization or window length for time-averaging, especially above spherical harmonic degree 9. We account for the influence of unresolved scales of the magnetic field via an augmented state approach (Gillet et al.2015a; Barrois et al.2017) where the small scales are assumed to be time-correlated (Gillet et al.2015a) and their impact is calculated iteratively using the induction equation. In Section 5 we present a detailed description of our new reference flow model, CoreFlo-LL.1, including its fit to the data, time-average structure and time-dependence. We find a steady background flow dominated by an eccentric planetary gyre similar to that found in previous studies and a time-varying flow that is strongly localized in the equatorial region and consists of vigorous non-zonal azimuthal flows. This flow explains well-observed SV and SA at mid and low latitudes. In Section 6 we discuss the origin of SA pulses, which we find are due to intense, alternating, bursts of non-zonal azimuthal flow acceleration. We find geomagnetic jerks can occur when there is an abrupt sign change or crossover in the azimuthal flow acceleration of sufficiently large scale. We go on to investigate variants of our flow model derived using different assumptions for the spatial regularization norm and regarding the equatorial symmetry. Connections are made to previous studies and implications for the dynamics of the core discussed, before we finish with concluding remarks in Section 7. 2 DATA 2.1 Ground observatory series Ground magnetic observatories provide high-quality vector observations of the geomagnetic field at a global network of locations on Earth’s surface (e.g. Matzka et al.2010). We make use of time-series of revised monthly means (Olsen et al.2014) based on quality-controlled (Macmillan & Olsen 2013) hourly mean magnetic field observations from the BGS database,1 version 0116, that was built from INTERMAGNET and WDC Edinburgh data as available in 2018 August. We use data spanning the interval 2000–2018. Revised monthly means aim to isolate as far as possible the core field signal and are constructed by removing from the hourly mean values (i) estimates of the large-scale magnetospheric field, taken in this case from the CHAOS-6-x7 field model (Finlay et al.2016); (ii) estimates of the ionospheric Sq field and its associated induced field [taken from the CM4 model (Sabaka et al.2004) driven by the latest F10.7 radio flux series], followed by taking a robust mean of the hourly values from all local times, using Huber weights to account for a long-tailed distribution of data. The resulting revised monthly means were then averaged over 4-month windows (to match the time sampling of the satellite-based virtual observatories described below) and annual differences taken to obtain time-series of SV at each observatory location. In this study, we choose to use only data collected at non-polar latitudes, equatorwards of ±60° geographic latitude, since we wish to focus on low-latitude flow dynamics and to avoid the complexities associated with magnetosphere–ionosphere coupling currents and related ionospheric currents in the polar regions. Our final data set then consists of data from 146 ground observatories; their locations are marked by squares in Fig. 1. Examples of SV time-series from low-latitude locations (MBO - M’Bour in Senegal, West Africa; KOU - Kourou in French Guinea, South America; GUA - Guam, USA, Western Pacific) are shown by the black dots in Fig. 2. Figure 1. Open in new tabDownload slide Maps of data locations showing the median absolute residuals from 2000 to 2018 for the three SV components at the virtual (circles) and ground (squares) observatories. Thick markers indicate locations of data shown in Fig. 2. Figure 1. Open in new tabDownload slide Maps of data locations showing the median absolute residuals from 2000 to 2018 for the three SV components at the virtual (circles) and ground (squares) observatories. Thick markers indicate locations of data shown in Fig. 2. Figure 2. Open in new tabDownload slide SV time-series for the ground observatories in GUA, MBO and KOU, and three virtual observatories, all showing the SV data (black), the CHAOS-6-x7 SV prediction (red) and the SV prediction from the CoreFlo-LL.1 model (green). Note that there is a data gap in the time-series of the virtual observatories due to lack of satellite data during that time. Figure 2. Open in new tabDownload slide SV time-series for the ground observatories in GUA, MBO and KOU, and three virtual observatories, all showing the SV data (black), the CHAOS-6-x7 SV prediction (red) and the SV prediction from the CoreFlo-LL.1 model (green). Note that there is a data gap in the time-series of the virtual observatories due to lack of satellite data during that time. Estimates of the data uncertainty for each component of each observatory were derived as follows. We calculated a three-by-three covariance matrix for each observatory from the revised monthly mean time-series of the three components after removing the predictions of the CHAOS-6-x7 field model and detrending. The uncertainty estimates were taken to be the square root of the diagonal elements of these covariance matrices. A similar procedure has been used to determine error estimates for observatory data in recent revisions of the CHAOS model and in the studies of core flow by Whaler et al. (2016) and Barrois et al. (2018b) that also used ground observatory monthly means. 2.2 Satellite-based virtual observatory series Low-Earth-orbit satellites, with dedicated systems for measuring the vector components of the geomagnetic field with absolute accuracy, provide global information on the geomagnetic field within a few weeks. However, since the position of the satellites are constantly changing, it is difficult to directly assess secular variation from the measurements. Here, we follow Mandea & Olsen (2006) and Olsen & Mandea (2007) and derive ‘virtual observatory’ field estimates at a grid of fixed locations based on measurements made on board the CHAMP and Swarm missions within 4-month time windows. Use of 4-month time windows helps reduce local time biases that are found when 1-month windows are used as in the originally proposed version of virtual observatories. Our virtual observatory time-series were derived using a modified version of the original procedure (for full details see Barrois et al.2018b). Here we give only a brief summary of the key steps. We used magnetic field measurements collected by the CHAMP vector field magnetometer (VFM) between 2000 July and 2010 September and by the Swarm VFMs, on board all three satellites (Alpha, Bravo and Charlie), between 2014 January and 2018 July. Starting from the CHAMP MAG-L3 and Swarm Level 1b MAG-L, version 0505, data products, we sub-sampled the data at 15 s intervals using the data in the VFM frame and rotating them into an Earth-Centered Earth-Fixed coordinate frame using the Euler rotation angles estimated in the CHAOS-6-x7 field model. We then selected data from dark regions during geomagnetically quiet times according to the criteria (i) the Sun was a maximum of 10° above the horizon; (ii) geomagnetic activity index Kp < 3°; (iii) the RC disturbance index (Olsen et al.2014) had |dRC/dt| < 3 nT/h; (iv) merging electric field at the magnetopause, Em ≤ 0.8 mV/m with |$E_m=0.33 v^{4/3} B_t^{2/3} \mathrm{sin}(\vert \Theta \vert /2)$|⁠, v the solar wind speed, Θ = arctan(By/Bz) and |$B_t=\sqrt{B_y^2+B_z^2}$| with By and Bz the components of the interplanetary magnetic field (IMF) calculated using 2 hr means of 1-min values of the IMF and solar wind extracted from the OMNI database2; (v) IMF Bz > 0 nT and IMF |By| < 10 nT, again based on 2 hr mean of 1-min values. Next, estimates of the magnetospheric and its related induced fields, as given by the CHAOS-6-x7 model, the ionospheric and its induced fields, as given by the CM4 model, as well as the static internal field for spherical harmonic degrees n > 20 given by the CHAOS-6-x7 model were removed from data. Robust inversions for a static local potential were then carried out in 4-month windows using all data within 700 km of a specified target point. As a pre-processing step and in order to aid the determination of robust data weights, the time-dependent internal field from CHAOS-6-x7 (for spherical harmonic degrees 1 to 20) was first removed from the data, then a robust least-squares inversion for the parameters of the local potential (up to cubic terms) was carried out. Barrois et al. (2018b) present the relevant expressions for the potential. The three components of the field at the target point were then calculated from the potential and added back on to the prediction of CHAOS-6-x7 at the target location. Note that this does not prevent our 4-monthly VO series from departing from CHAOS-6-x7; each estimate is free to depart from the CHAOS model as much as it is required by the data, and essentially identical results are obtained if IGRF-12 rather than CHAOS-6-x7 is used during this processing step. Rather than working directly with the satellite vector field data, we use sums and differences of the vector data, along track (using 15 s differences) and in the east–west direction (from Swarm Alpha and Charlie constructed by searching for data close in time with smallest differences in latitude). Use of sum and difference data has been found to be advantageous for studying the internal field (Olsen et al.2015; Sabaka et al.2015). In inverting for the potential, we use sums and differences of the matrices linking the vector data at each location to the potential, so no error is made due to the small change in the local coordinate system between the points being differenced. Time-series of such 4-month virtual observatory estimates were constructed on an approximately equal area grid of 300 locations, each separated by approximately 12°, using a recursive zonal equal area partitioning algorithm (Leopardi 2006). The virtual observatories derived from CHAMP data were constructed at 300 km altitude and those derived from Swarm data at 500 km altitude. This means, on average and only considering the increase in altitude, that the magnitude of the SV for Swarm is 7  nT yr−1 smaller than for CHAMP. We again chose to use only data series from latitudes equatorwards of ±60°, so we finally retained data from 258 virtual observatory locations at low and mid latitudes marked using the circles in Fig. 1. Secular variation estimates were computed by taking annual differences of the 4-monthly mean virtual observatory estimates. Examples of virtual observatory time-series at low latitudes are shown in the lower half of Fig. 2. Note the gap in the series between 2010, when the CHAMP mission ended, and 2014 when data was available from Swarm. Data uncertainty estimates for the virtual observatory series were derived by a procedure very similar to that used for the ground observatories. For each location, covariances were calculated between the time-series of the three components (after removing from each series the predictions of the CHAOS-6-x7 model and detrending), in order to obtain a three-by-three covariance matrix. A robust procedure for calculating the covariances was employed and the square root of the diagonal elements of the covariance matrices were taken to be the uncertainty estimates for each series. 3 MODEL FORMULATION 3.1 Flow in a rapidly rotating sphere In this section, we present details of our parametrization of the core flow. Our fundamental assumptions are (i) that rotation dominates core dynamics, and (ii) that the core geometry may be approximated by a sphere. The theory presented below follows closely the notation of Zhang & Liao (2017) where more details of the mode-based approach to studying flow in a rotating fluid can be found. We begin with the Navier–Stokes equations describing fluid motions in an incompressible electrically conducting fluid, uniformly rotating with angular frequency |$\boldsymbol {\Omega }$| (e.g. Gubbins & Roberts 1987) \begin{eqnarray*} \begin{aligned} \frac{\partial ^{}\mathbf {u}}{\partial t^{}} + \mathbf {u}\cdot \nabla \mathbf {u} \mathbf{ +2\boldsymbol {\Omega }}\times \mathbf {u}&=-\frac{1}{\rho _0}\nabla P+ \mathbf {F}\\ \nabla \cdot \mathbf {u}&=0, \end{aligned} \end{eqnarray*} (1) where |$\mathbf {F}= (\nu \nabla ^2\mathbf {u} +\rho ^{\prime }\mathbf {g}+\frac{1}{\rho _0}\mathbf {J}\times \mathbf {B}$|⁠) is a forcing/dissipation term, ρ′ is the relative departure from the reference density ρ0, |$\mathbf {g}$| is the acceleration due to gravity, ν is the coefficient of kinematic viscosity, |$\mathbf {J}$| is the current density and |$\mathbf {B}$| is the magnetic field. P is the departure from the pressure of the reference state that contains both the hydrostatic pressure and the contribution of the centrifugal force. In order to obtain fundamental solutions representing the basic structure of rotation-dominated flows, we set the forcing and dissipation term |$\mathbf {F}=0$|⁠, and further neglect the nonlinear term |$\mathbf {u}\cdot \nabla \mathbf {u}$| by assuming that the flow departs only slightly from solid body rotation. For a homogeneous and inviscid fluid confined to a spherical container of radius r0 that is uniformly rotating with |$\boldsymbol {\Omega }=\Omega \mathbf {\hat{z}}$|⁠, where |$\mathbf {\hat{z}}$| is the unit vector in the z-direction, and using the inverse of the angular frequency Ω−1 and the radius r0 as reference timescales and lengthscales, eq. (1) reduces in non-dimensional form to \begin{eqnarray*} \begin{aligned} \frac{\partial ^{}\mathbf {u}}{\partial t^{}} +2\mathbf {\hat{z}}\times \mathbf {u}&=-\frac{1}{\rho _0}\nabla P\\ \nabla \cdot \mathbf {u}&=0, \end{aligned} \end{eqnarray*} (2) while the required boundary condition is that the flow normal to the boundary (in direction |$\mathbf {\hat{r}}$|⁠) vanishes at the spherical container at r = 1, \begin{eqnarray*} \mathbf {\hat{r}}\cdot \mathbf {u}=0. \end{eqnarray*} (3) Eqs (2) and (3) together define a boundary value problem with solutions consisting of a single geostrophic mode, corresponding to a steady flow and an infinite number of time-dependent inertial modes (Zhang & Liao 2017). The total flow can then be conveniently expressed as a linear combination of these modes. Here, we briefly present the key parameters of these modes and refer the reader to Zhang & Liao (2017) for a complete description. 3.1.1 Geostrophic mode The geostrophic mode is the solution for a steady flow in a rotating incompressible fluid that is moving slowly with respect to rigid rotation and when viscosity can be neglected. It follows from the geostrophic balance equation which is found from eq. (2) by neglecting time-dependence \begin{eqnarray*} 2\mathbf {\hat{z}}\times \mathbf {u}=-\frac{1}{\rho _0}\nabla P. \end{eqnarray*} (4) Applying the curl leads to the Taylor–Proudman theorem which states that \begin{eqnarray*} \frac{\partial ^{}\mathbf {u}(\mathbf {r})}{\partial z^{}}=0. \end{eqnarray*} (5) Hence, the geostrophic solution describes a flow that is invariant along the rotation axis (i.e. the z-axis). In order to satisfy the boundary condition in a spherical container, eq. (3), the geostrophic mode must be axisymmetric and purely azimuthal. The geostrophic mode can, therefore, be written as an infinite sum of the so-called ‘geostrophic’ polynomials G2k − 1 for integers k ≥ 1 (Zhang & Liao 2017) of the form \begin{eqnarray*} {\bf u}^{\rm G}(r,\theta )=\sum_{k\ge 1}a^{\rm G}_k G_{2k-1}(r,\theta ){\hat{\boldsymbol {\Phi}}} \end{eqnarray*} (6) involving unknown coefficients |$a^\mathrm{G}_k$|⁠, the radius r, the colatitude θ and the unit vector |${\hat{\boldsymbol \Phi}}$| in azimuth (see Appendix A1 for details of G2k − 1). Although the geostrophic polynomials depend on the two coordinates r and θ, they take a single argument s ≡ rsin θ, of which they are odd functions. They are also orthogonal and can be normalized with respect to their mean value over the sphere, that is, given a second geostrophic polynomial G2l − 1 for an integer l, they satisfy \begin{eqnarray*} \frac{3}{4\pi }\int _{\mathcal V} G_{2k-1}G_{2l-1}{\rm d}{\mathcal V}=\delta _{kl}, \end{eqnarray*} (7) where δkl is the kronecker delta and |$\int _{\mathcal {V}}\mathrm{d}\mathcal {V}$| denotes the integration \begin{eqnarray*} \int _{\mathcal {V}}\mathrm{d}\mathcal {V}\equiv \int _{0}^{2\pi }\mathrm{d}\phi \int _{0}^{\pi }\sin \theta \mathrm{d}\theta \int _{0}^{1}r^2\mathrm{d}r, \end{eqnarray*} (8) with ϕ being the longitude. 3.1.2 Inertial modes Allowing for time-dependent flows, the full solution to eq. (2) is obtained using the superposition of an infinite number of the so-called inertial modes. Following Zhang & Liao (2017), these modes can be specified in terms of three complex-valued velocity components in spherical polar coordinates |$\mathbf {r}=(r,\theta ,\phi )$| as \begin{eqnarray*} \mathbf {u}(\mathbf {r},t)=[u_r(r,\theta ,\phi ),u_\theta (r,\theta ,\phi ),u_\phi (r,\theta ,\phi )]\mathrm{e}^{\mathrm{i}2\sigma t}, \end{eqnarray*} (9) where i is the imaginary number and |$\sigma =\tfrac{\omega }{2}$| a half-frequency, which is real-valued and bounded by 0 < |σ| < 1. Each inertial mode has a specific spatial structure that is associated with a specific value of σ. The inertial modes can be divided into two classes according to their symmetry with respect to the equatorial plane. The equatorially symmetric (ES) inertial modes satisfy \begin{eqnarray*} \big (u_r,u_\theta ,u_\phi \big )(r,\theta ,\phi )&=\big (u_r,-u_\theta ,u_\phi \big )(r,\pi -\theta ,\phi ), \end{eqnarray*} (10) whereas the equatorially antisymmetric (EA) inertial modes obey \begin{eqnarray*} \big (u_r,u_\theta ,u_\phi \big )(r,\theta ,\phi )&=\big (-u_r,u_\theta ,-u_\phi \big )(r,\pi -\theta ,\phi ). \end{eqnarray*} (11) We follow the notation of Zhang & Liao (2017) and specify the inertial modes by a triple index (m, n, k) with m, k ≥ 0 and n ≥ 1. The azimuthal wavenumber m controls the periodic structure in azimuth and is zero in case of axisymmetric modes. The index n determines the degree of complexity in the z direction, along the rotation axis, whereas k indicates the structure along the cylindrical radius perpendicular to the axis of rotation. The inertial modes, hereafter denoted |$\mathbf {u}_{mnk}$|⁠, are orthogonal when integrated over the full sphere and can therefore be normalized such that \begin{eqnarray*} \frac{3}{4\pi }\int _\mathcal {V}\mathbf {u}_{mnk}\cdot \mathbf {u}_{m^{\prime }n^{\prime }k^{\prime }}^*\mathrm{d}\mathcal {V}=\delta _{mm^{\prime }}\delta _{nn^{\prime }}\delta _{kk^{\prime }}, \end{eqnarray*} (12) where * denotes the complex conjugate. Since the expressions for the inertial modes are sinusoidal in azimuth, both the real and the imaginary parts are needed to correctly account for the phase in azimuth. Hence, there are always two coefficients associated with each inertial mode. We introduce ES and EA modes separately in the following. 3.1.2.1 Equatorially symmetric modes The ES inertial modes |$\mathbf {u}^\mathrm{S}_{mnk}$| can be divided into axisymmetric (m = 0) and non-axisymmetric (m ≥ 1) modes. The half-frequencies σ of the axisymmetric ES modes are the roots of the polynomial (Zhang & Liao 2017) \begin{eqnarray*} 0=\sum _{j=0}^{k-1}(-1)^{j}\frac{[2(2k-j)]!}{j!(2k-j)![2(k-j)-1]!}\sigma ^{2(k-j)}, \end{eqnarray*} (13) with k ≥ 2. Hence, for a given k, there are k − 1 axisymmetric ES modes corresponding to the k − 1 distinct roots. In case of the non-axisymmetric ES modes, the half-frequencies satisfy \begin{eqnarray*} 0=\sum _{j=0}^{k}(-1)^{j}\frac{[2(2k+m-j)]!}{j!(2k+m-j)![2(k-j)]!}\left[(m+2k-2j)-\frac{2(k-j)}{\sigma }\right]\sigma ^{2(k-j)}, \end{eqnarray*} (14) with k ≥ 1. There are 2k modes for a given k and m ≠ 0. The absolute value of the roots can be arranged in ascending order 0 < |σ1| < |σ2| < … < |σ2k| with the first being the longest period ES inertial mode for a given m and k (see Appendix A2 for explicit expressions). 3.1.2.2 Quasi-geostrophic modes A special subset of the ES inertial modes are those characterized by having the longest periods. Their force balance is closest to the geostrophic balance in eq. (4) that results in a flow that is almost invariant with respect to the rotation axis. Following Zhang & Liao (2017), modes of this subset of ES inertial modes will therefore be referred to as being QG. They correspond to the half-frequencies |$\sigma ^\mathrm{S}_{m1k}$| and are modes with the index n = 1 in accordance with the previously applied ordering. 3.1.2.3 Equatorially antisymmetric modes The EA inertial modes are denoted |$\mathbf {u}^\mathrm{A}_{mnk}$| and can be again separated into axisymmetric and non-axisymmetric modes. The half-frequencies of the axisymmetric EA modes satisfy (Zhang & Liao 2017) \begin{eqnarray*} 0=\sum _{j=0}^{k}(-1)^{j}\frac{[2(2k-j+1)]!}{j!(2k-j+1)![2(k-j)]!}\sigma ^{2(k-j)}, \end{eqnarray*} (15) with k ≥ 1, whereas the half-frequencies of the non-axisymmetric EA modes are the roots of \begin{eqnarray*} 0=\sum _{j=0}^{k}\frac{(-1)^{j}[2(2k+m-j+1)]!}{j!(2k+m-j+1)![2(k-j)+1]!} \left[(m+2k-2j+1)-\frac{2(k-j)+1}{\sigma }\right]\sigma ^{2(k-j)+1}, \end{eqnarray*} (16) with k ≥ 0. Hence, there are k axisymmetric and 2k + 1 non-axisymmetric EA modes for a given k (see Appendix A3 for explicit expressions). Fig. 3 shows examples of a geostrophic polynomial, ES and EA inertial modes, presenting the structure at the core surface and in meridional sections. Note that the inertial modes have little energy at small cylindrical radius and are well-suited for describing fluid motions in the equatorial region. Figure 3. Open in new tabDownload slide Examples of a geostrophic polynomial (a), and ES (b) and EA (c) inertial modes in a meridional plane (left-hand side of each panel) and on the core surface (right-hand side). The cyan line indicates the intersection of the meridional plane with the core surface. Note that the geostrophic polynomials are axisymmetric. Figure 3. Open in new tabDownload slide Examples of a geostrophic polynomial (a), and ES (b) and EA (c) inertial modes in a meridional plane (left-hand side of each panel) and on the core surface (right-hand side). The cyan line indicates the intersection of the meridional plane with the core surface. Note that the geostrophic polynomials are axisymmetric. 3.2 Decomposition of the flow into steady and time-dependent parts We model the core flow as a linear combination of a finite number of geostrophic polynomials and inertial modes. Further, we divide the flow into two parts, a steady background flow |$\mathbf {u}_0$| and a time-dependent flow |$\mathbf {u}_t$| \begin{eqnarray*} \mathbf {u}(\mathbf {r},t)=\mathbf {u}_0(\mathbf {r})+\mathbf {u}_t(\mathbf {r},t). \end{eqnarray*} (17) We allow the geostrophic mode, the QG modes (i.e. ES modes with n = 1) and the slowest EA modes (n = 1) to be time-dependent, since these are the modes expected to be excited by slow (relative to the rotation timescale) forcing processes in the core (e.g. Zhang & Liao 2017) \begin{eqnarray*} \begin{aligned} \mathbf {u}_t(\mathbf {r},t)=\sum _{k=1}^{K}a_k^\mathrm{G}(t)G_{2k-1}(\mathbf {r}){\hat{\boldsymbol \Phi }}+ \sum _{k=1}^{K}\sum _{m=1}^{M}a^\mathrm{S}_{m1k}(t)\Re \left[\mathbf {u}^\mathrm{S}_{m1k}(\mathbf {r})\right]+\sum _{k=0}^{K}\sum _{m=1}^{M}a^\mathrm{A}_{m1k}(t)\Re \left[\mathbf {u}^\mathrm{A}_{m1k}(\mathbf {r})\right] + \text{ imaginary parts}, \end{aligned} \end{eqnarray*} (18) where ℜ[ · ] denotes the real part. The steady background consists of the axisymmetric and non-axisymmetric inertial modes with n ≠ 1, since forcing processes on very long timescales, for example, from inhomogeneous mantle heat flow, is expected to include both equatorial symmetries \begin{eqnarray*} \begin{aligned} \mathbf {u}_0(\mathbf {r})=\sum _{k=1}^{K}\sum _{m=0}^{M}\sum _{n\ne 1}a^\mathrm{S}_{mnk}\Re \left[\mathbf {u}^\mathrm{S}_{mnk}(\mathbf {r})\right]+\sum _{k=0}^{K}\sum _{m=0}^{M}\sum _{n\ne 1}a^\mathrm{A}_{mnk}\Re \left[\mathbf {u}^\mathrm{A}_{mnk}(\mathbf {r})\right]+\text{ imaginary parts}. \end{aligned} \end{eqnarray*} (19) The upper summation limits K and M of the indexes k and m serve as the truncation degrees in this representation. The index n of the inertial modes is automatically bounded by the number of roots of the polynomials in eqs (13)–(16) for a given k. It is important to clarify at this point that we do not use the time-dependence of the inertial modes to parametrize the time-dependence of our core flows. We use only their spatial structure that we believe to be an efficient way of parametrizing the morphology of rotation-dominated motions at mid and low latitudes. As discussed in the next section, we solve for the structure of the flow independently at each epoch (i.e. every 4 months), except for a temporal regularization minimizing the flow acceleration. In effect, we are using the geostrophic and inertial modes to provide a kinematic description of the structure of the core flow epoch by epoch, and do not use them to describe dynamics. Forcing by buoyancy and magnetic effects, not accounted for in the above framework, will determine the time evolution of each mode, rather than the free oscillation periods of the modes. 4 FLOW ESTIMATION FROM GEOMAGNETIC OBSERVATIONS 4.1 Induction equation and unresolved scales The starting point for inferring the core flow from observations of the magnetic field is the radial component of the induction equation at the core surface, under the frozen-flux approximation (e.g. Bloxham & Jackson 1991) \begin{eqnarray*} \frac{\partial ^{}B_r}{\partial t^{}}=-\nabla _\mathrm{H}\cdot (\mathbf {u}B_r),\qquad \mbox{at } \quad r=1, \end{eqnarray*} (20) where |$\nabla _\mathrm{H}\equiv \nabla -\mathbf {\hat{r}} \frac{\partial ^{}}{\partial r^{}}$| is the horizontal gradient. At the core surface, only the large-scale structure of the core field and SV are resolved and can be used to study the core flow. Since all lengthscales of the core field and the flow can interact to give large-scale SV (Backus 1968), we decompose the radial core field as |$B_r=\overline{B}_r+\tilde{B}_r$|⁠, the sum of the resolved large-scale radial core field and an unresolved small-scale part. According to eq. (20) the large-scale radial SV at the core surface, denoted |$\frac{\partial ^{}\overline{B}_r}{\partial t^{}}$|⁠, is then (Gillet et al.2015a) \begin{eqnarray*} \frac{\partial ^{}\overline{B}_r}{\partial t^{}}=-\overline{\nabla _\mathrm{H}\cdot (\mathbf {u}\overline{B}_r)}+e, \end{eqnarray*} (21) where |$e=-\overline{\nabla _\mathrm{H}\cdot (\mathbf {u}\tilde{B}_r)}$| is the radial component of the so-called small-scale error that represents the contribution of the small-scale radial core field interacting with the flow. In order to use eq. (21) in an inversion scheme, we assume that the radial core field and SV can be derived from a scalar potential field that we expand into spherical harmonics (SH). Similarly, we express the flow at the core surface with toroidal and poloidal potentials in terms of an SH expansion. Inserting the SH expansions in eq. (21) and numerically integrating over a grid at the core surface (e.g. Lloyd & Gubbins 1990; Jackson 1997) results in the matrix equation \begin{eqnarray*} \mathbf {\dot{b}}=\mathbf {H}_b(\mathbf {b})\mathbf {x}+\mathbf {e}, \end{eqnarray*} (22) where the column vectors |$\mathbf {\dot{b}}$|⁠, |$\mathbf {b}$|⁠, |$\mathbf {x}$| and |$\mathbf {e}$| contain the SH expansions of the SV up to degree |$N_{\dot{b}}$|⁠, the resolved core field up to degree Nb, the toroidal–poloidal flow potentials at the core surface truncated at degree Nx and the small-scale error with truncation degree Ne, all evaluated at the core surface. Note that the matrix |$\mathbf {H}_b$| depends on the core field |$\mathbf {b}$|⁠, which we consider here to be known up to the truncation degree and given by a field model; this approximation is often made in such studies (e.g. Whaler & Beggan 2015). Eq. (22) is a global solution, however, here we use only observations from mid and low latitudes that are less contaminated by noise. Thanks to the form of the Green’s function for Laplace’s equation subject to Neumann boundary conditions, which links the core surface field to the field at the Earth’s surface and at satellite altitude, these observations do provide information at all latitudes. Nevertheless, the data constraint is weaker at polar latitudes and we thus concentrate on time-dependent flows at mid and low latitudes. 4.2 Forward problem According to eq. (22), the core surface flow and the small-scale error are connected to the large-scale radial SV at a single epoch. By making use of this equation at discrete times tp, p ∈ [1, P] spanning the time interval for which observational data is provided, we build a multi-epoch matrix equation to simultaneously compute the SV at all tp (e.g. Whaler et al.2016). For our formulation of the core flow in terms of modes, the forward problem is completed by including a matrix that allows the computation of the toroidal–poloidal flow expansion from a linear combination of geostrophic polynomials and inertial modes. Hereafter, the column vectors |$\mathbf {\dot{b}}_p$|⁠, |$\mathbf {b}_p$|⁠, |$\mathbf {x}_p$| and |$\mathbf {e}_p$| denote the SH expansions of the SV, the core field, the core surface flow and the small-scale error evaluated at time tp, respectively. We express the SH components of the time-dependent flow as linear combination of geostrophic polynomials, described by coefficients |$a^\mathrm{G}_k(t_p)$|⁠, QG modes, described by coefficients |$a^\mathrm{S}_{m1k}(t_p)$|⁠, and the slowest EA modes, described by the coefficients |$a^\mathrm{A}_{m1k}(t_p)$| in a time-dependent column vector |$\mathbf {a}_p\equiv \mathbf {a}(t_p)$|⁠. Similarly, we store the coefficients |$a^\mathrm{S}_{mnk}$| and |$a^\mathrm{A}_{mnk}$| with n ≠ 1 in a time-independent column vector |$\mathbf {a}_0$| that describes the remaining steady background flow. The toroidal–poloidal expansion of flow eq. (17) at time tp can therefore be written as \begin{eqnarray*} \mathbf {x}_p=\mathbf {M}_0\mathbf {a}_0+\mathbf {M}_t\mathbf {a}_{p}, \end{eqnarray*} (23) where |$\mathbf {M}_0$| and |$\mathbf {M}_t$| are fixed matrices that contain as columns the toroidal–poloidal expansion of the geostrophic and inertial modes matching the ordering of the coefficients in |$\mathbf {a}_0$| and |$\mathbf {a}_p$|⁠. Inserting this expression into eq. (22), the SH expansion of the SV on the core surface at time tp becomes \begin{eqnarray*} \mathbf {\dot{b}}_p=\mathbf {H}_b(\mathbf {b}_p)\mathbf {x}_p +\mathbf {e}_p=\mathbf {H}_b(\mathbf {b}_p)\bigg [\mathbf {M}_0\mathbf {a}_0+\mathbf {M}_t\mathbf {a}_p\bigg ] +\mathbf {e}_p. \end{eqnarray*} (24) In order to calculate from our model the three-vector components of the SV (i.e. ∂Br/∂t, ∂Bθ/∂t and ∂Bϕ/∂t) at the known locations of the ground and virtual observatories, listed in the form of a data vector |$\mathbf {d}_p$|⁠, we multiply the SV expansion |$\dot{\mathbf {b}}_p$| in eq. (24) with the matrix |$\mathbf {Y}_p$|⁠, containing the appropriate spatial derivatives of SH, such that \begin{eqnarray*} \mathbf {d}_p=\mathbf {Y}_p\mathbf {\dot{b}}_p=\mathbf {Y}_p\mathbf {H}_b(\mathbf {b}_p)\bigg [\mathbf {M}_0\mathbf {a}_0+\mathbf {M}_t\mathbf {a}_p\bigg ] +\mathbf {Y}_p\mathbf {e}_p. \end{eqnarray*} (25) In this step, it is assumed that the mantle conductivity is small and can be neglected. By adjusting the number of rows in |$\mathbf {Y}_p$|⁠, we can also handle the changing number of SV data over time. To summarize, we define our model vector to be |$\mathbf {m}=[\mathbf {a}_0^\mathrm{T},\mathbf {a}_1^\mathrm{T},\dots ,\mathbf {a}_P^\mathrm{T},\mathbf {e}_1^\mathrm{T},\dots ,\mathbf {e}_P^\mathrm{T}]^\mathrm{T}$| and combine the P single-epoch expressions in eq. (25) into one matrix equation of the form \begin{eqnarray*} \mathbf {d}=\mathbf {G}\mathbf {m}, \end{eqnarray*} (26) where the forward matrix is \begin{eqnarray*} \mathbf {G} = \left[\begin{array}{lcccccc} \mathbf {F}_1\mathbf {M}_0 & \mathbf {F}_1\mathbf {M}_t & & & \mathbf {Y}_1 & & \\ {\vdots } & & {\ddots } & & & {\ddots } & \\ \mathbf {F}_P\mathbf {M}_0 & & & \mathbf {F}_P\mathbf {M}_t & & & \mathbf {Y}_P \end{array}\right] \end{eqnarray*} (27) with |$\mathbf {F}_p=\mathbf {Y}_p\mathbf {H}_b(\mathbf {b}_p)$|⁠. This allows us to calculate model predictions |$\mathbf {d}=[\mathbf {d}_1^\mathrm{T},\dots ,\mathbf {d}_P^\mathrm{T}]^\mathrm{T}$| at the locations and times of the observational data. 4.3 Inverse problem In order to infer the core flow from the geomagnetic observations stored in |$\mathbf {d}^\mathrm{obs}$|⁠, we minimize a cost function of the model vector |$\mathbf {m}$| of the form \begin{eqnarray*} \Phi (\mathbf {m}) = (\mathbf {G}\mathbf {m}-\mathbf {d}^\mathrm{obs})^\mathrm{T}\mathbf {W}_d(\mathbf {G}\mathbf {m}-\mathbf {d}^\mathrm{obs}) +\lambda _0\Vert \mathbf {Q}_0\mathbf {a}_0 \Vert _1 + \lambda _t^\mathrm{S}\Vert \mathbf {Q}_t^\mathrm{S}\mathbf {a}_t \Vert _1 + \lambda _t^\mathrm{A}\Vert \mathbf {Q}_t^\mathrm{A}\mathbf {a}_t \Vert _1 + \lambda _a\Vert \mathbf {D}\mathbf {a}_t \Vert ^2_2 + \mathbf {e}_t^\mathrm{T}\mathbf {C}^{-1}_e\mathbf {e}_t, \end{eqnarray*} (28) where |$\Vert \mathbf {x} \Vert _p=(\sum _i\left|x_i\right|^p)^{1/p}$| denotes the lp norm of a vector |$\mathbf {x}$|⁠, |$\mathbf {a}_t\equiv [\mathbf {a}_{1}^\mathrm{T},\dots ,\mathbf {a}_{P}^\mathrm{T}]^\mathrm{T}$| is the vector containing all time-dependent flow coefficients and |$\mathbf {e}_t\equiv [\mathbf {e}_{1}^\mathrm{T},\dots ,\mathbf {e}_{P}^\mathrm{T}]^\mathrm{T}$| is the multi-epoch small-scale error. The strictly positive parameters λ0, |$\lambda _t^\mathrm{S}$|⁠, |$\lambda _t^\mathrm{A}$| and λa control the degree to which the steady background flow, the time-dependent equatorially symmetric and antisymmetric flow, and the flow acceleration are regularized. The first term on the right-hand side of eq. (28) measures the misfit to the ground and satellite SV data. It involves the diagonal matrix |$\mathbf {W}_d$| that contains the estimated data error variances |$\sigma _i^2$| modified following a Tukey biweight scheme in order to downweight outliers (e.g. Constable 1988), \begin{eqnarray*} (\mathbf {W}_d)_{ii} = \left\lbrace \begin{aligned} &\frac{\left(1-\big (\frac{r_i}{c\sigma _i}\big )^2\right)^2}{\sigma _i^2}, & &\frac{\left|r_i\right|}{\sigma _i} \le c\\ &0, & &\frac{\left|r_i\right|}{\sigma _i} \gt c \end{aligned} \right. \end{eqnarray*} (29) with breakpoint c = 4.685 and ri the residual of the data point |$(\mathbf {d}^\mathrm{obs})_i$| relative to the CHAOS-6-x7 geomagnetic field model. The choice of the data error variances |$\sigma _i^2$| is described in Section 2.2. We regularize the spatial structure of the flow by minimizing an l1 norm of the square root of the mode enstrophies calculated with the diagonal matrices |$\mathbf {Q}_0$| for the steady part |$\mathbf {u}_0$|⁠, |$\mathbf {Q}_t^\mathrm{S}$| and |$\mathbf {Q}_t^\mathrm{A}$| for the time-dependent equatorially symmetric and antisymmetric part of |$\mathbf {u}_t$|⁠. We regularize the time-dependence of the flow using an l2 norm of the flow acceleration calculated by applying a time finite-differencing operator |$\mathbf {D}$| to the time-dependent part of the flow in the model vector. The small-scale error is regularized using a quadratic form based on the inverse of its multi-epoch covariance matrix |$\mathbf {C}_e=\mathbf {C}_e(\mathbf {m})$|⁠, which depends on the model vector. We define the matrices |$\mathbf {D}$|⁠, |$\mathbf {Q}_0$|⁠, |$\mathbf {Q}_t^\mathrm{A}$|⁠, |$\mathbf {Q}_t^\mathrm{S}$| and |$\mathbf {C}_e$| in the next section. Regarding the minimization of the cost function in eq. (28), we note that this is a nonlinear problem as it involves l1 spatial norms of the flow model and the small-scale error covariance matrix that depends on the flow model. We therefore compute the model solution using an iterative procedure by solving the equation \begin{eqnarray*} \mathbf {m}_{k+1} = \left(\mathbf {G}^\mathrm{T}\mathbf {W}_d\mathbf {G}+\mathbf {R}(\mathbf {m}_k)\right)^{-1}\mathbf {G}^\mathrm{T}\mathbf {W}_d\mathbf {d}^\mathrm{obs} \end{eqnarray*} (30) at each iteration step k until the convergence criteria are met. We consider a computed model |$\mathbf {m}_{k+1}$| to be converged if the relative change between successive iterations of the data misfit |$(\mathbf {G}\mathbf {m}-\mathbf {d}^\mathrm{obs})^\mathrm{T}\mathbf {W}_d(\mathbf {G}\mathbf {m}-\mathbf {d}^\mathrm{obs})$|⁠, the l1 norm of the two parts of the flow |$\Vert \mathbf {Q}_t^\mathrm{S}\mathbf {a}_t \Vert _1+\Vert \mathbf {Q}_t^\mathrm{A}\mathbf {a}_t \Vert _1$| and |$\Vert \mathbf {Q}_0\mathbf {a}_0 \Vert _1$|⁠, and the small-scale error norm |$\mathbf {e}_t^\mathrm{T}\mathbf {C}^{-1}_e\mathbf {e}_t^\mathrm{\vphantom{T}}$| are each smaller than 10−2. The regularization matrix |$\mathbf {R}$| is block-diagonal and constructed from the model of the previous iterate |$\mathbf {m}_k$| according to \begin{eqnarray*} \mathbf {R}= \left[\begin{array}{lcc}\lambda _0\mathbf {W}_0 & & \\ & \lambda _t^\mathrm{S}\mathbf {W}_t^\mathrm{S}+\lambda _t^\mathrm{A}\mathbf {W}_t^\mathrm{A}+\lambda _a\mathbf {W}_a & \\ & & \mathbf {C}^{-1}_e \end{array}\right], \end{eqnarray*} (31) where the matrices |$\mathbf {W}_0$|⁠, |$\mathbf {W}_t^\mathrm{S}$| and |$\mathbf {W}_t^\mathrm{A}$| iteratively implement the l1 norm, and |$\mathbf {W}_a$| the l2 norm. With |$\mathbf {R}$| defined in this way, the model parameter estimation usually converged within 20 iterations. The definition of the matrices |$\mathbf {W}_0$|⁠, |$\mathbf {W}_t^\mathrm{S}$|⁠, |$\mathbf {W}_t^\mathrm{A}$|⁠, |$\mathbf {W}_a$| and |$\mathbf {C}_e$| is given in the following. 4.4 Regularization norms 4.4.1 Spatial regularization: l1 norm based on square root of mode enstrophies The flow coefficients |$\mathbf {a}_0$| and |$\mathbf {a}_t$| are spatially regularized using an l1 norm based on the Ekblom measure (Farquharson & Oldenburg 1998) \begin{eqnarray*} \mathbf {W}_0=\mathbf {Q}^\mathbf {T}_0\left(\frac{\delta _{ij}}{\sqrt{(\mathbf {Q}_0\mathbf {a}_0)_i^2+\epsilon ^2}}\right)\mathbf {Q}^\mathbf {\vphantom{T}}_0, \end{eqnarray*} (32) for the steady background part and similarly |$\mathbf {W}_t^\mathrm{S}$| and |$\mathbf {W}_t^\mathrm{A}$| for the time-dependent equatorially symmetric and antisymmetric part of the flow by replacing |$\mathbf {Q}_0$| and |$\mathbf {a}_0$| with |$\mathbf {Q}_t^\mathrm{S}$| and |$\mathbf {a}_t$|⁠, or |$\mathbf {Q}_t^\mathrm{A}$| and |$\mathbf {a}_t$|⁠, respectively. The constant ε > 0 is set to a value of 10−8 in order to prevent numerical problems in case the square root of the enstrophy of a mode approaches zero. In general, use of an l1 norm promotes a sparse model representation, by favouring a scenario where the model parameters are if possible zero, while reducing the amplitude of the remaining parameters. In order to also favour large-scale flows (to help the spatial power spectrum of our flow to converge), we apply the l1 norm with weights corresponding to the square root of the enstrophy |$Q(\mathbf {u})\ge 0$| for a flow |$\mathbf {u}$|⁠, with the enstrophy, the square of this quantity, defined as \begin{eqnarray*} Q^2(\mathbf {u})=\frac{3}{4\pi }\int _{\mathcal {V}}\left|\nabla \times \mathbf {u}\right|^2\mathrm{d}\mathcal {V}. \end{eqnarray*} (33) We compute the square root of the enstrophy for every individual inertial mode and geostrophic polynomial, and build the diagonal matrices |$\mathbf {Q}_0$| for the steady background, and |$\mathbf {Q}_t^\mathrm{S}$| and |$\mathbf {Q}_t^\mathrm{A}$| for the time-dependent equatorially symmetric and antisymmetric part of the flow. Expressions of the enstrophy of the modes can be found in the Appendix eqs (A4), (A8) and (A12). The parameters |$\lambda_t^\mathrm{S}$| and |$\lambda_t^\mathrm{A}$| control the extent to which the time-dependent part of the flow is regularized in space while λ0 controls the regularization of the steady background flow. 4.4.2 Temporal regularization: l2 norm of mode Acceleration We regularize the flow in time by minimizing the flow acceleration via an l2 norm based on the matrix |$\mathbf {W}_a=\mathbf {D}^\mathrm{T}\mathbf {D}$|⁠, where |$\mathbf {D}$| is the first-order forward time difference operator for equally spaced times tp \begin{eqnarray*} \mathbf {D}=\frac{1}{\Delta t} \left[\begin{array}{lccc}-\mathbf {I} & \quad\mathbf {I} \quad& \quad& \\ & \quad{\ddots } & \quad{\ddots } & \\ & \quad& -\mathbf {I} &\quad\mathbf {I} \end{array}\right], \end{eqnarray*} (34) where Δt = tp + 1 − tp is the step-size and |$\mathbf {I}$| is the unit matrix that is compatible with the size of |$\mathbf {a}_p$|⁠. The strength of the regularization of the flow acceleration is controlled by the parameter λa. 4.4.3 Regularization of the small-scale error We compute the multi-epoch covariance matrix |$\mathbf {C}_e$| of the small-scale error by projecting the covariance of the unresolved small-scale core field |$\mathbf {C}^\mathrm{\vphantom{T}}_{\tilde{b}}$| of SH degree |$N_b\lt n\le N_{\tilde{b}}$| onto large lengthscales of the SV, via the flow from the previous iteration \begin{eqnarray*} \mathbf {C}_e=\mathbf {H}^\mathrm{\vphantom{T}}\mathbf {C}^\mathrm{\vphantom{T}}_{\tilde{b}}\mathbf {H}^\mathrm{T}, \end{eqnarray*} (35) using the block-diagonal matrix \begin{eqnarray*} \mathbf {H}=\mathrm{diag}\big (\mathbf {H}_x(\mathbf {x}_1),\dots ,\mathbf {H}_x(\mathbf {x}_P)\big ). \end{eqnarray*} (36) The matrix |$\mathbf {H}_x(\mathbf {x}_p)$| is based on the induction equation eq. (20) and is similar to |$\mathbf {H}_b(\mathbf {b}_p)$| but takes the flow as the argument to connect the expansions of the magnetic field and the SV \begin{eqnarray*} \mathbf {\dot{b}}_p=\mathbf {H}_x(\mathbf {x}_p)\mathbf {b}_p. \end{eqnarray*} (37) We build the covariance matrix for the unresolved small-scale field |$\mathbf {C}_{\tilde{b}}$| as a P × P block matrix where each block |$(\mathbf {C}_{\tilde{b}})_{pq}$| for p, q ∈ [1, P] is equal to the covariance of the small-scale core field of epoch pairs tp and tq \begin{eqnarray*} (\mathbf {C}_{\tilde{b}})_{pq}=\mathrm{E}\left[(\tilde{\mathbf {b}}_p)^\mathrm{T}(\tilde{\mathbf {b}}_q)\right]=\underset{n}{\mathrm{diag}}\left(\sigma ^2_{\tilde{b}}(n)\left(1+\sqrt{3}\frac{\left|t_p-t_q\right|}{\tau _\mathrm{c}(n)}\right)\exp \left(-\sqrt{3}\frac{\left|t_p-t_q\right|}{\tau _\mathrm{c}(n)}\right)\right), \end{eqnarray*} (38) with E[ · ] the expectation, τc(n) the correlation time, |$\sigma ^2_{\tilde{b}}(n)$| the variance of the small-scale core field coefficient of degree n, which we have stored as elements in the vector |$\tilde{\mathbf {b}}_p$|⁠. Here, we use a Matérn class time correlation function with |$\nu =\frac{3}{2}$| (Rasmussen & Williams 2006). This has previously been used for a similar purpose by Gillet et al. (2015a); it describes an AR-2 process compatible with the observed secular variation spectrum (Lesur et al.2017) and allows abrupt slope changes in the SV (jerks). The parameters τc and |$\sigma ^2_{\tilde{b}}$| are estimated for each degree from the field and SV spectra of the CHAOS-6-x7 model in epoch 2015 (see Appendix B for details). In the computation, we use the time average of the flow over all considered epochs in |$\mathbf {H}$| since we expect the covariance of the small-scale error to be only weakly time-dependent when the flow is predominantly steady (Bloxham 1992; Gillet et al.2015a); this helps improve model convergence. 5 RESULTS In this section, we report results for our reference core flow model CoreFlo-LL.1. For this model the truncation degree of the geostrophic polynomials, the ES and EA inertial modes was chosen to be K = 10 and M = 20. To ensure an adequate representation of these modes in the computations an SH truncation degree Nx = 60 was used and the modes were evaluated to high precision using Matlab's symbolic toolbox. The multi-epoch inversion covers the period from 2000 September to 2018 January in 4-month steps, resulting in P = 53 epochs. The flow is a combination of one geostrophic mode (represented by 10 geostrophic polynomials) and 4720 inertial modes (200 QG modes and the 220 slowest EA modes in the time-dependent part, and the remaining 4300 ES and EA inertial modes in the steady background). The coefficients of the time-dependent part of the flow together with the small-scale error have to be estimated for every epoch. In total, we estimate 68 914 model coefficients: 2 × 4300 = 8600 determine the steady background flow, P × (10 + 2 × 200 + 2 × 220) = 45 050 the time-dependent part of the flow and P × 288 = 15 264 the SH coefficients of the small-scale error. The truncation degree of the small-scale error Ne = 16 is identical to the truncation degree |$N_{\dot{b}}$| of the computed SV. In order to compute the SV produced by a given core surface flow, we use the CHAOS-6-x7 model to estimate the core surface field only up to degree Nb = 14 since higher degrees are contaminated by the lithospheric field (Finlay et al.2016). In addition to our reference flow model CoreFlo-LL.1, we computed several other flows that will be discussed further in Section 6. Table 1 summarizes these flows and various statistics (defined below) derived from them. The additional flow models were designed to test the robustness of CoreFlo-LL.1 by changing aspects of the inversion scheme and the model parametrization. In model 1b, we constrained the time-dependent part of the flow to be purely equatorially symmetric (i.e. purely geostrophic and QG), whilst continuing to regularize in the same way as for CoreFlo-LL.1. In model 1c, we replaced our spatial flow regularization, based on the l1 norm of the square root of the enstrophy of the modes, with a more traditional approach using an l2 norm that scales with the third power of the spherical harmonic degree in the toroidal–poloidal expansion of the flow. Table 1. Summary of the computed flow models. Model 1 is our reference solution CoreFlo-LL.1, which is derived using an l1 norm of the square root of the mode enstrophies, and including time-dependent geostrophic, QG and the slowest EA modes (G+QG+EA). In model 1b, we use the same l1 regularization norm as for model 1, but restrict the time-dependence to the geostrophic mode and QG inertial modes (G+QG). For model 1c, which uses the same parametrization as model 1b, we use as spatial regularization an l2 norm of the toroidal–poloidal flow expansion that scales as n3 (Gillet et al.2009). Definitions of the diagnostics are given in the text. Model . Regularization norm . |$\lambda _t^\mathrm{S}$| . |$\lambda _t^\mathrm{A}$| . λ0 . λa . ΦSV . ΦSA . ΦLOD . ΦE . fS . ft . 1 l1 (G+QG+EA) 0.60 3.0 5.0 167 0.77 0.58 0.55 0.57 0.82 0.17 1b l1 (G+QG) 0.25 – 6.3 167 0.76 0.57 1.7 0.45 0.83 0.41 1c l2 (G+QG) 0.16 – 1.6 167 0.77 0.58 0.67 9.1 · 106 0.90 0.16 CHAOS-6-x7 – – – – – 0.80 0.64 – – – – Model . Regularization norm . |$\lambda _t^\mathrm{S}$| . |$\lambda _t^\mathrm{A}$| . λ0 . λa . ΦSV . ΦSA . ΦLOD . ΦE . fS . ft . 1 l1 (G+QG+EA) 0.60 3.0 5.0 167 0.77 0.58 0.55 0.57 0.82 0.17 1b l1 (G+QG) 0.25 – 6.3 167 0.76 0.57 1.7 0.45 0.83 0.41 1c l2 (G+QG) 0.16 – 1.6 167 0.77 0.58 0.67 9.1 · 106 0.90 0.16 CHAOS-6-x7 – – – – – 0.80 0.64 – – – – Open in new tab Table 1. Summary of the computed flow models. Model 1 is our reference solution CoreFlo-LL.1, which is derived using an l1 norm of the square root of the mode enstrophies, and including time-dependent geostrophic, QG and the slowest EA modes (G+QG+EA). In model 1b, we use the same l1 regularization norm as for model 1, but restrict the time-dependence to the geostrophic mode and QG inertial modes (G+QG). For model 1c, which uses the same parametrization as model 1b, we use as spatial regularization an l2 norm of the toroidal–poloidal flow expansion that scales as n3 (Gillet et al.2009). Definitions of the diagnostics are given in the text. Model . Regularization norm . |$\lambda _t^\mathrm{S}$| . |$\lambda _t^\mathrm{A}$| . λ0 . λa . ΦSV . ΦSA . ΦLOD . ΦE . fS . ft . 1 l1 (G+QG+EA) 0.60 3.0 5.0 167 0.77 0.58 0.55 0.57 0.82 0.17 1b l1 (G+QG) 0.25 – 6.3 167 0.76 0.57 1.7 0.45 0.83 0.41 1c l2 (G+QG) 0.16 – 1.6 167 0.77 0.58 0.67 9.1 · 106 0.90 0.16 CHAOS-6-x7 – – – – – 0.80 0.64 – – – – Model . Regularization norm . |$\lambda _t^\mathrm{S}$| . |$\lambda _t^\mathrm{A}$| . λ0 . λa . ΦSV . ΦSA . ΦLOD . ΦE . fS . ft . 1 l1 (G+QG+EA) 0.60 3.0 5.0 167 0.77 0.58 0.55 0.57 0.82 0.17 1b l1 (G+QG) 0.25 – 6.3 167 0.76 0.57 1.7 0.45 0.83 0.41 1c l2 (G+QG) 0.16 – 1.6 167 0.77 0.58 0.67 9.1 · 106 0.90 0.16 CHAOS-6-x7 – – – – – 0.80 0.64 – – – – Open in new tab We characterize and compare the derived flows with the help of the following statistics: normalized quadratic measure of the SV misfit between the observational data |$\mathbf {d}^\mathrm{obs}$| and our model predictions |$\mathbf {d}=\mathbf {G}\mathbf {m}$| \begin{eqnarray*} \Phi _\mathrm{SV} = \frac{1}{N_d}(\mathbf {d}-\mathbf {d}^\mathrm{obs})^\mathrm{T}\mathbf {W}_d(\mathbf {d}-\mathbf {d}^\mathrm{obs}), \end{eqnarray*} (39) where Nd is the number of data, normalized quadratic measure of the SA misfit by computing the central time difference |$\mathbf {\dot{d}}$| and |$\mathbf {\dot{d}}^\mathrm{obs}$| of the model predictions and data \begin{eqnarray*} \Phi _\mathrm{SA} = \frac{1}{N_{\dot{d}}}(\mathbf {\dot{d}}-\mathbf {\dot{d}}^\mathrm{obs})^\mathrm{T}\mathbf {W}_{\dot{d}}(\mathbf {\dot{d}}-\mathbf {\dot{d}}^\mathrm{obs}), \end{eqnarray*} (40) where |$N_{\dot{d}}$| is the number of the SA data, length-of-day (LOD) misfit between our model predictions and independent estimates |$\mathbf {d}_\gamma$| from Gillet et al. (2015b) \begin{eqnarray*} \Phi _\mathrm{LOD} = \frac{1}{\sigma _\gamma ^2N_\gamma }(\mathbf {G}_\gamma \mathbf {m}-\mathbf {d}_\gamma )^\mathrm{T}(\mathbf {G}_\gamma \mathbf {m}-\mathbf {d}_\gamma ), \end{eqnarray*} (41) where |$\sigma _\gamma =0.24\, \mathrm{ms}$| is the standard deviation of the LOD over the considered time period and |$\mathbf {G}_\gamma$| is the matrix to compute the LOD from our modelled core flow (e.g. Holme 2015), l1 norm of the square root of the mode enstrophy as a measure of spatial complexity \begin{eqnarray*} \Phi _\mathrm{E} = \frac{1}{N}(\Vert \mathbf {Q}_t^\mathrm{S}\mathbf {a}_t\Vert _1+\Vert \mathbf {Q}_t^\mathrm{A}\mathbf {a}_t\Vert _1+\Vert \mathbf {Q}_0\mathbf {a}_0\Vert _1), \end{eqnarray*} (42) where N is the total number of flow coefficients, fraction fS of the total time-averaged flow power that is contributed by the equatorially symmetric flow |$\mathbf {u}^\mathrm{S}$| (geostrophic mode and ES inertial modes) integrated over the core surface \begin{eqnarray*} f_\mathrm{S} = \frac{\left\langle \int _{\Omega }\left|\mathbf {u}^\mathrm{S}\right|^2\mathrm{d}\Omega \right\rangle }{\left\langle \int _{\Omega }\left|\mathbf {u}\right|^2\mathrm{d}\Omega \right\rangle }, \end{eqnarray*} (43) where 〈 · 〉 stands for the time-average, and fraction ft of the total time-averaged flow power that is contributed by the time-dependent part of the flow |$\mathbf {u}_t$| (geostrophic mode, QG inertial modes, and the slowest EA inertial modes for CoreFlo-LL.1), integrated over the core surface \begin{eqnarray*} f_\mathrm{t} = \frac{\left\langle \int _{\Omega }\left|\mathbf {u}_t\right|^2\mathrm{d}\Omega \right\rangle }{\left\langle \int _{\Omega }\left|\mathbf {u}\right|^2\mathrm{d}\Omega \right\rangle }. \end{eqnarray*} (44) 5.1 Fit to observations In Fig. 1 we present the median absolute deviation between CoreFlo-LL.1 model predictions and the observed SV at ground and virtual observatories as a function of location, for the radial, north–south and east–west SV components, respectively. The data at the majority of the sites is fit to better than 2 nT yr−1, with the exception of some ground stations, for example those in India, where the mean residuals are somewhat larger. There is however no obvious systematic trend in the amplitude of the residuals, and the data at both low and mid latitudes is overall fit well, suggesting that CoreFlo-LL.1 provides a reasonable explanation of field variations at these latitudes. Fig. 2 shows time-series of ground and SV data from example low-latitude sites, compared with the predictions from the CHAOS-6-x7 field model (red) and CoreFlo-LL.1 (green). CoreFlo-LL.1 predictions fit the data as well as CHAOS-6-x7 and it predicts remarkably similar patterns of SV. There is no evidence for systematic offsets in the CoreFlo-LL.1 series, and it succeeds in capturing ‘V’ shaped structures in the SV associated with geomagnetic jerks, see for example the ϕ-component at MBO in 2007 (e.g. Chulliat et al.2010) and the VO at latitude 5.964° and longitude −32.995°. The time-dependence in our SV predictions comes almost entirely from changes in the time-dependent part of our flow, with the small-scale error term simply adding an almost constant offset to the prediction of the steady part of the flow at a given location. Thus in CoreFlo-LL.1 changes of the SV (field accelerations) are largely explained by changes (accelerations) in the large-scale flow interacting with the large-scale core field. We come back to this point in Section 6.1. Histograms of the residuals between the SV data and CoreFlo-LL.1 predictions, over the entire time period, are presented in Fig. 4. We separate ground and virtual observatories as well as the corresponding vector components. The distributions are mostly symmetric but more long-tailed than for a Gaussian distribution of errors; this has been accounted for during the model estimation by use of a robust Tukey biweight scheme in the misfit function. The radial, north–south and east–west SV residuals have rms values of 2.07, 2.30 and 2.30 nT yr−1 at the ground observatories and 1.20, 1.42 and 1.21  nT yr−1 at the virtual observatories. These histograms provide further evidence that CoreFlo-LL.1 is describing well the observed SV at mid and low latitudes. Figure 4. Open in new tabDownload slide Histograms of the residuals between predictions of the CoreFlo-LL.1 model and the SV at ground and satellite altitude. For comparison, the dashed curves show the probability density function of Gaussian distributed residuals given the respective mean values and variances. Figure 4. Open in new tabDownload slide Histograms of the residuals between predictions of the CoreFlo-LL.1 model and the SV at ground and satellite altitude. For comparison, the dashed curves show the probability density function of Gaussian distributed residuals given the respective mean values and variances. In Fig. 5, we show the SV power spectrum, taking 2015 September as an example, a time well within a period of good data coverage provided by the Swarm satellites. Figure 5. Open in new tabDownload slide Power spectrum at the core surface in 2015 September of the SV computed from CoreFlo-LL.1 (green), the SV generated by the modelled flow alone, that is, excluding the small-scale error (blue), the small-scale error (black), the CHAOS-6-x7 model (red) and the difference between our prediction and the CHAOS-6-x7 model (grey). Figure 5. Open in new tabDownload slide Power spectrum at the core surface in 2015 September of the SV computed from CoreFlo-LL.1 (green), the SV generated by the modelled flow alone, that is, excluding the small-scale error (blue), the small-scale error (black), the CHAOS-6-x7 model (red) and the difference between our prediction and the CHAOS-6-x7 model (grey). The SV spectrum at the core surface estimated in the CHAOS-6-x7 model (red) is well reproduced by CoreFlo-LL.1 (green), up until SH degree 11 where the power of the difference (grey) becomes of the same order of magnitude as the SV itself. The power in the small-scale error (black) is consistently less than the power in the flow-generated SV (blue), which confirms that the majority of the observed SV is reproduced by the CoreFlo-LL.1 flow interacting with the large-scale core field. 5.2 Time-averaged flow In Fig. 6 we present the flow from CoreFlo-LL.1 time-averaged over the modelled time interval. The regions poleward of the tangent cylinder in this and subsequent plots have been masked in grey to emphasize that the determined flows are not reliable in this region. The basic time-averaged flow structure is that of a planetary-scale eccentric gyre familiar from previous studies (e.g. Pais & Jault 2008), with equatorward flow around 100°E, then strong and large-scale westward flow at mid and low latitudes under the Atlantic, stronger in the southern hemisphere in agreement with previous findings (Amit & Pais 2013; Baerenzung et al.2016) and poleward flow under North America. Note that because our static background flow is described by both ES and EA flow modes, departures from equatorial symmetry, including equator crossing, are allowed. The basic gyre structure, including its meridional flows, is recovered despite the fact that we only use SV data from mid and low latitudes. We find that |$82\, \hbox{ per cent}$| of the time-averaged power at the core surface results from equatorially symmetric flow, which is remarkably close to the value found by Baerenzung et al. (2014) via a probabilistic inversion scheme. Fig. 7 shows the toroidal–poloidal SH power spectrum of the time-averaged flow as well as several examples of the time-varying flow. The toroidal spectrum decreases gradually and the poloidal spectrum is relatively flat up to degree 12, both decrease more rapidly beyond degree 15. Figure 6. Open in new tabDownload slide Time-averaged flow at the core surface. The grey area masks the region inside the tangent cylinder where our flow parametrization is inappropriate. The projection is Hammer-Aitoff. Figure 6. Open in new tabDownload slide Time-averaged flow at the core surface. The grey area masks the region inside the tangent cylinder where our flow parametrization is inappropriate. The projection is Hammer-Aitoff. Figure 7. Open in new tabDownload slide Toroidal–poloidal power spectrum of the time-averaged flow (thick lines) and time-varying flow (thin lines, every 2-yr starting in 2000 September). Figure 7. Open in new tabDownload slide Toroidal–poloidal power spectrum of the time-averaged flow (thick lines) and time-varying flow (thin lines, every 2-yr starting in 2000 September). 5.3 Time-varying flow By removing the time-average from the instantaneous flow, we extract the time-varying part of the flow, which is by construction primarily equatorially symmetric in CoreFlo-LL.1. Looking at the amplitudes of the modes representing our time-dependent flow we find that the large-scale geostrophic polynomials (small k) and large-scale QG modes (small m and k) contribute most to the time-dependent part of the flow (93% of power of the time-varying flow is equatorially symmetric). More complex modes (high m and k) are in contrast less pronounced; this is at least partly a consequence of the applied spatial regularization which penalizes modes according to the square root of the enstrophy. In Fig. 8 snapshots of the time-varying flow at 3-yr intervals are presented to illustrate how the flow changes over the studied time interval. Figure 8. Open in new tabDownload slide Snapshots of the time-varying core surface flow (arrows) and the related horizontal divergence (contours) during the observation period from 2000 to 2018. The projection is Hammer-Aitoff and the grey area masks the region inside the tangent cylinder. Figure 8. Open in new tabDownload slide Snapshots of the time-varying core surface flow (arrows) and the related horizontal divergence (contours) during the observation period from 2000 to 2018. The projection is Hammer-Aitoff and the grey area masks the region inside the tangent cylinder. We find strong time-dependent, non-zonal, azimuthal flow in the equatorial region between latitudes 15°N and S, which is particularly enhanced under Indonesia (around 120°E) and under the northern part of South America (around 80°W). The overall amplitude of the time-varying flow apparently increases in the snapshots after 2014. This is a data-driven effect, since after 2014 is when high-quality data from the Swarm satellite mission are available. Also noteworthy is that the time alternating patterns of converging and diverging azimuthal flow around the equatorial region are associated with regions of downwelling and upwelling, respectively. Movies of the time-varying flow show that flow under northern South America (60 to 80°W) was converging in 2004, diverging in 2008, weakly converging in 2011 and then more strongly diverging in 2014 and strongly converging in 2017. See also Fig. 9 which is described in greater detail below. As an independent test of the time-variations of the geostrophic part of CoreFlo-LL.1, we computed the predicted change in the LOD and compared with geodetic observations3, corrected for known atmospheric and oceanic signals, and filtered using a 1-yr moving average (as in Gillet et al.2015b). We applied a moving average to the LOD changes predicted by our model to be able to compare with the observations. We find a trend of generally increasing LOD from 2003 to 2016 in agreement with the trend in the observations. About this trend our flow predicts a maximum in the change of LOD in 2007 and again in 2016, similar to the observations, but with lower amplitude especially in 2007. In order to more clearly visualize the interesting converging and diverging patterns of azimuthal flow at low latitude, we present the average of the azimuthal component of the time-varying flow over 15°N, S latitude in Fig. 9 as a function of longitude and time. Around 130°E, a prominent region of initially converging azimuthal flow changes sign to become an azimuthal flow divergence around 2011. The azimuthal flow under northern South America (60 to 90°W) changes its direction more frequently in time around a relatively fixed location on the equator. There is some indication of westward propagation of azimuthal flow features with alternating signs between 45°E and 60°W especially from 2005 to 2015. In the next section we discuss in detail how this time-varying flow gives rise to SA pulses and geomagnetic jerk events. Figure 9. Open in new tabDownload slide Time-varying part of the azimuthal flow, averaged over low latitudes, as a function of longitude and time. The azimuthal flow was averaged over latitudes 15°N, S of the equator in approx. 1° bins of longitude. Figure 9. Open in new tabDownload slide Time-varying part of the azimuthal flow, averaged over low latitudes, as a function of longitude and time. The azimuthal flow was averaged over latitudes 15°N, S of the equator in approx. 1° bins of longitude. 6 DISCUSSION 6.1 The origin of core surface field acceleration pulses The future evolution of the geomagnetic field is difficult to predict due to intermittent pulses of field acceleration that cause linear extrapolations to fail. It is therefore of significant interest to investigate whether our time-dependent flows can reproduce observed pulses of field acceleration, and if so by what process. In this respect it is important to recall that our flows were not derived from SV computed from spherical harmonic field models (e.g. CHAOS or GRIMM type models where the applied regularization plays an important role in the time-dependence of the core surface field acceleration), but directly from ground and satellite-based SV time-series. Fig. 10 presents the observed field accelerations at Earth’s surface and at satellite altitude, calculated by taking annual differences of the SV series from Fig. 2, together with the field accelerations predicted by the CoreFlo-LL.1 model and the CHAOS-6-x7 model. The CoreFlo-LL.1 predictions match the observations well with a weighted rms misfit of 2.75, 3.42 and 3.42 nT yr−2 to the accelerations in Br, Bθ and Bϕ, compared to 3.00, 3.49, 3.50 nT yr−2 for CHAOS-6-x7. Notable pulses of high-amplitude SA are seen in ∂2Bϕ/∂2t at MBO in 2006 and 2009, and in ∂2Br/∂2t in GUA near 2013. The SA observed at the satellite-based virtual observatories, including those at low latitudes, as shown for example in Fig. 10, are also generally well fit. From this we conclude that CoreFlo-LL.1 accounts well for mid- and low-latitude SA observed at ground and satellite altitude. Figure 10. Open in new tabDownload slide SA time-series for the ground and virtual observatories in Fig. 2 showing the observed SA (black), the CHAOS-6-x7 SA prediction (red) and our model’s SA prediction (green) both truncated at spherical harmonic degree 9. Figure 10. Open in new tabDownload slide SA time-series for the ground and virtual observatories in Fig. 2 showing the observed SA (black), the CHAOS-6-x7 SA prediction (red) and our model’s SA prediction (green) both truncated at spherical harmonic degree 9. Next, we examine the occurrence and form of the SA pulses at the core surface comparing those predicted by our flow model with those found in CHAOS-6-x7. Fig. 11 presents the time evolution of the SA power integrated over the core surface (e.g. Finlay et al.2015), calculated up to SH degree 9, from our modelled flow alone (i.e. without small-scale error contribution) and, for comparison, the same quantity from the CHAOS-6-x7 model. Peaks of the SA power around 2006 January, 2009 January, 2013 September and 2016 September correspond to well-known SA pulses at the core surface (Chulliat et al.2010; Chulliat & Maus 2014; Finlay et al.2015). There are, however, differences in the times of the pulse peaks of up to a year compared to CHAOS-6-x7. This may be, at least in part, due to the strong temporal regularization applied when constructing CHAOS-6-x7 (and similar models) which results in an averaging in time of the true field that depends on the spherical harmonic degree (Olsen et al.2009), and also because our flow model does not predict well the SA at high latitudes that will contribute to the SA power in the case of the field model. Figure 11. Open in new tabDownload slide SA power for SH degrees n ≤ 9 predicted by the CoreFlo-LL.1 flow at the core surface as a function of time (black), 1-yr moving average (green) and CHAOS-6-x7 prediction (red). The vertical lines indicate times of enhanced SA power, for which we plot model snapshots in Fig. 12. Figure 11. Open in new tabDownload slide SA power for SH degrees n ≤ 9 predicted by the CoreFlo-LL.1 flow at the core surface as a function of time (black), 1-yr moving average (green) and CHAOS-6-x7 prediction (red). The vertical lines indicate times of enhanced SA power, for which we plot model snapshots in Fig. 12. Fig. 12 presents the corresponding maps of the core surface SA plotted together with the core surface flow acceleration from CoreFlo-LL.1 at pulse times. The SA has again been truncated here at spherical harmonic degree 9 in order to facilitate comparisons with the CHAOS-6-x7 field model shown in the right-hand column. CoreFlo-LL.1 generates localized patches of strong SA in the equatorial region, for example, under northern South America, which change polarity over time in agreement with the patterns seen in CHAOS-6-x7. We find that pulses of SA power in Fig. 11 correspond to times of high-amplitude equatorial SA patches that alternate in sign and are located in specific longitudinal sectors. There are, however, differences in the patterns of core surface SA predicted by CoreFlo-LL.1 and those found in the CHAOS-6-x7 model, despite that both explain well the observed SA at Earth’s surface and at satellite altitude (recall Fig. 10). This suggests that care should perhaps be taken interpreting the small scale details of the core surface SA in field models, since a variety of possibilities exist for adequately explaining the observations, depending on the chosen regularization and other constraints imposed during model construction. Our time-dependent flow is by construction primarily equatorially symmetric—this leads to core surface SA patterns that are more equatorially symmetric than those seen in CHAOS-6-x7. Figure 12. Open in new tabDownload slide Snapshots of the flow acceleration and predicted SA from CoreFlo-LL.1 (left-hand side) during times of increased SA power as indicated in Fig. 11. The CHAOS-6-x7 prediction for the SA is also shown for comparison at the same times (right-hand side). SA has been truncated at spherical harmonic degree 9 in each case. The grey area masks the region inside the tangent cylinder for the flow predictions. Figure 12. Open in new tabDownload slide Snapshots of the flow acceleration and predicted SA from CoreFlo-LL.1 (left-hand side) during times of increased SA power as indicated in Fig. 11. The CHAOS-6-x7 prediction for the SA is also shown for comparison at the same times (right-hand side). SA has been truncated at spherical harmonic degree 9 in each case. The grey area masks the region inside the tangent cylinder for the flow predictions. Figs 12 and 13 illustrate how SA pulses in CoreFlo-LL.1 result from localized accelerations of the non-zonal azimuthal flow. In 2006 the strongest SA at low latitudes was situated under the northeastern corner of South America and under the equatorial Atlantic. At this time in CoreFlo-LL.1, there was a generally eastward acceleration under northern South America and westward flow acceleration under equatorial Africa. In 2009 the strongest SA feature was again under north eastern South America (around 45 °W), just slightly to the west of the location of the maximum amplitude feature in 2006 and opposite in sign. In CoreFlo-LL.1, this is associated with a strong westward acceleration under north eastern South America. The SA patterns under northern South America alternated in sign for the subsequent pulses in 2013 and 2016, as did the direction of the non-zonal azimuthal flow acceleration in this location. There is a clear divergence of azimuthal flow acceleration under northern South America in 2013. In 2016, when there is excellent Swarm data available for a year on either side, the flow accelerations in this region had reversed direction to consist of a strong convergence. Figure 13. Open in new tabDownload slide Flow acceleration at the core surface (top), on two meridional planes and the equatorial plane (middle) cutting across the core. The cyan lines at 87°W and 121°E indicate where the meridional planes intersect with the equatorial plane. The colour scale shows the strength of the azimuthal flow acceleration, while the grey patch masks the area of the inner core. We also plot the azimuthal flow acceleration averaged over 15°N, S latitude as a function of time and longitude (bottom). The green lines coincide with the time of increased SA as shown in Figs 11 and 12. Figure 13. Open in new tabDownload slide Flow acceleration at the core surface (top), on two meridional planes and the equatorial plane (middle) cutting across the core. The cyan lines at 87°W and 121°E indicate where the meridional planes intersect with the equatorial plane. The colour scale shows the strength of the azimuthal flow acceleration, while the grey patch masks the area of the inner core. We also plot the azimuthal flow acceleration averaged over 15°N, S latitude as a function of time and longitude (bottom). The green lines coincide with the time of increased SA as shown in Figs 11 and 12. The spatial structure of the non-zonal azimuthal flow acceleration at the core surface in 2013 is illustrated in greater detail in Fig. 13 where the azimuthal flow acceleration is contoured. The low-latitude flow acceleration is directed from regions of acceleration divergence towards regions of acceleration convergence. In order to satisfy the conservation of mass the former must be related with accelerations in upwelling and the latter to regions of acceleration in downwelling. The relevance of flow acceleration upwellings and downwellings to understanding rapid field accelerations has previously been highlighted by Olsen & Mandea (2008). The bottom panel in Fig. 13 shows a time longitude plot of the azimuthal flow accelerations in CoreFlo-LL.1, averaged 15°N and S of the equator. Times of SA acceleration pulses, marked by the green lines, correspond to times when there were strong azimuthal flow accelerations at specific locations. The longitude sectors from 60°W to 100°W, and from 80°E to 130°E, have experienced a succession of strong azimuthal flow acceleration features since 2000 with alternating sign. The Pacific region, from 130°E to 150°W, has also experienced strong azimuthal flow accelerations since 2012, with a change in sign from eastward to westward acceleration during 2014. Sign changes in the azimuthal flow acceleration are apparent at times and locations where geomagnetic jerks have been reported, for example, in the Pacific region during 2014 (Torta et al.2015). Sign changes in the azimuthal flow acceleration correspond to times when the azimuthal flow acceleration is changing rapidly; the time between strongly positive and negative accelerations can be short, of the order of 1 yr or less. Under the equatorial Pacific in 2014, as seen by comparing Figs 9 and 13, the sign change in the acceleration occurred as the time-varying azimuthal flow reached its highest amplitude, in this case a strong eastward flow. Not all zero crossings of the non-zonal azimuthal flow acceleration (or extrema in the time-varying, non-zonal, azimuthal flow) are associated with jerk-type events at Earth’s surface. In order for jerks to have a significant effect at the Earth’s surface, the flow acceleration sign changes must involve sufficiently large-scale flow structures, or similar flow acceleration sign changes must occur in several locations at about the same time. Following this logic, strong global jerk events seen in historical records (e.g. in 1969, 1978) could be the result of sign changes in the non-zonal azimuthal flow acceleration that are of sufficiently large amplitude and spatial extent. More localized jerk events can be due to changes in the sign of the non-zonal azimuthal flow acceleration structures with smaller lengthscale in specific locations. Overall, a spectrum of jerk events of different sizes and amplitudes is therefore expected, depending on the spatial extent and amplitude of the changes in the underlying non-zonal azimuthal flow acceleration. In Fig. 13 we also explore the 3-D structure of the flow accelerations as captured by CoreFlo-LL.1. Since the modes on which the model is based are 3-D and mode amplitudes have been inferred from the observations, we can present the flow and its accelerations throughout the core volume. At the core surface, we find the non-zonal azimuthal accelerations required to explain the observations are strongly localized in a belt within 15° of the equator and consist of successive regions of convergence and divergence, with enhanced amplitude at specific longitudes, for example, under the region west of northern South America in 2013. The flow acceleration structures of largest scale in 2013 are under the Pacific region, close to where there was a change in the sign of the acceleration in 2014. Low-latitude azimuthal flow acceleration is not present at all longitudes; it is weaker under the central Pacific and west of Africa in 2013. We find that the acceleration changes sign moving away from the equator. For example, under the Pacific region from 150°E to 180°E, the azimuthal acceleration is eastwards on the equator in 2013, but westwards at plus and minus 15° latitude. The time-dependent modes in CoreFlo-LL.1 are primarily equatorially symmetric and nearly axially invariant, hence meridional sections of the azimuthal flow acceleration show essentially columnar structures, with highest amplitudes close to the core surface at the equator. The sign change of the acceleration moving away from the equator at the core surface is seen in the meridional sections to be due to a change in sign of the columnar structures of the azimuthal flow acceleration in the cylindrical radial direction. Strong azimuthal flow accelerations are focused in a small part of the volume of the outer core, close to the core surface, as can be seen in the equatorial section. Considering a time-longitude plot of the azimuthal flow acceleration, alternation of the non-zonal azimuthal flow acceleration is evident, for example, close to 90°W and 120°E (cyan lines); local intensifications of the azimuthal flow acceleration are found at times of SA pulses (green lines). 6.2 Experiments with alternative flow regularizations and parametrizations In order to address the question of the robustness of the time-varying flow in CoreFlo-LL.1, and the associated accelerations, we have carried out a number of experiments with different choices of flow parametrization and model regularization, as summarized in Table 1, see the related discussion in Section 5 for a description of these experimental models. We examined in detail the azimuthal flow accelerations for these experimental flows. In equatorial time-longitude plots constructed as in Fig. 13., model 1band model 1c show many of the same flow acceleration patterns found for CoreFlo-LL.1. In particular, around 90°W we see unambiguous alternation in the direction of the azimuthal flow acceleration and we find that the intense eastward and westward accelerations again occur in similar pairs of converging and diverging accelerations. The amplitudes of the accelerations are higher for model 1b and there are more small-scale structures compared with CoreFlo-LL.1, as might be expected for a model constructed with equatorial symmetry of flow acceleration strictly imposed. Nevertheless, it consists of the same essential details: the non-zonal azimuthal flow acceleration alternates in time, for example around 90°W, or under the Pacific after 2012, and is organized around centres of flow acceleration convergence and divergence. In both model 1b and model 1c the azimuthal flow accelerations remain localized close to the core surface, and they vary little in the axial direction. Time-longitude plots made with models 1band1c also suggest a coherent westward movement of flow acceleration features from 90°W to 90°E between 2003 and 2014. Overall we find that flows regularized in different ways, or with different assumptions made regarding the equatorial symmetry of the time-varying flow, produce oscillating SA pulses in a similar manner to CoreFlo-LL.1, via alternating non-zonal azimuthal flow acceleration at low latitudes. Although the basic patterns of azimuthal flow acceleration are similar across the investigated flows, if one is willing to accept some time-dependent cross-equatorial flow (as in CoreFlo-LL.1) then less vigorous accelerations of flow convergence and divergence, and smaller amplitude flow accelerations, are required. Comparisons with numerical dynamo models in a regime appropriate for studying rapid dynamics are required to determine whether or not such transient EA flow accelerations are physically reasonable. 6.3 Relation to previous studies Patterns of time-varying core flow at low latitudes have been examined in a number of previous studies. Jackson (1997) mentioned the presence of oscillatory motions at low latitudes in his tangentially geostrophic flows, but did not analyse the motions in detail. Gillet et al. (2015b) inverted the COV-OBS field model for QG core flows and found vigorous azimuthal jets in the equatorial belt below 10° latitude; these were particularly intense after 2000 when satellite observations contributed to the field model. Gillet et al. (2015b) argued that the dynamics in the equatorial region was different to that occurring at mid latitudes, though coupled via meridional flows arriving or leaving at longitudes of 130°E and 60°W. They found interannual variations of the equatorial jets, for example at 85°W. Finlay et al. (2016) carried out a similar inversion of the higher resolution CHAOS-6 model, which included Swarm data, and found similar oscillating, non-zonal, azimuthal jets in the equatorial region, for example at 40°W. They noted that pulses of SA coincided with times of large azimuthal flow acceleration, between maxima and minima of the time-varying azimuthal flow. Our results are in good accord with these previous QG studies, despite our flows being parametrized and regularized in a very different way, and being based directly on ground and observatory data rather than on field models. Our flow model 1c is most similar to that of Gillet et al. (2015b) and Finlay et al. (2016) since it has purely equatorially symmetric time-dependent flows and a spatial regularization applied to its poloidal and toroidal flow components. Barrois et al. (2017) recently introduced a new Kalman filter approach to core flow estimation, taking output from a numerical dynamo simulation as prior information, and the COV-OBS.x1 field model as input. Barrois et al. (2018b) have extended this approach to work with ground and satellite data sets similar to those employed here. We have compared our flow results with those of Barrois et al. (2018b) as well as the corrected flows recently presented by Barrois et al. (2018a), in particular by producing the same time–longitude plots of azimuthal flow acceleration; we find Barrois’s flows show similar large-scale patterns of azimuthal flow acceleration to CoreFlo-LL.1, with regions of azimuthal flow acceleration convergence and divergence alternating on interannual timescales. The amplitudes of the accelerations in the models of Barrois are, however, smaller than those found in CoreFlo-LL.1; this may be in part due to their reliance on a dynamo prior which may lack realistic rapid dynamics at low latitudes (Aubert 2018). It is possible that magnetic diffusion, due to field gradients being built up by the time-varying core flow, could also make a contribution to the observed SA. This could conceivably be accounted for following the approach described by Amit & Christensen (2008) where diffusion is correlated with horizontal divergence at the core surface. Note, however, that in this case diffusion is enslaved to the flow and cannot drive field changes on its own. Another approach to studying the recent time-varying core flow has been pursued by Whaler et al. (2016). They inverted ground observatory data alone, applying a strong norm on the spatial complexity of the flow (varying with n5, forcing the flow to be large scale), but without enforcing any symmetry constraint or any requirement that the flow be quasi or tangentially geostrophic. Like us, they also minimized a measure of the flow acceleration. Whaler et al. (2016) obtained flows that were primarily toroidal, equatorially symmetric and tangentially geostrophic, with the time-averaged part explaining much of the variance, and being dominated by an eccentric planetary gyre. The details of their time-varying flows were sensitive to the choice of flow parametrization and amount of applied temporal regularization. Nonetheless, they reported that their time-dependent flow had a larger component of poloidal energy than the time-averaged flow, as is the case for the largest-scale part of our flows (see Fig. 7). They also highlighted an anticlockwise rotation of the time-varying flow under the equatorial western Pacific between 1997 and 2015. In CoreFlo-LL.1 we find that the azimuthal flow direction changes from westward to eastward around 2012 in this region, and this occurs via an anticlockwise rotation. Whaler et al. (2016) also found peaks of poloidal flow acceleration in the equatorial region under the Western Atlantic and the eastern Americas with opposite signs in 2006 and 2009; we find similar features in our flows, although the accelerations are generally stronger. They discussed a sign change in their non-zonal azimuthal flow at 85°W around 2006; we find the same type of change in our time-varying azimuthal flow from eastwards to westwards at this location in 2007 (see Fig. 9). To summarize, the azimuthal flow variations found by Whaler et al. (2016) are generally consistent with our results, but our flows have more energy at smaller scales, and the amplitude of our time-varying flows are larger. Given the restriction of our observations to mid and low latitudes, and because we base our flow parametrization on inertial modes in a full sphere that are not suitable for studying flow close to and inside the tangent cylinder, we are unable to comment on the acceleration of high-latitude flows recently discussed by Livermore et al. (2017) and Barrois et al. (2018b). The data employed, the choice of flow parametrization and the applied regularization all affect the details of the inferred time-varying core flow. Nevertheless, there now seems to be some agreement across the above studies that non-zonal azimuthal flow acceleration at low latitudes is required by the observations. Localized peaks in the azimuthal flow acceleration, associated with flow acceleration convergence/divergence, give rise to SA pulses with alternating polarity. The meridional part of the time-varying low-latitude flow is less well constrained by the data alone. Purely ES time-varying flow can account for the secular variation observations, but allowing cross-equatorial flow enables a similar fit to the data with a less energetic, larger scale, flow. 6.4 Possible implications for core dynamics CoreFlo-LL.1 supports the hypothesis that the equatorial region beneath the core–mantle boundary is a place of vigorous localized fluid motions, with oscillations on interannual timescales involving strong acceleration of the non-zonal azimuthal flow. Chulliat et al. (2015) proposed that these may be related to waves that could propagate in a stratified layer beneath the core–mantle boundary at low latitudes. Knezek & Buffett (2018) have developed detailed numerical models of how non-zonal MAC waves can be trapped near the equator, for plausible configurations of the background magnetic field. Such waves propagate in the azimuthal direction, and Knezek & Buffett (2018) find that they travel predominantly eastward for plausible field strengths. We do not find compelling evidence in our flow models for eastward propagation of waves close to the core surface. Rather we find strong oscillation in the non-zonal flow at particular locations, and perhaps tentative indications of westward propagation of disturbances under the Atlantic sector, see Fig. 13 as previously noted by Chulliat et al. (2015). An alternative source for low-latitude flow accelerations could be the arrival at the core–mantle boundary of QG Alfvén waves triggered by buoyancy release deep in the core (Aubert 2018). The energy of these waves is focused by the geometry of the magnetic field within the core, and their arrival at the core surface is characterized by non-zonal azimuthal flow accelerations of alternating sign, with a sub-decadal timescale related to the period of the underlying Alfvén waves. This scenario seems to be in good agreement with many of the details of our time-varying flows. The apparent westward motion of azimuthal flow acceleration features under the Atlantic may in this case perhaps be a signature of a QG Alfvén wave front arriving at the core surface first under western Africa and then progressively later at locations to the west. Aubert & Finlay (2019) have argued that such waves may provide an explanation for geomagnetic jerks. As in our flows, jerk events are then associated with a rapid change in the sign of the non-zonal azimuthal flow acceleration pattern at the core surface. 7 CONCLUSIONS We have derived a new model describing the time-varying flow at low latitudes in the Earth’s outer core, called CoreFlo-LL.1. This model, together with related movies, is available online at http://www.spacecenter.dk/files/magnetic-models/COREFLO-LL/. It is based on ground and satellite magnetic observations and assumes that the time-varying flow can be parametrized primarily in terms of geostrophic and QG modes. We find evidence for strong non-zonal azimuthal flow accelerations in the equatorial region that change sign on an interannual timescale. SA pulses at the core surface are associated with localized extrema of these azimuthal flow accelerations. We find that geomagnetic jerks may occur when the non-zonal azimuthal flow acceleration rapidly changes sign, for example, if the time-varying azimuthal flow reaches a peak in its amplitude. The 2014 jerk, that was seen most clearly in the Pacific region, involved a large-scale non-zonal azimuthal flow acceleration structure that changed sign in this fashion. One explanation for such an event could be a QG Alfvén wave, triggered deep within the core, arriving at the core surface (Aubert & Finlay 2019). We find that flows based on geostrophic and inertial mode spatial structures are able to provide an adequate description of geomagnetic observations, with time-averaged flows dominated by a planetary scale eccentric gyre. Compared to previous studies, CoreFlo-LL.1 shows more intense time-varying flows, and small lengthscale structures, localized in the equatorial region. It is remarkable that field accelerations at both mid and low latitudes at Earth’s surface, and at satellite altitude, can be well fit by time-varying flows focused in the equatorial region close to the core surface. The meridional component of our time-varying flows are however dependent on a priori assumptions made concerning the flow structure, in particular whether or not time-dependent cross-equatorial flow is permissible. In this work we have neglected the role of magnetic diffusion. Though not expected to drive field acceleration, it should be included in future studies, for example, following the approach of Amit & Christensen (2008). Moving forward it will also be important to use more complete data error covariance descriptions, for example, taking into account the correlations between errors on neighbouring ground and virtual observatories. A fundamental limitation of this study is that it uses the free oscillation modes of a rotating sphere, rather than a spherical shell, motivated by our focus on low latitudes and the existence of analytic expression for the full sphere modes. One could also, in principle, solve a numerical eigenvalue problem and obtain correct modes for a spherical shell and take into account finite dissipation; these could then instead be used as the basis for an inversion. Hypotheses testing, for example, exploring different choices of the background density profile and core magnetic field structure could then be undertaken. A more complete understanding of the physical processes underlying time-varying low-latitude core flows will ultimately require fitting of observations by full dynamic models, rather than the kinematic descriptions we have employed epoch-by-epoch. As an initial step, such models should describe non-zonal QG Alfvén waves, as these appear to be of interest in explaining SA pulses and geomagnetic jerks. The lengthening high-resolution data series being generated by the Swarm satellite mission will provide increasingly powerful observational tests of such ideas. ACKNOWLEDGEMENTS We thank ESA for providing L1b magnetic field data from the Swarm satellite mission that were crucial to this study. The staff of the geomagnetic observatories and INTERMAGNET are thanked for supplying high-quality observatory data, and BGS are thanked for providing us with checked and corrected observatory hourly mean values. The support of the CHAMP mission by the German Aerospace Center (DLR) and the Federal Ministry of Education and Research is gratefully acknowledged. Nicolas Gillet, Olivier Barrois and Julien Aubert are thanked for helpful discussions and suggestions. CK and CCF were partly supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 772561). Ciaran Beggan and an anonymous reviewer are thanked for constructive comments that improved the manuscript. Footnotes 1 ftp://ftp.nerc-murchison.ac.uk/geomag/Swarm/AUX_OBS 2 http://omniweb.gsfc.nasa.gov 3 https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html 4 The normalization is different here from the one given in Zhang & Liao (2017), where +3 is +1 in the first denominator with the double factorial. We believe this is a misprint and should be +3 as given here. REFERENCES Amit H. , Christensen U.R., 2008 . Accounting for magnetic diffusion in core flow inversions from geomagnetic secular variation , Geophys. J. Int. , 175 ( 3 ), 913 – 924 . Google Scholar Crossref Search ADS WorldCat Amit H. , Pais M.A., 2013 . Differences between tangential geostrophy and columnar flow , Geophys. J. Int. , 194 ( 1 ), 145 – 157 . Google Scholar Crossref Search ADS WorldCat Amit H. , Olson P., Christensen U., 2007 . Tests of core flow imaging methods with numerical dynamos , Geophys. J. Int. , 168 ( 1 ), 27 – 39 . Google Scholar Crossref Search ADS WorldCat Aubert J. , 2012 . Flow throughout the earth’s core inverted from geomagnetic observations and numerical dynamo models , Geophys. J. Int. , 192 ( 2 ), 537 – 556 . Google Scholar Crossref Search ADS WorldCat Aubert J. , 2018 . Geomagnetic acceleration and rapid hydromagnetic wave dynamics in advanced numerical simulations of the geodynamo , Geophys. J. Int. , 214 , 531 – 547 . Google Scholar Crossref Search ADS WorldCat Aubert J. , Finlay C.C., 2019 . Geomagnetic jerks and rapid hydromagnetic waves focusing at earth’s core surface , Nature Geosci. , (under review) . Google Scholar OpenURL Placeholder Text WorldCat Aubert J. , Gastine T., Fournier A., 2017 . Spherical convective dynamos in the rapidly rotating asymptotic regime , J. Fluid Mech. , 813 , 558 – 593 . Google Scholar Crossref Search ADS WorldCat Backus G. , Rieutord M., 2017 . Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid , Phys. Rev. E. , 95 ( 5 ), 053116 , doi:10.1103/PhysRevE.95.053116 . Google Scholar Crossref Search ADS PubMed WorldCat Backus G.E. , 1968 . Kinematics of geomagnetic secular variation in a perfectly conducting core , Phil. Trans. R. Soc. Lond. A , 263 ( 1141 ), 239 – 266 . Google Scholar Crossref Search ADS WorldCat Baerenzung J. , Holschneider M., Lesur V., 2014 . Bayesian inversion for the filtered flow at the earth’s core-mantle boundary , J. geophys. Res. , 119 ( 4 ), 2695 – 2720 . Google Scholar Crossref Search ADS WorldCat Baerenzung J. , Holschneider M., Lesur V., 2016 . The flow at the earth’s core-mantle boundary under weak prior constraints , J. geophys. Res. , 121 ( 3 ), 1343 – 1364 . Google Scholar Crossref Search ADS WorldCat Baerenzung J. , , Holschneider M., Wicht J., Sanchez S., Lesur V., 2018 . Modeling and predicting the short term evolution of the geomagnetic field , J. geophys Res. , 123 , 4539 – 4560 . Google Scholar Crossref Search ADS WorldCat Bardsley O.P. , Davidson P.A., 2016 . Inertial-alfvén waves as columnar helices in planetary cores , J. Fluid Mech. , 805 , R2 , doi:10.1017/jfm.2016.577 . Google Scholar Crossref Search ADS WorldCat Barrois O. , Gillet N., Aubert J., 2017 . Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics , Geophys. J. Int. , 211 ( 1 ), 50 – 68 . Google Scholar Crossref Search ADS WorldCat Barrois O. , Gillet N., Aubert J., Hammer M.D., Finlay C.C., Martin Y., 2018a . Erratum: ‘contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics’ and ‘assimilation of ground and satellite magnetic measurements: inference of core surface magnetic and velocity field changes’ , Geophys. J. Int. , doi:10.1093/gji/ggy47. Google Scholar OpenURL Placeholder Text WorldCat Barrois O. , Hammer M.D., Finlay C.C., Martin Y., Gillet N., 2018b . Assimilation of ground and satellite magnetic measurements: inference of core surface magnetic and velocity field changes , Geophys. J. Int. , 215 , 695 – 712 . Google Scholar Crossref Search ADS WorldCat Bloxham J. , 1992 . The steady part of the secular variation of the earth’s magnetic field , J. geophys. Res. , 97 , 19 565 – 19 579 . Google Scholar Crossref Search ADS WorldCat Bloxham J. , Jackson A., 1991 . Fluid flow near the surface of earth’s outer core , Rev. Geophys. , 29 ( 1 ), 97 – 120 . Google Scholar Crossref Search ADS WorldCat Bloxham J. , Zatman S., Dumberry M., 2002 . The origin of geomagnetic jerks , Nature , 420 ( 6911 ), 65 – 68 . Google Scholar Crossref Search ADS PubMed WorldCat Brown W. , Mound J., Livermore P., 2013 . Jerks abound: an analysis of geomagnetic observatory data from 1957 to 2008 , Phys. Earth planet. Inter. , 223 , 62 – 76 . Google Scholar Crossref Search ADS WorldCat Busse F.H. , 1970 . Thermal instabilities in rapidly rotating systems , J. Fluid Mech. , 44 ( 3 ), 441 – 460 . Google Scholar Crossref Search ADS WorldCat Busse F.H. , Zhang K., Liao X., 2005 . On slow inertial waves in the solar convection zone , Astrophys. J. , 631 , L171 – L174 . Google Scholar Crossref Search ADS WorldCat Chulliat A. , Maus S., 2014 . Geomagnetic secular acceleration, jerks, and a localized standing wave at the core surface from 2000 to 2010 , J. geophys. Res. , 119 ( 3 ), 1531 – 1543 . Google Scholar Crossref Search ADS WorldCat Chulliat A. , Thebault E., Hulot G., 2010 . Core field acceleration pulse as a common cause of the 2003 and 2007 geomagnetic jerks , Geophys. Res. Lett. , 37 ( 7 ), doi:10.1029/2009GL042019 . Google Scholar OpenURL Placeholder Text WorldCat Chulliat A. , Alken P., Maus S., 2015 . Fast equatorial waves propagating at the top of the earth’s core , Geophys. Res. Lett. , 42 ( 9 ), 3321 – 3329 . Google Scholar Crossref Search ADS WorldCat Constable C. G. , 1988 . Parameter estimation in non-Gaussian noise , Geophys. J. Int. , 94 ( 1 ), 131 – 142 . Google Scholar Crossref Search ADS WorldCat Courtillot V. , Le Mouël J., 1984 . Geomagnetic secular variation impulses , Nature , 311 ( 5988 ), 709 – 716 . Google Scholar Crossref Search ADS WorldCat Cui Z. , Zhang K., Liao X., 2014 . On the completeness of inertial wave modes in rotating annular channels , Geophys. astrophys. Fluid. Dyn. , 108 ( 1 ), 44 – 59 . Google Scholar Crossref Search ADS WorldCat Eymin C. , Hulot G., 2005 . On core surface flows inferred from satellite magnetic data , Phys. Earth planet. Inter. , 152 ( 3 ), 200 – 220 . Google Scholar Crossref Search ADS WorldCat Farquharson C.G. , Oldenburg D.W., 1998 . Non-linear inversion using general measures of data misfit and model structure , Geophys. J. Int. , 134 ( 1 ), 213 – 227 . Google Scholar Crossref Search ADS WorldCat Finlay C.C. , Olsen N., Tøffner-Clausen L., 2015 . DTU candidate field models for IGRF-12 and the CHAOS-5 geomagnetic field model , Earth Planets Space , 67 ( 1 ), doi:10.1186/s40623-015-0274-3 . Google Scholar OpenURL Placeholder Text WorldCat Finlay C.C. , Olsen N., Kotsiaros S., Gillet N., Tøffner-Clausen L., 2016 . Recent geomagnetic secular variation from swarm and ground observatories as estimated in the CHAOS-6 geomagnetic field model , Earth Planets Space , 68 ( 1 ), doi:10.1186/s40623-016-0486-1 . Google Scholar OpenURL Placeholder Text WorldCat Gillet N. , Pais M.A., Jault D., 2009 . Ensemble inversion of time-dependent core flow models , Geochem. Geophys. Geosyst. , 10 ( 6 ), doi:10.1029/2008GC002290 . Google Scholar OpenURL Placeholder Text WorldCat Gillet N. , Schaeffer N., Jault D., 2012 . Rationale and geophysical evidence for quasi-geostrophic rapid dynamics within the earth’s outer core , Phys. Earth planet. Inter. , 202 , 78 – 88 . Google Scholar Crossref Search ADS WorldCat Gillet N. , Jault D., Finlay C.C., Olsen N., 2013 . Stochastic modeling of the earth’s magnetic field: Inversion for covariances over the observatory era , Geochem. Geophys. Geosyst. , 14 ( 4 ), 766 – 786 . Google Scholar Crossref Search ADS WorldCat Gillet N. , Barrois O., Finlay C.C., 2015a . Stochastic forecasting of the geomagnetic field from the COV-OBS.x1 geomagnetic field model, and candidate models for IGRF-12 , Earth Planets Space , 67 ( 1 ), doi:10.1186/s40623-015-0225-z . Google Scholar OpenURL Placeholder Text WorldCat Gillet N. , Jault D., Finlay C.C., 2015b . Planetary gyre, time-dependent eddies, torsional waves, and equatorial jets at the earth’s core surface , J. geophys. Res. , 120 ( 6 ), 3991 – 4013 . Google Scholar Crossref Search ADS WorldCat Gubbins D. , 1982 . Finding core motions from magnetic observations , Phil. Trans. R. Soc. Lond. A , 306 ( 1492 ), 247 – 254 . Google Scholar Crossref Search ADS WorldCat Gubbins D. , Roberts P., 1987 . Magnetohydrodynamics of the earth’s core , Geomagnetism , 2 , 1 – 183 . Google Scholar OpenURL Placeholder Text WorldCat Hide R. , 1966 . Free hydromagnetic oscillations of the earth’s core and the theory of the geomagnetic secular variation , Phil. Trans. R. Soc. Lond. A , 259 ( 1107 ), 615 – 647 . Google Scholar Crossref Search ADS WorldCat Holme R. , 2015 . Large-Scale Flow in the Core , in: Treatise on Geophysics , vol. 8 , ed. Olson P., Elsevier , 2nd edn., pp. 91 – 113 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Holme R. , Olsen N., 2006 . Core surface flow modelling from high-resolution secular variation , Geophys. J. Int. , 166 ( 2 ), 518 – 528 . Google Scholar Crossref Search ADS WorldCat Hulot G. , Mouël J.L., Wahr J., 1992 . Taking into account truncation problems and geomagnetic model accuracy in assessing computed flows at the core-mantle boundary , Geophys. J. Int. , 108 ( 1 ), 224 – 246 . Google Scholar Crossref Search ADS WorldCat Ivers D.J. , Jackson A., Winch D., 2015 . Enumeration, orthogonality and completeness of the incompressible coriolis modes in a sphere , J. Fluid Mech. , 766 , 468 – 498 . Google Scholar Crossref Search ADS WorldCat Jackson A. , 1997 . Time-dependency of tangentially geostrophic core surface motions , Phys. Earth planet. Inter. , 103 ( 3–4 ), 293 – 311 . Google Scholar Crossref Search ADS WorldCat Jackson A. , Bloxham J., Gubbins D., 1993 . Time-dependent flow at the core surface and conservation of angular momentum in the coupled core-mantle system , in Dynamics of Earth’s Deep Interior and Earth Rotation , vo. 72 , eds Le Mouël J.‐L., Smylie D.E., Herring T., AGU Publ , pp. 97 – 107 . Google Scholar OpenURL Placeholder Text WorldCat Jault D. , 2008 . Axial invariance of rapidly varying diffusionless motions in the earth’s core interior , Phys. Earth planet. Inter. , 166 ( 1–2 ), 67 – 76 . Google Scholar Crossref Search ADS WorldCat Jault D. , Gire C., Le Mouël J., 1988 . Westward drift, core motions and exchanges of angular momentum between core and mantle , Nature , 333 ( 6171 ), 353 . Google Scholar Crossref Search ADS WorldCat Kahle A.B. , Vestine E., Ball R., 1967 . Estimated surface motions of the earth’s core , J. geophys. Res. , 72 ( 3 ), 1095 – 1108 . Google Scholar Crossref Search ADS WorldCat Knezek N. , Buffett B., 2018 . Influence of magnetic field configuration on magnetohydrodynamic waves in earth’s core , Phys. Earth planet. Inter. , 277 , 1 – 9 . Google Scholar Crossref Search ADS WorldCat Le Huy M. , Mandea M., Le Mouël J.-L., Pais A., 2000 . Time evolution of the fluid flow at the top of the core. geomagnetic jerks , Earth Planets Space , 52 ( 3 ), 163 – 173 . Google Scholar Crossref Search ADS WorldCat Le Mouël J. , 1984 . Outer-core geostrophic flow and secular variation of earth’s geomagnetic field , Nature , 311 ( 5988 ), 734 – 735 . Google Scholar Crossref Search ADS WorldCat Leopardi P. , 2006 . A partition of the unit sphere into regions of equal area and small diameter , Electron. Trans. Numer. Anal. , 25 , 309 – 327 . Google Scholar OpenURL Placeholder Text WorldCat Lesur V. , Wardinski I., Asari S., Minchev B., Mandea M., 2010a . Modelling the earth’s core magnetic field under flow constraints , Earth Planets Space , 62 ( 6 ), 503 – 516 . Google Scholar Crossref Search ADS WorldCat Lesur V. , Wardinski I., Hamoudi M., Rother M., 2010b . The second generation of the gfz reference internal magnetic model: Grimm-2 , Earth Planets Space , 62 ( 10 ), 6 , doi:10.5047/eps.2010.07.007 . Google Scholar Crossref Search ADS WorldCat Lesur V. , Whaler K., Wardinski I., 2015 . Are geomagnetic data consistent with stably stratified flow at the core–mantle boundary? , Geophys. J. Int. , 201 ( 2 ), 929 – 946 . Google Scholar Crossref Search ADS WorldCat Lesur V. , Wardinski I., Baerenzung J., Holschneider M., 2018 . On the frequency spectra of the core magnetic field Gauss coefficients , Phys. Earth planet. Inter. , 276 , 145 – 158 . Google Scholar Crossref Search ADS WorldCat Liao X. , Zhang K., 2010 . Asymptotic and numerical solutions of the initial value problem in rotating planetary fluid cores , Geophys. J. Int. , 180 ( 1 ), 181 – 192 . Google Scholar Crossref Search ADS WorldCat Liao X. , Zhang K., Li L., 2012 . Asymptotic solutions of differential rotation driven by convection in rapidly rotating fluid spheres with the non-slip boundary condition , Geophys. astrophys. Fluid Dyn. , 106 , 643 – 659 . Google Scholar Crossref Search ADS WorldCat Livermore P.W. , Hollerbach R., Finlay C.C., 2017 . An accelerating high-latitude jet in earth’s core , Nature Geosci. , 10 , 62 – 68 . Google Scholar Crossref Search ADS WorldCat Lloyd D. , Gubbins D., 1990 . Toroidal fluid motion at the top of the earth’s core , Geophys. J. Int. , 100 ( 3 ), 455 – 467 . Google Scholar Crossref Search ADS WorldCat Macmillan S. , Olsen N., 2013 . Observatory data and the swarm mission , Earth Planets Space , 65 , 1355 – 1362 . Google Scholar Crossref Search ADS WorldCat Maffei S. , Jackson A., Livermore P.W., 2017 . Characterization of columnar inertial modes in rapidly rotating spheres and spheroids , Proc. R. Soc. A , 473 ( 2204 ), 20170181 , doi:10.1098/rspa.2017.0181 . Google Scholar Crossref Search ADS WorldCat Mandea M. , Olsen N., 2006 . A new approach to directly determine the secular variation from magnetic satellite observations , Geophys. Res. Lett. , 33 ( 15) , doi:10.1029/2006GL026616. Google Scholar OpenURL Placeholder Text WorldCat Mandea M. , Holme R., Pais A., Pinheiro K., Jackson A., Verbanac G., 2010 . Geomagnetic jerks: rapid core field variations and core dynamics , Space Sci. Rev. , 155 ( 1–4 ), 147 – 175 . Google Scholar Crossref Search ADS WorldCat Matzka J. , Chulliat A., Mandea M., Finlay C.C., Qamili E., 2010 . Geomagnetic observations for main field studies: from ground to space , Space Sci. Rev. , 155 , 29 – 64 . Google Scholar Crossref Search ADS WorldCat Olsen N. , Mandea M., 2007 . Investigation of a secular variation impulse using satellite data: the 2003 geomagnetic jerk , Earth planet. Sci. Lett. , 255 ( 1–2 ), 94 – 105 . Google Scholar Crossref Search ADS WorldCat Olsen N. , Mandea M., 2008 . Rapidly changing flows in the earth’s core , Nature Geosci. , 1 ( 6 ), 390 , doi:10.1038/ngeo203 . Google Scholar Crossref Search ADS WorldCat Olsen N. , Mandea M., Sabaka T.J., Tøffner-Clausen L., 2009 . Chaos-2—a geomagnetic field model derived from one decade of continuous satellite data , Geophys. J. Int. , 179 ( 3 ), 1477 – 1487 . Google Scholar Crossref Search ADS WorldCat Olsen N. , Lühr H., Finlay C.C., Sabaka T.J., Michaelis I., Rauberg J., Tøffner-Clausen L., 2014 . The chaos-4 geomagnetic field model , Geophys. J. Int. , 197 , 815 – 827 . Google Scholar Crossref Search ADS WorldCat Olsen N. et al. , 2015 . The swarm initial field model for the 2014 geomagnetic field , Geophys. Res. Lett. , 42 ( 4 ), 1092 – 1098 . Google Scholar Crossref Search ADS WorldCat Pais M. , Jault D., 2008 . Quasi-geostrophic flows responsible for the secular variation of the earth’s magnetic field , Geophys. J. Int. , 173 ( 2 ), 421 – 443 . Google Scholar Crossref Search ADS WorldCat Pais M.A. , Morozova A.L., Schaeffer N., 2014 . Variability modes in core flows inverted from geomagnetic field models , Geophys. J. Int. , 200 ( 1 ), 402 – 420 . Google Scholar Crossref Search ADS WorldCat Rasmussen C.E. , Williams C. K.I., 2006 . Gaussian Process for Machine Learning , MIT press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rau S. , Christensen U., Jackson A., Wicht J., 2000 . Core flow inversion tested with numerical dynamo models , Geophys. J. Int. , 141 ( 2 ), 485 – 497 . Google Scholar Crossref Search ADS WorldCat Roberts P. , Scott S., 1965 . On analysis of the secular variation , J. Geomag. Geoelectr. , 17 ( 2 ), 137 – 151 . Google Scholar Crossref Search ADS WorldCat Roberts P.H. , Glatzmaier G.A., 2000 . A test of the frozen-flux approximation using a new geodynamo model , Phil. Trans. R. Soc. Lond., A. , 358 ( 1768 ), 1109 – 1121 . Google Scholar Crossref Search ADS WorldCat Sabaka T.J. , Olsen N., Purucker M.E., 2004 . Extending comprehensive models of the earth’s magnetic field with ørsted and champ data , Geophys. J. Int. , 159 ( 2 ), 521 – 547 . Google Scholar Crossref Search ADS WorldCat Sabaka T.J. , Olsen N., Tyler R.H., Kuvshinov A., 2015 . Cm5, a pre-swarm comprehensive geomagnetic field model derived from over 12 yr of champ, ørsted, sac-c and observatory data , Geophys. J. Int. , 200 ( 3 ), 1596 – 1626 . Google Scholar Crossref Search ADS WorldCat Schaeffer N. , Cardin P., 2006 . Quasi-geostrophic kinematic dynamos at low magnetic Prandtl number , Earth planet. Sci. Lett. , 245 ( 3–4 ), 595 – 604 . Google Scholar Crossref Search ADS WorldCat Schaeffer N. , Pais M.A., 2011 . On symmetry and anisotropy of earth-core flows , Geophys. Res. Lett. , 38 ( 10 , doi:10.1029/2011GL046888) , Google Scholar OpenURL Placeholder Text WorldCat Schaeffer N. , Jault D., Nataf H.-C., Fournier A., 2017 . Turbulent geodynamo simulations: a leap towards earth’s core , Geophys. J. Int. , 211 ( 1 ), 1 – 29 . Google Scholar Crossref Search ADS WorldCat Silva L. , Hulot G., 2012 . Investigating the 2003 geomagnetic jerk by simultaneous inversion of the secular variation and acceleration for both the core flow and its acceleration , Phys. Earth planet. Inter. , 198 , 28 – 50 . Google Scholar Crossref Search ADS WorldCat Torta J.M. , Pavón-Carrasco F.J., Marsal S., Finlay C.C., 2015 . Evidence for a new geomagnetic jerk in 2014 , Geophys. Res. Lett. , 42 ( 19 ), 7933 – 7940 . Google Scholar Crossref Search ADS WorldCat Wardinski I. , Lesur V., 2012 . An extended version of the c3fm geomagnetic field model: application of a continuous frozen-flux constraint , Geophys. J. Int. , 189 ( 3 ), 1409 – 1429 . Google Scholar Crossref Search ADS WorldCat Wardinski I. , Holme R., Asari S., Mandea M., 2008 . The 2003 geomagnetic jerk and its relation to the core surface flows , Earth planet. Sci. Lett. , 267 ( 3–4 ), 468 – 481 . Google Scholar Crossref Search ADS WorldCat Whaler K. , 1986 . Geomagnetic evidence for fluid upwelling at the core-mantle boundary , Geophys. J. Int. , 86 ( 2 ), 563 – 588 . Google Scholar Crossref Search ADS WorldCat Whaler K. , Beggan C., 2015 . Derivation and use of core surface flows for forecasting secular variation , J. geophys. Res. , 120 ( 3 ), 1400 – 1414 . Google Scholar Crossref Search ADS WorldCat Whaler K.A. , Olsen N., Finlay C.C., 2016 . Decadal variability in core surface flows deduced from geomagnetic observatory monthly means , Geophys. J. Int. , 207 ( 1 ), 228 – 243 . Google Scholar Crossref Search ADS WorldCat Zhang K. , 1993 . On equatorially trapped boundary inertial waves , J. Fluid Mech. , 248 , 203 – 217 . Google Scholar Crossref Search ADS WorldCat Zhang K. , Liao X., 2004 . A new asymptotic method for the analysis of convection in a rapidly rotating sphere , J. Fluid Mech. , 518 , 319 – 346 . Google Scholar Crossref Search ADS WorldCat Zhang K. , Liao X., 2017 . Theory and Modeling of Rotating Fluids , Cambridge monographs on mechanics , Cambridge Univ. Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Zhang K. , Earnshaw P., Liao X., Busse F., 2001 . On inertial waves in a rotating fluid sphere , J. Fluid Mech. , 437 , 103 – 119 . Google Scholar Crossref Search ADS WorldCat Zhang K. , Liao X., Busse F., 2007 . Asymptotic solutions of convection in rapidly rotating non-slip spheres , J. Fluid Mech. , 578 , 371 – 380 . Google Scholar Crossref Search ADS WorldCat APPENDIX A: EXPLICIT EXPRESSIONS OF MODES A1 Geostrophic mode The unnormalized geostrophic polynomials with k ≥ 1 are defined by (Zhang & Liao 2017) \begin{eqnarray*} G_{2k-1}(s)=\sum _{j=1}^{k} C^\mathrm{G}_{k;j}s^{2j-1}, \end{eqnarray*} (A1) with s = rsin θ and \begin{eqnarray*} C^\mathrm{G}_{k;j} = \frac{(-1)^{k-j}[2(k+j)-1]!!}{2^{k-1}(k-j)!(j-1)!(2j)!!}. \end{eqnarray*} (A2) In this case, the mean over the full sphere of the squared absolute value becomes \begin{eqnarray*} \frac{3}{4\pi }\int _\mathcal {V}\left|G_{2k-1}\right|^2\mathrm{d}\mathcal {V}=\frac{3(2k+1)!!(2k-1)!!}{(4k+1)(2k)!!(2k-2)!!}, \end{eqnarray*} (A3) which is used to normalize the polynomials. The enstrophy as defined in eq. (33) of the unnormalized geostrophic polynomials is \begin{eqnarray*} Q^2\left( G_{2k-1}{\hat{\boldsymbol \Phi }}\right) =3\sum _{jl}^kC^\mathrm{G}_{k;j}C^\mathrm{G}_{k;l}2^{j+l}jl\frac{[j+l-2]!}{[2(j+l)-1]!!}. \end{eqnarray*} (A4) A2 Symmetric inertial modes Following closely Zhang et al. (2001) and Zhang & Liao (2017), the unnormalized and complex-valued three components of both non-axisymmetric and axisymmetric modes are given in spherical coordinates by \begin{eqnarray*} \begin{split} \mathbf {\hat{r}}\cdot \mathbf {u}_{mnk}^\mathrm{S} = -\frac{\mathrm{i}}{2}\sum _{i=0}^{k}\sum _{j=0}^{k-i}C^\mathrm{S}_{mk;ij}\sigma ^{2i-1}(1-\sigma ^2)^{j-1} r^{m+2(i+j)-1}\sin ^{m+2j}\theta \cos ^{2i}\theta \mathrm{e}^{\mathrm{i}m\phi } [\sigma (m + m\sigma + 2j\sigma )-2i(1-\sigma ^2)] \end{split} \end{eqnarray*} (A5a) \begin{eqnarray*} \begin{split} {\hat{\boldsymbol \theta }}\cdot \mathbf {u}_{mnk}^\mathrm{S}= -\frac{\mathrm{i}}{2}\sum _{i=0}^{k}\sum _{j=0}^{k-i}C^\mathrm{S}_{mk;ij}\sigma ^{2i-1}(1-\sigma ^2)^{j-1} r^{m+2(i+j)-1}\sin ^{m+2j-1}\theta \cos ^{2i-1}\theta \mathrm{e}^{\mathrm{i}m\phi } [\sigma (m + m\sigma + 2j\sigma )\cos ^2\theta +2i(1-\sigma ^2)\sin ^2\theta ] \end{split} \end{eqnarray*} (A5b) \begin{eqnarray*} \begin{split} {\hat{\boldsymbol \Phi }}\cdot \mathbf {u}_{mnk}^\mathrm{S}= \frac{1}{2}\sum _{i=0}^{k}\sum _{j=0}^{k-i} C^\mathrm{S}_{mk;ij}\sigma ^{2i}(1-\sigma ^2)^{j-1} r^{m+2(i+j)-1}\sin ^{m+2j-1}\theta \cos ^{2i}\theta \mathrm{e}^{\mathrm{i}m\phi } (m + m\sigma + 2j), \end{split} \end{eqnarray*} (A5c) where the coefficients |$C^\mathrm{S}_{mk;ij}$| are defined as \begin{eqnarray*} C^\mathrm{S}_{mk;ij}=\frac{(-1)^{i+j}[2(k+m+i+j)-1]!!}{2^{j+1}(2i-1)!!(k-i-j)!i!j!(m+j)!}. \end{eqnarray*} (A6) The subscripts of σmnk have been omitted to make the equations easier to read. The normalization of the modes is given by \begin{eqnarray*} \begin{aligned} \frac{3}{4\pi }\int _{\mathcal {V}}\left|\mathbf {u}^\mathrm{S}_{mnk}\right|^2\mathrm{d}\mathcal {V}=\sum _{i=0}^{k}\sum _{j=0}^{k-i}\sum _{q=0}^{k}\sum _{l=0}^{k-q}&C^\mathrm{S}_{mk;ij}C^\mathrm{S}_{mk;ql} \frac{3\cdot 2^{m+j+l-3}}{[2(m+i+j+q+l)+1]!!}\sigma ^{2(i+q)}(1-\sigma ^2)^{j+l}\\ &\times \bigg (\big [(m+m\sigma +2j)(m+m\sigma +2l) +(m+m\sigma +2j\sigma )(m+m\sigma +2l\sigma )\big ]\\ &\phantom{\bigg (}\times (m+j+l-1)!\frac{[2(i+q)-1]!!}{(1-\sigma ^2)^2} +8iq(m+j+l)!\frac{[2(i+q)-3]!!}{\sigma ^2}\bigg ). \end{aligned} \end{eqnarray*} (A7) Given the components in eq. (A5), one can also compute the inner product of the vorticity of two unnormalized ES modes |$\mathbf {u}_\alpha ^\mathrm{S}$| and |$\mathbf {u}_\beta ^\mathrm{S}$| with \begin{eqnarray*} \begin{aligned} \frac{3}{4\pi }\int _\mathcal {V}(\nabla \times \mathbf {u}_\alpha ^\mathrm{S})\cdot (\nabla &\times \mathbf {u}_\beta ^\mathrm{S})^*\mathrm{d}\mathcal {V}=\\ \sum _{i=0}^{k_\alpha }\sum _{j=0}^{k_\alpha -i}\sum _{q=0}^{k_\beta }\sum _{l=0}^{k_\beta -q}& C^\mathrm{S}_{mk_\alpha ;ij} C^\mathrm{S}_{mk_\beta ;ql}\frac{3\cdot 2^{j+l+m-1}iq}{[2(i+j+q+l+m)-1]!!} \sigma _\alpha ^{2i-2}\sigma _\beta ^{2q-2}(1-\sigma _\alpha ^2)^{j-1}(1-\sigma _\beta ^2 )^{l-1}\\ &\times \bigg (\sigma _\alpha \sigma _\beta [l+j+m-1]![2(i+q)-3]!\\ &\phantom{\bigg (}\times \big [(2j\sigma _\alpha +m+m\sigma _\alpha )(2l\sigma _\beta +m+m\sigma _\beta )+(2j+m+m\sigma _\alpha )(2l+m+m\sigma _\beta )\big ]\\ &\phantom{\times \bigg (}+2(2i-1)(2q-1)(1-\sigma _\alpha ^2)(1-\sigma _\beta ^2)[l+j+m]![2(i+q)-5]!!\bigg ), \end{aligned} \end{eqnarray*} (A8) for mα = mβ = m and zero otherwise. The asterisk denotes the complex conjugate. In this paper, we only used terms with α = β. A3 Antisymmetric inertial modes Again following closely Zhang et al. (2001) and Zhang & Liao (2017) the three unnormalized and complex-valued velocity components in spherical polar coordinates are \begin{eqnarray*} \begin{split} \mathbf {\hat{r}}\cdot \mathbf {u}^\mathrm{A}_{mnk}=-\frac{\mathrm{i}}{2}\sum _{i=0}^{k}\sum _{j=0}^{k-i}C^\mathrm{A}_{mk;ij}\sigma ^{2i-1}(1-\sigma ^2)^{j-1} r^{m+2(i+j)}\sin ^{m+2j}\theta \cos ^{2i+1}\theta \mathrm{e}^{\mathrm{i}m\phi } [\sigma (m + m\sigma + 2j\sigma )-(2i+1)(1-\sigma ^2)] \end{split} \end{eqnarray*} (A9a) \begin{eqnarray*} \begin{split} {\hat{\boldsymbol \theta}}\cdot {\mathbf u}^{\mathrm A}_{mnk}=-\frac{\mathrm{i}}{2}\sum_{i=0}^{k}\sum_{j=0}^{k-i} C^{\mathrm A}_{mk;ij}\sigma^{2i-1}(1-\sigma^2)^{j-1} r^{m+2(i+j)}\sin^{m+2j-1}\theta \cos^{2i}\theta {\mathrm e}^{\mathrm{i}m\phi}\\ \qquad \times [\sigma (m + m\sigma + 2j\sigma )\cos^2\theta +(2i+1)(1-\sigma^2)\sin^2\theta ] \end{split} \end{eqnarray*} (A9b) \begin{eqnarray*} \begin{split} {\hat{\boldsymbol \Phi }}\cdot \mathbf {u}^\mathrm{A}_{mnk}=\frac{1}{2}\sum _{i=0}^{k}\sum _{j=0}^{k-i}C^\mathrm{A}_{mk;ij}\sigma ^{2i}(1-\sigma ^2)^{j-1} r^{m+2(i+j)}\sin ^{m+2j-1}\theta \cos ^{2i+1}\theta \mathrm{e}^{\mathrm{i}m\phi } (m + m\sigma + 2j), \end{split} \end{eqnarray*} (A9c) where the coefficients |$C^\mathrm{A}_{mk;ij}$| are defined as \begin{eqnarray*} C^\mathrm{A}_{ij;km}=\frac{(-1)^{i+j}[2(k+m+i+j)+1]!!}{2^{j+1}(2i+1)!!(k-i-j)!i!j!(m+j)!}. \end{eqnarray*} (A10) The normalization4 of the components in this definition is \begin{eqnarray*} \begin{aligned} \frac{3}{4\pi }\int _{\mathcal {V}}\left|\mathbf {u}^\mathrm{A}_{mnk}\right|^2\mathrm{d}\mathcal {V}=\sum _{i=0}^{k}\sum _{j=0}^{k-i}\sum _{q=0}^{k}\sum _{l=0}^{k-q}&C^\mathrm{A}_{mk;ij}C^\mathrm{A}_{mk;ql} \frac{3\cdot 2^{m+j+l-3}}{[2(m+i+j+q+l)+3]!!}\sigma ^{2(i+q)}(1-\sigma ^2)^{j+l}\\ &\times \bigg ([(m+m\sigma +2j)(m+m\sigma +2l)+(m+m\sigma +2j\sigma )(m+m\sigma +2l\sigma )]\\ &\phantom{\bigg (}\times (m+j+l-1)!\frac{[2(i+q)+1]!!}{(1-\sigma ^2)^2}+2(2i+1)(2q+1)(m+j+l)!\\ &\phantom{\bigg (}\times \frac{[2(i+q)-1]!!}{\sigma ^2}\bigg ). \end{aligned} \end{eqnarray*} (A11) The inner product of the vorticity of two unnormalized modes becomes \begin{eqnarray*} \begin{aligned} \frac{3}{4\pi }\int _\mathcal {V}(\nabla \times \mathbf {u}_\alpha ^\mathrm{A})\cdot (\nabla &\times \mathbf {u}_\beta ^\mathrm{A})^*\mathrm{d}\mathcal {V}=\\ \sum _{i=0}^{k_\alpha }\sum _{j=0}^{k_\alpha -i}\sum _{q=0}^{k_\beta }\sum _{l=0}^{k_\beta -q}& C^\mathrm{A}_{mk_\alpha ;ij} C^\mathrm{A}_{mk_\beta ;ql}\frac{3\cdot 2^{j+l+m-3}(2i+1)(2q+1)}{[2(i+j+q+l+m)+1]!!} \sigma _\alpha ^{2i-2}\sigma _\beta ^{2q-2}(1-\sigma _\alpha ^2)^{j-1}(1-\sigma _\beta ^2 )^{l-1}\\ &\times \bigg (\sigma _\alpha \sigma _\beta [l+j+m-1]![2(i+q)-1]!\bigg .\\ &\phantom{\bigg (}\times \big [(2j\sigma _\alpha +m+m\sigma _\alpha )(2l\sigma _\beta +m+m\sigma _\beta )+(2j+m+m\sigma _\alpha )(2l+m+m\sigma _\beta )\big ]\\ &\phantom{\times \bigg (}+8iq(1-\sigma _\alpha ^2)(1-\sigma _\beta ^2)[l+j+m]![2(i+q)-3]!!\bigg ), \end{aligned} \end{eqnarray*} (A12) for mα = mβ = m and zero otherwise. In this paper, we only used terms with α = β. APPENDIX B: DETAILS ON SMALL-SCALE FIELD COVARIANCE For |$\mathbf {C}^\mathrm{\vphantom{T}}_{\tilde{b}}$|⁠, we assume at all epochs that the expansion of the unresolved small-scale core field |$\mathbf {\tilde{b}}_p$| has zero mean and variance |$\sigma ^2_{\tilde{b}}$|⁠, and are correlated in time according to a correlation function ρ that only depends on the SH degree n. We estimate the variance by extending the spectrum of the CHAOS-6-x7 model at the surface up to degree |$N_{\tilde{b}} = 30$| in the year 2015 \begin{eqnarray*} W_b(n)=(n+1)\sum _{m=1}^n\big (g_n^m(t)^2+h_n^m(t)^2\big )_{t=2015} \end{eqnarray*} (B1) with an exponential fit to the spectrum for degrees n ∈ [2, 13]. The variance of the SH expansion, which we assume to be independent of order m and suitable for all epochs, is then \begin{eqnarray*} \sigma ^2_{\tilde{b}}(n)=\frac{1.10\cdot 10^9}{(n+1)(2n+1)}\exp (-1.26n), \end{eqnarray*} (B2) with |$N_b\lt n\le N_{\tilde{b}}$|⁠. The correlation time in years can also be derived from the CHAOS-6-x7 model via (Gillet et al.2013) \begin{eqnarray*} \tau _\mathrm{c}(n)=\sqrt{3\frac{W_b(n)}{W_{\dot{b}}(n)}}, \end{eqnarray*} (B3) where |$W_{\dot{b}}$| is the power spectrum of the SV and given by the expression in eq. (B1) by replacing the SH coefficients with the respective time derivatives |$\dot{g}_n^m$| and |$\dot{h}_n^m$|⁠. An extension towards the degrees |$N_b\lt n\le N_{\tilde{b}}$| is achieved by fitting an exponential model for degrees n ∈ [2, 10] \begin{eqnarray*} \tau _\mathrm{c}(n) = 5.00\cdot 10^2\exp (-0.23n). \end{eqnarray*} (B4) © The Author(s) 2019. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © The Author(s) 2019. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Time-dependent low-latitude core flow and geomagnetic field acceleration pulses
JO - Geophysical Journal International
DO - 10.1093/gji/ggy545
DA - 2019-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/time-dependent-low-latitude-core-flow-and-geomagnetic-field-TD9BWwj26e
SP - 140
EP - 168
VL - 217
IS - 1
DP - DeepDyve
ER -