Forecasting yearly geomagnetic variation through sequential estimation of core flow and magnetic diffusion

Earth’s internal magnetic field is generated through motion of the electrically conductive iron-alloy fluid comprising its outer core. Temporal variability of this magnetic field, termed secular variation (SV), results from two processes: one is the interaction between core fluid motion and the magnetic field, the other is magnetic diffusion. As diffusion is widely thought to take place over relatively long, millennial time scales, it is common to disregard it when considering yearly to decadal field changes; in this frozen-flux approximation, core fluid motion may be inferred on the core–mantle boundary (CMB) using observations of SV at Earth’s surface. Such flow models have been used to forecast variation in the magnetic field. However, recent work suggests that diffusion may also contribute significantly to SV on short time scales provided that the radial length scale of the magnetic field structure within the core is sufficiently short. In this work, we introduce a hybrid method to forecast field evolution that considers a model based on both a steady flow and diffusion, in which we adopt a two-step process: first fitting the SV to a steady flow, and then fitting the residual by magnetic diffusion. We assess this approach by hindcasting the evolution for 2010–2015, based on fitting the models to CHAOS-6 using time windows prior to 2010. We find that including diffusion yields a reduction of up to 25% in the global hindcast error at Earth’s surface; at the CMB this error reduction can be in excess of 77%. We show that fitting the model over the shortest window that we consider, 2009–2010, yields the lowest hindcast error. Based on our hindcast tests, we present a candidate model for the SV over 2020–2025 for IGRF-13, fit over the time window 2018.3–2019.3. Our forecasts indicate that over the next decade the axial dipole will continue to decay, reversed-flux patches will increase in both area and intensity, and the north magnetic (dip) pole will continue to migrate towards Siberia.


Introduction
Earth's time-dependent internal magnetic field, commonly referred to as the core or main field, is generated through turbulent motion of the fluid and electrically conducting iron alloy comprising its outer core. Although this field is generated within the core, it permeates the core-mantle boundary (CMB) and the overlying solid mantle, allowing it to be observed at Earth's surface and above. Downward continuation of historical and recent observations allows the construction of models of the core field at the edge of the source region, which describe its spatial variability as well as its time-dependence over the past decades to centuries (e.g., Jackson et al. 2000;Lesur et al. 2008;Sabaka et al. 2015;Gillet et al. 2015;Finlay et al. 2016b). While temporal variation of the field, referred to as secular variation (SV), is ultimately the result of core fluid motion and magnetic diffusion (e.g., Jackson and Finlay 2015;Holme 2015), the latter contribution is often neglected under the frozen-flux assumption (Roberts and Scott 1965). With additional constraints, this approximation has allowed various workers to use geomagnetic observations to infer core fluid motion along the CMB (e.g., Vestine et al. 1967;Whaler 1980;Le Mouël 1984;Bloxham 1988;Maus et al. 2008;Lesur et al. 2010;Kloss and Finlay 2019). However, in reality, the flow described in all such models remains poorly resolved due to the effects of limited resolution and measurement errors (e.g., Rau et al. 2000;Bärenzung et al. 2016).
Such flow models have been used to predict field evolution on yearly to decadal time scales, for example by calculating the field change that is expected when these fluid motions persist in time (Beggan and Whaler 2009;Whaler and Beggan 2015). In contrast, Aubert (2015) utilised simulations of the geodynamo to forecast field evolution, an approach that is computationally more expensive, but can directly account for both core fluid flow and magnetic diffusion. Another forecasting method is given by Barrois et al. (2017), who do not employ geodynamo simulations directly, but instead predict field change using the statistical properties of the so-called Coupled Earth simulation (Aubert 2013). Various authors have also implemented the use of (ensemble) Kalman filtering for geomagnetic forecasting (Beggan and Whaler 2009;Fournier et al. 2015; Barrois et al. 2017Barrois et al. , 2018Bärenzung et al. 2018;Beggan and Whaler 2018), an approach that also enables a prediction of the uncertainties associated with the core field, and for the forecast to be updated rather easily when new data become available.
There are several benefits to creating models that can estimate the future state of the core-generated field. Firstly, we can test the plausibility of specific mechanisms by retrospectively comparing field forecasts with geomagnetic observations (referred to as hindcasting). Secondly, more practical advantages relate to the creation of forecasts applicable to the widespread use of global magnetic maps for orientation and directional information by mobile devices and industry. These publically available maps of the present day magnetic field (e.g., the International Reference Geomagnetic Field (IGRF) or the World Magnetic Model (WMM)) are updated only once every 5 years, and intervening epochs must be approximated by forecasts. The more accurate the forecast the better; SV that is not well predicted may lead to localised error, mitigated only by the release of an out-of-cycle update as occurred with the WMM in 2019 (Chulliat et al. 2019). In addition, forecasts allow the tracking of regions hazardous to satellites, such as the South Atlantic Anomaly (SAA), where low-orbit spacecraft are more prone to damage from solar radiation (Heirtzler 2002). The SAA is a region on Earth's surface with anomalously low intensity (typically less than 30,000 nT) (e.g., Hartmann and Pacca 2009), and currently covers parts of South Brazil and Uruguay .
In this work, we introduce a novel hybrid method for forecasting field evolution on yearly time scales, by combining the work of Whaler and Beggan (2015) and Metman et al. (2019). While magnetic diffusion is often considered negligible on short time intervals, we model it alongside core fluid motion, motivated by the fact that diffusion (like core fluid motion) can explain field evolution over several decades (Metman et al. 2019), and that a description of both physical processes governing SV (i.e. fluid flow and diffusion) could produce forecasts of increased accuracy. Several authors have already forecast field evolution by accounting for magnetic diffusion (Aubert 2015;Barrois et al. 2017). However, these models rely on parametrisations or computational models run in a non-Earth-like regime (e.g., Christensen et al. 2010), introducing additional sources of uncertainty in the outcomes.
"Methods" section below contains our forecasting strategy in which we first consider the special case of frozen-flux (negligible magnetic diffusion), which we then extend by including magnetic diffusion. "Results" section contains several hindcasting tests, and demonstrations of how diffusion improves hindcast performance. We then apply the method to produce hybrid geomagnetic forecasts for 2019-2025, and describe our IGRF SV candidate model. In "Discussion" section, we reflect on our results and their implications.

