ArchKalmag14k: A Kalman‐Filter Based Global Geomagnetic Model for the Holocene

We propose a global geomagnetic field model for the last 14 thousand years, based on thermoremanent records. We call the model ArchKalmag14k. ArchKalmag14k is constructed by modifying recently proposed algorithms, based on space‐time correlations. Due to the amount of data and complexity of the model, the full Bayesian posterior is numerically intractable. To tackle this, we sequentialize the inversion by implementing a Kalman‐filter with a fixed time step. Every step consists of a prediction, based on a degree dependent temporal covariance, and a correction via Gaussian process regression. Dating errors are treated via a noisy input formulation. Cross correlations are reintroduced by a smoothing algorithm and model parameters are inferred from the data. Due to the specific statistical nature of the proposed algorithms, the model comes with space and time‐dependent uncertainty estimates. The new model ArchKalmag14k shows less variation in the large‐scale degrees than comparable models. Local predictions represent the underlying data and agree with comparable models, if the location is sampled well. Uncertainties are bigger for earlier times and in regions of sparse data coverage. We also use ArchKalmag14k to analyze the appearance and evolution of the South Atlantic anomaly together with reverse flux patches at the core‐mantle boundary, considering the model uncertainties. While we find good agreement with earlier models for recent times, our model suggests a different evolution of intensity minima prior to 1650 CE. In general, our results suggest that prior to 6000 BCE the data is not sufficient to support global models.

2 of 17 all (Constable et al., 2000;Jackson et al., 2000;Korte & Constable, 2003) or relied on 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). In contrast, Nilsson and Suttie (2021) (and earlier Hellio et al. (2014) for local field models) used a Bayesian formulation of the proposed Gaussian process (GP) approach, to estimate uncertainties based on the posterior distribution. Holschneider et al. (2016) extended the GP approach to the spatial domain, to also reflect uncertainties resulting from the data distribution, and in two recent studies this method was adapted to paleomagnetic records Schanner et al., 2021). The major challenge with the modeling strategies proposed there is related to the inversion of large-scale matrices, and the methods were found computationally unfeasible for the number of records available for the Holocene. In the area of modeling the recent field, this challenge was overcome by applying sequentialization by means of a Kalman-filter (Kalman, 1960) to the inversion problem Ropp et al., 2020). This way, models from a way higher number of satellite observations have been constructed, while retaining the strategies proposed by Holschneider et al. (2016). In this study we apply sequentialization to the earlier developed strategy (Schanner et al., 2021, in the following referred to as SMKH21) and propose a new global geomagnetic field model for the Holocene.
Global geomagnetic field models on archeological scales are inferred from two classes of data: Data from materials with thermoremanent magnetization, such as volcanic rocks, bricks or burnt clay fragments from archeologic sites, and data from marine or lacustrine sediments with embedded magnetic particles. As in earlier studies, we focus on the former class and loosely refer to it as archeomagnetic data in this paper. The extension to sediments poses several additional challenges, some of which are addressed and discussed by Nilsson and Suttie (2021). The a priori model that results from the sequentialization of SMKH21 is similar to the one proposed by Nilsson and Suttie (2021). Besides a focus on a different and smaller data set, the main difference lies in the inversion procedure: While Nilsson and Suttie (2021) employ a Markov Chain Monte-Carlo (MCMC) based strategy, we rely on Kalman-filter based inversion. This method relies mostly on linear algebra and may therefore be called deterministic, in contrast to the sample based and therefore probabilistic MCMC approach.
The rest of this article is structured as follows: In Section 2, we discuss prior assumptions, showcase the modeling method and introduce the data set. Section 3 contains a brief validation section, using synthetic data, but mainly focuses on the description of features of the new model, which are discussed in Section 4. We conclude in Section 5 by reconsidering possible extensions and shortcomings of the method, as well as an outlook to future work.

