CIR Process Age Inference Algorithm v1.0: scalable and consistent sedimentation rate modeling for ocean sediment cores via the Cox-Ingersoll-Ross process

Lee, Taehee; Lisiecki, Lorraine E.; Newall, Samuel R.; Zhou, Yuxin; Lawrence, Charles E.

doi:10.5194/egusphere-2026-2299

Preprints

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

Preprints

24 Jun 2026

| 24 Jun 2026

CIR Process Age Inference Algorithm v1.0: scalable and consistent sedimentation rate modeling for ocean sediment cores via the Cox-Ingersoll-Ross process

Taehee Lee, Lorraine E. Lisiecki, Samuel R. Newall, Yuxin Zhou, and Charles E. Lawrence

Abstract. Deep-sea sediment cores provide a key archive of past climate, with geochemical measurements of microfossils providing both environmental information and age constraints; however, constructing continuous age-depth models is challenging due to sparse and uncertain direct dating (e.g., radiocarbon up to 50 kyr BP) and the need to integrate more densely sampled but indirect proxies such as benthic δ¹⁸O, which requires pattern alignment to reference stacks while maintaining physically plausible sediment rates. Bayesian frameworks have therefore become standard, with approaches such as BACON (Blaauw and Christen, 2011) that models inverse sedimentation rates via an autoregressive gamma process and BIGMACS (Lee et al., 2023) that integrates both direct and indirect constraints through probabilistic alignment and empirically informed priors. Despite their practical utilities, these method exhibit key limitations: BACON relies heavily on user-specified hyperparameters that are not statistically inferred from the data, while BIGMACS could employ a sedimentation-rate model defined on uneven depth grids, potentially leading to inconsistent smoothness and sensitivity to proxy resolution. Here we propose the Cox-Ingersoll-Ross (CIR) Process Age Inference Algorithm to tackle the aforementioned limitations. This multi-layer Bayesian hierarchical framework employs the CIR process as a prior on the inverse sedimentation rates to guarantee a consistent smoothness over depths to address the drawback of BIGMACS, and it allows estimation of the CIR process model parameters via the Expectation-maximization (EM) algorithm. To validate our framework, we first estimated the model parameters from a carefully curated dataset of 79 radiocarbon records, and then applied the algorithm to four other radiocarbon-dated benthic δ¹⁸O sediment core records to compare the performance to BIGMACS. The resulting age models not only show greater consistency and robustness but also preserve smoothness of posterior sedimentation rates over depths and successfully avoid alignment artifacts.

Received: 22 Apr 2026 – Discussion started: 24 Jun 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

Taehee Lee, Lorraine E. Lisiecki, Samuel R. Newall, Yuxin Zhou, and Charles E. Lawrence

Status: final response (author comments only)

RC1:
'Comment on egusphere-2026-2299', Anonymous Referee #1, 08 Sep 2026
In this article the authors present a new age modeling algorithm that uses a multi-layer Bayesian framework to infer sedimentation rates (and then age-depth models) in marine cores using δ18O and 14C data. The key advances over past work is more realistic smoothness and process-based fitting of sedimentation rates which was a limitation in BIGMACS, the previous approach to this work.
The manuscript is well written and justified, and the work is important. I have several suggestions to revision to the manuscript and codebases, but believe the paper should be published following these revisions.
I have two major suggestions for the manuscript itself.
The structure and logic of the paper is sound, and the sets up the need for the CIR Process algorithm by discussing the utility and limitations of BIGMACS, the direct predecessor of this algorithm and Bacon, a more general utility age modeling algorithm mostly used for radiocarbon dated sediment sequences. The background and methodological comparisons are solid, and the authors return to the discussion of Bacon in the conclusion. However, Bacon completely disappears from the results and discussion. I recognize that it serves a slightly different purpose, but if possible, including a 1:1 comparison of the method to Bacon, and discussing the similarities and differences, would be an improvement. If it's not possible to do so, consider removing Bacon from the motivation and conclusions, if it's not similar enough to be useful for any of the examples.

Figures 2-6 carry nearly all of the visual results, and could used some improvement. Structurally, there are few problems, mostly that the axis ticks and labels are too small to comfortably read, and the captions need to be greatly expanded in terms of the detail and description of the figures. I know the content is in the text, but it should also (or instead?) be presented in at least the first caption. In addition, I'd suggest considering an alternative presentation. Rather than replicating the same format for each approach, it might be more effective to plot both results on the same axes, simplifying the density clouds. For panels (C), it's stack, it's currently quite difficult to see the differences, as the overall structure is quite similar. It would probably be better to show one, that has both, and perhaps the second panel could be anomalies relative to the stack for both. A similar structure might be helpful for age depth and log sed rate as well. Ultimately, given that the discussion asks the reader to compare the two approaches, the current structure makes the reader work quite hard for that comparison.

Beyond the manuscript, I also reviewed the Matlab and python code shared as zenodo repositories. Given that a primary goal of the study is to enable others to use the methodology presented here to create age models, the functionality and useability of the code is in many ways as important as the manuscript itself. First, I tested reproducing an age model for MD95-2042 using both the Matlab and python approaches, and it worked well. The median ages were very similar using different random seeds and all of the posterior age paths are strictly monotone in depth, and the credible intervals brackets the medians. This all supports the science presented in the manuscript.
I was able to get the Matlab code to run on without statistics or parallel toolboxes in Matlab 2026a by only replacing a single function (normrnd), and could replace that with a very simple (4-line) function. The authors might consider making a similar change to support users who don't have or don't want to install those toolboxes. After making this change, the δ¹⁸O scale/shift and sedimentation rate all agree across this Matlab implementation, the Matlab results that came in zenodo, and python across all four shared cores. This all suggests that the code and results are robust, which is great.
However, a couple minor issues in the code prevented me from running them out of the box:
Trainer_MATLAB/draw_Samples.m doesn't work as documented. The manual says to use draw_Samples('DATA_02282025_K_4_INT_0.5', ...), where the first argument names a file, but that gives an error. Fixed by setting the name as the third argument.

iters_for_R is documented as a flag but used as a number. main.m's header comment and the manual both describe it as yes/no, but AccMODEL.m treats it as an update period. This makes the documented example silently do the opposite of what the manual says. Passing a number works correctly, and I ended up using Inf.

Both of these are simple fixes that need to be done so that people can replicate your results without fixing code.
A couple other thoughts about the code:
The aligners and trainers build the grids slightly differently. Is that intentional?

It would probably be good to add seeding options for reproducibility.

The python code is much slower than the matlab code. It seems like this is mostly because python is building the N x N0 x num_samples tensors inside every SVGD/HMC call, whereas Matlab builds them once per record outside the loop. Unless there's a good reason not to, I'd make that change, because as is the documented examples in python would take a very long time to run on most computers.