Methods
Here, we present our hybrid method of forecasting core field evolution over yearly timescales. The evolution of radial magnetic field is governed by where η is the magnetic diffusivity (assumed uniform), ∇ H = (∇ −r∂ r ) is the horizontal gradient, u is the horizontal core fluid velocity and r is radius. Here we set η = 18 km 2 /yr , about 0.6 m 2 s −1 (Pozzo et al. 2012), consistent with the work of Metman et al. (2019). We assume that B r and its time derivative are known up to degree 14 on the CMB from the CHAOS-6-x9 model (Finlay et al. 2016a). The challenge is to find a large-scale flow and a radial structure of B r within the core, both of which are themselves unobservable, which lead to the observed SV. Although this is a linear equation, diffusion and core flow induction are interlinked; such dynamics may, for example, govern reversed-flux patches, with a combination of advection and diffusion controlling the location and mechanism of magnetic flux expulsion from the core. (1) Here, we consider a simplified two-step model, in which the effects of core flow and diffusion are treated sequentially and separately as optimised initial value problems for some time T. We first choose a time window, hereafter referred to as the model fitting period [T 1 , T ] with T 1 < T , over which both models are optimised for a residual computed on a regular temporal grid spacing of 1 month. The models are then allowed to evolve independently over a new time window, [T , T 2 ] with T 2 > T ; summing the contributions from core flow and diffusion then produces our estimated forecast for the geomagnetic field. In many of our examples, T is chosen as 2010.0 in order to assess forecasting techniques against the true field evolution of 2010.0-2015.0. Several model fitting periods are used, their length ranging from 1 to 9 years.
Step 1: modelling core flow only In this section, we present our determination of core surface flow, using an approach largely based on Whaler and Beggan (2015). We neglect magnetic diffusion and describe the rate of change of the magnetic field for a perfectly conducting fluid immediately below the CMB as As there is only one equation (for B r ) and two unknowns (the horizontal components of flow) any inference for core flow is non-unique (see e.g., Holme (2015)); to infer a single structure of flow we further assume that the flow is both large-scale and steady over both the model fitting and forecast periods. We fit the flow using B r and Ḃ r as prescribed by the CHAOS-6-x9 field model (Finlay et al. 2016a), which is expanded up to spherical harmonic degree L = 14 . Additionally, u is assumed solenoidal (describing an incompressible fluid) allowing its partitioning into toroidal and poloidal parts, each of which is also expanded in spherical harmonics up to degree L u (see e.g., Whaler 1986). We adopt the 'strong' norm of Bloxham (1988), penalising flow complexity with a damping parameter γ.
Our approach to core flow modelling begins by minimising a least-squares (L2-norm) damped quadratic cost function by taking where the matrix A corresponds to the equations of condition generated from the Gaunt-Elsasser matrix representation of the toroidal and poloidal flow, C is the data covariance matrix, ġ is the vector containing the time derivative of all Gauss coefficients, each evaluated on the temporal grid within the model fitting period, and D is the regularisation matrix scaled with a damping parameter γ = 5 × 10 −5 (see for more details, Beggan and Whaler 2008;Whaler and Beggan 2015). For simplicity, we set the values of the diagonal of C to 1 (nT/yr) 2 , though we could include more realistic estimates of the error on each Gauss coefficient (e.g., Lowes and Olsen 2004;Fournier et al. 2015). We point out that we explicitly ignore errors and uncertainties due to truncation or small scale flows ( L > 14 ) and assume that their effect will be captured within the diffusion part of the modelling process. However, as we show later, the fit of the model to the data is acceptable in any case.
The residuals of the L2-norm fit are then used to further refine the fit of the flow model to the input Gauss coefficients with an iteratively reweighted L1-norm method . The vector of residuals ( e =ġ − Aû) are used to create a matrix, R , with diagonal elements R ii = √ 2/|e i | . For each iteration, k, the solution is The same value of γ is used for all iterations. The superscript hat, applied here to the flow, denotes a best-fit estimate. The solution is iterated fifteen times to ensure convergence (which typically occurs within ten iterations). This produces a L1-norm best-fit flow model for a single time epoch.
To compute the steady flow over the entire time period, the matrix solution for the final iteration from each epoch are summed and then solved in a least-squares inversion. At this point, there is an additional damping applied using γ = n · 5 × 10 −3 , where n is the number of epochs. The steady flow model ( û ) can be used to compute the end-member case of purely advective secular variation (described only by fluid flow), denoted by Bu r (Eq. 2); although the flow is steady, any SV induced is not steady since the magnetic field varies with time. Note that no other geophysical constraints (e.g., tangentially geostrophic or helical) have been placed on the flow.
An example of a steady flow model, derived using 13 months of SV data from 2018.3 to 2019.3, and the sum of the square of the residuals (SSR) are shown in Fig. 1. The left panel shows the modelled steady flow on the coremantle boundary. The right panel illustrates the SSR values for a single epoch (2018.3) from the initial L2-norm inversion and the improvement of the fit as the model is iteratively updated via the L1-norm algorithm for 15 iterations. Over 13 months, the weighted mean fit to the SV Gauss coefficients is 0.01 nT, while the normalised standard deviation is 0.31. The histograms of the residuals are very peaked and much more Laplacian than Gaussian, as noted in Walker and Jackson (2000). (4) Using these fluid flow models, core field forecasts are realised using a first-order Taylor series of B r under the assumption that fluid flow modelled over the fitting period remains steady over the forecasting period as with the integer k an index corresponding to the time discretisation of the forecast period. This description allows a field forecast to be computed iteratively, using the field state and therefore the flow contribution to SV at t k−1 . For compatibility with the end of the modelling window, the initial field state at the start of the forecast period at time T (i.e. for k = 0 ) is prescribed by CHAOS-6-x9. For simplicity, the time step over the forecasting period t k − t k−1 is once more set at 1 month for all k.
As the parameter L u is increased, we might expect the fit of the steady flow over the model fitting period to improve, and thereby possibly improve the forecasts. However, as we will show later, this happens only up to L u = 8 , thereafter any improvement is small. This result might be attributed to our disregard of processes such as flow acceleration and radial flow, or to our means of regularisation; however, we consider below how accounting for magnetic diffusion can improve the model fit.

