Time-dependent low latitude core flow and geomagnetic field acceleration pulses

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 {\it Swarm} and CHAMP satellites and at ground magnetic observatories. The model, called {\it 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. We find that the equatorial region beneath the core-mantle boundary is a place of vigorous, localised, 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 acceleration sign changes can occur within a year or less, and when the structures involved are of large spatial scale they can give rise to geomagnetic jerks at the Earth's surface.


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;Brown et al. 2013;Mandea et al. 2010). 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;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 & Maus 2014;Torta et al. 2015).
An important step towards understanding the underlying physical mechanism producing this 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 years, 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 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% 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 they possessed significant power in periods between 4 and 9.5 years, it was argued that they 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. (2017Barrois et al. ( , 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 equatorially region. In a follow up study Barrois et al. (2018b,a), applying the same technique to ground and satellite SV data from the past twenty years, found low latitude flow accelerations on interannual time scales, 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 present 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, that 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 they 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 parameterize 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 parameterization 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 secular variation 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 secular variation (SV) and secular acceleration (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 localised 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 which occur when there is an abrupt sign change or cross-over 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.

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 , version 0113, that was built from INTER-MAGNET and WDC Edinburgh data as available in January 2018. We use data spanning the interval 2000-2018. Revised monthly means are aimed 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 from the CHAOS-6-x7 field model ) (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 times series of secular variation (SV) at each observatory location.
In this study, we choose to use only data collected at non-polar latitudes, equatorward 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 dataset 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; and GUA -Guam, USA, Western Pacific) are shown by the black dots in Fig. 2.
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 de-trending. 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.

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 between July 2000 and September 2010 and by the Swarm vector field magnetometers, on board all three satellites (Alpha, Bravo, Charlie), between January 2014 and July 2018. Starting from the CHAMP MAG-L3 and Swarm Level 1b MAG-L, version 0501, data products, we sub-sampled the data at 15s intervals using the data in the vector field magnetometer (VFM) frame and rotating them into an Earth-Centered Earth-Fixed (ECEF) 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 o ; (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 Em = 0.33v 4/3 B 2/3 t sin(|Θ|/2), v the solar wind speed, Θ = arctan(By/Bz) and Bt = B 2 y + B 2 z with By and Bz the components of the interplanetary magnetic field (IMF), calculated using 2 hourly means of 1-min values of the IMF and solar wind extracted from the OMNI database † ; (v) IMF Bz > 0 nT and IMF |By| < 10 nT, again based on 2 hourly mean of 1 minute 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 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 was then calculated from the potential, adding this back on to the prediction of CHAOS-6-x7 at the target location. Note 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.  , the CHAOS-6-x7 SV prediction (red), and the SV prediction from the CoreFlo-LL.1 model (green). Note that there is data gap in the time series of the virtual observatories due to lack of satellite data during that time.
Rather than working directly with the satellite vector field data, we use sums and differences of the vector data, along track (using 15 second 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). 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.
Times series of such 4-monthly 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, 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 smaller than for CHAMP. We again chose to use only data series from latitudes equatorward 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 times 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 de-trending), 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.

Flow in a rapidly rotating sphere
In this section, we present details of our parameterization 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 Ω (e.g. Gubbins & Roberts 1987 where F = (ν∇ 2 u + ρ g + 1 ρ 0 J × B) is a forcing/dissipation term, ρ is the relative departure from the reference density ρ0, g is the acceleration due to gravity, ν is the coefficient of kinematic viscosity, J is the current density and B is the magnetic field. P is the departure from the pressure of the reference state which 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 F = 0, and further neglect the non-linear term u · ∇u by assuming that the flow departs only slightly from rigid body rotation. For a homogeneous and inviscid fluid confined to a spherical container of radius r0 that is uniformly rotating with Ω = Ωẑ, whereẑ 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 while the required boundary condition is that the flow normal to the boundary (in directionr) vanishes at the spherical container at r = 1, Eqns.
(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 be then 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.

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- Applying the curl leads to the Taylor-Proudman theorem, which states that 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 so-called 'geostrophic' polynomials G 2k−1 for integers k ≥ 1 (Zhang & Liao 2017) of the form involving unknown coefficients a G k , the radius r, the colatitude θ and the unit vectorφ in azimuth (see Appendix. A1 for details of G 2k−1 ). Although the geostrophic polynomials depend on the two coordinates r and θ, they take a single argument s ≡ r sin θ, of which they are odd functions. They are also orthogonal and can be normalized with respect to their mean value over the sphere, i.e. given a second geostrophic polynomial G 2l−1 for an integer l, they satisfy where δ kl is the kronecker delta and V dV denotes the integration with φ being the longitude.

Inertial modes
Allowing for time-dependent flows, the full solution to Eq. (2) is obtained using the superposition of an infinite number of 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 r = (r, θ, φ) as where i is the imaginary number and σ = ω 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 (E S ) inertial modes satisfy whereas the equatorially antisymmetric (E A ) inertial modes obey 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 u mnk , are orthogonal when integrated over the full sphere and can therefore be normalized such that 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 E S and E A modes separately in the following.

Equatorially Symmetric Modes
The E S inertial modes u S mnk , can be divided into axisymmetric (m = 0) and non-axisymmetric (m ≥ 1) modes. The half-frequencies σ of the axisymmetric E S modes are the roots of the polynomial (Zhang & Liao 2017 with k ≥ 2. Hence, for a given k, there are k − 1 axisymmetric E S modes corresponding to the k − 1 distinct roots. In case of the nonaxisymmetric E S modes, the half-frequencies satisfy 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 E S inertial mode for a given m and k (see Appendix A2 for explicit expressions).

Quasi-Geostrophic Modes
A special subset of the E S inertial modes are those characterized by having the longest periods. Their force balance is closest to the geostrophic balance in Eq. (4) which results in a flow that is almost invariant with respect to the rotation axis. Following Zhang & Liao (2017), modes of this subset of E S inertial modes will therefore be referred to as being quasi-geostrophic (QG). They correspond to the half-frequencies σ S m1k , and are modes with the index n = 1 in accordance with the previously applied ordering.

Equatorially Antisymmetric Modes
The E A inertial modes are denoted u A mnk and can be again separated into axisymmetric and non-axisymmetric modes. The half-frequencies of the axisymmetric E A modes satisfy (Zhang & Liao 2017) with k ≥ 1, whereas the half-frequencies of the non-axisymmetric E A modes are the roots of with k ≥ 0. Hence, there are k axisymmetric and 2k + 1 non-axisymmetric E A modes for a given k (see Appendix A3 for explicit expressions). Fig. 3 shows examples of a Geostrophic polynomial, E S and E A inertial modes, presenting the structure at the core surface and in meridional sections. Note that these modes have little energy at small cylindrical radius and are well-suited for describing fluid motions in the equatorial region.

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 u0 and a time-dependent flow ut We allow the geostrophic mode, the QG modes (i.e. E S modes with n = 1) and the slowest E A modes (n = 1) to be time-dependent, since these are the modes expected to be excited by the very slow (relative to the rotation time scale) forcing processes in the core (e.g. Zhang & Liao 2017) where [·] denotes the real part. The steady background consists of the axisymmetric and non-axisymmetric inertial modes with n = 1, since forcing processes on long timescales, for example inhomogeneous mantle heat flow, is expected to include both equatorial symmetries 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 parameterize the time-dependence of our core flows. We use only their spatial structure, which we believe to be an efficient way of parameterizing the morphology of rotationdominated 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.

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 where ∇H ≡ ∇ −r ∂ ∂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 Br = Br +Br, 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 ∂ Br ∂t , is then (Gillet et al. 2015a) where e = −∇H · (uBr) is the radial component of the so-called small-scale error, which 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 which 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 equatioṅ where the column vectorsḃ, b, x and e contain the SH expansions of the SV up to degree N˙b, the resolved core field up to degree N b , 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 H b depends on the core field 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, that links the core surface field to the field at 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 based on inertial modes that are most prominent at mid and low latitudes.

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 which allows the computation of the toroidal-poloidal flow expansion from a linear combination of geostrophic polynomials and inertial modes. Hereafter, the column vectorsḃp, bp, xp and ep 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 G k (tp), QG modes, described by coefficients a S m1k (tp), and the slowest E A modes, described by the coefficients a A m1k (tp) in a time-dependent column vector ap ≡ a(tp). Similarly, we store the coefficients a S mnk and a A mnk with n = 1 in a timeindependent column vector a0 that describes the steady background flow. The toroidal-poloidal expansion of flow Eq. (17) at time tp can therefore be written as where M0 and Mt are fixed matrices which contain as columns the toroidal-poloidal expansion of the geostrophic and inertial modes matching the ordering of the coefficients in a0 and ap. Inserting this expression into Eq. (22), the SH expansion of the SV on the core surface at time tp becomesḃ 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 dp, we multiply the SV expansionḃp in Eq. (24) with the matrix Yp, containing the appropriate spatial derivatives of spherical harmonics, such that In this step, it is assumed that the mantle conductivity is small and can be neglected. By adjusting the number of rows in Yp, we can also handle the changing number of SV data over time.
To summarise, we define our model vector to be m = [a T 0 , a T 1 , . . . , a T P , e T 1 , . . . , e T P ] T and combine the P single-epoch expressions in Eq. (25) into one matrix equation of the form where the forward matrix is with Fp = YpH b (bp). This allows us to calculate model predictions d = [d T 1 , . . . , d T P ] T at the locations and times of the observational data.

Inverse problem
In order to infer the core flow from the geomagnetic observations stored in d obs , we minimize a cost function of the model vector m of the form where x p = ( i |xi| p ) 1/p denotes the lp norm of a vector x, at ≡ [a T 1 , . . . , a T P ] T is the vector containing all time-dependent flow coefficients, and et ≡ [e T 1 , . . . , e T P ] T is the multi-epoch small-scale error. The strictly positive parameters λ0, λ S t , λ A t 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 side of Eq. (28) measures the misfit to the ground and satellite SV data. It involves the diagonal matrix W d which contains the estimated data error variance σ 2 i modified following a Tukey biweight scheme in order to downweight outliers (e.g. Constable 1988).
with breakpoint c = 4.685 and ri the residual of the data point (d obs )i relative to the CHAOS-6-x7 geomagnetic field model. The choice of the data error variances σ 2 i is described in section 2.2. We regularize the spatial structure of flows by minimizing an l1 norm of the square root of the mode enstrophies calculated with the diagonal matrices Q0 for the steady part u0, Q S t and Q A t for the time-dependent equatorially symmetric and antisymmetric part of ut. We regularize the time-dependence of the flow using an l2 norm of the flow acceleration calculated by applying a time finite-differencing operator D to the time-dependent part of the model vector. The small-scale error is regularized using a quadratic form based on the inverse of its multi-epoch covariance matrix Ce = Ce(m), which depends on the model vector. We define the matrices D, Q0, Q A t , Q S t and Ce in the next section.
Regarding the minimization of the cost function in Eq. (28), we note that this is a non-linear 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 at each iteration step k until convergence criteria are met. We consider a computed model m k+1 to be converged if the relative change between successive iterations of the data misfit (Gm − d obs ) T W d (Gm − d obs ), the l1 norm of the two parts of the flow Q S t at 1 + Q A t at 1 and Q0a0 1 , and the small-scale error norm e T t C −1 e e t are each smaller than 10 −2 . The regularization matrix R is block-diagonal and constructed from the model of the previous iterate m k according to where the matrices W0, W S t and W A t iteratively implement the l1 norm, and Wa the l2 norm. With R defined in this way, the model parameter estimation usually converged within 20 iterations. The definition of the matrices W0, W S t , W A t , Wa and Ce is given in the following.

Regularization norms
4.4.1 Spatial regularization: l1 norm based on square root of mode enstrophies The flow coefficients a0 and at are spatially regularized using an l1 norm based on the Ekblom measure (Farquharson & Oldenburg 1998) for the steady background part and similarly W S t and W A t for the time-dependent equatorially symmetric and antisymmetric part of the flow by replacing Q0 and a0 with Q S t and at, or Q A t and at, 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. Use of an l1 norm promotes a sparse mode representation of the flow, by favouring a scenario where the mode amplitudes are if possible zero, and reducing the amplitude of the remaining modes. In order to also weakly 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(u) ≥ 0 for a flow u, with the enstrophy, the square of this quantity, defined as We compute the square root of the enstrophy for every individual inertial mode and geostrophic polynomial, and build the diagonal matrices Q0 for the steady background, and Q S t and Q A t 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 Eq. .
where ∆t = tp+1 − tp is the step-size and I is the unit matrix that is compatible with the size of ap. The strength of the regularization of the flow acceleration is controlled by the parameter λa.

Regularization of the small-scale error
We compute the multi-epoch covariance matrix Ce of the small-scale error by projecting the covariance of the unresolved small-scale core field Cb of SH degree N b < n ≤ Nb onto large lengthscales of the SV, via the flow at the previous iteration using the block-diagonal matrix The matrix Hx(xp) is based on the induction equation Eq. (20) and is similar to H b (bp) but takes the flow as the argument to connect the expansions of the magnetic field and the SVḃ We build the covariance matrix for the unresolved small-scale field Cb as a P × P block matrix where each block (Cb)pq for p, q ∈ [1, P ] is equal to the covariance of the small-scale core field of epoch pairs tp and tq with E[·] the expectation, τc(n) the correlation time, σ 2 b (n) the variance of the small-scale core field coefficient of degree n, which we have stored as elements in the vectorbp. Here, we use a Matérn class time correlation function with ν = 3 2 (Rasmussen & Williams 2006). This has previous been used for a similar purpose by Gillet et al. (2015a), it describes an AR-2 process compatible with the observed secular variation spectrum ) and allows abrupt slope changes in the SV (jerks). The parameters τc and σ 2 b are estimated for each degree from the field and SV spectra of the CHAOS-6-x7 model in epoch 2015 (see App. B for details). In the computation, we use the time average of the flow over all considered epochs in H since we expect the covariance of the small-scale error to be only weakly time-dependent when the flow is predominantly steady (Gillet et al. 2015a;Bloxham 1992), this helps improve model convergence.

RESULTS
In this section, we report the results for our reference flow model CoreFlo-LL.1. For this model the truncation degree of the geostrophic polynomials, the E S and E A inertial modes was chosen to be K = 10 and M = 20. To ensure an adequate representation of these modes in the computations a SH truncation degree Nx = 60 was used. The multi-epoch inversion covers the period from September 2000 to January 2018 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 E A modes in the time-dependent part, and the remaining 4300 E S and E A 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 68914 model coefficients: 2 · 4300 = 8600 determine the background flow, P · (10 + 2 · 200 + 2 · 220) = 45050 the time-dependent part of the flow and P · 288 = 15264 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˙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 N b = 14 since higher degrees are contaminated by the lithospheric field .
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 other flow models were designed to test the robustness of CoreFlo-LL.1 by changing aspects of the inversion scheme and the model parameterization. In model 1b, we constrained the time-dependent part of the flow to be purely equatorially symmetric, whilst continuing to regularize in the same way as CoreFlo-LL.1. For 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.
We characterize and compare the derived flows with the help of the following statistics: (1) normalized quadratic measure of the SV misfit between the observational data d obs and our model predictions d = Gm where N d is the number of data, (2) normalized quadratic measure of the SA misfit by computing the central time differenceḋ andḋ obs of the model predictions and data where Nḋ is the number of the SA data,  Table 1. Summary of the computed flow models. Model 1 is our reference solution CoreFlo-LL.1, which is derived using an l 1 norm of the square root of the mode enstrophies, and including time-dependent geostrophic, QG and the slowest E A modes (G+QG+E A ). In model 1b, we use the same l 1 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 parameterization as model 1b, we use as spatial regularization an l 2 norm of the toroidal-poloidal flow expansion that scales as n 3 (Gillet et al. 2009). Definitions of the diagnostics are given in the text. where σγ = 0.24 ms is the standard deviation of the LOD over the considered time period and Gγ is the matrix to compute the LOD from our modeled core flow (e.g. Holme 2015), (4) l1 norm of the square root of the mode enstrophy as a measure of spatial complexity where N is the total number of flow coefficients, (5) fraction fS of the total time-averaged flow power that is contributed by the equatorially symmetric flow u S (geostrophic mode and E S inertial modes) integrated over the core surface where · stands for the time-average, and (6) fraction ft of the total time-averaged flow power that is contributed by the time-dependent part of the flow ut (geostrophic mode, QG inertial modes, and the slowest E A inertial modes for CoreFlo-LL.1), integrated over the core surface.

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, 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 fits the data as well as CHAOS-6-x7, and shows remarkably similar patterns of SV. There is no evidence for systematic offsets in the CoreFlo-LL.1 series, and it successed in capturing 'V' shaped structures in the SV associated with geomagnetic jerks, see for example the φ-component at MBO in 2007 (e.g.  and the VO at latitude 5.964 • degrees, 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 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 at the ground observatories and 1.20, 1.42 and 1.21 nT/yr at the virtual observatories. These histograms provide further evidence that CoreFlo-LL.1 is describing well the observed SV at mid and low latitudes.
In Fig. 5, we show the SV power spectrum, taking September 2015 as an example, a time well within a period of good data coverage provided by the Swarm satellites. 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 Units are (km/yr) −2 . 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.

Time-averaged flow
In Fig. 6 we present the flow from CoreFlo-LL.1 time-averaged over the modeled time interval. The regions poleward of the tangent cylinder in this and subsequent plots have been masked in grey to emphasise that the determined flows are not reliable in this region. The basic timeaveraged 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, slightly 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 background flow is described by both E S and E A flow modes, departures from equatorial symmetry, including equator crossing, are allowed. Nonetheless the basic gyre structure, including its meridional flows, are recovered despite the fact that we only use SV data from mid and low latitudes. We find that 82 % of the time-averaged power at the core surface results from equatorially symmetric flow, which is remarkably close to the value of 82 % 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 14.

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 (0.93 of its time-averaged power due to equatorially symmetric flow). More complex modes (high m and k) are in contrast less pronounced; this is at least partly a consequence of the applied spatial regularization.
In Fig. 8 snapshots of the time-varying flow at three year intervals are presented to illustrate how the flow changes over the studied time interval. We find strong non-zonal azimuthal flow in the equatorial region between latitudes 15 • N and S, that is particularly enhanced around Next, 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 observations ‡ , that have been corrected for known atmospheric and oceanic signals filtered using a one year 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 2006 in agreement with the trend in the observations. ‡ https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html  In the next section we discuss in detail how this time varying flow gives rise to secular acceleration pulses and geomagnetic jerk events.

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 timedependence 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 ∂ 2 B φ /∂ 2 t at MBO in 2006 and 2009, and in ∂ 2 Br/∂ 2 t 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.
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 modeled 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 January 2006, January 2009, September 2013 and September 2016 correspond to well-known SA pulses at the core surface Finlay et al. 2015;Chulliat & Maus 2014). 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), or because our flow does not predict well the SA at high latitudes that will contribute to the SA power in the case of the field model. 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 CHAOS-6-x7 field model, shown in the right column. CoreFlo-LL.1 generates localised patches of strong SA in the equatorial region, for example around 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 care should perhaps be taken interpreting the 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.  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 regions of strong convergence.
The spatial structure of the non-zonal azimuthal flow acceleration at the core surface in 2013 is more clearly seen 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 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 order one year or less. Under the equatorial Pacific in 2014, as seen by comparing Figs. 13 and 9, the sign change in the acceleration occurred as the timevarying 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. 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 scale. More localised jerk events can in this case be due to changes in the sign of the non-zonal azimuthal flow acceleration of smaller lengthscale. Overall, a spectrum of jerk events of different sizes and amplitudes is therefore expected, depending on the size and amplitude of the changes in the underlying non-zonal azimuthal flow acceleration.
In Fig. 13 we also explore the 3D structure of the flow accelerations, as captured by CoreFlo-LL.1. Since the modes on which the model is based are three dimensional, 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, the non-zonal azimuthal accelerations required to explain the observations are found to Figure 10. 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. We find that the acceleration changes sign as we move 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 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 very 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, oscillation in the polarity of the non-zonal azimuthal flow acceleration are seen for example close to 90 • W and 120 • E (cyan lines); local intensifications of the azimuthal flow acceleration are also evident at times of SA pulses (green lines).

Experiments with alternative flow regularizations and parameterizations
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 parameterization and model regularization, as summarized in Table 1, see the related discussion in section 5 for a description of the experimental models. We examined in detail results for the azimuthal flow accelerations in 2013 for these experimental flows, similar in form to Fig. 13. Regarding the time-longitude plots, for both model 1b and model 1c, we can identify many similar features as in Fig. 13. In particular, around 90 • W we see unambiguous alternation in the direction of the azimuthal flow acceleration, and find that the intense eastward and westward accelerations again occur in similar pairs of converging and diverging accelerations. The amplitudes of the accelerations are however 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 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 centers of flow acceleration convergence and divergence. In both model 1b and model 1c the azimuthal flow accelerations remain localised close to the core surface, and they vary little in the axial direction. Their time-longitude plots also suggests 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 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. Further comparisons with numerical dynamo models in a regime appropriate for studying rapid dynamics are required to determine whether or not such transient E A flow accelerations are physically reasonable.

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. Figure 12. Snapshots of the flow acceleration and predicted SA from CoreFlo-LL.1 (left) during increased SA power as indicated in Fig. 11. CHAOS-6-x7 prediction of the SA is also shown for comparison at the same time (right). SA has been truncated at spherical harmonic degree 9 in each case. The grey area masks the region inside the tangent cylinder. (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 flow 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, that 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 parameterized 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 datasets similar to those employed here. We have compared our flow results with those of Barrois et al. (2018b) as well as the corrected flows reported by Barrois et al. (2018a), 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 secular acceleration. 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 independently 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 n 5 , 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 are however sensitive to the choice of flow parameterization 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 anti-clockwise rotation of the time-varying flow under the equatorial western Pacific between 1997 and 2015. In CoreFlo-LL.1 we find the azimuthal flow direction changes from westward to eastward around 2012 in this region, and this occurs via an anti-clockwise 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 eastward to westward 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 greater power at small scales, and the amplitude of our time-varying flows are much larger. Given the restriction of our observations to mid and low latitudes, and because we base our flow parameterization 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 parameterization, 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 are required by the observations. Localised peaks in the azimuthal flow acceleration, associated with flow acceleration convergence/divergences, give rise to SA pulses with alternating polarity. The meridional part of the time-varying low latitude flows is less well constrained by the data alone. Purely E S time-varying flows can account for the secular variation observations, but allowing cross-equatorial flow enables a similar fit to the data with simpler flows.

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 localised 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 could 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 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 Figs. 13 as previously noted by Chulliat et al. (2015).
An alternative source for low latitude flow accelerations is 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 that of 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 first under western Africa and then progressively later at locations to the west. Aubert & Finlay (2019) have argued that such waves may indeed provide an explanation for geomagnetic jerks. In our flows, jerk events are associated with a rapid change in the sign of the low latitude, non-zonal, azimuthal flow acceleration pattern.

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. It is based on ground and satellite magnetic observations and assumes that the time-varying flow can be parameterized primarily in terms of geostrophic and QG modes. We find evidence for strong non-zonal azimuthal flow accelerations in the equatorial region, that change polarity on an interannual timescale. SA pulses at the core surface are associated with localised 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, localised 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 at the equatorial region on the core surface. The meridional component of our time-varying flows are however strongly 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 driving 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 taking into account finite dissipation; these could then instead be used these 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.
(A.11) § 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.

APPENDIX B: DETAILS ON SMALL-SCALE FIELD COVARIANCE
For Cb, we assume at all epochs that the expansion of the unresolved small-scale core fieldbp has zero mean and variance σ 2 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 Nb = 30 in the year 2015 W b (n) = (n + 1) n m=1 g m n (t) 2 + h m n (t) 2 t=2015 (B.1) 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 σ 2 b (n) = 1.10 · 10 9 (n + 1)(2n + 1) exp(−1.26n), (B.2) with N b < n ≤ Nb. The correlation time in years can also be derived from the CHAOS-6 model with (Gillet et al. 2013) where W˙b is the power spectrum of the SV and given by the expression in Eq. (B.1) by replacing the SH coefficients with the respective time derivativesġ m n andḣ m n . An extension towards the degrees N b < n ≤ Nb is achieved by fitting an exponential model for degrees n ∈ [2, 10] τc(n) = 5.00 · 10 2 exp(−0.23n). (B.4)