Gaussian Process-Based Modeling
In the 1980s, Constable and Parker (1988) proposed using GPs to model the Earth's magnetic field (EMF). The technique was later applied by Gillet et al. (2013) and extended by Holschneider et al. (2016). A GP is a stochastic process that is uniquely characterized by a mean function ̄ and a covariance function K B ∼  (̄, ) . (1) Gaussian process-based modeling is a Bayesian approach, where a GP is used as a prior and an update is given by some normal likelihood. The posterior is then a GP as well, so that the model is also uniquely characterized by a mean function and a covariance function (Rasmussen & Williams, 2006). The main difficulty in applying this technique to paleomagnetic records lies in constructing the normal likelihood, as archeomagnetic observations are nonlinearly related to the magnetic field.

Data Model
To apply GP based modeling, one has to construct a normal likelihood, relating observations to model predictions of the magnetic field. In paleomagnetism, the observations are the field directions (declination D and inclination I) and intensity F. At locations x and times t, the data model can then be formulated as: ≈̃⊤ . (5)̃̃ , and ̃ indicate the point of expansion (POE) and we summarize the linearized expressions as H lin. (see also Mauerberger et al., 2020). The observation errors E are also non-Gaussian, as the directional errors are given by a Fisher-von Mises distribution. We approximate this two dimensional distribution with 95% confidence cone (α 95 ) by two centered normal distributions with standard deviations (Piper, 1989; where o I is the observed inclination. We label these approximate errors E prox. . Next, we consider dating uncertainties as suggested in SMKH21. The precise times t at which the archeomagnetic specimen received their magnetization are unknown. Instead, a corrupted date t o = t + e t is reported, and we consider e t to be a centered normal error, even though dating errors for archeological artifacts may have a non-Gaussian distribution. This error in the inputs is handled by another linearization, as proposed by McHutchon and Rasmussen (2011, the noisy input Gaussian process (NIGP)). As the errors are centered, the a priori mean is not affected by this procedure. However, via linearization the dating uncertainties are translated into observation uncertainties, and the covariance gets an additional term Σ tt′ is the covariance matrix of the dating errors and • is the Hadamard product (see also Schanner et al., 2021, Section 2.4). The effect of the NIGP model is thus the inclusion of dating errors as contributions to the data covariance, similar to measurement errors. To realize this, dating errors are weighed by the second order time derivative of the kernel. The idea is related, but not equal, to the approach of estimating the contribution of dating uncertainties by using the secular variation (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. Finally, a residual term ξP, where ∼  (0, ) is added to cover any effects that are not modeled, like crustal field or ellipticity of the Earth (see also Mauerberger et al., 2020;Schanner et al., 2021). This way, the data model reads

A Priori Process
We consider the common spherical harmonics (SH) expansion of the geomagnetic potential Φ, which is valid outside of the Earth's conducting core, assuming an insulating mantle: is the unit vector x/|x| and refers to the real valued and Schmidt seminormalized SH of degree ℓ and order m with related Gauss coefficient . Here, the Gauss coefficients are defined at the reference radius R. Note that in the figures below, we give the extrapolated coefficients at the Earth's surface. From this, the Earth's magnetic field is given as the gradient and mean and covariance function of the EMF can be derived from assumptions about correlations of the Gauss coefficients. A priori we assume all Gauss coefficients except for the axial dipole to be of zero mean. The mean function for the axial dipole is assumed constant, with value 0 1 . We assume all coefficients to be uncorrelated and identically distributed at a reference radius R = 2,800 km (within the Earth's core). This is the "virtual" source region where the spectrum is flat and the field has no direct physical meaning. The magnetic field given by this assumption is only a valid representation of the actual field above the core-mantle boundary (CMB). Inside of the core it can be seen as an artificial connection of the physical field at the CMB to the virtual sources inside of the core. We assume two different a priori variances, one for the dipole coefficients α DP and one for all higher degrees α ND . For each coefficient we assume a temporal correlation in the form of an AR(2)-process, as proposed by Gillet et al. (2013) and employed also by others Hellio & Gillet, 2018;Nilsson & Suttie, 2021;Ropp et al., 2020). This way, the temporal correlation of each coefficient is given by Similar to Baerenzung et al. (2020), we assume one correlation time τ DP for the dipole and a relation for all higher degrees ℓ ≥ 2 The posterior may be smoother or more detailed than these timescales, depending on the data.

