`#> Loading required package: ape`

This vignette explains how to estimate age-depth models (ADMs) from
sedimentation rates and tie points using the function
`sedrate_to_multiadm`

The `sedrate_to_multiadm`

function estimates age-depth
models from sedimentation rates and position and timing of tie points.
It takes the following inputs that encode user knowledge:

`h_tp`

: a function encoding stratigraphic positions of tie points`t_tp`

: a function encoding times of the tie points`sed_rate_gen`

: a function encoding how sedimentation rates change with stratigraphic positions. This information can for example be derived from cyclostratigraphic analyses.

In addition, it takes the following inputs that specify the estimation procedure:

`h`

: vector of heights where the ADM is determined`no_of_rep`

: integer, number of runs

Additional parameters that determine the numeric behavior of the
integration method used are `subdivisions`

and
`stop.on.error`

. The parameters `T_unit`

and
`L_unit`

can be used to associate time and length units with
the generated age-depth model.

We construct an age depth model for a section of 10 m thickness where upper and lower bounds on sedimentation rates are available.

We start by defining the section as well as lower and upper bounds on sedimentation rates

```
h_min = 2 # lower boundary of the section
h_max = 10 # upper boundary of the section
T_unit = "Myr"
L_unit = "m"
```

We want to know the age-depth model every 10 cm, so we define

We assume there is one tie point in the section at 5 m height, and its mean age is 66 years with a standard deviation of 0.25 Myr

Then the tie point timing is given by

and the tie point height is given by

Every time `t_tp`

is evaluated, it returns one possible
time of the tie points. Similarly, every time `h_tp`

is
evaluated, it returns the stratigraphic position of the tie point (which
is deterministic in this case):

We assume the following upper and lower limits for sedimentation rates:

```
sedrate_max_y = c(2,5,8,5)
sedrate_max_x = c(1,4,6,10)
sedrate_min_y = c(1,1,7,0.5)
sedrate_min_x = sedrate_max_x
```

Here, `sedrate_max_y[i]`

is the upper limit on
sedimentation rate at `sedrate_max_x[i]`

(mutatis mutandis
for `sedrate_min`

). Between these points, we assume linear
interpolation. This is done by the function
`sed_rate_gen_from_bounds`

:

```
sedrate = sed_rate_gen_from_bounds(h_l = sedrate_min_x,
s_l = sedrate_min_y,
h_u = sedrate_max_x,
s_u = sedrate_max_y,
rate = 1)
```

Because the sedimentation rates are uncertain, functions returned by
`sedrate`

will differ each time the function is evaluated. As
an example, we plot three different sample paths (realizations) of the
sedimentation rate through the section:

```
plot(NULL,
xlim = range(h),
ylim = c(0, max(c(sedrate_max_y))),
xlab = "Height [m]",
ylab = "Sedimentation Rate [m/Myr]")
no_sedrates = 3
cols = c("red", "blue", "black")
for (i in seq_len(no_sedrates)){
sedrate_sample = sedrate()
lines(h, sedrate_sample(h), lwd = 3, col = cols[i])
}
```

These sedimentation rates assume sedimentation rates are drawn from a
uniform distribution between the upper and lower limits of sedimentation
rate provided by the user. This is done at random locations determined
according to a Poisson point process with rate `rate`

.

Times and heights of tie points are coded via the functions
`t_tp`

(timing) and `h_tp`

(height) that take no
inputs. They serve as wrappers around user-defined procedures that
reflect uncertainties around tie points. Every time `t_tp`

and `h_tp`

are evaluated, they return possible values for the
tie points. Conceptually, both `t_tp`

and `h_tp`

are user implemented random number generators that draw from the
distributions of tie points. Writing these functions requires some
effort, but it allows the user to hand over arbitrarily complex
uncertainties of the tie points to the
`strat_cont_to_multiadm`

function.

Multiple wrappers are available to simplify coding tie points:

`tp_height_det`

for specifying deterministic stratigraphic heights`tp_time_floating_scale`

to encode time tie points for floating time scale`tp_time_norm`

for normally distributed tie points in time

Both `t_tp`

and `h_tp`

must return strictly
ordered numeric vectors of times/heights. This means that it is the
users responsibility to avoid inversions of times/heights.

