Modeling and predicting the short term evolution of the Geomagnetic field
The coupled evolution of the magnetic field and the flow at the Earth’s core mantle boundary is modeled within the time period. To constraint the dynamical behavior of the system with a core field model deriving from direct measurements of the Earth’s magnetic field we used an Ensemble Kalman filter algorithm. By simulating an ensemble of possible states, access to the complete statistical properties of the considered fields is available. Furthermore, the method enables to provide predictions and to assess their reliability.
In this study, we could highlight the cohabitation of two distinct flow regimes. One associated with the large scale part of the eccentric gyre, which evolves slowly in time and posses a very long memory of its past, and a faster one associated with the small scale velocity field. We show that the latter can exhibit rapid variations in localized areas. The combination of the two regimes can predict quite well the decadal variations in length of day, but it can also explain the discrepancies between the physically predicted and the observed trend in these variations.
Hindcast tests demonstrate that the model is well balanced and that it can provide accurate short term predictions of a mean state and its associated uncertainties. However, magnetic field predictions are limited by the high randomization rate of the different velocity field scales, and after approximately years of forecast, no reliable information on the core field can be recovered.
BAERENZUNG, HOLSCHNEIDER, WICHT, SANCHEZ, LESUR
In the Earth’s outer core, turbulent motions of the electrically conducting fluid sustains the geomagnetic field through dynamo action. Part of this field, the poloidal one, crosses the mantle and can be observed at the Earth’s surface and above. Because of the low conductivity of the mantle (see Velímský (2010); Jault (2015)), once measured and modeled, the poloidal field can be estimated everywhere outside and at the outer boundary of the core. Since at the core mantle boundary (CMB) the flow is advecting the magnetic field, a close examination of the geomagnetic field variations at the Earth’s surface, make it possible to infer the velocity field responsible for them.
Within the outer core, the evolution of the magnetic field is prescribed by the induction equation. At the CMB this equation can be simplified. Under the assumption that the mantle is a perfect electrical insulator, the toroidal field, which interacts with the poloidal field inside the outer core, vanishes at the CMB. In addition, since the fluid cannot penetrate the mantle, its associated velocity field is purely two dimensional. And finally, on short timescales, diffusion effects can be considered as negligible in comparison to advection effects as shown in Holme (2007). So all in all, the induction equation expressed at the CMB can be simplified into the so-called Frozen Flux (FF) approximation (see Backus et al. (1996)). By inverting this equation, which couples the velocity field to the radial component of the magnetic field and the secular variation, fluid motions at the CMB can be recovered.
However, since the velocity field has two components for one equation, and since any flow scale can interact with the magnetic field to generate the large scale, observable, secular variation (SV), the inverse problem is ill-posed. To reduce the non-uniqueness of the velocity field, physical assumptions decreasing the dimension of possible solutions can be made. Most of the constraints generally used in core flow inversions, such as quasi-geostrophy, columnar, tangential-geostrophy or purely toroidal, are derived in Holme (2007); Finlay et al. (2010). Nevertheless, to obtain a unique solution, additional constraints on the velocity field need to be imposed. Typically, one enforces the energy associated with the small scale velocity field to rapidly decay, based on the so-called large scale assumption (see Holme (2007); Finlay et al. (2010)). However, Baerenzung et al. (2016) has shown that although the flow is dominant at large scales, its total kinetic energy spectrum does not exhibit a strong decaying slope. Recently, other strategies have been developed to bypass the issues raised by the nonuniqueness of the velocity field. In particular, Aubert (2014) proposed to use the statistical properties of an Earth like geodynamo simulation (the Coupled-Earth model of Aubert (2013)) to characterize the flow and the magnetic field a priori. A major advantage of such an approach is that the correlations between the fields at the CMB and the fields within the outer core are available allowing for imaging the entire outer core state.
Constraining a priori the temporal dependency of the velocity field is a more delicate operation as constraining it spatially. Optimally, one should account for the dynamics of the outer core fluid and magnetic field prescribed by the magnetohydrodynamic equations. Approaches such as variational data assimilation (see Canet et al. (2009); Li et al. (2014)) allow for physical modeling within an inversion framework. To do so, the method search for the optimal initial conditions of the system in order for the deterministic trajectories of the different fields to explain at best the observations. The drawback of the method is that all data are treated simultaneously, so whenever the dimension of data or their amount are large, the algorithm become computationally expensive. An alternative to avoid such a block inversion is to operate sequentially. Kuang and Tangborn (2008) were the first to adopt a sequential assimilation algorithm in the context of geomagnetic modeling. The optimal interpolation algorithm they used proceeds recursively in two steps. In the first one, the state variables are propagated in space and time with a given physical model, in the case of Kuang and Tangborn (2008), a three-dimensional geodynamo simulation. Once observations become available, the state variables are corrected accordingly, and then the prediction restart. Since with this approach uncertainties are not modeled, they have to be specified in an ad-hoc manner. Algorithms which permit to simulate a mean model and its associated uncertainties exist, the Kalman filter (KF) approach (see Kalman (1960); Talagrand (1997); Cohn (1997); Evensen (2003)) is one them. As for the optimal interpolation algorithm the KF proceeds sequentially with a forecast and an analysis steps. The main difference is that the evolution of errors is also predicted, and whenever data become available, these errors are taken into account for the Bayesian update of the state variables. The KF can only be applied to systems exhibiting a linear dynamics. When the dynamics of the system is nonlinear, as it is the case here for the geomagnetic field, propagation of errors cannot be analytically derived. However, it can be either approximated by linearization with the extended version of the KF, or represented through an ensemble of possible solutions in the ensemble Kalman filter (EnKF). As for Barrois et al. (2017) the latter option is the one we chose in this study.
In the EnKF, the different fields of interest are represented through an ensemble of possible states. For the prediction step, the dynamical model of the system prescribes the spatio temporal evolution of each individual member of the ensemble. At the analysis, covariances deriving from the forecasted fields and data are combined to correct the state of the ensemble predicted at observation time. Due to the limitations in available computational power, a balance between complexity of the dynamical model and size of the ensemble has to be found. Here we decided to favor accuracy in statistical representation. Therefore, the evolution of the core magnetic field and velocity field are only modeled at the level of the core-mantle boundary, through respectively the frozen flux equation and a first order auto regressive process. Extending the approach of Baerenzung et al. (2016) to the time domain, the parameters of the auto regressive process, are assumed to derive from scale dependent power laws, and are directly estimated with the COV-OBS.x1 core magnetic field secular variation model of Gillet et al. (2013, 2015a).
2 Models and parameters
2.1 Quantities of interest and notations
In this study three fields are of particular interest, the radial component of the magnetic field , the secular variation , and the velocity field all expressed at the Earth’s core mantle boundary.
The spectral counter parts of and are respectively given by the spherical harmonics coefficients and such as:
with the Schmidt semi-normalized spherical harmonics (SH) of degree and order .
The velocity field at the CMB is decomposed into a poloidal and toroidal scalar field such as:
where corresponds to the horizontal divergence operator. In spectral space, the poloidal and toroidal fields are respectively deriving from the coefficients and through the formulation:
According to the spherical harmonics expansion of , and , we define the magnetic field and secular variation energy spectra as following:
and the velocity field poloidal and toroidal energy spectra respectively as:
In the following, we will also adopt two types of notation, one with normal characters and another one with bold characters. Normal characters will be associated with single epoch quantities. This includes , and , which will respectively contain the spherical harmonics coefficients of the secular variation , the magnetic field , and both the poloidal and toroidal fields and , but also any scalar, vector or matrix associated with a quantity expressed at a given time.
Bold characters will be used for quantities depending on both space and time. As an example, a vector will contain spatial vectors at different epochs, such as:
Finally, a constant formalism for some statistical quantities will be used all over the manuscript. The mean and maximum value of a distribution will respectively be written with an over bar and a hat such as:
and the covariance associated with a random variable or between a random variable and a random variable will respectively be expressed as:
Note that the use of bold characters for space time dependent variables also applies to the latter statistical quantities.
2.2 Characterization of the magnetic field dynamical behavior
As mentioned in the introduction, under the assumption that the mantle is a perfect electrical insulator, and supposing that the observed secular variation is mainly induced by advection of the magnetic field at the CMB, the dynamical evolution for the radial component of the magnetic field at the CMB is given by the frozen flux equation which reads:
Following the notations given in section 2.1, this equation can be written in spectral space as:
where the linear operators and and the third order tensor allow us to calculate the SH coefficients associated with the advection term when they are respectively applied to , and .
2.3 Characterization of the flow spatio-temporal behavior
The dynamical evolution of the fluid within the Earth’s outer core is prescribed by the Navier-Stokes equations. Numerically solving this equation is not only numerically expensive, it remains out of reach at the regime the Earth is exhibiting. Since on short timescales, the observable magnetic field and secular variation only depend on the velocity field at the CMB, we propose to model the outer core flow only at this location. Following Gillet et al. (2015b) we chose a first order autoregressive process to do so. With such a model, the velocity field , depends on the velocity field through the relation:
where is the matrix associated with the parameters of the autoregressive model, we refer it as memory term, and is a white noise characterized by the Gaussian distribution with the covariance of the noise given a certain parametrization of the AR process. If each singular value of is smaller than , the process is stationary, a state that the flow at the CMB is certainly fulfilling. Under such conditions, the time-averaged spatial covariance of the velocity field is related to the covariance of the white noise through the relation:
Defining the vector containing the velocity fields at different epochs, its full covariance is given by:
If is a scalar, this parametrization of the flow spatio-temporal behavior is similar to the one used by Gillet et al. (2015b). Indeed, in their study they chose a temporal correlation of the velocity field depending on the characteristic time such as , with the time stepping of the auto regressive process.
2.4 Parametrization of the auto regressive process
To fully describe the autoregressive process, its two main parameters, the spatial covariance and the memory term of the process , have to be characterized. Following the developments of Baerenzung et al. (2016), is chosen to derive from the poloidal and toroidal stationary spectra of the flow, the latter being assumed to behave as power laws with different spectral ranges. , which contains the information on the temporal correlations of the velocity field, is assumed to be scale dependent, therefore we choose it to derive from power laws presenting the same ranges than the flow energy spectra. Under such assumptions, the poloidal and toroidal stationary spectra of the flow and memory terms of the AR process are given by:
where the ’s are the magnitudes of the energy spectra and the memory terms, and the ’s are their slopes within the spherical harmonics ranges ’s. The constants ’s are given by:
being the spherical harmonics degrees where transitions in slope occur.
According to this parametrization, whereas the matrix is diagonal and contains both and , the spatial covariance derives from the poloidal and toroidal energy spectra such as:
The total parameters of the autoregressive process, referred as , can be divided into two categories, one containing the parametrization of the spatial covariances, namely , and the other associated with the memory terms of the process . and are respectively given by:
where the index and are associated with the different poloidal and toroidal spherical harmonics ranges. Once is given, the prior distribution of the velocity field can be expressed and reads:
2.5 Posterior distribution the autoregressive parameters
To forecast the evolution of the velocity field in the ensemble Kalman filter algorithm, the parameters of the autoregressive process have to be known. Their posterior distribution, , can be expressed through , the joint posterior distribution of the AR parameters, the core mantle boundary velocity field and magnetic field, such as:
where is the observed secular variation. For this study, is taken from the COV-OBS.x1 model of Gillet et al. (2013, 2015a), and the covariance matrix of the secular variation , is derived from the ensemble members provided by the model. Because of the singular nature of the covariance matrix, only its diagonal part is kept. Under such conditions, the likelihood distribution reads:
The prior distribution of the magnetic field is decomposed into two parts. The first part provides the statistical properties of the large scale field , whereas the second part describes our prior knowledge on the small scale field . The large scale magnetic field is characterized by the prior distribution:
where and are respectively the COV-OBS.x1 magnetic field and covariance matrix. As for the secular variation, is evaluated with the ensemble members of the model, and only its diagonal part of kept.
The small scale magnetic field is chosen to be at any time isotropically distributed, with a mean and a covariance deriving from the extrapolation of the large scale field spectrum . Here we chose the formulation proposed by Buffett and Christensen (2007) to characterize the magnetic field spectrum at the CMB. It reads:
where . To determine the constant , we used the COV-OBS.x1 magnetic field sampled every two years between and , and performed a weighted least square fit of the associated energy spectra between SH degree and . We obtained that nT. Note that another type of extrapolation has also been tried, assuming an exponential decay of the magnetic field spectrum. Although we do not show the results associated with this modeling, we observed that such an assumption would provide insufficient levels of energy at small scales, leading to suboptimal predictions of the magnetic field evolution. From the extrapolation given in equation (34) we construct the covariance of the small scale magnetic field at a given time through the relation:
Neglecting a priori the temporal correlations between the small scales of the magnetic field, the full covariance is simply a bloc diagonal matrix where every block are identical and given by .
As expressed in the previous section, the prior distribution of the velocity field conditioned by the AR parameters is given by equation (30).
The last distribution entering equation (31), is the prior distribution of the AR parameters . These parameters are depending on the magnitudes ’s, the slopes ’s, and the spherical harmonics ranges ’s of the flow stationary spectra and the AR memory terms. Whereas the ranges ’s will be a priori imposed, and therefore considered as known, the magnitudes and slopes are completely undetermined. To reflect this lack of knowledge we characterize them by uniform distributions such as:
The full prior distribution of is simply the product of the prior distributions of each individual AR parameter.
Finally, following the development of Baerenzung et al. (2016), the posterior distribution of the AR parameters given the secular variation of equation (31) is approximated by the following distribution:
where is the dimension of the secular variation vector. To construct the matrix , the covariance between at a time and at a time , with respect to the distribution is calculated for every combination of epochs considered. The component at a row index and a column index of the resulting covariance matrix reads:
where we recall that the third order tensor is defined such as .
2.6 Sequential assimilation of the core secular variation and magnetic field
To combine our dynamical model for the magnetic field and the velocity field at the CMB to a magnetic field and secular variation model derived from geomagnetic data, we implemented the Ensemble Kalman filter approach proposed by Evensen (2003). This method proceeds in two steps. In the first one, referred as the prediction step, the spatio temporal evolution of the state variables (the velocity and magnetic fields) represented through an ensemble is forecasted until data become available. Then, the second step, called the analysis, is initiated, and each member of the ensemble is corrected accordingly to the data.
As mentioned previously, while the dynamical behavior of the velocity field is prescribed by a first order auto regressive process, the evolution of the magnetic field is constrained by the frozen flux equation. Nevertheless, directly solving the FF equation is numerically unstable. Indeed, since in this equation has no diffusion mechanism, cascading magnetic energy will have a tendency to accumulate on the smallest simulated scales, and slowly contaminate the entire field through non linear interactions. To counter this latter effect an extra hyperdiffusion term is added to the FF equation. Under such conditions the prediction step of the EnKF for each member of the ensemble of velocity and magnetic fields is given by:
where the memory term and the random noise of the AR process are scaled accordingly to the time step . The hypperdiffusivity is set to km.yr. If the magnetic field was only diffused over time with such an hyperdiffusion, in years the latter would loose , and of energy at respectively spherical harmonics degrees , and . The numerical resolution of the FF equation is performed through an Euler scheme for the first iteration, and a second order Adams-Bashforth scheme for the following ones. The time step has been set to half a year.
For the analysis step of the EnKF, both the observed magnetic field and secular variation are assimilated simultaneously. To do so, the model state has to be augmented in order to take the secular variation into account. Therefore, from each pair of the forecasted ensemble, a prediction for the observables is build accordingly to the following relations:
where the linear operator simply truncates the forecasted magnetic field at the level of the observed one. Normally equation (42) should contain the hyperdiffusion term of equation (41). However, the effects of this hyperdiffusion are so weak on the range of scales describing the observed secular variation (the spectral expansion of the COV-OBS.x1 SV does not excess SH degree in this study), that they are neglected.
From the augmented ensemble the covariances necessary for the analysis step of the EnKF are calculated. Recalling that the covariance of a field is referred as and the covariance between a field and a field is expressed as , each updated pair of velocity field and magnetic field is given by:
where and are random realizations from the distributions of the COV-OBS.x1 model for respectively the secular variation and the magnetic field .
3 Geophysical application
3.1 Numerical set up
The time period considered in this study is , and the magnetic field and secular variation data are taken from the COV-OBS.x1 model with a year sampling rate, corresponding to the knots of the model’s B-spline expansion. Every simulation is performed through a pseudo spectral approach on the Gauss-Legendre grid provided by Schaeffer (2013). Both the poloidal and toroidal parts of the velocity field are expanded up to SH degree , and the radial component of the magnetic field is expressed up to SH degree in order for the field to possess a large enough diffusion range. Whereas the COV-OBS.x1 magnetic field is always taken up to SH degree , the expansion of the COV-OBS.x1 secular variation depends on the variance level associated with each scale. If globally at a certain scale the standard deviation of the secular variation is larger than the absolute value of the mean field, the total field is truncated at this scale. Under such a condition, the COV-OBS.x1 secular variation is taken up SH degrees , , and for the respective time windows , , and . Finally, the state of the system is characterized by pairs of magnetic field and velocity field at the core mantle boundary.
3.2 Estimation of the flow optimal auto regressive parameters
To simulate the spatio temporal evolution of the flow at the CMB, the parameters of the auto regressive process have to be estimated. We recall that their posterior distribution, , is expressed in equation (38). By maximizing this distribution, one should get the optimal parameters for the AR process. However, instead of estimating both temporal and spatial parameters simultaneously, we proceed in two steps. First, only the spatial covariance of the velocity field is evaluated following the method proposed and tested by Baerenzung et al. (2016). This approach consists in maximizing the distribution in which only the block diagonal part of the covariance matrix is kept. Once the spatial covariance is determined, it is assumed to be known, and the maximum of is calculated.
Following the developments of section 2.4, the AR parameters are decomposed into scale dependent power laws exhibiting different spectral ranges. Baerenzung et al. (2016) showed that if the stationary spectra of the flow are decomposed into two spectral ranges, the optimal scales where transition in slope occurs are and for respectively the toroidal and poloidal energy spectra. Here, whereas we keep the same decomposition for the spectrum associated with the poloidal field, more degrees of freedom are allowed for the toroidal field spectrum. Since we wish to accurately determine the spatio-temporal evolution of the eccentric gyre, toroidal field component at SH degree and , the main components of the gyre, are free to exhibit any variance level and characteristic time. Similarly, toroidal field components at SH degree are also assumed to be unconstrained by surrounding velocity field scales. This choice is motivated by the particular low level of energy that these scales are exhibiting over recent epochs (see Baerenzung et al. (2016); Whaler et al. (2016)). Finally, one spectral range is used to characterize the toroidal field spatial variance and memory effects between SH degrees and . So all in all, the AR parameters associated with the toroidal field exhibit the four respective spectral ranges, , , and .
As mentioned in the beginning of the section, the estimation of the stationary energy spectra parameters is performed between and , taking the COV-OBS.x1 magnetic field and secular variation every years. On figure 1 the resulting power law spectra are displayed with crosses. As already observed in Baerenzung et al. (2016), the toroidal field (in black), and in particular its large scales (SH degree and ), exhibits a much larger energetic level than the poloidal field (in gray). Nevertheless the toroidal energy spectrum also presents a strong increase of energy towards its smallest scales. This effect, which is in contradiction with the results of Baerenzung et al. (2016), is very likely to be attributed to a slight underestimation of the COV-OBS.x1 secular variation uncertainties. Therefore, in order to better estimate the small scale energy spectra of the velocity field, we performed different estimations of the stationary spectra parameters by varying the time window in which the evaluation is computed. We found that the largest period where the spectra did not exhibit an anomalous behavior was . The resulting prior spectra are shown in figure 1 with circles. Combining the small scale spectra of the evaluation to the large scale ones of the estimation, we get the final prior spectra for both the toroidal and poloidal field displayed in figure 1 with solid lines. The values associated with the spectra parameters are given in table LABEL:spectraTable.
The spatial covariances of the AR process being characterized, the evaluation of the memory terms is now performed by maximizing the distribution within the time window. The results, expressed through the scale dependent characteristic time are shown on figure 2 and the parameters of the memory terms are given in table LABEL:tableParameters. The most striking feature one can observe, is the very long memory time (of the order of thousand years) associated with the main component of the eccentric gyre. This indicates that this structure is very persistent over time. Nevertheless, these values should be taken with care since their evaluation is performed on a comparatively short time window of years. In contrast with the large scale field, the toroidal field components at spherical harmonics degree varying from to exhibit lower characteristic memory times, with a decaying behavior where yr and yr. Note that this limiting time is similar to the e-folding time of the Geodynamo as calculated by Hulot et al. (2010); Lhuillier et al. (2011). The characteristic times associated with the poloidal field indicate that the latter presents a slower dynamical behavior at large scales than at small scales with varying from yr to yr.
3.3 General properties of the flow at the core mantle boundary
The AR process for the flow being parametrized, estimation of the velocity field and magnetic field at the CMB through the EnKF algorithm can be performed. To initialize the fields in , we applied the Gibbs sampling algorithm proposed by Baerenzung et al. (2016) and which is detailed in the appendix A. However, instead of sampling the joint posterior distribution of the flow and the magnetic field at the epoch only, we sampled the distribution characterizing simultaneously the fields in , and , in order to constraint the initial state with recent observations.
Because of the sequential nature of the EnKF algorithm, accuracy of the estimated fields is not constant over time, but increases whenever new data are assimilated. This effect is well illustrated by figure 3 displaying at two different epochs, (grey) and (black), the energy spectra of the toroidal (top) and poloidal (bottom) mean velocity fields (thick solid lines) and uncertainty fields (circles). Whereas in the mean toroidal field exhibit a level of energy larger or of the same order than the variance of the field up to SH degree , in reliable information become available up to SH degree . Above these scales the posterior variance of the flow rapidly reaches its prior level. One can also notice that the larger the scale, the stronger the variance reduction over time. The latter observation, which is also valid for the poloidal field, is linked to the behavior of the characteristic times associated with the different flow scales (see figure 2). The fact that is a strictly decaying function of , implies that small scales velocity field will exhibit a higher randomization rate than large scales, and so between two analysis step, the prior variance will increase faster at small than at large scales.
In physical space, the gain of flow accuracy over time is particularly stricking for the toroidal part of the velocity field as shown on figure 4. In this figure, the toroidal (left) and poloidal (right) mean velocity fields are displayed with black arrows for three different epochs, (top), (middle) and (bottom). Color maps, representing the confidence interval on the velocity field orientation, provide information on locations where the mean flow direction can be reliably estimated (violet and blue) or not (red). In , very little parts of the eccentric gyre can be confidently estimated. Only the westward flow below Africa and the Atlantic ocean, the Southern branches of the gyre, and the Northern circulation around and partially inside the tangent cylinder (the cylinder tangent to the inner core and aligned with the axis of rotation of the Earth), appear as reliable patterns. At later times the gyre is well defined, and many of its small scale structures become visible. Globally, uncertainties on the toroidal part of the flow are decreasing with time. This does not seem to be the case for the poloidal field where in reliable patterns are covering a larger surface of the CMB than in . Nevertheless, the r.m.s velocity of the poloidal field and associated standard deviation of in and in , indicate that the global uncertainty level of the poloidal field as well as its magnitude have been decreasing between and .
Contrary to flow models constrained to be geostrophic (see Bloxham and Jackson (1991); Amit and Pais (2013)), up and downwelling fluid motions, which converts into poloidal field when meeting the CMB, are not particularly located around the equator. Instead, a strong and persistent poloidal structure evolves below the Indian ocean and South Africa. According to Bloxham (1986), such a poloidal field could be at the origin, through the expulsion of magnetic flux from the outer core, of the intense reversed flux patch located there. Although the frozen flux equation cannot model the transport of magnetic structures between the core and its outer boundary, poloidal field spreading or concentrating magnetic patches around location of intense up or downwelling motions can be detected.
Other specific features of the velocity field are in apparent contradiction with a possible geostrophic state of the outer core flow. This includes the reliable part of the toroidal field penetrating the tangent cylinder, or its component crossing the equator below India and South America as already reported in other studies (see Barrois et al. (2017)). Particularly striking, is the apparent violation by the eccentric gyre of the equatorial symmetry condition imposed by quasi-geostrophy (see Amit and Olson (2004)). Indeed, the flow responsible for the westward drift together with the circulations around the tangent cylinder, exhibit different levels of intensity in the Northern hemisphere than in the Southern one. These visual observations are confirmed in figure 5, where the r.m.s. velocity and associated standard deviation of the toroidal part of the flow, are measured in different locations of the CMB. One can observe on the top right of figure 5, that the flow evolving below Africa and the Atlantic ocean within the areas shown on bottom right of the figure, is at any time more intense in the South than in the North. For the circulations around the tangent cylinder, both southern and northern part exhibit similar levels of energy between and as shown on the bottom left of figure 5. However, in the flow below Alaska and the Eastern part of Siberia starts to accelerate, and intensifies almost continuously over the last decades to reach a r.m.s. velocity around km.yr in . This acceleration has already been observed during the satellite era by Livermore et al. (2017). Here we observe that it takes its origin at quite early times. Although the toroidal part of the flow exhibits some clear deviation from geostrophy, comparisons of its r.m.s. velocity over the entire CMB (black symbols on the top left of figure 5) with the r.m.s. velocity of its equatorial symmetric part (gray symbols), shows that the latter component remains at any time dominant.
3.4 Variations in length of day
Many factors influence the Earth’s rotation and thus the length of the day (LOD) . Variations on short timescales are mostly caused by an exchange of angular momentum between Earth’s solid body and atmospheric and oceanic currents. Jault et al. (1988) have shown that the coupling between core and mantle induced decadal LOD changes. Nevertheless, modification of the Earth’s oblateness () through the melting of glacier and ice caps and through the global sea level rise is also thought to have a non negligible influence on the LOD over such timescales (see Munk (2002)). Over centenial and millennial timescales, two principle mechanisms are modifying the rotation rate of the Earth. The first one is the tidal friction. By deforming the Earth’s surface, tidal forces are inducing a dissipation of energy in the Earth-Moon system. As a consequence, the rotation rate of the Earth is decreasing, and the LOD increases with a rate of ms/cy as estimated by Williams and Boggs (2016). The other phenomenon modifying the long term variations in LOD is the glacial isostatic adjustment (GIA). When accumulating on the polar caps over the last glaciation period, ice was compressing the mantle, inducing an increase of the Earth’s oblateness. When the ice melted, the mantle tend to regain its initial shape at a rate depending on its viscosity profile. Peltier (2015) has shown that the associated decrease in oblateness would decrease the LOD at a rate of ms/cy.
The sum of the latter effects ms/cy explains very well the millennial trend of increase in LOD of ms/cy observed by Stephenson et al. (2016). Nevertheless, if one accounts for the influence of the global sea level rise and ice melting over the last century (see Church and White (2011); Hay et al. (2015)), the trend in the observed variations in LOD () should be slightly larger than ms/cy. More precisely, through measurements of the Earth’s oblateness with satellite laser ranging, Cheng et al. (2013) have shown that the oblateness was decreasing before the 1990’s, but not sufficiently to be explained by the GIA alone, and increasing after. Therefore, the trend in should lie between ms/cy and ms/cy before the 1990’s and it should be larger than ms/cy after.
As highlighted by Munk (2002), these estimations are in contradiction with the trend of ms/cy observed over the almost last years. Therefore, either a mechanism increasing the rotation rate of the Earth is missing, or GIA effects are underestimated. Mitrovica et al. (2015) showed that combining a GIA model exhibiting a lower rate than the one of Peltier (2015), to lower estimates of the global sea level rise over the last century, could explain the discrepancies between observed and predicted trend in . Nevertheless, for this model to be compatible with the millennial trend, a mechanism slowing down the rotation of the Earth over very long periods of time is necessary. They suggested that the outer core flow estimated by Dumberry and Bloxham (2006) could be responsible for this. This means that the oscillations of the core angular momentum highlighted by Dumberry and Bloxham (2006) would also be accompanied by a global decay.
So because of the uncertain nature of the trend in , we estimated the optimal value it should take accordingly to our ensembles of velocity fields. To compute the variations in LOD deriving from the flow at the core mantle boundary, we used the formula of Jault (2015) which reads:
Based on the previous arguments, we decompose the LOD into the core contribution and an additional long-term linear trend with rate modeling the effects discussed above such as:
where is related to a reference observation at : . Instead of prescribing and we search for their optimal values accordingly to our core flow model in a Bayesian approach. Assuming that and are a priori unknown (uniform prior over infinite ranges), the posterior distribution of can be expressed as:
Because of the abrupt change in the Earth’s oblateness during the ’s, we restrict the analysis to pre epochs. Observed LOD variations are taken from Gross (2001) who also provides uncertainty estimates . The likelihood distribution is approximated by a Gaussian distribution such as
The prior distribution of is also assumed to be Gaussian with a mean and a covariance deriving from the ensemble of time series calculated with equation (45). The distribution is thus also a Gaussian distribution. Our computation suggest a mean of ms/cy and a large standard deviation of ms/cy which embraces the values discussed above.
The mean estimated trend lies well between the ms/cy and ms/cy interval prescribed by tidal forces and by the measurements of the Earth’s oblateness. Furthermore, once the optimal trend is removed from , comparisons with the core flow contribution exhibit a good agreement on decadal time scales as illustrated by figure 6. Although the large uncertainty levels associated with and , forbids any precise conclusions on the impact of the last century sea level rise and ice melting on the variations in length of day, our results suggest that the scenario proposed by Mitrovica et al. (2015) is not the most likely (otherwise the optimal trend would be close to ms/cy). Instead, estimates of given by Peltier (2015) are probably more appropriate. This means that discrepancies between predicted and observed trend in over the last centuries could be compensated by a recent increase in core angular momentum, and that the outer core flow would not slow down the rotation of the Earth over very long periods of time.
The ability of a model to successfully predict the system evolution not only suggests that the model correctly captures the dynamics on the considered time scales but also points towards useful applications. We test our model with so-called hindcast simulations, which means that the analysis steps in the assimilation are only carried out until a time , and then integrated as a free model run up to a time . The free run corresponds to the forecast, or prediction, which can be compared with the data.
Here we use six different values, , , , , , and , and compare predictions for with the respective epoch in the CHAOS-6 magnetic field model by Finlay et al. (2016). This means that we attempt predictions over periods ranging from to years. Figure 7 presents the results in terms of energy spectra at the Earth’s core mantle boundary. Thick black lines show the CHAOS-6 reference field while thin gray lines show the predicted ensemble mean . Crosses illustrate the prediction error which can be compared with thick gray line that depicts the predicted error, which is simply the standard deviation of the ensemble members. As expected, the errors increase with the prediction period (decreasing ). The predicted error provides a good estimate for the prediction error with the exception of degree in the year period. This likely owed to the fact that the CHAOS-6 core field contains the crustal field that starts to contribute more significantly at this scale. The larger prediction error outweighs this deficiency for longer prediction periods.
The comparison with two other more trivial prediction methods further allows to judge the advantage of our more sophisticated approach. ’No cast’ refers to the assumption that the field remains identical to the field at . ’Linear extrapolation’ uses the secular variation at to linearly extrapolate the field from to : . The prediction errors of these two trivial methods are shown as triangles and circles in figure 7. Linear extrapolation and assimilation prediction errors remain similar for the two shortest prediction periods. However, for predictions beyond the yr horizon, our assimilation formalism starts to pay off. For the longest prediction period of years, the errors in the two trivial methods already exceed the spectral energy at degree while the assimilation predictions remains appropriate until degree .
Working with an ensemble also allows the evaluation of statistical properties of quantities like inclination or declination that depend nonlinearly on the state variables. Figure 8 compares the inclination (left) and declination (right) for (black lines) and (red lines) with the ensemble mean predictions using (yellow lines). The prediction errors are quantified by the absolute local difference in degrees and are shown as color maps in the top two panels. Color maps in the bottom two panels show the confidence interval of the ensemble prediction.
As expected, inclination errors are large for the small or vanishing values around the equator. This is illustrated in the top panel of the right column and also well captured by the larger variance in the ensemble used to predict the error shown in the lower right panel. We calculated the area where the prediction error remains within the % confidence interval defined by the prediction ensemble, which amounts to % of the total surface. This indicates that the inclination uncertainties provide a good estimate of the prediction error.
Since the declination becomes undefined at the North and South dip poles, the surrounding regions show larger prediction errors and ensemble variances, as is demonstrated in the left column of figure 8. The area where prediction errors remain within the confidence interval amounts to only of the surface, which suggests a somewhat inferior error prediction likely related to the dip poles. The error prediction seems more reliable in the southern hemisphere alone where the relative area increases to .
Although not clearly visible in figure 8 because the color has been saturated at , the prediction error for the location of the North Magnetic Dip Pole (NMDP) is also particularly large. Its real and predicted position for 2015 have been marked by gray circles and triangles, respectively. According to Chulliat et al. (2010), the rapid acceleration of the NMDP drift during the 1990’s can be explained by the expulsion of magnetic flux below the New Siberian Islands. This may be the reason for the larger prediction error, since such expulsions are not modeled in the frozen flux approximation used here. Note that the NMDP drift is much better captured when the hindcast test starts at .
Finally, we note that our model accurately predicts the evolution of the inclination and declination associated with South Atlantic anomaly despite the significant changes between and . This highlight the potential usefulness of the method for forecasting core’s magnetic field features.
In order to quantify the different sources of forecast errors we analyze the secular variation which is responsible for advancing the field from to : . Under the frozen flux approximation, depends on the flow and the magnetic field, through the relation . When separating into the observable part and the non observable small scale part we can distinguish three error sources:
Here the primed quantities on the right hand side indicate deviations from the ensemble expectation value, for example .
Figure 9 compares the spectra of the different error contributions to the secular variation spectrum at the Earth’s surface for (left) and (right). The last analysis step performed at directly constraints and leads to a small related variance and thus a small error contribution that lies about two orders of magnitude below the and related error. At the end of the forecast at , however, the related error has grown significantly but remains the smallest contribution. Since is mostly constrained by the apriorily assumed statistics at all times, the related error changes only mildly. While the related error is smaller than the error at , it becomes the dominant contribution at due to the increase in flow dispersion.
Thus neither the observational error in the large scale magnetic field nor the lack of knowledge on the small scale contributions is the limiting factor for the predictions but rather randomization of the different velocity field scales over time. We also tested the predictability range for the magnetic field. Starting a forecast in , we let the system evolve until the scale per scale variance of the magnetic field exceeds the mean predicted energy. Like in the forecast from 1990 to 2015 illustrated in the middle left panel of figure 7, after yrs only, the variance in contributions beyond degree exceed the mean energy. The predictability limit further decreases to , , and after yrs, yrs and yrs repsectively. After yrs, even the dipole energy is exceeded by the respective variance level.
We have employed a sequential data assimilation framework to model the dynamics of the geomagnetic field and the flow at the top of Earth’s core in the 20th century using the COV-OBS.x1 model of Gillet et al. (2015b) as observations. The method extends the approach in Baerenzung et al. (2016) to the time domain, as a sequential propagation in time of flow and field inversions under weak prior constraints. The prior is a dynamical model that combines the induction equation in the frozen flux approximation with a simple AR1 process describing the flow evolution. The latter comprises a memory term and stochastic forcing which are both constrained by the secular variation observations following the ideas presented in Baerenzung et al. (2016).
We use an ensemble approach, the EnKF (Evensen (2003)), where the dynamical model uncertainties are characterized by statistically sound covariances. Using the AR1 process falls short of integrating a proper Navier-Stokes equation but allows us to forward a large ensemble of 40,000 members in time in order to characterize the errors and prior covariances.
The optimal parameters characterizing the flow prior spatial and temporal properties points to particularly long time scales of several centuries to millennia for the toroidal field at SH degrees and . This flow contribution can be attributed to a large scale slowly evolving gyre which has also been identified in core flow inversions (Pais and Hulot (2000); Gillet et al. (2015b)) and numerical simulations (Aubert et al. (2013); Schaeffer et al. (2017)). The most prominent feature of the gyre is the well documented pronounced westward drift at low and mid-latitudes of the Atlantic hemisphere. Smaller scale flows have characteristic time scales in the decadal range that are consistent with previous estimates (Christensen and Tilgner (2004); Hulot et al. (2010)). Typical related features are local modifications of the gyre in the southern hemisphere or the acceleration of the westward flow at and around the tangent cylinder underneath Alaska and the eastern part of Siberia, which has already been reported by Livermore et al. (2017). This points to an important contribution of ageostrophic motions to the dominantly geostrophic overall core flow.
Predictions of decadal length of day variations (LOD) from changes in angular momentum of our core flow model yield remarkable resemblance to corresponding independent observations. Furthermore, the recent increase in core angular momentum, that we attribute to an acceleration of the geostrophic contribution of the gyre, enables to compensate the difference between the recently observed trend in LOD changes and the expected one highlighted by Munk (2002).
We further tested the capability of our model to forecast the evolution of magnetic and flow fields through hindcast experiments. Comparisons of the magnetic field evolution with linear extrapolations and no casts (in which the field is assumed static) show that our more sophisticated model significantly improves predictions beyond 10 or 15 years. Moreover, we inspected the reliability of the forecast errors predicted by the ensemble dispersion, which resulted in good agreement with the hindcast errors. Such match reveals that the characteristic times estimated for the auto regressive process of the flow correspond to a realistic randomization rate of the fields.
However, it is the dispersion of the velocity field itself which seems to dominate the uncertainties in secular variation estimations, what limits therefore the predictability of the geomagnetic field. Within this limitation, the scale dependent predictability corresponds to 1950, 640, 400, 160 and 75 years for degrees , 2, 3, 5 and 7, respectively. A more realistic dynamical model such as a geodynamo simulation would possibly extend the predictability horizon. Nevertheless, the enormous numerical power required to perform dynamo simulations at more extreme parameters would preclude the type of ensemble Kalman filter approach followed here. Even with compromises in dynamo model parameters and the ensemble size, the computational costs would still increase by orders of magnitude. Moreover, since dynamo simulations are strongly nonlinear, the system bears an intrinsic sensitivity to initial perturbations. This amounts to an important e-folding time, which is estimated to be about 30 years with a temporal rescaling within the secular variation time scale (Hulot et al. (2010); Lhuillier et al. (2011)). Since this characteristic time is not so different from the one associated with the small length scales of our flow model, we expect that the rate at which information is lost in the system will be somewhat equivalent in both modeling strategies.
Appendix A Gibbs sampling
The Gibbs sampling algorithm permits to randomly draw an ensemble characterizing statistically a given joint distribution, by recursively sampling conditional probability distributions deriving from it. In our case, the joint distribution of interest is . Therefore, at a step , the algorithm samples alternatively the two following distributions:
In this study, three epochs were considered simultaneously, , and . Since over long periods of time the real temporal correlations of the slow varying components of the flow probably differ from the one induced by the auto regressive process, the characteristic times associated with the degree and of the toroidal field and of the poloidal field were reestimated with a year time step.
Acknowledgements.This work has been supported by the German Research Foundation (DFG) within the Priority Program SPP1788 “Dynamic Earth”. The data for this paper are available by contacting the corresponding author at email@example.com.
- affiliationtext: Institute for Mathematics, University of Potsdam, Potsdam, Germany
- affiliationtext: Max Planck Institute for solar system research, Göttingen, Germany
- affiliationtext: Institut de Physique du Globe de Paris, Paris, France
- Amit, H., and P. Olson (2004), Helical core flow from geomagnetic secular variation, Physics of the Earth and Planetary Interiors, 147, 1–25, doi:10.1016/j.pepi.2004.02.006.
- Amit, H., and M. A. Pais (2013), Differences between tangential geostrophy and columnar flow, Geophysical Journal International, 194, 145–157, doi:10.1093/gji/ggt077.
- Aubert, J. (2013), Flow throughout the Earth’s core inverted from geomagnetic observations and numerical dynamo models, Geophysical Journal International, 192, 537–556, doi:10.1093/gji/ggs051.
- Aubert, J. (2014), Earth’s core internal dynamics 1840-2010 imaged by inverse geodynamo modelling, Geophysical Journal International, 197, 1321–1334, doi:10.1093/gji/ggu064.
- Aubert, J., C. Christopher, and A. Fournier (2013), Bottom-up control of geomagnetic secular variation by the Earth/’s inner core, Nature, 502, 219–223, doi:10.1038/nature12574.
- Backus, G., R. Parker, and C. Constable (1996), Foundations of Geomagnetism.
- Baerenzung, J., M. Holschneider, and V. Lesur (2016), The flow at the earth’s core-mantle boundary under weak prior constraints, Journal of Geophysical Research: Solid Earth, 121(3), 1343–1364, doi:10.1002/2015JB012464, 2015JB012464.
- Barrois, O., N. Gillet, and J. Aubert (2017), Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics, Geophysical Journal International, 211(1), 50–68, doi:10.1093/gji/ggx280.
- Bloxham, J. (1986), The expulsion of magnetic flux from the Earth’s core., Geophysical Journal, 87, 669–678, doi:10.1111/j.1365-246X.1986.tb06643.x.
- Bloxham, J., and A. Jackson (1991), Fluid flow near the surface of earth’s outer core, Reviews of Geophysics, 29, 97–120, doi:10.1029/90RG02470.
- Buffett, B. A., and U. R. Christensen (2007), Magnetic and viscous coupling at the core-mantle boundary: inferences from observations of the Earth’s nutations, Geophysical Journal International, 171, 145–152, doi:10.1111/j.1365-246X.2007.03543.x.
- Canet, E., A. Fournier, and D. Jault (2009), Forward and adjoint quasi-geostrophic models of the geomagnetic secular variation, Journal of Geophysical Research (Solid Earth), 114(13), B11101, doi:10.1029/2008JB006189.
- Cheng, M., B. D. Tapley, and J. C. Ries (2013), Deceleration in the Earth’s oblateness, Journal of Geophysical Research (Solid Earth), 118, 740–747, doi:10.1002/jgrb.50058.
- Christensen, U. R., and A. Tilgner (2004), Power requirement of the geodynamo from ohmic losses in numerical and laboratory dynamos, nature, 429, 169–171, doi:10.1038/nature02508.
- Chulliat, A., G. Hulot, and L. R. Newitt (2010), Magnetic flux expulsion from the core as a possible cause of the unusually large acceleration of the north magnetic pole during the 1990s, Journal of Geophysical Research (Solid Earth), 115, B07101, doi:10.1029/2009JB007143.
- Church, J. A., and N. J. White (2011), Sea-Level Rise from the Late 19th to the Early 21st Century, Surveys in Geophysics, 32, 585–602, doi:10.1007/s10712-011-9119-1.
- Cohn, S. E. (1997), An Introduction to Estimation Theory, Journal of the Meteorological Society of Japan. Ser. II, 75, 257–288, doi:10.2151/jmsj1965.75.1B257.
- Dumberry, M., and J. Bloxham (2006), Azimuthal flows in the Earth’s core and changes in length of day at millennial timescales, Geophysical Journal International, 165, 32–46, doi:10.1111/j.1365-246X.2006.02903.x.
- Evensen, G. (2003), The Ensemble Kalman Filter: theoretical formulation and practical implementation, Ocean Dynamics, 53, 343–367, doi:10.1007/s10236-003-0036-9.
- Finlay, C. C., M. Dumberry, A. Chulliat, and M. A. Pais (2010), Short Timescale Core Dynamics: Theory and Observations, Space Science Reviews, 155, 177–218, doi:10.1007/s11214-010-9691-6.
- Finlay, C. C., N. Olsen, S. Kotsiaros, N. Gillet, and L. Tøffner-Clausen (2016), Recent geomagnetic secular variation from Swarm and ground observatories as estimated in the CHAOS-6 geomagnetic field model, Earth, Planets, and Space, 68, 112, doi:10.1186/s40623-016-0486-1.
- Gillet, N., D. Jault, C. C. Finlay, and N. Olsen (2013), Stochastic modeling of the Earth’s magnetic field: Inversion for covariances over the observatory era, Geochemistry, Geophysics, Geosystems, 14, 766–786, doi:10.1002/ggge.20041.
- Gillet, N., O. Barrois, and C. C. Finlay (2015a), Stochastic forecasting of the geomagnetic field from the COV-OBS.x1 geomagnetic field model, and candidate models for IGRF-12, Earth, Planets, and Space, 67, 71, doi:10.1186/s40623-015-0225-z.
- Gillet, N., D. Jault, and C. C. Finlay (2015b), Planetary gyre, time-dependent eddies, torsional waves, and equatorial jets at the Earth’s core surface, Journal of Geophysical Research (Solid Earth), 120, 3991–4013, doi:10.1002/2014JB011786.
- Gross, R. S. (2001), A combined length-of-day series spanning 1832â1997: Lunar97, Physics of the Earth and Planetary Interiors, 123(1), 65 – 76, doi:http://dx.doi.org/10.1016/S0031-9201(00)00217-X.
- Hay, C. C., E. Morrow, R. E. Kopp, and J. X. Mitrovica (2015), Probabilistic reanalysis of twentieth-century sea-level rise, nature, 517, 481–484, doi:10.1038/nature14093.
- Holme, R. (2007), Large-Scale Flow in the Core, in Treatise on Geophysics, edited by E. in Chief:Gerald Schubert, pp. 107 – 130, Elsevier, Amsterdam, doi:http://dx.doi.org/10.1016/B978-044452748-6.00127-9.
- Hulot, G., F. Lhuillier, and J. Aubert (2010), Earth’s dynamo limit of predictability, \grl, 37, L06305, doi:10.1029/2009GL041869.
- Jault, D. (2015), Illuminating the electrical conductivity of the lowermost mantle from below, Geophysical Journal International, 202, 482–496, doi:10.1093/gji/ggv152.
- Jault, D., C. Gire, and J. L. Le Mouel (1988), Westward drift, core motions and exchanges of angular momentum between core and mantle, Nature, 333, 353–356, doi:10.1038/333353a0.
- Kalman, R. E. (1960), A New Approach to Linear Filtering and Prediction Problems, Journal of Basic Engineering, 82, 35–45, doi:10.1115/1.3662552.
- Kuang, W., and A. Tangborn (2008), MoSSTDAS: The first generation geomagnetic data assimilation framework, Communications in Computational Physics, 3, 85–108.
- Lhuillier, F., J. Aubert, and G. Hulot (2011), Earth’s dynamo limit of predictability controlled by magnetic dissipation, Geophysical Journal International, 186, 492–508, doi:10.1111/j.1365-246X.2011.05081.x.
- Li, K., A. Jackson, and P. W. Livermore (2014), Variational data assimilation for a forced, inertia-free magnetohydrodynamic dynamo model, Geophysical Journal International, 199, 1662–1676, doi:10.1093/gji/ggu260.
- Livermore, P. W., R. Hollerbach, and C. C. Finlay (2017), An accelerating high-latitude jet in Earth’s core, Nature Geoscience, 10, 62–68, doi:10.1038/ngeo2859.
- Mitrovica, J., C. C. Hay, E. Morrow, R. E. Kopp, M. Dumberry, and S. Stanley (2015), Reconciling past changes in Earthâs rotation with 20th century global sea-level rise: Resolving Munkâs enigma, Science Advances, 1(11), doi:10.1126/sciadv.1500679.
- Munk, W. (2002), Twentieth century sea level: An enigma, Proceedings of the National Academy of Science, 99, 6550–6555, doi:10.1073/pnas.092704599.
- Pais, A., and G. Hulot (2000), Length of day decade variations, torsional oscillations and inner core superrotation: evidence from recovered core surface zonal flows, Physics of the Earth and Planetary Interiors, 118, 291–316, doi:10.1016/S0031-9201(99)00161-2.
- Peltier, W. (2015), The history of the earthâs rotation: impacts of deep earth physicsand surface climate variability., Treatise on Geophysics, second ed., 9, 221â279.
- Schaeffer, N. (2013), Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations, Geochemistry, Geophysics, Geosystems, 14, 751–758, doi:10.1002/ggge.20071.
- Schaeffer, N., D. Jault, H.-C. Nataf, and A. Fournier (2017), Turbulent geodynamo simulations: a leap towards Earth’s core, Geophysical Journal International, 211, 1–29, doi:10.1093/gji/ggx265.
- Stephenson, F. R., L. V. Morrison, and C. Y. Hohenkerk (2016), Measurement of the Earth’s rotation: 720 BC to AD 2015, Proceedings of the Royal Society of London Series A, 472, 20160404, doi:10.1098/rspa.2016.0404.
- Talagrand, O. (1997), Assimilation of Observations, an Introduction, Journal of the Meteorological Society of Japan. Ser. II, 75, 191–209, doi:10.2151/jmsj1965.75.1B191.
- Velímský, J. (2010), Electrical conductivity in the lower mantle: Constraints from CHAMP satellite data by time-domain EM induction modelling, Physics of the Earth and Planetary Interiors, 180, 111–117, doi:10.1016/j.pepi.2010.02.007.
- Whaler, K. A., N. Olsen, and C. C. Finlay (2016), Decadal variability in core surface flows deduced from geomagnetic observatory monthly means, Geophysical Journal International, 207, 228–243, doi:10.1093/gji/ggw268.
- Williams, J. G., and D. H. Boggs (2016), Secular tidal changes in lunar orbit and Earth rotation, Celestial Mechanics and Dynamical Astronomy, 126, 89–129, doi:10.1007/s10569-016-9702-3.