The other major issue is the documentation. There is no readme, license, requirements or environment files. I assume these all come from a github repository (one is mentioned on zenodo), but it is dead link (I suspect it's private). Having this code in a well described public git repository is a necessary improvement to help users find and use the code, as well as report issues and make contributions. I believe that the authors plan to add this to github, as it's referenced in a few places, and I encourage them to do this sooner than later.
There are other more minor issues with the documentation:
Aligner_MATLAB supports a user-supplied calibration curve that's not mentioned in the documentation, and as far as I can tell the python version has no equivalent.

The Python manual says C14_individual/ is used by the visualization, but provides no function to generate them for a new core. This functionality does exist in Matlab, and it would be good if it was in python too.

Both MATLAB manuals say you must copy settings files into the input folder, but it works without doing that. Worth clarifying in the instructions.

The trainer manual reports outputs going to a different path than they actually do.

getFigures runs and produces figures, but doesn't save to file, whereas the python version does write out figures.

Python and matlab ask for different input layouts. Better to standardize to let folks share input files without reformatting.

The python manual describes an input folder that doesn't actually seemed to be used, as the notebook takes explicit paths. The calibration_models paths are relative, so the notebook must be run from the root.

Aligner_MATLAB/Defaults/stacks/ is empty
Citation: https://doi.org/10.5194/egusphere-2026-2299-RC1
- AC1:
  'Reply on RC1', Taehee Lee, 06 Oct 2026
  In this article the authors present a new age modeling algorithm that uses a multi-layer Bayesian framework to infer sedimentation rates (and then age-depth models) in marine cores using δ18O and 14C data. The key advances over past work is more realistic smoothness and process-based fitting of sedimentation rates which was a limitation in BIGMACS, the previous approach to this work.
  The manuscript is well written and justified, and the work is important. I have several suggestions to revision to the manuscript and codebases, but believe the paper should be published following these revisions.
  We thank the reviewer for the positive assessment of the manuscript and for the detailed examination of both software implementations. We particularly appreciate the effort to reproduce the MD95-2042 chronology, examine results across random seeds, and compare the MATLAB and Python outputs. These checks provide valuable independent evidence of reproducibility and consistency across implementations.
  
  The comments also identify practical obstacles that users should not have to resolve themselves. We agree that usable, documented software is an important component of this contribution. Below, we describe the changes we will make to the manuscript, software, and documentation in the revision.
  
  I have two major suggestions for the manuscript itself.
  The structure and logic of the paper is sound, and the sets up the need for the CIR Process algorithm by discussing the utility and limitations of BIGMACS, the direct predecessor of this algorithm and Bacon, a more general utility age modeling algorithm mostly used for radiocarbon dated sediment sequences. The background and methodological comparisons are solid, and the authors return to the discussion of Bacon in the conclusion. However, Bacon completely disappears from the results and discussion. I recognize that it serves a slightly different purpose, but if possible, including a 1:1 comparison of the method to Bacon, and discussing the similarities and differences, would be an improvement. If it's not possible to do so, consider removing Bacon from the motivation and conclusions, if it's not similar enough to be useful for any of the examples.
  
  We agree that the motivation, comparisons, and conclusions should be aligned more closely. The primary objective of this study is to improve the age-inference component of BIGMACS in a setting that combines radiocarbon constraints with probabilistic alignment of benthic δ¹⁸O observations to a reference stack. Our numerical comparisons address that objective.
  
  We will therefore follow the reviewer’s second suggestion. We will revise the Abstract, Introduction, and Conclusion to focus the problem statement and claims of improvement on BIGMACS. We will remove wording that presents the method as resolving limitations of BACON or implies comparative superiority over BACON without a corresponding evaluation.
  
  BACON will remain in the relevant background discussion as an established age-depth modeling approach. If needed, we will also retain a concise explanation of the differences between its accumulation rate formulation and the CIR formulation. These references will provide methodological context rather than motivate unsupported comparative claims.
  
  We will clarify that BIGMACS was selected as the benchmark because it accommodates the same combination of radiocarbon constraints and benthic δ¹⁸O stack alignment used in our applications. A comparison under radiocarbon-only conditions would address a different, worthwhile question, but it is not the objective of the present evaluation.
  
  Figures 2-6 carry nearly all of the visual results, and could used some improvement. Structurally, there are few problems, mostly that the axis ticks and labels are too small to comfortably read, and the captions need to be greatly expanded in terms of the detail and description of the figures. I know the content is in the text, but it should also (or instead?) be presented in at least the first caption. In addition, I'd suggest considering an alternative presentation. Rather than replicating the same format for each approach, it might be more effective to plot both results on the same axes, simplifying the density clouds. For panels (C), it's stack, it's currently quite difficult to see the differences, as the overall structure is quite similar. It would probably be better to show one, that has both, and perhaps the second panel could be anomalies relative to the stack for both. A similar structure might be helpful for age depth and log sed rate as well. Ultimately, given that the discussion asks the reader to compare the two approaches, the current structure makes the reader work quite hard for that comparison.
  
  We agree. The figures should allow readers to assess the comparison without repeatedly consulting the main text or moving between separate panels.
  
  We will increase the sizes of axis labels, tick labels, and legends, and check their readability at the intended publication size. We will also expand the captions to identify the plotted quantities, uncertainty intervals, observational constraints, and reference-stack information. The first comparison caption will explain the common plotting conventions in detail, and subsequent captions will identify the relevant core and any figure-specific features.
  
  We will redesign the main comparisons to place the two methods on shared axes where this improves readability. Posterior medians and 95% credible bands can provide a clearer comparison than overlapping clouds of individual posterior paths. Illustrations of sampled paths can be retained separately where they convey information about trajectory behavior.
  
  For the benthic δ¹⁸O panels, we will evaluate a combined alignment display and a complementary residual display relative to the stack. Any residual panel will explicitly define the residual and explain how age uncertainty is represented. We will select the final arrangement based on whether it makes the differences between methods clearer without obscuring their uncertainty.
  
  We will also add a panel showing the difference between the two methods’ posterior median ages at common depths, alongside the age-depth comparison showing each method’s credible bands.
  
  Beyond the manuscript, I also reviewed the Matlab and python code shared as zenodo repositories. Given that a primary goal of the study is to enable others to use the methodology presented here to create age models, the functionality and useability of the code is in many ways as important as the manuscript itself. First, I tested reproducing an age model for MD95-2042 using both the Matlab and python approaches, and it worked well. The median ages were very similar using different random seeds and all of the posterior age paths are strictly monotone in depth, and the credible intervals brackets the medians. This all supports the science presented in the manuscript.
  I was able to get the Matlab code to run on without statistics or parallel toolboxes in Matlab 2026a by only replacing a single function (normrnd), and could replace that with a very simple (4-line) function. The authors might consider making a similar change to support users who don't have or don't want to install those toolboxes. After making this change, the δ¹⁸O scale/shift and sedimentation rate all agree across this Matlab implementation, the Matlab results that came in zenodo, and python across all four shared cores. This all suggests that the code and results are robust, which is great.
  We appreciate these independent checks, including the agreement of median ages across random seeds, the monotonicity of sampled age-depth paths, and the consistency of results between MATLAB and Python.
  
  We will replace calls to normrnd that can be expressed using MATLAB’s built-in randn function with the corresponding location-and-scale transformation. We will also check the documented example workflows for other toolbox dependencies and revise the requirements accordingly. This will allow us to state the dependencies based on the revised code rather than impose unnecessary toolbox requirements.
  
  However, a couple minor issues in the code prevented me from running them out of the box:
  Trainer_MATLAB/draw_Samples.m doesn't work as documented. The manual says to use draw_Samples('DATA_02282025_K_4_INT_0.5', ...), where the first argument names a file, but that gives an error. Fixed by setting the name as the third argument.
  
  Thank you for reporting the error and identifying this discrepancy. We will correct the documented call so that the filename is supplied in the position required by the function. We will check the function header, manual, and example calls for consistency and run the corrected example as written.
  
  iters_for_R is documented as a flag but used as a number. main.m's header comment and the manual both describe it as yes/no, but AccMODEL.m treats it as an update period. This makes the documented example silently do the opposite of what the manual says. Passing a number works correctly, and I ended up using Inf.
  
  Thank you for identifying this documentation error. iters_for_R specifies the number of iterations between updates of the record-specific standardization parameter. Setting it to inf disables these updates and keeps the parameter fixed at its initially specified value, as the reviewer correctly observed.
  
  The description of this argument as a ‘yes’ or ‘no’ flag in the main.m header comment was incorrect. We will correct the header comment to explain the numerical update interval and the use of inf.
  
  Both of these are simple fixes that need to be done so that people can replicate your results without fixing code.
  A couple other thoughts about the code:
  The aligners and trainers build the grids slightly differently. Is that intentional?
  
  Yes, the difference is intentional and reflects the different purposes of the trainer and aligner.
  
  In the trainer, the grid is constructed to maintain a specific interval length in standardized depth coordinates. This allows us to assess the sensitivity of hyperparameter estimation to grid resolution under controlled standardized spacing, as reported in Section 5.1.
  
  In the aligner, the grid is held fixed in the original depth coordinates throughout inference, as described in line 341-344. Reconstructing the grid as the standardization parameter changes would alter the discretized latent representation and could require additional HMC equilibration. Keeping the grid fixed avoids this computational overhead. Although updates to the standardization parameter change the effective grid spacing in standardized coordinates, these changes remain small when the parameter changes modestly.
  
  The resolution experiments support this choice within the range examined. Section 5.1 shows that, for each candidate a, the estimates of b and r remain stable across sufficiently fine standardized grid spacings. More directly, the sensitivity analysis at the end of Section 5.2 shows that the inferred age models remain consistent across the tested grid resolutions. These results support using a fixed grid in the aligner, provided its effective standardized spacing remains within the range for which stability was demonstrated.
  
  We will explain this intentional distinction and its computational rationale in the revised documentation.
  
  It would probably be good to add seeding options for reproducibility.
  
  We agree and will add a documented random-seed option to both implementations. We will explain how to reproduce a run within a given implementation and software environment.
  
  We will also clarify that using the same numerical seed in MATLAB and Python does not imply identical random draws or identical posterior samples across languages. Cross-implementation comparisons should concern the inferred distributions and summaries rather than sample-by-sample identity.
  
  The python code is much slower than the matlab code. It seems like this is mostly because python is building the N x N0 x num_samples tensors inside every SVGD/HMC call, whereas Matlab builds them once per record outside the loop. Unless there's a good reason not to, I'd make that change, because as is the documented examples in python would take a very long time to run on most computers.
  
  Thank you for highlighting the runtime difference and suggesting a possible optimization. On inspecting the Python file currently under review, we find that the integration weights are represented by an N1 x N0 matrix and N1 x num_samples correction arrays, rather than an N1 x N0 x num_samples tensor. These arrays are constructed once per SVGD or HMC function invocation, outside the internal iteration and leapfrog loops.
  
  There is nevertheless scope to reduce repeated setup work across function calls. Some quantities depend on the standardization parameter and sample count, so they cannot necessarily be computed only once per record and reused unchanged throughout inference. We will examine which quantities can safely be cached and refreshed only when their dependencies change. We will also check repeated construction of interpolation objects for the fixed calibration curves and reference stack.
  
  We have not yet established whether this setup work is the principal cause of the reported MATLAB and Python runtime difference. We will profile representative example runs before attributing the difference to a particular operation, and use those measurements to guide optimizations while checking that the underlying calculations are preserved. We will also check the correspondence between the archived release and the inspected file so that any implementation differences are clear.
  
  The other major issue is the documentation. There is no readme, license, requirements or environment files. I assume these all come from a github repository (one is mentioned on zenodo), but it is dead link (I suspect it's private). Having this code in a well described public git repository is a necessary improvement to help users find and use the code, as well as report issues and make contributions. I believe that the authors plan to add this to github, as it's referenced in a few places, and I encourage them to do this sooner than later.
  Thank you for bringing this to our attention. We had prepared the MATLAB and Python GitHub repositories before publication of the preprint and intended to make them public when the preprint became available. Unfortunately, they inadvertently remained private. Both repositories were made public on September 17^th, 2026. We apologize for this oversight and the resulting difficulty accessing the code.
  
  The GitHub repositories can be found in https://github.com/eilion/CIR_Process_Aligner for Python and https://github.com/eilion/CIR_Process_Aligner_MATLAB for MATLAB.
  
  There are other more minor issues with the documentation:
  Aligner_MATLAB supports a user-supplied calibration curve that's not mentioned in the documentation, and as far as I can tell the python version has no equivalent.
  
  Thank you for identifying this inconsistency. Support for user-supplied calibration curve in Aligner_MATLAB was inherited from BIGMACS and was inadvertently retained, although we had intended to exclude this functionality from the current release. Its absence from the Python implementation and documentation therefore reflects the intended scope.
  
  In the revision, we will remove the legacy user-supplied calibration curve option from the MATLAB implementation and clarify the supported calibration options in the documentation.
  
  The Python manual says C14_individual/ is used by the visualization, but provides no function to generate them for a new core. This functionality does exist in Matlab, and it would be good if it was in python too.
  
  Thank you for identifying this documentation ambiguity. The MATLAB function for individual radiocarbon-age inference was inherited from BIGMACS. This auxiliary functionality is outside the scope of the CIR-based inference method presented in this paper was therefore intentionally omitted from the Python implementation.
  
  The individual radiocarbon-age results are used only to optionally display radiocarbon-only age constraints alongside the inferred chronology. Neither CIR-based age inference nor the Python visualization requires these files; the visualization runs without them.
  
  In fact, the Jupyter Notebook already mentioned it as ‘optional’, and we do not plan to include this additional functionality to the Python code.
  
  Both MATLAB manuals say you must copy settings files into the input folder, but it works without doing that. Worth clarifying in the instructions.
  
  Thank you for pointing this out. When no settings file is supplied in the input folder, the MATLAB code uses the default settings file. Copying a settings file into the input folder is therefore necessary only when users wish to override the defaults.
  
  We will revise both MATLAB manuals to explain this fallback behavior and clarify that copying the settings file is optional when the default settings are appropriate.
  
  The trainer manual reports outputs going to a different path than they actually do.
  
  We will correct the trainer documentation to match the output paths used by the code and verify these paths by running the documented example. The revised instructions will identify the files produced and their locations.
  
  getFigures runs and produces figures, but doesn't save to file, whereas the python version does write out figures.
  
  Thank you for this suggestion. We will modify MATLAB’s getFigures function to save the generated figures automatically and document the output directory and file formats.
  
  Python and matlab ask for different input layouts. Better to standardize to let folks share input files without reformatting.
  
  We appreciate the suggestion and recognize that a common input layout would facilitate sharing datasets between implementations. However, the two implementations were developed with different practical aims. The MATLAB code was released primarily to facilitate reproduction of the paper’s results, while the Python implementation was developed to make the method accessible through a freely available programming language, rather than as a direct translation of the MATLAB code.
  
  Our priority is consistency in the underlying inference method and results, rather than identical input interfaces. We will therefore retain the existing input layouts and clarify the required formats and example inputs in the respective manuals. We will also explain the intended roles of the MATLAB and Python implementations in the manuscript so that users can choose the implementation appropriate for their needs.
  
  The python manual describes an input folder that doesn't actually seemed to be used, as the notebook takes explicit paths. The calibration_models paths are relative, so the notebook must be run from the root.
  
  Thank you for pointing out this ambiguity. The input folder described in the manual is a recommended arrangement for organizing user data, rather than a required location; the notebook accepts explicitly specified input paths. We will revise the manual to make this distinction clear.
  
  We will also state explicitly in the Jupyter notebook instructions that the working directory must be the project root folder so that the relative paths to calibration_models resolve correctly.
  
  Aligner_MATLAB/Defaults/stacks/ is empty.
  
  Thank you for pointing this out. The empty Aligner_MATLAB/Defaults/stacks/ folder is a legacy component inherited from BIGMACS and is not referenced in the current manual. We will remove it from the revised distribution to avoid ambiguity.
  
  Citation: https://doi.org/10.5194/egusphere-2026-2299-AC1
RC2:
'Comment on egusphere-2026-2299', Anonymous Referee #2, 09 Sep 2026

This article presents an advancement of the author(s) previous age-depth modelling algorithm "BIGMACS" (Lee et al 2023). The new work is a significant advance, it uses an alternative to the AR(1) gamma process used by Bacon and BIGMACS - the Cox-Ingersoll-Ross (CIR) process - which allows for consistent smoothness along a modelled record - in contract to BIGMACS, in which smoothness varied with the local density of data points.
Like BIGMACS, the new approach can combine both direct (e.g. 14C) and indirect (e.g. matching to a reference stack) age control data, but in a more consistent way. The manuscript focusses on modelling sediment cores with radiocarbon dates, supplemented by matching to a d18O stack.
I was able to get the Python code running, but it is not really packaged sufficiently well for general use. The Github repo referenced from Zenodo was not reachable. In the Zenodo snapshot there was no licensing info and limited documentation.
I have tested the referenced Python code, and it seems to work well for a set of difficult to fit cores. As implied in the introduction sections, the CIR based algorithm is much less sensitive to user chosen meta parameters than is Bacon, particularly the resolution of the model (thickness of modelled sections in Bacon, number of inducing depth in CIR).
I have personal experience using Bacon, but not BIGMACS. In the manuscript, the performance of the new approach is compared only to the author(s) previous age-modelling algorithm "BIGMACS", perhaps for diplomatic reasons. However, I think the most widely used age-depth modelling software for radiocarbon is still Bacon. If they wish to convince people to adopt their new method it would help to have some direct comparisons.
The tunable parameters of the algorithm have been calibrated against a large set of sediment cores - but this set was restricted to marine sediment cores with relatively high sedimentation rates (> 8 cm / kyr). The four example cores also have relatively high and similar sedimentation rates. So it is unclear how well the default parametrisation would work for low sedimentation rate environments or terrestrial cores.
One further concern I have is that d18O values are supplied to the algorithm with no uncertainties, this might put too much weight on the d18O matching given that d18O values can be very uncertain.
I consider this work suitable for publication in GMD as long as the issues with the documentation of the code are corrected. I would recommend adding direct comparisons to Bacon for the work to have more impact in the field.

Detailed comments on the text.
Figures and tables
I found the figure and table legends too brief.
For figs 2,3,4,5 it is difficult to see the differences between CIR and BIGMACS - could the difference be plotted directly?
Abstract
I like the terminology "direct-" and "indirect-proxies" for age to refer to e.g. radiocarbon and d18O values that will be compared to a reference stack
line 65
"Query sedimentation records may then be standardised and integrated into the framework for age estimation."
I guess this means "new" sedimentation records for which an age-depth model is wanted?
line 147
Although technically it is an hierarchical model, I think using this term might be misleading. The description of the model is a closer to a latent process model: a latent stochastic process, deterministically transformed, then observed via a likelihood. It is not hierarchical in the sense of multilevel regression, with regularisation / shrinkage / partial pooling of the parameters via hierarchical levels.
193
"is a non-informative improper prior on the bias random variable B" - should this be biased random variable B?
271-274
Is the ability to estimate these shifts part of the existing code?
284-285
This was for the training data, but can the model handle variable reservoir ages at the inference stage?
323
typo "inducing"
326-327
"Since the record-specifcc standardisation parameter is defined as the median of the average sedimentation rates, it is not surprising that the estimated values of β are close to the corresponding α."
is this because the mean of a gamma is shape / rate and the mean has been set to ~1 by the rescaling, so alpha must be approximately beta by construction?
328-331
"Although V(0.5) is nearly identical to – albeit marginally greater than – V(2.0) for certain interval lengths, the corresponding estimate of ρ is close to zero. This would imply that sedimentation rates exhibit virtually no autocorrelation with respect to depth, a conclusion which contradicts established scientific understanding."
By this you mean you can exclude the apparent good fit at alpha = 0.5 because the fitted rho is unrealistic - why is a good fit achieved in this parameter space?
358-359
I think you mean the lower part of the core (top of the figure, but the deepest part of the core)
361-367
It is not really possible to see the differences between CIR and BIGMACS from these figures. Can you plot the difference directly?
388-394
It would be helpful to have a summarising figure showing some metric of "fit" as a function of grid resolution.

Citation: https://doi.org/10.5194/egusphere-2026-2299-RC2
- AC2:
  'Reply on RC2', Taehee Lee, 06 Oct 2026
  This article presents an advancement of the author(s) previous age-depth modelling algorithm "BIGMACS" (Lee et al 2023). The new work is a significant advance, it uses an alternative to the AR(1) gamma process used by Bacon and BIGMACS - the Cox-Ingersoll-Ross (CIR) process - which allows for consistent smoothness along a modelled record - in contract to BIGMACS, in which smoothness varied with the local density of data points.
  Like BIGMACS, the new approach can combine both direct (e.g. 14C) and indirect (e.g. matching to a reference stack) age control data, but in a more consistent way. The manuscript focusses on modelling sediment cores with radiocarbon dates, supplemented by matching to a d18O stack.
  We thank the reviewer for the positive assessment of the study and for testing the Python implementation on challenging sediment records. We appreciate the comments on the method’s practical behavior, the applicability of the learned parameters, and the presentation of uncertainty. We also agree that the software needs clearer documentation and packaging for general use.
  
  A central revision will be to align the motivation and conclusions more closely with the scope of our evaluation. The study addresses the age-inference component of BIGMACS, combining radiocarbon constraints with probabilistic alignment of benthic δ¹⁸O observations to a reference stack. We will clarify this focus and remove statements that imply superiority over BACON without a direct comparison.
  
  We respond to the general and detailed comments below.
  
  I was able to get the Python code running, but it is not really packaged sufficiently well for general use. The Github repo referenced from Zenodo was not reachable. In the Zenodo snapshot there was no licensing info and limited documentation.
  I have tested the referenced Python code, and it seems to work well for a set of difficult to fit cores. As implied in the introduction sections, the CIR based algorithm is much less sensitive to user chosen meta parameters than is Bacon, particularly the resolution of the model (thickness of modelled sections in Bacon, number of inducing depth in CIR).
  We thank the reviewer for testing the Python implementation on challenging cores and for sharing these encouraging observations.
  
  The MATLAB and Python GitHub repositories had been prepared before publication of the preprint, and we intended to make them public when the preprint became available. Unfortunately, they inadvertently remained private. Both repositories were made public on September 17^th, 2026. We apologize for this oversight and the resulting difficulty accessing the code.
  
  We appreciate the reviewer’s observations regarding sensitivity to model resolution. In the revised manuscript, we will distinguish the grid resolution stability demonstrated in our experiments from broader comparisons with BACON, which were not directly evaluated in this study.
  
  I have personal experience using Bacon, but not BIGMACS. In the manuscript, the performance of the new approach is compared only to the author(s) previous age-modelling algorithm "BIGMACS", perhaps for diplomatic reasons. However, I think the most widely used age-depth modelling software for radiocarbon is still Bacon. If they wish to convince people to adopt their new method it would help to have some direct comparisons.
  We appreciate the reviewer’s practical observations and agree that BACON is an important reference point for readers interested in age–depth modeling. However, our present evaluation addresses a specific task: improving the BIGMACS framework for combining direct chronological constraints with indirect proxy alignment.
  
  BIGMACS was selected as the benchmark because it uses the same types of observations and reference-stack information as the proposed method. This enables a comparison within the multiproxy setting that motivated the study. BACON supports radiocarbon and other chronological constraints, but its standard formulation does not provide the benthic δ18O–stack alignment likelihood used here. A radiocarbon-only comparison of the two methods would therefore address a related but different question.
  
  We recognize that parts of the current motivation and conclusions extend beyond this evaluation. Rather than adding a separate BACON benchmark, we will revise the Abstract, Introduction, and Conclusion to remove claims or implications of superiority over BACON. We will focus the stated advances on the BIGMACS comparison actually undertaken.
  
  BACON will remain in the background discussion where it provides relevant context for age–depth modeling and accumulation-rate priors. We will distinguish that context from the methods evaluated experimentally. We will also limit statements about sensitivity to model resolution to the grid-resolution results demonstrated in this study, rather than generalizing them to comparative performance against BACON.
  
  The tunable parameters of the algorithm have been calibrated against a large set of sediment cores - but this set was restricted to marine sediment cores with relatively high sedimentation rates (> 8 cm / kyr). The four example cores also have relatively high and similar sedimentation rates. So it is unclear how well the default parametrisation would work for low sedimentation rate environments or terrestrial cores.
  We agree that the applicability of the model formulation must be distinguished from the transferability of the learned default parameters.
  
  The training dataset was deliberately restricted to selected marine records, including the sedimentation-rate threshold described in Section 4. Consequently, the estimated parameters characterize the training collection under the assumptions of our model. The four application cores do not establish that these defaults are appropriate for low-sedimentation-rate or terrestrial environments.
  
  Depth standardization makes it possible to express different records on a common scale, but it does not establish that their sedimentation dynamics, depositional discontinuities, or bioturbation effects are equivalent. We will make this distinction explicit and qualify statements suggesting globally applicable parameter values.
  
  The revised Discussion will state that the default parameterization has been evaluated for the marine applications presented here. Applications outside that setting may require sensitivity analysis, re-estimation using representative training records, and appropriate calibration and observation models. We will present this as a limitation of the current validation rather than claim that standardization alone guarantees transferability.
  
  One further concern I have is that d18O values are supplied to the algorithm with no uncertainties, this might put too much weight on the d18O matching given that d18O values can be very uncertain.
  Thank you for raising this distinction. Although the current input does not supply a separate uncertainty for each benthic δ¹⁸O observation, the alignment likelihood does not treat the observations as exact matches to the stack. Equations (8) and (9) incorporate the stack variance, with Equation (9) providing a Student’s t likelihood for a more robust treatment.
  
  Nevertheless, stack uncertainty and observation-specific uncertainty are not interchangeable. The current formulation does not separately represent every source of analytical uncertainty or record-specific discrepancy, and the use of a robust likelihood does not remove the need to distinguish these sources.
  
  We will clarify precisely what uncertainty is represented by the stack variance and the chosen likelihood, and what is not represented separately. We will also discuss the potential consequences for the relative influence of δ¹⁸O alignment when additional observation or record-specific uncertainties are substantial. We will avoid suggesting that the existing likelihood necessarily captures all uncertainty in individual δ¹⁸O measurements.
  
  I consider this work suitable for publication in GMD as long as the issues with the documentation of the code are corrected. I would recommend adding direct comparisons to Bacon for the work to have more impact in the field.
  We thank the reviewer again for the constructive suggestions and independent software testing. The planned revisions will sharpen the scope of the paper, clarify the assumptions and uncertainty represented by the model, improve the comparative figures, and make the software more accessible for practical use.
  
  Detailed comments on the text.
  Figures and tables
  I found the figure and table legends too brief.
  We agree and will expand the legends so that the figures and tables can be understood without relying extensively on the main text. The revised captions will identify the quantities displayed, the posterior summaries and credible intervals, the observational constraints, and the meaning of colors and symbols. Table captions will explain parameter definitions, units, and the configurations being compared.
  
  For figs 2,3,4,5 it is difficult to see the differences between CIR and BIGMACS - could the difference be plotted directly?
  We agree that the differences should be easier to assess. We will enlarge the axis labels, tick labels, and legends, expand the captions, and display the two methods on shared axes where appropriate, using posterior medians and 95% credible bands to reduce visual clutter.
  
  We will also add a panel showing the difference between the two methods’ posterior median ages at common depths, while retaining the credible bands for each method in the accompanying age–depth comparison. The caption will clearly identify the difference as a comparison of posterior medians.
  
  For the benthic δ¹⁸O panels, we will evaluate a combined alignment display and a complementary residual display relative to the reference stack. We will define any residual explicitly and select a presentation that makes the differences clearer without obscuring uncertainty.
  
  Abstract
  I like the terminology "direct-" and "indirect-proxies" for age to refer to e.g. radiocarbon and d18O values that will be compared to a reference stack
  We thank the reviewer for this positive feedback on the terminology.
  
  line 65
  "Query sedimentation records may then be standardised and integrated into the framework for age estimation."
  I guess this means "new" sedimentation records for which an age-depth model is wanted?
  Yes, that is the intended meaning. “Query sedimentation records” refers to records for which age-depth models are to be inferred using the learned CIR parameters. We will replace the sentence with “Sediment records for which age-depth models are sought can then be standardized and analyzed using the learned CIR parameters.”
  
  line 147
  Although technically it is an hierarchical model, I think using this term might be misleading. The description of the model is a closer to a latent process model: a latent stochastic process, deterministically transformed, then observed via a likelihood. It is not hierarchical in the sense of multilevel regression, with regularisation / shrinkage / partial pooling of the parameters via hierarchical levels.
  We agree that “Bayesian latent-process model” describes our formulation more precisely. Our use of “hierarchical” referred to the successive layers of the model - a prior on inverse sedimentation rates, a deterministic transformation into an age–depth function, and an observation likelihood - rather than a multilevel structure that induces partial pooling of record-specific parameters.
  
  We will revise the terminology throughout the manuscript to avoid this ambiguity while retaining the explanation of these model components. In particular, we will replace the sentence at line 147 with “Our CIR process age inference method is formulated as a Bayesian latent-process model, in which a CIR prior on inverse sedimentation rates is linked to proxy observations through a deterministic transformation into an age-depth function and an observation likelihood.”
  
  193
  "is a non-informative improper prior on the bias random variable B" - should this be biased random variable B?
  “Biased” is not intended here. We used “bias” to mean an additive offset, but agree that this wording is ambiguous. B represents the age at the reference depth from which the age–depth function is constructed.
  
  We will replace the phrase with “an improper uniform prior on the reference-depth age B, subject to the lower bound B ≥ τ₀.” This states both the interpretation of B and the prior specification directly.
  
  271-274
  Is the ability to estimate these shifts part of the existing code?
  The aligner implementations provided on Zenodo and GitHub already include estimation of the record-specific benthic δ¹⁸O shift and scale parameters. Users can choose whether to estimate these parameters during age inference or keep them fixed at specified values.
  
  We will clarify this functionality in the manuscript and identify the corresponding options in the software documentation.
  
  284-285
  This was for the training data, but can the model handle variable reservoir ages at the inference stage?
  The inference procedure uses the reservoir corrections supplied with the radiocarbon observations. These corrections, and their associated uncertainties, can differ between observations; the model does not require a single correction for the entire core. However, the supplied corrections are treated as fixed inputs rather than estimated during age inference. The algorithm therefore does not jointly infer an unknown time-varying reservoir-age history.
  
  The exclusion described in lines 284–285 concerns the training dataset, where we sought to reduce the influence of uncertain temporal reservoir-age variability on estimation of the shared CIR parameters. Those lines will be changed to say, “We also exclude cores from > 40°N in the Atlantic Ocean due to uncertainty in the temporal variability of marine reservoir ages there.” We will also clarify the distinction between accommodating externally specified, observation-specific corrections and inferring reservoir-age variability.
  
  323
  typo "inducing"
  Thank you. We will correct this typographical error.
  
  326-327
  "Since the record-specifcc standardisation parameter is defined as the median of the average sedimentation rates, it is not surprising that the estimated values of β are close to the corresponding α."
  is this because the mean of a gamma is shape / rate and the mean has been set to ~1 by the rescaling, so alpha must be approximately beta by construction?
  Yes, this is the intended interpretation, with a distinction between approximate normalization and an imposed parameter constraint. Under the shape-rate parameterization, the stationary Gamma mean is α/β.
  
  This normalization makes a fitted stationary mean α/β near one, and hence β near α, unsurprising. However, it does not impose α = β: the path-wise average is distinct from the stationary ensemble mean, and β is subsequently estimated from the data rather than fixed to α. We will clarify this distinction in the revised text.
  
  328-331
  "Although V(0.5) is nearly identical to – albeit marginally greater than – V(2.0) for certain interval lengths, the corresponding estimate of ρ is close to zero. This would imply that sedimentation rates exhibit virtually no autocorrelation with respect to depth, a conclusion which contradicts established scientific understanding."
  By this you mean you can exclude the apparent good fit at alpha = 0.5 because the fitted rho is unrealistic - why is a good fit achieved in this parameter space?
  Thank you for raising this point. Our original wording overstated what the selection criterion establishes and did not sufficiently distinguish that criterion from our interpretation of the inferred dependence structure.
  
  V(α) is a moment-consistency criterion, not a measure of fit to the proxy observations. It compares the specified shape parameter with a shape value calculated from the pooled posterior mean and variance of inverse sedimentation rates. It does not directly assess their dependence across depth. Consequently, candidates with substantially different autocorrelation structures can yield similar values of V, and a small V(0.5) should not be interpreted as evidence of equally good observational or predictive fit.
  
  In the reported calculations, V(2.0) is smaller than V(0.5), so α = 2.0 is selected by the stated rule. However, the small difference indicates that this criterion alone does not strongly distinguish the two candidates. Our preference for α = 2.0 was also supported by the dependence structure: the near-zero fitted ρ for α = 0.5 implies a short correlation length in standardized depth, which is less consistent with the persistence we intended the sedimentation-rate prior to represent. This is a modeling consideration, rather than evidence that α = 0.5 has been statistically excluded.
  
  Therefore, we will revise the text to acknowledge the limited discrimination of the moment-based criterion and remove the categorical statement that the alternative “contradicts established scientific understanding.” We will distinguish the numerical selection of α = 2.0 from the additional scientific rationale for preferring it.
  
  To be specific, we will replace the last two sentences (lines 328-331) of Section 5.1 with the following: “Although V(0.5) is only marginally greater than V(2.0) for some grid spacings, the corresponding estimate of ρ is close to zero (about 0.02), implying a short correlation length for inverse sedimentation rates in standardized depth. The criterion V assesses marginal moment consistency rather than dependence across depth and therefore provides limited discrimination between these candidates. Our preference for α = 2.0 is additionally supported by evidence and modeling arguments for persistence in sedimentation rates. For instance, Figure 3 of Huybers and Wunsch (2004) found that a model with autocorrelated accumulation rates reproduced observed age variances across marine cores under both on- and two-age-control-point configurations more adequately than an uncorrelated model. Such persistence is also represented in the autoregressive model of Blaauw and Christen (2011), and Hercman and Pawlak (2012) motivate their use of LOESS regression in MOD-AGE by arguing that closely spaced layers are more likely to have similar deposition rates than more distant layers. These studies support retaining a dependence structure, although they do not determine the appropriate correlation length for our records or independently exclude the α = 0.5 candidate.
  
  358-359
  I think you mean the lower part of the core (top of the figure, but the deepest part of the core)
  Thank you for pointing out this ambiguity. We will replace “the upper portion of the core” with an explicit description of the depth interval shown in Figure 3, which extends to 17.2 m, to avoid confusion between position within the core and position in the figure.
  
  361-367
  It is not really possible to see the differences between CIR and BIGMACS from these figures. Can you plot the difference directly?
  We agree. This will be addressed through the revised comparison figures described above, including a direct display of posterior median age differences at common depths. We will also revise the discussion so that the cited differences can be identified readily in the corresponding panels.
  
  388-394
  It would be helpful to have a summarising figure showing some metric of "fit" as a function of grid resolution.
  We agree that a quantitative summary would make the grid-sensitivity results easier to assess. We will add a figure showing, for each core, the root-mean-square difference between posterior median ages at each grid resolution and those obtained with the finest tested grid (approximately 0.125), evaluated at common depths.
  
  This metric will summarize sensitivity to discretization rather than goodness of fit to a known chronology; the finest grid will serve as a numerical reference, not as ground truth. The existing Appendix figures will retain the comparisons of posterior uncertainty.
  
  Citation: https://doi.org/10.5194/egusphere-2026-2299-AC2
EC1:
'Comment on egusphere-2026-2299', Heather Kim, 15 Sep 2026

Dear Authors,
Please see the following review from the third reviewer and address the comments by the deadline provided. Best regards, Heather Kim
----
The paper discusses an alternative chronology building inference process in paleo-studies. In a similar approach to other available methods, but using a CIS AR Gamma process, the authors use an informal step-wise hybrid inference process to fit the data. The paper is interesting and the examples are illustrative. However, I have major concerns regarding the overall use of the methodology, the authors conceptual understanding of the problem and the paper presentation.
1. Note that this is a journal article, not some lecture notes. I do not see the need to explain 1) Gaussian Processes, 2) mid point rule, 3) HMC and 4) SVGD. These are standard tools that only need to be mentioned and cited. A supplementary material could be used for a more pedagogical explanation.
2. The model in (10) is a piece-wise linear model, adding the mid point for additional smoothness, and the accumulations rates are modeled with a (CIS) Gamma process. This is as in Blaauw and Christen (2011), but changing the gamma process to the CIS, and the additional knot for integration. The CIS process is governed by two parameters, β and ρ, and these need to be set a priori (see l. 198). An interactive scheme is used in which the chronology is sampled from the posterior and then these parameters are updated. This hybrid scheme runs without justification, with no clear objective funtion in sight, nor a posterior distribution to sample from. What are the authors doing? This is surely not Bayesian; but indeed, “our algorithm returns the age model in the form of samples, not their single estimates only” (l. 203), but samples of what?
3. Why not run an MCMC (HMC) on the full posterior, as in Blaauw and Christen (2011)?
4. Was it worth the great added complication of the CIS model? Using the OU −Γ process is quite simpler, as in Blaauw and Christen (2011). Did the authors make a comparison?
5. The authors are worried about prior specification, but point “estimates” of hyperparameters, in some ad hoc procedure, do not assure us that priors are now estimated from data, and therefore are “better”. Input of prior knowledge of accumulation rates is critical here, given the few radiocarbon, and others, calendar markers along cores. That knowledge is simply not in data, and thus an informative, and committed, prior needs to be used to fill in the gaps.

