Correlation Based Time Evolution of the Archeomagnetic Field

In a previous study, a new snapshot modeling concept for the archeomagnetic field was introduced (Mauerberger et al., 2020, https://doi.org/10.1093/gji/ggaa336). By assuming a Gaussian process for the geomagnetic potential, a correlation‐based algorithm was presented, which incorporates a closed‐form spatial correlation function. This work extends the suggested modeling strategy to the temporal domain. A space‐time correlation kernel is constructed from the tensor product of the closed‐form spatial correlation kernel with a squared exponential kernel in time. Dating uncertainties are incorporated into the modeling concept using a noisy input Gaussian process. All but one modeling hyperparameters are marginalized, to reduce their influence on the outcome and to translate their variability to the posterior variance. The resulting distribution incorporates uncertainties related to dating, measurement and modeling process. Results from application to archeomagnetic data show less variation in the dipole than comparable models, but are in general agreement with previous findings.

in the low degrees, and to exploit the data to its fullest, we suggest a Bayesian modeling approach based on Gaussian processes (GPs), both in space and time. With this already in mind, we implemented a closed-form covariance function for the spatial domain in a previous study , hereafter referred to as MSKH20). The present work extends this study to the temporal domain. We again employ the closedform correlation kernel, introduced by Holschneider et al. (2016), and extend it to a space-time kernel using a squared exponential (SQE) kernel. Knowing that such a kernel is unphysical, we abstain from suggesting a new geomagnetic field model. Instead, the aim of this work is to show the potential of the proposed modeling approach and to lay out the modeling strategy and its implications in detail.
The application of a spatiotemporal GP in a Bayesian framework includes the natural availability of well quantified uncertainties via the posterior standard deviation. While early models (C. G. Constable et al., 2000;Jackson et al., 2000;Korte & Constable, 2003) do not provide uncertainty estimates, more recent field models use ensemble techniques to quantify (modeling related) errors (Hellio & Gillet, 2018;Korte et al., 2009;Licht et al., 2013;Pavón-Carrasco et al., 2014;Senftleben, 2019). Within the space-time correlation framework that we suggest, uncertainties arising from the uneven data distribution, from inaccurate dating and from the modeling itself can be accounted for in a well defined, statistically sound manner (McHutchon & Rasmussen, 2011;Rasmussen & Williams, 2006). The inversion scheme is embedded in a functional analytic frame of non-parametric modeling. The result is a distribution over functions, in this case in both space and time. This distribution is characterized by a mean function, which gives the most likely field model, and a two-point covariance function, describing the variability of the field. This article is structured as follows: The rest of this section covers some basic introduction into magnetic field theory and GP inversion. We use those paragraphs mainly to introduce our notation. In Section 2 we discuss our prior assumptions, construct the correlation kernel and describe the full modeling algorithm. Section 3 contains a brief validation section, using synthetic data, as well as a case study to showcase the capabilities of our method. We conclude with a discussion in Section 4. The appendix provides further insight into the mathematical footing of the introduced methods.

Magnetic Field Theory
Assuming an insulating mantle, outside of the conducting core, the EMF B can be approximated by the gradient of the geomagnetic potential Φ (Backus et al., 1996): Φ is a scalar potential, satisfying Laplace's equation ∇ 2 Φ = 0. Assuming the sources of the potential lie at some reference sphere with radius R, at locations x outside of this sphere |x| > R the field can be represented using spherical harmonics (SH) (1) x is the unit vector x/|x| and  m Y refers to the real valued and Schmidt semi-normalized SH of degree ℓ and order m with related Gauss coefficient  m g . Similar to MSKH20, we do not consider the Earth's ellipticity.
The dependence of  m g on a reference radius R is not explicitly expressed. The time dependence of the field is typically encoded in the Gauss coefficients  ( ) m g t . We use upright letters x = (x, t) to distinguish spacetime inputs from purely spatial inputs. Often the Gauss coefficients are expressed in form of a spline model (Bloxham & Jackson, 1992) where θ, ϕ, and r are colatitude, longitude and radius of a field location x.
Paleomagnetic records of the EMF are provided as declination D, inclination I and intensity F, which relate to the field vector in a non-linear fashion: The horizontal intensity   is an auxiliary quantity. H is called observation functional.

