Satellite magnetic data reveal interannual waves in Earth’s core

We have inverted core surface flow motions for 1999–2021 (Section 4) from models of the geomagnetic field constructed using ground-based and satellite data (11, 12). The calculated flows, obtained under the spatiotemporal constraints of a geodynamo simulation, are mostly compatible with the QG hypothesis (SI Appendix, Fig. S1) and can thus, for a large part, be continued into the core interior. QG flows with relatively large latitudinal length scales ${\ell }_{\lambda }$ at the core surface project onto flows, with small length scales ${\ell }_s$ in the cylindrical radial direction near the equator. Then, the shear in the equatorial plane is the largest in the cylindrical radial direction. Mass conservation implies $u_s\sim u_{\varphi }{\ell }_s/{\ell }_{\varphi }$, with ${\ell }_{\varphi }$ the azimuthal length scale, so that the azimuthal flow component $u_{\varphi }$ dominates. We find recurring quasi-periodic patterns in the equatorial belt, of period $T\simeq 7$ y, propagating westward at a phase velocity $V_{\varphi }\sim 1,500$ km/y or ${25}^{°}/$ y (Fig. 1 A, Lower, and SI Appendix, Fig. S2). They are carried by low azimuthal wave numbers m, predominantly m = 2, and, together with torsional waves, they account for the interannual magnetic signal. Time–latitude diagrams (Fig. 1 A, Upper, and SI Appendix, Figs. S2 and S3) indicate outward propagation at a speed $V_s\sim 200$ km/y along the cylindrical radius s, as well as several zero-crossings in latitude. The latitude of the first one is ${\lambda }^{\mathrm{\star }}\simeq \pm {10}^{°}$. It corresponds to a horizontal length scale at the core surface ${\ell }_{\lambda }=2r_0{\lambda }^{\mathrm{\star }}\sim 1,200$ km, where $r_0=3,485$ km is the outer core radius. Once projected onto the equatorial plane, it gives ${\ell }_s=r_0(1-\cos ({\lambda }^{\mathrm{\star }}))\sim 50$ km.

Assuming that the mantle is electrically insulating, we focus the discussion of our results on the equatorial region. At the core surface, magnetic fluctuations b have to match a potential field so that $\overrightarrow{b}=-\nabla \mathrm{\Phi }$, where $\mathrm{\Phi }$ is the magnetic potential. We have $\left|b_r\right|\sim (l+1)\left|\mathrm{\Phi }\right|/r_0$ and $\left|b_{\varphi }\right|\sim m\left|\mathrm{\Phi }\right|/r_0$, leading to $\left|b_r\right|\sim (l+1)\left|b_{\varphi }\right|/m$. Therefore, the two components have similar amplitude at r = r₀. Conversely, in the interior, the azimuthal component is much stronger than the radial component because ${\ell }_s\ll {\ell }_{\varphi }$. For a solenoidal field, we have near the equator $\left|b_r\right|\simeq \left|b_s\right|\sim \left|b_{\varphi }\right|{\ell }_s/{\ell }_{\varphi }\ll \left|b_{\varphi }\right|$. The relative magnitude of the magnetic and kinetic parts can be estimated in the interior from the linearized induction equationwith η the magnetic diffusivity and ${\mathbf{B}}_0$ the slowly evolving background field. Assuming diffusion does not affect much the rapid changes in the magnetic field, as it is the case in Earth-like geodynamo simulations (13), the magnetic field perturbation in the interior scales as $\left|b_{\varphi }\right|\sim \left|B_{0,r}\right|\left|u_{\varphi }\right|/\omega {\ell }_s$, with $\omega =2\pi /T$. We can transform this result into a magnetic to kinetic energy ratio in the interior

with $\rho \simeq {10}^4$ kg/m³ the core density, $\mu =4\pi \times {10}^{-7}$ H/m its magnetic permeability, and $\left|B_{0,r}\right|=5\times {10}^{-4}$ T an estimate of the radial magnetic field at the core surface in the equatorial region (14). This result agrees with the classification of the waves we have detected as MC waves, whose energy is mostly magnetic. From Eq. 1 below the equator, we get $\left|b_r\right|\sim m\left|B_{0,r}\right|\left|u_{\varphi }\right|/(\pi r_0\omega )$, with m the azimuthal wave number of the magnetic perturbation. As a result, using the continuity of b_r and of $u_{\varphi }$ across the boundary layer, the magnetic to kinetic energy ratio at the core–mantle boundary behaves as

