%% This is a LaTeX template file for EPS
%%
%% modified by M.M. on 17 March 2014
%% modified by Y.O. on 7 September 2014
%% modified by Y.O. on 12 December, 2018
%% modified by Y.O. and M.M. on 6 December, 2019
%% modified by Y.O. on 28 March, 2020
%% modified by T.S. on 2 February, 2022

\documentclass{EPS}

\usepackage{mathtools, amsfonts, physics, amsmath}
\usepackage{lineno}
\linenumbers*[1]
\setcounter{page}{1}

% The lineno package requires math environments to be wrapped
\let\oldequation\equation
\let\oldendequation\endequation
\renewenvironment{equation}
  {\linenomathNonumbers\oldequation}
  {\oldendequation\endlinenomath}
\let\oldalign\align
\let\oldendalign\endalign
\renewenvironment{align}
  {\linenomathNonumbers\oldalign}
  {\oldendalign\endlinenomath}

\title{The global geomagnetic field over the historical era: What can we learn from ship-log declinations?}
\author{
Maximilian Schanner, Institute of applied mathematics,
Potsdam University, 14467 Potsdam, Germany,
arthus@gfz-potsdam.de\\
Lukas Bohsung, German Research Centre for Geosciences (GFZ) Potsdam, lbohsung@gfz-potsdam.de\\
Clara Fischer, GFZ Potsdam, clara\_fischer@gmx.de\\
Monika Korte, GFZ Potsdam, monika@gfz-potsdam.de\\
Matthias Holschneider, Potsdam University, hols@math.uni-potsdam.de
}

\abstract{
Modern geomagnetic field models are constructed from satellite and observatory data, while models on the millenial timescale are constructed from indirect records of thermoremanent and sedimentary origin.
An intermediate period, spanning the last four centuries, is covered by historical survey data and ship-logs, which is strongly dominated by geomagnetic declination information.
We apply a sequentialized, Gaussian process based modeling technique to this dataset and propose a new field model for this era.
In order to investigate the information gained from declination records from ship-logs, we seperate the dataset and construct a second model, where unpaired declination records are removed.
Instead of uncovering a more detailed field, the availability of more records leads to better resolution of the global field structure.
The availability of more records helps notably to constrain global field properties like the dipole moment.
It also allows to resovle some detailed field structures more accurately.
Based on the model constructed from the full dataset, we perfom an analysis of the South Atlantic Anomaly and regions of low field intensity in general.
We extend a recent analysis of center of mass movement and area evolution of the South Atlantic Anomaly further back in time and confirm the findings of its non-monotonous growth.
}

\keywords{Geomagnetic field, Statistical modeling, Gaussian processes, Kalman filter, South Atlantic Anomaly, Bayesian inversion}

\begin{document}

\maketitle

\section{Introduction}

Studies of the geomagnetic field can shed light on the Earth's interior processes.
By projecting global field models to the core-mantle boundary, the evolution of the geomagnetic field can be related to the core flow \citep{Bloxham1992, Hulot2002, Gillet2019}.
Both the geomagnetic field and the dynamo process in the liquid core exhibit a rich dynamic on multiple timescales \citep[e.g.][]{Constable2015}.
A good worldwide data coverage is required to provide the full global view of geomagnetic field evolution.
During recent times, satellite data is available and allows for the inversion of high-resolution models on the decadal timescale \citep[e.g.][]{Finlay2020, IGRF12, Baerenzung2022}.
When going back in time, data becomes sparse and the model resolution is limited \citep{Hellio2018, Schanner2022}.
Still, archeomagnetic records allow for the reconstruction of the global geomagnetic field at least until 6000 BCE and together with sediment records, large scale features of the field can be recovered on millenial and longer timescales \citep{Constable2015}.
An intermediate timescale, going back several hundred years, is covered by historical records from land surveys and ship-logs.
While the former include directional and intensity measurements, the latter almost exclusively comprises declination records, that have been measured in large numbers for navigational purposes.
Even though the magnetic declination contains limited information about the field vector, the ship-log dataset may provide useful information on global and regional field structure, due to its dense coverage of the oceans.


