Estimating post‐Depositional Detrital Remanent Magnetization (pDRM) Effects: A Flexible Lock‐In Function Approach

The primary data sources for reconstructing the geomagnetic field of the past millennia are archeomagnetic and sedimentary paleomagnetic data. Sediment records, in particular, are crucial in extending the temporal and spatial coverage of global geomagnetic field models, especially when archeomagnetic data are sparse. The exact process on how sediment data acquire magnetization including post‐depositional detrital remanent magnetization is still poorly understood. However, it is widely accepted that these effects lead to a smoothing of the magnetic signal and offsets with respect to the sediment age. They impede the direct inclusion of sediment records in global geomagnetic field modeling. As a first step, we model these effects for a single sediment core using a new class of flexible parameterized lock‐in functions. The parameters of the lock‐in function are estimated by the maximum likelihood method using archeomagnetic data as a reference. The effectiveness of the proposed method is evaluated through synthetic tests. Our results demonstrate that the proposed method is capable of estimating the parameters associated with the distortion caused by the lock‐in process.


10.1029/2023JB027373
2 of 17 While the lock-in process in the TRM occurs on short time scales (hours to weeks), the lock-in time of magnetic moments in sediment records can be much longer (years to centuries).The magnetization in sediments is called detrital remanent magnetization (DRM), which was first measured by McNish and Johnson (1938).During the sedimentation process, magnetic particles are deposited in such a way that their magnetic moments tend to point in the direction of the geomagnetic field while interaction with other particles and the ongoing solidification increasingly impede the particles to fully align.Additional sediment particles lead to a consolidation of the underlying layers and thus to a mechanical lock-in of the magnetic particles.The magnetization in sediments is affected by the interaction of the magnetic particles with the substrate at the sediment-water interface and by dewatering of the sediment (Irving, 1957).The terminology and classification of these effects are not completely consistent in the literature.In the following, we will be using the terminology recommended in the review paper by Verosub (1977).According to Verosub (1977), the term DRM refers to the remanent magnetization found in sediments.By depositional DRM (dDRM) we describe the magnetization acquired by the interaction of the particles with the substrate at the sediment/water interface.The term post-depositional DRM (pDRM) refers to the longer timescale and describes any magnetization that is acquired after the particle settles on the sediment/ water interface.
There are various effects that are summarized in the term dDRM.One example is the inclination error, which occurs when non-spherical particles settle flat on the sediment/water interface.This leads to a distortion of the inclination to smaller values (R. King, 1955).Another distortion of the inclination can occur when aligned particles roll into the nearest depression of the sediment/water interface (Griffiths et al., 1960).
In this paper, we will focus on the investigation of the post-depositional DRM.The traditional pDRM model, based on decades of research (e.g., Hamano, 1980;Irving, 1957;Irving & Major, 1964;Kent, 1973;Otofuji & Sasajima, 1981), can be described as follows.In general, only coarse-grained fractions are mechanically fixed more or less immediately after deposition.Smaller particles which are embedded in water-filled voids or pore spaces of the sediment can move freely for a longer period of time.With progressive consolidation and dewatering of the sediment, also these particles become locked in.However, alternative theories have challenged the classical pDRM acquisition concept, suggesting that sediment flocculation restricts significant post-depositional grain movement within pore spaces (Katari et al., 2000).Also, the effects and potential roles of bioturbation have been investigated and resulted in alternative sediment mixing models (e.g., Egli & Zhao, 2015).
To sum up, how exactly the pDRM process works and what effects are more or less important is still not fully understood.Nonetheless, decades of investigation, several experiments, and modeling approaches lead to the conclusion that the pDRM process results in a delayed and smoothed signal.In other words, the magnetic moment of an entire layer is given by the weighted sum of the geomagnetic field over the lock-in period.The weights are given by so called lock-in functions (Roberts & Winklhofer, 2004;Suganuma et al., 2011).
Figure 1 provides a practical illustration of this effect.We show declination and inclination data obtained from a real sediment record from Sweden, together with measurement errors.Additionally, we show the predictions from the ArchKalmag14k.rmodel (Schanner et al., 2022).Note that ArchKalmag14k.r exclusively relies on archeomagnetic data.The visual representation clearly demonstrates a noticeable delay and mild smoothing in the sediment data compared to the predictions from the ArchKalmag14k.rmodel.It's important to highlight that Over the last decades, many lock-in functions have been suggested.Exponential lock-in functions (e.g., Kent & Schneider, 1995;Løvlie, 1976), constant (e.g., Bleil & Von Dobeneck, 1999), linear (e.g., Meynadier & Valet, 1996), cubic (e.g., Roberts & Winklhofer, 2004), Gaussian (e.g., Suganuma et al., 2011) and parameterized lock-in functions that can cover multiple classes (e.g., Nilsson et al., 2018).
In this study, we present a new class of parameterized lock-in functions designed for the comprehensive modeling of the offset and smoothing effects caused by the pDRM process.These lock-in functions are characterized by four parameters giving them the flexibility to approximate a wide range of previously suggested lock-in functions and therefore a wide range of possible lock-in behaviors.
The cornerstone of our methodology lies in the estimation of the lock-in function parameters for a given core sample.The aim is to identify the most suitable lock-in function for a given core sample, one that accomplishes the tasks of shifting the data adequately and removing the right degree of smoothing.The determination of this optimal lock-in function necessitates leveraging external geomagnetic field information that is not affected by distortions associated with the lock-in process.This is where archeomagnetic data plays a pivotal role.Archeomagnetic data acquire magnetization through TRM and are therefore independent of any lock-in distortion.
Our model uses maximum likelihood methods and global archeomagnetic data to estimate the lock-in function parameters.To illustrate this concept, Figure 1 provides a visual representation.It is essential to emphasize that our model does not rely on the pre-existing Archkalmag14k.r model or any similar models.Instead, it exclusively utilizes the original data used in the creation of the Archkalmag14k.rmodel (Schanner et al., 2022).
In Section 2, we first briefly outline the geomagnetic field modeling method and then introduce the new class of lock-in functions.Subsequently, we delve into two technical subsections (Sections 2.3 and 2.4).In Section 2.3, we describe the relationship between the sediment observations and the geomagnetic field resulting in a rigorously defined data model for the sediment observations.Within Section 2.4, we describe the adoption and adjustments of the Kalman Filter-based method employed in Schanner et al. (2022) to suit our context.While these two subsections are not essential for a basic grasp of the results, they are crucial for a deeper understanding of the methodology and for making it easier to replicate the research.The newly developed method undergoes extensive testing on multiple synthetic data sets in Section 3. We discuss the findings and give an outlook on future work, including the prospective application of this methodology to real sediment records in Section 4 before ending with a summary, conclusions, and outlook on future research.

