Overcoming multimodalities in age&ndash;depth model posteriors using On&ndash;the&ndash;fly Probability Enhanced Sampling and Parallel Tempering

Kjellson, Hanna; Nilsson, Andreas; Irbäck, Anders

doi:10.5194/egusphere-2026-3766

Preprints

https://doi.org/10.5194/egusphere-2026-3766

Preprints

15 Jul 2026

| 15 Jul 2026

Status: this preprint is open for discussion and under review for Geochronology (GChron).

Overcoming multimodalities in age–depth model posteriors using On–the–fly Probability Enhanced Sampling and Parallel Tempering

Hanna Kjellson, Andreas Nilsson, and Anders Irbäck

Abstract. Bayesian age-depth models are central to paleoclimate research, linking depth in natural archives to calendar age. When synchronizing variable data from different sites, inference in these models requires sampling from high-dimensional, often multimodal posteriors. Consequently, standard samplers become trapped in local optima, miss plausible chronologies, and bias downstream analyses. Here, we investigate two enhanced sampling approaches, On-the-fly Probability Enhanced Sampling (OPES) and Parallel Tempering (PT), and find that tempering only the likelihood components that drive multimodality can improve exploration. Our results suggest that OPES and PT combined with targeted tempering provide robust age-depth models with implications for a wide range of applications in Earth sciences, including paleoclimate and paleomagnetic reconstructions.

Received: 30 Jun 2026 – Discussion started: 15 Jul 2026

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Download & links

Hanna Kjellson, Andreas Nilsson, and Anders Irbäck

Status: open (until 03 Oct 2026)

Post a comment Subscribe to comment alert

RC1:
'Comment on egusphere-2026-3766', Anonymous Referee #1, 22 Aug 2026 reply

I found the manuscript interesting in that it considers alternative sampling approaches for complex Bayesian age–depth models. The idea of exploring alternative sampling strategies for complex age-depth models is certainly interesting.
However, I have substantial concerns with the motivation, experimental design, and statistical validation of the proposed approach. In its current form, I do not think the manuscript establishes what problem is being solved, nor does it demonstrate that the proposed methods improve the computational viability of solving this problems, Bayesian inference or uncertainty quantification relative to existing alternatives. In its current form the manuscript appears to lack of have a concrete objective.
Major comments:

