# Large-scale palaeoflow at the top of Earth's core

Neil Suttie, Andreas Nilsson, Nicolas Gillet, Mathieu Dumberry

► To cite this version:

Neil Suttie, Andreas Nilsson, Nicolas Gillet, Mathieu Dumberry. Large-scale palaeoflow at the top of Earth's core. Earth and Planetary Science Letters, 2025, 652, pp.119185. 10.1016/j.epsl.2024.119185. hal-04938345

HAL Id: hal-04938345

https://hal.science/hal-04938345v1

Submitted on 10 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Large-scale palaeoflow at the top of Earth's core

Neil Suttie<sup>a, *, ID</sup>, Andreas Nilsson<sup>a, ID</sup>, Nicolas Gillet<sup>b</sup>, Mathieu Dumberry<sup>c</sup>

## ARTICLE INFO

Editor: DR J Badro

Keywords:
Core flow
Geomagnetism
Length-of-day

## ABSTRACT

A number of recent studies have made use of statistics derived from geodynamo simulations to constrain the flow at the top of the Earth's core. Here we adapt these methods to infer possible core surface flow solutions over the past 9000 years, using a low degree geomagnetic field model, based on archaeomagnetic and sedimentary archives. Despite the low spatial resolution of the field model, we find that the largest scale aspects of the recovered flow agree well with those derived from historical field models. In particular we see the growth of an eccentric planetary gyre over the past 400 years and find that a similar feature appears on bimillennial timescales, concurrent with pulses in the zonal flow. In the flow model vigorous gyre structures are not persistent but appear as transient features with durations of a few hundred years. The zonal flow is westward over the entire period with a mean value of 0.09°/yr. Finally, we consider the variations in length-of-day that are implied by the flow model under the assumption of geostrophy, and compare the total clock error with that inferred from records of historical solar eclipses. Good agreement is found between the two series of mantle rotation for at least the past 1200 years, with anomalously slow mantle rotation coupled with reduced westward zonal flow. Further back in time a rigorous assessment is hard to make, owing to uncertainty surrounding the exact locations of eclipses.

## 1. Introduction

Since the advent of continuous satellite monitoring of the geomagnetic field, core flow inversions under a variety of dynamical assumptions have revealed the existence of an eccentric planetary gyre in the Earth's outer core, that can explain much of the secular variation on decadal timescales e.g. Finlay et al. (2023), Gillet et al. (2019, 2015b), Pais and Jault (2008). The westward flow of the gyre reaches low latitudes beneath the Atlantic and is pushed to high latitudes in the opposite hemisphere. Since its existence was first proposed, the questions of its longevity and stability have naturally arisen. A number of mechanisms may account for persistent differences in the flow beneath the Atlantic and Pacific hemispheres. Variations in the lower mantle, which could persist over geological timescales, could influence the geodynamo through lateral variations in heat flux (Gubbins et al., 2007; Mound and Davies, 2017, 2023) or electrical conductivity (Dumberry and More, 2020). Gyre structures are also seen in numerical dynamo simulations (Aubert et al., 2013), in the presence of gravitational coupling between the inner core and the mantle (Pichon et al., 2016). On the other hand, there is the possibility that the gyre is a temporary feature, possibly

arising from some type of wave (Hori et al., 2022, 2020) and such transient gyre-like structures have been produced in numerical dynamos with isotropic forcing (Schaeffer et al., 2017). More generally, the question as to whether the present secular variation is typical of its long-term behaviour is of great interest.

The problem of recovering flow within the Earth's core from magnetic observations is seriously underdetermined; even with complete knowledge of the field and its secular variation at the surface of the Earth, no unique solution exists for the flow at the top of the core (Roberts and Scott, 1965). By imposing some dynamical constraints on the type of flow that is expected, it is possible to derive flow solutions, although the inherent non-uniqueness cannot generally be eliminated (Holme, 2015). Even when the large-scale flow is so prescribed as to reduce the non-uniqueness, small-scale flows may contribute to the observed secular variation. Such flows are commonly damped through regularisation, the so-called "large-scale approximation", but the justification for this has been called into question (Baerenzung et al., 2016; Gillet et al., 2015b). (See Schwaiger et al. 2023 for a fuller discussion of these issues).

In order to answer questions regarding the longevity of flow

<footnote><sup>a</sup> Department of Geology, University of Lund, Sweden</footnote>

<footnote><sup>b</sup> Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, Université Gustave Eiffel, ISTerre, France</footnote>

<footnote><sup>c</sup> Department of Physics, University of Alberta, Canada</footnote>