Traditionally, global geomagnetic field models are represented in the spherical harmonics expansion with a B-spline model for the time evolution.
The series coefficients are determined by a regularized least squares procedure \citep{Bloxham1992, Korte2003}.
During recent times, statistical methods have been suggested, allowing for assessment of model uncertainties \citep{Leonhardt2007, Holschneider2016, Mauerberger2020}.
While several bootstrapping approaches exist \citep{Korte2011, Hellio2018, Arneitz2019}, more rigorous algorithms have been proposed recently \citep{Schanner2021, Nilsson2021}.
Based on Gaussian process regression, the non-linear relation of directional and intensity records is either tackled via sampling or linearization.
Both methods are limited by high numerical costs when considering large datasets.


In this study we apply a sequentialized Bayesian inversion procedure \citep{Schanner2022}  to a combined dataset of archeomagnetic and historical records to construct a global field model for the past 1000 years.
The main difficulty is the estimation of model parameters, as this step requires inverting the dataset multiple times.
We argue, that using hyperparameters which are estimated from a longer timescale archeomagnetic dataset is reasonable and enables the application of the rigorous inversion to the larger combined dataset.
By separating the ship-log declinations and performing two inversions, we illustrate what information can be inferred from declination records on a global scale.

\citet{Arneitz2021} present a model for a similar era, based on a similar database.
Major differences to their work are in data selection and the inversion procedure.
We compare our results to their BIGMUDIh.1 model whenever appropriate.

\section{Methods}
We construct two new global geomagnetic field models, covering the last 1000 years.
The modeling procedure is described in detail in \citet{Schanner2022}.
The central idea is modeling the global geomagnetic field as a Gaussian process, i.e.
\begin{equation}
\boldsymbol{B} \sim \mathcal{GP}\big(\bar{\boldsymbol{B}},~K_{\boldsymbol{B}}\big)~.\label{eq:GP}
\end{equation}
As historical and archeomagnetic records are given as directions and intensities, the field vector $\boldsymbol{B}$ is non-linearly related to the observations.
To be able to access the posterior, the observation functionals are linearized around a proxy model, resulting in a normal likelihood and therefore a normal approximation of the posterior distribution, thanks to a Gaussian prior.
The historical dataset also contains measurements of the horizontal intensity $F_H$.
These are handled analog to the intensity $F$, i.e.
\begin{equation}
F_H \approx
    \frac{1}{\tilde F_H} \qty(\tilde{B}_N,\, \tilde{B}_E,\,0)^\top\cdot\boldsymbol B ~.
\end{equation}
Due to the amount of data, a direct Gaussian process regression is unfeasible.
Instead, the inversion is sequentialized by means of a Kalman-filter \citep{Kalman1960, Baerenzung2020}.
Cross correlations are reintroduced by a smoothing step.

Dating errors are handled similar to the ArchKalmag14k-approach by using a noisy input Gaussian process \citep{McHutchon2011}.
However, due to the shorter interval covered by the model, we do not expect very large field changes over the interval spanned by the dating uncertainties.
We therefore set upper bounds for the translated errors, in constrast to the original approach, where translated errors could become arbitrary large.
The upper bounds are 30$^\circ$, 15$^\circ$, 5$\mu$T and 3.5$\mu$T for declination, inclination, intensity and horizontal intensity respectively.

We start the modeling process at 1950 CE and go back until 1000 CE.
New data is incorporated every year and the output stored every ten years.
Similar to ArchKalmag14k, the cutoff degree (in a spherical harmonics expansion) is set to 20.
This leads to an output of 440 field coefficients and 440 secular variation coefficients, together with the corresponding covariance matrix, at 96 knot points.

With the Kalman filter approach, it is straightforward to connect the model to existing ones.
We choose the initial mean and covariance to agree with the Kalmag model output for 1950 CE \citep[Kalmag spans the interval 1900 to 2022;][]{Baerenzung2022}.
Besides constraining the model by incorporating knowledge from the satellite era, this has the advantage of improving the initial linearization point.


