practical.qmd

---
title: "Age-depth modeling in R"
format:
  html:
    embed-resources: true
    toc: true
    toc-depth: 3
author: 
  - Niklas Hohmann
  - David De Vleeschouwer
  - Christian Zeeden
bibliography: references.bib
execute: 
  cache: true
---

# Age-depth modeling in R: Practical

These are the materials a workshop on **Age-depth modeling in R**, held at the [CycloNet meeting in Utrecht (June 26th 2025)](https://www.uu.nl/en/research/department-of-earth-sciences/cyclonet-meeting-utrecht) and the [GESEP School 2026](https://dfg-spp-icdp.de/de/forschung/kolloquium/icdp-iodp3-kolloquium-2026/gesep-school-2026) as part of the ICDP-IODP3 Kolloquium in Münster (March 8th to 12th 2026).

### Learning goals

The learning goals for this practical are

-   Getting to know the available software packages for deep-time age-depth modeling in R
-   Understand the syntax and model assumptions of the software packages
-   Determine ages of horizons and time contained within a stratigraphic interval
-   Transform proxy records using these age-depth models

We cover three packages: modifiedBChron [@trayler2019] , astroBayes [@trayler2023] and admtools [@hohmann_nonparametric_2024].

-   `modifiedBChron` use radioisotope data to construct age-depth models, building on the sediment accumulation model of BChron [@parnell2008] originally developed for radiocarbon data

-   `astroBayes` combines radioisotope data with orbital fitting, and assumes piecewise constant sedimentation rates

-   `admtools` provides two methods: ICON and FAM. ICON combines tie point data with estimates on sedimentation rates (e.g. from orbital matching or semiquantitiative data) to construct age-depth models, while FAM matches assumed tracer fluxes with observations on tracer values in the section

### Package installation

The code below assumes you have all three packages installed, plus the `dplyr` package for some data wrangling and `StratPal` for generating synthetic proxy records. To install them, run the following code:

```{r}
#| eval: false
install.packages("remotes")

install.packages("admtools")
install.packages("StratPal")
install.packages("dplyr")

remotes::install_github("robintrayler/modifiedBChron")
remotes::install_github("robintrayler/astroBayes")
```

### Workshop structure

This part of the workshop is a practical, meaning the participants work through the materials on their own pace, and the trainers are available for questions and discussion. Estimated time for the practical is 2 to 3 h. You can work through the materials in any order, so feel free to start with a package that seems most relevant to you. To familiarize yourself with the syntax and handling of the packages, a list of tasks is provided at the end of each section. They are not mandatory, but feel free to work through them to deepen your understanding of the functionality.

20 minutes before workshop end we will come together to discuss insights you had and problems you ran into. For this part, please think about

1.  What model assumptions are justified for your specific time scale and depositional environment?
2.  What functionality of age-depth modeling software is most useful for you? What is missing in the ecosystem?

### Further reading

If you want to dive deeper into the mechanics of age-depth modeling, we recommend working through the materials in @david_de_vleeschouwer_2021_4749027. They provide excellent examples on the problems that can arise in age-depth modeling, and step-by-step instructions to build complex age-depth models.

## modifiedBChron

### Introduction

The package `modifiedBChron` [@trayler2019] is a modification of `BChron` [@parnell2008] for deep-time usage. `BChron` was originally developed for usage with radiocarbon records, and assumes that each possible chrononlogy is piecewise linear, with a random number of points between dated horizons.

### Usage

We begin by loading the package and the data from @trayler2019. This code is modified from [this repository](https://github.com/robintrayler/modifiedBChron).

```{R}
library(modifiedBChron)
data("SCF") # load data from the Santa Cruz formation
```

The age-depth model is estimated using the function `ageModel`. This might take a few minutes to run.

```{r}
#| results: hide
#| warning: false
age_model = modifiedBChron::ageModel(ages = SCF$age,
                                      ageSds = SCF$ageSds,
                                      positions = SCF$position,
                                      ids = SCF$ids,
                                      positionThicknesses = SCF$thickness,
                                      distTypes = SCF$distType,
                                      predictPositions = seq(48.25, 193.50, by = 1),
                                      MC = 10000,  # how many iterations 
                                      burn = 2000) # how many iterations to discard
```

It takes the following data:

-   `ages`: mean age

-   `ageSds`: age uncertainty - standard deviation for normal distributions and half range for uniform distribution

-   `id`: sample ID. Entries with identical ID will be combined into one probability distribution

-   `positions`: stratigraphic position of the sample above base of the section

-   `positionThickness`: stratigraphic uncertainty of positions as half thickness

-   `distType` : `U` or `G` for a uniform distribution or a normal distribution

The age-depth model can be plotted using

```{R}
modifiedBChron::modelPlot(model = age_model,
                          scale = 5)
```

### Model assumptions

First, `modifiedBChron` combines multiple radiometric dates specified by a normal or uniform distribution into tie point uncertainties. The tie points can also have stratigraphic uncertainty (as specified by `positionThickness`). Then, it uses the model by @parnell2008 to estimate the age uncertainty between the tie points. In this model, the chronology between tie points is approximated using a monotonous, piecewise linear function. Here, number of accumulation events is determined by a Poisson distribution and the amount of sediment accumulated is drawn from a gamma distribution.

### Extracting ages and transforming proxy records

Ages of new stratigraphic positions can be estimated using `agePredict`:

```{r}
#| results: hide
# estimate ages every 50 cm, without height uncertainty
new_positions = seq(from = 49, to = 190, by = 0.5)
age_predictions = modifiedBChron::agePredict(model = age_model, 
                                            newPositions = new_positions, # height of new positions
                                            newPositionThicknesses = rep(0, length(new_positions))) # thicknesses
```

Now we hand these estimated median ages to the `admtools` package:

```{r}
library(admtools)
adm = admtools::tp_to_adm(t = age_predictions$HDI$`0.5`,
                          h = age_predictions$HDI$newPositions,
                          T_unit = "Myr",
                          L_unit = "m")
plot(adm)
admtools::T_axis_lab(label = "Age")
admtools::L_axis_lab()
title(main = "Median Age in SCF")
```

As an example proxy record, we simulate a random walk using the `StratPal` package [@Hohmann2025_stratpal]:

```{R}
library(StratPal)
set.seed(42)
proxy_record = StratPal::random_walk(t = seq(from = min_time(adm),
                                             to = max_time(adm),
                                             by = 0.01),# 10 kyr resolutions
                                     sigma = 1,
                                     mu = 1) 
plot(proxy_record,
     type = "l",
     xlab = "Age [Myr]",
     ylab = "Proxy value",
     main = "Synthetic proxy record in the time domain")
```

This proxy record is currently in the time domain, and can be transformed into the stratigraphic domain using `time_to_strat`:

```{R}
proxy_record_height = admtools::time_to_strat(obj = proxy_record,
                                              x = adm)
plot(proxy_record_height,
     orientation = "du", # plot up/down
     type = "l",
     xlab = "Proxy value",
     ylab = "Height [m]",
     main = "Synthetic proxy record in the straitgraphic domain")
```

Here the distorting effects of the age-depth model are clearly visible: the low sedimentation rate at around 80 m (17.2 Ma) results in many short oscillations, while the interval from 140 to 100 meter is stratigraphically diluted, resulting in fewer oscillations than in the time domain.

We can get the sedimentation rate in the stratigraphic domain from the median age-depth model using `sed_rate_l`:

```{R}
plot(x = new_positions,
     y = abs(admtools::sed_rate_l(adm, new_positions)),
     type = "l",
     xlab = "Height [m]",
     ylab = "Sedimentation rate [m / Myr]",
     ylim = c(0, 310),
     main = "Sed. rate from median ADM")
```

Let's say we're interested in the amount of time recorded between 100 and 120 m height. We can determine this using

```{r}
#| results: hide
age_predictions2 = modifiedBChron::agePredict(model = age_model, 
                                              newPositions = c(100, 120), # height of new positions
                                              newPositionThicknesses = rep(0, 2)) # thicknesses
```

```{R}
diff = age_predictions2$raw$V1 - age_predictions2$raw$V2
hist(diff,
     xlab = "Time [Myr]",
     freq = TRUE,
     main = "Time between 100 and 120 m")
```

### Tasks

1. Plot the condensation in the section (tip: use `admtools::condensation` )
2. Add another tie point at around 110 m, and examine how the uncertainties and the downstream analyses change


### Solutions
::: {.callout-note collapse="true"}
#### Solution task 1

Task: Plot the condensation in the section


To plot the condensation in the section, we use `admtools::condensaton`:
```{R}
plot(x = new_positions,
     y = abs(admtools::condensation(adm, new_positions)),
     type = "l",
     xlab = "Height [m]",
     ylab = "Condensation [Myr / m]",
     ylim = c(0, 0.03),
     main = "Condensation from median ADM")
```
:::

::: {.callout-note collapse="true"}
#### Solution task 2
Task: Add another tie point at around 110 m, and examine how the uncertainties and the downstream analyses change

For the tie point at 110 m, we assume a normally distributed age ("G" for Gaussian) with a mean age of 17.1 Ma and a standard deviation of 0.8 Ma and no height uncertainty:

```{r}
#| message: false
#| output: false
levels(SCF$ids) = c(levels(SCF$ids), "NewID")
age_model_mod = modifiedBChron::ageModel(ages = c(SCF$age, 17.1),
                                         ageSds = c(SCF$ageSds, 0.8),
                                         positions = c(SCF$position, 110),
                                         ids = factor(c(as.character(SCF$ids), "newID")),
                                         positionThicknesses = c(SCF$thickness, 0),
                                         distTypes = factor(c(as.character(SCF$distType), "G")),
                                         predictPositions = seq(48.25, 193.50, by = 1),
                                         MC = 10000,  # how many iterations 
                                         burn = 2000) # how many iterations to discard
```

```{r}
modifiedBChron::modelPlot(model = age_model_mod,
                          scale = 5)
```

We focus on the newly introduced height 110 m

```{r}
#| message: false
#| output: false
new_positions = 110
age_predictions = modifiedBChron::agePredict(model = age_model_mod, 
                                              newPositions = new_positions, # height of new positions
                                              newPositionThicknesses = rep(0, length(new_positions))) # thicknesses


```

```{r}
age_predictions$HDI
```
Here you can see that the age uncertainty at the newly introduced tie point is strongly reduced compared to the assumed large standard deviation.
:::

## astroBayes: Astrochronology and radioisotopic geochrononlogy

`astroBayes` is an R package for combining astrochronology and radioisotopic geochronology published by @trayler2023. It creates age-depth models by simultaneously matching a proxy signal to target frequencies and absolute radiometric dates.

### Usage

This code is modified from [this repository](https://github.com/robintrayler/astroBayes).

First, let's load the package and the associated data:

```{r}
#| warning: false
library(astroBayes)
# target frequencies for matching
# taken from LA04 (prec & obliquity) and LA10d (LA10d)
data("target_frequencies") 
# absolute age constraints
# is available as `dates` in the workspace
data("radioisotopic_dates")
# cyclostratigraphic data
data("cyclostratigraphic_data")
# boundaries of the stratigraphic layers
data("layer_boundaries")
```

To construct the age-depth model, pass the data to `astro_bayes_model`. Warning: the estimation can take a few minutes.

```{r}
#| message: false
age_model = astroBayes::astro_bayes_model(geochron_data = dates,
                                          cyclostrat_data = cyclostrat,
                                          target_frequency = target_frequencies,
                                          layer_boundaries = layer_boundaries,
                                          iterations = 10000,
                                          burn = 1000)
```

The age-depth model can be plotted using

```{R}
plot(age_model,type = "age_depth")
```

Estimates for the sedimentation rates can be shown using

```{R}
plot(age_model, type = "sed_rate")
```

To check the match of the age-depth model to the target frequencies, use

```{r}
#| warning: false
plot(age_model, type = "periodogram")
```

This will show a periodogram of the data, with the target frequencies shown as dashed vertical lines.

### Model Assumptions

`astroBayes` fits a proxy record (specified by the variable `cyclostrat_data`) to a set of target frequencies (specified by `target_frequency`) and radioisotopic data (specified by `geochron_data`). It assumes constant sedimentation rate (prior specified by `layer_boundaries$sed_min` and `layer_boundaries$sed_max`) over stratigraphic intervals specified by `layer_boundaries$position`, and uses the radioisotopic dates to anchor this age-depth model. The age-depth model can also include hiatus surfaces. For a visual overview of these parameters, see Figure 1 in @trayler2023

### Extracting ages

To get new ages, use `predict` with a suitable formatted data frame

```{R}
new_positions = data.frame(position = c(5, 10, 15), 
                           thickness = c(1, 0, 1), 
                           id = c("x", "y", "z"))
predictions = predict(age_model, new_positions)
plot(predictions)
```

### Example: CIP 2.0 case study

As example for the application of `astroBayes`, we reproduce the CIP case study number 1 [@SINNESAEL2019] , code is based on [this repository](https://github.com/robintrayler/astroBayes_manuscript).

```{R}
library(dplyr)
# load data
cyclostrat        = read.csv(file = "./data/CIP2/cyclostrat_data.csv")
target_frequency  = read.csv(file = "./data/CIP2/target_frequency.csv")
true_data         = read.csv(file = "./data/CIP2/true_age.csv")
layer_boundaries  = read.csv(file = "./data/CIP2/layer_boundaries.csv") 

# get true ages of 4 positions
date_positions = true_data[c(125, 350, 700, 950), ]
  
geochron_data  = # assemble into synthetic radiometric dates
  data.frame(age = date_positions$age,
             age_sd = date_positions$age * 0.015,
             position = date_positions$position,
             thickness = 0,
             id = letters[1:4])  |> 
  arrange(position)
```

Now we can estimate the age-depth model (again, this can take a few minutes)

```{r}
#| message: false
 model = astroBayes::astro_bayes_model(geochron_data = geochron_data,
                                       cyclostrat_data = cyclostrat,
                                       target_frequency = target_frequency,
                                       layer_boundaries = layer_boundaries,
                                       iterations = 10000, 
                                       burn = 1000, 
                                       method = "malinverno")
```

The plot immediately shows the hiatus surface at around 6 m:

```{R}
plot(model, type = "age_depth")
```

Now we extract the median age-depth model

```{R}
new_positions = data.frame(position = cyclostrat$position, 
                           thickness = rep(0, length(cyclostrat$position)), 
                           id = seq_along(cyclostrat$position))
predictions = predict(model, new_positions)
```

and plot it:

```{R}
adm = admtools::tp_to_adm(t = predictions$CI$median,
                          h = predictions$CI$position,
                          L_unit = "m",
                          T_unit = "Myr")
plot(adm)
admtools::T_axis_lab(label = "Age")
admtools::L_axis_lab(label = "Depth")

```

Now we can convert the observed cyclostratigraphic signal back into the time domain;

```{R}
cyc_signal = list(h = cyclostrat$position,
                  y = cyclostrat$value)
cyc_signal_time = admtools::strat_to_time(cyc_signal, adm)
plot(cyc_signal_time,
     type = "l",
     xlab = "Age [Myr]",
     ylab = "Value")
```

The destructive effect of the hiatus is clear at around 1.1 Myr

### Tasks

1. Modify the prior on the sedimentation rate (e.g., broaden or shift the interval). How does that change the estimated age-depth model? What happens if you provide unrealistically high/low sedimentation rates?
2. Remove the information on the hiatus from the data. How do your estimates change? Compare the periodograms of the data with and without the hiatus. How do your age estimates for the age of 2 m depth change?

### Solutions 

::: {.callout-tip collapse="true"}
#### Solution task 1

Task: Modify the prior on the sedimentation rate (e.g., broaden or shift the interval). How does that change the estimated age-depth model? What happens if you provide unrealistically high/low sedimentation rates?

We construct two cases, one over and the other underestimating the sedimentation rates:

```{r}
# sed rate priors from 1 to 4
layer_boundaries_under = layer_boundaries
layer_boundaries_under$sed_min = layer_boundaries_under$sed_min - 4
layer_boundaries_under$sed_max = layer_boundaries_under$sed_max - 16
# sed rate priors from 15 to 30
layer_boundaries_over = layer_boundaries
layer_boundaries_over$sed_min = layer_boundaries_over$sed_min + 10
layer_boundaries_over$sed_max = layer_boundaries_over$sed_max + 10
```

Underestimating the sedimentation rates:
```{r}
#| message: false
 age_model_under = astroBayes::astro_bayes_model(geochron_data = geochron_data,
                                                  cyclostrat_data = cyclostrat,
                                                  target_frequency = target_frequency,
                                                  layer_boundaries = layer_boundaries_under,
                                                  iterations = 10000, 
                                                  burn = 1000, 
                                                  method = "malinverno")
```

Plotting this, we get

```{r}
plot(age_model_under, type = "age_depth")
```

and 
```{r}
plot(age_model_under, type = "sed_rate")
```


Whereas with overestimating the sedimentation rates, we get
```{r}
#| message: false
 age_model_over = astroBayes::astro_bayes_model(geochron_data = geochron_data,
                                                cyclostrat_data = cyclostrat,
                                                target_frequency = target_frequency,
                                                layer_boundaries = layer_boundaries_over,
                                                iterations = 10000, 
                                                burn = 1000, 
                                                method = "malinverno")
```

Plotting this, we get

```{r}
plot(age_model_over, type = "age_depth")
```

and 

```{r}
plot(age_model_over, type = "sed_rate")
```

:::

::: {.callout-tip collapse="true"}

#### Solution task 2

Task: Remove the information on the hiatus from the data. How do your estimates change? Compare the periodograms of the data with and without the hiatus. How do your age estimates for the age of 2 m depth change?

We keep all the layers but remove the info that one of the layer boundaries is a hiatus:

```{r}
layer_boundaries_no_hiat = layer_boundaries
layer_boundaries_no_hiat$hiatus_boundary[3] = FALSE
```

Estimating the age-depth model results in

```{r}
#| message: false
 age_model_no_hiat = astroBayes::astro_bayes_model(geochron_data = geochron_data,
                                                  cyclostrat_data = cyclostrat,
                                                  target_frequency = target_frequency,
                                                  layer_boundaries = layer_boundaries_no_hiat,
                                                  iterations = 10000, 
                                                  burn = 1000, 
                                                  method = "malinverno")
```


Plotting this, we get

```{r}
plot(age_model_no_hiat, type = "age_depth")
```

Comparing the age estimates at 2 m we get

```{r}
#| message: false
poi = data.frame(position = 2,
                 thickness = 0,
                 id = as.character(1))
predictions_w_hiat = predict(model, poi)
predictions_no_hiat = predict(age_model_no_hiat, poi)
```

```{r}
# no hiatus
predictions_no_hiat$CI$CI_2.5
predictions_no_hiat$CI$CI_97.5
# with hiatus
predictions_w_hiat$CI$CI_2.5
predictions_w_hiat$CI$CI_97.5
```
Surprisingly, removing the hiatus has upcore effects even on intervals that are not adjacent to the hiatus.
:::


## admtools: Age-depth models from sedimentation rates

The `admtools`package [@hohmann_nonparametric_2024] provides two different methods to estimate age-depth models. Here, we focus on the estimation of age-depth models from information on sedimentation rates (e.g., from evolutive methods, or qualitative constraints from sequence stratigraphy). The second methods is discussed in detail in the section below.

### Usage

We demonstrate the usage of the admtools package by reproducing the example from @hohmann_nonparametric_2024 , which builds on @dasilva2020 (code available in @dasilva2024). In this example, first sedimentation rates are estimated using eTimeOpt @meyers2019. Then, uncertainties on sedimentation rates are propagated into the age-depth model, which is anchored using dates from @percival2018

First, we load the required data from @dasilva2020

```{R}
data=read.csv("./data/raw/SbS_XRF_forfactor3.csv",header = T, sep=";")
```

First, let's clean that data and define some constants for the analysis:

```{r}
#| message: false
#| warning: false
Height=data[,c(1)]
mydata=data[,c(2:23)]
mydata[is.na(mydata)] = 0
# select magnetic susceptibility
MS=data[,c(1,23)]
MS=na.omit(MS)
MSi=astrochron::linterp(dat = MS,
                        dt = 0.02,
                        verbose = FALSE,
                        genplot = FALSE)

# target frequencies for the matching
targetE=c(130.7, 123.8, 98.9, 94.9)
astrochron::bergerPeriods(372, genplot = FALSE)
targetP=c(19.98, 16.87)

# ages and uncertainties from Percival
t_ashbed_mean = 372.360 * 1000
uncertainty_radiometric = 0.053 # [Ma], 2 sigma around the age
t_ashbed_sd = uncertainty_radiometric * 1000 / 2 # Percival uses 2 sigma

# Heights of the ash beds, lower and upper kellwasser events, and Frasnian-Famennian boundary
h_ashbed = 1.55 
h_top_lkw = 1.15
h_bottom_ukw = 3.6
h_top_ukw = 4.05
h_f_f_bdry = 4.05
```

Here, we focs only on the case with precession amplitude modulation - for short eccentricity amplitude modulation see @hohmann_nonparametric_2024. We start by estimating the sedimentation rates using eTimeOpt

```{r}
etimeOptMS_prec=astrochron::eTimeOpt(MSi,
                                     win=0.02*100,
                                     step=0.02*10,
                                     sedmin=0.1,
                                     sedmax=0.6, 
                                     numsed=100,
                                     linLog=1,
                                     limit=T,
                                     fit=1,
                                     fitModPwr=T,
                                     flow=NULL,
                                     fhigh=NULL,
                                     roll=NULL,
                                     targetE=targetE,
                                     targetP=targetP,
                                     detrend=T,
                                     ydir=1,
                                     output=1,
                                     genplot=T,
                                     check=T,
                                     verbose=0)
```

This provides us with estimates of r\^2 envelope, power, and their product as a function of sedimentation rate and height. We focus on estimating the age-depth model from the product (third plot)

```{r}
# extract r^2_opt from eTimeOpt
fa_prec = admtools::get_data_from_eTimeOpt(etimeOptMS_prec, index = 3)
# generate sed rate generator
rate = 3 # rate parameter: no. of changes in sed. rate per stratigraphic unit

se_prec = admtools::sed_rate_from_matrix(height = fa_prec$heights,
                                         sedrate = fa_prec$sed_rate / 100, # convert cm/kyr to m/kyr
                                         matrix = fa_prec$results,
                                         mode = "poisson",
                                         rate = rate)
```

`se_prec` is a function that upon each evaluation returns one possible sedimentation rate in the section. Using a Monte Carlo method, this can be converted into an age-depth model.

Next, we define the stratigraphic tie points.

```{r}
# stratigraphic tie points: height of the tie point
h_tp = function(){
  return(h_ashbed)
}

# absolute age tie point with age from Percival et al 2018
t_tp_absolute = function(){
  t = rnorm(1, mean =  - t_ashbed_mean, sd = t_ashbed_sd)
  return(t)
}
```

With the information on sedimentation rates and tie points at hand, we can now estimate the age-depth model by calling the function

```{r}
no_of_rep = 1000
h = seq(1.05, 4.25, by = 0.05) # height to estimate ages - every 5 cm
adm_prec_abs = admtools::sedrate_to_multiadm(h_tp , 
                                             t_tp_absolute,
                                             sed_rate_gen = se_prec,
                                             h, 
                                             no_of_rep = no_of_rep )
```

```{r}
plot(adm_prec_abs)
mtext("Time BP", side = 1, line = 3)
mtext("Stratigraphic Position", side = 2, line = 3)
lines(x = c(-10^9, 0), y = c(h_ashbed, h_ashbed))
```

Here the line is the stratigraphic position of the ash layer.

### Model assumptions

A priori, `admtools` assumes that the law of superposition holds, and that there are no age reversals between tie points. However, here we import data from `eTimeOpt` using the `sed_rate_from_matrix` function, which assumes sedimentation rates follow a probability distribution specified by r\^2_opt, and that there are breakpoints following a Poisson process according to the specified `rate` parameter.

### Determining ages

As an application, we determine the age of the Frasnian-Famennian boundary using `get_time`:

```{r}
ages = admtools::get_time(adm_prec_abs, h = h_f_f_bdry)
hist(abs(unlist(ages))/1000,
     xlab = "Age [Ma]",
     freq = FALSE,
     main = "Age of F-F boundary")
```

In addition, we determine the duration of the Upper Kellwasser event

```{r}
dur = admtools::get_time(adm_prec_abs, h = c(h_top_ukw, h_bottom_ukw))
duration = sapply(dur, function(x) x[1] - x[2])
hist(duration,
     xlab = "Duration [kyr]",
     freq = FALSE,
     main = "Duration Upper Kellwasser Event")
```

### Tasks

1.   Modify the time tie point to have no uncertainty. How does this change the shape of the age-depth model? How does this propagate into the downstream analyses (e.g., duration of the Kellwasser events, age of the F-F boundary)?

### Solutions

::: {.callout-tip collapse="true"}
#### Solution task 1
Task: Modify the time tie point to have no uncertainty. How does this change the shape of the age-depth model? How does this propagate into the downstream analyses (e.g., duration of the Kellwasser events, age of the F-F boundary)?

We keep the same mean age, but remove all uncertainty:

```{r}
t_tp_no_uncertainty = function(){
  t = rnorm(1, mean =  - t_ashbed_mean, sd = t_ashbed_sd)
  return(t)
}
```

Re-estimating the age-depth model, we get
```{r}
adm_tp_no_uncertainty = admtools::sedrate_to_multiadm(h_tp , 
                                             t_tp_no_uncertainty,
                                             sed_rate_gen = se_prec,
                                             h, 
                                             no_of_rep = no_of_rep )
```


Plotting this
```{r}
plot(adm_tp_no_uncertainty)
mtext("Time BP", side = 1, line = 3)
mtext("Stratigraphic Position", side = 2, line = 3)
lines(x = c(-10^9, 0), y = c(h_ashbed, h_ashbed))
```

The reduction of tie point uncertainty does reduce the uncertainty of absolute ages, with the error now being only due to uncertainties on sedimentation rates:
```{r}
ages = admtools::get_time(adm_tp_no_uncertainty, h = h_f_f_bdry)
hist(abs(unlist(ages))/1000,
     xlab = "Age [Ma]",
     freq = FALSE,
     main = "Age of F-F boundary")
```

In contrast, there is no influence on  the estimation of durations, as uncertainties cancel out:

```{r}
dur = admtools::get_time(adm_tp_no_uncertainty, h = c(h_top_ukw, h_bottom_ukw))
duration = sapply(dur, function(x) x[1] - x[2])
hist(duration,
     xlab = "Duration [kyr]",
     freq = FALSE,
     main = "Duration Upper Kellwasser Event")
```

:::

## admtools: Age-depth models from tracer values

The `admtools` package [@hohmann_nonparametric_2024] provides a second way to estimate age-depth models, which is based on the comparison of observed tracer values (e.g., Pollen, extraterrestrial Helium 3, 210Pb) with assumptions of their flux into the sediment. This procedure is called FAM for Flux - Assumption Matching.

### Usage

To demonstrate this approach, we reproduce the age-depth model from @murphy_extraterrestrial_2010 for the Paleocene-Eocene Thermal Maximum (PETM) from Walvis Ridge, IODP Site 1266. This example is taken from @hohmann_nonparametric_2024, code is available under @hohmann_2025_15489276

We begin by importing data from the publication:

```{R}
# load data from Murphy et al. https://doi.org/10.1016/j.gca.2010.03.039
data = read.csv(file = "./data/raw/murphy_et_al_2010_1-s2.0-S0016703710003108-mmc3.csv",
                header = TRUE,
                sep = "\t")
```

For the analysis, we need some constants:

```{R}
# heights
h = data$Depth..mcd.

# base of clay layer, from Murphy et al., Table 1
base_clay = 306.78
# heights where the ages are determined - cm resolution
h_eval = seq(from = 303.5,
             to = base_clay,
             by = 0.01)
# 3He flux mean and uncertainty             
mean_3He_flux = 0.48 # mu based on supplementary materials
twosigma_3H3_flux = 0.08 # 2 sigma based on supplementary materials from Murphy et al.
```

Next, we define the stratigraphic tie points. As we do not have any absolute ages, we construct a floating age-depth model, where time is measured relative to the age of the base of the clay layer. For details on how more complex tie points can be constructed see the documentation in the `admtools` package [webpage](https://mindthegap-erc.github.io/admtools/articles/adm_from_trace_cont.html).

```{R}
# stratigraphic tie point at the base of the clay layer
h_tp = function(){
  return(base_clay)
}
# time tie point: measure time relative to the deposition of the base of the clay layer
t_tp = function(){
  return(0)
}
```

Next we define our assumed and observed 3He fluxes. The assumption made is that the tracer flux in the time domain is constant (but uncertain), so increased observed fluxed in the stratigraphic domain reflect low sedimentation rates. Accumulating this information can used to construct an age-depth model.

```{R}
# assumed flux in the time domain
# flux is constant, with an error according to the flux values determined by Murphy et al. (2010)
time_const_gen_supp = function(){
  eps = 0.0001
  r = max(eps,rnorm(1, mean = mean_3He_flux, sd = twosigma_3H3_flux/2)) # murphy et al 2010, pcc/cm^-2/kyr
  f = approxfun(x = c(-1000, 1000), y = rep(r, 2), rule = 2)
  return(f)
}
```

```{R}
# 3He flux observed in the stratigraphic domain 
# Based on data from Murphy et al. 2010
strat_cont_gen_det = function(){
  f = approxfun(x = data$Depth..mcd.,
                y = data$X3HeET..pcc.g...20. * data$DBD..g.cm3. * 100,
                rule = 2)
  return(f)
}
plot(x = h, 
     y = strat_cont_gen_det()(h),
     xlab = "Depth [m]",
     ylab = "3He flux [pcc/cm^-2/kyr]",
     main = "Observed 3He flux",
     type = "l")
```

This already indicates that something interesting is happening: If 3He flux is constant in time, the peaks in observed tracer flux correspond to stratigraphic condensation (low sedimentation rates).

Now we can estimate the age-depth model by calling `strat_cont_to_multiadm`. This can again take a few minutes to calculate.

```{R}
# height where the adm is estimated
h_est = seq(304, base_clay, by = 0.05)         
subdiv = 1000 # numeric options for integration
adm = admtools::strat_cont_to_multiadm(h_tp = h_tp,
                               t_tp = t_tp,
                               strat_cont_gen = strat_cont_gen_det,
                               time_cont_gen = time_const_gen_supp,
                               h = h_est,
                               subdivisions = subdiv,
                               stop.on.error = FALSE,
                               no_of_rep = 500)
```

This age-depth model can be plotted by calling the standard `plot` function.

```{R}
plot(adm)
mtext("Time before clay layer [kyr]", side = 1, line = 3)
mtext("Depth [mbsf]", side = 2, line = 3)
```

Here, the red line is the median age, and the blue lines are the 95 % coverage interval.

### Model assumptions

A priori, `strat_cont_to_multiadm` makes no model assumptions other than that the law of superposition holds and the tracer values are positive. Reversals between tie points must be resolved by the user, which is why this method is most powerful when few tie points are present. The procedure is based on a simple "pattern matching" approach: If it is know how much tracer is deposited in a certain amount of time (assumption), then this can be used to constrain how much time is recorded in the corresponding stratigraphic interval (observation). This can in turn be used to construct age-depth models.

### Extracting ages

The `adm` object contains a number of realizations of the age-depth model, each of which is one potential chronology. From these chronologies, a variety of data can be extracted.

#### Transforming data

First, we extract the median age-depth model:

```{R}
adm_m = admtools::median_adm(adm, h = h_est)
plot(adm_m)
mtext("Time before clay layer [kyr]", side = 1, line = 3)
mtext("Depth [mbsf]", side = 2, line = 3)
title(main = "Age-depth model for the PETM at ODP site 1266")
```

We can get the sedimentation rate in the stratigraphic domain using `sed_rate_l`:

```{R}
plot(x = h_est, 
     y = admtools::sed_rate_l(adm_m, h_est) * 100,
     xlab = "Time [kyr]",
     ylab = "Sedimentation rate [cm/kyr]",
     type = "l",
     main = "Sedimentation rate at OPD site 1266")
```

Note the low sedimentation rate at the base of the clay layer (around 306.8 m) which coincides with the PETM main interval.

To study how this condensation affects proxy records, we simulate a synthetic proxy record over 200 kyr using a random walk model.

```{R}
set.seed(42)
proxy_record = StratPal::random_walk(t = seq(from = 0,
                                             to = 250,
                                             by = 1),# 10 kyr resolutions
                                     sigma = 0.1,
                                     mu = 0.01) 
proxy_record$t = - proxy_record$t # reverse time direction
plot(proxy_record,
     type = "l",
     xlab = "Time to clay layer [kyr]",
     ylab = "Proxy value")
```

Transforming this record into the stratigraphic domain yields

```{R}
proxy_record_height = admtools::time_to_strat(proxy_record, adm_m)
plot(proxy_record_height,
     type = "l",
     xlab = "Proxy value",
     ylab = "Height [m]")
```

The effect of the condensation at the PETM main interval coinciding with the clay layer is clearly visible, with the increased volatility of the proxy record at the top of the section (which is not causally linked to the PETM).

#### Duration of the PETM

Now we extract some durations from the section. The section was subdivided by into multiple intervals, we follow the intervals given by @murphy_extraterrestrial_2010

```{R}
clay_layer_top = 306.15
base_recovery = 306.4
recovery_1_top = 306.15 # “Shoulder” δ13C inflection point F
recovery_2_top = 304.7 # δ13C inflection point G
recovery_3_top = 304.19 # End of anomalously high carbonate sedimentation
```

We get the durations of the recovery interval by using `get_time`:

```{R}
recovery_times = admtools::get_time(adm, h = rev(c(base_recovery, recovery_3_top)))
diff = sapply(recovery_times, diff)
hist(diff,
     main = "Duration PETM recovery interval",
     xlab = "Duration [kyr]",
     freq = FALSE)

```

This gives us a median duration of the PETM recovery interval of 100 kyr at ODP site 1266.

### Tasks

1.   What is the duration of the three recovery intervals?
2.   Double the assumed 3He flux. How does that change the estimates of the PETM duration?
3.   Assume the onset of the clay layer is at exactly 55.9 Ma, with an uncertainty of 0.2 Myr (2 sigma). Add this as a tie point to the analysis, and re-estimate the age-depth model.
4.   How does this affect the downstream analyses (e.g., sedimentation rates, median age-depth models, durations)? Which ones are affected, which ones are not? Why?

### Solutions

::: {.callout-tip collapse="true"}
#### Solution task 1

Task: What is the duration of the three recovery intervals?

To calculate the duration of the the three recovery intervals, we reuse the code above and modify the stratigraphic interval examined:

```{r}
# 1st recover interval
recovery_interval_1 = c(base_recovery, recovery_1_top)
recovery_times = admtools::get_time(adm, h = rev(recovery_interval_1))
diff = sapply(recovery_times, diff)
hist(diff,
     main = "Duration PETM recovery interval 1",
     xlab = "Duration [kyr]",
     freq = FALSE)

```
```{r}
# 2nd recover interval
recovery_interval_2 = c(recovery_1_top, recovery_2_top)
recovery_times = admtools::get_time(adm, h = rev(recovery_interval_2))
diff = sapply(recovery_times, diff)
hist(diff,
     main = "Duration PETM recovery interval 2",
     xlab = "Duration [kyr]",
     freq = FALSE)

```

```{r}
# 3rd recovery interval
recovery_interval_3 = c(recovery_2_top, recovery_3_top)
recovery_times = admtools::get_time(adm, h = rev(recovery_interval_3))
diff = sapply(recovery_times, diff)
hist(diff,
     main = "Duration PETM recovery interval 3",
     xlab = "Duration [kyr]",
     freq = FALSE)

```

:::

::: {.callout-tip collapse="true"}

#### Solution task 2

Task: Double the assumed 3He flux. How does that change the estimates of the PETM duration?

To model this, we multiply the initially modeled 3He flux by 2, and rerun the age-depth model:
```{r}
double_mean_3He_flux = 2 * mean_3He_flux
time_const_gen_supp_doubled = function(){
  eps = 0.0001
  r = max(eps,rnorm(1, mean = double_mean_3He_flux, sd = twosigma_3H3_flux/2)) # murphy et al 2010, pcc/cm^-2/kyr
  f = approxfun(x = c(-1000, 1000), y = rep(r, 2), rule = 2)
  return(f)
}

adm_doubled_flux = admtools::strat_cont_to_multiadm(h_tp = h_tp,
                               t_tp = t_tp,
                               strat_cont_gen = strat_cont_gen_det,
                               time_cont_gen = time_const_gen_supp_doubled,
                               h = h_est,
                               subdivisions = subdiv,
                               stop.on.error = FALSE,
                               no_of_rep = 500)
```

Plotting it again, it becomes immediately obvious that duration of the covered age-depth model has halved:

```{R}
plot(adm_doubled_flux)
mtext("Time before clay layer [kyr]", side = 1, line = 3)
mtext("Depth [mbsf]", side = 2, line = 3)
```
Extracting the durations, we see that the estimated recovery duration has also halved:

```{R}
recovery_times = admtools::get_time(adm_doubled_flux, h = rev(c(base_recovery, recovery_3_top)))
diff = sapply(recovery_times, diff)
hist(diff,
     main = "Duration PETM recovery interval, halved 3He flux",
     xlab = "Duration [kyr]",
     freq = FALSE)

```

:::


::: {.callout-tip collapse="true"}
#### Solution task 3

Task: Assume the onset of the clay layer is at exactly 55.9 Ma, with an uncertainty of 0.2 Myr (2 sigma). Add this as a tie point to the analysis, and re-estimate the age-depth model.

For this, we modify the tie point age to match the age assumptions:

```{r}
t_tp_mod = function(){
  x = rnorm(1, mean = -55900, sd = 100) # unit conversion into kyr & 1 sigma!
  return(x)
}
```


Now we can re-estimate the age-depth model

```{R}
# height where the adm is estimated
h_est = seq(304, base_clay, by = 0.05)         
subdiv = 1000 # numeric options for integration
adm_t_tp_mod = admtools::strat_cont_to_multiadm(h_tp = h_tp,
                               t_tp = t_tp_mod,
                               strat_cont_gen = strat_cont_gen_det,
                               time_cont_gen = time_const_gen_supp,
                               h = h_est,
                               subdivisions = subdiv,
                               stop.on.error = FALSE,
                               no_of_rep = 500)
```

and plot it:

```{R}
plot(adm_t_tp_mod)
mtext("Time before present [kyr]", side = 1, line = 3)
mtext("Depth [mbsf]", side = 2, line = 3)
title(main = "ADM in absolute age")
```
Not the high uncertainty in the age-depth model now. This uncertainty is dominated by the uncertainty at the tie point rather the uncertainty introduced by the 3He measurements or the fluctuations in 3He flux.
However, calculating the duration of the PETM recovery interval shows that it remains identical: The uncertainty of the tie point does not propagate into durations.

```{r}
recovery_times = admtools::get_time(adm_t_tp_mod, h = rev(c(base_recovery, recovery_3_top)))
diff = sapply(recovery_times, diff)
hist(diff,
     main = "Duration PETM recovery interval",
     xlab = "Duration [kyr]",
     freq = FALSE)
```
:::

::: {.callout-tip collapse="true"}
#### Solution task 4

Task: How does this affect the downstream analyses (e.g., sedimentation rates, median age-depth models, durations)? Which ones are affected, which ones are not? Why?

The uncertainty of the tie point only matters for the absolute age of events, e.g., when the second recovery interval started.

The median age-depth model is only shifted by the median absolute age, and all analyses that are based on durations (e.g., duration of the recovery intervals) are not affected because the uncertainty of the tie point cancels out.
:::

## References