Step 2: inclusion of magnetic diffusion
For the hybrid model, we now include diffusion so that where Ḃ d r = η r ∇ 2 (rB r ). We first apply step 1 to optimise a steady flow û . However, this flow model will in general never describe the evolving magnetic field over the model fitting period perfectly, and therefore has an associated residual both in terms of the field itself ( B r −B u r ), and its time derivative ( Ḃ d r =Ḃ r −B u r ). To model magnetic diffusion, we find a radial structure of magnetic field whose diffusion matches as closely as possible the residual, thereby correcting the model from step 1. Magnetic diffusion is modelled using the formalism given by Metman et al. (2019), see specifically Eqs. (8) and (15) of their work, who provide a forward solution to the initial-value diffusion problem in terms of Galerkin polynomials (Chen et al. 2018;Li et al. 2018). While their study considers fitting the time-dependent B r using this Galerkin solution, we note that here we would need to evolve the approximate Eq. 5, unavoidably introducing an error compared with the true frozen-flux evolution. Instead we use the SV residual (which is known exactly) and for which the matrix exponential notation of Metman et al. (2019) can be directly applied. Given an initial structure of magnetic field, our forward solution of the first time-derivative of magnetic field is then where Y α is the real-valued Schmidt semi-normalised spherical harmonic of degree 1 ≤ l α ≤ L , and order 0 ≤ m α ≤ l α , which has either azimuthal sine or cosine dependence, and where L is the degree of truncation of the core field (here we set L = 14 ). The angles θ and φ denote, respectively, colatitude and longitude, the coefficients q α describe the initial structure of magnetic field throughout the core, and expm[·] is the matrix exponential. The reader is referred to Metman et al. (2019) for a detailed description of the radial basis functions ξ α (r) (whose dimension N corresponds to the number of radial modes used to describe the radial profile of magnetic field), and that of the associated matrix H α relating to the radially dependent part of the Laplacian in Eq. (1).
Our aim is to best explain the residual with diffusion, therefore the optimised coefficients q α are obtained through the minimisation of the quadratic  2019), who penalise large amplitudes of the magnetic field. Here, because we are ultimately penalising the residuals of the first time derivatives, it is consistent to regularise based on the same quantity, the square of the first time derivative (assuming only magnetic diffusion) of the initial field at time T integrated over all space. The norm can also be interpreted as a penalisation of small spatial scales in the radial profile, similar to the norm used by Bloxham (1988) to minimise flow heterogeneity. The cost function is expressed in terms of our model coefficients, reducing to with q = (q 1 ,q 2 , . . . ,q L(L+2) ) T , W a diagonal weighting matrix related to Simpson's rule for the numerical integration in time (see Appendix B in Metman et al. 2019), a blockwise-diagonal forward mapping describing purely diffusive SV, where each block is defined as and t k is the kth nodal point in time. While in this notation a factor η 2 would appear in the regularisation term in Eq. (9) (a result of taking the time derivative of B 0 ), this scalar has been absorbed within the damping parameter . The least-squares cost function is then minimised as where ġ d ≡ġ −ĝ u , and ġ and ĝ u are the secular variation Gauss coefficients from CHAOS-6 and the best-fit estimate of the flow model of step 1, respectively. Having estimated both the fluid flow and magnetic diffusion over the modelling period, we are now in a position to make hybrid forecasts of the core field using a first-order Taylor series where the contributions Bu r and Bd r are computed through forward continuation of their respective descriptions obtained for the model fitting period.