Citation: https://doi.org/10.5194/egusphere-2026-2299-EC1
- AC3:
  'Reply on EC1', Taehee Lee, 06 Oct 2026
  The paper discusses an alternative chronology building inference process in paleo-studies. In a similar approach to other available methods, but using a CIS AR Gamma process, the authors use an informal step-wise hybrid inference process to fit the data. The paper is interesting and the examples are illustrative. However, I have major concerns regarding the overall use of the methodology, the authors conceptual understanding of the problem and the paper presentation.
  We thank the reviewer for recognizing the interest of the study and the illustrative value of the examples. The comments raise questions about the presentation, the distinction between hyperparameter estimation and posterior inference, and the choice of comparison methods. We address these separately below.
  
  A central distinction is that our study separates learning a shared sedimentation-rate prior from inferring the chronology of an application core. We estimate the prior hyperparameters using 79 training cores and subsequently apply the learned prior to four separate cores. In the application stage, HMC targets the posterior conditional on the learned hyperparameters. This is an empirical Bayes approach, rather than a fully Bayesian analysis that assigns hyperpriors to all shared hyperparameters. We will make this distinction more explicit in the manuscript.
  
  Note that this is a journal article, not some lecture notes. I do not see the need to explain 1) Gaussian Processes, 2) mid point rule, 3) HMC and 4) SVGD. These are standard tools that only need to be mentioned and cited. A supplementary material could be used for a more pedagogical explanation.
  
  We agree that background material should be concise and directly relevant to the proposed method. However, merely mentioning and citing these concepts might not necessarily be sufficient for the intended readership. The paper is written for geoscientists, including readers familiar with age-depth modeling who may not have specialist training in stochastic processes or computational Bayesian inference.
  
  Sections 2.2-2.5 occupy approximately one and a half pages and support specific components of the method. The Gaussian-process discussion introduces the representation of CIR as a sum of squared OU Gaussian processes, which underlies our inference implementation. The midpoint rule explains how nonnegative inverse sedimentation rates are converted into monotone age-depth paths. The HMC and SVGD descriptions explain the distinct roles of posterior sampling and initialization.
  
  These explanations are intended to make the subsequent model formulation accessible, rather than to provide a general tutorial. We therefore propose to retain concise background explanations in the main text while strengthening their connections to the proposed method and trimming general discussion that is not needed.
  
  The model in (10) is a piece-wise linear model, adding the mid point for additional smoothness, and the accumulations rates are modeled with a (CIS) Gamma process. This is as in Blaauw and Christen (2011), but changing the gamma process to the CIS, and the additional knot for integration. The CIS process is governed by two parameters, β and ρ, and these need to be set a priori (see l. 198). An interactive scheme is used in which the chronology is sampled from the posterior and then these parameters are updated. This hybrid scheme runs without justification, with no clear objective funtion in sight, nor a posterior distribution to sample from. What are the authors doing? This is surely not Bayesian; but indeed, “our algorithm returns the age model in the form of samples, not their single estimates only” (l. 203), but samples of what?
  
  We understand “CIS” repeatedly appearing in the comments to refer to the Cox-Ingersoll-Ross (CIR) process.
  
  We respectfully disagree with the characterization that the manuscript provides neither an optimization objective nor a posterior sampling target. Both have been clearly specified in lines 210-230 in Section 3.2.2 and in lines 468-501 in Appendix A.
  
  First, Equation (10) provides a piecewise-linear numerical representation of the age-depth function. The midpoint rule converts inverse sedimentation rates into ages; it does not itself introduce additional smoothness or dependence between rates. Those properties are specified by the Cox-Ingersoll-Ross (CIR) prior.
  
  Second, our intended approach separates learning a shared prior from applying that prior to individual records. Rather than assigning hyperpriors to β and ρ and sampling them jointly with all latent variables, we estimate them from 79 curated training cores, conditional on each candidate α. The selected CIR hyperparameters are then used for age inference in four separate application cores. This is an empirical-Bayes strategy, rather than a fully Bayesian treatment of the shared hyperparameters.
  
  The manuscript specifies the joint model in Equation (5), the augmented formulation in Section 3.2.2, and the multiple-record formulation in Appendix Equations (A1) and (A2). Equation (16) and Appendix Equation (A3) define the EM Q-function, the accompanying M-step specifies its maximization, and Appendix Equation (A5) gives its Monte Carlo approximation. Posterior sampling in this training procedure is used to approximate the expectations required for hyperparameter estimation; the parameter updates are optimization steps, not posterior draws of β and ρ. The separate criterion for selecting α is described in Section 3.2.3 and is not part of that EM maximization.
  
  Regarding “samples of what?”, the sentence, “we draw samples from the posterior of the age function by the HMC algorithm, given the model parameters”, immediately preceding the quoted statement specifies that the samples are drawn from the posterior of the age function, conditional on the model parameters. Section 3.2.4 gives the corresponding mathematical description. More explicitly, HMC samples the latent Ornstein-Uhlenbeck (OU) Gaussian process values and initial age variable conditional on the observations and the selected CIR hyperparameters. These samples determine inverse sedimentation rate paths and, through the integration rule, age-depth paths. The returned age samples therefore represent the conditional posterior distribution of the discretized age-depth function. SVGD is used for initialization, rather than as a replacement for the subsequent HMC sampling.
  
  Treating the shared hyperparameters as estimated quantities does not remove the Bayesian character of this conditional inference: the CIR prior and observation likelihood define the posterior being sampled. We agree that this differs from a fully Bayesian analysis that also propagates hyperparameter uncertainty. The reported age-depth credible intervals are conditional on the selected hyperparameters, and we will state this limitation explicitly.
  
  We will revise the manuscript to distinguish more clearly the training-stage optimization, the separate selection of α, and the application-stage posterior sampling in Section 3.2. We will also replace the potentially confusing wording “specified a priori” with wording that states explicitly that the shared CIR hyperparameters are learned from the training records before application-stage age inference.
  
  Why not run an MCMC (HMC) on the full posterior, as in Blaauw and Christen (2011)?
  
  Our objective is not only to introduce an alternative stochastic model for inverse sedimentation rates, but also to learn a shared, empirically informed prior from a curated collection of 79 marine sediment records and apply it to separate records requiring age inference.
  
  Accordingly, we adopt an empirical-Bayes strategy: β and ρ are estimated from the training records conditional on candidate values of α, and α is selected using the criterion described in Section 3.2.3. The selected hyperparameters then define the inverse sedimentation rate prior used for the application cores. HMC samples the latent variables conditional on these hyperparameters and the application-core observations; the corresponding inverse-sedimentation-rate and age-depth paths represent their conditional posterior distributions.
  
  Thus, HMC is already used for posterior inference. The distinction is that the shared CIR hyperparameters are learned in a separate training stage rather than assigned hyperpriors and sampled jointly with the application-core chronology. This separation allows information from the training collection to inform the prior for an individual core with sparse chronological constraints.
  
  A fully Bayesian analysis that propagates uncertainty in the shared hyperparameters would be a valid alternative. However, it is not necessary for defining or sampling the conditional posterior targeted here. Our choice reflects the intended separation between learning a reusable prior and applying it to new records. We will clarify this rationale and explicitly acknowledge that the resulting credible intervals are conditional on the selected hyperparameters and do not incorporate their estimation uncertainty.
  
  Was it worth the great added complication of the CIS model? Using the OU −Γ process is quite simpler, as in Blaauw and Christen (2011). Did the authors make a comparison?
  
  We did not conduct a direct comparison with BACON. Our principal objective is to address limitations of the age-inference component of BIGMACS in a setting that combines radiocarbon constraints with probabilistic alignment of benthic δ¹⁸O observations to a reference stack. BIGMACS is therefore the directly relevant benchmark, allowing comparison using the same proxy observations and reference stack.
  
  BACON accommodates radiocarbon and other chronological constraints, but its standard formulation does not provide the benthic δ¹⁸O stack alignment likelihood used in this study. Comparing our multiproxy results against a radiocarbon only BACON analysis would confound differences between methods with differences in the information supplied. A radiocarbon-only comparison of both methods would address a separate question. Isolating the effect of CIR versus an alternative process prior would instead require a comparison under a common likelihood and data configuration.
  
  The rationale for choosing CIR is discussed in Section 2. CIR and OU-Γ are distinct stochastic processes despite sharing Gamma marginal distributions under the stationary formulations considered. In particular, CIR provides continuous inverse sedimentation rate paths. Our choice reflects this modeling property, while the OU Gaussian-process representation provides a means of implementing inference.
  
  We agree that the choice of CIR must be distinguished from claims of comparative superiority. The results demonstrate the behavior of the proposed method relative to BIGMACS in the applications examined; they do not establish that CIR is universally preferable to OU-Γ in BACON. We will clarify this scope.
  
  The authors are worried about prior specification, but point “estimates” of hyperparameters, in some ad hoc procedure, do not assure us that priors are now estimated from data, and therefore are “better”. Input of prior knowledge of accumulation rates is critical here, given the few radiocarbon, and others, calendar markers along cores. That knowledge is simply not in data, and thus an informative, and committed, prior needs to be used to fill in the gaps.
  
  We agree that informative priors are important when chronological observations are sparse. This is precisely the role of the CIR prior in our model. The issue is therefore how that prior is specified, rather than whether an informative prior is needed.
  
  We respectfully disagree with characterizing the entire estimation procedure as “ad hoc.” Conditional on each candidate α, β and ρ are estimated through the explicitly formulated EM procedure: Equation (16) and Appendix Equation (A3) define the Q-function, the M-step specifies its maximization, and Appendix Equation (A5) gives the Monte Carlo approximation. These estimates use observations from 79 curated training cores, rather than relying solely on the sparse chronological constraints within an individual application core.
  
  The selection of α is a separate, moment-based procedure described in Section 3.2.3. We will distinguish it explicitly from the EM estimation of β and ρ and clarify its limitations, including its limited discrimination between some candidates. We do not claim that all three hyperparameters are selected through a single likelihood-maximization procedure.
  
  The distinction between the training collection and an application core is central to our approach. Although a few dates within one core may provide limited information about sedimentation rate variability, this does not establish that a collection of dated cores contains no information relevant to specifying a shared prior. Our aim is to use that information under an explicit common process assumption, while retaining an informative prior to constrain age inference between observations in a new core.
  
  Nevertheless, being estimated from data does not by itself make a prior superior. Its suitability depends on the information available in the training collection and the extent to which its assumptions apply to the target record. Point estimation also does not propagate hyperparameter uncertainty into the conditional age-depth posterior. We will state these qualifications explicitly in the Conclusion section.
  
  We will therefore revise the motivation in the Abstract and Introduction and discussion in the Conclusion section to remove suggestions that the proposed prior is generally superior to BACON’s prior specification. BACON will remain as relevant methodological background, while comparative claims will be restricted to the BIGMACS experiments actually presented. Specifically, the Conclusion will clarify that our contribution is an empirically informed CIR prior and an associated multiproxy age-inference framework, not a demonstration that learned priors universally outperform expert-specified priors or replace geological judgment.
  
  Citation: https://doi.org/10.5194/egusphere-2026-2299-AC3