1. The objective of the manuscript.
I found it difficult to determine what the principal objective of the manuscript is. The authors argue that OPES and PT, particularly with partial tempering, perform better than HMC for the examples considered. However, it is not clear what “better” means in this context.
Is the objective to improve computational speed? Effective sample size? Exploration of the posterior? Recovery of the correct posterior distribution? Uncertainty quantification? Robustness to initial conditions? Ease of implementation?
These are different objectives and require different comparisons and validation criteria. Once one of these objective is being well define, and a comparison structure has being decided on, the authors has to establish a comparison with current age-depth methods.
In its current form, the manuscript moves between computational arguments and claims concerning more robust age–depth modelling without establishing a clear quantity against which the methods should be evaluated. Before a comparison between sampling methods can be meaningful, the authors need to define what constitutes successful inference.
2. Multimodality is not itself a problem to be “overcome”
I have a conceptual concern with how multimodality is framed throughout the manuscript, starting with the title. In Bayesian inference, multimodality is not necessarily an undesirable property of the posterior. In age–depth modelling, several substantially different chronologies may be supported by the available information, and such multimodality can therefore represent an important component of the chronological uncertainty that the Bayesian analysis is intended to quantify.
The computational difficulty arises when a sampling algorithm cannot adequately explore this posterior, for example because it becomes trapped within one region of posterior mass or moves only very rarely between separated regions. This distinction is well recognised in the MCMC literature. Medina-Aguayo and Christen (2022), for example, explicitly note that multimodality can arise naturally in well-defined Bayesian problems, including age–depth models, while the computational problem occurs when the sampler becomes trapped in only a subset of the modes.
A successful sampling method should therefore not aim to “overcome” multimodality itself, but rather to overcome the sampling barriers that prevent an adequate representation of a multimodal posterior. This includes not only identifying the relevant regions of posterior support, but also representing their relative posterior probability correctly.
I therefore think that the terminology and framing of the manuscript should be reconsidered. A more precise description of the objective would be to improve exploration of, or overcome sampling barriers within, multimodal posterior distributions rather than to overcome multimodality itself.
This distinction is particularly important for uncertainty quantification. The posterior need not have a simple form: it may contain several modes, asymmetric structures, narrow and broad regions, ridges, or other complex geometry. If these features are supported by the model and the data, they form part of the inferential uncertainty and should be reproduced by the sampling procedure rather than treated as a pathology to be removed.
3. The manuscript does not demonstrate that the proposed algorithms recover the posterior distribution correctly
This follows directly from the previous point.
Demonstrating that a sampler can reach several high-probability regions is not sufficient to demonstrate correct posterior recovery. Medina-Aguayo and Christen (2022) explicitly note that samples obtained from separate modes cannot simply be combined without accounting for the probability mass associated with each region. This issue is also relevant for tempering methods: Tawn et al. (2020) show that standard power tempering can substantially alter modal masses at intermediate temperatures, potentially impairing mixing between modes. Although PT targets the original posterior when $T=1$, adequate finite-sample exploration and recovery of the relative posterior mass of the different regions still need to be demonstrated.
This is particularly important here because the, perceived, motivation of the manuscript is that standard samplers may fail to represent a complicated posterior. However, the main comparison in Figure 4 is based on the mean and 95% quantiles of the potential energy. Similar potential-energy distributions do not demonstrate that two algorithms are sampling the same posterior, since they may visit regions of similar density while allocating very different probability across parameter space.
I therefore do not think Figure 4 provides sufficient validation of the central claim. The authors have a simulation setting in which the underlying age–depth relationship is known, which could be used for stronger validation through repeated simulations assessing coverage and bias.
At present, the results demonstrate improved exploration relative to the HMC implementation considered, but not correct posterior recovery.
4. The comparison with HMC is insufficient to support the broader claims concerning existing sampling approaches
I do not think the comparison with HMC is entirely fair in the context of the manuscript's motivation. The authors introduce BChron, OxCal, and Bacon as established age–depth models, all of which rely on MCMC-based posterior computation in their implementations. Bacon, in particular, uses the self-adjusting t-walk sampler (Christen & Fox, 2010). A more natural benchmark would therefore be against sampling approaches that are already used successfully in age–depth modelling, rather than only against a fixed-step HMC implementation that the authors themselves acknowledge could be improved with NUTS and adaptive step sizes.
This is also relevant to the suggestion that more complex chronological models necessarily require more elaborate sampling schemes. Aquino-López et al. (2018), for example, implemented the PLUM model using the t-walk and reported satisfactory MCMC performance without further tuning, despite jointly estimating several components of the chronological model. Thus, complex age–depth inference has previously been tackled successfully with a general-purpose MCMC sampler.
This does not imply that the t-walk will necessarily solve the particular posterior geometry considered here. However, if the purpose of the manuscript is to establish an advantage of OPES or PT for complex age–depth inference, I believe a fair comparison should include established MCMC approaches already used for these models.
5. The age–depth model itself makes the comparison difficult to interpret
The model is described as being based on Bacon, but an important component of the Bacon model is deliberately removed. In particular, the autoregressive memory structure between accumulation rates is omitted because, according to the authors, it can make the posterior easier to sample.
I find this choice problematic in the context of the main objective of the paper.
If the purpose is to demonstrate that realistic age–depth models generate posterior distributions that cannot be adequately sampled using existing methods, deliberately removing part of an established model because it simplifies sampling may artificially increase the difficulty of the benchmark.
The resulting experiment therefore does not show whether the proposed methods are needed for Bacon itself, or even for a realistic Bacon-type formulation. It shows that they are useful for the modified model constructed here.
I believe that a much more informative analysis would first reproduce the established model and then compare sampling algorithms while keeping the statistical model fixed. Otherwise, the model and sampling algorithm are being changed simultaneously.
At minimum, the authors should investigate whether their conclusions remain when the autoregressive accumulation structure is included.
6. The construction of the accumulation-rate prior needs substantial clarification and justification
Around Eq. (2), the authors state that $\mu$ and $\sigma$ are chosen to obtain an approximately correct mean accumulation rate based on the upper- and lowermost dates. I find this prior construction problematic.
If these same dates are subsequently included in the likelihood, then the observations are being used both to determine the prior and to update that prior through the likelihood. This is effectively an empirical-Bayes or data-informed prior construction, and the uncertainty associated with estimating $\mu$ and $\sigma$ does not appear to be propagated.
More importantly, I do not think such a construction should be the default choice here. Accumulation-rate priors should, where possible, be based on information external to the observations being analysed, for example previous studies of comparable archives, accumulation rates from the same study region, or other independent geological knowledge. If genuinely weak prior information is intended, this should instead be specified explicitly and its sensitivity investigated.
This issue is particularly important in this manuscript because the prior helps determine the shape and multimodality of the posterior that the sampling algorithms are being asked to explore. Using the same data to construct both the prior and likelihood may therefore influence the very computational difficulty that the paper seeks to study.
7. The computational comparison should account for computational cost
If one of the advantages claimed for these methods is computational performance, comparisons based simply on the number of iterations are not sufficient.
PT uses 20--30 temperature levels, each involving HMC calculations, whereas the reference HMC algorithm has a very different computational cost per iteration.
A meaningful computational comparison should therefore use a common measure such as run time, number of likelihood or gradient evaluations, or effective sample size per unit computational effort.
Without such a comparison it is difficult to determine whether the proposed approaches are actually more computationally efficient or simply use substantially more computation.
Conclusion
Taken together, these concerns go beyond issues that could be resolved through revision of the current manuscript. In my view, the paper requires a fundamental reconsideration of its objective, the design of the methodological comparison, the choice of appropriate benchmarks, and the criteria used to evaluate posterior recovery and computational performance.
The topic is interesting, but the present study does not yet provide a sufficiently clear or rigorous basis for its main claims. Addressing these points would require a substantial rethinking and restructuring of the work rather than a conventional revision.