\subsection{Hyperparameters}
With the established algorithm, the model depends on several hyperparameters, that define the Gaussian process kernel and prior mean.
These consist of two timescales, that give the temporal smoothness of the process, two variances, which control the scale of the field dynamics, a mean value for the axial dipole and a residual term, that is related to contributions in the measurements, which are not reflected in the model (e.g.\ the crustal field).
To estimate these parameters, two central strategies have been pursued in the past:
\citet{Hellio2018} and \citet{Nilsson2021} used reference models from the satellite era to constrain the parameters, while \citet{Schanner2022} estimated the parameters from the dataset directly.
Applying the latter approach to the historical dataset is unfeasible, due to the large amount of data (about 140.000 declination records) and the resulting longer inversion times.
Inversion of the full archeomagnetic dataset, used to construct the ArchKalmag14k model, takes about 30 seconds, but the inversion of the full combined dataset takes about two and a half hours.
In order to find the optimal set of hyperparameters, around ten thousand inversions have to be performed, rendering the strategy unfeasible.
Another hindrance is the timescale on which the dipole changes.
Recently, \citet{Nilsson2022} found recurrent signals in the dipole moment with a period of 650 years.
To reliably estimate a correlation time corresponding to these signals, multiple periods should be covered by the model.
With the 1000 years covered by the dataset of this study, this is barely the case and a longer timespan will facilitate estimation of the a\,priori timescale.
We therefore choose to set the a\,priori hyperparamters to the values estimated from the ArchKalmag14k.r dataset, as given in Table 1 of \citet{SchannerPrep}.
In order to check for faster signals in the dipole, we performed an inversion with a white noise kernel (for the dipole), which amounts to neglect the time correlation in the prior for the dipole.
Analysis of the posterior mean and empirical autocorrelation function indicates that the ArchKalmag14k.r prior is able to capture the variations present in the data.
We further inverted the data using the two parameter kernel proposed by \citet{Bouligand2016}.
As this type of kernel contains two timescales with different spectral behaviors, it might be able to reflect both the millenial and the faster, centennial dynamics of the dipole.
Two sets of parameters were chosen:
For one model, all hyperparameters were kept from ArchKalmag14k.r and the dipole correlation times were chosen as proposed by \citet{Hellio2018}.
The resulting model is not able to resolve the variations present in the data, most likely due to the large value of the long timescale ($T_s$, c.f. \citet{Hellio2018}, above paragraph 3).
For the second model, the hyperparameters were estimated from 10 percent of the data, that have been chosen randomly.
The obtained paramaters for the non-dipole agree well with ArchKalmag14k.r.
The dipole parameters lead to a model that reproduces the posterior features of the one we get when inverting with the ArchKalmag14k.r-parameters.
This confirms our concern that 1000 years is not enough to determine the suitable dipole hyperparameters and supports our decision to use the a\,priori hyperparamters estimated from the longer archeomagnetic dataset.


\section{Data}
There are two classes of data considered in this study.
Indirect measurements of the geomagnetic field are taken from the GEOMAGIA database \citep{Brown2015}.
The dataset consists of declination, inclination and intensity measurements from volcanic rocks and archeologic artefacts.
Direct measurements stem from the HISTMAG database \citep{Arneitz2017}.
In this database, records from land surveys, ship voyages and observatories are compiled.
In addition to declination, inclination and intensity records, some historical horizontal intensity values are included.
In general, the dataset is quite similar to the one used for constructing the BIGMUDIh.1 model \citep{Arneitz2021}.
The HISTMAG database does not contain error estimates for most of the records.
We assign the rather conservative values of 4.5$^\circ$ to directions and 8.25 $\mu$T to intensities.