<footnote>* Corresponding author.
E-mail address: neil.suttie@geol.lu.se (N. Suttie).</footnote>
structures in the core, it is necessary to consider changes in the geomagnetic field over centennial to millennial timescales, which can only be accessed through geological and archaeological archives that are capable of recording elements of the past geomagnetic field. Models derived from these archives can only capture the coarsest features of the field, so the problems of non-uniqueness are further exacerbated if they are inverted for core flow. Despite this difficulty, some simple models of the flow at the core's surface on millennial timescales have been proposed. Dumberry and Bloxham (2006) used the 3000 year archaeological field model CALS3K.2 (Korte and Constable, 2005) to investigate changes in the core's angular momentum by constraining the surface flow to be time invariant in all but its zonal equatorially symmetric components. By assuming that the flows are surface expressions of axially invariant geostrophic motions, they calculated the variation in mantle rotation required to compensate the changing angular momentum of the core. They found that while the amplitude and timescales of the changes were similar to length-of-day (LOD) variations inferred from the locations of historical eclipses (Morrison and Stephenson, 2001), the overall fit was poor. They proposed that axial shear arising from thermal and magnetic winds could account for the difference. This hypothesis remains controversial, and in the 71 % path dynamo simulation (Aubert and Gillet, 2021) the angular momentum of the core is well approximated by the surface flow on interannual to millennial timescales, under the assumption of geostrophy (Schwaiger et al., 2024).

Taking the CALS7k.2 (Korte and Constable, 2005) model, Wardinski and Korte (2008) obtained surface flow solutions for the past 7000 years at 50 year intervals by requiring the instantaneous solutions to be either toroidal or tangentially geostrophic. In these solutions the zonal flow switched from westward to eastwards quasi-periodically, with periods of 1050 and 540 years detected in both sets of solutions. Some insight into flow at the top of the core can be afforded by the study of time-longitude plots of the radial field at a given latitude, after removal of the axisymmetric field, which is unaffected by azimuthal motions. Using this method to trace the movement of features in CALS7k.2, Dumberry and Finlay (2007) again found evidence for periods of both eastward and westward movement of magnetic field features over the past 3000 years, most clearly seen at mid-latitudes in the Northern hemisphere. However, a similar study performed on a variety of field models found that the motion of features at 60°N could be best explained by two different rates of westward drift, and saw no clear eastwardly progression of features (Nilsson et al., 2020). They interpreted the slower rate (0.09°/yr) as representative of the mean zonal flow, with the faster rate of 0.25°/yr possibly the expression of a planetary wave.

Recently, another method of constraining coreflow inversions has been developed, relying on kinematic statistics from numerical simulations of the geodynamo (Aubert, 2014; Barrois et al., 2017; Fournier et al., 2011; Gillet et al., 2019). An important realisation was that the largescale, long period features of these simulations remain invariant along a certain path through parameter space, as the path advances towards ever more realistic conditions pertaining to the Earth's core (Aubert, 2019; Aubert et al., 2017). As it is just these large scale, long period features that we might hope to investigate using archaeological field models, the use of a numerical dynamo to constrain the flow is appealing. Already, this approach has established the persistence of the gyre since at least the latter 19th century (Aubert, 2014; Gillet et al., 2019).

Furthermore, the development of the stochastic modelling package pygeodyn (Huder et al., 2019), allows for a simple implementation of this approach. The algorithm uses an augmented state Kalman filter (Barrois et al., 2017) to account for both sub-grid errors and diffusion, which can be significant even on short timescales with high-degree field models. The approach was further developed to incorporate temporal statistics derived from numerical dynamo simulations by Gillet et al. (2019), who produced a model of the flow at the surface of the core extending back to 1880.

It is the possibility of recovering large scale features of the flow over

millennial timescales that we address here. In particular we consider two of the outstanding questions regarding the temporal evolution of the flow at the top of the Earth's core. Does a persistent gyre structure account for the observed secular variation over centennial to millennial timescales, and do the flow solutions agree with LOD as inferred from historical eclipses, under the assumption that the angular momentum of the core can be determined from the flow at its surface?

## 2. Methods

A model of flow at the surface of the core was derived using the freely available modelling package pygeodyn 2.0.0, a geomagnetic data assimilation package (Huder et al., 2019). The inputs are a geomagnetic field model, in this case pfm9k.2 (Nilsson et al., 2022) and the magnetic field and surface flow from a numerical dynamo simulation, which provides the spatio-temporal statistics from which a prior stochastic model is derived. In essence, the model is an attempt to explain the secular variation of the geomagnetic field model by flows that are statistically consistent with the numerical dynamo simulation. It must be stressed that the results that follow depend strongly on the numerical simulation used to define the prior statistics.

The modelling method is explained in detail by Gillet et al. (2019) and only a brief summary is given here. The data assimilation algorithm uses a Kalman filter described by Barrois et al. (2017) based on the radial induction equation at the core mantle boundary for the large-scale secular variation which can be written in the spectral domain as

$$ \dot{\mathbf{b}} = \mathbf{A}(\mathbf{b})\mathbf{u} + \mathbf{e} \quad (1) $$