Taehee Lee, Lorraine E. Lisiecki, Samuel R. Newall, Yuxin Zhou, and Charles E. Lawrence

Viewed

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

HTML	PDF	XML	Total	BibTeX	EndNote
223	99	31	353	13	13

HTML: 223
PDF: 99
XML: 31
Total: 353
BibTeX: 13
EndNote: 13

Views and downloads (calculated since 24 Jun 2026)

Month	HTML	PDF	XML	Total
Jun 2026	1	1	0	2
Jul 2026	48	33	4	85
Aug 2026	34	11	3	48
Sep 2026	100	40	13	153
Oct 2026	40	14	11	65

Cumulative views and downloads (calculated since 24 Jun 2026)

Month	HTML	PDF	XML	Total
Jun 2026	1	1	0	2
Jul 2026	48	33	4	85
Aug 2026	34	11	3	48
Sep 2026	100	40	13	153
Oct 2026	40	14	11	65

Viewed (geographical distribution)

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

Country	#	Views	%

Latest update: 14 Oct 2026

Short summary

We propose a new method, the CIR Process Age Inference Algorithm, to improve age estimates from deep-sea sediment cores. It uses a Cox-Ingersoll-Ross process as a prior on inverse sedimentation rates, ensuring smooth, realistic changes with depth. The model parameters are learned from 79 sediment records. Tested on real sediment core data, our method produces more consistent, robust age estimates than the software BIGMACS, avoids alignment errors, and maintains smooth sedimentation patterns.


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