The indirect measurements are distributed equally over the time interval under consideration.
Spatially, they are unevenly distributed, with clusters in Europe, East Asia and America (see figure S1 in the supplementary material).
From the direct measurements, we seperate records that consist of unpaired declinations.
The resulting historical dataset with single declinations removed covers the last 400 years, with the majority of records in the nineteenth and twentieth century.
The spatial coverage is more even than in the archeomagnetic dataset.
Still, a bias towards the northern hemisphere remains.
An illustration of the dataset without declinations, showing composition, spatial and temporal distribution, is given in figure \ref{fig:data}.
The declination dataset goes further back in time, until the 16th century.
Except from Antarctica and some inland regions, it covers the globe densely, as declination was measured both from shipboard and on land surveys.
Further illustrations of the datasets, similar to figure \ref{fig:data}, are provided with the supplementary material.

Following \citet{Schanner2022}, we identify outliers by means of a Naive Bayes classifier \citep[e.g.][]{Berrar2018}.
Records that have already been identified as outliers during construction of the ArchKalmag14k.r model are removed first.
Then, an inversion is performed.
For every record, the probability to be generated from the model distribution or from a flat (``noise'') distribution is caluclated.
Records that are more likely to stem from the flat distribution are discarded as outliers.
This leads to the exclusion of 169 records.
The small number of rejected records is likely due to the conservative error values assigned to the historical records, but also demonstrates a good internal consistency of the large historical dataset.
The distribution of outliers mostly reflects the distribution of the data.
The majority of rejected records are declinations, for which seemingly a minus sign has been lost during transcription.
Comparison of the model with the full database and models without outliers shows only minor differences, especially on a global scale.
We still reject the records, as we believe the approach with the naive Bayes classifier to be more objective than previous approaches, where data was rejected according to mostly subjective quality criteria \citep[e.g.][]{Arneitz2021}.
We give the number of records ending up in our dataset in table \ref{tab:data}.
For comparison purposes, also the number of indirect records is given.

\begin{table}[h!bt]
    \centering
    \caption{Number of records of individual field components in several datasets. ``Indir.'' refers to indirect measurements and ``no\_D'' to the dataset with single declinations removed. ``Full'' refers to the full dataset, which the HistKalmag model is built from.}
 	\label{tab:data}
    \begin{tabular}{r|ccccc}
    Name & \# $D$ & \# $I$ & \# $F$ & \# $H$ & \# tot\\\hline
    indir.\  &1,945 & 3,102 & 1,477 & 0 & 6,524\\
    no\_D &15,564 & 22,934 & 5,413 & 11,935 & 55,846\\
    full &156,471 & 22,934 & 5,413 & 11,935 & 196,753\\
    \end{tabular}
\end{table}


\section{Results and Discussion}
In analogy to the Kalmag and ArchKalmag14k models, we call the presented models HistKalmag and Histkalmag.no\_D, where the latter refers to the version without declinations.
Figure \ref{fig:DPvND} shows the dipole and non-dipole energy of several models at the Earth's surface.
Both versions of HistKalmag show less variation than the comparison model BIGMUDIh.1 in the dipole, while more variation is present in the non-dipole coefficients.
Before 1800 CE, when no survey data is available, the non-dipole energy of HistKalmag.no\_D quickly rises to a level comparable to the ArchKalmag14k.r model.
In contrast, the non-dipole contribution in the HistKalmag model rises more uniformly and reaches a similar level at around 1500 CE, when the database comprises archeological and volcanic records only.
An explanation for this may be the global information contained in the declination records.
While both ArchKalmag14k.r and HistKalmag.no\_D contain limited global information and thus explain observations locally, by higher spherical harmonic degrees, the dense coverage of the oceans in the 16th century leads to a more reliable recovery of the dipole in HistKalmag and thus less energy in the higher order degrees.
The deviance between ArchKalmag14k.r and the HistKalmag models prior to 1400 CE, where the database is the same, is due to the slight modification of the noisy input Gaussian process (i.e.\ the constraint translated dating errors) mentioned above.

Figure \ref{fig:regional} shows the global field intensity and standard deviation for the epoch 1700 CE.
The maps are centered at the Pacific, as this region is only covered densely by the ship-log declination data.
Clearly, uncertainties are lowest in the HistKalmag model, which is based on the full dataset.
This model also shows lower intensity around the South Pole and a more constrained South Atlantic Anomaly.
Another difference is the low intensity field patch over South East Asia.
Location, shape and intensity of this patch are quite different among the three models.
In the ArchKalmag14k.r model it is merged with the South Atlantic Anomaly.
BIGMUDIh.1 also shows this patch, but it is merged with the South Atlantic Anomaly as well, separating later, around 1760 CE.


