diff --git a/.Rbuildignore b/.Rbuildignore
index 850ff3a..a8b5434 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -7,3 +7,5 @@
 ^pkgdown$
 
 ^CONTRIBUTING.md$
+^doc$
+^Meta$
diff --git a/.gitignore b/.gitignore
index a992cae..b06d6a0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,5 @@
 .Rproj.user
 inst/doc
 docs
+/doc/
+/Meta/
diff --git a/NAMESPACE b/NAMESPACE
index d115956..fe1ebb4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(apply_niche_pref)
+export(apply_taphonomy)
 export(cnd)
 export(ornstein_uhlenbeck)
 export(p3)
diff --git a/R/apply_niche_pref.R b/R/apply_niche_pref.R
index 7868c8c..d0a2eab 100644
--- a/R/apply_niche_pref.R
+++ b/R/apply_niche_pref.R
@@ -38,7 +38,11 @@ apply_niche_pref = function(x, niche_def, gc){
     #'
     #' }
     #'
-  niche_with_time = function(t) niche_def(gc(t))
-  r = thin(x, niche_with_time)
+  #' @seealso [apply_taphonomy()] to model taphonomic effects based on the sampe principle, [thin()] for the underlying mathematical procedure
+
+  # function that returns niche preference as a function of y (typically time)
+  change_in_niche = function(y) niche_def(gc(y))
+  # thin events based on niche preference
+  r = thin(x, change_in_niche)
   return(r)
   }