The vectors **b** and **u** contain the Schmidt semi-normalized spherical harmonic coefficients of the radial magnetic field and surface core flow respectively, and the upper dot denotes the time derivative. The term involving matrix **A** accounts for the secular variation induced by the large-scale flow, while the errors of representativeness, **e**, also given as spherical harmonic coefficients, accounts for the contribution from sub-grid processes and diffusion. The Gauss coefficients for the magnetic field (and errors of representativeness) are given their conventional labels $s_i^m$ and $h_i^m$, with *l* and *m* spherical harmonic degree and order. Coefficients for the toroidal surface flow are labelled $t_i^n$ as in Holme (2015) with *c* and *s*, subscripts denoting the azimuthal phase, as in Barrois et al. (2017). The vector **u** also contains coefficients for the poloidal flow but this will not be discussed in detail. For a given truncation of field and flow **e** is calculated simply as the difference between the SV of the numerical dynamo and that calculated from the action of the truncated flow on the truncated field at each step of the dynamo run, by Eq. (1). The time series of the dynamo run, $\mathbf{u}'(t)$ and $\mathbf{e}'(t)$, provide covariance matrices for each e.g. $\text{cov}(\mathbf{u}')$ and also temporal cross-covariance matrices e.g. $\text{cov}(\mathbf{u}'(t), \mathbf{u}'(t+dt))$, where $dt$ is the sampling interval, here set to 30 years. These matrices are used to build “drift matrices”, which describe the evolution of **u** and **e** as stochastic processes.

The algorithm operates on an ensemble of realizations of the field, in this case 300 samples taken from the field model pfm9k.2, and proceeds in two stages: forecasts which are performed every 10 years and analysis steps every 50 years. The forecast step takes each member of the ensemble separately and, from an initial value for **u** and **e**, calculates new values based on the drift matrices and the addition of some noise. The field at this new epoch is then calculated via Eq. (1), giving a new “core state” consisting of **b**, **u** and **e**. The analysis step interpolates between the ensemble of forecasts for **b** and the magnetic field model, with the interpolation weighted towards the field model so that the magnetic field and SV recovered at this step are constrained by the field model and its uncertainties.

The model pfm9k.2, used as the observations, is a 9000 year field model consisting of an ensemble of models, each expressed as a degree 5 spherical harmonic expansion in terms of Gauss coefficients, given at 50 year intervals. 300 realizations of this model were used up to 1800.
From 1850 onwards 300 realisations of COV-OBS.x2 (Huder et al., 2020) were taken at 50 yearly intervals, truncated at spherical harmonic degree 5. Each Gauss coefficient of each ensemble member of pfm9k.2 is described by a Gaussian process. We take the derivative of the mean function of each as the SV at each 50 year timestep (Riihimäki and Vehtari, 2010).

For the simulation used to derive the prior statistics, we chose to use the 71 % path dynamo (Aubert and Gillet, 2021), which provides both a long enough time-series (10 ka) and short enough time-steps to generate the covariance matrices. Other available path dynamos (Aubert, 2023; Aubert et al., 2017) are less suitable; 0 % path has too long a timestep (90 years) and 100 % path has too short a run. Another suitable option would be 50 % path and using this as the prior was not found to greatly change the results (see Supplementary Fig. S1). While this is to be expected, using a dynamo with different boundary conditions could certainly change the inferred flow and we reiterate that the results that follow are contingent on the choice of geodynamo simulation.

The choice of how much spatial complexity to allow the flow, which we impose by simple truncation of the spherical harmonic basis, is somewhat arbitrary. A lower truncation will see more of the observed secular variation incorporated into the errors of representativeness. A higher level of truncation will generate flow complexity that is not required to reproduce the secular variation, but reflects the numerical geodynamo prior. Within the dynamo simulation, certain features such as polar vortices require degree 7 and 9 zonal flows, so we opt for truncating the flow at degree 10, to preserve these features in the prior. Following Gillet et al. (2019), the errors of representativeness were calculated for the specific case of a degree 5 magnetic field model, and a degree 10 core flow model. We found that a more severe truncation of the flow at degree 5 did not change the low degree flows, on which this study focusses, significantly (see Supplementary Fig. S1). We stress that the higher degree flows are retained so that each individual member of the ensemble appears consistent with the types of flow seen in the dynamo simulation; this does not imply that they can be resolved.

The spatial complexity of the flow is further limited by using the first 80 principal components, reducing the number of independent flow components by one third. The Graphical LASSO (Istas et al., 2023) was implemented to reduce spurious terms in the empirical forecast covariance matrices.

To summarize, the modelled secular variation results from two stochastic processes, **u** and **e** with spatial and temporal covariance matrices

matching the dynamo simulation. These processes are generated in such a way that the total SV matches the SV of the archaeomagnetic field model, within its uncertainties. The errors of representativeness (ER) are necessary to account for the SV generated by inaccessible processes: contributions from sub-grid processes and diffusion. In the dynamo simulation, the contribution of diffusion to the secular variation on decadal to millennial timescale is found to be quite small, accounting for <10 % of the power in the secular acceleration (Aubert and Gillet, 2021). A much stronger contribution of diffusion over these timescales would make the method used here unworkable, as ER would account for the majority of the SV. The results presented below include the diffusive part of the SV as a process with properties derived from the dynamo simulation. They are not definitive solutions for the flow at the top of the core over the past 9000 years, but are consistent with both the archaeomagnetic field model and the specific dynamo simulation that was used to generate the prior statistics.

## 3. Results

### 3.1. Contributions to secular variation

As explained above, the secular variation (SV) in the model is the sum of two processes: SV that is induced from the flow, and a contribution from the errors of representativeness (ER). In Fig. 1, the model SV is shown (red) along with these two contributions, shown in black for the Gauss coefficients $h_1^1$ and $g_2^2$. The green line represents the mean SV of the original archaeomagnetic field model and the black line represents, in (a) and (c), the SV induced by the inverted large-scale flow, while in (b) and (d) it shows the contribution from ER. One prominent feature is the rapid secular acceleration in $h_1^1$ between 800 BCE and 200 BCE and this is almost entirely induced by the large-scale flow (Fig. 1a). This change in $h_1^1$ is linked to the dispersion of an intense flux patch below western Asia, which is discussed later. The variation in $g_2^2$ illustrates how, in this model, ER (black line, Fig. 1d) typically accounts for a significant fraction of the observed SV on centennial time-scales, while most of millennial SV changes are the consequence of variations in the large scale flow (black line, Fig. 1c).

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Modelled secular variation (red) with contribution from induced SV (a,c) and ER (b,d) in black for two Gauss coefficients. Envelopes show the standard deviation of the ensemble. The green line shows the mean SV from the pfm9k.2 field model with the standard deviation shown as green shading in a and c.</figcaption></figure>
## 3.2. Large-scale flows

As we aim to study the core angular momentum and the evidence for a persistent gyre, the most important flow coefficients are the equatorial zonal flow ($t_1^0$), and the toroidal flows $t_{c2}^1$ and $t_{s2}^1$ so it is the evolution of these that we focus on. The three coefficient time series are shown in Fig. 2. Also shown are the mean values from the flow model of Gillet et al. (2019), based on the COV-OBS.x1 field model (Gillet et al., 2015a) at the analysis timesteps of 1900, 1950 and CE 2000. In that model, the strengthening gyre can chiefly be characterised by the equatorial zonal flow, $t_1^0$, increasing from 4 to 9 km/yr westwardly, and a $t_{c2}^1$ component enhancing westward equatorial flow beneath the Atlantic and opposing it beneath the Pacific. The general agreement between the models is notable; although pfm9k2 is tied to COV-OBS.x2 for the recent period, it is still truncated at degree 5, rather than degree 14 as are the COV-OBS models. It is apparent then, that similar large-scale flows can be recovered despite the highly truncated field model being used.

The mean value of $t_1^0$ is seen to be westward (negative) for the whole duration of the model, although the uncertainties encompass the possibility of easterly equatorial flow, and is somewhat slower than the prior distribution, shown in dashed/dotted lines. The difference in the total zonal flow speed between the model and prior is apparent in the time average. In the model, the average of the total zonal flow outside the tangent cylinder over latitude and time is 0.09°/yr westward, whereas in the dynamo prior the same calculation yields 0.15°/yr, equivalent to 9km/yr at the equator. As mentioned in the introduction, a westward drift rate of 0.09°/yr, was detected by Nilsson et al. (2020), who suggested that it represented a mean background flow. That the same value is recovered here, using a different field model strengthens this assertion, and shows that the magnetic field model is informative enough to significantly modulate the prior. The period between CE 600 and CE 1600 is marked by much reduced zonal flow; the average outside the tangent cylinder is <0.05°/yr over this interval. At CE 1600, $t_1^0$ starts to increase in magnitude along with the $t_{c2}^1$ component, both becoming more negative (Fig. 2) and it is around this time that a feature similar to today's gyre becomes apparent. As well as the recent increase in westward flow, there are seen to be four previous periods of enhanced westerly flow at 6400, 3700, 1700 and 200 BCE. At 1700 BCE the configuration of the proposed flow is similar to that of the present day, with $t_1^0$ and $t_{c2}^1$ combining to produce a gyre extending beneath the Atlantic (Fig. 3b) The strong equatorial flow at 200 BCE is not associated with large scale meridional flow (Fig. 3a), but is accompanied by a strong polar vortex. In the geomagnetic field model, the 7th to 2nd centuries BCE see an intense flux patch under western Asia being rapidly

dispersed westwards (Nilsson et al., 2022), an event captured as a rapid swing in declination in Northern hemisphere lake sediment records. The increasing magnitude of $t_1^0$ from 350 BCE to 200 BCE (Fig. 2) seems to be associated with this event and Fig. 1a and c show that much of the SV at this time is attributable to the modelled flow.

## 3.3. Influence of the geodynamo prior

The 71 % path numerical dynamo simulation (Aubert and Gillet, 2021) that provides the spatio-temporal statistics at the heart of the inversions, produces eccentric gyres preferentially at certain longitudes. By rotating the prior model azimuthally, it is possible to examine its influence on the location of the inferred gyre. We rotated the prior azimuthally by angles of 90°, 180° and 270° and ran the pygeodyn algorithm with all other parameters unchanged. For each of these, the mean value of $t_1^0$ was not significantly altered, with each new flow model showing the same peaks in equatorial flow as the original (Fig. 4). After a rotation of 90° E, the inferred gyre at 1700 BCE is pulled slightly to the East (Fig. 5a), but is broadly in the same position as the original, suggesting that its position under the Atlantic hinges on the geomagnetic field model and is data driven. However, after a rotation of 180°, the flow at 1700 BCE does not take the form of an eccentric gyre, and is instead chiefly characterised by equatorial zonal flow (Fig. 5b). This leads us to speculate that the period of enhanced zonal flow at 200 BCE in the original model may in fact represent a gyre extending under the Pacific, a feature which is not supported by the spatial statistics derived from the geodynamo prior. Some support for this notion is found by removing the time average from the original flow model, which leaves a gyre-like structure extending beneath the Pacific hemisphere (Fig. 6b). Furthermore, after a rotation of the prior 180°, a gyre structure extending under the Pacific at 200 BCE becomes apparent (Fig. 6a). This suggestion is admittedly speculative, and emphasises that the results given here are strongly influenced by the prior statistics.

## 3.4. Variations in length-of-day

Under the assumption of geostrophy, the angular momentum of the core can be calculated from the zonal surface flow coefficients, along with the attendant changes in mantle rotation (Jault and Finlay, 2015). The variations in length-of-day (ΔLOD) inferred from our flow model are almost entirely governed by $t_1^0$, as can be seen by comparing ΔLOD in Fig. 7a with Fig. 2. The contribution from $t_3^0$, although used in the calculation, was found to be small, as it only varies slightly throughout the model, staying a more or less constant –1km/yr, and probably

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Evolution of three toroidal flow components over 9000 years, showing the mean and standard deviation of the ensemble. Black dots are from the mean model of Gillet et al. (2019) at the analysis timesteps. Dashed and dotted black lines show the mean and standard deviation of the flow components in the dynamo prior.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Orthographic projections of the northern hemisphere showing the mean of the ensemble of core flow solutions at 200 BCE (a) and 1700 BCE (b). White dashed line shows the projection of the inner core tangent cylinder onto the CMB. Coastlines are shown for reference only.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Evolution of $\epsilon_1^0$ over the past 9000 years using azimuthally rotated geodynamo simulations for the prior distributions.</figcaption></figure>

beyond the resolution of the model. Also shown is the variation in LOD predicted by the 3ka flow model of Dumberry and Bloxham (2006). While there is good overall agreement, despite very different modelling strategies, there is one important difference; there is no long-term trend in LOD in our model over the past 3000 years, with the value at 1000 BCE equal to that at CE 1900. Conversely, the flow inversion of Dumberry and Bloxham (2006) implies a change in LOD due to core-mantle coupling of around 12 ms over the same period, which on integration has a significant effect on the long-term secular deceleration (Mitrovica et al., 2015).

To compare the implied variations in mantle rotation with observations of historical eclipses, it is insightful to integrate ΔLOD derived from core flow back in time. The site at which an eclipse was observed provides the angle through which the Earth has rotated since the observation. This is usually expressed as clock error, or ΔT, the accumulated difference between time measured by standard days of 86400 s, and the variable solar day which depends on Earth's rotation i.e. the difference between fixed Terrestrial Time and variable Universal Time (Stephenson, 2007). For each of the ensemble of 300 core flow models, we sum the derived change in LOD back in time from 1800, close to

when the long-term trend in LOD approached 24 h. As only the change in LOD is predicted from the flow, there is freedom to choose the absolute value of LOD at 1800, although it should be close to zero. From telescopic observation of lunar occultation, LOD appears to be slightly shorter than 24 h at 1800 (Stephenson et al., 2016), and we set the mean LOD of our ensemble at 1800 to -0.5 ms. In Fig. 7c we show the resulting change in ΔT along with the spline fit and untimed eclipse data from Stephenson et al. (2016) after both have had the long-term trend of 1.78 ms /century (Stephenson et al., 2016) subtracted. Also shown is the similarly detrended ΔT from an account of an eclipse possibly observed at Athens in 484 (Hayakawa et al., 2022; Stephenson, 2007), which was not used in the spline fit. Note that the error bars denote absolute limits in ΔT, representing the band of totality: if a total eclipse was observed at a specific site, there are absolute limits on ΔT that would place the site within the path of the umbra.

## 4. Discussion

It is unsurprising, given the coarseness of the geomagnetic field model, and the influence of the prior dynamo simulation, that the results
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Orthographic northern hemisphere projections of the flow obtained after rotating the numerical dynamo prior distributions around the rotation axis by 90° E (a) and 180° (b) at 1700 BCE. White dashed line shows the projection of the inner core tangent cylinder onto the CMB. Coastlines are shown for reference only.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Flow solutions at 200 BCE after rotating the prior numerical dynamo by 180° around the rotation axis (a) and using the original prior statistics, but with the time-average flow removed (b). White dashed line shows the projection of the inner core tangent cylinder onto the CMB. Coastlines are shown for reference only.</figcaption></figure>

presented here remain equivocal. Perhaps the most striking feature of the 9kyr flow model is the occurrence of pulses of westward zonal flow at 3700, 1700, 200 BCE, and the present day, with the best resolved ancient ones at 1700 and 200 BCE having a duration of approximately 250 years. Whether both of these pulses are associated with an eccentric gyre remains uncertain; at 1700 BCE there appears to be a gyre in a similar location as the present, with its location changing only slightly after a 90° eastwards rotation of the prior, rather than solely determined by the spatial statistics inherited from the dynamo simulation. At 200 BCE there is little large-scale meridional flow, but removal of the time average flow hints at the possibility of a gyre structure extending under the Pacific. However, this hypothesis awaits further testing, possibly by using a numerical dynamo with homogenous boundary conditions. In our model, the present day gyre becomes evident around CE 1700 and is seen to strengthen over the past 300 years. Further back in time there is no clear gyre structure under the Atlantic; in this model strong gyre

structures appear as transient features with durations of a few hundred years, supporting previous suggestions that the present asymmetry in secular variation may be unusual (Mound and Davies, 2023; Nilsson et al., 2022). However, given the influence of the prior, evident in Figs. 5 and 6, these results remain uncertain.

Under the geostrophic assumption, variations in LOD required by the zonal flow can be integrated over time to estimate the accumulated clock error. This fits reasonably well with historical observations of eclipses as far back as CE 700. Further back, the flow model diverges rapidly from the spline curve which is forced through the eclipse at CE 454, shown in green in Fig. 7. It is this one data point that is largely responsible for the disagreement between ΔT derived from the modelled flow and the spline fit to the eclipse data. We note that the observation is assumed to have been made at Jiankang although the observation site was never explicitly recorded. The arguments for assigning it to the royal observatory at Jiankang are to be found in Stephenson et al. (2018).
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** (a). Change in length-of-day (ΔLOD) inferred from the flow model (blue) with envelope representing the standard deviation of the ensemble of 300 solutions. The black line shows ΔLOD from the model of Dumberry and Bloxham (2006). (b) Contribution to clock error (ΔT) calculated by integrating ΔLOD back from 1800. (c) Clock error (ΔT) as in (b) compared with eclipse data, after removal of the long term deceleration of 1.78 ms /century (Stephenson and Morrison 2016). Black and coloured bars are absolute limits inferred from total solar eclipses, and the red line the spline fit (all from Stephenson and Morrison 2016). The green bar highlights the 454 eclipse discussed in the text. Blue line shows ΔT obtained from core flow model by integrating ΔLOD, with the envelope representing the standard deviation of the ensemble of 300 solutions. Purple bar shows ΔT (after adjustment for long term trend) for a total solar eclipse in CE 484 if witnessed in Athens (Stephenson, 2007; Hayakawa et al., 2022).</figcaption></figure>

Conversely, around the same epoch, the predicted ΔT from the flow model agrees well with a total eclipse being witnessed in Athens in CE 484; for a discussion of this event and its recording see Stephenson (2007) and Hayakawa et al. (2022). The 454 eclipse timing is seemingly inconsistent with the 484 eclipse being observed at Athens; to reconcile the two events would require a change in ΔT of in excess of 4 min in 30 years, with LOD at least 25 ms shorter than normal, at the very bounds of what may be feasible as a result of core mantle interactions. To validate the core flow model using records of ancient eclipses requires a choice as to which observation can be reliably attributed to a specific site, especially prior to CE 700, when the records are sparse. Stephenson et al. (2018) has developed a convincing argument that eclipses reported in the treatises of Chinese dynasties can usually be placed at the capital, in this case Jiankang (modern day Nanjing), with only a few exceptions. Yet the fact remains that the observation site was not explicitly recorded. It should be noted that the value of ΔT implied by integrating LOD from the mean flow model would place the edge of totality 39 km to the north of Jiankang (see Supplementary Fig. S3), so a report of totality does not seem unfeasible.

While it will always be problematic to confidently derive, on the
strength of historical eclipse records, the small amplitude variations in
ΔT which might arise from core-mantle interactions, the long-term
secular trend is well-established. This trend, a constant deceleration of
1.7ms/cy is some 0.7 ms/cy less than the 2.4 ms/cy which would be
expected by tidal friction alone (Williams and Boggs, 2016). The dif-
ference might be attributed to the postglacial rebound of the polar re-
gions reducing the oblateness of the Earth, and so countering the
deceleration. Mitrovica et al. (2015) showed that the incomplete aver-
aging of millennial core flow fluctuations over the past 3000, based on
the flow model of Dumberry and Bloxham (2006), amounts to an

accumulation of nearly 2 h of clock error. This corresponds to an
equivalent (positive) LOD trend of approximately 0.4 to 0.6 ms/yr,
which implies a larger countering trend from postglacial rebound in the
range of 1.1 to 1.4 ms/cy. The flow model computed in the present study
suggests, in contrast, a much smaller contribution to accumulated clock
error, limited to approximately 15 min. Accordingly, core flows can only
account for a small fraction of the observed LOD trend, implying a
postglacial rebound induced (negative) LOD rate close to 0.7 ms/cy.
This is in keeping with the recent study of Kiani Shahvandi et al. (2024)
which argues for a post-glacial rebound trend of -0.8 ± 0.1 ms/cy and
that the total contribution of core flow to the secular trend of Earth's
rotation over the past three millennia is minimal.

**5. Conclusions**

Despite the low spatial resolution of the geomagnetic field model,
some confidence in the flow solutions is afforded by both the similarity
of large scale flows to those derived by Gillet et al. (2019) and a
reasonable agreement with changes in clock error derived from histor-
ical eclipses over the past 1200 years, when both the field model and the
eclipse records are most trustworthy. The basic picture is one of low
westward zonal flow corresponding to slow mantle rotation from about
CE 800 until CE 1700. Further back in time, there is some uncertainty
surrounding the reliability of eclipse observations and variations in
mantle rotation. Given the difficulty of attributing accounts of eclipses
from ancient China to specific sites, it is perhaps unlikely that ancient
eclipse timings will ever provide a thoroughly reliable test with which to
validate core flow models, or indeed, to argue the case for or against the
dominance of geostrophy within the core, on millennial timescales.

Approximately bimillennial pulses in zonal flow are a distinct feature
of the model. These may either be associated with hemispherical gyres, similar to the present day, or with an intense polar vortex around 200 BCE, however the influence of the prior statistical model make it impossible to make any strong claims regarding such features. The long-term average of the zonal flow, averaged outside of the inner-core tangent cylinder is 0.09° westward, the same magnitude as the westward drift rate previously observed in palaeomagnetic field models spanning the past few thousand years (Nilsson et al., 2020).

### CRediT authorship contribution statement

Neil Suttie: Writing – original draft, Investigation, Formal analysis. Andreas Nilsson: Supervision, Investigation, Funding acquisition. Nicolas Gillet: Visualization, Software, Conceptualization. Mathieu Dumberry: Software.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The research was funded by the European Union (ERC, PALEOCORE, 101125394). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. NS and AN also acknowledge funding from the Swedish Research Council (DNR2020-04813). MD is funded by a Discovery Grant from NSERC/CRSNG of Canada. NG has been funded by ESA in the framework of EO Science for Society, through contract 4000127193/19/NL/IA (Swarm + 4D Deep Earth: Core). He has also been supported by the French Space Agency (CNES) in the context of the *Swarm* mission of ESA. NG is part of Labex OSUG@2020 (ANR10 LABX56).

### Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.epsl.2024.119185.

### Data availability

Data will be made available on request.

### References

* Aubert, J., 2014. Earth's core internal dynamics 1840–2010 imaged by inverse geodynamo modelling. Geophys. J. Int. 197, 1321–1334.

* Aubert, J., 2019. Approaching Earth's core conditions in high-resolution geodynamo simulations. Geophys. J. Int. 219, S137–S151.

* Aubert, J., 2023. State and evolution of the geodynamo from numerical models reaching the physical conditions of Earth's core. Geophys. J. Int. 235, 468–487.

* Aubert, J., Finlay, C.C., Fournier, A., 2013. Bottom-up control of geomagnetic secular variation by the Earth's inner core. Nature 502, 219–223.

* Aubert, J., Gastine, T., Fournier, A., 2017. Spherical convective dynamos in the rapidly rotating asymptotic regime. J. Fluid. Mech. 813, 558–593.

* Aubert, J., Gillet, N., 2021. The interplay of fast waves and slow convection in geodynamo simulations nearing Earth's core conditions. Geophys. J. Int. 225, 1854–1873.

* Baerenzung, J., Holschneider, M., Lesur, V., 2016. The flow at the Earth's core-mantle boundary under weak prior constraints. J. Geophys. Res. Solid Earth 121, 1343–1364.

* Barrois, O., Gillet, N., Aubert, J., 2017. Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics. Geophys. J. Int. 211, 50–68.

* Dumberry, M., Bloxham, J., 2006. Azimuthal flows in the Earth's core and changes in length of day at millennial timescales. Geophys. J. Int. 165, 32–46.

* Dumberry, M., Finlay, C.C., 2007. Eastward and westward drift of the Earth's magnetic field for the last three millennia. Earth Planet. Sci. Lett. 254, 146–157.

* Dumberry, M., More, C., 2020. Weak magnetic field changes over the Pacific due to high conductance in lowermost mantle. Nat. Geosci. 13, 516–520.

* Finlay, C.C., Gillet, N., Aubert, J., Livermore, P.W., Jault, D., 2023. Gyres, jets and waves in the Earth's core. Nature Rev. Earth Environ. 4, 377–392.

* Fournier, A., Aubert, J., Thébault, E., 2011. Inference on core surface flow from observations and 3-D dynamo modelling. Geophys. J. Int. 186, 118–136.

* Gillet, N., Barrois, O., Finlay, C.C., 2015a. Stochastic forecasting of the geomagnetic field from the COV-OBS.x1 geomagnetic field model, and candidate models for IGRF-12. Earth Planets Space 67, 71.

* Gillet, N., Huder, L., Aubert, J., 2019. A reduced stochastic model of core surface dynamics based on geodynamo simulations. Geophys. J. Int. 219, 522–539.

* Gillet, N., Jault, D., Finlay, C.C., 2015b. Planetary gyre, time-dependent eddies, torsional waves, and equatorial jets at the Earth's core surface. J. Geophys. Res. Solid Earth 120, 3991–4013.

* Gubbins, D., Willis, A.P., Sreenivasan, B., 2007. Correlation of Earth's magnetic field with lower mantle thermal and seismic structure. Phys. Earth Planet. Interiors 162, 256–260.

* Hayakawa, H., Murata, K., Sōma, M., 2022. The variable Earth's rotation in the 4th–7th centuries: new ΔT constraints from byzantine eclipse records. Publ. Astronom. Soc. Pac. 134, 094401.

* Holme, R., 2015. 8.04 - Large-scale flow in the core. In: Schubert, G. (Ed.), Treatise on Geophysics. Elsevier, Oxford, pp. 91–113. Second Edition.

* Hori, K., Nilsson, A., Tobias, S.M., 2022. Waves in planetary dynamos. Rev. Mod. Plasma Phys. 7, 5.

* Hori, K., Tobias, S.M., Jones, C.A., 2020. Solitary magnetostrophic Rossby waves in spherical shells. J. Fluid. Mech. 904, R3.

* Huder, L., Gillet, N., Finlay, C.C., Hammer, M.D., Tchoungui, H., 2020. COV-OBS.x2: 180 years of geomagnetic field evolution from ground-based and satellite observations. Earth Planets Space 72, 160.

* Huder, L., Gillet, N., Thollard, F., 2019. pygeodyn 1.1.0: a Python package for geomagnetic data assimilation. Geosci. Model Dev. 12, 3795–3803.

* Istas, M., Gillet, N., Finlay, C.C., Hammer, M.D., Huder, L., 2023. Transient core surface dynamics from ground and satellite geomagnetic data. Geophys. J. Int. 233, 1890–1915.

* Jault, D., Finlay, C.C., 2015. 8.09 - Waves in the core and mechanical core-mantle interactions. In: Schubert, G. (Ed.), Treatise on Geophysics (Second Edition). Elsevier, Oxford, pp. 225–244.

* Kiani Shahvandi, M., Adhikari, S., Dumberry, M., Mishra, S., Soja, B., 2024. The increasingly dominant role of climate change on length of day variations. Proc. Natl. Acad. Sci. 121, e2406930121.

* Korte, M., Constable, C.G., 2005. Continuous geomagnetic field models for the past 7 millennia: 2. CALS7K. Geochem. Geophys. Geosyst. 6, Q02H16. https://doi.org/10.1029/2004GC000801.

* Mitrovica, J.X., Hay, C.C., Morrow, E., Kopp, R.E., Dumberry, M., Stanley, S., 2015. Reconciling past changes in Earth's rotation with 20th century global sea-level rise: resolving Munk's enigma. Sci. Adv. 1, e1500679.

* Morrison, L.V., Stephenson, F.R., 2001. Historical eclipses and the variability of the Earth's rotation. J. Geodyn. 32, 247–265.

* Mound, J.E., Davies, C.J., 2017. Heat transfer in rapidly rotating convection with heterogeneous thermal boundary conditions. J. Fluid. Mech. 828, 601–629.

* Mound, J.E., Davies, C.J., 2023. Longitudinal structure of Earth's magnetic field controlled by lower mantle heat flow. Nat. Geosci. 16, 380–385.

* Nilsson, A., Suttie, N., Korte, M., Holme, R., Hill, M., 2020. Persistent westward drift of the geomagnetic field at the core-mantle boundary linked to recurrent high-latitude weak/reverse flux patches. Geophys. J. Int. 222, 1423–1432.

* Nilsson, A., Suttie, N., Stoner, J.S., Muscheler, R., 2022. Recurrent ancient geomagnetic field anomalies shed light on future evolution of the South Atlantic Anomaly. Proc. Natl. Acad. Sci. 119, e2200749119.

* Pais, M.A., Jault, D., 2008. Quasi-geostrophic flows responsible for the secular variation of the Earth's magnetic field. Geophys. J. Int. 173, 421–443.

* Pichon, G., Aubert, J., Fournier, A., 2016. Coupled dynamics of Earth's geomagnetic westward drift and inner core super-rotation. Earth Planet. Sci. Lett. 437, 114–126.

* Riihimäki, J., Vehtari, A., 2010. Gaussian processes with monotonicity information. In: Yee Whye, T., Mike, T. (Eds.), Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pp. 645–652. PMLR, Proceedings of Machine Learning Research.

* Roberts, P.H., Scott, S., 1965. On analysis of the secular variation 1. A hydromagnetic constraint: theory. J. Geomagn. Geoelectr. 17, 137–151.

* Schaeffer, N., Jault, D., Nataf, H.-C., Fournier, A., 2017. Turbulent geodynamo simulations: a leap towards Earth's core. Geophys. J. Int. 211, 1–29.

* Schwaiger, T., Gillet, N., Jault, D., Istas, M., Mandea, M., 2024. Wave-like motions and torques in Earth's core as inferred from geomagnetic data: a synthetic study. Phys. Earth Planet. Interiors 346, 107104.

* Schwaiger, T., Jault, D., Gillet, N., Schaeffer, N., Mandea, M., 2023. Local estimation of quasi-geostrophic flows in Earth's core. Geophys. J. Int. 234, 494–511.

* Stephenson, F.R., 2007. Variations in the earth's clock error ΔT between AD 300 and 800 as deduced from observations of solar and lunar eclipses. J. Astronom. Hist. Heritage 10, 211–220.

* Stephenson, F.R., Morrison, L.V., Hohenkerk, C.Y., 2016. Measurement of the Earth's rotation: 720 BC to AD 2015. Proc. R. Soc. A Math. Phys. Eng. Sci. 472, 20160404.
Stephenson, F.R., Morrison, L.V., Hohenkerk, C.Y., 2018. The provenance of early Chinese records of large solar eclipses and the determination of the Earth's rotation. J. Hist. Astronomy 49, 425-471.

Wardinski, I., Korte, M., 2008. The evolution of the core-surface flow over the last seven thousands years. J. Geophys. Res. 113, B05101. https://doi.org/10.1029/2007JB005024.

Williams, J.G., Boggs, D.H., 2016. Secular tidal changes in lunar orbit and Earth rotation. Celest. Mech. Dyn. Astron. 126, 89-129.