We show local model predictions, together with local data, in figures \ref{fig:india} and \ref{fig:rapa}.
Another location is provided with the supplementary material.
In the 17th century, the Indian ocean was crossed by many ship voyages.
The collected declination data leads to significant differences in the local predictions from the HistKalmag.no\_D and the HistKalmag model.
Interestingly, this difference is not in the declination itself, but in the local predictions of the field intensity (figure \ref{fig:india}, top row).
We believe this is again a consequence of global vs.\ local information.
The HistKalmag.no\_D model contains survey data from India as well as inclination data from some ship voyages across the Indian ocean.
These records contain intensity variations that are resolved locally by the model.
Only the global coverage, introduced by the ship-log declination dataset, relates these features to global field structure, that also results in a locally different field intensity.
Similar differences are visible in local predictions at the Azores (supporting information, figure S3).
However, the predictions there show some deviance for the declination as well.
This may be due to the Azores being closer to Europe, which is densely covered by survey- and other data from the 18th century, and by data from ship voyages across the atlantic. 
The different global structure is also evident from figure \ref{fig:DPvND} and figure \ref{fig:regional}.

Even though no data is present in the direct surrounding, predictions at Rapa Iti, a small island in the Pacific, show a declination change around 1400 CE.
This change is related to the appearance of a low intensity field patch in the Pacific around this time (figure \ref{fig:1400patch}).
The signal of this patch is recorded in intensity records from the Pacific and North America and declination records from South America.


\subsection{South Atlantic Anomaly}
The South Atlantic Anomaly, a low intensity field region in the South Atlantic region, has been investigated by several studies \cite[e.g.][]{Hartmann2009, PavonCarrasco2016, TerraNova2017, Campuzano2019, Finlay2020}.
Recently, \citet{Hamit2021} proposed a novel definition of the South Atlantic Anomaly region and center of mass, taking global field changes into account.
We consider the proposed definition and estimate area and center of mass of the South Atlantic Anomaly accordingly, extending the analysis of \citet{Hamit2021} further back in time.

The HistKalmag mean model shows a low intensity field region over South East Asia in the beginning. This region splits, with one part quickly moving westward and the other decreasing and disappearing around 1100 CE.
The westward moving part stops its movement slightly north of South America and simliarly splits, with one part decreasing and the other moving eastward.
In the mean model, today's South Atlantic Anomaly appears as a merging of the decreasing part noth of South America, that slightly moves eastward, and a low intensity field patch emerging close to Madagascar around 1200 CE.
However, individual realisations from the Histkalmag distribution show vastly different field configurations before ca. 1350 CE.
The short lived low intesity field patch depicted in figure \ref{fig:1400patch} is the first feature that is consistent within the ensemble.
Figure \ref{fig:SAA_area} shows the area of the low intensity field region as defined by \citet{Hamit2021} (using their terminology, herafter to referred as S1 region).
The mentioned differing field configurations also reflect in the S1 area, as evident from the spread in the sample curves depicted.
After 1600 CE, the ensemble shows less variability around the mean, due to the increase in the number of data.
This is the reason why we end the plot of the S1 center of mass (figure \ref{fig:SAA_com}) at 1500 CE.
Another reason is the aforementioned low intensity field patch in the western Pacific.
With multiple patches, the center of mass of the whole region jumps abruptly and is not comparable to earlier epochs.
To consistently track the center of mass, one would have to isolate the different patches, which is beyond the scope of this article.
Evolution of both area and center of mass of the S1 region of the HistKalmag model are quite similar to BIGMUDIh.1, disagreeing most significantly before 1600 CE.
The area decreases after 1500 CE, with two minima at 1680 CE and 1770 CE.
The decrease is linked to the disappearance of a reverse flux patch at the core mantle boundary at the southern tip of Africa.
The rapid westward movement of the center of mass after 1700 CE is likely caused by the emergence of reverse flux patches in the western and central Pacific.
Earlier, the center of mass is more stable as the associated flux patch at the core mantle boundary moves only slightly.
The described features are present in the mean model and most of the ensemble members, while fluctuations around the mean model are higher at the core mantel boundary than at the surface.