Gaussian Process Regression
In the eighties C. G. Constable and Parker (1988) already proposed using GPs to model the EMF. A GP is a stochastic process, characterized by a mean function B and a covariance function Given observations o(y) of B at locations and times y = (y, s) with Gaussian measurement errors, characterized by a covariance matrix Σ o , the posterior of B is again a GP. Its (conditional) mean and covariance functions read (Rasmussen & Williams, 2006) x = (x, t) refers to location and time of interest. Note, that herein already a difference to previous GP based models is visible: The covariance function is defined both in space and time. While in principle the truncated spherical harmonics and the B-spline basis may also be used to construct a covariance function, the language and formalism of GP regression have so far only been applied to either temporal correlations (Gillet et al., 2013;Hellio & Gillet, 2018) or spatial correlations S. Sanchez et al., 2016).

Modeling Concept
We propose a fully Bayesian modeling concept, embedded in a functional analytic setting. Therefore, GP based techniques are employed. One key ingredient to GP regression is the a priori covariance function, also called the (correlation) kernel. In this section, we formulate the covariance function we employ, based on our a priori assumptions. Additionally, we formulate the paleomagnetic data model and discuss approximations that are necessary to apply the GP regression scheme.

A Priori Process
Translating the uninformative dipole prior from MSKH20 to a time-dynamic realm presents a challenge, as temporal correlations vanish in the limit of zero a priori precision and thus cannot easily be recovered in the posterior. Instead, the a priori mean function of the EMF is assumed to be constant in time, with only axial dipole contribution: The strength  0 1 of the a priori dipole is a free parameter, that will later be marginalized.
We suggest building the space-time covariance using a sum of tensor products. The proposed closed-form covariance function for the spatial correlations of the non-dipole part K Φ,S,ND includes all SH degrees. Holschneider et al. (2016) describe how to construct this kernel, and in MSKH20 we describe in detail how to adapt it for paleomagnetic applications. In Appendix A we layout this procedure. The strategy is to translate an idea about correlations amongst Gauss coefficients to the potential, using the SH representation. The field covariance function then follows from the gradient. Our main a priori assumption for the covariance is that at some reference sphere, close to where the core field is generated, the geomagnetic Gauss coefficients are uncorrelated. Assuming a flat spectrum at this reference sphere's radius R, it is possible to derive a closed-form for the potential covariance function. This closed form is called Legendre kernel and reads (Holschneider et al., 2016, Equation 53) where b = x ⋅ x′/R 2 and a = |x‖x′|/R 2 .
Temporal correlations are incorporated by a tensor product of this spatial kernel with a squared exponential (SQE) kernel: τ is the a priori correlation time, that gives an idea about the timescale on which the dynamics of the process happen. Note that the posterior curve may be smoother or more detailed than this scale, depending on the data. Similar to MSKH20, we split the kernel into dipole and non-dipole part, as the statistical properties of the dipole are known to differ from the higher field degrees (C. G. Constable & Parker, 1988). Each part is coupled to its own temporal correlation kernel and thus has its own correlation time: The index S stands for spatial and α • are the a priori variances of the dipole and non-dipole part. They can be interpreted as the expected standard deviation in the dipole ( 1 m g ) and non-dipole coefficients. See Appendix A for the explicit forms and further details. The kernel implements a single, constant correlation time τ ND for all degrees ℓ ≥ 2. We are aware, that previous work indicates different behavior (Bouligand et al., 2016). However, implementing the SQE kernel as suggested is straightforward and sufficient for the conceptual work we present here. The field's covariance function reads

Linearization
Paleomagnetic observations are reported as declination, inclination and intensity. With measurement errors E, the data model reads SCHANNER ET AL.
Clearly, the relationship to the field vector is non-linear (Equation 4). Handling non-linear transformations in the framework of GP regression is technically demanding and often analytically impossible, as the transformed random variables are no longer Gaussian. While more sophisticated methods exist (e.g., Snelson et al., 2003), the standard approach is to linearize the observation functional by means of a Taylor approximation of first order. For declination, inclination and intensity, the approximate, linear functionals read Here,    , , D I F and  B indicate the point of expansion (POE). We implement this approximate transformation, to have a linear relation between the observations and the modeled quantity. Linear transformations preserve normality and thus the standard GP formalism is applicable. The proxy data model reads H lin. refers to the linearized observation functionals Equations 14-16.
Previous works implemented an axial dipole as POE (e.g., Hellio & Gillet, 2018). MSKH20 shows that the performance of inversion for the POE can be enhanced if we separate the data into two disjoint groups. One group consists of records with full vector information (complete) and the other of records with at least one component missing (incomplete). In a first step, only complete records are considered and the posterior distribution for these records is calculated. This first step posterior then serves as the prior and POE for the second step, where the remaining, incomplete records are treated.