It is $\mathrm{\lesssim }{10}^{-2}$ for observable length scales. Using $\left|u_{\varphi }\right|\sim 4$ km/y (Fig. 1), we find $b_r\sim 1 \text{μ}$T for m = 5, in agreement with interannual field changes observed at the core surface (3, 7). The difference between the kinetic to magnetic energy ratio at the surface and in the interior is one reason why a detection of such waves is difficult from the magnetic field alone, as has been the case for torsional waves (1). Eq. 3 tends to indicate that the direct extraction of waves from magnetic field models is easier for their high m part. As a matter of fact, previous investigations of equatorial waves emphasized magnetic structures with wave number $m=5-6$ (3), whereas the present study, relying on core surface flows, focuses on m = 2 components.

Field changes induced by the recovered flow at interannual periods account well for the oscillations in the rate of change of the magnetic field observed above the Earth’s surface, in particular near the equator. This is well illustrated by low latitudes observatory series (Fig. 2). Still, subgrid processes play a significant role (SI Appendix, Figs. S4 and S5). Maps of the interannual induction at the core surface show a clear correlation between these two sources of magnetic field changes (SI Appendix, Fig. S6). This is a property shared with geodynamo simulations (SI Appendix, Fig. S7).

We have computed linear solutions to the QG model in the presence of a nonaxisymmetric poloidal background field ${\mathbf{B}}_0$, accounting for weak magnetic diffusion as in the Earth (Section 4). The computed eigenmodes include QG inertial (or Rossby) modes at periods shorter than T_A, torsional Alfvén modes at periods around T_A, and QG MC modes with periods longer than T_A. In Fig. 1B, we show time–latitude (Upper) and time–longitude (Lower) diagrams for the surface azimuthal velocity of a QG MC mode with a period $T\simeq 7.16$ y and a much longer decay time of 92.56 y. It is readily seen that the numerical QG MC mode compares well with the flow inferred from observations (Fig. 1). The mode is characterized by small length scales in the latitudinal direction and larger length scales in longitude. The complexity observed in the time–longitude diagram for the flow reflects the mixing of azimuthal orders in the eigenmode, in link with azimuthal variations in ${\mathbf{B}}_0$ (see SI Appendix, Fig. S10 and the local dispersion relation below). Because ${\mathbf{B}}_0$ is nonaxisymmetric, eigensolutions mix azimuthal orders m. Still, motions for interannual eigenmodes are dominated by low m patterns. For our particular choice of background field, it happens that eigenmodes separate between modes whose flow presents only even or odd orders m (Section 4). Among the whole set of calculated MC eigenmodes, the one presented here is dominated by a m = 2 flow structure and presents nine zero crossings in the radial direction. The azimuthal complexity of the magnetic perturbation is higher (Fig. 3), as a result of the magnetic boundary condition (Section 2A).