As an example, I assume the stratigraphic positions of the tie points are known without uncertainty, and are at 10 and 20 m stratigraphic height.

```
h_min = 10 # stratigraphic height of lower tie point [m]
h_max = 20 # stratigraphic height of upper tie point [m]
```

`h_tp`

is then implemented as follows:

When evaluated, this function returns the stratigraphic positions of the tie points:

For a more complex example, I assume that the timing of the first tie
point follows a normal distribution with mean 0 and standard deviation
0.5. For the second tie point, only maximum and minimum time is
available. Due to the lack of information, I assume a uniform
distribution between the minimum (9) and the maximum (11). This is
implemented as follows in `t_tp`

:

```
t_tp = function() {
repeat{
# timing first tie point
t1 = rnorm(n = 1, mean = 0, sd = 0.5)
# timing second tie point
t2 = runif(n = 1, min = 9, max = 11)
if (t1 < t2){ # if order is correct, return values
return(c(t1, t2))
}
}
}
```

```
t_tp() # evaluating the function returns a random pair of times drawn from the specified distribution
#> [1] -0.2982402 10.1211065
```

Using Myr as time unit, the distribution of times for the tie points is as follows:

Mathematically, sedimentation rates are assumed to be stochastic processes. With each iteration of the estimation procedure, a sample path is generated from the stochastic processes. This sample path reflects one possible change of sedimentation rate in the section, given our uncertainties about it.

Computationally, this is implemented using *function
factories*, which are functions that return functions. A function
factory defines a stochastic process, and each function \(f\) generated by a function factory is a
sample path. In turn, \(f(x)\) returns
the value of the sample path at \(x\)
.

Here, function factories are used as complex random number generators: Instead of returning one or multiple random numbers, they return a random function.

Available wrappers to define sedimentation rates are

`sed_rate_gen_from_bounds`

: generate sed. rate from upper and lower bounds on the sedimentation rate (see above)`sed_rate_from_matrix`

: specify sedimentation rate based on matrix, to be used in conjunction with`get_data_from_eTimeOpt`

Sedimentation rates must be coded as function factories, i.e., functions that return functions. They must be able to take vector inputs and return a vector of the same length as output, and always return strictly positive values.

As example, I use a simple sedimentation rate model, where only upper and lower bounds on sedimentation rates in the section are known. Between these limits, I assume a uniform distribution.

```
h_min = 10
h_max = 90
# limits on sed. rates
lower_limit = c(0.1,2,0.1,10)
upper_limit = c(0.2,3,2,12)
# strat intervals where sed rates are defined
s = c(h_min, 30,65, 80, h_max)
```

Based on these parameters, the sedimentation rate function factory is defined as follows:

```
# define function factory
sed_rate_fun = function(){
# draw sed rates from uniform distribution
aa = runif(n = length(lower_limit), min = lower_limit, max = upper_limit)
# define sed rate "realization" based on samples from uniform distribution
sed_rate_fun = approxfun(x = s,
y = c(aa, aa[length(aa)]),
method = "constant",
rule = 2,
f = 1)
return(sed_rate_fun)
}
```

Note that the inner function is a function of one variable (height),
while the outer function takes no arguments - it simply returns the
inner function. To visualize this, let’s plot three sedimentation rates
generated by the “sedimentation rate function factory”
`sed_rate_fun`

:

```
plot(NULL,
xlim = c(h_min, h_max),
ylim = c(0, max(upper_limit)),
xlab = "Stratigraphic Height [m]",
ylab = "Sedimentation Rate")
no_of_sedrates = 3 # no. of sed rates displayed
h = seq(h_min,h_max, by = 0.1) # strat. positions where sed rates are plotted
cols = c("red", "blue", "black")
for (i in seq_len(no_of_sedrates)){
# generate sed rate from the factory
sed_rate_sample = sed_rate_fun()
# plot sed rate in the section
lines(h, sed_rate_sample(h), col = cols[i])
}
```

All sedimentation rates generated by `sed_rate_fun`

will
be different, because they are determined by random numbers.

For information on estimating age-depth models from tracer contents of rocks and sediments, see

For an overview of the structure of the `admtools`

package
and the classes used therein see

For an overview over all available vignettes for the
`admtools`

package use