Sequentialization
In previous studies Schanner et al., 2021), we aimed at performing standard GP regression in the introduced setting. However, as determining the hyperparameters of the model requires this regression to be performed many times, this proved to be computationally unfeasible. To overcome this, we perform a sequentialized inversion, in form of a Kalman-filter Kalman, 1960). Starting at an initial time, the Kalman-filter consists of a series of steps, each consisting of a prediction based on the current model and a correction, which updates the model if data is available. In contrast to the previous study SMKH21, this requires us to define a cutoff degree ℓ max , so that the model can be characterized by a finite vector of coefficients and their derivatives = ( ,̇) . The prediction equations from step i to i + 1 are given by. where is the forward operator of the AR(2)-process and ̃= − ⊤ with the a priori correlations Σ. The correction step consists of a Bayesian GP inversion, as described in detail in SMKH21. The linearization is performed around the current model, beginning with the prior. We run the Kalman-filter "backwards," that is, from modern times to the past, as the data distribution is sparser toward earlier years. We expect the bigger amount of data in the beginning of the filtering to constrain the model and improve the POE for earlier times. We choose a cutoff degree of ℓ max = 20 and a step size of Δt = 10 years. Both choices are believed to allow for a way higher resolution 10.1029/2021JB023166 5 of 17 than present in the data, so that every dynamic present in the data can be captured by the model. After running the Kalman-filter we run a smoothing algorithm, following the formulation of Rauch et al. (1965) [see also Baerenzung et al. (2020)]. This way, cross correlations that are not present in the Kalman-filter are reintroduced to the posterior.
We store a set of coefficients every 50 years, so that the output of a sequentialized inversion consists of 281 sets of 440 main field coefficients, 440 secular variation coefficients and the respective covariances.