The mode features propagation toward the equator, similar to the flow recovered from observations. The velocity peaks near the equator, with the first zero crossing at ${\lambda }^{\mathrm{\star }}\simeq \pm {10}^{°}$. The time–longitude diagram shows a clear westward propagation of the mode. Notwithstanding their small cylindrical radial length scale, QG flows project on surface flows with large length scales in the latitudinal direction in the equatorial region of the core surface. This is also seen in Fig. 3 A and B, which illustrate, respectively, the surface flow and magnetic field perturbations of the mode for a snapshot in time. In order to compare the magnetic and the kinetic energies, we transform the magnetic field into a quantity that has the dimension of a velocity: $\mathbf{b}/\sqrt{\rho \mu }$. Its radial component $b_r/\sqrt{\rho \mu }$ at the core surface is a factor six times smaller than the flow, in qualitative agreement with the observations at interannual periods. QG eigenmodes present a strong cylindrical radial shear of the z-invariant azimuthal velocity, in particular in the vicinity of the equator. This translates into a much stronger azimuthal magnetic field perturbation in the bulk of the core, as we expect from the scale analysis outlined in the previous section. Deep in the core, the magnetic field perturbation dominates the velocity (see equatorial slices in Fig. 3 C and D), with the magnetic energy $|\mathbf{b}{|}^2/(2\mu )$ larger than the kinetic energy $\rho |\mathbf{u}{|}^2/2$ by a factor of $\approx 36$ (averaged over the volume). This is in agreement with Eq. 2 and the energy ratio anticipated for MC eigenmodes (10). The distinction between the strength of the perturbation at the surface and in the bulk is more easily seen in the azimuthal rms values, $\langle f\rangle (s,z)=\sqrt{\frac{1}{2\pi }{\displaystyle \oint f}{(s,\varphi ,z)}^2\hspace{0.17em}\text{d}\varphi }$, of the velocity and the magnetic field perturbation, evaluated at the core surface ($z=H=\sqrt{r_0^2-s^2}$) and in the equatorial plane (z = 0), as shown in Fig. 4. In the equatorial plane, the dominant component is along the azimuth, with $\langle b_{\varphi }\rangle (s,0)$ dominating over $\langle b_r\rangle (s,0)$ for $s/r_0\mathrm{\gtrsim }0.5$. Except in the vicinity of the equator, where $\langle b_r\rangle (s,0)$ and $\langle b_r\rangle (s,H)$ tend to superimpose, the magnetic perturbation is much stronger in the bulk than at the surface. Deep in the core, the velocity component of the mode is weaker than the magnetic one everywhere, with $\langle u_{\varphi }\rangle (s)$ less than $\langle b_{\varphi }\rangle (s,0)/\sqrt{\rho \mu }$. Only close to the surface (and in particular the equator), the kinetic part of the mode dominates. At the core surface, the magnetic to kinetic energy ratio of the torsional Alfvén and the QG MC eigenmodes decreases with their frequency (SI Appendix, Figs. S8 and S9).

We use as observations Gauss coefficients from geomagnetic field models. We consider both CHAOS-7 (11) (1999–2021) and COV-OBS-x2 (12) (1840–2020). These are built from ground-based and satellite data, the latter being included in COV-OBS-x2 by means of geomagnetic virtual observatories (31). Both models are projected in time on splines: of order 4 with knots spacing 2 y for COV-OBS-x2 and of order 6 with knot spacing 6 mo for CHAOS-7. Another important difference between the two models is associated with the prior information used to reduce nonuniqueness. CHAOS-7 is a regularized field model that penalizes the third time derivative plus the second time derivative near the endpoints, while in COV-OBS-x2, a stochastic prior coherent with the occurrence of geomagnetic jerks is imposed using temporal cross-covariances. Also note that near zonal coefficients are more heavily damped in CHAOS-7, in order to reduce the leakage of unmodeled signals at high latitudes.

The flow models have been inverted with the pygeodyn data-assimilation algorithm (32). This tool uses an augmented state-ensemble Kalman filter that forecasts a stochastic model of the core surface dynamics (22). The latter is anchored to spatiotemporal statistics derived from the 71%-path geodynamo simulation (13). The core surface flow is inverted from the radial component of the induction equation at the core surface,where overlines stand for the projection onto spherical harmonic degrees $l\le L_b=13$. The flow is truncated at degree L_u = 18. The last term e_r in Eq. 7 stands for errors of representativeness, which incorporate all contributions to the rate of change of the field that involve small length-scale fields (degrees $l>L_u$ and $l>L_b$ for, respectively, the flow and the magnetic field), plus diffusion. The separate contributions from e_r and u to geomagnetic field changes are illustrated in SI Appendix, Fig. S5 with some examples of Gauss coefficients and in Fig. 2 and SI Appendix, Fig. S4 with predictions to several ground-based observatories.

A first inversion has been performed using as observations Gauss coefficients that describe the rate of change of magnetic field (or secular variation) from the COV-OBS-x2 model over 1880–2020. This initial flow model has been used as the starting state for the inversion of a second model over 1999–2021 from the CHAOS-7 field model. The obtained solutions thus extend previous models (22) to more recent epochs. In both cases, secular variation data uncertainties are estimated from the spread within the ensemble of COV-OBS-x2 models and depend on time. Ensembles of 50 realizations have been calculated. The two decades covered by satellite data limit the period range where flow models are best resolved to time scales shorter than 10 y. Within this subdecadal range, a spectral line at $T\approx 7$ y has been isolated in ground-based series recorded in the vicinity of the equator (5). To focus on the source of this interannual signal, we band-pass filter the obtained flow models for periods between 4 and 9.5 y, with an order 4 causal Butterworth filter.