\section{Conclusions}
The HistKalmag model constitutes a bridge between the longer timescale ArchKalmag14k model and the satellite- and observatory era based Kalmag model.
The large database of direct, non-linear observations of the magnetic field from surveys and ship-logs provides global field information over the last five centuries.
In order to assess what information can be extracted from declination records, an alternative model, HistKalmag.no\_D, was constructed, where single declination records were disregarded.
The model shows a higher non-dipole energy for the era covered by the declination dataset than the HistKalmag model, which is built from the full dataset.
The availability of more records helps notably to constrain global field properties like the dipole moment.
A similar behavior is observed over longer timescales in the ArchKalmag14k model.
More records also allow to resovle some detailed field structures more accurately.
Still, the general field structure is already captured well by the archeomagnetic database.
This is likely due to the field sources lying deep inside the Earth and the resulting suppression of higher spherical harmonics.
When looking at maps of the radial component at the core mantle boundary (supplemenatry material, fig. S4), one might get the impression that for recent times, when more data is available, the field shows more small scale features.
This is true for the mean field, however, samples from the posterior show a constant resolution (i.e.~characteristic size of the posterior fluctuations) over the whole model period.
The spread in the distribution is bigger for earlier times, leading to smoother (and thus bigger) patches in the mean.
This highlights the importance of considering uncertainties when discussing model properties and features.

For the period covered by direct observations, we linked the South Atlantic Anomaly and another low intensity field patch to reverse flux at the core mantle boundary.
For earlier times, uncertainties are too high to reliably locate features at the core mantle boundary.
Similar to \citet{Hamit2021}, we find that the growth of low intensity field patches is related to growth of reverse flux patches at the core mantle boundary.
The reverse flux patches at the core mantle boundary are relatively static and instead of following their movement, the low intensity field regions at the Earth's surface move from one patch to another
(most notably, the South Atlantic Anomaly is first located above a reverse flux patch at the southern tip of Africa and then ``dragged'' westward, by reverse flux appearing east to the coast of South America).
By estimating the area of the South Atlantic Anomaly further back in time, we could extend the analysis of \citet{Hamit2021} and confirm their finding of a non-monotonous growth.
More extremly, the South Atlantic Anomaly shrinks from 1500 CE to ca.\ 1680 CE and starts to grow rapidly after 1770 CE, with milder size changes in between.
For earlier times, the area of low intensity is also affected by a low intensity patch in the western Pacific and a detailed seperation of the South Atlantic Anomaly from other low intensity field regions remains open for these times.

\section{Declarations}