diff --git a/R/apply_taphonomy.R b/R/apply_taphonomy.R
new file mode 100644
index 0000000..3d500df
--- /dev/null
+++ b/R/apply_taphonomy.R
@@ -0,0 +1,22 @@
+apply_taphonomy = function(x, pres_potential, ctc){
+  #' @export
+  #'
+  #' @title model taphonomy
+  #'
+  #' @param x numeric vector. event type data, e.g. fossil occurrences
+  #' @param pres_potential function. takes taphonomic conditions as input and returns the preservation potential (a number between 0 and 1)
+  #' @param ctc function, change in taphonomic conditions (ctc) with time.
+  #'
+  #' @description
+    #' models taphonomy by combining the change in taphonomic conditions with the preservation potential as a function of taphonomic conditions to determine how preservation potential changes. This is then used to systematically remove (thin) the event data using `thin`
+    #'
+    #'
+  #' @seealso [apply_niche_pref()] for modeling niche preferences based on the same principle, [thin()] for the underlying mathematical procedure
+  #'
+
+  # function that returns preservation potential as a function of input (e.g. time or position)
+  change_pres_pot = function(y) pres_potential(ctc(y))
+  # thin events based on preservation potential
+  r = thin(x, change_pres_pot)
+  return(r)
+}
diff --git a/R/ornstein_uhlenbeck.R b/R/ornstein_uhlenbeck.R
index 7442a6d..1d40c9b 100644
--- a/R/ornstein_uhlenbeck.R
+++ b/R/ornstein_uhlenbeck.R
@@ -10,7 +10,7 @@ ornstein_uhlenbeck = function(t, mu = 0, theta = 1, sigma = 1, y0 = 0){
   #' @param y0 scalar. initial value (value of process at the first entry of t)
   #'
   #' @description
-    #' simulates ornstein-uhlenbeck using the Euler-Maruyame method. The process is simulated on a scale of 0.25 * min(diff(t)) and then interpolated to the values of `t`.
+    #' simulates the Ornstein-Uhlenbeck process using the Euler-Maruyame method. The process is simulated on a scale of 0.25 * min(diff(t)) and then interpolated to the values of `t`.
     #'
   #' @returns a list with two elements: `t` and `y`. `t` is a duplicate of the input `t`, `y` are the values of the OU process at these times. Outputs are of class `timelist` and can thus be plotted directly using `plot`, see `?plot.timelist`
   #'
diff --git a/R/stasis.R b/R/stasis.R
index b4da51f..36ca1cc 100644
--- a/R/stasis.R
+++ b/R/stasis.R
@@ -8,7 +8,7 @@ stasis = function(t, mean = 0, sd = 1){
   #' @param sd standard deviation
   #'
   #' @description
-    #' simulates stasis as independent, normally distributed random variables with mean `mean` and standard deviatin `sd`
+    #' simulates stasis as independent, normally distributed random variables with mean `mean` and standard deviation `sd`
     #'
   #' @examples
     #' \dontrun{
diff --git a/R/thin.R b/R/thin.R
index 4f2cb4c..165c494 100644
--- a/R/thin.R
+++ b/R/thin.R
@@ -11,6 +11,8 @@ thin = function(x, thin){
   #' @description
     #' Thins a vector of events using the function thin, meaning the probability that the ith event in x is preserved is given by _thin (x(i))_. Values of
     #' `thin` below 0 and above 1 are ignored
+    #' Is used to model niche preferences in apply_niche_pref, where events are thinned based on how close they are to their preferred niche
+    #' Can also be used to model taphonomic effects. In this case, `thin` describes the preservation potential.
     #'
   #' @examples
     #' x = p3(rate = 100, from = 0, to = 3 * pi) # simulate Poisson point process
@@ -19,10 +21,11 @@ thin = function(x, thin){
     #' yy = thin(x, function(x) 5 * sin(x))
     #' hist(yy) # note how values of 5 * sin above 1 are not affecting the thinning
     #'
+  #' @seealso [apply_niche_pref()] and [apply_taphonomy()] for use cases with biological meaning
 
   p = pmax(pmin(thin(x), rep(1, length(x))), rep(0, length(x))) # cut off at 0 and 1
-  ind = stats::rbinom(n = length(x), size = 1, prob = p)
-  y = x[as.logical(ind)]
+  ind = stats::rbinom(n = length(x), size = 1, prob = p) # determine if preserved or not
+  y = x[as.logical(ind)] # select preserved events
   return(y)
 
 }
diff --git a/man/apply_niche_pref.Rd b/man/apply_niche_pref.Rd
index 9e86590..3f74b66 100644
--- a/man/apply_niche_pref.Rd
+++ b/man/apply_niche_pref.Rd
@@ -45,3 +45,6 @@ Combines the functions niche_def and gc ("gradient change") to determine how the
 }
 
 }
+\seealso{
+\code{\link[=apply_taphonomy]{apply_taphonomy()}} to model taphonomic effects based on the sampe principle, \code{\link[=thin]{thin()}} for the underlying mathematical procedure
+}
diff --git a/man/apply_taphonomy.Rd b/man/apply_taphonomy.Rd
new file mode 100644
index 0000000..c2a63ec
--- /dev/null
+++ b/man/apply_taphonomy.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/apply_taphonomy.R
+\name{apply_taphonomy}
+\alias{apply_taphonomy}
+\title{model taphonomy}
+\usage{
+apply_taphonomy(x, pres_potential, ctc)
+}
+\arguments{
+\item{x}{numeric vector. event type data, e.g. fossil occurrences}
+
+\item{pres_potential}{function. takes taphonomic conditions as input and returns the preservation potential (a number between 0 and 1)}
+
+\item{ctc}{function, change in taphonomic conditions (ctc) with time.}
+}
+\description{
+models taphonomy by combining the change in taphonomic conditions with the preservation potential as a function of taphonomic conditions to determine how preservation potential changes. This is then used to systematically remove (thin) the event data using \code{thin}
+}
+\seealso{
+\code{\link[=apply_niche_pref]{apply_niche_pref()}} for modeling niche preferences based on the same principle, \code{\link[=thin]{thin()}} for the underlying mathematical procedure
+}
diff --git a/man/ornstein_uhlenbeck.Rd b/man/ornstein_uhlenbeck.Rd
index 18202bf..0fb21d6 100644
--- a/man/ornstein_uhlenbeck.Rd
+++ b/man/ornstein_uhlenbeck.Rd
@@ -21,7 +21,7 @@ ornstein_uhlenbeck(t, mu = 0, theta = 1, sigma = 1, y0 = 0)
 a list with two elements: \code{t} and \code{y}. \code{t} is a duplicate of the input \code{t}, \code{y} are the values of the OU process at these times. Outputs are of class \code{timelist} and can thus be plotted directly using \code{plot}, see \code{?plot.timelist}
 }
 \description{
-simulates ornstein-uhlenbeck using the Euler-Maruyame method. The process is simulated on a scale of 0.25 * min(diff(t)) and then interpolated to the values of \code{t}.
+simulates the Ornstein-Uhlenbeck process using the Euler-Maruyame method. The process is simulated on a scale of 0.25 * min(diff(t)) and then interpolated to the values of \code{t}.
 }
 \examples{
 \dontrun{
diff --git a/man/stasis.Rd b/man/stasis.Rd
index 4120723..ff0d74d 100644
--- a/man/stasis.Rd
+++ b/man/stasis.Rd
@@ -14,7 +14,7 @@ stasis(t, mean = 0, sd = 1)
 \item{sd}{standard deviation}
 }
 \description{
-simulates stasis as independent, normally distributed random variables with mean \code{mean} and standard deviatin \code{sd}
+simulates stasis as independent, normally distributed random variables with mean \code{mean} and standard deviation \code{sd}
 }
 \examples{
 \dontrun{
diff --git a/man/thin.Rd b/man/thin.Rd
index e194486..060232b 100644
--- a/man/thin.Rd
+++ b/man/thin.Rd
@@ -14,6 +14,8 @@ thin(x, thin)
 \description{
 Thins a vector of events using the function thin, meaning the probability that the ith event in x is preserved is given by \emph{thin (x(i))}. Values of
 \code{thin} below 0 and above 1 are ignored
+Is used to model niche preferences in apply_niche_pref, where events are thinned based on how close they are to their preferred niche
+Can also be used to model taphonomic effects. In this case, \code{thin} describes the preservation potential.
 }
 \examples{
 x = p3(rate = 100, from = 0, to = 3 * pi) # simulate Poisson point process
@@ -23,3 +25,6 @@ yy = thin(x, function(x) 5 * sin(x))
 hist(yy) # note how values of 5 * sin above 1 are not affecting the thinning
 
 }
+\seealso{
+\code{\link[=apply_niche_pref]{apply_niche_pref()}} and \code{\link[=apply_taphonomy]{apply_taphonomy()}} for use cases with biological meaning
+}
diff --git a/vignettes/advanced_functionality.Rmd b/vignettes/advanced_functionality.Rmd
index 18c5897..fc50b3e 100644
--- a/vignettes/advanced_functionality.Rmd
+++ b/vignettes/advanced_functionality.Rmd
@@ -18,4 +18,65 @@ knitr::opts_chunk$set(
 library(StratPal)
 ```
 
-TBA
+Here we go through how the presented workflows can be generalized and expanded upon.
+
+## Different forward models
+
+In the examples, we used age-depth models from a carbonate platform simulated with CarboCAT lite. This can easily be expanded to any other sedimentary forward model (be it marine or terrestrial, siliciclastic, mixed, or carbonate system).
+
+For this, age-depth information needs to be extracted from the forward model. The vectors of elapsed model time vs. accumulated height can then be handed over to `tp_to_adm` to define age-depth model objects. These can then be used in the pipelines as in the examples.
+
+If the forward model included erosion, you need to first account for time intervals where sediment is first deposited and later removed. For this, first define a sediment accumulation curve using `tp_to_sac`, and pass it to `sac_to_adm` to turn it into an age-depth model. You can then use the resulting age-depth model for your analysis pipeline.
+
+## Empirical age-depth models
+
+Of course you can use empirical age-depth models in your pipeline. You can create them from tie points using `tp_to_adm` just as in the case with forward models.
+
+The workflow in `StratPal` does not work with `multiadm` objects as defined by `admtools`. These objects store multiple age-depth models to account for uncertainty. You can however reduce their complexity, e.g. via `admtools::median_adm`, and use the resulting age-depth model for the modeling pipeline.
+
+## Other ecological gradients and niche definitions
+
+In the examples we used water depth as ecological gradient. This was because water depth was stored by the forward model used. The approach for niche modeling used in `StratPal` works for all types of niches. Potentially relevant niches where information can be extracted from forward models are e.g. substrate consistency, temperature, water energy, and many others.
+
+We used a simple niche definition based on a probability density function of a normal distribution. This can be expanded to arbitrary niche definitions along a gradient.
+
+For niche modeling, two components are needed:
+
+1.  A function `niche_def` . It takes as input values of the gradient, and returns numbers returns numbers between 0 and 1 as output (where 1 corresponds to the taxon being in its preferred niche, and 0 to it being completely out of its preferred niche). It must take vectorized input, and return a vector of equal length
+2.  A function `gc` that describes how the gradient changes with time. It must take vectorized inputs, and return a vector of equal length.
+
+As long as these conditions are met, niche preferences along any gradient can be modeled using `apply_niche_pref`.
+
+## Different models of phenotypic evolution
+
+Other types of phenotypic evolution in the time domain can easily be incorporated. They should stick to the following conventions to be seamlessly integrated:
+
+-   They should take a vector of times as first argument
+
+-   The implementation should be for a continuous time version of the mode of interest, as the time domain is generally sampled at irregular times (see Hohmann et al. 2024, methods section, for details)
+
+-   They should return a list with two elements: One named `t` that is a duplicate of the first argument handed to the function. The second should be named `y` and should contain the simulated trait values.
+
+-   The output list `l` should be assigned both the class "list" and "timelist" via the command `class(l) = c("timelist", "list")` , see `random_walk` as example
+
+The last two points make sure plotting via `plot.timelist` and `plot.stratlist` works seamlessly.
+
+## Different models of event type data
+
+We used (constant and variable rate) Poisson point processes to simulate event type data. Any model of event-type data can be used (e.g. with temporal correlation), as long as it returns a vector with the timing of the events
+
+## Taphonomy for event type data
+
+Taphonomic effects can be incorporated using using the `apply_taphonomy` function. This is based on the same principle as niche modeling, and makes use of the `thin` function, but was not shown in the examples.
+
+## Transform different types of data
+
+`strat_to_time` and `time_to_strat` are generic functions of the `admtools` package that can transform different types of data. Currently they transform phylogenetic trees (`phylo` objects), lists (`list` class), and numeric data such as the event type data (class `numeric`).
+
+In general, any type temporal (resp. stratigraphic) data can be transformed. For this, you need define an `S3` class for the object you would like to transform, and then implement the transformation of this object for `time_to_strat` and `strat_to_time`.
+
+If you would like to add transformations for new S3 classes to the `admtools` package, please see the contribution guidelines of the `admtools` package for details on how to contribute code.
+
+## References
+
+\* Niklas Hohmann, Joël R Koelewijn, Peter Burgess, Emilia Jarochowska. Identification of the mode of evolution in incomplete carbonate successions.bioRxiv 2023.12.18.572098; doi: [10.1101/2023.12.18.572098]([https://doi.org/10.1101/2023.12.18.572098),](https://doi.org/10.1101/2023.12.18.572098),) Data published under a [CC-BY 4.0]([https://creativecommons.org/licenses/by/4.0/)](https://creativecommons.org/licenses/by/4.0/)) license.
diff --git a/vignettes/event_data.Rmd b/vignettes/event_data.Rmd
index 67af56b..4a7a3b0 100644
--- a/vignettes/event_data.Rmd
+++ b/vignettes/event_data.Rmd
@@ -14,43 +14,48 @@ knitr::opts_chunk$set(
 )
 ```
 
+
+## Introduction
+
+This vignette provides an introduction to stratigraphic paleobiology applied to "_event type data_" such as location and age of specimens, and first/last occurrences of taxa.
+
+First, let's load the required packages:
+
 ```{r setup}
 library(StratPal)
 library(admtools)
 ```
 
-
-## Introduction
-
-This vignette provides an introduction to stratigraphic paleobiology applied to fossil occurrences and first/last taxon occurrences.
-
 ## Modeling fossil occurrences and first/last fossil occurrences
 
-Fossil occurrences and first/last occurrences can be considered distinct events that occur at specific points in time or stratigraphic position. We model them as point processes. In the `StratPal` package, we provide two functions to simulate point processes:
+Location and age of fossil specimens and first/last occurrences of taxa can be considered distinct events that occur at specific points in time or stratigraphic positions. Up to a few notable exceptions, the same principles apply to this "_event type data_", so we treat them together.
+
+We model events as _point processes_. In the `StratPal` package, we provide two functions to simulate point processes:
 
 * `p3`, short for Poisson point process, to simulate events that occur at a constant rate
 
 * `p3_var_rate` for variable rate Poisson point processes to simulate events that occur at variable rates (e.g. extinction events etc.)
 
-The usage of `p3` is straightforward:
+The usage of `p3` is straightforward.
 
 ```{r}
 # simulate fossil occurrences over one Myr with an average of 8 occurrences per Myr
-p3(rate = 8, from = 0, to = 1) |>
+p3(rate = 15, from = 0, to = 1) |>
   hist(main = "Fossil occurrences",
-       xlab = "Time [Myr]")
+       xlab = "Time [Myr]",
+       ylab = "# Fossils")
 ```
 
-For `p3_var_rate`, you can pass either a function that specifies the rate:
+For `p3_var_rate`, you can pass either a function that specifies the rate to `x`:
 
 ```{r}
 # return 100 occurrences by setting n parameter
-p3_var_rate(sin, from = 0, to = 9, n = 100) |>
+p3_var_rate(x = sin, from = 0, to = 9, n = 100) |>
   hist(xlab = "Time [Myr]", main = "Fossil occurrences")
 # note that negative rates (where sin < 0) are ignored
 ```
 
-Alternatively, you can pass it two vectors `x` and `y`. The rate is then determined by linear interpolation between these vectors. This is equivalent to using the function `approxfun(x,y)`
+Alternatively, you can pass it two vectors `x` and `y`. The rate is then determined by linear interpolation between these vectors. This is equivalent to using the function `approxfun(x,y, rule = 2)` as input for `x`
 
 ```{r}
 # decline in fossil occurrences from 50 to 0 over 1 Myr
@@ -58,69 +63,76 @@ p3_var_rate(x = c(0,1), y = c(50, 0), from = 0, to = 1, f_max = 50) |>
   hist(xlab = "Time [Myr]", main = "Fossil occurrences")
 ```
 
-Note that it is only a matter of interpretation whether the "events" produced by `p3` and `p3_var_rate` are fossil occurrences or first/last occurrences. The conceptual framework holds nonetheless.
+Keep in mind that it is only a matter of interpretation whether the "events" produced by `p3` and `p3_var_rate` are age/timing or locations of specimens from one taxon or of first/last occurrences of multiple taxa. The conceptual framework holds for all event type data.
+
+**Task:** Model an "extinction event" with a background rate of 3 species disappearing per Myr, and a peak in extinction rate where the rate is 10-fold increased over an interval of 200 kyr. How distinct is the peak in the histogram, and how much do the results differe due to randomness when you run the analysis multiple times?
 
 ## Age-depth models
 
-Throughout the examples, we use the same two age-depth model, which is from 2 km and 12 km offshore in scenario A. We define them
+Throughout this examples, we the age-depth models from 2 km and 12 km offshore in scenario A.
 
 ```{r, figures-side, fig.show="hold", out.width="50%", fig.align='left'}
-# define ADM 2 km from shore
-adm_2km = tp_to_adm(t = scenarioA$t_myr,
+adm_2km = tp_to_adm(t = scenarioA$t_myr,   # 2 km from shore
                 h = scenarioA$h_m[,"2km"],
                 T_unit = "Myr",
                 L_unit = "m")
 
-adm_12km = tp_to_adm(t = scenarioA$t_myr,
+adm_12km = tp_to_adm(t = scenarioA$t_myr,  # 12 km from shore
                 h = scenarioA$h_m[,"12km"],
                 T_unit = "Myr",
                 L_unit = "m")
 
-# plot age-depth model
-plot(adm_2km,
-     lwd_acc = 2,   # plot thicker lines for accumulative intervals (lwd = line width)
-     lty_destr = 0) # don't plot destructive intervals (lty = line type)
-T_axis_lab() # add time axis label
-L_axis_lab() # add length axis label
-plot(adm_12km,
-     lwd_acc = 2,   # plot thicker lines for accumulative intervals (lwd = line width)
-     lty_destr = 0) # don't plot destructive intervals (lty = line type)
-T_axis_lab() # add time axis label
-L_axis_lab() # add length axis label
+plot(adm_2km,       # plot age-depth model 2 km from shore
+     lwd_acc = 2,   # plot thicker lines for intervals with sediment accumulation (lwd = line width)
+     lty_destr = 0) # don't plot destructive intervals/gaps (lty = line type)
+T_axis_lab()        # add time axis label
+L_axis_lab()        # add length axis label
+title("Age-depth model 2 km from shore")
+plot(adm_12km,      # plot age-depth model 12 km from shore
+     lwd_acc = 2,   
+     lty_destr = 0)
+T_axis_lab() 
+L_axis_lab() 
+title("Age-depth model 12 km from shore")
 ```
 
+**Task:** Plot the sea level curve for scenarion A, and the water depth curves at the two locations. How is this connected to the age-depth models, specifically to the presence and duration of gaps and low/high sedimentation rates? You can use `plot_sed_rate_t`, `get_hiat_list` and `get_hiat_duration` from `admtools`.
 
 ## Stratigraphic distortions of "event" type data
 
-### Fossil occurrences
+Any event type of data can directly be transformed using an age-depth model via `time_to_strat`. Depending on how the events are interpreted, we need to choose different options. 
 
-Any "event" type data can directly transformed using `time_to_strat`:
+### Fossil specimens
+
+Assuming the events are fossil specimens of one taxon, we can examine where specimens appear in the stratigraphic column:
 
 ```{r}
 p3(rate = 200, from = min_time(adm_2km), to = max_time(adm_2km)) |>
   time_to_strat(adm_2km, destructive = TRUE) |>
   hist(xlab = "Stratigraphic height [m]",
-       main = "Fossil occurrences 2 km offshore",
+       main = "Fossil abundance 2 km offshore",
+       ylab = "# Fossils",
        breaks = seq(from = min_height(adm_2km), to = max_height(adm_2km), length.out = 20))
 ```
 
-Note the option `destructive = TRUE` in `time_to_strat`. It makes sure any fossil occurrences that coincide with a hiatus are destroyed, and do not appear in the stratigraphic domain.
+Note the option `destructive = TRUE` in `time_to_strat`. This is because we assume specimens that lived during a hiatus are destroyed.
 
-The pattern looks very different 2 km offshore:
+Assuming the same rate of fossil specimens, we get a very different pattern 12 km from shore:
 
 ```{r}
 p3(rate = 200, from = min_time(adm_12km), to = max_time(adm_12km)) |>
   time_to_strat(adm_12km, destructive = TRUE) |>
   hist(xlab = "Stratigraphic height [m]",
-       main = "Fossil occurrences 12 km offshore",
+       main = "Fossil abundance 12 km offshore",
+       ylab = "# Fossils",
        breaks = seq(from = min_height(adm_12km), to = max_height(adm_12km), length.out = 20))
 ```
 
-Not only is the section much shorter (not even 20 m), but ther are 2 distinct peaks at the bottom of the section and around 13 m. Comparing with the age-depth model shows that these locations correspond to times of extreme stratigraphic condensation (low sedimentation rates, indicated by a near-horizontal age-depth model). As a result, more time is represented per stratigraphic increment, leading to accumulations of fossil occurrences.
+Not only is the section much shorter (not even 20 m), but there are 2 distinct peaks at the bottom of the section and around 13 m. Comparing with the age-depth model shows that these locations correspond to times of extreme stratigraphic condensation (low sedimentation rates, indicated by a near-horizontal age-depth model). As a result, more time is represented per stratigraphic increment, leading to accumulations of fossil occurrences (Hohmann 2021).
 
-### Last occurrences of common taxa
+### First/Last occurrences of common taxa
 
-Last (and first) occurrences can also directly be transformed using `time_to_strat` under the condition that the taxa are very common. This assumption will be explored further below.
+In the fossil record, we can not directly observe origination and extinction rates of taxa, but only their first and last occurrences. This type of event data can also directly be transformed using `time_to_strat` when taxa are very common. This assumption will be explored further below.
 
 ```{r}
 p3(rate = 200, from = min_time(adm_2km), to = max_time(adm_2km)) |>
@@ -130,7 +142,7 @@ p3(rate = 200, from = min_time(adm_2km), to = max_time(adm_2km)) |>
        breaks = seq(from = min_height(adm_2km), to = max_height(adm_2km), length.out = 20))
 ```
 
-This would be the pattern of last occurrences in a section 2 km offshore when the extinction rate is constant! The peaks coincide with the hiatuses, as the last occurrences of the species that disappear during the hiatus cluster at the hiatus surface. Note here that we used the option `destructive = FALSE` for `time_to_strat` to make sure the last occurrences appear at the location of the hiatus. This is why the assumption of abundant taxa is important: If taxa are rare, their last occurrence does not match their actual time of extinction (Signor-Lipps effect), and the last occurrences can be found further below the hiats.
+This would be the pattern of last occurrences in the section 2 km offshore when the rate of last occurrences is constant! The peaks coincide with the hiatuses, as the last occurrences of the species that disappear during the hiatus cluster at the hiatus surface. Note here that we used the option `destructive = FALSE` for `time_to_strat` to make sure the last occurrences appear at the location of the hiatus. This is why the assumption of abundant taxa is important: If taxa are rare, their last occurrence does not match their actual time of extinction (Signor-Lipps effect), and the last occurrences can be found further below the hiatus.
 
 The patter 12 km offshore is very different:
 
@@ -258,3 +270,25 @@ list("t" = t, "y" = wd) |>
   time_to_strat(adm_2km) |> 
   plot(orientation = "lr", type = "l", xlab = "Stratigraphic position [m]", ylab = "Water depth [m]")
 ```
+
+## Further reading
+
+Go to
+
+```{r, eval=FALSE}
+vignette("phenotypic_evolution")
+```
+
+for details on how to model stratigraphic paleobiology of phenotypic evolution, or explore the vignette online under ([mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution](https://mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution.html)).
+
+See also
+
+```{r, eval=FALSE}
+vignette("advanced_functionality")
+```
+
+for details on how to expand on the modeling pipelines described here, or explore the vignette online under ([mindthegap-erc.github.io/StratPal/articles/advanced_functionality](https://mindthegap-erc.github.io/StratPal/articles/advanced_functionality.html)).
+
+## Reference
+
+* NIKLAS HOHMANN; INCORPORATING INFORMATION ON VARYING SEDIMENTATION RATES INTO PALEONTOLOGICAL ANALYSES. PALAIOS 2021;; 36 (2): 53–67. doi: https://doi.org/10.2110/palo.2020.038
diff --git a/vignettes/phenotypic_evolution.Rmd b/vignettes/phenotypic_evolution.Rmd
index f29f668..54f372b 100644
--- a/vignettes/phenotypic_evolution.Rmd
+++ b/vignettes/phenotypic_evolution.Rmd
@@ -29,11 +29,13 @@ library(admtools)
 
 The `StratPal` package provides three functions for simulating different modes of evolution in the time domain:
 
--   `stasis` simulates evolutionary stasis as independent, normally distributed random variables with mean `mean` and standard deviation `sd`.
+-   `stasis` simulates evolutionary stasis as independent, normally distributed random variables with mean `mean` and standard deviation `sd` (Hunt 2006).
 
--   `random_walk` simulates trait evolution following a (potentially biased) random walk with variability `sigma`, directionality `mu`, and initial trait value `y0`. Setting `mu` to a negative (positive) value will make the random walk decrease (increase) over time, `sigma` determines the effect of randomness on the trait values.
+-   `random_walk` simulates trait evolution following a (potentially biased) random walk with variability `sigma`, directionality `mu`, and initial trait value `y0` (e.g. Bookstein 1987). Setting `mu` to a negative (positive) value will make the value of the random walk decrease (increase) over time, `sigma` determines the effect of randomness on the trait values.
 
--   `ornstein_uhlenbeck` corresponds to convergence to a phenotypic optimum, where `mu` is the optimal value (long term mean), `theta` determines how fast `mu` is approached, `sigma` the noise level, and `y0` the initial trait value.
+-   `ornstein_uhlenbeck` corresponds to convergence to a phenotypic optimum, where `mu` is the optimal value (long term mean), `theta` determines how fast `mu` is approached, `sigma` the noise level, and `y0` the initial trait value (Lande 1976).
+
+Use `vignette("advanced_functionality")` for guidelines to implement your own modes of evolution.
 
 ### Visualization
 
@@ -55,31 +57,35 @@ We are interested in how phenotypic evolution is preserved in the stratigraphic
 
 ### Age-depth models
 
-As an example, we compare the the preservation of trait evolution 2 km and 12 km offshore in the carbonate platform in scenarion A (see `?scenarioA`). For this, first define the age-depth models
+As an example, we compare the the preservation of trait evolution 2 km and 12 km offshore in the carbonate platform in scenario A. See `?scenarioA` and `vignette("StratPal")` for details on scenario A.
+
+First define the age-depth models:
 
 ```{r, figures-side, fig.show="hold", out.width="50%", fig.align='left'}
-# define ADM 2 km from shore
-adm_2km = tp_to_adm(t = scenarioA$t_myr,
+adm_2km = tp_to_adm(t = scenarioA$t_myr,    # 2 km from shore
                 h = scenarioA$h_m[,"2km"],
                 T_unit = "Myr",
                 L_unit = "m")
 
-adm_12km = tp_to_adm(t = scenarioA$t_myr,
+adm_12km = tp_to_adm(t = scenarioA$t_myr,   # 12 km from shore
                 h = scenarioA$h_m[,"12km"],
                 T_unit = "Myr",
                 L_unit = "m")
 
-# plot age-depth model
-plot(adm_2km,
-     lwd_acc = 2,   # plot thicker lines for accumulative intervals (lwd = line width)
-     lty_destr = 0) # don't plot destructive intervals (lty = line type)
-T_axis_lab() # add time axis label
-L_axis_lab() # add length axis label
-plot(adm_12km,
-     lwd_acc = 2,   # plot thicker lines for accumulative intervals (lwd = line width)
-     lty_destr = 0) # don't plot destructive intervals (lty = line type)
-T_axis_lab() # add time axis label
-L_axis_lab() # add length axis label
+
+plot(adm_2km,       # plot age-depth model 2 km from shore
+     lwd_acc = 2,   # plot thicker lines for intervals with sediment accumulation (lwd = line width)
+     lty_destr = 0) # don't plot destructive intervals/gaps (lty = line type)
+T_axis_lab()        # add time axis label
+L_axis_lab()        # add length axis label
+title("Age-depth model 2 km from shore")
+
+plot(adm_12km,      # plot age-depth model 12 km from shore
+     lwd_acc = 2,   
+     lty_destr = 0) 
+T_axis_lab() 
+L_axis_lab() 
+title("Age-depth model 12 km from shore")
 ```
 
 ### Stratigraphic expression of phenotypic evolution
@@ -93,12 +99,12 @@ seq(from = min_time(adm_2km), to = max_time(adm_2km), by = 0.01) |> # sample eve
   time_to_strat(adm_2km, destructive = FALSE) |>                    # transform data from time to strat domain
   plot(type = "l",                                                  # plot
        orientation = "lr",
-       xlab = "Stratigraphic height [m]",
+       xlab = paste0("Stratigraphic height [", get_L_unit(adm_2km), "]"),
        ylab = "Trait value")
 
 ```
 
-This is what a biased random walk in the time domain would look like if it was observed 2 km offshore in the simulated carbonate platform. The large jumps in traits correspond to long hiatuses caused by prolonged drops in relative sea level. Compare this figure to the random walk in the time domain (without stratigraphic distortions.)
+This is what a biased random walk in the time domain would look like if it was observed 2 km offshore in the simulated carbonate platform. The large jumps in traits correspond to long hiatuses caused by prolonged drops in relative sea level (Hohmann et al. 2024). Compare this figure to the random walk in the time domain (without stratigraphic distortions).
 
 12 km offshore, the preservation is very different:
 
@@ -108,11 +114,11 @@ seq(from = min_time(adm_12km), to = max_time(adm_12km), by = 0.01) |> # sample e
   time_to_strat(adm_12km, destructive = FALSE) |>                     # transform data from time to strat domain
   plot(type = "l",                                                    # plot results
        orientation = "lr",
-       xlab = "Stratigraphic height [m]",
+       xlab = paste0("Stratigraphic height [", get_L_unit(adm_12km), "]"),
        ylab = "Trait value")
 ```
 
-You can also see two jumps, but the first is at a different location compared to the section 2 km offshore. Looking at the age-depth models, we can see why this is: 2 km offshore, the jumps are caused by prolonged hiatuses that are associated with the massive drops in sea level around 0.5 Myr and 1.5 Myr. 12 km offshore, the age-depth model shows fewer gaps, but prolonged intervals with low sedimentation rates at the beginning of the simulation and at around 1.5 Myr. This leads to stratigraphically condensed intervals, at which phenotypic evolution appears to be accelerated. You can see that not all artefactual jumps in traits observable in the stratigraphic record are caused by gaps.
+You can also see two jumps, but the first is at a different location compared to the section 2 km offshore. Looking at the age-depth models, we can see why this is: 2 km offshore, the jumps are caused by prolonged hiatuses that are associated with the massive drops in sea level around 0.5 Myr and 1.5 Myr. 12 km offshore, the age-depth model shows fewer gaps, but prolonged intervals with low sedimentation rates at the beginning of the simulation and at around 1.5 Myr. This leads to intervals where the rate phenotypic evolution appears to be accelerated because of stratigraphic condensation (Hohmann 2021). You can see that not all artefactual jumps in traits observable in the stratigraphic record are caused by gaps.
 
 **Task:** How does this effect vary between different modes of evolution (with different parameter choices), and at different locations in the platform (e.g., along an onshore-offshore gradient)? Can you make any general statements about where you see the strongest effects?
 
@@ -145,22 +151,22 @@ sampling_loc_m |>                                  # sampling locations
   plot(orientation = "lr",                         # plot stratigrahic data
        type = "l",        
        ylab = "Trait Value",
-       xlab = "Stratigraphic Height [m]")
+       xlab = paste0("Stratigraphic height [", get_L_unit(adm_2km), "]"))
 ```
 
-There we have it! This is what the lineage would look like if it would be sampled every 2 meters in a sections 2 km offshore.
+There we have it! This is what the lineage would look like if it was sampled every 2 meters in a sections 2 km from shore.
 
 **Task:** How does prescribing a sampling strategy change the interpretations of your last task? Can you draw any conclusions about which sampling strategy is best suited for a given environment or mode of evolution?
 
 ## Further reading
 
-Go to
+Run
 
 ```{r, eval=FALSE}
-vignette("phenotypic_evolution")
+vignette("event_data")
 ```
 
-for details on how to model stratigraphic paleobiology of phenotypic evolution, or explore the vignette online on the package webpage ([mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution](https://mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution.html)).
+for details on how to model stratigraphic paleobiology of event data such as individual fossils and first/last occurrences of taxa, or explore the vignette online under [mindthegap-erc.github.io/StratPal/articles/event_data](https://mindthegap-erc.github.io/StratPal/articles/event_data.html).
 
 See also
 
@@ -168,6 +174,19 @@ See also
 vignette("advanced_functionality")
 ```
 
-for details on how to expand on the modeling pipelines described here, or explore the vignette online on the package webpage ([mindthegap-erc.github.io/StratPal/articles/advanced_functionality](https://mindthegap-erc.github.io/StratPal/articles/advanced_functionality.html)).
+for details on how to expand on the modeling pipelines described here, or explore the vignette online under [mindthegap-erc.github.io/StratPal/articles/advanced_functionality](https://mindthegap-erc.github.io/StratPal/articles/advanced_functionality.html).
+
+## References
+
+* Bookstein, Fred L. 1987. "Random walk and the existence of evolutionary rates." Paleobiology. https://doi.org/10.1017/S0094837300009039.
+
+* Hohmann, Niklas. 2021. "Incorporating information on varying sedimentation rates into paleontological analyses." PALAIOS. https://doi.org/10.2110/palo.2020.038.
+
+* Hohmann, Niklas; Koelewijn, Joël R.; Burgess, Peter; Jarochowska, Emilia. 2024. "Identification of the mode of evolution in incomplete carbonate successions". BMC Ecology and Evolution, In Press. https://doi.org/10.1101/2023.12.18.572098.
+
+* Hunt, Gene. 2006. “Fitting and Comparing Models of Phyletic Evolution: Random Walks and Beyond.” Paleobiology. https://doi.org/10.1666/05070.1.
+
+* Lande, Russell. 1976. "Natural selection and random gentic drift in phenotypic evolution." Evolution. https://doi.org/10.1111/j.1558-5646.1976.tb00911.x
+