The flow perturbation u over a mean state $[{\mathbf{u}}_0=\mathbf{0},{\mathbf{B}}_0]$ of an inviscid and electrically conducting planetary core under rapid rotation can be described by the linearized momentum equationcomplemented by the induction Eq. 1. To compute eigenmodes of this system for nonidealized mean fields and sufficient spatial complexity, we consider further simplifications. We exploit the columnarity of magnetohydrodynamic flows under rapid rotation (33). Then, the flow can be described as QG, so that

with $\psi =\psi (s,\varphi ,t)$ the stream function and $H=\sqrt{r_0^2-s^2}$ the half height of the fluid column. By projecting the three -dimensional (3D) momentum Eq. 8 onto the subset of QG velocities [9], a reduced model for a fully liquid core is obtained, where the velocity is described by the evolution of ψ in the equatorial plane (8). This reduced momentum equation is combined with a 3D magnetic field that matches a potential field at the core–mantle boundary, i.e., the mantle is assumed to be insulating (10). We choose orthonormal velocity and magnetic field bases for the projection to obtain a standard, instead of a generalized, eigen problem, substantially improving convergence of the eigen solver. For the velocity basis, we use the orthogonal QG inertial eigenmodes, which can be written explicitly (34). The magnetic field basis presented in Gerick et al. (10) is orthonormalized by using a modified Gram–Schmidt method.

We impose a poloidal background magnetic field ${\mathbf{B}}_0$ designed so that it satisfies geophysical constraints associated with the phase velocity $\sqrt{\langle V_A^2(s,\varphi )\rangle }$ of (geostrophic) torsional eigenmodes (1, 14):

The background field is defined aswhere ${\widehat{\mathbf{B}}}_p(l,m,n)$ are the basis vectors given by equation A21 in Gerick et al. (10), normalized to have unit energy. Here, l and m correspond to the spherical harmonic degree and order, respectively, and n is the radial degree. The coefficients ${\beta }_{l,m,n}$ are given by ${\beta }_{1,0,0}\approx -2.3,\hspace{0.17em}{\beta }_{1,0,1}\approx 1.45,\hspace{0.17em}{\beta }_{1,-1,0}\approx 2.04$, and ${\beta }_{1,-1,1}\approx -1.77$, so that ${\left({\displaystyle \int {\mathbf{B}}_0^2}\hspace{0.17em}\text{d}V\right)}^{1/2}\approx 12.7\hspace{0.17em}\text{mT}$. The choice for the coefficients ${\beta }_{l,m,n}$, based on the requirements in Eq. 10, results in a period for the gravest torsional Alfvén mode about 6 y (1). A nonaxisymmetric ${\mathbf{B}}_0$ with m = –1 is chosen, so that $\langle V_A^2\rangle (s)\ne 0$ everywhere, in particular at s = 0. This requirement is crucial for the existence of diffusionless torsional Alfvén modes (35). It seems to be relevant as well for obtaining numerically converged MC eigenmodes in the interannual period range, with or without diffusion (although $\eta \ne 0$ slightly helps for numerical convergence). The background field in Eq. 11 includes the two largest Gauss coefficients in the Earth’s magnetic field exterior to the core [here, the time-averaged CHAOS-7 model (11)], with the correct ratio of their amplitude, but nine times larger so that requirement Eq. 10b is met.

Our background field presents point antisymmetry with respect to $\mathbf{r}=\mathbf{0}$: ${\mathbf{B}}_0(\mathbf{r})=-{\mathbf{B}}_0(-\mathbf{r})$. As a consequence, the QG modes separate between modes for which the flow is either symmetrical or antisymmetrical with respect to $\mathbf{r}=\mathbf{0}$ (and vice versa for the magnetic part of the mode). QG flows are furthermore symmetrical with respect to the equatorial plane [${\mathbf{u}}_e(z)={\mathbf{u}}_e(-z),u_z(z)=-u_z(z)$, where ${\mathbf{u}}_e$ denotes the equatorial flow component]. As a result, ${\mathbf{u}}_e$ is either symmetrical with respect to a rotation of angle π about the axis (${\mathbf{u}}_e(s,\varphi )={\mathbf{u}}_e(s,\varphi +\pi )$) or antisymmetrical (${\mathbf{u}}_e(s,\varphi )=-{\mathbf{u}}_e(s,\varphi +\pi )$). In the first (resp., second) instance, the mode presents only odd (resp., even) m. The same reasoning does not hold for the magnetic part of the mode because it is not equatorially symmetrical. The flow for the selected eigenmode of period 7.16 y has only even m.