\section{Availability of data and materials}
The dataset is a combination of the GEOMAGIA \citep{Brown2015} and HISTMAG \citep{Arneitz2017} databases, which are freely available. 
A list of records identified as outliers is provided upon request. 
The modeling software is an extension of paleokalmag, wich is available via a git-repository \texttt{https://sec23.git-pages.gfz-potsdam.de/korte/paleokalmag/}.
The HistKalmag model is also distributed via a dedicated website at \texttt{https://ionocovar.agnld.uni-potsdam.de/Kalmag/Histo/}.

\section{Competing interests}

The authors declare that they have no competing interests.

\section{Funding}
This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), grant 388291411.

\section{Authors' contributions}

M. Schanner, M. Korte and C. Fischer designed this research.
M. Korte and C. Fischer performed data preprocessing and selection.
M. Schanner performed theoretical work, analysis and software development.
L. Bohsung assisted in analysing the impact of the hyperparameters on the model.
M. Schanner assembled the manuscript, with contributions from all co-authors.
M. Korte and M. Holschneider supervised the findings of this work.

\section{Acknowledgements}
The authors thank J. Baerenzung for his assistance with making the model available via the mentioned website.

\bibliographystyle{plainnat}
\bibliography{refs}

\section{Figures}
\begin{figure}[h!tb]
    \centering
    \includegraphics[width=\textwidth]{data_no_D.pdf}
    \caption{Composition of the dataset with single declinations removed, together with spatial and temporal distribution.}
    \label{fig:data}
\end{figure}

\begin{figure}[h!tb]
    \centering
    \includegraphics[width=\textwidth]{DP_v_ND.pdf}
    \caption{Dipole (top) and non-Dipole (bottom) energy at the Earth's surface. For ArchKalmag14k.r and the HistKalmag models, samples from the posterior are drawn transparently in the background to illustrate the uncertainties. The thick lines give the ensemble mean, while the dashed lines represent the energy calculated from the mean model directly.}
    \label{fig:DPvND}
\end{figure}

\begin{figure}[h!tb]
    \centering
    \includegraphics[width=\textwidth]{regional_1700.pdf}
    \caption{Geomagnetic field intensity (top) and standard deviation (bottom) for the epoch 1700 CE. Depicted are the three models ArchKalmag14k, HistKalmag.no\_D and HistKalmag. The maps are centered at the Pacific, as this region is covered densely only by the ship-log declination data.}
    \label{fig:regional}
\end{figure}


\begin{figure}[h!tb]
    \centering
    \includegraphics[width=\textwidth]{indian_ocean.pdf}
    \caption{
        Local predictions of three different models in the Indian ocean (-6.5$^\circ$, 73$^\circ$), 
        together with spatial and temporal distribution of the surrounding data.
        The upper right panel contains all records from the spatial distribution, while only data from a 500 km radius (depicted in orange in the top left panel) is shown together with the local predictions.
        The errorbars reflect one standard deviation.
        Declination and intensity are translated along the corresponding axial diple \citep{Merrill1996}.
        Predictions from the posterior distribution are drawn as transparent lines in the background, to illustrate the model uncertainties.
        }
    \label{fig:india}
\end{figure}

\begin{figure}[h!tb]
    \centering
    \includegraphics[width=\textwidth]{rapa_iti.pdf}
    \caption{
        Local predictions of three different models at Rapa Iti (-27.605556$^\circ$, -144.344444$^\circ$), 
        together with spatial and temporal distribution of the surrounding data.
        The upper right panel contains all records from the spatial distribution, while only data from a 500 km radius (depicted in orange in the top left panel) is shown together with the local predictions.
        The errorbars reflect one standard deviation.
        Declination and intensity are translated along the corresponding axial diple \citep{Merrill1996}.
        Predictions from the posterior distribution are drawn as transparent lines in the background, to illustrate the model uncertainties.
        }
    \label{fig:rapa}
\end{figure}

\begin{figure}[h!tb]
    \centering
    \includegraphics[width=\textwidth]{minimum_map.pdf}
    \caption{Evolution of a short lived low intensity field patch in the western Pacific around the year 1400 CE, calculated from the HistKalmag model. The top row shows the field intensity at the Earth's surface. The second row the field intensity standard deviation. The third row depicts the radial field component at the core mantle boundary and the bottom row the corresponding standard deviation. The yellow contour line indicates the S1 low field intensity region, according to the definition by \citet{Hamit2021}.}
    \label{fig:1400patch}
\end{figure}

\begin{figure}[h!tb]
    \centering
    \includegraphics[width=\textwidth]{SAA_area.pdf}
    \caption{Area of the S1 low intensity field region over time, according to the definition by \citet{Hamit2021}. The thick blue line corresponds to the HistKalmag mean model, the transparent lines in the background show the area of 1000 samples from the posterior distribution.}
    \label{fig:SAA_area}
\end{figure}

\begin{figure}[h!tb]
    \centering
    \includegraphics[width=\textwidth]{SAA_paths.pdf}
    \caption{Location of the South Atlantic Anomaly center of mass, according to the definition by \citet{Hamit2021}. The thick blue line corresponds to the HistKalmag mean model, the transparent lines in the background show the location of 1000 samples from the posterior distribution.}
    \label{fig:SAA_com}
\end{figure}


\end{document}