Measurement Errors
In order to apply the GP regression formalism (Equation 5), the full data model has to be Gaussian. Therefore, linearizing the observation functional as described in the previous section is not sufficient, but a normal proxy error model has to be constructed as well. Intensity records often provide the error as standard deviations of a normal distribution, and thus linearizing the observation functional is sufficient for the intensities. Records of the archeomagnetic directions (declination and inclination) on the other hand, are reported together with the 95% confidence cone (α 95 ) of a von Mises-Fisher distribution. Thus for the archeomagnetic directions, we construct a Gaussian proxy, using (Suttie & Nilsson,  Additionally, similar to MSKH20, we implement a scaling factor ϵ to compensate possible false error estimates, and a residual term P with scaling factor ρ, to address modeling related errors (e.g., observational bias due to crustal field contributions). This way, the data model reads where E prox. are the approximate errors, constructed from Equation 18.

Dating Uncertainties
To a large amount, archeomagnetic specimen are dated using either radiocarbon dating or archeologic age estimation. Both methods suffer from uncertainties, as the former depends on carbon models of the atmosphere and the latter on contextual knowledge. To incorporate these uncertainties, and to represent them in the resulting models, previous studies mostly relied on sampling strategies (e.g., Hellio & Gillet, 2018;Korte et al., 2009;Senftleben, 2019). Hellio et al. (2014) used a normal error model for the dates, and applied Markov Chain Monte-Carlo (MCMC) methods to estimate the posterior distribution. We pursue a similar, hierarchical approach, but instead of MCMC methods, we perform analytic approximations. The dating uncertainties are translated to measurement errors, as presented below, by weighing them with the temporal derivatives of the kernel. This is related, but not equal, to the idea of using the secular variation to estimate the contribution of dating uncertainties (see e.g., Korte et al., 2005). Due to the GP structure of the proposed model, the covariance structure for the secular variation is available a priori. This covariance mediates the dating uncertainties to the measurement errors.
Summarizing the errors from the previous section as ɛ for readability, the data model is However, one does not know the precise time t at which the specimen received its magnetization, but a corrupted date is a normal error. Plugging this into the data model gives This is known as the noisy input Gaussian process (NIGP) (McHutchon & Rasmussen, 2011). Due to the random variable e t appearing at the inputs of the GP B, this data model is non-Gaussian again. To tackle it, McHutchon and Rasmussen (2011) suggest once more a linearization. This gives The first term is normal and the second term allows for easy construction of a moment matching proxy. With this modifications, Equation 23 can be used for GP regression in the usual way. Since the error e t is centered, the a priori mean is not affected by the dating uncertainties. However, the covariance gets an additional term Here Σ tt′ is the dating error covariance matrix and • is the Hadamard product, that is, element wise multiplication along the t direction. To this end, K B (x, x′) is considered as a matrix consisting of 3 × 3 blocks. The effect of the NIGP model is thus the inclusion of dating errors as contributions to the data covariance, similar to measurement errors. The translation is realized by weighing the dating uncertainties by the second order time derivative of the kernel.
10.1029/2020JB021548 6 of 22 In Figure 1 we present a synthetic, one-dimensional example, to compare the proposed NIGP strategy to the standard GP inversion and inference via MCMC. Data was generated from a one-dimensional SQE kernel and assigned large input uncertainties and small errors, to mimic the situation of large dating uncertainties. The standard GP regression shows the typical constrictions at the input points, while the NIGP shows a larger standard deviation, especially at the input points. We believe that an MCMC approach gives a better estimate of the actual posterior, though in a realistic setting this is computationally unfeasible. However, as can be seen from the bottom panel of Figure 1, the NIGP gives a reasonable proxy to the MCMC result at immensely reduced computational cost.

Hyperparameters
The model we constructed throughout this section consists of several parameters, which are a priori unknown. Most obvious are the a priori dipole strength  0 1 , the variances α DP and α ND and the correlation times τ DP and τ ND . Additionally, there are the two scalings, ϵ and ρ, for the measurement errors and the residual, respectively. The least obvious is the kernel's reference radius R. R basically controls the slope of the prior power spectrum. As suggested in our MSKH20 (Figure 2), we take R = 2,800 km, which gives a slope similar to the IGRF power spectra (Thébault et al., 2015). Note, that R gives the radius of a sphere of "virtual" sources and has no direct physical meaning. The other parameters are marginalized, so that the outcome of the modeling procedure is a compound distribution Here η is the set of hyperparameters { , , , , , , }       In the top panel, the standard GP inversion, which can not incorporate input uncertainty, is compared to the noisy input Gaussian process (NIGP) approach. The NIGP approach translates uncertainties in the inputs to uncertainties in the posterior, in stark contrast to the standard GP regression, which shows the typical constrictions at the inputs. In the bottom panel, an MCMC based approach is compared to the NIGP.
One can see, that the NIGP gives a reasonable proxy to the MCMC posterior, which we believe to be a good estimate for the actual posterior. For this small data set, the runtime for the MCMC was 43 s, in distinction to 0.3 ms for the NIGP. This factor of roughly 10 4 makes using MCMC infeasible for the later steps of our proposed modeling procedure. The compound distribution is no longer Gaussian and includes modeling uncertainties, resulting from the a priori lack of knowledge about the hyperparameters, but does not depend on these parameters. This distribution is the central result of the suggested modeling strategy. It includes a most probable field model, the mean of the distribution, as well as uncertainty estimates, resulting both from modeling and measurement process. Similar distributions can be obtained for other quantities of interest, such as the Gauss coefficients.
To actually evaluate these expressions, numerical approximations have to be employed, as described in the following section. The integral is approximated by a sum, which results in a Gaussian mixture distribution. Moments for this mixture are easily obtained.

Application
In this section we demonstrate the potential of the suggested modeling scheme, by first applying it to synthetic test data and finally conducting a case study based on actual archeomagnetic records. The major task to this end is the implementation of the covariance matrices from the kernel, the linearization and the twostep strategy. As this process is described in detail in MSKH20, we outsource it to Appendix B. However, two points are to be discussed here. One is the explicit second derivative of the temporal kernel, appearing in Equation 24. In this study the correlation kernels for dipole and non-dipole contributions are each considered tensor product kernels, so that the time derivative only affects the temporal part. The full derivative of each contribution is then given by the product of the respective spatial part and the second order derivative of the SQE-kernel. For the SQE-kernel, the derivative is straightforward to calculate and reads When using a more realistic kernel, especially one with different correlation times for different degrees ℓ, the kernel can no longer be constructed as a tensor product, and calculating this derivatives analytically may pose a challenge, so that numerical methods have to be employed. This is one reason why in this conceptual study we chose the SQE kernel over a more realistic one e.g., the one proposed by Gillet et al. (2013).
The second point is the marginalization integral in Equation 25. As the proposed model contains seven parameters and as the data is incorporated all at once, instead of in bins as in MSKH20, the brute-force parameter space exploration and integration suggested in MSKH20 are now computationally unfeasible. This is due to the cost for fixed point integration growing exponentially with the number of dimensions and the high cost of matrix inversion ( ).
Instead, we perform numerical integration similar to the strategy suggested by Rue et al. (2009, Section 6.5).
The idea is to center the integration around the maximum a posteriori probability estimator (MAP) of the marginal posterior. Collocation points are added according to central composite design (CCD) (S. M. Sanchez & Sanchez, 2005), in order to capture the bulk of the uncertainties in the hyperparameters (see Figure 2). In seven dimensions, the integration is approximated by a sum over 79 collocation points. To find the MAP of the marginal posterior, we use the LIPO-TR global optimization algorithm (King, 2009(King, , 2017. The parameters are assigned box priors, as is required by most global optimization algorithms. We choose as upper and lower bonds for the hyperparameters  where • stands for DP and ND. The poles induced by Jeffery's priors do not cause trouble, as the box constraints are far enough from zero. By numerically approximating the integral in Equation 25, the compound distribution is approximated by a Gaussian mixture. Details can again be found in Appendix B.

Synthetic Tests
To validate the proposed algorithm, we test it on synthetic data. As inputs we choose dates and locations from the archeomagnetic and volcanic data in GEOMAGIA v3.3 (Brown et al., 2015) for the interval from 800 to 2,000 (cf., Figure 5). At these locations we generate data from the ARCH10k.1 model (C. Constable et al., 2016) and corrupt it by artificial noise. For the directions we use a von Mises-Fisher distribution.
Intensities are corrupted by gamma distributed noise and the dates by normal noise. The error levels are taken from GEOMAGIA as well. Additionally, to simulate effects of other magnetic field sources and from the measurement process, we add random contributions with a constant standard deviation of 2.5 μT to the generated data.
The resulting MAP for the hyperparameters is . Scaling the non-dipole variance α ND is not possible, due to the degree dependent scaling and the closed-form kernel. The error level scalings are as expected. As the reported error levels have been used to corrupt the data, they are accurate (ϵ ≈ 100%). The residual term reflects the random contributions that have been added.
In Figure 3 we compare the dipole of the reference model to our findings from running the model on the synthetic data. One can see that only the long-term behavior of the reference model can be recovered. Nonetheless, we believe that the results are meaningful and related to the data distribution. To illustrate this, we present a simplified example. Consider a one-dimensional time series that we want to model by a Gaussian process. We generate data from some reference model (here we chose the axial dipole of the ARCH10k.1 model) and corrupt it with artificial noise. For simplicity and to emphasize the effect, we neglect dating uncertainties. We synthesize two data sets, one with many observations and small uncertainties and one with SCHANNER ET AL.
10.1029/2020JB021548 9 of 22 Comparison of the reference model to the one recovered from synthetic test data. The shaded area represents one standard deviation. The long-term behavior could be recovered, while data is too sparsely distributed and uncertainties are too large to recover the short timescale features. See the text for additional discussion. few observations and large errors. The "data poor" situation was designed to mimic the actual situation, based on findings in MSKH20: Every one hundred years an observation was generated and the error level chosen to agree with the one found in MSKH20.
Then for both data sets, a separate correlation time is estimated by maximizing the marginal posterior. In the situation with many observations, the proposed modeling algorithm gives a correlation time corresponding to the short term variability of the reference model and the reference model is recovered well and detailed (Figure 4, left panel). In the situation with few observations; however, one can only recover the long-term variability of the reference model (Figure 4, right panel). The actual situation is more complicated, but we believe that the large correlation times and smooth curves come from a similar mechanism. While locally, the data distribution is much denser than in the "data poor" case, global information (e.g., about the axial dipole) is sparse and unreliable, as can be seen from the results of MSKH20. We performed a similar experiment for synthetic field data and found similar results, that is, only the long-term information can be recovered.

Case Study
Here we present the results from applying the proposed modeling strategy to actual archeomagnetic data, taken from GEOMAGIA v3.3 (Brown et al., 2015). The data covers the interval (800, 2,000) AD and consists of 7,801 records, of which 3.9% are complete vector triples. Figure 5 shows the spatial and temporal distribution of the data. Similar to MSKH20, we do not consider the Earth's ellipticity. We use the originally reported error estimates and assign α 95 = 4.5° as directional errors and σ F = 8.25 μT as intensity errors to the 8.4% of data, where no error is reported. As described above, the dating uncertainties are considered as standard deviations of independent normal distributions. When different values for upper and lower temporal error are reported, we use the bigger value. The 0.7% of data for which no dating uncertainty is reported are assigned a standard deviation σ t = 100 years. Note that for the recent times fewer records are available from the archeomagnetic data set. This results in bigger uncertainties toward recent times, as can be seen for example in the inclination series in Figure 8.
We compare our findings to two existing magnetic field models as well as to the results of MSKH20. The models are ARCH10k.1 (C. Constable et al., 2016) and COV-ARCH (Hellio & Gillet, 2018). They are considered reasonably comparable, as they are based on similar data compilations. Both models report Gauss coefficients up to SH degree ℓ = 10, while the actual spatial resolution is determined by regularization in ARCH10k.1 and by cross-covariances based on prior assumptions in COV-ARCH and lies around ℓ = 5 in both cases. ARCH10k.1 does not report uncertainties, while COV-ARCH provides an ensemble of 50 models from which uncertainties are constructed by calculating sample mean and sample standard deviation. In SCHANNER ET AL.  . Synthetic one-dimensional example illustrating how the lack of data influences the information one can recover from the reference model. The shaded area represents one standard deviation. In the data rich situation on the left, the reference model can be recovered quite well and posterior uncertainties are so small, they are barely visible. In the data poor situation on the right only the long-term behavior can be recovered.
contrast to the results presented by Hellio and Gillet (2018), the publicly available model (https://geodyn. univ-grenoble-alpes.fr/) is not time continuous but reports coarsely binned coefficients for every hundred years in the interval.
We want to stress again, that the presented results stem from a conceptual design. Especially the common temporal correlation time for all degrees ℓ ≥ 2 should be reconsidered, when building an actual model from the proposed strategy. The hyperparameter-MAP for the actual data is quite similar to the one in the synthetic data test (see Section 3.1): T. All values are in a reasonable order of magnitude. Surprisingly the residual scaling did not decrease in comparison to MSKH20. The proximity of both correlation times may be explained by the dominance of the quadrupole over the larger degrees. For further insight, we provide profiled distributions together with the modeling software .

Field Predictions
Predicting on the EMF's vector components is straightforward and given by wrapping Equations 6 and 7, into the marginalization procedure described in the appendix. However, to get a reasonable spatial resolution too many design points have to be included to store the covariance matrices for all integration points, which is necessary to calculate the mixture distribution. Instead, similar to MSKH20, we resort to calculating the moment matching Gaussian proxy. The top row of Figure 6 shows mean and standard deviation for the down component of the EMF for the epoch 1700. Similarly, predictions at the core-mantle boundary (CMB) can be performed, by translating the design points accordingly. The results for the down component are shown for the epoch 1700 in Figure 7. Inferring the observables (declination, inclination, and intensity) is hindered by the non-linear relation to the field. Utilizing again a linearization, mean and covariance can be constructed, similar to MSKH20 and Hellio et al. (2014, Appendix A). The bottom row of Figure 6 shows a prediction of the EMF's intensity for the epoch 1700.
Compared to MSKH20 for snapshots in time, the new results show a slightly lower strength of the down component and lower field intensity. Moreover, the new results have a reduced standard deviation, which can be attributed to additional constraints due to the temporal correlations. The overall field structure is similar, showing features such as the South Atlantic Anomaly. The reconstruction at the CMB (Figure 7) reveals a region of lower field strength at the southern tip of Africa, which is at the limit of significance, but was not present in the snapshot model.  . Spatial and temporal distribution of the data. We separate the data only into complete and incomplete without indicating declination, inclination and intensity separately, to avoid overloading the plots. The turquoise dots indicate the locations at which we present the time series (Paris, Figure 8 and Pacific, Figure 9). ters to infer is again too big when considering a reasonable temporal resolution. Figures 8 and 9 show time series at two distinct locations. Figure 8 presents time series for Paris, together with comparison models and data, while the series in Figure 9 are for a location in the Pacific, where no data is present in the surroundings. The data in Figure 8 stem from a surrounding of 250 km.
Inclinations and intensities are translated to the coordinates of Paris along the corresponding axial dipoles (Merrill et al., 1996). Declinations are taken as reported. For Paris, the different models mostly agree, with larger deviations toward the recent epoch, where the database is thinning out. Our inferred series shows less variation than comparable models, which probably can be attributed to our somewhat unrealistic temporal kernel with equal correlation time for all spatial wavelengths from ℓ = 2 on. As expected, the uncertainties are bigger for the location in the Pacific, and the comparing model series show larger deviations.

Gauss Coefficients
Although the proposed model is inherently non-parametric in both space and time, predictions on Gauss coefficients can be performed. As they are linearly related to the field, the procedure is straightforward (see Equation 21 in Holschneider et al., 2016). When predicting coefficients for a specific epoch, the full mixture distribution is accessible. When predicting on coefficient time series; however, the number of parameters one has to include in the prediction to get a reasonable temporal resolution is too memory intense to store SCHANNER ET AL.     Differences to ARCH10k.1 and COV-ARCH might partly be due to some differences in the underlying data compilation.

Spectra
Power spectra are considered to condense the information contained in the Gauss coefficients (e.g., Backus et al., 1996). Using sampling techniques, mean and percentiles of the geomagnetic power spectrum distribution are available (for further details consider MSKH20, Section 5.6). Figure 11 shows the geomagnetic SCHANNER ET AL.  power spectrum for the epoch 1700, together with 16-and 84-percentiles. Within the reported uncertainties, the power in dipole and quadrupole agrees to the comparison models. For ℓ = 3 the reconstruction reports less power than COV-ARCH and MSKH20. Noteworthy is the faster power decay for degrees ℓ = 4 … 7 when compared to MSKH20, which also indicates larger deviations from the prior. This may be due to temporal correlations increasing the information or the long correlation time damping small scale structures. Implementing a separate correlation time for each coefficient may provide further insight.
Similar to the geomagnetic power spectrum, the spectrum of the secular variation can be calculated (Alldredge, 1984). Therefore, one has to predict on the derivatives of Gauss coefficients, also called the secular variations. As the derivative is a linear operator, this is straightforward. To explain the basic concept, consider the simplified example of direct observations of the EMF: , Since the a priori mean in our model is constant, the first term vanishes. Thus, to predict on the secular variations one only has to calculate the derivative of the correlation between coefficients and observations, which in the suggested model reduces to the first derivative of the SQE kernel: All other quantities are known from inferring the Gauss coefficients themselves. The covariance of the secular variation translates analogously. Mean and percentiles of the secular variation spectrum are then again available via sampling.
As MSKH20 consists of snapshot models, the respective secular variation spectrum is not accessible. Similarly, as the publicly available version of COV-ARCH only reports values every 100 years, no secular variation can be calculated and the spectrum is missing in Figure 12. Note, that we choose an earlier epoch for the secular variation spectrum, as ARCH10k.1 is constrained to gufm1 (Jackson et al., 2000) for the recent times and therefore shows higher than average secular variation for these centuries. The secular variation spectra for the two models are fairly similar, with a very good agreement for the dipole and slightly higher values for the higher degrees in our new model.

Dipole
Finally we present the dynamics of the EMFs dipole. Figure 13 shows time series of the dipole moment magnitude. The magnitude is higher than the ones reported by comparison models, while the "outliers" SCHANNER ET AL.
10.1029/2020JB021548 15 of 22 from MSKH20 are not present in the new results. From 1840 on ARCH10k.1 is constrained by the gufm1 model, which in turn is constrained by a large amount of direct observations and can be considered to represent the dipole moment quite reliably from that time on. The deviation of our model from ARCH10k.1 during the last century is likely caused by a lack of archeomagnetic data for these epochs. Figure 14 shows the movement of the geomagnetic north pole. The mean curve of our model (black line) is calculated via sampling. For a given epoch the full density is available analytically (cf., Mauerberger et al., 2020, Equation 96). To not overload the plot, we only show mean and one-sigma ellipses every century. The stereographic projection is responsible for the crescent-shaped distortion of the ellipses. The rapid movement of the dipole for earlier epochs suggested by other models is not found by our new reconstruction and for more recent times, the path lies further to the west. Deviance for the most recent epochs is again caused by a lack of data.

Conclusions
The presented work extends the Bayesian strategy for correlation-based modeling of the archeomagnetic field introduced in MSKH20 to the temporal domain. In Section 2, all necessary modifications are discussed, together with a novel approach to include dating uncertainties. In contrast to previous works (Hellio & Gillet, 2018;Hellio et al., 2014;Korte et al., 2009;Licht et al., 2013;Senftleben, 2019), using a NIGP (McHutchon & Rasmussen, 2011) to incorporate dating uncertainties does not rely on sampling techniques. The a priori model is again constructed with the aim of being as objective as possible. The uninformative dipole prior from MSKH20 cannot easily be translated to the time-dynamic realm, as temporal correlations vanish in the limit of zero a priori precision and thus cannot easily be recovered in the posterior. Instead, we assume a priori a constant axial dipole with the dipole strength being a free parameter. Together with all but one other model parameters, the dipole strength is marginalized so that the model does not depend on the specific value. This marginalization presents another challenge, as numerical integration in a seven dimensional space has to be performed. The Riemann sum approach from MSKH20 is not applicable, due to the curse of dimensionality (i.e., unfeasible computation time). As a practicable alternative to the brute force integration we implement a CCD SCHANNER ET AL.  (S. M. Sanchez & Sanchez, 2005) based integration, as proposed by Rue et al. (2009). The major challenge in implementing this strategy consists of finding the MAP of the marginal posterior. Running the LIPO-TR algorithm (King, 2009(King, , 2017 for the validation and case study data sets took around 25 h each on a regular workstation. Once the MAP is found, the set of integration points consists of 79 hyperparameter combinations and the marginalization takes between half an hour and 5 h, depending on the quantity one predicts on, as some quantities require sampling, which is more computationally demanding. With the marginalization performed, the model depends only on the a priori choice of the Gauss coefficients covariance structure at the reference radius, the value of the reference radius and the temporal covariance structure. For the conceptual work presented, we chose an unphysical SQE kernel. This has the advantage of being easy to implement, but the main point why we use this kernel instead of a more reasonable one is the necessity to calculate temporal derivatives to implement the NIGP. For the SQE kernel this is straightforward. We have shown by means of a synthetic test and a case study on real data from 800 AD to 2000 AD that even with the simplified kernel the results compare well with previous models. Notably the "outliers" for the years 1100 and 1300, present in MSKH20, do not appear in the present work. This may be explained by the new model considering temporal errors, and thus covering a false binning, or by the smoothness imposed by the SQE-kernel, which suppresses the influence of single records. The results of the case study show smoother curves than the models we compare to, in particular for the axial dipole contribution. While partly this can be attributed to the SQE-kernel, the proposed method is able to resolve faster variations locally, if the data supports these variations. Due to the closed form, such variations can be reflected by higher order coefficients and do not introduce global variations. Implementing a more realistic kernel, such as the one proposed by Gillet et al. (2013), will be the direction of future work. Together with a Bayesian framework for data selection this will allow the construction and proposition of a new correlation based field model. In MSKH20, the expansion of the database by records from ship logs was discussed. Incorporating uncertainties arising from imprecise locations may be performed by the use of the proposed NIGP. Instead of temporal derivatives, the spatial gradients of the kernel are used to translate the input uncertainties to contributions to the measurement errors. The challenges to scale relative intensities and preserve stratification (Nilsson et al., 2014) persist, so that sediment records require a different approach than the application of the NIGP proposed here.
We again developed a python framework to save the effort of implementing the proposed algorithm . Together with extensive documentation, the software source code provides further insight into the modeling algorithm. It is available at https://sec23.git-pages.gfz-potsdam.de/korte/ corbam/.

Appendix A: Constructing the Spatial Covariance Kernel
The Lagrange kernel (Equation 9) is constructed from the Gauss coefficient correlations as follows: Consider the covariance of the magnetic potential in SH decomposition. Then, for a potential of internal origin Assuming that at some reference sphere the Gauss coefficients are uncorrelated with a flat spectrum, that is, Following Holschneider et al. (2016), evaluating the sums gives the kernel Equation 9 where b = x ⋅ x′/R 2 and a = |x‖x′|/R 2 . The dipole kernel can be extracted by setting ℓ = 1. This yields Thus the non-dipole kernel reads where the last term excludes the monopole (ℓ = 0).

Appendix B: Detailed Modeling Algorithm
The modeling algorithm consists of two stages. The first one deals with finding the MAP of the hyperparameters. The MAP then serves as a center point for marginalizing the hyperparameters in the second step. We begin this section by laying out the inversion process. From quantities that are calculated along the way, the marginal posterior can be constructed. Using both procedures, the full algorithm can presented in a compact way.

B1. Inversion
The inversion closely follows the modeling concept described in MSKH20. To provide insight into the mathematical background of the inversion, we lay out the full inversion process for the field B at locations of interest y. Inverting for other quantities, such as the Gauss coefficients or the field's intensity, is straightforward.
To keep the equations concise, we use the following notation for matrices: ( ) x B is the a priori (mean) field at locations x   ( ) | x H B refers to the gradient of the observation functionals (cf., Equations 14-16), evaluated at the POE  B at locations x ( ) x H refers to the linearized, transformed prior field at locations x The transformed prior field serves as a mean proxy to observations at locations x.
Σ yy refers to the a priori covariance of the field at locations y. This is This is a matrix, composed of 3 × 3 blocks, containing correlations at each point.
is short hand for the posterior mean of the field at locations y, given observations o, that is, Σ oo refers to the covariance amongst observations: The residual Σ p is an identity matrix of the dimension of number of observations and Σ e is the (typically diagonal) matrix of approximate measurement errors, see Section 2.3. Σ T is the correction term for dating uncertainties, see Section 2.4: As the linearization is tackled by means of a two-step strategy, at first the data is partitioned into complete c and incomplete i records: The first step in the two-step strategy only deals with complete records. However, as the posterior mean from the first step serves as the point of expansion (POE) in the second step, predictions at locations of incomplete records x i have to be included as well. The posterior mean and covariance from the first step read In the second step, the posterior mean of the first step is used as POE: SCHANNER ET AL.

10.1029/2020JB021548
19 of 22 Exploration The first stage consists of finding the maximum a posteriori probability estimator (MAP) ˆ of the marginal posterior. Therefore, the log marginal posterior is optimized using global optimization techniques Integration With the MAP as center, a set of integration points   is constructed as described by Rue et al. (2009) and S. M. Sanchez and Sanchez (2005). With weights Δ η , the integral for the compound distribution Equation 25 is approximated by a sum This way the compound distribution is approximated by a Gaussian mixture. Similar expressions exist for the compound distributions of all quantities of interest, such as Gauss coefficients or observables like the linearized intensity F.

Data Availability Statement
There is no new data involved in this publication. The data used in the case study is available via GEOMA-GIA v3.3 (Brown et al., 2015). All results were produced using a python implementation of the discussed algorithm, which is publicly available .