Hyperparameters
The a priori model depends on several parameters, that have to be inferred before the actual inversion can be performed. One approach (e.g., Hellio & Gillet, 2018;Nilsson & Suttie, 2021) is to infer these parameters from outside knowledge, for example, from models based on observatory and satellite data. We followed this approach in selecting the reference radius R, which effectively controls the slope of the a priori spectrum, by comparison to the IGRF models. For the other parameters we suggest a more self-consistent strategy and estimate them based on a maximum likelihood procedure. This strategy did not work for the reference radius, most likely because the sparse data in earlier years do not constrain it well enough.
Consider the forward log-marginal likelihood with observations in the ith step o i and their observation covariance Σ o,i . N refers to the number of steps in the Kalman-filter. The forward likelihood depends on the hyperparameters and is considered a measure for how good a choice of hyperparameters describes the data. We maximize this expression using LIPO-TR (King, 2009(King, , 2017 and use the maximum estimator for the parameters in the inference. The search region is specified by lower and upper bounds for the hyperparameters, these are as follows: where • stands for DP and ND. Note that 0 1 and the α • are given at the reference radius R = 2,800 km.

Data Set
The data set is a slight variation of all records from the archeological and volcanic database from GEOMAGIA v3.4 (Brown et al., 2015) with ages between 12000 BCE and 2000 CE. Part of the records from Mexico contain wrong age and dating uncertainty estimates (Mahgoub, personal communication). Some of them are too old by several thousand years and have been removed, while for others updated 14 C ages have been published for the lava flows that they originate from. These updates have been used (see Table S1 in Supporting Information S1) and will be included in GEOMAGIA in the future. To identify other records that deviate from the rest, we use a Naive Bayes classifier (e.g., Berrar, 2018). This procedure is integrated into the Kalman-filter as follows: When a step i + 1 contains new data, we evaluate the probability of every record to either come from a normal distribution with standard deviation of the size of the reported error or from a flat distribution of larger variance ((100°) 2 for declination, (50°) 2 for inclination and (100 µT) 2 for intensities). Records that are more likely to stem from the flat distribution are considered outliers. In comparison to the standard approach of rejecting all data that deviates by a specific amount from the model, this procedure is more flexible and allows larger deviations, especially if the current model reports high uncertainties. By this procedure 276 records are identified and removed from the data set. The final data set contains 18,735 records from 11,637 locations. It consists of 5,611 declinations, 7,028 inclinations, and 6,096 intensities.
In the geomagnetic community, it is common to use L1 or Huber norms and employ reweighing techniques to address outliers (Olsen, 2002;. Hellio et al. (2014) and Nilsson and Suttie (2021) implement longer tailed error distributions and include all records in their models. In contrast to these approaches, we resort to outlier rejection for two practical reasons: first, the precise influence of the error distribution on the model is a question that is yet to be addressed by the community. Records which are accepted by the naive Bayes classifier are likely to be described well by a Gaussian error model, and by rejecting the rest we postpone a detailed analysis. Second, by considering a Gaussian error model only, the inversion is feasible analytically and can be addressed by numerical linear algebra, without resorting to sampling techniques. We acknowledge that the error model is a point that could be addressed more thoroughly in the future (Figure 1).

Validation
In order to validate the proposed modeling method, we performed a test inversion on synthetic data. We therefore set up a model with fixed hyperparameters and sampled coefficients from the prior distribution, which serve as reference. From these coefficients we generated data at the same input locations and times as the ones in the data set described in Section 2.6. The data was then corrupted by artificial noise from a Gamma distribution for the intensity and a von Mises-Fisher distribution for the directions and by normal noise in the ages. The error levels reported in the database were used. Table 1 shows the fixed hyperparameters and the inferred ones. Apart from one parameter they agree reasonably well. The deviance in the nondipole correlation time is likely due to the data distribution. We believe that the variations that are present in the data can be resolved with the larger a priori correlation time and shorter variations cannot be recovered. No additional contributions (white noise) were added to the synthetic data set and the algorithm chooses the lowest possible value for the residual scaling accordingly. Figure 2 shows generated and inferred axial dipole and quadrupole. The quadrupole behavior appears to be better recovered than the dipole behavior in the earlier times (12000 BCE to 2000 BCE). This may be due to the data distribution, as indicated by Figure 3. We there show the relative covariance gain (the diagonal of the last term on the right hand side of Equation 7 in Schanner et al. (2021), relative to the prior covariance) per coefficient for a stationary inversion. This quantity only depends on the data distribution. With the actual distribution, the gain is stronger in the quadrupole. With a symmetrized data set, large-scale degrees are resolved more evenly. A similar behavior is also visible in the power spectra of the test inversion, which are provided with Figure S2 in Model Note. 0 1 is the constant a priori axial dipole, α DP and α ND give the a priori scaling of the dipole and nondipole covariance kernel, respectively. τ DP and τ ND give the corresponding a priori correlation times. ξ is the scaling factor of the residual term. 0 1 and α • are given at the reference radius. Supporting Information S1. In general, Figure 2 shows a promising agreement, although some variation in the dipole, prominently between 10000 and 8000 BCE, is not present in the inferred model. This already hints at the data not containing enough information to recover global features during early times. Further figures from the validation process, showing the other dipole and some higher order coefficients, are available with Figure S1 in Supporting Information S1.

ArchKalmag14k
In the following, we propose and describe a new global geomagnetic field model, based on archeomagnetic records. It covers the last 14,000 years and we call it ArchKalmag14k, as it is based on methods similar to the Kalmag model by Baerenzung et al. (2020). The hyperparameters that maximize the marginal likelihood and define the prior used for constructing the model are given in Table 2. Table S2 in Supporting Information S1 compares the a priori correlation times calculated from τ DP and τ ND to the ones used by Hellio and Gillet (2018) and Nilsson and Suttie (2021). The general order of magnitude agrees. In general, the correlation times we find are slightly lower than what is calculated for the other studies. The largest deviance is observed in the octupole correlation time, where the other studies give a value that is twice the one we find. Below we compare ArchKalmag14k to the models ARCH10k.1 (Constable et al., 2016) and SHA.DIF.14k (Pavón-Carrasco et al., 2014), as both rest on a similar database and cover a similar timespan.  Running the inversion as described in Section 2 gives 281 sets of 440 main field and 440 secular variation coefficients together with the respective covariances, one set every 50 years. A comparison of the model coefficients to the prior is given with Figure S3 in Supporting Information S1. Figure 4 shows the dipole and axial quadrupole and octopole coefficients together with 95%-uncertainties and comparison models. The proposed model ArchKalmag14k shows less variation in the dipole degrees than comparable models, especially during earlier times when data is sparse. More variation is present in the quadrupole and octopole, with variation decreasing toward earlier times.
This behavior is also reflected in the power spectra. Figure 5 shows the spatial (top row) and secular variation (bottom row) spectra for two selected epochs, one with dense (1000 CE) and one with sparse (6000 BCE) data coverage. The blue lines show the power spectrum as a random variable (i.e., a quantity nonlinearly derived from the posterior GP, see also Mauerberger et al., 2020, Section 5.6), together with the corresponding prior as a light blue dashed line. These curves represent the nonlinear transformations of the prior and posterior distribution. We also plot the power spectrum of the mean model (gray lines), that is, the power spectrum directly inferred from the mean coefficients. The random variable gives higher values than the mean and comparison models, as it also includes the variance of the coefficients. The random variable can be compared to the prior, to determine the model resolution, while the power spectrum of the mean is better suited for comparison to existing models. For the recent epoch, the spectrum lies between the one for ARCH10k.1 (orange) and SHA.DIF.14k (green). For the earlier epoch, more power is present in degrees 2 and 3 and a more rapid decrease in power is observed for the higher degrees, than in the comparison models. For the secular variation the prior is reproduced from degree three on at both epochs. For the earlier epoch, the dipole secular variation power is also close to the prior. The mean model shows less secular variation in the dipole than the comparison models, with more power in degrees 2-4. For the recent epoch, more variation is observed in the higher degrees with a more rapid decrease in power for the earlier epoch, similar to the spatial spectrum.   Figures 6 and 7 show local curves for Paris and Hawaii, respectively. Data from a surrounding of 250 km around the respective location is included with the prediction. Inclination and intensity are translated to the location of prediction making the simplifying assumption that the field is an axial dipole (Merrill et al., 1996). Declinations are taken as reported. The two locations were chosen because they have very different data coverage: Paris is covered well during recent times with a decrease in data from 1000 BCE on and virtually no data for epochs earlier than 6000 BCE. This is reflected in the prediction curves, which show less variation and increasing uncertainties for times with low data coverage. Hawaii is not as densely covered during recent times, but due to the volcanic area, records are available over the whole timespan of the model. Consequently, the predictions show variations during earlier times and the reported uncertainties are smaller. The comparison models agree within the reported 95%-intervals for both locations. For Paris, the SHA.DIF.14k model shows more variation during times earlier than 5000 BCE and most prominently from 12000 to 8000 BCE. For Hawaii, all models show a similar amount of variation, with SHA.DIF.14k varying slightly more and ARCH10k.1 slightly less, especially in the intensity. Two additional local predictions, for the Indian ocean and New Zealand, are provided with Figures S4 and S5 in Supporting Information S1.
We investigate the misfit of the model in Table 3. The sum of residuals squared, divided by the variance is a χ 2 -distributed random variable: It is evident, that  lies below the 95% confidence interval of the corresponding χ 2 distribution in almost all cases. Also  is below 1 for all subsets. Both results might indicate, that the model is overfitting the data, but when looking at the relatively big mean absolute error (MAE), this does not seem to be the case. Instead, we think the reason for the low misfit and  values lies in the large contributions from the dating uncertainties. These may be overestimated, as we use the reported uncertainties as standard deviations for normal error distributions. If in fact they come from a different distribution, for example, a uniform one, this procedure gives errors that are too large. The large errors also result in a low impact of such records on the model. The large MAE may be caused by those  (Merrill et al., 1996) to the location of prediction (orange dot) and shown as gray dots. Horizontal and vertical Gy bars indicate the two-sigma temporal and field component data uncertainties, respectively. The temporal distribution (top right) includes all data visible in the top left plot. records as well. As discussed above, a more thorough treatment of the error model (and with this also the dating uncertainties) in general may be necessary to address this.

Dipole Moment and Location
During the Holocene, the geomagnetic field is dipole dominated. Therefore, it is of special interest to infer the dynamics of the dipole. Figure 8 shows the evolution of the dipole moment. To access the dipole moment mean and standard deviation, sampling techniques are employed. The proposed model ArchKalmag14k shows significantly less variation in the dipole moment than comparable models. We observe some rapid variations from 1000 BCE to today, but for earlier times no rapid variations are found. Interestingly we observe a higher dipole moment than the comparison models for the interval 6000 to 2000 BCE and also from 12000 to 8000 BCE.  (Merrill et al., 1996) to the location of prediction (orange dot) and shown as gray dots. Horizontal and vertical Gy bars indicate the two-sigma temporal and field component data uncertainties, respectively. The temporal distribution (top right) includes all data visible in the top left plot. Figure 9 shows the latitude and longitude of the dipole location, together with the angular standard deviation (Butler, 2004). The latter is inferred via sampling. In earlier studies Schanner et al., 2021) we analyzed the statistics of the dipole axis coordinates directly. Here we analyze the projection of the dipole onto the sphere instead. The corresponding distribution is approximated by a von Mises-Fisher distribution and we report the latitude and longitude of its location parameter, instead of the mean of the marginal distributions. The advantage of performing statistics on the sphere instead of considering the marginal distribution is that there is no critical point (resp. meridian). The disadvantage is that the distribution is not available in closed form and that uncertainties cannot easily be translated to latitude and longitude, as approximations become unreliable when close to the pole (singularity in Equation 6). Similar to the dipole moment, the proposed model shows less variation during earlier times. The dipole latitude shows a trend opposite to the SHA.DIF.14k model for the interval 12000 to 6000 BCE, with the geomagnetic pole being very close to the geographic one in the beginning and a decrease in latitude toward recent times, in contrast to an increase present in the SHA.DIF.14k model. The angular standard deviation ( Figure 9, bottom row) increases toward earlier times, as is expected from the thinning data distribution.

South Atlantic Anomaly
To conclude the results, we present investigations of the South Atlantic Anomaly (SAA). The SAA is a region of low field intensity, that has been linked to reverse flux patches at the CMB during recent times (e.g., Terra-Nova et al., 2017). We compare the appearance and evolution of the SAA as predicted by ArchKalmag14k to other studies (Campuzano et al., 2019;Hartmann & Pacca, 2009). We do not follow the kernel-based approach of Terra-Nova et al. (2017), but investigate maps of the magnetic fields radial component at the CMB. In general, due to the projection into the Earth's interior, uncertainties at the CMB are so large that reverse flux in the mean model is not resolved reliably and more data and future work are required to confirm these findings. We consider the projections qualitatively nevertheless.
We find a region of field intensity lower than 32 μT emerging close to the tip of Brazil at 1200 CE (bottom right in Figure 11). Reverse flux is present to the north and a patch of reverse flux is located directly south of the region. Together with this patch, the region of low intensity rapidly moves south-eastward to the coast of today's Namibia, where it is located in 1300 CE (Figure 10b). This contrasts the findings of Campuzano et al. (2019), where the low intensity region emerges approximately 100 years earlier close to Madagascar, although an earlier emergence is within the uncertainties of our model. The SAA then extends to the West and slightly to the East, with the center drifting westward until 1500 CE, back to the origin of the region. From there it moves East and constricts at the coast of today's Namibia, almost disappearing at 1650 CE. This dynamic is also not present in SHA.WQ.2k by Campuzano et al. (2019), where the SAA persists at the coast of Namibia and does not decrease in size. The described evolution precedes the dynamics found by Hartmann and Pacca (2009). The subsequent westward drift of the low intensity region generally agrees with their findings and the findings of Campuzano et al. (2019) within the uncertainties.
Further, we find a low field intensity region emerging in 250 BCE west of today's Peru (bottom left in Figure 11). It drifts south-eastward and in 500 CE merges with a second low field intensity region that emerges around 400 CE North-East of Madagascar. Both anomalies are accompanied by reverse flux in the Southern hemisphere. The joint low intensity region continues to drift eastward and shrinks, persisting until 900 CE. Campuzano et al. (2019) find a low intensity field region emerging at the coast of Namibia at 175 CE. In their findings the Note. N refers to the number of records in the data subset.  gives the sum of the residuals squared, divided by the variances. This is a χ 2 -distributed random variable and we give the 95%-intervals of the corresponding distribution as 2 .  is the normalized misfit, that is, √  ∕ . MAE refers to the mean absolute error. The early subset consists of all records before −6000 BCE, the middle subset covers −6000 BCE to 0 CE and the recent subset includes all records after 0 CE. A map of the localized subsets is given with Figure S6 in Supporting Information S1.

Table 3
Data Misfit for Several Data Subsets earlier anomaly is static and grows until 500 CE. It then shrinks and disappears at 700 CE, earlier than in our findings.
Low intensity regions around the equator are present from the beginning of the model timespan on, but uncertainties are too large to reliably interpret their appearance. First reliable hints on a low intensity field region in the Indian ocean are present around 3000 BCE, with the region drifting eastward (top left in Figure 11) and a second low intensity region appearing over the Northern part of South America at 2600 BCE (top right in Figure 11). The anomaly in the Indian ocean disappears at 2200 BCE. The one above South America is accompanied by pronounced reverse flux, although during these epochs uncertainties at the CMB are even higher than during  Overall, the model shows low field intensity anomalies, accompanied by reverse flux, emerging and vanishing regularly, with a cycle in the order of 1000 years. An animation of the field at the Earth's surface and the CMB can be found with Supporting Information S1.

Discussion
In the preceding section, we proposed the new global geomagnetic field model ArchKalmag14k and presented its features. The local predictions give a reasonable representation of the underlying archeomagnetic data and agree with comparison models within the uncertainties. If no data is present, local curves show significantly less variation than the compared models. Low order, global scale degrees are only resolved if a sufficient amount of data is present. In this case, local predictions for remote locations also show rapid variations and uncertainties are relatively small (see the local predictions for the Indian ocean in Figure S4 in Supporting Information S1). If the data cannot resolve the global scales, the prior is reproduced, which is evident from local curves with no data coverage ( Figure 6) and the analysis of the dipole itself (Figures 8 and 9). For times earlier than 6000 BCE, the axial dipole varies only slightly around the prior mean value of approximately −36.19 μT (Figure 4, top row; see also Figure S3 in Supporting Information S1). Nevertheless, local variations are resolved, if supported by the data (Figure 7, especially the dip in declination at 11000 BCE). Spatial power spectra provide insight on the resolution of the model on global scales. From a comparison of the spectra to the respective prior it is evident, that for recent times information up to degree six is obtained, while for the earlier times the prior is reproduced already at degree 3 ( Figure 5, top row). An investigation of low intensity field regions reproduces the emergence and evolution of the South Atlantic Anomaly (SAA) in recent times (from 1600 CE on), while the preceding dynamics differ from other studies (Campuzano et al., 2019). Low intensity field regions can be resolved from 3000 BCE on. Although uncertainties at the CMB are large, hints for reverse flux patches associated with these field anomalies are found. A detailed evaluation relating these patches to the anomalies, for example, based on kernels (Terra-Nova et al., 2017) remains to be done and more data are needed to reduce the uncertainties.
In contrast to other recently proposed Bayesian models (Hellio & Gillet, 2018;Nilsson & Suttie, 2021), most prior parameters of ArchKalmag14k are inferred from the data via maximization of the log-marginal likelihood. As the marginal likelihood drops off quickly around the maximum, we did not perform an integration as proposed in the last study (Schanner et al., 2021). The a priori assumption of a constant axial dipole may lead to an underestimation of uncertainties in the dipole degrees, moment and location, as the prior mean is constrained well by data from recent times and variations during earlier times are considered around this fixed, constant value. Using Figure 11. Paths of the minima of low intensity field regions described in the text. Dots are drawn every hundred years. The locations of the minima have been obtained using grid search, which leads to several discrete jumps in the curves. only part of the recent records to create a data set that is more homogeneous in time may improve this, but leads to other complications as hyperparameters become less constrained and harder to determine, when fewer records are available. Artificially increasing the a priori dipole variance leads to more variation around the constant mean during earlier times, but also to higher posterior uncertainties and the model we propose lies well within these. Two scenarios are reasonable, to explain the absence of variations during earlier times in our model. Either the statistical properties (and thus the underlying processes) of the EMF changed during the Holocene, sometime around 3000 BCE. This is supported by a visual inspection of the top row in Figures 4 and 8. Or the data do not contain enough information to recover the global dynamics of the field, which is supported by the findings of the validation section. Additional data, for example, from sediments may help recovering the actual field dynamics, but require significant adaptation of the modeling method. Additionally, and similar to Hellio and Gillet (2018), a two parameter AR2-process may be considered for temporal correlations of the (axial) dipole.

Conclusions
This study proposes a new global geomagnetic model for the Holocene, called ArchKalmag14k. We modified the algorithms suggested in earlier works Schanner et al., 2021) to be applicable to the archeomagnetic database. The inversion is sequentialized by means of a Kalman-filter Kalman, 1960). The resulting model consists of sets of Gauss coefficients, secular variations and covariances, stored every 50 years. The model can be reproduced by code that is publicly available (https://sec23.git-pages. gfz-potsdam.de/korte/paleokalmag/) or is provided upon request. ArchKalmag14k can be imported by pymagglobal , so that feature analysis is straight-forward. Together with the software, we provide a Jupyter notebook, that illustrates how to reproduce ArchKalmag14k. The model is also available via a website: https://ionocovar.agnld.uni-potsdam.de/Kalmag/Archeo/.
The central result of this study is that for times earlier than 6000 BCE the current database of thermoremanent records alone does not contain enough information to construct global models. For times earlier than 6000 BCE, ArchKalmag14k reproduces the prior on a global scale and only local variations are resolved. Existing models may further overconfidently report variations during times later than 6000 BCE, as local variations that are resolved by higher degrees in ArchKalmag14k result in variations of the large-scale dipole in existing models. The next step is to extend and adapt the modeling framework to incorporate sediment records. As the recent study by Nilsson & Suttie (2021) shows this requires significant modifications due to aspects of the sedimentation process and the respective statistical implications.