Aquino-López, M. A., Blaauw, M., Christen, J. A., & Sanderson, N. K. (2018). Bayesian analysis of $^{210}$Pb dating. Journal of Agricultural, Biological and Environmental Statistics, 23, 317–333. https://doi.org/10.1007/s13253-018-0328-7
Christen, J. A., & Fox, C. (2010). A general purpose sampling algorithm for continuous distributions (the t-walk). Bayesian Analysis, 5(2), 263–282. https://doi.org/10.1214/10-BA603
Medina-Aguayo, F. J., & Christen, J. A. (2022). Penalised t-walk MCMC. Journal of Statistical Planning and Inference, 221, 230–247. https://doi.org/10.1016/j.jspi.2022.04.008

Tawn, N. G., Roberts, G. O., & Rosenthal, J. S. (2020). Weight-preserving simulated tempering. Statistics and Computing, 30, 27–41. https://doi.org/10.1007/s11222-019-09863-3.

Reply

Citation: https://doi.org/10.5194/egusphere-2026-3766-RC1
- AC1: 'Reply on RC1', Hanna Kjellson, 27 Aug 2026 reply
  
  We would like to thank the reviewer for the constructive comments. We have taken several aspects of the feedback on board and we will revise the manuscript accordingly.
  Fundamentally, we believe that the main concerns raised by the reviewer are due to a misunderstanding of the work presented. A clear example of this is the second comment, where the reviewer stresses the importance of determining the relative weights of different modes, which is exactly what we have successfully done! Therefore, we believe that these concerns can be easily resolved by clarifying the objectives and revising the presentation of our results.
  Below are detailed answers to each comment by the reviewer.
  
  1. The objective of the manuscript.
  The objective of the manuscript is to develop sampling methods able to deal with multimodal posterior distributions encountered in age-depth models based on measurement data (e.g. d18O) that are synchronized to a reference timeseries. This problem represents a frequently encountered and largely unsolved challenge in Earth Sciences.
  The goal is for the methods to correctly recover multimodal posterior distributions. We validate our conclusions by comparing results from extensive calculations with two different methods, On-the-fly Probability Enhanced Sampling (OPES) and Parallel Tempering (PT). We find that OPES and PT, with targeted tempering, identify the same modes and assign nearly identical relative weights to the modes. In both cases, our starting point is Hamiltonian Monte Carlo (HMC), to which we add tempering techniques, either OPES or PT. For comparison, we also conduct simulations with HMC only. Not unexpectedly, pure HMC tends to get stuck in one mode and fails to identify all major modes.
  We note that the objective could be stated more clearly. We will revise the introduction to make the objective clear.
  
  2. Multimodality is not itself a problem to be “overcome”
  The title could indeed be rephrased. A suggestion for a new and more representative title is "Overcoming barriers in multimodal age-depth model posteriors using targeted tempering methods".
  With a rephrased title and clarified objective in the introduction, we believe that the general goal of finding all modes and recovering the correct posterior distribution will be made clear. For instance, we show several modes in Fig. 2 and Fig. A1 that are recovered with the same weights for both PT and OPES. The very close agreement between the PT and OPES results strongly indicates that the main modes and their relative weights can be reliably recovered. Again, we emphasize that these results have been confirmed through multiple independent runs for each method.
  
  3. The manuscript does not demonstrate that the proposed algorithms recover the posterior distribution correctly
  We believe that this comment is due to a misunderstanding of what we want to convey with Fig. 4. This figure was never intended to be viewed in isolation. The main purpose of Fig. 4 is to show that the sampling methods without targeted tempering fail. The recovery of the correct modes when targeted tempering is applied is shown by the fact that both PT and OPES yield identical results. We thank the reviewer for making us aware of this issue. In the revised manuscript we will extend Fig. 4 to also illustrate the reproducibility of the results with OPES and PT, as already discussed in 1 and 2.
  
  4. The comparison with HMC is insufficient to support the broader claims concerning existing sampling approaches
  We first want to clarify that the age-depth modelling that we do differs from what is done in Bacon, OxCal, and BChron because of the additional synchronization. This synchronization introduces serious multimodality, see for example Fig. S3 in Nilsson et al. (2022), with modes that are much narrower and more separated than those seen in, e.g., Fig. 7 in Medina-Aguayo and Christen (2022).
  Since this type of Bayesian synchronization approach is relatively new, with only a few studies attempting to implement similar methods (e.g. Nilsson and Suttie, 2021, Lee et al., 2023, Edmonsond and Dyer. 2024 and Eichenseeer et.al., 2025), we argue that there is no clear benchmark sampler. We choose to use HMC as it is a well-established general-purpose sampler for high-dimensional targets, successfully used for synchronization of slower-varying data in Nilsson et al. (2022). To overcome barriers in multimodal posterior landscapes we add tempering techniques, PT and OPES, while HMC remains the basic simulation engine. Therefore, to evaluate the usefulness of PT and OPES, a natural choice is to compare to a pure HMC sampler.
  The reviewer suggests that "complex age–depth inference has previously been tackled successfully with a general-purpose MCMC sampler”. However, we would like to point out that the posteriors in the provided reference (Aquino-López et al., 2018) show no clear multimodality (as far as we can tell) and, from a pure sampling stand-point, are not as challenging as the posteriors in our systems.
  Regarding the suggestion that a fair comparison should include established MCMC approaches used for similar models, such as the t-walk used by Aquino-López et al. (2018), we would like to re-iterate that the models referred to are not similar to ours. We note, however, that the t-walk sampler has previously been implemented for age-depth models including synchronization (Nilsson et al., 2018), but due to slow convergence it was abandoned in favour of HMC in subsequent studies (Nilsson and Suttie, 2021; Nilsson et al., 2022).
  
  5. The age–depth model itself makes the comparison difficult to interpret
  First, it is important to emphasize that, as clearly stated in the current manuscript, the goal of this study is not to create the best age-depth model, but to develop sampling methods to overcome barriers in multimodal posterior landscapes using realistic examples. We argue that, despite removing the memory term, our age-depth models are realistic (see below for more details). More importantly, since we are dealing exclusively with synthetic data and the "true" age-depth model is generated from the prior, any question with regards to the appropriateness of the prior becomes irrelevant.
  As discussed in Trachsel and Telford (2016), there are currently no age-depth models that work perfectly, and we thus choose a simple one that allows for a broad range of solutions, in particular excluding the memory term because it can sometimes artificially exclude possible solutions (Trachsel and Telford, 2016). This is particularly true for stalagmites, where large variations in growth are fairly common and the inclusion of a memory term may promote unrealistically smooth age-depth models.
  Changing the age-depth model could indeed change the posterior landscape slightly (see for example discussion around Fig. S3 in Nilsson et al., 2022), but if chosen reasonably, it would not change the fact that the synchronized data introduces narrow and separated modes.
  The motivation for the exclusion of the memory term in the age-depth models will be further clarified in the revised manuscript.
  
  6. The construction of the accumulation-rate prior needs substantial clarification and justification
  We would first like to re-iterate that the goal of this study is not to create the best age-depth model, but to develop sampling methods to overcome barriers in multimodal posterior landscapes using realistic examples. Since we are dealing exclusively with synthetic data and the "true" age-depth model is generated from the prior, any question with regards to the appropriateness of the prior becomes irrelevant.
  Nevertheless, with regard to the accumulation-rate prior, our choice does not affect the shape and multimodality of the posterior to the extent that it matters for the sampling. We find that the model is rather insensitive to the prior mean accumulation rate (within reasonable choices), which is true even when only the absolute ages are included. Also, the data are synthetic, and a similar prior could have been obtained in a real scenario, when considering external information.
  
  7. The computational comparison should account for computational cost
  To achieve proper sampling of our multimodal systems with HMC is prohibitively time-consuming (within our resources), whereas we do obtain satisfactory results with targeted PT and OPES. A quantification of the speed-up gained with PT or OPES is therefore not possible.
  For PT (a.k.a. replica exchange), the computational cost is proportional to the number of replicas (or temperatures). It is worth noting, however, that the wall-clock time may increase much more slowly with the number of replicas, since the algorithm is extremely well suited for parallel computations.
  OPES requires the determination of a biasing potential, which is done in an adaptive fashion. For this initial step, we used several replicas. However, unlike PT, OPES is fundamentally a single-replica method. The cost per iteration is essentially the same as for HMC.
  Our main point is not to show that PT and OPES are fast or faster than other samplers, but that if we adapt them slightly (using targeted tempering), we can sample from very difficult posteriors that we have not been able to sample from before.
  
  Conclusion
  To conclude, we believe that the suggestions made by the reviewer regarding clarifying the objective of the study will strengthen this work. However, we strongly disagree with the conclusion that these concerns go beyond what can be resolved through revision.
  
  References
  Aquino-López, M. A., Blaauw, M., Christen, J. A., & Sanderson, N. K.: Bayesian analysis of $^{210}$Pb dating. Journal of Agricultural, Biological and Environmental Statistics, 23, 317–333, https://doi.org/10.1007/s13253-018-0328-7, 2018.
  Edmonsond, S. and Dyer, B.: A Bayesian framework for inferring regional and global change from stratigraphic proxy records (StratMCv1.0), Geoscientific Model Development, 18, 4759–4788, https://doi.org/10.5194/gmd-18-4759-2025, 2025.
  Eichenseer, K., Sinnesael, M., Smith, M. R., and Millard, A. R.: StratoBayes: a Bayesian method for automated stratigraphic correlation and age modelling, Geochronology, 7, 545–570, https://doi.org/10.5194/gchron-7-545-2025, 2025.
  Lee, T., Rand, D., Lisiecki, L. E., Gebbie, G., and Lawrence, C.: Bayesian age models and stacks: combining age inferences from radiocarbon and benthic δ18O stratigraphic alignment, Climate of the Past, 19, 1993–2012, https://doi.org/10.5194/cp-19-1993-2023, 2023.
  Medina-Aguayo, F. J., & Christen, J. A.: Penalised t-walk MCMC. Journal of Statistical Planning and Inference, 221, 230–247, https://doi.org/10.1016/j.jspi.2022.04.008, 2022.
  Nilsson, A. and Suttie, N.: Probabilistic approach to geomagnetic field modelling of data with age uncertainties and post–depositional magnetisations, Physics of the Earth and Planetary Interiors, 317, 106 737, https://doi.org/https://doi.org/10.1016/j.pepi.2021.106737, 2021.
  Nilsson, A., Suttie, N., Hill, M.J.: Short-Term Magnetic Field Variations From the Post-depositional Remanence of Lake Sediments, Frontiers in Earth Science, 6, https://www.frontiersin.org/journals/earth-science/articles/10.3389/feart.2018.00039, 2018.
  Nilsson, A., Suttie, N., Stoner, J., and Muscheler, R.: Recurrent ancient geomagnetic field anomalies shed light on future evolution of the South Atlantic Anomaly, Proceedings of the National Academy of Sciences, 119, e2200749 119, https://doi.org/10.1073/pnas.2200749119, 2022.
  Trachsel, M. and Telford, R.: All age-depth models are wrong, but are getting better, The Holocene, 27, https://doi.org/10.1177/0959683616675939, 2016.
  
  Reply
  
  Citation: https://doi.org/10.5194/egusphere-2026-3766-AC1

Hanna Kjellson, Andreas Nilsson, and Anders Irbäck

Viewed

Total article views: 139 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
99	30	10	139	12	10

HTML: 99
PDF: 30
XML: 10
Total: 139
BibTeX: 12
EndNote: 10

Views and downloads (calculated since 15 Jul 2026)

Month	HTML	PDF	XML	Total
Jul 2026	83	25	7	115
Aug 2026	16	5	3	24

Cumulative views and downloads (calculated since 15 Jul 2026)

Month	HTML	PDF	XML	Total
Jul 2026	83	25	7	115
Aug 2026	16	5	3	24

Viewed (geographical distribution)

Total article views: 83 (including HTML, PDF, and XML) Thereof 83 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 28 Aug 2026

Short summary

Natural records provide valuable information about past climate variability. However, age uncertainties limit our ability to resolve the timing of climate events and identify leads and lags within the climate system. In this project, we successfully combine Bayesian modeling with problem-specific enhanced sampling techniques to reduce these uncertainties on four datasets. Applied to complete climate reconstructions, this could improve our understanding of past and future climate change.


Total:	0
HTML:	0
PDF:	0
XML:	0