Higher orders of the background magnetic field are not included here, to ensure numerical convergence of the eigensolutions. The latter is verified by requiring a decay in the kinetic and magnetic energy of the eigenmodes as a function of polynomial degree, shown in SI Appendix, Fig. S11 for the modes discussed above. In SI Appendix, Fig. S10, we show profiles of $\sqrt{\rho \mu }{V}_A(s,\varphi )$ at several longitudes and a map of $B_{0,r}(\theta ,\varphi )$ at the core surface. Defining $B_0={\left({\displaystyle \int {\mathbf{B}}_0^2}\hspace{0.17em}\text{d}V\right)}^{1/2}$, the Lehnert number associated to this background magnetic field is $\text{Le}=B_0/(\mathrm{\Omega }{r}_0\sqrt{{\mu }_0\rho })\approx 4.3\times {10}^{-4}$. We choose $\eta \approx 3.76\hspace{0.17em}{\text{m}}^2/\text{s}$, within the range of acceptable geophysical values (27, 36), so that the Lundquist number (the number of Alfvén times per magnetic diffusion time $r_0^2/\eta$) is $\text{Lu}=B_0{r}_0/(\eta \sqrt{{\mu }_0\rho })\simeq {10}^5$. As a consequence, the calculated eigenmodes are only secondarily affected by magnetic damping, with decay times much longer than their period.

We focus on several eigenmodes near the target period of 7 y using the shift-invert Arnoldi method. Many other modes at longer and shorter periods exist, but are not of interest in the context discussed here. The calculated QG MC eigenmodes with interannual periods show quality factors (i.e., the ratio of their frequency to their decay rate) larger than 10 and therefore propagate before being damped by Ohmic dissipation.

We write the streamfunction as $\psi (s,\varphi )=H^3{\psi }_1(s,\varphi )$, where ψ₁ is regular at s = 1 in order to ensure regularity of u at the equator, where H = 0. We assume that radial length scales are much shorter than horizontal length scales, or $k{r}_0\gg m/(2\pi )$. In order to derive the evolution equation for ψ, we project the momentum equation onto trial functions (8):

Here, we have kept only the highest derivatives in s in the expression of the magnetic force, in a configuration where the azimuthal field is not significantly larger than the radial one, as suggested by geodynamo simulations outside of the geostrophic cylinder attached to the inner core (37). The curly brackets $\left\{\right\}$ denote the z-integral (38),and the operator ${\mathcal{L}}^2$ is defined as

For a spherical core $dH/ds=-s/H$, so that the momentum Eq. 12, expressed for ψ₁, becomesunder the condition $kH\ll |dH/ds|$. The time evolution of the quantity $\left\{B_s{B}_{\varphi }\right\}$ is given by the linearized induction equation, where we keep, again, only the highest derivative in s:

Here, we have neglected magnetic dissipation, as it plays a secondary role for decadal to interannual geomagnetic field changes (13) and because the calculated eigenmodes present decay times much longer than their periods. Finally, by coupling [15] and [16], an evolution equation for ψ₁ is obtained. Assuming that the radial wavelength of the perturbation is much shorter than the length scale over which the medium ($H,s,\left\{B_{0,s}^2\right\}$) evolves, we getwhere V_A is given by Eq. 5. Under the plane-wave ansatz ${\psi }_1\propto \exp (\text{i}(ks+m\varphi -\omega t))$, it gives

leading to the dispersion relation

This latter can also be written aswith k₀ given by Eq. 4. In the small wavelength limit ($k\gg k_0$), we have at the first order in ${\left(k_0/k\right)}^3$

Then, the wave frequency ω only weakly departs from the Alfvén wave frequency. These waves have been found in dynamo simulations, where they have been called QG Alfvén waves (6, 13).

In the opposite limit $k\ll k_0$, from [19], Rossby waves and MC waves are clearly separated, with frequencies, respectively,where $\text{Le}=V_A/(H\mathrm{\Omega })$ compares the rotation frequency and the Alfvén wave frequency. The period of MC waves is distinct from that of Rossby or QG Alfvén waves only for large length scales ($k\le k_0$). They have frequencies similar to the torsional waves frequency kV_A if their wave number k does not differ too much from k₀.