Results
In this section, we use the methods described above applied over several choices of model fitting window (2009-2010, 2007-2010, 2005-2010, 2003-2010, 2001-2010) to the time-dependent CHAOS-6-x9 field model (Finlay et al. 2016a). By hindcasting for the period 2010-2015, we assess whether the addition of diffusion improves the ability of frozen-flux core flow models to forecast beyond 2010 by comparing with the actual evolution of the field according to CHAOS-6-x9. We quantify the accuracy of the hindcast by evaluating the global RMS residual in |B| , which we denote R , averaged over Earth's surface and over the forecasting period 2010-2015: where for our test case where g m l and h m l are the time-dependent Gauss coefficients from CHAOS-6-x9, ĝ m l (t) and ĥ m l (t) are calculated from our models, and the factor of 5 arises from the 5-year period.
We point out that the CHAOS-6-x9 model is derived from geomagnetic measurements up to April 2019, an epoch well beyond our model fitting periods. Over the model fitting periods, time derivatives of the field are therefore known to higher precision than when only data up to 2010 (i.e. the end point of our model fitting periods) would be employed, possibly resulting in overestimated forecasting accuracy over the period 2010-2015.
To assess the magnitude of this bias, we also compare our results with a second set of forecasts obtained with the older CHAOS-3 field model (Olsen et al. 2010) (which is derived from data up to 2010); we use the metrics above to calculate the accuracy of these forecasts with respect to the newer CHAOS-6-x9 model. The global RMS error of our purely steady-flow-driven hindcasts (represented by N = 0 , i.e. no diffusion modes) is given as a function of time in Fig. 2. These forecasts have been derived with 1, 3, 5, 7, or 9 years of SV data, with L u = 8 (a and c) and L u = 14 (b and d), and are evaluated at either Earth's surface (a and b) or at the CMB (c and d). Also indicated are the forecasts computed with the linear extrapolation given by IGRF-11 (thick dashed curves) and the one computed with CHAOS-3 (thick dot-dashed curves; for the model fitting window 2009-2010). Note that the IGRF-11 only provides the SV Gauss coefficients up to degree 8, and therefore we have set higher-degree Gauss coefficients to their value in 2010 over the forecasting period (a so-called no-cast).
Other line styles correspond to the length of model fitting window. For all forecasts, it can readily be seen that the global error increases monotonically with time, where this error grows in a much more linear fashion when evaluated at the CMB (Fig. 2c, d) than at Earth's surface (Fig. 2a, b). Moreover, we find that the choice of model fitting window affects the forecast accuracy more strongly when it is evaluated at Earth's surface, and that a shorter modelling period generally yields more accurate hindcasts. In all cases, the best hindcasts are produced by using the smallest model fitting window at the CMB all flow-based forecasts are more accurate than the IGRF-11 prediction. Additionally, by comparing panels a, b with c, d of Fig. 2 it can be seen there is no significant error reduction when the spatial flow complexity is increased by changing L u = 8 to L u = 14 . Lastly, we note that the use of CHAOS-3 during the model fitting step results in slightly higher forecast errors, although clearly these forecasts still outperform those obtained with IGRF-11. The RMS forecast errors at 2015.0 for the (CHAOS-6-x9) flow-only predictions are summarised in Tables 1 and 2.  Fig. 3b shows that at l = 5 , the model with L u = 4 has a lower error than the model with L u = 8 . Figure 3a, b both shows that the IGRF-11 prediction is comparable with most of our frozen-flux hindcasts for 8 ≤ l ≤ 14 , although interestingly for the case of Fig. 3a the simple IGRF-11 extrapolation gives a better fit for l < 8 . Finally, while the CHAOS-3 and CHAOS-6 error spectra are generally similar, the former shows increased errors at degrees 2, 6, and 7 for the 2009.0-2010.0 period (Fig. 3b).
Hybrid hindcasts for 2010.0-2015.0 As more spatially complex steady flows do not necessarily provide more accurate forecasts (increasing L u beyond 8 does not significantly reduce the RMS hindcast error), we now assess whether including diffusion reduces the global error. We use the hybrid forecasting scheme detailed above, for which we have set damping parameters between 10 −9 and 10 −8 (depending on the model fitting window) to regularise the diffusive solution (Eq. 11). These values ensure we invert matrices with a (2-norm) condition number < 10 8 , such that our solutions are computationally tractable when using double precision variables. We need to choose the number of radial modes, N, that describe the initial structure of magnetic field that will subsequently diffuse. Here we aim for parity between the number of parameters describing both core flow and diffusion, in order that both methods are treated on an equal footing. The number of core flow parameters is 2L u (L u + 2) , which is 448 when L u = 14 , whereas the number of diffusion parameters is NL(L + 2) . When L = 14 , choosing N = 2 then gives parity; we also test N < 2 to see the effect of reducing the resolution of the diffusive model. We also note that the choice N = 2 represents only large spatial scales within the core, unavoidably we then neglect diffusion of radially fine-scaled features possibly degrading the diffusion model fit. However, we find that larger N are generally associated with rapid temporal variation of the core field predictions, associated with large forecast residuals. A conservative choice for N is therefore in line with our primary objective of improving geomagnetic forecasts. Figure 4 shows the RMS hindcast error for the period 2010.0-2015.0, evaluated at Earth's surface or at the CMB (a, b and c, d, respectively), and obtained with L u = 8 (a, c) or L u = 14 (b, d); in all cases N = 2 . Compared to the purely flow-based predictions (Fig. 2), we find the hybrid scheme has a smaller global error at all times. As before, the best hindcasts are achieved by using the shortest model fitting window, 2009-2010. Most of Table 2 As Table 1   As with the flow-only method, we find that increasing L u from 8 to 14 within the hybrid scheme generally yields a negligible reduction in global hindcast error as the relative improvement is typically 1%. While the CHAOS-3 forecasts are characterised by larger prediction errors than the CHAOS-6 hindcasts, in particular at the CMB, these nevertheless outperform IGRF-11. The purely flow-based and hybrid RMS hindcast errors at 2015.0 (for CHAOS-6-x9) are listed in Table 1 (at Earth's surface) and Table 2 (at the CMB) for several model fitting windows.
Error spectra of the global RMS error at Earth's surface are given in Fig. 5 In particular, the increased errors at degree 2, 6, and 7 due to the use of CHAOS-3 (Fig. 3)  period, which is associated with a geomagnetic jerk in 2014 (Torta et al. 2015), i.e. a step change in the secular acceleration. The evolution of g 0 14 as given by CHAOS-6-x9 is not well matched by the purely flow-driven and IGRF-11 forecasts, possibly because l = 14 is at the truncation limit of the model, whereas the hybrid forecast is markedly more accurate. The corresponding CHAOS-3 forecast is poorer, although in contrast to the steady flow and IGRF forecasts it still captures the observed growth of g 0 14 . The fit from IGRF-11 is poor as only degrees 8 or less have a non-constant extrapolation over the hindcast period.
Lastly, we comment on the spatial distribution of hindcast errors compared to CHAOS-6. Figure 7 shows the unsigned error in B r at the CMB and at selected epochs for the hindcasts calculated with L u = 14 , the 2009.0-2010.0 model fitting period, and with either N = 0 (a and b) or N = 2 (c and d). It is evident that the hybrid hindcast produces much reduced errors everywhere (compare Fig. 7a, b with 7c, d), where particularly large error reductions are found below the west of Chile, the South Atlantic, and Central Asia. For both core-flow and hybrid models, the largest error amplitudes are found along the − 90 • and 90 • meridians; the smallest hindcast errors are generally in the Pacific Hemisphere. Similarly, Fig. 8 shows for the same hindcast the unsigned B r error at Earth's surface. Clearly, the error reduction acquired with the hybrid approach is not as significant as at the CMB since the improvements are only at high degree.
At Earth's surface, we find a hemispherical asymmetry in error distribution, with the largest errors within the Indo-Pacific Hemisphere and the smallest in the Atlantic.  Figure 9 shows predicted time series of selected Gauss coefficients. For g 0 1 and g 0 14 we predict a continuation of their behaviour from 2010.0-2015.0 (Fig. 6), namely decay (in absolute value). In contrast, g 0 7 has strengthened over 2010.0-2015.0 (Fig. 6), whereas we predict an accelerating decay up to 2025.0 for this coefficient. Faster axial dipole decay (in absolute value) is forecast with longer model fitting windows, while for the evolution  0 (a, b) or N = 2 (c, d) of g 0 7 and g 0 14 we find no discernible correlation with the length of model fitting windows used.
We have also calculated the forecast change in several global quantities describing the manifestation of reversed flux on the CMB (Fig. 10), motivated by the reduced hybrid hindcast errors in the South Atlantic (Fig. 7), where these features are particularly prevalent (Gubbins and Roberts 1987; Terra-Nova et al. 2015; Metman where S R is the combined reversed patch section of the CMB. These three quantities have all been computed with a degree-3 magnetic equator (see for more details Metman et al. 2018). For the period 2019.33-2025.0, we forecast an overall growth and intensification of reversed-flux patches (Fig. 10a, b), and predict these features to migrate towards the equator (Fig. 10c). This predicted movement contrasts with the overall poleward migration of reversed flux over the twentieth century, which has been considered the main cause for the axial | cos θ| dS, dipole decay over that period (Finlay et al. 2016a;Metman et al. 2018). As such, the forecast future decay of the axial dipole (Fig. 9), can be attributed only to the alternative mechanism of growth and intensification of the reversed field. One feature of the magnetic field on Earth's surface of broad societal interest is the location of the north magnetic dip pole, where the horizontal field is zero and the field points radially downwards. Although historically centred over the Canadian Arctic, over the last few decades the pole has accelerated towards Siberia (Mandea and Dormy 2003;Chulliat et al. 2010), crossing the international date line in 2017. The location of the pole is currently governed by the relative magnitudes of two patches of intense reversed field on the CMB centred under Canada and Siberia, with a weakening of the surface signature of the Canadian patch responsible for the recent acceleration (Livermore et al. 2020). Figure 11 shows the predicted position of the pole according to our hybrid forecast with L u = 14 , N = 2 and the 2018.33-2019.33 fitting window. Our prediction of the pole moving further towards Siberia agrees with the selection of models shown in Livermore et al. (2020), here with a time-averaged speed of approximately 46 km /yr. Finally, we present our IGRF-13 candidate model for secular variation. Using our favoured 1-year model fitting period, we fit a steady flow and subsequently the residual with magnetic diffusion ( N = 2 ) over the period 2018.33-2019.33 using Gauss coefficients up to degree 14. We then forward evolve our model to predict the Gauss coefficients in 2020 and 2025. Dividing the difference by five gives the annual SV coefficients for our candidate model over this 5-year period. Although we calculate the SV Gauss coefficients up to degree 14, only those up to degree 8 are relevant for our candidate model, which can be found in Additional file 1.

Discussion
We have presented a new, computationally inexpensive method to forecast the short-term evolution of the core magnetic field. With our framework, we combine the work of Whaler and Beggan (2015) and Metman et al. (2019), by first inverting a time-dependent field model (here chosen to be CHAOS-6-x9) for steady core fluid motion, and subsequently fitting a model of magnetic diffusion to the residual generated by this flow. Our methodology is designed specifically to improve forecasting; our separate treatment of core flow induction and magnetic diffusion prevents us from accurately modelling the underlying physics in which the effects are intertwined. It is worth noting that we have opted to first fit a core flow and then used diffusion to model the residual. In an alternative model, we could have chosen to fit diffusion first and then correct the model with core flow induction. Such a model would be dominated by diffusion, and indeed it is possible that diffusion plays a much more significant role than is traditionally believed, as recent work has shown that magnetic diffusion alone can account almost entirely for core field evolution even over short time windows (Metman et al. 2019). Nevertheless, we consider it more physically realistic to allow core flow to explain most of the SV, as a purely diffusive field evolution is not self-sustainable (e.g., Gubbins and Roberts 1987), and diffusive mechanisms always depend on fluid flow to generate the concentrations of the magnetic field that subsequently diffuse (a good example of this is flux expulsion, see e.g., Bloxham 1986). In our model, the diffusion term acts primarily as a correction for high spherical harmonic degrees. This behaviour is certainly physically reasonable as diffusion should operate most rapidly on small magnetic length scales. In terms of hindcast tests, we find limited reduction in global hindcast accuracy when the spatial complexity of the steady flow is increased beyond L u = 8 , whereas accounting for diffusion even with coarse radial resolution can result in a significant hindcast improvement on yearly timescales, in particular for features of high spherical harmonic degree. Therefore, since attenuation of the field as a function of radial distance above the CMB increases with with increasing spherical harmonic degree, the hybrid hindcasts (accounting for fluid flow and diffusion) perform comparatively better at the CMB than at Earth's surface (compare, for example, Tables 1, 2). Moreover, by correcting for diffusion we find relatively large hindcast error reductions in the South Atlantic, demonstrating that the evolution of patches of reversed flux are captured better with this hybrid scheme than a steady flow alone. Reversed-flux patches are forecast to contribute to axial dipole decay up to 2025.0 through their proliferation and intensification.
Our flow models have little effect on the prediction of field features characterised by a spherical harmonic degree higher than 8. For these degrees, the purely flowbased hindcast error spectra (Fig. 3) are independent of the spatial complexity of the flow. The part of these error spectra above degree 8 bears a striking resemblance to those of the IGRF SV hindcast, highlighting the need to introduce diffusion to outperform a simple linear extrapolation at high degree. It is possible that this behaviour of the flow models relates to the regularisation imposed in the flow inversion, which may damp small-wavelength flow to such an extent that it provides a negligible contribution to global field evolution. Further investigation is required to determine if this is a robust feature of our flow-based hindcasts. For example, the effect of different damping parameter values could be considered, as could different means of regularisation (e.g., the use of a kinetic energy norm (Madden and Le Mouël 1982)) in place of the 'strong' norm (Bloxham 1988) used here, which would reduce the penalisation of energetic flow at higher degrees. Nonetheless, the ambiguities of the recovered flow allow a small amount of diffusion to accommodate the misfit between the model and SV data, as observed in Lesur et al. (2015). The global hindcasting accuracy achieved with our hybrid method is comparable to, and occasionally better than, that reported in previous works. For example, for the model fitting periods 2001.0-2010.0, 2005.0-2010.0, and 2007.0-2010.0 our global RMS hindcast error at Earth's surface in 2015.0 is smaller than for the hindcasts presented by Whaler and Beggan (2015) and Beggan and Whaler (2018) using, respectively, time-dependent flows and ensemble Kalman filters (see, respectively, Table 2 and 1 in those works). When we compare our results with the best performing hindcast of Whaler and Beggan (2015) (i.e. that obtained with the 2007.0-2010.0 model fitting window) the relative difference in the RMS surface error at 2014.5 is slight, approximately 5%. Similarly, Bärenzung et al. (2018) report an RMS hindcast error at Earth's surface of 66 nT at the end of a 2010.0-2015.0 hindcast period, which is roughly 6% higher than our best performing hybrid hindcast (Table 1). Considering these rather small reductions in hindcast error obtained with our hybrid scheme, we emphasise once more that accounting for magnetic diffusion appears useful predominantly for hindcasting field evolution at the CMB instead of at Earth's surface, with the CMB error reductions being proportionally much larger (Table 2). Additionally, we note that we optimise for a global RMS hindcast error, and other methods may still outperform our hybrid approach regionally or locally. Hence, there are other conceivable combinations of flow and diffusion which may give a physically or mathematically better forecast model.
There is a clear positive correlation between the length of the model fitting period and the global hindcast error (Figs. 2, 4), which is likely due to the steady nature of our modelled core flows. Within a relatively long model fitting period (e.g., 2001.0-2010.0), there is likely significant time-dependence in the flow, whose signature in the SV is not reproducible with a steady flow (Waddington et al. 1995). Indeed, the reduction in the RMS surface error achieved with the hybrid scheme compared to core flow alone is relatively large for long modelling periods (Table 1), which may be attributed to diffusion allowing a more accurate description over such timescales.
A better description of short-period SV could be accomplished by including higher-order terms in eq. (12), along with observations of secular acceleration associated with fluid flow acceleration (Whaler and Beggan 2015) and/or by calculating the secular acceleration generated by diffusion (e.g., by taking a time derivative of Eq. (7)). The error of truncation associated with Eq. (12) would then be O (t k − t k−1 ) 3 (for sufficiently short time intervals); therefore, such an approach could produce forecasts of increased accuracy, although the method will still be hampered by the separate treatment of core flow and diffusion.
Nevertheless, our hybrid forecasts generally perform better than those based on fluid flow and flow acceleration alone, suggesting that accounting for diffusion is more advantageous than including more spherical harmonic degrees or higher-order terms in Eq. (12).

Conclusions
We have described a new hybrid scheme for producing forecasts of the geomagnetic field, in which core flow and diffusion are both represented. Hindcast tests showed that fitting the models over the shortest period (here, 1 year) gave the most accurate results, likely because the steady-flow assumption is less valid as the modelling windows increase. We presented our candidate model for SV over 2020-2025 for IGRF-13, and discussed the change in geomagnetic field that it forecasts.