Geomagnetic Field Model
We will model the geomagnetic field by using a Bayesian approach based on Gaussian Processes.Every Gaussian Process is uniquely defined by a mean and a covariance function (Rasmussen, 2004).
We describe the geomagnetic field as the realization of a Gaussian Process denotes the standard two-sphere associated to the space variable.Therefore, the knowledge about the geomagnetic field and its uncertainty is a distribution of functions   ∶  2 × ℝ → ℝ 3 .In the following, we will model the lock-in process for a single sediment core sample and treat the space variable as a constant, that is, we will consider   as a Gaussian process of time only.We follow the a priori assumptions of Schanner et al. (2022) and use the estimated hyperparameters given in Table 2 of Schanner et al. (2022).Hence, we assume that all Gauss coefficients are a priori uncorrelated at a reference radius   = 2,800 km with zero mean except for the axial dipole.For the axial dipole, we assume a constant mean value of   0 1 = −38 μT (at the Earth's surface).Further, we assume an a priori variance  DP = 39 μT for the dipole and an a priori variance  ND = 118.22μT for all higher degrees (at the reference radius).The temporal correlation of the Gauss coefficients is assumed to be .
These parameters reflect the statistics of the archeomagnetic data and a physical interpretation is not obvious.For example, the value of −38 μT for   0 1 is the optimal value when fitting an axial dipole to the data directly.The correlation times may relate to physical processes, but are derived from the variability that is resolved in the data.Using other prior parameters is straightforward and an investigation on the influence of those parameters will be done in a future study.

Lock-In Process
The pDRM or lock-in process is affected by several effects which are still not completely understood.However, it is widely accepted that the pDRM process leads to a delayed and smoothed signal in the magnetic moment of a layer.Roberts and Winklhofer (2004) describe this behavior by the following convolution or weighted average where  () describes the magnetic moment of the layer at depth   ,   describes the geomagnetic field (see Section 2.1, note that we use the depth-series here derived from an age-depth model) and   gives the lock-in depth, that is, the depth where all particles are completely consolidated.The weights are given by the lock-in function   .
The lock-in function does not depend on the depth   of the layer.The underlying assumption here is that the sedimentation material does not change significantly over the absolute depth of the sediment record.Alternatively, we could formulate our model in time instead of depth and assume a constant sedimentation rate, but this is a much stricter assumption.
For a given core sample, we are interested in the weights, given by the shape of the lock-in function   .We propose a new class of piecewise linear parameterized lock-in functions Depending on the four parameters  1, 2, 3, 4 ∈ ℝ≥0 with 4 can model the offset as well as the smoothing related to the lock-in process.The offset is given by the half lock-in depth and the smoothing by the width and height of the lock-in function (see Section 3 for a detailed discussion).
For a given core sample, we estimate the parameters  1-4 using a Bayesian approach based on Gaussian Processes and involving archeomagnetic data as a reference.While leveraging all archeomagnetic data components (declination, inclination, and intensities) we exclusively utilize the directional components from sediment data.This choice comes from the recognition that the lock-in function governing directional data may substantially differ from that governing intensities.A discussion of this point can be found in Section 4. Importantly, our method has the flexibility to be applied equally to intensities or even to all three components simultaneously.Also, largely unknown intensities in paleomagnetic records are challenging and motivated us to focus on directional components (e.g., Roberts et al., 2013).

Data Model
In this section, we will derive the data model which describes the relation between the measured signal in the sedimentary records and the geomagnetic field variations.While our primary focus here is on the sedimentary records, we need the information from archeomagnetic records at a later stage.The model of archeomagnetic data is outlined in Schanner et al. (2022).
The first functional   , used to describe the data model, is associated with the offset and smoothing caused by the lock-in process and given by . Later we will apply this function to the geomagnetic field which we model as such a function.The linearity of the functional   follows directly from the linearity of the integral.
Besides the natural smoothing caused by the lock-in process, there is a smoothing effect caused by the way the magnetization in a sediment core sample is measured.When investigating sediment core samples, cubes or u-channels of different sizes are taken from the core.Afterward, the magnetization in the extracted cube (or points in the u-channel) is measured.The resulting measurement is then an average of the actual magnetization across the width of the cube (or determined by the response function of the magnetometer in the case of u-channels).
We assume that the size of the extracted cubes (or the response function) does not change within a core sample.Therefore, we can define the size of the extracted cubes for one core sample as This results in a second smoothing and can be described by the following measurement smoothing functional The linearity of the functional   again follows from the linearity of the integral.The quantities that are measured in laboratory experiments are not provided in Cartesian field vector components (North [N], East [E], Down [Z]) but as two angles, declination (D) and inclination (I), and intensity (F).The nonlinear relationships between these components can be described by three observation functionals where the three components of  () are given by In the following, we will apply these functionals to the Gaussian Process associated with the geomagnetic field.Note that the lock-in function is defined in depth.Therefore, we can not directly use the time-dependent Gaussian Process given in Equation 1.By switching from time to depth we end up with a new Gaussian Process where the mean function coincides with the mean function of the Gaussian Process given in Equation 1.This is because the mean function is assumed to be constant.The kernel function follows directly by applying the age-depth model (a function that maps time to depth) to the kernel function of the Gaussian Process given in Equation 1.
By applying the functional   to the Gaussian process  B , we get, for all  ∈ ℝ , the first part of our data model Since Note that  2 is also a Gaussian Process.Since both functionals lead to a smoothing we absorb the smoothing associated with the measurement functional into the lock-in function and approximate The nonlinearity results in a data model that is not Gaussian anymore.However, as described in Schanner et al. (2022), these functionals can be linearized by a first-order Taylor expansion.Following Hellio et al. (2014) the linearization results in three functionals where By using these linearized functionals we approximate the data model  3 by a Gaussian Process.In conclusion, the components of our final data model are given by where   = ( E E E )⊤ ∈ ℝ 3 indicates the vector of measurement errors.

Parameter Estimation
In order to estimate the four lock-in function parameters  1-4 , we perform a maximum likelihood estimation (type-II MLE, Rasmussen, 2004).For a Gaussian process, the marginal likelihood is available in closed form (see e.g., Chapter 2.2, Equation 2.30 in Rasmussen, 2004).However, due to the amount of archeomagnetic data used as a reference in our approach, calculating the marginal likelihood is numerically costly (a single function call may take minutes) and hampers numerical optimization.The number of required function calls for global function optimization with four parameters is in the order of thousands.This makes estimation via the closed form unfeasible. Instead, we perform a sequentialization of the marginal likelihood evaluation, similar to Baerenzung et al. (2020) and Schanner et al. (2022).
The idea is to replace the closed form Gaussian process regression by a Kalman filter (Kalman, 1960).The closed form marginal likelihood is approximated by a sum over the marginal likelihood values calculated for the individual Kalman filter steps (see Equation 24in Baerenzung et al., 2020).The resulting expression provides a measure of how well a set of lock-in function parameters describes the pDRM process in a single sediment core, given the global set of archeomagnetic data and the respective sediment data.In other words, we use the archeomagnetic data to estimate the shape of the lock-in function for a single core.Due to the temporal distribution of the archeomagnetic data, we limit the estimation to the last eight thousand years.
A crucial difference between the existing implementations in Baerenzung et al. (2020) and Schanner et al. (2022) is the modified observation function discussed in Section 2.3.The convolution integral can be interpreted as a delay (and smoothing) in the measurements, resulting in cross-correlations between the Kalman filter steps.To respect these cross-correlations, the Kalman filter state vector is modified to contain recent time steps.The number of time steps   has to be big enough, so that the corresponding time period covers the assumed maximal lock-in depth   (Choi et al., 2009).The augmentation of the state vector leads, for each  ∈ [0,  ] , to a modified forward operator given as where   = (max, Δ) is the forward operator defined in Baerenzung et al. (2020) The recursive equations for the update step are given by To formulate the backward recursion equations assume that the recursion starts from the last time step   .We set     =   and The backward recursion equations are given as A detailed derivation of these formulas can be found in Appendix A.
Expanding the Kalman filter state vector by previous time steps leads to an increased memory demand and higher computational cost.Choosing a spherical harmonics cutoff degree of  max = 8 and a Kalman filter step size of  Δ = 40 yrs gives a sweet-spot in the trade-off between estimation accuracy and computation time.Several tests showed that lower time steps and higher cutoff degrees have no significant influence on the estimation, while increasing the computational time dramatically.

Log-Marginal Likelihood Optimization
The optimization of the log-marginal likelihood was conducted using dlib's LIPO-TR function optimization algorithm (D.E. King, 2009;Malherbe & Vayatis, 2017).From a mathematical point of view, this algorithm does not guarantee convergence.To prevent the algorithm from running indefinitely, two approaches have been employed.The first method is to limit the maximum number of function calls and using the best result.We started using this method with a maximum number of function calls of 3,500.However, after finding an optimum several times, the optimization algorithm switches to a random search.This random search continues until the maximum number of function calls is exhausted.In the cases, we investigated the random search found the optima after five hundred to one thousand steps and the following random search did not yield any new optima.Therefore, we decided to use the second method where we consider the algorithm as converged when the obtained optimum remains unchanged for several iterations.
The second method significantly reduces the computation time enabling us to perform 50 optimizations for each of the synthetic test cases.We consider two optima  o1, o2 to be the same when where  = 10 −7 denotes the relative tolerance and  = 10 −14 denotes the absolute tolerance.We set the upper bounds for the parameter estimation to 100 cm.
Note that the optimization of such a problem is not straightforward.Other optimization algorithms might lead to better results or better performance.In the future, we will optimize the procedure by including derivatives.

Synthetic Data
All synthetic data points are based on the same reference geomagnetic field time series drawn from the prior described in Section 2.1.Three synthetic data sets were generated from this reference time series.The first data set represents the archeomagnetic data with input locations and times being the same as in the archeomagnetic data used in Schanner et al. (2022).In addition, two synthetic sediment data sets were generated.The first synthetic data set corresponds to a core sample located in Sweden (60°9′3.6″N, 13°3′18″ E) and the other one, to one located on Rapa Iti (27°36′57.6″S, 144°16′58.8″W).Both have the same temporal distribution shown in the lower panel on the left side of Figure 2. The age-depth model used for both synthetic sediment data sets was derived from the age-depth model of the lake sediment core KLK used in Nilsson et al. (2022).
We then applied six different lock-in functions (orange functions in the first column of Figures 4 and 5) to each of the two sediment data sets.The synthetic data from Sweden distorted with the orange lock-in function in (A) of Figure 4 is shown in Figure 3.The reference process is shown in green.Declination and inclination of the synthetic data with measurement uncertainties are shown as blue dots with error bars.Note that we added noise to the data after the distortion associated with the lock-in process.We see an obvious offset as well as a smoothing of the sediment data compared to the reference process.Visualizations of all other synthetic sediment data used in this study as well as notebooks to generate additional synthetic data can be found on our website at https://sec23.git-pages.gfz-potsdam.de/korte/pdrm/.In conclusion, we ended up with six synthetic tests in Sweden and six in Rapa Iti.The reason why we chose these two locations is that Sweden is an area with many paleomagnetic sediment records and a decent coverage of archeomagnetic data while the coverage of archeomagnetic data around Rapa Iti is sparse.

Results
In this section, we will assess the proposed method by conducting synthetic tests.The data utilized in this section, along with the method's implementation, can be found on our website under https://sec23.git-pages.gfz-potsdam.de/korte/pdrm/ and in the corresponding GitLab repository (Bohsung & Schanner, 2023).Moreover, we provide scripts for generating synthetic data, enabling further testing.
Initially, we compared the estimated parameters  1, . . ., 4 with their counterparts in the true lock-in function.The results are visualized in Figure 4 for the synthetic test located in Sweden and in Figure 5 for the synthetic test located in Rapa Iti.The rows (A-F) correspond to the six synthetic test cases.The lock-in function used for the distortion (orange) and the lock-in functions of the 50 optimization runs are plotted in the first column.The colors of the 50 estimated lock-in function correspond to the associated log-marginal-likelihood (log-ml) value.The color ranges from blue for lock-in functions with a high log-ml value (better estimations) to red for lock-in functions with a low log-ml value (worse estimations).The distributions of the log-ml values are shown in the last column.
The visualization of the estimated lock-in functions reveals that the parameters  1, . . ., 4 exhibit some variability in their determination.This variability could potentially be attributed to the optimizer used.Consequently, we proceeded to explore further aspects of the estimated lock-in functions.First, we focused on the half lock-in depth denoted as  0.5 and given as the depth where half of the lock-in process is completed, that is,  0.5 fulfills the equation The second column of Figures 4 and 5 illustrates the distributions of this parameter across the six cases.Subsequently, the remaining two columns present the distributions of the lock-in function heights and widths given as Note that the parameters are weighted with the associated log-ml values.The results from Sweden (see Figure 4) and Rapa Iti (see Figure 5) show similar variance in the estimated parameters, implying that our method is robust enough for solid parameter estimation even for locations where archeomagnetic data is sparse.
As one can see especially the half lock-in depth,  0.5 , is well determined.Its significance lies in its interpretation as the horizontal offset within the observed data.Furthermore, the height and width parameters correspond to the smoothing in the observed sediment records.While not as precisely determined as the half lock-in depth, these parameters remain reliable indicators.In the first moment, the variance in the estimated lock-in functions might seem problematic (see first row of Figures 4 and 5).However, in Figure 6, we compare the data points generated using the real lock-in function (orange), the estimated lock-in function with the best log-ml value (blue), and the estimated lock-in function with the worst log-ml value (red).Obviously, there is not a big difference Note that the differences between the data points are also influenced by the artificially added noise.We compared the distorted sediment observations since the deconvolution of the distorted sediment observations is not straightforward.In a future study, we will work on this challenge.

Discussion
The presented class of parameterized lock-in functions is capable of modeling the offset as well as the smoothing effects associated with the pDRM process.The parametrization with four parameters delivers high flexibility to approximate a wide variety of lock-in behaviors.Functions of higher degree or with more interpolation points could possibly yield better approximations, but would also increase the number of hyperparameters.Decreasing the number of hyperparameters, for example, to one parameter handling the offset and one for the smoothing would lead to less flexibility.We think the presented class of lock-in functions with four parameters gives a good sweet-spot.
In our study, we focus on the directional components (declination and inclination) of sediment records.As mentioned above, this choice comes from the recognition that the lock-in function for the directional components may substantially differ from that for the intensities.
To illustrate this point, consider a scenario where the intensity of the geomagnetic field is constant, while the direction varies.After completion of the lock-in process, the direction recorded in the sediment layer emerges as a weighted average of geomagnetic field directions.Conversely, the sediment layer's intensity appears to decreased in comparison to the beginning of the lock-in process.This is because in the beginning all particles align with the field, but during the lock-in an increasing proportion of particles deviates from this direction, so that their magnetic moments can partly cancel out.The resulting intensity of the whole layer is the sum of the particle magnetic moments and consequently might be decreased.In such cases, the lock-in function for directional components must be distinguished from that of intensity.In fact, the intensity signal captured in the sediment depends not only on the lock-in process, but also on the directional dynamics of the geomagnetic field.This becomes even more complicated if we consider a realistic scenario, where both the field's direction and intensity The orange data points were generated using the real lock-in function in row (A) from (see Figure 4).For the blue data points, we used the estimated lock-in function with the best log-ml value.And for the red data points the estimated lock-in function with the worst log-ml value.
BOHSUNG ET AL. 10.1029/2023JB027373 14 of 17 change over time.Consequently, considering relative paleointensities from sediments is much more complicated than incorporating directions only.
Additionally to the described effect, assuming one lock-in function with identical parameters throughout the sedimentary record remains a reasonable assumption for directional data (at least for the Holocene).It might however be too strict for intensities.Investigation of these challenges will be done in future research.
The synthetic tests presented in Section 3 emphasize the importance of correctly interpreting the estimated lock-in functions.The four estimated parameters  1-4 can show strong variations.Our investigations reveal certain features of the estimated lock-in function that are reliably determined.Particularly significant among these features is the half lock-in depth,  0.5 , which directly corresponds to the offset induced by the pDRM process.The smoothing effects caused by the pDRM process are associated with the height and width of the estimated lock-in function.Although these two parameters are less precisely determined compared to the half lock-in depth, they still yield valuable insights.Especially the comparison of data points distorted with best, worst, and real lock-in function (see Figure 6) shows the accuracy of our estimation even if the parameters  1 to  4 show high variance.The comparison of outcomes for the geographical regions of Sweden and Rapa Iti demonstrates the robustness of our methodology, even in locations with limited coverage by archeomagnetic data.

Conclusion
We present a new class of parameterized lock-in functions.These functions are capable of modeling the offset and smoothing effects associated with the pDRM process.The four parameters deliver high flexibility to approximate a wide variety of lock-in behaviors.Extensive testing on synthetic data sets has demonstrated that our method is highly effective in estimating the parameters associated with the lock-in process, even in areas where archeomagnetic data is sparse.
Transitioning from theoretical developments and synthetic testing to the practical world of real sediment data presents a crucial next step.There are many challenges to this endeavor, including not only the complex pDRM process, but also ancillary effects such as inclination shallowing.Furthermore, real data often provides declination values in relative terms, necessitating the estimation of the declination offset parameter.
Currently, our methodology focuses on the estimation of lock-function parameters.However, the goal is the deconvolution of sedimentary records using the estimated lock-in functions.This endeavor is poised to yield a dedicated sediment preprocessing software, enhancing the reliability of sediment data for geomagnetic field modeling.
Another approach we will work on is the expanding our modeling technique to accommodate multiple sediment cores and potentially incorporating field parameters like correlation length as additional hyperparameters.The advantage of this method is that we can potentially learn from the cross-correlations between different sediment cores.
A comparative analysis of these two methods, simultaneous estimation of the parameters versus parameter estimation as a part of the preprocessing procedure will be the final step.

Appendix A
In this section, we present the derivations of the formulas used in Section 2.4.

𝚺𝚺
15 of 17 The backward recursion equations are derived as follows where the inverse of the matrix   − +1 is derived as follows.We define The inverse of   − +1 is then given by where (⋆) follows since and (⋆⋆) follows since

Figure 1 .
Figure 1.Declination and inclination of a sediment record from Sweden together with measurement errors (blue dots with error bars) are compared to predictions of ArchKalmag14k.r (relying exclusively on archeomagnetic data).An offset and smoothing in the signal of the sediment data compared to the signal derived from ArchKalmag14k.r can be observed.

Figure 2 .
Figure 2. Spatial and temporal distribution of synthetic data.The upper panel on the left side shows the temporal distribution of the synthetic archeomagnetic data while the temporal distribution for the synthetic sediment data is shown in the lower panel.Spatial distributions for synthetic archeomagnetic data (blue dots) and the two synthetic sediment locations (red stars) are shown on the right side.

Figure 3 .
Figure 3.The figure shows the declination and inclination of the noise synthetic sediment data (orange dots) distorted using the orange lock-in function in row (A) of Figure 4. Measurement errors are shown as orange error bars.The reference process for these data is shown in green.

Figure 4 .
Figure 4.The figure shows the results of the 50 parameter estimations for the six synthetic tests located in Sweden.Rows (A-F) correspond to the six cases.The first column shows the true lock-in function (orange) and 50 estimations.The color of the estimated lock-in functions depends on the associated log-ml values ranging from red (low log-ml) to blue (high log-ml).Columns two to three show the distributions of the half lock-in depth, lock-in function height, and lock-in function width.The values of the true parameters are visualized in orange.All parameters are weighted with the log-ml values.The distributions of the log-ml values are visualized in the last column.

Figure 5 .
Figure 5.The figure shows the results of the 50 parameter estimations for the six synthetic tests located on Rapa Iti.See Figure 4 for details.

Figure 6 .
Figure6.The figure shows the comparison of synthetic data generated using three different lock-in functions.The orange data points were generated using the real lock-in function in row (A) from (see Figure4).For the blue data points, we used the estimated lock-in function with the best log-ml value.And for the red data points the estimated lock-in function with the worst log-ml value.
function defined in Section 2.2 and    0 is the lock-in depth, that is, the relative depth where the last particle of the layer at depth is fully locked in.Note that   is a functional that maps functions in   ( ℝ, ℝ 3 ) (set of all continuous functions from  ℝ → ℝ 3 ) to a vector in  ℝ 3 .In the notation above the function   is the argument of the functional   , that is, In conclusion, our lock-in function models not only the offset and smoothing associated with the pDRM process but also the measurement associated smoothing.
2() ≈ 1() .As a heuristic explanation consider two extreme cases.If the lock-in depth   is significantly larger than the size of the sample cube   , the measurement smoothing is negligible.If   is close to zero we can approximate the measurement functional   by the lock-in functional   by choosing the lock-in function with a tableau of width   and setting  =  .

,
max is the spherical harmonics cutoff degree, and  Δ is the step size.Due to the constant step size  Δ , the forward operator does not depend on the Kalman filter step, that is,  ] , the Bayesian filtering equations are recursively defined as = −1 +   =  +  where   is the measurement,   ∼  (, ) the process noise and   ∼  (, ) the measurement noise.The matrix   is the operator that projects the model to the data.The matrix  Σ = 0 − 0 = For  ∈ [1,  ] and with the modified forward operator, the recursive equations of the prediction step are given by