\[\newcommand{unit}[1]{{\rm #1}}\]

CarboKitten.jl
an open source toolkit for carbonate stratigraphic modeling

Johan Hidding1, Emilia Jarochowska2, Niklas Hohmann2, Xianyi Liu2, Peter Burgess3, Hanno Spreeuw1

  1. Netherlands eScience Center, Science Park 402 (Matrix THREE), 1098 XH Amsterdam, the Netherlands
  2. Utrecht University, Faculty of Geosciences, Princetonlaan 8a, 3584 CB Utrecht, The Netherlands
  3. University of Liverpool, School of Environmental Sciences, 4 Brownlow Street, Liverpool, L69 3GP, United Kingdom

This paper was originally published in Geoscientific Model Development GMD, 19, 6207–6229, 2026. The present version contains all the source code we used to generate the presented results. These scripts were run with CarboKitten version 0.6.1. This site is best enjoyed in a light color theme.

Abstract

Stratigraphic forward modeling is a powerful tool for testing hypotheses about the geological record and conducting numerical experiments in stratigraphy at timescales not accessible to direct human observation. Open Source software for stratigraphic modeling available so far has previously focused on siliciclastic or terrestrial depositional environments. CarboKitten is a stratigraphic forward modeling toolkit for carbonate platforms. With performance and accessibility in mind, CarboKitten is implemented in Julia, using the literate programming approach.

CarboKitten integrates three components: the carbonate production model of Bosscher and Schlager (1992), the cellular automaton for spatial heterogeneity introduced by Burgess (2013), and a novel finite difference transport model inspired by Paola et al. (1992). The model simulates carbonate production through multiple biological factories (typically euphotic, oligophotic and aphotic), accounts for ecological processes that create spatial facies patterns through cellular automaton rules, and implements sediment transport via an active layer approach where material moves along paths of steepest descent.

Key features include support for different boundary conditions, variable sea level and insolation inputs, wave-induced transport capabilities, and visualization tools aiming at beautiful plots. The software exports data in the interoperable HDF5 format and includes functions for creating stratigraphic cross-sections, chronostratigraphic (Wheeler) diagrams, topographic maps, and sediment accumulation curves. Performance benchmarks demonstrate linear scaling with grid size and time steps, enabling efficient execution on consumer hardware.

CarboKitten addresses a gap in available carbonate modeling tools by providing an accessible, well-documented, and modifiable toolkit for hypothesis testing in carbonate stratigraphy. The model operates on timescales from centuries to millions of years and can simulate various scenarios including orbital forcing, sea level change, and biological succession patterns. CarboKitten’s accessibility should encourage broader adoption of stratigraphic forward modeling in carbonate research and education, supporting hypothesis-driven approaches to understanding the structure of the geological record and reconstructing the history of the Earth from carbonate strata.

1 Introduction

Stratigraphic forward modelling is well established as a means of examining our understanding of the formation of stratal architectures (Granjeon and Joseph 1999; Burgess et al. 2001; Paterson et al. 2006; Schlager and Warrlich 2009; Ding et al. 2019; Jean Borgomano et al. 2020; Liu et al. 2022), prediction, correlation and imputation of architectures from incomplete data (Warrlich et al. 2008; Masiero et al. 2021), and testing hypotheses on the structure of the geological record (e.g., Kemp et al. 2018; Masiero et al. 2020; Liu and Liu 2021) and the preservation of proxies (Curtis et al. 2025), fossils (Holland 2000; Hannisdal 2006; Hohmann et al. 2024), or forcing mechanisms (Kemp et al. 2016; Kemp and Van Manen 2019; Burgess et al. 2019). Owing to their economic interest, most such models are proprietary to exploration companies and their availability to researchers is limited. Some older models developed by researchers share the fate of many other research software packages and their maintenance ceases, e.g. when a project ends (Warrlich 2000). It is not always possible to resuscitate such models, especially if documentation or license are lacking or code has not been shared (e.g., Strobel et al. 1989; Demicco 1998; Barrett and Webster 2017). As a result, the choice of stratigraphic forward models available to researchers at the moment is narrow and shifted towards siliciclastic (e.g., Hutton and Syvitski 2008; Sylvester et al. 2024) or fluvial depositional systems (Wild et al. 2019; Falivene et al. 2019), to the point that researchers may resort to these models to create simulations of carbonate sections (Zimmt et al. 2021).

Modeling carbonate depositional systems requires not only accounting for aqueous and atmospheric processes, but also for the biological character of sediment production and dispersal. Ecological processes, such as facilitation, competition and dispersal, may on one hand confound the relationships between sediment composition and water depth (e.g. Granjeon and Joseph 1999; Dyer et al. 2018; Salles, Pall, et al. 2018; Pastier et al. 2019; Weij et al. 2019) and, on the other hand, lead to creation of complex facies patterns under stable sea level conditions (Drummond and Dugan 1999; Purkis et al. 2016; Xi and Burgess 2022). Complex models accounting for this have been mostly developed for exploration, e.g. Carbonate 3D (Warrlich et al. 2002; Warrlich et al. 2008), DIONISOS (Granjeon and Joseph 1999) and Carbonate GPM (Hill et al. 2009). Of research-driven models operating in more than one dimension, two include a wider range of depositional environments with carbonate production modules: CARB3D+ (Paterson et al. 2006), SedSimple (Tetzlaff 2023) and Badlands (Salles 2016), including its Python interface pyBadlands (Salles, Ding, et al. 2018), but due to their general focus these models do not account for the spatial heterogeneity driven by biological processes. Finally, CarboCAT (Burgess 2013) is a research-driven 3D model dedicated to stratigraphic forward modeling of carbonate platforms, which includes a cellular automaton that approximates the spatial heterogeneity formed through ecological interactions between carbonate-producing organisms. CarboCAT has been used in multiple studies (e.g. Masiero et al. 2020; Xi and Burgess 2022; Hohmann et al. 2024), but having been written in Matlab, it is not accessible to contributions from the entire scientific community. Based on the successful applications of CarboCAT, we set out to develop a new generation model with the following specifications:

  1. it should be Open Source and it should be easy for researchers to understand the algorithm, which is a prerequisite to being able to contribute to it or modify it to one’s needs,
  2. it should allow for spatial heterogeneity of carbonate facies,
  3. it should include a sediment transport algorithm operating on different carbonate facies and produces realistic results without decreasing the model’s performance substantially,
  4. it should allow exporting and plotting multiple types of data users may need, including slices through the model grid, age-depth models, sediment accumulation curves, and stratigraphic columns,
  5. it should be performant, easy to parallelize, and platform-independent,
  6. it should be well documented and easy to use at a level accessible to a geosciences student.

The above prerequisites led us to re-designing the original architecture of CarboCAT and implementing its successor in Julia. In this article we present CarboKitten, an efficient and accessible Open Source model for stratigraphic forward simulations of carbonate platforms.

2 Model

CarboKitten combines the carbonate production model by Bosscher and Schlager (1992), the cellular automaton from Burgess (2013), and a custom finite difference transport model inspired on an approach by Paola et al. (1992). We describe each of these components in detail in the following sections.

2.1 Quantities

Since the model describes the accumulation of sediment under a range of variable conditions, a short discussion of different measures in the vertical column is in order.

Subsidence rate

Quantified as a rate \(\sigma\) in units of \(\unit{m/Myr}\). The growth of sediment is only sustainable in scenarios where there is a steady subsidence. In our models we use a default value of \(50\ \unit{m/Myr}\) (or \(0.5\ \unit{mm/kyr}\)). This parameter can be set by the users.

Initial topography

The model starts at an initial topography \(\eta_0(x) = \eta(x, t_0)\), consisting of impenetrable bedrock. A more complex topography can be provided as an input array, e.g. by running a previous model and extracting the height of sediment.

Topography

The present topography \(\eta(x, t)\) is given as the initial topography plus any amount of sediment accumulated over time. In our definition of \(\eta\) we do not correct for subsidence (see also the definition for water depth below), so it should be considered relative to a bedrock reference. This definition matches how coordinates are handled in CarboKitten internally.

Relative sea level

The relative sea level \(R(t)\) is usually a function of time, given as an input parameter of the model.

Water depth

The water depth is computed from the current topography, relative sea level and subsidence rate,

\[w(x, t) = R(t) - \eta(x, t) + \int_{t_0}^{t} \sigma \textrm{d}t.\]

2.2 Carbonate Production

The general form of our production model follows that of Bosscher and Schlager (1992). This model finds the sediment accumulation curve by integrating an ODE that outside of the model parameters only depends on the initial topography.

\[\frac{\partial \eta}{\partial t} = P(\eta),\]

where \(P\) is the sediment production in \(\textrm{m/Myr}\),

\[P(w) = \begin{cases} g_m \tanh\left(\frac{I_0}{I_k}\ e^{-kw}\right) & \text{if $w \ge 0$}\\ 0 & \text{if $w < 0$} \end{cases},\](1)

where \(I_0\) is the insolation in units of energy flux, \(I_k\) is the saturation intensity and should be provided in the same units as \(I_0\), \(k\) the extinction coefficient in \(\unit{m^{-1}}\) and \(g_m\) the maximum growth rate in \(m/Myr\).

This model encapsulates both the exponential extinction of sun light as water depth increases, and the idea that the growth of organisms interpolates between no growth at great depth and saturated growth in shallow waters (i.e. solar input is not the limiting factor at those depths).

Here we parametrize \(P\) as a function of \(w\) for \(w > 0\). Note that \(\nabla w = - \nabla \eta\), but otherwise we will use \(w\) and \(\eta\) wherever one or the other is more convenient.

Following Burgess (2013), we extend the BS92 model by introducing multiple facies that each have their own growth characteristics (except for insolation \(I_0\), which is a global input variable).

\[P(w) = \sum_f P_f(w)\]

Table 1: Parameters for the production model of the three default carbonate factories.
Factory \(g_m\) \([\unit{m/Myr}]\) \(I_k\) \([\unit{W/m^2}]\) \(k\) \([\unit{m^{-1}}]\)
Euphotic 500.0 60.0 0.8
Oligophotic 400.0 60.0 0.1
Aphotic 100.0 60.0 0.005

In the examples presented here we set default values of factory parameters that corrspond to three biological facies based on sediment produced by three carbonate factories: the euphotic (E), oligophotic (O) and aphotic (A) factories. These default values are shown in Table 1, and the resulting production curves shown in Figure 1.

We also provide the option of specifying pelagic production curves. Here the production is computed as the integral of Equation 1 over the entire water column.

Production curves

Figure 1: Production curves for our three default carbonate factories. Additionaly, we show the production curve for a pelagic facies with the same extinction coefficient and saturation intensity as the oligophotic facies, and a maximum growth rate of \(5\ \unit{Myr^{-1}}\). All of these production curves were computed with an insolation of \(400\ \unit{Wm^{-2}}.\)

runs/production-curves.jl 
#| file: "runs/production-curves.jl"
#| classes: ["task"]
#| creates: ["md/fig/production-curves.pdf"]
#| collect: figures
module Script

using CarboKitten
using CarboKitten.Production
using CairoMakie

@kwdef struct Input <: CarboKitten.AbstractInput
    insolation = 400.0u"W/m^2"
end

inch = 96
pt = 4/3
cm = inch / 2.54

function main()
    water_depth = (0.01:0.1:50.0)u"m"
    fig = Figure(size=(8.3cm, 9.0cm), fontsize=8pt)
    input = Input()

    ax = Axis(fig[1, 1], title="production at 400.0 W m⁻²",
        yreversed=true, ylabel="depth [m]", xlabel="production [m/Myr]")
    for (k, prod) in pairs(Production.EXAMPLE)
        f = production_profile(input, prod)
        p = water_depth .|> (w -> f(input.insolation, w))
        lines!(ax, p |> in_units_of(u"m/Myr"),
            water_depth |> in_units_of(u"m"), label = string(k))
    end
    axislegend(ax, "factories", valign = :bottom)
    save("md/fig/production-curves.pdf", fig)
end

end

Script.main()

2.3 Cellular Automaton

Representing spatial heterogeneity resulting from positive and negative biological interactions in a computationally simple model requires meeting the following conditions:

  1. infinite heterogeneity in space and time, i.e. no convergence on stable patterns;
  2. authigenic variation that does not require external drivers;
  3. the approach must be scalable to n facies, even if we focus the examples here on 3;
  4. adjustable temporal persistence: it must be possible to set the turnover frequency.

These requirements are met by Celullar Automata (CA), as proposed by Drummond and Dugan (1999), Burgess and Emery (2004) and Burgess (2013), without substanial computational costs. For this reason we adopt this modeling approach in CarboKitten, directly reimplementing the automaton described by Burgess (2013) in his package CarboCAT. Cellular automata are commonly used to generate spatial heterogeneity in forward modeling, with some models serving as discrete approximations of partial differential equations that can generate complex spatial dynamics such as Turing patterns Drummond and Dugan (1999).

The CA emulates the biological succession of species by following a set of simple rules. If conditions are right, a species will multiply and occupy neighbouring territory. However, when there are too many of the same kind, the species will die from over population. For each cell in the grid a centered neighbourhood of \(5\times 5\) pixels is considered. We count the number of neighbouring cells of the same species. Then we consider two ranges: the activation range (default \(6 \le n \le 10\)) and viability range (default \(4 \le n \le 10\)). If the number of live neighbours is in the viability range, the cell stays alive. If the cell was dead, but the number of live neighbours is in the activation range, the cell becomes alive. The neighborhood size and the rules represent a case of Larger than Life family of two-dimensional cellular automata (Evans 1996), but this particular set of rules has been proposed specifically for CarboCAT (Burgess 2013) and, to our knowledge, does not correspond to any documented Larger than Life rules. We examined other sets of rules (Hidding 2026), but most lead to rapid stabilization of spatial patterns. In Figure 2 we illustrate the long-term and short-term behaviour of our default CA rule set.

The initial CA grid is randomized. A dead cell may qualify to become alive for different carbonate factories at the same time. To resolve this priority collision, facies priority for occupying a dead cell is rotated every iteration using a deterministic cyclic shift (fixed round-robin pattern), which ensures that there is no priority given to any facies.

In the default configuration we emulate three species, corresponding to the factory species discussed in the section on carbonate production. The state of the CA determines which carbonate factory is switched on for each cell in the grid.

The \(5\times 5\) neighborhood strikes a balance between small-scale heterogeneity and computational cost. A smaller neighborhood would result in finer-grained spatial patterns, whereas a larger one in spatial smoothing and larger, more coherent patches of each facies. However, it would slow the model down. In a real depositional system, different carbonate producer guilds have their own length scales at which they disperse and interact and one-size-fits-all neihgborhood is clearly a simplification. The size of the neighborhood is fixed in the current version of CarboKitten, but adjustment of granularity of facies distribution generated by the CA can be achieved by users by changing the grid dimensions and the CA interval.

code 
#| classes: ["task"]
#| creates: ["md/fig/ca-first-steps.pdf"]
#| collect: figures

module Script
using CarboKitten
using CarboKitten.Components: CellularAutomaton as CA
using CairoMakie

function main()
  input = CA.Input(
      box = CarboKitten.Box{Periodic{2}}(
        grid_size=(50, 50), phys_scale=1.0u"m"),
      facies = fill(CA.Facies(), 3)
  )

  state = CA.initial_state(input)
  step! = CA.step!(input)

  fig = Figure(size=(1000, 500))
  axes_indices = Iterators.flatten(eachrow(CartesianIndices((2, 4))))
  xaxis, yaxis = box_axes(input.box)
  for (i, idx) in enumerate(axes_indices)
    ax = Axis(fig[Tuple(idx)...], aspect=AxisAspect(1), title="step $(i)")
    if idx[1] == 2
      ax.xlabel = "x"
    end
    if idx[2] == 1
      ax.ylabel = "y"
    end
    heatmap!(ax, xaxis/u"m", yaxis/u"m", state.ca)
    step!(state)
  end
  save("md/fig/ca-first-steps.pdf", fig)
end

end

Script.main()
#| classes: ["task"]
#| creates: ["md/fig/ca-long-term.png"]
#| collect: figures
module Script
using CarboKitten
using CarboKitten.Components: CellularAutomaton as CA
using CairoMakie

function main()
  input = CA.Input(
      box = CarboKitten.Box{Periodic{2}}(
        grid_size=(50, 50), phys_scale=1.0u"m"),
      facies = fill(CA.Facies(), 3)
  )

  state = CA.initial_state(input)
  step! = CA.step!(input)

  for _ in 1:1000
    step!(state)
  end

  fig = Figure(size=(1000, 500))
  axes_indices = Iterators.flatten(eachrow(CartesianIndices((2, 4))))
  xaxis, yaxis = box_axes(input.box)
  i = 1000
  for row in 1:2
    for col in 1:4
      ax = Axis(fig[row, col], aspect=AxisAspect(1), title="step $(i)")

      if row == 2
        ax.xlabel = "x"
      end
      if col == 1
        ax.ylabel = "y"
      end

      heatmap!(ax, xaxis/u"m", yaxis/u"m", state.ca)
      step!(state)
      i += 1
    end
    for _ in 1:996
      step!(state)
      i += 1
    end
  end
  save("md/fig/ca-long-term.pdf", fig)
end

end

Script.main()
CA

Figure 2: Iterations of the CA, as described by Burgess (2013), on a periodic grid of \(50\times50\). Each color represents a carbonate facies, but since the CA process is not affected by facies identity we leave out any designation. Starting with random noise, we first iterate 1000 times to get into a typical state. The top row shows iterations 1000 to 1003, the bottom row 2000 to 2003. This shows that the patterns keep reasonably stable on the short term, while evolving more extensively over the long term.

2.4 Transport

Our transport model is borrowed from other similar approaches in siliclastic (river bed) modeling (see Paola et al. 1992; James et al. 2010), where it is made plausible that this approach is viable for models that work on long time scales. Because our transport model is novel (at least for modelling carbonate platforms), we discuss the full model in a separate section. Here, we discuss how transport is embedded in the larger model.

We consider all sediment transport to happen in an active layer close to the sea floor. This layer has a certain amount of sediment \(C_f\) (in units of \(\unit{m}\)) that travels along a path of steepest descent. We say that this material is entrained. Every time step the active layer is supplied with sediment produced in the production step as well as older sediment through disintegration. After transport a fraction of the entrained sediment is deposited on the sea floor as the transported version of the original (donor) facies in a process that we refer to as lithification, being the process of turning loose sediment into rock. Although in reality sediment might not be mobile for a while before lithification sets in, for the purpose of our model, we chose the term to represent the immobilisation of sediment as a whole, see Figure 3.

Diagram showing concepts of production, lithification and disintegration

Figure 3: Diagram showing concepts of production, lithification and disintegration. Every time step newly produced sediment and older disintegrated material (configured as a disintegration rate) is added to the active layer. After transport, a set fraction of the sediment (configured as a lithification half-life time) is lithified, becoming the sea floor.

The actual transport is computed using a finite difference approach that is further discussed in Section 3.

2.5 Composed model

Diagram showing flow control of the model

Figure 4: Flow chart of model control. The shaded areas show where the internal states of the CA, sediment buffer and active layer are affected. We can follow the model, starting from the time step node. The time step implies a change in external control factors like sea level, subsidence and insolation. Together with the sediment buffer, those factors determine water depth which in turn affects both production and transport. The CA evolves by itself (there is an option for feedback from the production step, but we left it out of the diagram for simplicity), but determines which facies are produced where. We then produce and disintegrate sediment that is added to the active layer for transport. The transport step is adaptive, so it can loop several times before we lithify a fraction of the active layer. After lithification, we increment the time step, completing the loop.

Putting everything together, we evaluate the model as follows each iteration:

  1. Advance the cellular automaton.
  2. Compute the production \(P_f\).
  3. Disintegrate sediment at the rate \(r_d\).
  4. Transport entrained sediment \(C_f\).
  5. Lithify (cement) sediment with half-time of \(t_l\).

Advancing the CA can be configured to happen one-in-\(n\) iterations to slow it down. Transporting the sediment can be computed on smaller time steps if required for numeric stability. The complete loop is illustrated in Figure 4.

2.6 Input parameters

CarboKitten has many input parameters: box geometry, time parameters, a list of facies properties, transport model intrinsics and external conditions such as initial topography, relative sea level and insolation. We have already discussed the facies properties in Section 2.2, and the transport model is discussed in Section 3. That leaves us the external conditions that should be considered the driving forces of carbonate platform formation.

The initial topography, sea level and insolation can all be entered in three different ways: a given constant, a Julia function or an array exactly matching the box size or number of time steps.

In Section 6 we provide two examples where we use external sources to drive the sea level and insolation curves. A full list of input parameters is available in the Appendix.

2.7 Visualisations

CarboKitten generates data in the accessible, binary HDF5 format, thus output can be visualised with most common tools, e.g. imported into R or a Jupyter notebook. Nevertheless, we provide some routines based on Makie (Danisch and Krumbiegel 2021) for creating cross-sections, chronostratigrahic (Wheeler) diagrams and topographic overviews. Some of the most common plot types have been collected into a summary plot, which is shown in Figure 5.

Summary plot

Figure 5: Overview of different visualizations supported by CarboKitten. Panel (a) shows a stratigraphic cross-section, including an indication for unconformities, (b) a topographic overview including two intermediate time steps, (c) the production curves used, (d) sedimentation rate as a function of time, (e) Wheeler diagram with preserved dominant facies as a function of time, (f) the sea-level curve given as input. The combined plot is arranged such that spatial data is on the top row, while time-dependent information is shown at the bottom with matching y-axes.

runs/standard_example_run.jl 
#| file: runs/standard_example_run.jl
#| classes: ["task"]
#| creates: data/alcap-example.h5
module StandardExample
using CarboKitten
using CarboKitten.Models: ALCAP as M

function main()
    CarboKitten.init()
    run_model(Model{M}, M.Example.INPUT, "data/alcap-example.h5")
end
end

StandardExample.main()
#| file: runs/standard_example_plot.jl
#| classes: ["task"]
#| requires: data/alcap-example.h5
#| creates: md/fig/summary-plot.png
#| collect: figures
module StandardExamplePlot
using CairoMakie
using CarboKitten.Visualization: summary_plot

function add_panel_labels!(fig::Figure;
                          labels::Vector{String},
                          fontsize::Int64,
                          offset1::Tuple{Int64, Int64},
                          offset2::Tuple{Int64, Int64})

    panel_positions = [
        (1, 1), # a
        (1, 3), # b
        (2, 3), # c
        (3, 1), # d
        (3, 2), # e
        (3, 3)  # f
    ]

    for (i, (row, col)) in enumerate(panel_positions)
        if i <= 3
            Label(fig[row, col], labels[i],
                  tellwidth=false,
                  tellheight=false,
                  halign=:left,
                  valign=:top,
                  padding=(offset1[1], offset1[1], -offset1[2], offset1[1]),
                  fontsize=fontsize)
        else
            Label(fig[row, col], labels[i],
                    tellwidth=false,
                    tellheight=false,
                    halign=:left,
                    valign=:top,
                    padding=(offset2[1], offset2[1], -offset2[2], offset2[1]),
                    fontsize=fontsize)
        end
    end

    return fig
end


function main()
    fig = summary_plot("data/alcap-example.h5")
    add_panel_labels!(fig, labels = ["a", "b", "c", "d", "e", "f"], fontsize = 22, offset1 = (-25,-5), offset2 = (-25,10))
    save("md/fig/summary-plot.png", fig)
end
end

StandardExamplePlot.main()

3 Transport

The transport module is designed with timesteps of the order of centuries in mind and therefore it aims to lead to a similar average sediment redistribution as would have resulted from individual short-term physical processes (storms, waves, tidal currents, grain-size sorting), without attempting to resolve them. CarboKitten currently does not resolve grain-size-dependent transport, storm-driven episodic redistribution, longshore transport, or bioturbation contributing to sediment disintegration. However, such extensions are principally possible, should modelling on shorter time scales become the focus of the model’s users.

In Section 2.4 we discussed how transport is embedded in the larger model. Our transport model supposes that all entrained sediment resides in a layer of constant thickness just above the sea floor, also known as the active layer. The concentration of sediment \(C_f\) (in units of \(\unit{m}\)) is considered separately for each facies (as with all quantities with the \(f\) subscript). Each iteration of the larger model we supply the active layer with freshly produced (autochthonous) sediment as well as disintegrated older (allochthonous) sediment. We then compute transport of the active layer for as many sub-iterations as is deemed needed for the solver to remain stable. After that, a percentage of the contents in the active layer is deposited on the sea floor. The lithification percentage depends both on the time step taken and the given lithification time \(t_l\), which is configured in terms of a half-life time. There are many ways to compute sediment transport in the active layer. We have opted for a finite difference strategy inspired by Paola et al. (1992).

We assume a local sediment flux proportional to the local gradient,

\[\vec{q}_f = - C_f d_f \vec{\nabla} \eta + C_f \vec{v}_f(w),\]

where \(d_f\) is a facies-dependent transport coefficient (\(\unit{m Myr^{-1}}\)), and \(v_f(w)\) is a chosen additional velocity as a function of water depth (\(\unit{m Myr^{-1}}\)). Optionally, we use \(v_f(w)\) to model wave induced sediment transport (for an example see Section 6.3). The mass balance (continuity equation) is then,

\[\frac{\partial C_f}{\partial t} = -\vec{\nabla} \cdot \vec{q}_f\]

This gives us an advection equation for the sediment concentration \(C_f\). We also express everything in terms of water depth, having \(\nabla w = -\nabla \eta\), arriving at

\[\frac{\partial C_f}{\partial t} = -(d_f \vec{\nabla} w + \vec{v}_f(w)) \cdot \vec{\nabla}C_f - (d_f \nabla^2 w + \vec{s}_f(w) \cdot \vec{\nabla} w) C_f,\](2)

where \(\vec{s}_f(w) = \vec{v}_f'(w)\) is the velocity shear, or the derivative of the velocity with respect to water depth. We solve this PDE using a finite difference method-of-lines approach with an explicit solver (forward Euler and \(4^{th}\) order Runge-Kuta are supported).

Note that, although the transport equation is an advection equation in \(C_f\), we need to consider that \(C_f\) acts as a proxy for \(\eta\). The continuing cycle of disintegration, transport, and lithification makes sure that the change in sediment concentration due to transport and the change in topography are strongly coupled. Then, what seems like an innocent reaction term in Equation 2, turns out to behave as a diffusion equation in \(\eta\). Any attempt at modelling sediment transport where there is an effective down-slope flux combined with some form of disintegration will yield diffusive behaviour.

3.1 Other approaches

In the critical angle approach developed by Warrlich (2000), sediment is transported from unstable slopes to the nearest down-slope stable region. Stability is defined separately for different grain sizes. This method is motivated by the empirical relationship between grain composition and maximum slope angle (Kenter 1990).

The problem with this critical angle-based method of transport is that production across an unstable region is deposited on a small strip, where slopes are below the critical angle. It becomes unclear how to interpret these models from a physics point of view, as results depend heavily on the time-step that is chosen. Contrasting to that, both our transport model and production model (with the exception of the cellular automaton) are discretisations of otherwise continuous processes. This means that, at least assymptotically (i.e. as long as the time step is small enough), our implementation is independent of the chosen time step.

One aspect of critical angle theory that we do use is that we can modulate the disintegration rate \(r_d\) (and therefore the amount of entrained material) with the magnitude of the slope \(|\nabla \eta|\). If we only disintegrate material where the slope is supercritical, the net effect is that sediment is transported from supercritical to stable areas. The difference is that we have a much better control over the physics, and there is no need to convert back and forth between gridded values and a particle representation used in the critical angle approach (e.g. Warrlich 2000). A different approach has been used in the early model CARBPLAT by Bosscher and Southam (1992), which took empirically observed carbonate slopes Adams and Schlager (2000) and defined a slope function that returned slope parameters bounded by the limits of the angle of repose. In the study by Bosscher and Southam (1992) an exponential slope function was assumed, although it should be noted that there is literature debate on the distribution of slope shapes of carbonate platforms (e.g., Schlager and Camber 1986; Kenter 1990; Adams and Schlager 2000). This modelling approach is agnostic with respect to sediment properties and transport mechanisms and optimises the similarity to observed shapes, allowing the user to choose the parameter that produces the best result. However, it does not allow modelling a mixture of sediment types with different properties and requires an a priori assumption on the expected slope shape. It had not been adapted in subsequent models.

3.2 Parameter choices

Our transport model is based on the elementary assumption that sediment flux is proportional to the slope of the sea floor. Nevertheless, we are extrapolating this idea to time scales on which it is hard to reason or otherwise measure the parameters to our model. Especially the combination of transport coefficient \(d_f\), disintegration rate \(r_d\) and lithification time \(t_l\) can be pivotal in acquiring a set of physical outcomes, while we have no good way to estimate acceptable ranges of values for them, other than trying them out and see if the results are plausible.

That being said, by considering some artificial scenarios we can gain more insight into the behaviour of our main parameters.

3.2.1 Disintegration versus subsidence

The chosen disintegration rate will determine wheter our model is on average erosive or accumulative. In the case of a platform morphology, the potential production exceeds the subsidence, meaning that the subsidence rate sets the pace at which the platform grows. At the edge of the platform, there is a transitional region where the maximum production is at pace with the subsidence. If the distintegration rate is higher than the subsidence rate, produced sediment will be immediately disintegrated, stay in the active layer for much longer, and be transported down slope. If the disintegration rate is much lower than the subsidence rate, produced (autochthonous) sediment can accumulate in situ, generating steeper morphologies.

3.2.2 Equilibrium Concentration

The model parametrizes sediment disintegration (i.e. activation or entrainment of older sediments) by a global constant disintegration rate \(r_d\). Entrained sediment is transported by the mechanism described above, and then (re)lithifies by a given percentage every time step. The lithification time is specified as a half-life time \(t_l\). In absence of production, and with infinite available sediment for disintegration, we can see the amount of entrained sediment \(C\) reaching an equilibrium:

\[C(t + \Delta t) = C(t) 2^{-\Delta t / t_l} + r\Delta t,\]

and taking the limit \(\Delta t \to 0\), the equilibrium is reached at,

\[\langle C\rangle = \frac{1}{\ln 2}\ r_d\ t_l.\]

This equilibrium (having units of \(\unit{m}\)) can be useful when estimating the effects of choosing the disintegration rate and lithification time.

3.2.3 Disintegration versus lithification

Both the disintegration rate and the lithification time modulate how long sediment resides in the active layer. By carefully scaling one or the other, the effective diffusion of material can be controlled without changing the specific diffusivity. However, choosing a high lithification time (thus a slow lithification) over a high disintegration rate can help in transporting only freshly produced sediments.

Note that not setting the lithification time (which would amount to immediately depositing all of the active layer on every iteration) results in models that depend heavily on a chosen time step.

To understand the relative effects of choosing a certain lithification time and/or disintegration rate, we ran a one-dimensional model where sediment is produced in a central patch. Then we can study the rate at which sediment is dispersed, either by direct transport before lithification happens, or by subsequent disintegration and re-deposition. By carefully choosing the parameters, we can make a slow lithification process look very similar to a high disintegration rate, as shown in Figure 6.

Comparison between lithification and disintegration

Figure 6: Comparison between lithification and disintegration. The four panes show different combinations of parameters for a one-dimensional model. We have enabled a production of \(100\ \unit{m/Myr}\) for a \(4\ \unit{km}\) wide patch in the middle of the box, and chose a runtime of \(1\ \unit{Myr}\) with a time step of \(100\ \unit{yr}\) (the sharp edges in the production profile induce fast transport, requiring small time steps), and the transport coefficient was set to \(10\ \unit{m/Myr}\). Panels \((a)\) and \((b)\) have a short lithification time \(t_l\) (\(100\ \unit{yr}\)), while panels \((c)\) and \((d)\) have a long lithification time (\(1000\ \unit{yr}\)). On the columns, \((a)\) and \((c)\) have a low disintegration rate \(r_d\) (\(10\ \unit{m/Myr}\)), while \((b)\) and \((d)\) have a high disintegration rate (\(500\ \unit{m/Myr}\)). Values were chosen to have a similar net effect on the dispersion of produced sediment.

runs/ParameterScan.jl 
#| file: runs/ParameterScan.jl
module ParameterScan

function cartesian_product(; pars...)
    ks = keys(pars)
    vs = Iterators.product(values(pars)...)
    collect(ks .=> v for v in vs)
end

kwsplat(f) = d -> f(;d...)

end
#| file: runs/TransportPlots.jl
module TransportPlots

using CarboKitten
using CairoMakie
using Unitful
using Printf

function plot_topography(result)
    fig = Figure()
    ax = Axis(fig[1, 1])
    plot_topography!(ax, result)
    fig
end

function plot_topography!(ax, result)
    diffusivity::typeof(1.0u"m/yr") = result.header.attributes["diffusivity"]
    disintegration_rate::typeof(1.0u"m/Myr") = result.header.attributes["disintegration_rate"]

    dt = result.header.Δt

    ax.xlabel = "x [km]"
    ax.ylabel = "h [m]"

    h0 = result.header.initial_topography[:,1]
    x_axis = result.header.axes.x |> in_units_of(u"km")
    slice = result.data_volumes[:all][:, 1]
    t_axis = result.header.axes.t[1:slice.write_interval:end] |> in_units_of(u"Myr")
    time_steps = div(result.header.time_steps, slice.write_interval)

    ixs = [div(time_steps, 10) + 1, div(time_steps, 2) + 1, time_steps]
    for ix in ixs
        h = (h0 .+ slice.sediment_thickness[:, ix]) |> in_units_of(u"m")
        lines!(ax, x_axis, h, label=@sprintf("%.1f", t_axis[ix]))
    end
end

function plot_matrix(plot!, results, row_names, col_names; fig_pars...)
    fig = Figure(;fig_pars...)
    nrows, ncols = size(results)
    @assert nrows == length(row_names)
    @assert ncols == length(col_names)

    for (i, c) in enumerate(col_names)
        Label(fig[i, 0], c, rotation = pi/2, tellheight=false, tellwidth=true)
    end
    for (i, r) in enumerate(row_names)
        Label(fig[0, i], r, tellheight=true, tellwidth=false)
    end
    axes = [Axis(fig[reverse(Tuple(i))...], title="($('`'+l))")
        for (l,i) in enumerate(eachindex(IndexCartesian(), results))]

    for i in eachindex(IndexCartesian(), results)
        ax = axes[i]
        if i != CartesianIndex(1, 1)
            linkyaxes!(axes[1, 1], ax)
        end
        plot!(ax, results[i])
        if i[2] != ncols
            ax.xlabel = ""
        end
        if i[1] != 1
            ax.ylabel = ""
        end
    end

    Legend(fig[:, nrows+1], axes[1, 1], "time steps [Myr]")

    return fig
end

end
#| file: runs/CustomProductionModel.jl
using CarboKitten
using CarboKitten.Components.Common
using CarboKitten.Components:
    TimeIntegration, Boxes, FaciesBase, SedimentBuffer, WaterDepth,
    Tag, ActiveLayer, Output, Diagnostics
using CarboKitten.Output: Frame
using ModuleMixins: @compose, @for_each

@compose module CustomProduction
@mixin Tag, ActiveLayer, Output, Diagnostics

using CarboKitten.Output: Frame

@kwdef struct Input <: AbstractInput
    production    # a function of (x, y, wd)
end

function write_header(input::AbstractInput, output::AbstractOutput)
    @for_each(P -> P.write_header(input, output), PARENTS)
end

function initial_frame(input::AbstractInput)
    nf = n_facies(input)
    return Frame(
        production=zeros(Sediment, nf, input.box.grid_size...),
        disintegration=zeros(Sediment, nf, input.box.grid_size...),
        deposition=zeros(Sediment, nf, input.box.grid_size...))
end

function initial_state(input::AbstractInput)
    sediment_height = zeros(Height, input.box.grid_size...)
    sediment_buffer = zeros(Float64, input.sediment_buffer_size, n_facies(input), input.box.grid_size...)
    active_layer = zeros(Amount, n_facies(input), input.box.grid_size...)

    state = State(
        step=0, sediment_height=sediment_height,
        sediment_buffer=sediment_buffer,
        active_layer = active_layer)

    return state
end

function step!(input::Input)
    disintegrate! = ActiveLayer.disintegrator(input)
    transport! = ActiveLayer.transporter(input)
    local_water_depth = water_depth(input)
    x, y = box_axes(input.box)
    na = [CartesianIndex()]
    produce(_, wd) = input.production.(x[:,na], y[na,:], wd)[na,:,:]
    pf = lithification_factor(input)

    function (state::State)
        wd = local_water_depth(state)
        p = produce(state, wd)
        d = disintegrate!(state)

        state.active_layer .+= p
        state.active_layer .+= d
        transport!(state)

        deposit = pf .* state.active_layer
        push_sediment!(state.sediment_buffer, deposit ./ input.depositional_resolution .|> NoUnits)
        state.active_layer .-= deposit
        state.sediment_height .+= sum(deposit; dims=1)[1,:,:]
        state.step += 1

        return Frame(
            production = p,
            disintegration = d,
            deposition = deposit)
    end
end

end
#| file: runs/TransportTest.jl
module TransportTest

include("CustomProductionModel.jl")

using CarboKitten
using CairoMakie
using CarboKitten: Box
using CarboKitten: set_attribute
using .CustomProduction: CustomProduction as M

const Time = typeof(1.0u"Myr")

function run_with(;dt, diffusivity, disintegration_rate, lithification_time, patch_width = 2.0u"km")
    facies = [
        M.Facies(
            transport_coefficient=diffusivity)  # 10u"m/yr"
    ]

    box = Box{Periodic{2}}(
        grid_size=(500, 1), phys_scale=30.0u"m")

    t_end = 1.0u"Myr"
    time = TimeProperties(
        Δt = dt,
        steps = (t_end / dt) |> round |> Int)

    @info "Running at Δt = $(time.Δt) and steps = $(time.steps)"
    @info "Production in a single step: $(100.0u"m/Myr" * time.Δt)"


    centre = box.grid_size[1] * box.phys_scale / 2.0
    production(x, y, w) = abs(x - centre) < patch_width ?
        100.0u"m/Myr" * time.Δt :
        0.0u"m"

    write_interval = div(time.steps, 100)

    input = M.Input(
        box=box,
        time=time,
        output = Dict(
            :all => OutputSpec(write_interval=write_interval)),
        initial_topography=(_, _) -> -100.0u"m",
        sea_level=t -> 0.0u"m",
        subsidence_rate=0.0u"m/Myr",
        disintegration_rate=disintegration_rate,
        sediment_buffer_size=50,
        depositional_resolution=0.5u"m",
        lithification_time=lithification_time,
        transport_solver=Val{:forward_euler},
        facies=facies,

        production=production)

    result = run_model(Model{M}, input, MemoryOutput(input))
    set_attribute(result, "diffusivity", diffusivity)
    set_attribute(result, "disintegration_rate", disintegration_rate)
    return result
end

end
#| file: runs/disintegration-vs-lithification.jl
#| classes: ["task"]
#| creates:
#|   - md/fig/disintegration-vs-lithification.pdf
#| requires:
#|   - runs/TransportTest.jl
#|   - runs/TransportPlots.jl
#| collect: figures
module DisintegrationVsCementation

include("TransportPlots.jl")
include("ParameterScan.jl")
include("TransportTest.jl")

using CarboKitten
using CairoMakie
using .ParameterScan: cartesian_product
using .TransportTest: run_with
using .TransportPlots: plot_matrix, plot_topography!
using LaTeXStrings

function main()
    CarboKitten.init()

    pars = (
        diffusivity = [ 10.0 ] * u"m/yr",
        disintegration_rate = [ 10.0, 500.0 ] * u"m/Myr",
        dt = [ 100.0 ] * u"yr",
        lithification_time = [ 100.0, 1000.0 ] * u"yr" )
    cp = cartesian_product(;pars...)

    result = Array{Union{Missing, MemoryOutput}}(missing, size(cp)...)
    Threads.@threads for i in eachindex(cp)
        result[i] = run_with(;cp[i]...)
    end

    fig = plot_matrix(result[1,:,1,:],
            [latexstring("r_d = $(d.val) m/Myr") for d in pars.disintegration_rate],
            [latexstring("t_l = $(d.val) yr") for d in pars.lithification_time];
            fontsize = 10) do ax, result
        plot_topography!(ax, result)
    end

    save("md/fig/disintegration-vs-lithification.pdf", fig)
end

end

DisintegrationVsCementation.main()

3.2.4 Facies-specific transport coefficient \(d_f\)

The facies-specific transport coefficient used in CarboKitten is expressed in \(\unit{m/Myr}\), because it is derived from the parameter \(\nu_f\), transport velocity, that is expressed per unit slope. A diffusion coefficient \(d_f\) appears in \(\partial_t \eta = d_f \nabla^2 \eta\) and have units \(\unit{m^2/Myr}\). In CarboKitten’s formulation \(\nu_f\) is one dimension of length smaller because the active-layer concentration \(C_f\) [m] is already present in the flux. This approach presents two challenges: 1) setting transport coefficients that yield realistic results for carbonate facies modeled at the timescales at which CarboKitten is run, 2) should empirically justified diffusion coefficients for carbonate sediment be available, converting these diffusion coefficients to values of the transport coefficient used in the model. Because advection-diffusion is a modeling approach in carbonate transport and not a direct representation of the physical process of sediment transport, prior empirical diffusion coefficient values are limited. Sultana et al. (2022) reviewed published values, which lie in the range of \(10^5 \unit{m^2/Myr}\) to \(7 \times 10^9 \unit{m^2/Myr}\) (Bosence et al. 1994; Mitchell et al. 1996). Modeling studies differ on the orders of magnitude, which is partly a matter of what processes are accounted for in effective diffusivity, and partly reflects different timescales of measurement. In Dionisos simulations, Sultana et al. (2022) used values ranging from \(1.25 \times 10^6 \unit{m^2/Myr}\) for the sand fraction in the photozoan factory to \(50 \times 10^6 \unit{m^2/Myr}\) for the mud produced by the heterozoan factory and identified 2500 \(10^5\ \unit{m^2/Myr}\) as the upper limit, beyond which no sediment accumulation took place. In a different model, values many orders of magnitude lower have been proposed for the effective diffusion coefficient that implicitly accounts for lithification and depth-dependent wave velocity, introduced by Kaufman et al. (1991):

\[D_{\textrm{Kaufman}}(W) = C_0 \exp(-C_1 W)\]

where \(C_0 = 0.005\ \unit{m^2/Myr}\) for carbonates and \(C_1\) values considered are 0.05 and 0.1 \(\unit{m^{-1}}\), resulting in maximum \(D_{\textrm{Kaufman}}\) values of 0.005 \(\unit{m^2/Myr}\), i.e. much lower than the empirical ones.

Effective sediment diffusion coefficient \(D_{\textrm{eff}}\) values in CarboKitten runs can be estimated from the dispersal of a sediment pulse under any scenario with a given value of transport coefficient, lithification time and disintegration rate. Diffusivity values obtained from runs with a transport coefficient of \(5\ \unit{m/Myr}\) lie in the range of \(3.7\ \times 10^5 \unit{m^2/Myr}\) to \(8.8\ \times 10^6 \unit{m^2/Myr}\) {2}, i.e. well within those reported empirically and overlapping with those used by Sultana et al. (2022) to obtain realistic platform morphologies. Effective \(D_{\textrm{eff}}\) values obtained using this estimate scale linearly with the input transport coefficient \(d_f\).

runs/diffusivity_estimation.jl 
#| file: runs/diffusivity_estimation.jl

# Estimate location and width of the production peak
# using moments of the distribution, based on:
# https://github.com/MindTheGap-ERC/CarboKitten.jl/issues/110#issuecomment-2754962558
# https://gist.github.com/jhidding/1c1f0248de722fb26b356e00aa20ff5a

module DiffusivityEstimation

using CarboKitten.Export: read_slice, Header, DataSlice
using CarboKitten: MemoryOutput
using CarboKitten.Utility: in_units_of
using Unitful

trapezoid(t, y) = (t[2:end] .- t[1:end-1]) .* (y[2:end] .+ y[1:end-1]) ./ 2 |> sum

moment(g, t, y) = trapezoid(t, y .* g.(t))

"""
Compute the location (μ) and width (σ) of a sediment peak from
the spatial profile `h(x)` using moments.

- `x`: position vector (with units)
- `h`: sediment thickness vector (with units)

Returns `(μ, σ)` — the center of mass and standard deviation.
"""
function peak_moments(x, h)
    s = moment(t -> 1, x, h)           # zeroth moment (total mass)
    μ = moment(t -> t, x, h ./ s)      # first moment (mean position)
    v = moment(t -> (t - μ)^2, x, h ./ s)  # second central moment (variance)
    σ = √(v)
    return (μ=μ, σ=σ)
end

"""
    peak_evolution(header::Header, data::DataSlice; use_active_layer=true)

Compute peak location and width at each saved time step
from `Header` and `DataSlice`.

Returns `(t, μ, σ)` vectors.
"""
function peak_evolution(header::Header, data::DataSlice; use_active_layer=true)
    x = header.axes.x
    t = header.axes.t[1:data.write_interval:end]

    n_t = length(t)
    Length = typeof(1.0u"m")
    μ = Vector{Length}(undef, n_t)
    σ = Vector{Length}(undef, n_t)

    if use_active_layer && data.active_layer !== nothing
        # active_layer has shape [n_facies, n_x, n_t] — sum over facies
        al = dropdims(sum(data.active_layer; dims=1); dims=1)
    else
        al = nothing
    end

    for i in 1:n_t
        h = al !== nothing ? al[:, i] : data.sediment_thickness[:, i]
        if all(h .<= 0.0u"m")
            μ[i] = NaN * u"m"
            σ[i] = NaN * u"m"
        else
            m = peak_moments(x, h)
            μ[i] = m.μ
            σ[i] = m.σ
        end
    end

    return (t=t, μ=μ, σ=σ)
end

"""
    peak_evolution(filename, group=:profile; use_active_layer=true)

Read an HDF5 file and compute peak location and width at each saved time step.

Returns `(t, μ, σ)` vectors.
"""
function peak_evolution(filename::AbstractString, group=:profile; use_active_layer=true)
    header, data = read_slice(filename, group)
    peak_evolution(header, data; use_active_layer)
end

"""
    peak_evolution(output::MemoryOutput, group=:profile; use_active_layer=true)

Compute peak location and width at each saved time step
from MemoryOutput.

Returns `(t, μ, σ)` vectors.
"""
function peak_evolution(output::MemoryOutput, group::Symbol=:profile; use_active_layer=true)
    data = output.data_slices[group]
    peak_evolution(output.header, data; use_active_layer)
end

"""
    linear_fit(x, y)

Ordinary least-squares fit of y = a + b·x.
Returns `(a, b, R²)`.
"""
function linear_fit(x, y)
    n = length(x)
    sx = sum(x)
    sy = sum(y)
    sxx = sum(x .* x)
    sxy = sum(x .* y)
    syy = sum(y .* y)
    denom = n * sxx - sx^2
    b = (n * sxy - sx * sy) / denom
    a = (sy - b * sx) / n
    ss_res = syy - a * sy - b * sxy
    ss_tot = syy - sy^2 / n
    R² = 1 - ss_res / ss_tot
    return (a=a, b=b, R²=R²)
end

"""
    _fit_diffusivity(t, μ, σ; R²_threshold=0.99)

Internal: fit the diffusion coefficient from peak evolution data.
"""
function _fit_diffusivity(t, μ, σ; R²_threshold=0.99)

    # Filter out NaN entries (e.g. t=0 with no sediment)
    valid = .!isnan.(σ)
    t_v = t[valid]
    σ_v = σ[valid]
    μ_v = μ[valid]

    # σ² as a function of time
    σ² = σ_v .^ 2

    # Fit σ²(t) = σ₀² + 2D·t, trimming saturated tail
    n = length(t_v)
    n_fit = n
    fit = linear_fit(t_v, σ²)

    while fit.R² < R²_threshold && n_fit > 3
        n_fit -= 1
        fit = linear_fit(t_v[1:n_fit], σ²[1:n_fit])
    end

    D = fit.b / 2

    @info "Estimated effective diffusion coefficient: $D"
    @info "Linear fit R² = $(fit.R²), using $(n_fit)/$(n) valid time steps"
    @info "Fit range: t = $(t_v[1]) to $(t_v[n_fit])"
    @info "Peak center at t_first: $(μ_v[1]), t_last: $(μ_v[n_fit])"
    @info "Peak width σ at t_first: $(σ_v[1]), t_last: $(σ_v[n_fit])"

    return (D=D, R²=fit.R², n_fit=n_fit, t=t, μ=μ, σ=σ)
end

"""
    estimate_diffusivity(filename, group=:profile; use_active_layer=true, R²_threshold=0.99)

Estimate the effective diffusion coefficient from an HDF5 file.
For pure diffusion, σ²(t) = σ₀² + 2Dt,
so D = slope / 2 from a linear fit of σ²(t).

The fit automatically trims late time steps where σ² saturates
(e.g. due to boundary effects) by progressively removing points
from the end until R² exceeds `R²_threshold`.

Returns a named tuple with `D`, `R²`, `n_fit`, `t`, `μ`, `σ`.
"""
function estimate_diffusivity(filename::AbstractString, group=:profile; use_active_layer=true, R²_threshold=0.99)
    t, μ, σ = peak_evolution(filename, group; use_active_layer)
    _fit_diffusivity(t, μ, σ; R²_threshold)
end

"""
    estimate_diffusivity(output::MemoryOutput, group=:profile; use_active_layer=true, R²_threshold=0.99)

Estimate the effective diffusion coefficient from MemoryOutput.
See the `AbstractString` method for details.
"""
function estimate_diffusivity(output::MemoryOutput, group::Symbol=:profile; use_active_layer=true, R²_threshold=0.99)
    t, μ, σ = peak_evolution(output, group; use_active_layer)
    _fit_diffusivity(t, μ, σ; R²_threshold)
end

end # module
#| file: runs/diffusivity_scan.jl
using CarboKitten
using CarboKitten.Models: ALCAP as M
using CairoMakie
using Unitful: ustrip, unit, NoUnits
using DelimitedFiles

include("diffusivity_estimation.jl")

box = CarboKitten.Box{Periodic{2}}(grid_size=(100, 1), phys_scale=50.0u"m")

t_end = 1.0u"Myr"
Δt = 50.0u"yr"
t_steps = t_end / Δt |> ceil |> Int
peak_centre = box.phys_scale * box.grid_size[1] ÷ 2
peak_width = 200.0u"m"
peak_height = 10.0u"m"
write_interval = max(1, t_steps ÷ 1000)

facies1 = M.Facies(
    viability_range = (0, 0),
    activation_range = (0, 0),
    production = BenthicProduction(
        maximum_growth_rate = 0u"m/Myr",
        extinction_coefficient = 0u"m^-1",
        saturation_intensity = 60u"W/m^2"),
    transport_coefficient = 5u"m/yr",
    wave_velocity = _ -> (Vec2(0.0u"m/Myr", 0.0u"m/Myr"), Vec2(0.0u"1/Myr", 0.0u"1/Myr")),
    initial_sediment = (x, _) -> peak_height * exp(-(x - peak_centre)^2/(2 * peak_width^2)),
)

function make_input(; cementation_time, disintegration_rate)
    M.Input(
        facies = [facies1],
        box = box,
        time = TimeProperties(Δt=Δt, steps=t_steps),
        output = Dict(:profile => OutputSpec(slice=(:, 1), write_interval=write_interval)),
        lithification_time = cementation_time,
        disintegration_rate = disintegration_rate,
        subsidence_rate = 0.0u"m/Myr",
        sea_level = _ -> 0.0u"m",
        initial_topography = (_, _) -> -100.0u"m",
        insolation = 0.0u"W/m^2",
        transport_solver = Val{:forward_euler},
        depositional_resolution = 1.0u"km",
        sediment_buffer_size = 2,
    )
end

parameters = Dict(
    :cementation_time    => [1000, 2500, 5000] .* u"yr",
    :disintegration_rate => [5.0, 10.0, 20.0, 50.0] .* u"m/Myr",
)

function cartesian_product(pars::Dict{Key,Vector}) where {Key}
    if isempty(pars)
        return [ Dict{Key, Any}() ]
    end
    pars = copy(pars)
    result = []
    k, vs = first(pairs(pars))
    for item in cartesian_product(delete!(pars, k))
        for v in vs
            push!(result, merge(item, Dict(k => v)))
        end
    end
    return result
end

ct_values = parameters[:cementation_time]
dr_values = parameters[:disintegration_rate]
D_values   = Matrix{Any}(undef, length(ct_values), length(dr_values))

for (i_ct, ct) in enumerate(ct_values)
    for (i_dr, dr) in enumerate(dr_values)
        inp    = make_input(cementation_time=ct, disintegration_rate=dr)

        output = run_model(Model{M}, inp, MemoryOutput(inp))

        fig = Figure()
        ax  = Axis(fig[1, 1])
        x   = output.header.axes.x |> in_units_of(u"km")
        y   = output.data_slices[:profile].sediment_thickness |> in_units_of(u"m")
        for i in [1, 100, 1000]
            lines!(ax, x, y[:, i])
        end

        ct_val = round(Int, ustrip(u"yr", ct))
        dr_val = round(ustrip(u"m/Myr", dr), digits=2)
        save("data/diffusivity_scan/ct_$(ct_val)yr_dr$(dr_val)mMyr.png", fig)

        est = DiffusivityEstimation.estimate_diffusivity(output, use_active_layer=false)
        D_values[i_ct, i_dr] = est.D
    end
end

# strip units for plotting
D_unit   = unit(first(D_values))
D_matrix = ustrip.(D_values)
ct_axis  = ustrip.(u"yr",   ct_values)
dr_axis  = ustrip.(u"m/Myr", dr_values)

# save D_matrix to CSV: rows = cementation_time, cols = disintegration_rate
open("data/diffusivity_scan/D_matrix.csv", "w") do io
    write(io, "cementation_time_yr," * join(["dr_$(d)mMyr" for d in dr_axis], ",") * "\n")
    for (i, ct) in enumerate(ct_axis)
        write(io, "$(ct)," * join(D_matrix[i, :], ",") * "\n")
    end
end

fig_summary = Figure()
ax = Axis(fig_summary[1, 1],
    xlabel = "cementation time [yr]",
    ylabel = "disintegration rate [m/Myr]",
    title  = "Estimated diffusivity for transport coefficient = $(facies1.transport_coefficient)",
    aspect = length(ct_axis) / length(dr_axis))
hm = heatmap!(ax, ct_axis, dr_axis, D_matrix)
Colorbar(fig_summary[1, 2], hm, label = "D [$D_unit]")
save("data/diffusivity_scan/D_summary.png", fig_summary)
Table 2: Estimated effective diffusivity \(D_{\textrm{eff}}\ [\unit{m^2 Myr^{-1}}]\) for combinations of lithification time \(t_l\) and disintegration rate \(r_d\) at facies transport coefficient \(d_f\) equal to 5 \(\unit{m/yr}\).
Lithification time \(t_l\) \(r_d = 5\ \unit{m/Myr}\) \(r_d = 10\ \unit{m/Myr}\) \(r_d = 20\ \unit{m/Myr}\) \(r_d = 50\ \unit{m/Myr}\)
1000 yr 36,565 72,237 137,595 299,711
2500 yr 84,136 165,496 301,921 540,862
5000 yr 156,904 300,574 503,548 879,985

3.3 Implementation and limitations

Our implementation of the transport model first computes the gradient of the sea floor (or equivalently the water depth) \(\vec{\nabla} w\) using central differences. From this gradient we can compute the advection coefficients in Equation 2, \(\vec{c}_{\textrm{adv}} = d_f \vec{\nabla} \eta + \vec{v}_f(w)\). The maximum advection coefficient sets the Courant number and determines how many time steps we need to take to solve Equation 2. For an advection equation integrated with the forward Euler method, we need

\[|c_{\textrm{adv}}|_{\infty} \frac{\Delta t}{\Delta x} \le 1.\]

This states that we cannot move matter further than a single pixel distance in one iteration, or our computation becomes unstable.

In practice, we compute the transport coefficients and the maximum resulting advective Courant number. Then we integrate the transport equation adaptively, choosing the minimum number of subdivided time steps that keeps the advection stable. Since we computed the transport coefficients in advance, it is relatively cheap to apply multiple iterations of the advection solver, for which we use a first-order upwind scheme.

Now consider our transport model in the context of the larger carbonate platform model. Each time we disintegrate some matter which gets entrained and transported as part of the sediment concentration \(C_f\), after which a fraction is cemented, increasing \(\eta\). If we consider \(\partial{\eta}/\partial{t} \sim \partial{C}/\partial{t}\), then part of the transport equation is the diffusion equation \(\partial{\eta}/\partial{t} \sim d_f C_f \nabla^2 \eta\). This leaves our implementation vulnerable to instabilities when the global time step is taken too large. For just this diffusion term the CFL limit is

\[d_f C_f \frac{\Delta t}{(\Delta x)^2} \le 1.\]

This means that increasing the resolution of a model by a factor two may need a time step four times smaller for the integration to remain stable. CarboKitten has a diagnostic mode where this condition is checked against, allowing the user to make informed changes to the input parameters. Because the sediment concentration is not known in advance, it is not possible to make this check in advance.

4 Validation

We validated the model by attempting to replicate the general architecture of an isolated carbonate platform in Central Luconia (part of the Sarawak Basin), offshore Malaysia (Cycle IV to Cycle Lower V). We chose this area is because this platform has been extensively studied (Henglai et al. 2024) and the interpreted seismic profiles and reconstructed sea-level curves were publicly reported. We informed model runs with the sea-level data obtained from Figure 2 of Henglai et al. (2024). In addition, the tectonic setting of the region in the studied period (around \(15.48\ \unit{Ma}\) to \(11\ \unit{Ma}\)) is believed to be uncomplicated (Fathiyah Jamaludin 2025) and thus possible to represent with a constant subsidence rate. Four wells were drilled into the atoll, providing discontinous cores and, thus, geological facies information.

This platform has a pinnacle shape, with the length of approximately \(5\ \unit{km}\) in NE-SW direction, suggesting it was drowned in during this time period. The cores indicate that the top of carbonate platform is dominated by coral boundstones. Although the sampling in the lower part of the cores is discontinous, the cores suggest that the facies are dominated by pack- to wackestones and interbeded by coral boundstones, bioclastic boundstones and argillaceous carbonate. We related the geological facies to the carbonate producers via Table 3.

Table 3: Facies used in validation. We used three producing facies in modelling our validation case, and a fourth facies type to track transported coral.
Geological Facies Producers Model Producer Production rate
Coral boundstones Euphotic coral (T factory) Facies 1 \(1800\ \unit{m\ Myr^{-1}}\)
Pack- to wackestone Oligophotic algae (M factory) Facies 2 \(800\ \unit{m\ Myr^{-1}}\)
Argillaceous carbonate Pelagic (P factory) Facies 3 \(8\ \unit{Myr^{-1}}\)
Bioclastic boundstone Transported coral Facies 4 n.a.

Three carbonate producers are fed into the model, corresponding to three geological facies (producers). Their respective production rates are listed in Table 3. A fourth facies is the transported facies of euphotic corals. These values are within the production rates estimated by (Lopez-Gamundi et al. 2025). A pre-run was set in order to generate an atoll-like initial topography for the actual run. The total duration of modelling covers the entire Cycle IV (from \(15.8\ \unit{Ma}\) to \(11.8\ \unit{Ma}\), having a duration of \(4\ \unit{Myr}\)).

The criteria for comparison between the modelled results and the interpreted seismic profiles were:

  1. We obtain five major zones within the platform.
  2. The platform should show a pinnacle shape.
  3. The total thickness of carbonate stratigraphy is approximately (\(600\ \unit{m}\)).
  4. The platform has a slope indicating presence of off-shore sediments transport.

The results are listed in Figure 7. We compared the features in the figure and the interpreted seimic profile:

  1. The interpreted seismic profiles suggested five zones, and this is captured in the resultant figure (five parasequences).
  2. The results displayed a ‘pinacle’ shape, showing upper parasequences bearing smaller horizonal dimensions than the lower parasequences.
  3. The total thickness of carbonate (\(\sim 550\ \unit{m}\)) is comparable to the data suggested by well penetration (\(\sim 600\ \unit{m}\)).
  4. The slope is presented in the modelled results, as a mixture of facies 3 and 4.

The model yields nearly identical - i.e., qualitatively the same - results for different timesteps. Users should note that changing the timestep requires adjusting the CA interval, which has been applied for the tested versions of the code: a timestep of 100, 250 and 250 yr. Results for a time step of 400 yr are shown in 7.

Case_validation

Figure 7: Validation case. A: This shows a cross-section of our model of the platform in Central Luconia, offshore Malaysia. The white lines indicate unconformities, and the black lines are coeval lines at regular intervals, every solid black line separating one million years. The colours indicate the dominant facies type of the modelled sediment, blue for facies 1 (T factory), yellow for facies 2 (M factory), green for facies 3 (P factory) and pink represents facies 4 (transported coral, i.e. formerly facies 1). B: Seismic interpretation of the platform, redrawn from Figure 9 of Henglai et al. (2024).

Therefore, we claim CarboKitten is able to reproduce major architectural features of a real-life carbonate platform.

5 Software design

5.1 Box topology

runs/topology_coast.jl 
#| file: runs/topology_coast.jl
#| classes: ["task"]
#| creates: data/topology_coast.h5
module TopologyCoast

using CarboKitten
using CarboKitten.Production
using CarboKitten.Models: ALCAP as M

function main()
    CarboKitten.init()

    facies = [
        M.Facies(
            production=Production.EXAMPLE[:euphotic],
            transport_coefficient=10.0u"m/yr"),
        M.Facies(
            production=Production.EXAMPLE[:oligophotic],
            transport_coefficient=10.0u"m/yr"),
        M.Facies(
            production=Production.EXAMPLE[:aphotic],
            transport_coefficient=10.0u"m/yr")
    ]

    sea_level(t) =
        10.0u"m" * sin(2π * t / 123456.0u"yr") +
         5.0u"m" * sin(2π * t /  23456.0u"yr")

    kwargs = (
        sea_level = sea_level,
        subsidence_rate=50.0u"m/Myr",
        disintegration_rate=50.0u"m/Myr",
        insolation=400.0u"W/m^2",
        sediment_buffer_size=50,
        depositional_resolution=0.5u"m",
        transport_solver=Val{:forward_euler},
    )

    coast_input = M.Input(
        time = TimeProperties(
            Δt=0.0002u"Myr",
            steps=5000),
        box = Box{Coast}(
            grid_size=(250, 50),
            phys_scale=60.0u"m"),
        facies = facies,
        initial_topography = (x, y) -> let x_prime = x - 3.0u"km"
            x_prime > 0.0u"km" ? - x_prime / 300.0 : - x_prime / 30.0
        end,
        output = Dict(
            :topography => OutputSpec(write_interval = 1000),
            :profile => OutputSpec(slice = (:, 25)));
        kwargs...)

    run_model(Model{M}, coast_input, "data/topology_coast.h5")
end

end

TopologyCoast.main()
#| file: runs/Noise.jl
module Noise
using FFTW

"""
    make_noise(box::Box, n, s, σ)

Make some Gaussian Random Noise with a power spectrum of

``P(k) = (s k)^n \\exp(-2 π^2 k^2 σ^2)``

The `box` should have periodic boundaries, `n` be a unitless
number (usually between -2.0 and 2.0), `s` is a scaling in meters
which only affects the amplitude and is needed to make `s * k` unitless,
and `σ` is the standard deviation of the Gaussian filter to reduce
small scale noise.

The noise is first generated as Gaussian white noise, then convolved in
Fourier space by multiplying with the square root of the power spectrum.

Example:

    box = Box{Periodic{2}}(grid_size=(100, 100), phys_scale=300.0u"m")
    cnoise = make_noise(box, -1.5, 50.0u"m", 500.0u"m")

"""
function make_noise(box, n, s, σ)
    white_noise = randn(box.grid_size...)
    P(k) = (k * s)^n * exp(-π^2 * k^2 * 2 * σ^2)
    kx = FFTW.rfftfreq(box.grid_size[1], 1/box.phys_scale)
    ky = FFTW.fftfreq(box.grid_size[2], 1/box.phys_scale)
    kabs = sqrt.(kx.^2 .+ ky'.^2)

    fy = FFTW.rfft(white_noise)
    p  = P.(kabs)
    p[1] = 0.0
    fy .*= sqrt.(p)
    FFTW.irfft(fy, box.grid_size[1])
end
end
#| file: runs/topology_periodic.jl
#| classes: ["task"]
#| creates:
#|   - data/topology_periodic.h5
#| requires:
#|   - runs/Noise.jl
include("Noise.jl")

module TopologyPeriodic

using CarboKitten
using CarboKitten.Production
using CarboKitten.Models: ALCAP as M
using ..Noise: make_noise

function main()
    CarboKitten.init()

    facies = [
        M.Facies(
            production=Production.EXAMPLE[:euphotic],
            transport_coefficient=10.0u"m/yr"),
        M.Facies(
            production=Production.EXAMPLE[:oligophotic],
            transport_coefficient=10.0u"m/yr"),
        M.Facies(
            production=Production.EXAMPLE[:aphotic],
            transport_coefficient=10.0u"m/yr")
    ]

    sea_level(t) =
        10.0u"m" * sin(2π * t / 123456.0u"yr") +
         5.0u"m" * sin(2π * t /  23456.0u"yr")

    kwargs = (
        sea_level = sea_level,
        subsidence_rate=50.0u"m/Myr",
        disintegration_rate=50.0u"m/Myr",
        insolation=400.0u"W/m^2",
        sediment_buffer_size=50,
        depositional_resolution=0.5u"m",
        transport_solver=Val{:forward_euler},
    )

    box = Box{Periodic{2}}(grid_size=(256, 256), phys_scale=60.0u"m")
    initial_topography = make_noise(box, -1.5, 5.0u"m", 1.0u"km") .* 1.0u"m" .- 30.0u"m"

    coast_input = M.Input(
        time = TimeProperties(
            Δt=0.0002u"Myr",
            steps=5000),
        box = box,
        facies = facies,
        initial_topography = initial_topography,
        output = Dict(
            :topography => OutputSpec(write_interval = 500),
            :profile_x => OutputSpec(slice = (:, 125)),
            :profile_y => OutputSpec(slice = (125, :)));
        kwargs...)

    run_model(Model{M}, coast_input, "data/topology_periodic.h5")
end

end

TopologyPeriodic.main()
#| file: runs/topology_plot.jl
#| classes: ["task"]
#| creates:
#|   - md/fig/topologies.png
#| requires:
#|   - data/topology_coast.h5
#|   - data/topology_periodic.h5
#| collect: figures
module TopologyPlot

using CairoMakie
using CarboKitten.Export: read_volume
using CarboKitten.Visualization: glamour_view!

function main()
    coastal_header, coastal_data =
        read_volume("data/topology_coast.h5", :topography)
    periodic_header, periodic_data =
        read_volume("data/topology_periodic.h5", :topography)

    fig = Figure(size=(800, 500))

    glamour_view!(
        Axis3(fig[1, 1], title="(a) periodic boundaries",
              width=250, height=250),
        periodic_header, periodic_data,
        colormap=Reverse(:GnBu))

    glamour_view!(
        Axis3(fig[1, 2], title="(b) coastal boundaries"),
        coastal_header, coastal_data,
        colormap=Reverse(:GnBu))

    ax1 = Axis(fig[2, 1], title="(c)", aspect=1.0)
    hlines!(ax1, [0.0, 15.0])
    vlines!(ax1, [0.0, 15.0])
    arrows2d!(ax1,
              [Point(7.5, 0.0), Point(7.5, 15.0),
               Point(0.0, 7.5), Point(15.0, 7.5)],
              [Vec(0.0, 2.0), Vec(0.0, 2.0),
               Vec(2.0, 0.0), Vec(2.0, 0.0)],
              color=[:red, :red, :blue, :blue])

    ax2 = Axis(fig[2, 2], title="(d)", aspect=DataAspect(), height=100)
    hlines!(ax2, [0.0, 5.0])
    vlines!(ax2, [0.0, 15.0])
    arrows2d!(ax2,
              [Point(7.5, 0.0), Point(7.5, 5.0),
               Point(0.0, 2.5), Point(0.0, 2.5),
               Point(15.0, 2.5), Point(15.0, 2.5)],
              [Vec(0.0, 0.8), Vec(0.0, 0.8),
               Vec(1.0, 0.0), Vec(-1.0, 0.0),
               Vec(1.0, 0.0), Vec(-1.0, 0.0)],
              color=[:red, :red, :blue, :blue, :green, :green])

    save("md/fig/topologies.png", fig)
    return fig
end

end

TopologyPlot.main()

CarboKitten needs to work with different choices for box topology, i.e. how the boundaries of a model box connect to each other. For example, when we simulate a small strip of coastline it is best to have one axis (in this case the \(x\)-axis) reflect onto itself, while the other axis is periodic, leaving fewer edge effects.

In another case, where we want to simulate an entire island, or even an archipelago, it is more convenient to use fully periodic coordinates. We illustrate these choices in Figure 8.

Different topologies, 3d view and boundaries

Figure 8: Model topologies. CarboKitten allows the user to choose different topologies for the spatial modelling. In panel (a) we see a group of reef islands that were modelled on a fully periodic grid of size \(250 \times 250\), using a randomly generated initial topography. A more common use case is shown in panel (b), where the \(x\) coordinate is reflected at the boundaries, while the \(y\) coordinate is periodic, thus modelling a small strip of coastline. Here the grid size is \(250 \times 50\), and the initial topography is a linearly declining slope of \(0.3\%\) (with the exception of the shore, which is steeper). The superimposed surfaces represent different equidistant time steps in the same runs, while the colors indicate the elevation along with the values on the \(z\)-axis, here provided to guide the viewer. Panels (c) and (d) schematically illustrate these same box topologies using coloured arrows.

5.2 The sediment buffer

In our models of sediment transport and denudation it is important to remember the sedimentation history for all produced facies for some time into the past. We keep a three-dimensional fixed-size buffer, where two dimensions represent the \(x\) and \(y\) coordinates of the system, and the third dimension discretizes the amount of deposited material. Each cell in the buffer represents a parcel of sediment, where we store the relative fractions of each contributing facies. The data-structure acts like a stack, working on a First-In-First-Out (FIFO) basis.

We emphasise that this buffer is only used to determine the facies composition of disintegrated sediment. The sediment output of the overall model is written to disk at each iteration for post-analysis. This means that the model output can be much more precise than the depositional resolution of the buffer.

The user can set the size of the buffer as well as the amount of sediment represented by each cell. The depth of the buffer is controlled by an input parameter sediment_buffer_size, while the resolution is set through the depositional_resolution parameter in units of \(\unit{m}\) of sediment.

The rows in the buffer represent a constant amount of sediment. An alternative approach is to have rows that represent time slices. In that case, when we need to disintegrate an amount of sediment, we need to search the buffer back in time until enough sediment is collected. This can be very slow, and it also means that one needs to have the full sedimentation history in memory. A buffer that is discretized on depth however does not have those requirements, at the expense of a small amount of facies mixing.

While the sediment buffer is allocated as a single 4-dimensional array (depth, facies, \(x\), \(y\)), it is best to explain its functioning from the perspective of a single cell in our model. We are left with two dimensions: depth (rows) and facies (columns).

We choose to have the head of our sediment stack always be at the first row. When sediment out-grows the buffer, the deepest layers are dropped from memory. The head can contain an incomplete amount of sediment, while all rows below the head are either full or empty. When sediment is pushed to the stack and the head row overflows, all rows are copied down one row and the surplus is assigned to the now empty head row. The inverse happens when removing (popping) material from the stack. This process is illustrated below in Figure 9.

Sediment buffer diagram

Figure 9: Above we see a buffer. Suppose we first produce some sediment amounting to \(3/4\) of the buffer resolution, and after that an amount of \(1/2\) gets disintegrated for further transport. In the beginning, the top cell of our buffer is \(1/2\) full. First we push a parcel of (relative) size \(3/4\). The remaining \(1/4\) overflows to the next cell. Then we pop an amount of \(1/2\). The first \(1/4\) of this amount is retrieved from the top cell, and the second \(1/4\) from the cell below that. This popped parcel will have different fractions from the pushed one, since it also draws from the half filled row that was in the stack before pushing. In this sense, a small amount of facies mixing will take place, depending on the depositional resolution chosen.

5.3 User interface

The user interfaces CarboKitten by writing a Julia script that defines the relevant model parameters and runs the chosen model. Effectively, very little Julia needs to be known to take an example input and modify parameters. Output is written to HDF5 files for post-processing and visualization.

CarboKitten ships with routines for visualisation and data extraction into CSV files. This makes it easier for novice users to use results from CarboKitten in further processing pipelines that rely on other programming languages. Data extracted includes sediment accumulation curves, age-depth models, water depth, and stratigraphic columns with facies code, allowing to test a wide range of hypotheses. These include, but are not limited to, testing hypotheses on orderedness of strata (Burgess 2016), preservation orbital forcing (Kemp et al. 2016), proxy records (Curtis et al. 2025), or preservation of biotic information such as patterns of origination and extinction, biostratigraphic precision, and evolutionary change (Hohmann and Jarochowska 2025; Hohmann et al. 2024; Holland and Patzkowsky 2002).

5.4 Performance

Since CarboKitten is written in Julia with performance in mind, it should be efficient to run, even on consumer grade hardware, i.e. an average laptop. We are yet to substantiate this claim. Since Julia is a just-in-time compiled language, the first execution of any code in a new session always takes a bit longer than subsequent runs. Measurements presented in this section do not include this initial overhead.

5.4.1 Baseline

Our baseline model is the example included in CarboKitten, grid size \(100 \times 50\) with 5000 time steps of 200 years each (results shown in Figure 5). This model runs in 27 seconds on a Intel Core i7 at \(3.0\ \unit{GHz}\).

With regards to memory consumption, CarboKitten allocates a fixed amount of memory at the start of a model run, which scales linearly with the size of the grid. The most significant fraction of the memory is occupied by the sediment buffer. In the example run we have a buffer size of 50. With three facies types being stored this results in an array size of \(100 \times 50 \times 50 \times 3\), stored in double precision gives a mere \(6 \unit{MB}\). However, for a \(300 \times 300\) sized grid this already increases to \(108 \unit{MB}\).

5.4.2 Scaling

The run-time and memory consumption of CarboKitten should scale linearly with the number of pixels in the grid, with two complicating factors. Firstly, for smaller models the run-time can become limited by many smaller writes to HDF5. For those cases we provide a method of running models entirely in-memory. The second complication is the transport model. Here run times may vary due to the number of integration steps required for stability reasons. Increasing the resolution of a model also means increasing the number of transport integration time steps required by the same factor (considering the CFL condition for advective transport). Transport efficiency is also affected by the local topography: increasing the slope also increases the number of integration steps required. Carbonate platforms have the tendency to generate steep slopes due to exponential sedimentation rates in the production model. These steep slopes can be mitigated by setting the transport coefficient \(d_f\): higher values will lead to smoothed-out and gentler slopes. On the other hand, modelling on-shore transport due to wave transport can induce steeper slopes, again requiring smaller integration time steps. Note that we are speaking of integration steps of the transport model, which can be any integer fraction of a full model time step. When the transport model needs too many steps for every model step, we can start to question the accuracy of the model as a whole, and the user should try decreasing the time-step of the full model to compensate.

5.4.3 Benchmark

To further quantify these complications in our estimated run-times, we run a model of a single atoll at three resolutions (\(200, 100\), and \(50\ \unit{m}\), corresponding to grid sizes of \(75^2, 150^2, 300^2\)) with three step sizes (\(400, 200\), and \(100\ \unit{yr}\), corresponding to 2500, 5000, and 10000 steps), for a total of nine benchmark cases. We set the interval of the cellular automaton to compensate for the number of time steps. This way, runs with the same grid size should have very similar output. The results are shown in Figure 10.

The combination of 2500 time steps with a \(300^2\) grid size yields instabilities in the transport model and is left out of the results. Other than that, CarboKitten scaled as predicted from our previous considerations.

Benchmark plots

Figure 10: Benchmark with respect to number of time steps and grid size. Panel (a) shows the run-time dependency on the number of time-steps, while panel (b) shows the dependency on the number of grid cells on each axis, both on a log-log scale. This scaling follows the predicted behaviour: linear in both the number of time-steps and total number of grid cells (on this plot being the grid size squared). Note that the run with 2500 time steps and \(300^2\) grid size is left out, since the transport model was unstable for that configuration. These numbers were consistent throughout multiple runs.

5.4.4 Validation

We may validate our benchmark by looking at the results of the runs with grid size \(150^2\). This is shown in Figure 11. These results show that, when time steps are taken small enough, CarboKitten converges to a consistent result that does not depend on the size of the time step.

Benchmark validation

Figure 11: Benchmark validation. This shows a cross-section of the runs with a grid size of \(150^2\). Looking at the first output, using only 2500 time steps, we see a wave like pattern even where the deep sea facies dominate. These waves are not physical, but a result from taking the time step too large. When we look at the results from 5000 and 10000 time steps, they look so similar that we can conclude that in this case 5000 steps was enough to get accurate results.

runs/benchmarks.jl 
#| file: runs/benchmarks.jl
module Benchmarks

using BenchmarkTools
using CarboKitten
using GeometryBasics
using DataFrames
using CSV

# A constant homogeneous wave velocity
v_const(v_max) = _ -> (Vec2(v_max, 0.0u"m/yr"), Vec2(0.0u"1/yr", 0.0u"1/yr"))

initial_topography(x, y) =
    - sqrt((x - 7.5u"km")^2 + (y - 7.5u"km")^2) / 100.0

const FACIES = [
    ALCAP.Facies(
        viability_range=(4, 10),
        activation_range=(6, 10),
        production = BenthicProduction(
            maximum_growth_rate=500u"m/Myr",
            extinction_coefficient=0.8u"m^-1",
            saturation_intensity=60u"W/m^2"
        ),
        transport_coefficient=50.0u"m/yr",
        wave_velocity=v_const(-2.0u"m/yr")),
    ALCAP.Facies(
        viability_range=(4, 10),
        activation_range=(6, 10),
        production = BenthicProduction(
            maximum_growth_rate=400u"m/Myr",
            extinction_coefficient=0.1u"m^-1",
            saturation_intensity=60u"W/m^2"
        ),
        transport_coefficient=10.0u"m/yr",
        wave_velocity=v_const(-0.5u"m/yr")),
    ALCAP.Facies(
        viability_range=(4, 10),
        activation_range=(6, 10),
        production = BenthicProduction(
            maximum_growth_rate=100u"m/Myr",
            extinction_coefficient=0.005u"m^-1",
            saturation_intensity=60u"W/m^2"
        ),
        transport_coefficient=10.0u"m/yr",
        wave_velocity=v_const(-2.0u"m/yr"))
]

box(res) = Box{Periodic{2}}(
    grid_size=(res, res),
    phys_scale=15.0u"km" / res)

time(steps) = TimeProperties(
    Δt=1.0u"Myr" / steps,
    steps=steps)

output(res, steps) = Dict(
    :topography => OutputSpec(write_interval = max(1, div(steps, 50))),
    :profile    => OutputSpec(slice = (:, div(res, 2))))

ca_interval(steps) = max(div(steps, 5000), 1)

sea_level(t) =
    10.0u"m" * sin(2π * t / 123456.0u"yr") +
     5.0u"m" * sin(2π * t /  23456.0u"yr")

input(res, steps) = ALCAP.Input(
    tag="atoll_$(res)_$(steps)",
    box=box(res),
    time=time(steps),
    output=output(res, steps),
    ca_interval=ca_interval(steps),
    initial_topography=initial_topography,
    sea_level=sea_level,
    subsidence_rate=50.0u"m/Myr",
    disintegration_rate=50.0u"m/Myr",
    lithification_time=100.0u"yr",
    insolation=400.0u"W/m^2",
    sediment_buffer_size=50,
    depositional_resolution=0.5u"m",
    transport_solver=Val{:forward_euler},
    facies=FACIES)


function cartesian_product(pars::Dict{Key,Vector{Elem}}) where {Key, Elem}
    if isempty(pars)
        return [ Dict{Key, Any}() ]
    end

    pars = copy(pars)
    result = []

    k, vs = first(pairs(pars))

    for item in cartesian_product(delete!(pars, k))
        for v in vs
            push!(result, merge(item, Dict(k => v)))
end
    end

    return result
end

const BM = @NamedTuple{value::String, time::Float64, bytes::Int64, alloc::Int64, gctime::Float64}

function run_benchmark(; res, steps)
    output_file = "data/bench_$(res)_$(steps).h5"
    try
        bm = @btimed run_model(Model{ALCAP}, $(input(res, steps)), $output_file)
        return (res=res, steps=steps, bm...)
    catch e
        return (res=res, steps=steps, value="error", time=NaN, bytes=0, alloc=0, gctime=NaN)
    end
end

function main()
    CarboKitten.init()
    run_model(Model{ALCAP}, input(50, 1000), "data/bench_init.h5")

    pars = Dict(
        :res => [75, 150, 300],
        :steps => [2500, 5000, 10000] )

    results = DataFrame(res=Int[], steps=Int[], value=String[], time=Float64[], bytes=Int64[], alloc=Int64[], gctime=Float64[])

    for p in cartesian_product(pars)
        push!(results, run_benchmark(; p...))
    end

    CSV.write("data/benchmark.csv", results)
end

end

Benchmarks.main()
#| file: runs/benchmark_plot.jl
#| classes: ["task"]
#| requires:
#|   - data/benchmark.csv
#| creates:
#|   - md/fig/benchmark.pdf
#| collect: figures
module BenchmarkPlot

using CairoMakie
using AlgebraOfGraphics
using DataFrames
using CSV

function main()
    df = CSV.read("data/benchmark.csv", DataFrame)
    fig = Figure()

    layer1 = data(df) * mapping(:steps => "time steps", :time => "run time [s]", color=:res => string => "grid size") * visual(ScatterLines)
    layer2 = data(df) * mapping(:res => "grid size", :time => "run time [s]", color=:steps => string => "time steps") * visual(ScatterLines)

    fig = Figure(size=(800, 400))
    fg1 = draw!(fig[1, 1], layer1, axis=(xscale=Makie.pseudolog10, yscale=log10))
    legend_args = (tellwidth=false, tellheight=false, halign=:left, valign=:top, margin=(10, 10, 10, 10))
    legend!(fig[1, 1], fg1; legend_args...)
    fg2 = draw!(fig[1, 2], layer2, axis=(xscale=Makie.pseudolog10, yscale=log10))
    legend!(fig[1, 2], fg2; legend_args...)

    Label(fig[1, 1, TopLeft()], "a", halign=:left, fontsize=20)
    Label(fig[1, 2, TopLeft()], "b", halign=:left, fontsize=20)

    save("md/fig/benchmark.pdf", fig)

    fig
end

end

BenchmarkPlot.main()
#| file: runs/benchmark_validation.jl
module BenchmarkValidation

using CairoMakie
using CarboKitten.Visualization: sediment_profile!
using CarboKitten.Export: read_slice

function main()
    h1, d1 = read_slice("data/bench_150_2500.h5", :profile)
    h2, d2 = read_slice("data/bench_150_5000.h5", :profile)
    h3, d3 = read_slice("data/bench_150_10000.h5", :profile)

    fig = Figure(size = (1200, 400))
    sediment_profile!(Axis(fig[1, 1], title="2500 steps"), h1, d1)
    sediment_profile!(Axis(fig[1, 2], title="5000 steps"), h2, d2)
    sediment_profile!(Axis(fig[1, 3], title="10000 steps"), h3, d3)

    fig.content[1].title = "2500 steps"
    fig.content[2].title = "5000 steps"
    fig.content[3].title = "10000 steps"

    save("md/fig/benchmark_validation.png", fig)

    fig
end

end

BenchmarkValidation.main()

5.4.5 Potential for GPU optimisation

At the time of writing, CarboKitten is a single threaded CPU code. However, the structure of the model, is highly ammenable to optimisation on a GPU, which would drastically improve run-times further. Going through the steps of the composed model in Section 2.5: the cellular automaton is a stencil operation, production a map, disintegration a stencil, transport is implemented as an iterated stencil, and deposition is a map. Both stencil and map operations are highly localized in memory and are ideal for implementation on a GPU.

5.5 Documentation

CarboKitten is written entirely using literate programming (Knuth 1984). This means that the implementation of CarboKitten is written as an integral part of its own documentation, using a system called Entangled (Hidding 2023). The aim is that interested readers have a direct reference to the code implementing the methods that are explained in the documentation.

6 Examples

We provide two examples of typical tasks users are likely to undertake: creating a simulation using an empirical, externally provided sea level curve and one with explicit forcing by an insolation curve. These examples are supported by the code used to generate this executable manuscript (FIXME ref to the code) and users are encourage to use that code as a starting point for modifying these examples for their needs.

The third example serves to illustrate the details of how the wave-induced transport is modelled and how modelling decisions and parameter choices affect the outcomes obtained using this feature.

6.1 Sea level

Variables external to the model, which modulate the output the most, are the sea level and insolation. The sea level, together with subsidence, result in the relative sea level, which translates into water depth at any given position in the basin. The sea level must be specified as a function of time. It can be a constant, a continuous function or an empirical dataset. Empirical datasets can be read in as text files and need to be interpolated to equidistant intervals corresponding to the time step with which the model is run. The example here uses the sea level curve by Lisiecki and Raymo (2005), reproduced in the compilation by Miller et al. (2005). The dataset of relative sea level records derived from foraminifer \(\delta^{18}O\) extracted from this compilation is included in CarboKitten to facilitate simulations of the most typical sea-level scenarios. In this example we start the model at \(2\ \unit{Ma}\) and build the platform until \(134.54\ \unit{ka}\), i.e. until the end of the record by Lisiecki and Raymo (2005), using a time step of 200 y (Figure 12).

variable_SL 
#| id: variable_SL
#| file: runs/variable_sl.jl
#| creates: md/fig/variable-sl.png
module VariableSL

using CarboKitten
using DelimitedFiles: readdlm
using Unitful
using DataFrames
using Interpolations
using CategoricalArrays
using CarboKitten.DataSets: artifact_dir
using CairoMakie
using CarboKitten.Visualization: sediment_profile!
using CarboKitten.Export: read_slice

function miller_2020()
    dir = artifact_dir()
    filename = joinpath(dir, "Miller2020", "Cenozoic_sea_level_reconstruction.tab")

    data, header = readdlm(filename, '\t', header=true)
    return DataFrame(
        time=-data[:,4] * u"kyr",
        sealevel=data[:,7] * u"m",
        refkey=categorical(data[:,2]),
        reference=categorical(data[:,3]))
end

function sea_level()
    df = miller_2020()
    lisiecki_df = df[df.refkey .== "846 Lisiecki", :]
    lisiecki_df = filter(row -> -2.0u"Myr" <= row.time, lisiecki_df)
    sort!(lisiecki_df, [:time])

    return linear_interpolation(
        lisiecki_df.time,
        lisiecki_df.sealevel)
end

const TIME_PROPERTIES = TimeProperties(
    t0 = -1999.7u"kyr",
    Δt = 100.0u"yr",
    steps = 18650
)

const TAG = "lisiecki-sea-level"

const FACIES = [
    ALCAP.Facies(
        production=BenthicProduction(
            maximum_growth_rate=200u"m/Myr",
            extinction_coefficient=0.8u"m^-1",
            saturation_intensity=60u"W/m^2"),
        transport_coefficient=20.0u"m/yr",
        name="euphotic"),
    ALCAP.Facies(
        production=BenthicProduction(
            maximum_growth_rate=500u"m/Myr",
            extinction_coefficient=0.1u"m^-1",
            saturation_intensity=60u"W/m^2"),
        transport_coefficient=10.0u"m/yr",
        name="oligophotic"),
    ALCAP.Facies(
        production=BenthicProduction(
            maximum_growth_rate=100u"m/Myr",
            extinction_coefficient=0.005u"m^-1",
            saturation_intensity=60u"W/m^2"),
        transport_coefficient=50.0u"m/yr",
        name="aphotic")
]

const INPUT = ALCAP.Input(
    tag="$TAG",
    box=CarboKitten.Box{Coast}(grid_size=(200, 50), phys_scale=150.0u"m"),
    time=TIME_PROPERTIES,
    ca_interval=1,
    initial_topography=(x, y) -> -x / 200.0 + 20.0u"m",
    sea_level=sea_level(),
    output=Dict(
        :profile => OutputSpec(slice = (:, 25), write_interval = 1)),
    subsidence_rate=5.0u"m/Myr",
    disintegration_rate=50.0u"m/Myr",
    insolation=500.0u"W/m^2",
    sediment_buffer_size=50,
    depositional_resolution=0.5u"m",
    facies=FACIES)

function main()
    CarboKitten.init()
    run_model(Model{ALCAP}, INPUT, "data/variable-sl.h5")
end

function plot(result_file)
    header, result_profile = read_slice(result_file, :profile)

    inch = 96
    pt = 4/3
    cm = inch / 2.54

    fig = Figure(size=(20cm, 12cm), fontsize=8pt)
    ax_left = Axis(fig[1, 1])
    ax_right = Axis(fig[1, 2])
    colsize!(fig.layout, 1, Relative(0.2))
    Label(fig[1, 1, TopLeft()], "(a)")
    Label(fig[1, 2, TopLeft()], "(b)")

    sl_fn = sea_level()
    times = collect(time_axis(TIME_PROPERTIES))
    sl_values = [ustrip(u"m", sl_fn(t)) for t in times]
    times_myr = ustrip.(u"Myr", times)

    lines!(ax_left, sl_values, times_myr)
    ax_left.xlabel = "Sea level [m]"
    ax_left.ylabel = "Time [Myr]"

    sediment_profile!(ax_right, header, result_profile, show_unconformities=50)

    save("md/fig/variable-sl.png", fig, px_per_unit=300/inch)
end

end

# result = VariableSL.main()
result = "data/variable-sl.h5"
VariableSL.plot(result)
variable-sl

Figure 12: Platform generated using the sea level curve of Lisiecki et al. (2005). Panel (a) shows the sea level curve. Panel (b) the corresponding output stratigraphy.

6.2 Insolation

The relationship between production and insolation can be modified with user-provided parameters. It may be confusing that the extinction coefficient \(k\) is, in CarboKitten, a property of the carbonate factory and the facies it deposits and not of the basin or position in it. In reality extinction coefficient varies for different wavelengths of the sunlight spectrum, but the set of its values across the spectrum is constant for a given water body. While different carbonate factories exploit (or ignore, in the case of the aphotic factory) different parts of the light spectrum, the model is agnostic to it and allows users to set \(k\) to values that may represent an average across different producers using different wavelengths.

As default, we use insolation of \(400\ \unit{W/m^2}\), which is approximately equivalent to the \(2000\ \unit{\mu E m^{-2} s^{-1}}\) used by Bosscher and Schlager (1992). This is representative of insolation on the sea surface at midday in the tropics. However, insolation varies with the position of the Earth with respect to the Sun and over geological timescales this variation may affect the patterns of sediment production. Incoming solar radiation can be used as an input vector to modulate production. CarboKitten is agnostic with respect to the source of this information. As an example in Figure 13, here we use the daily mean insolation on June solstice, calculated using the astronomical solution by Laskar (2004), obtained through the R package palinsol (Crucifix 2023). Here we obtain it for the coming million year (starting in 1950, which is when the astronomical solution starts) at the 25° N latitude and use the total solar irradiance value of 1361 \(\unit{kW m^{-2}}\). Variation in solar irradiance is so small that it would hardly manifest itself if linearly propagated to the sea level curve. A universal transfer function describing the relationship between insolation and sea level does not exist. For the purpose of illustrating the functionality of the model, we calculate the sea level as an amplified insolation value. The amplification is chosen arbitrarily as the square of the insolation anomaly, with the anomaly being the deviation from mean irradiation.

runs/extract_insolation.R 
#| file: runs/extract_insolation.R

if (!require("palinsol")) {
    install.packages("palinsol", repos = "https://cran.r-project.org")
}

library(palinsol)
time_end <- 500000
time_start <- 0
time_step <- 200
times <- seq(time_start, time_end, time_step)
param_la04 = t(sapply(times, function(t) astro(t, solution = la04, degree = TRUE)))
orbit <- list()
insolation <- list()
lat_degree = 25

for (t in 1:length(times)) {
  orbit[[t]] <- list(
    eps = param_la04[t,1] * pi / 180,
    ecc = param_la04[t,2],
    varpi = (param_la04[t,3] - 180) * pi / 180
  )

  insolation[[t]] <- Insol(
    orbit[[t]],
    long = pi / 2,
    lat = lat_degree * pi / 180,
    S0 = 1361,
    H = NULL
  )
}

insolation = inso_values <- unlist(insolation)
write.csv(insolation, file="data/insolation.csv", row.names = FALSE)

The insolation file can be read into a CarboKitten script defining the model to be run. The alternative is calling R directly from Julia using RCall.jl.

runs/insolation_run.jl 
#| file: runs/insolation_run.jl
module Insolation

using CarboKitten
using DelimitedFiles: readdlm
using Unitful
using Interpolations
using CairoMakie
using CarboKitten.Visualization: sediment_profile, sediment_profile!
using Statistics

function import_insolation(file::String)
    dir = "data"
    filename = joinpath(dir, file)
    insolation = readdlm(filename, '\t', header=false, skipstart=1)
    vec(insolation) .|> Float64
end

const TIME_PROPERTIES = TimeProperties(
    t0 = 0u"Myr",
    Δt = 200.0u"yr",
    steps = length(import_insolation("insolation.csv"))-1
)

function get_insolation(times::Vector, insolation::Vector)
    interpolator = linear_interpolation(times, insolation)
    return t -> interpolator(ustrip(u"yr", t)) * u"W/m^2"
end

const TAG = "insolation-future"

const FACIES = [
    ALCAP.Facies(
        production=BenthicProduction(
            maximum_growth_rate=200u"m/Myr",
            extinction_coefficient=0.8u"m^-1",
            saturation_intensity=60u"W/m^2"),
        transport_coefficient=50.0u"m/yr",
        name="euphotic"),
    ALCAP.Facies(
        production=BenthicProduction(
            maximum_growth_rate=500u"m/Myr",
            extinction_coefficient=0.1u"m^-1",
            saturation_intensity=60u"W/m^2"),
        transport_coefficient=25.0u"m/yr",
        name="oligophotic"),
    ALCAP.Facies(
        production=BenthicProduction(
            maximum_growth_rate=100u"m/Myr",
            extinction_coefficient=0.005u"m^-1",
            saturation_intensity=60u"W/m^2"),
        transport_coefficient=12.5u"m/yr",
        name="aphotic")
]

const time_vector = collect(time_axis(TIME_PROPERTIES)) / u"yr" .|> NoUnits
const insolation_vector = import_insolation("insolation.csv")

function get_sea_level(times::Vector, insolation::Vector)
    insolation_anomaly = (insolation .- mean(insolation)) ./ mean(insolation)
    sea_level_anomaly = -(100 .* (insolation_anomaly)) .^ 2
    sea_level_values = -100.0 .+ sea_level_anomaly
    interpolator = linear_interpolation(times, sea_level_values)
    return t -> interpolator(ustrip(u"yr", t)) * u"m"
end

const INPUT = ALCAP.Input(
    tag="$TAG",
    box=CarboKitten.Box{Coast}(grid_size=(100, 50), phys_scale=150.0u"m"),
    time=TIME_PROPERTIES,
    ca_interval=1,
    initial_topography=(x, y) -> -x / 200.0 - 100.0u"m",
    sea_level = get_sea_level(time_vector, insolation_vector),
        output=Dict(
        :profile => OutputSpec(slice=(:, 25), write_interval=1)),
    subsidence_rate=50.0u"m/Myr",
    disintegration_rate=50.0u"m/Myr",
    insolation= get_insolation(time_vector, insolation_vector),
    sediment_buffer_size=50,
    depositional_resolution=0.5u"m",
    facies=FACIES)

    function main()
        run_model(Model{ALCAP}, INPUT, MemoryOutput(INPUT))
    end

    function plot(result::MemoryOutput)
        inch = 96
        pt = 4/3
        cm = inch / 2.54

        fig = Figure(size=(20cm, 12cm), fontsize=8pt)
        ax_left = Axis(fig[1, 1])
        ax_right = Axis(fig[1, 2])
        colsize!(fig.layout, 1, Relative(0.2))
        Label(fig[1, 1, TopLeft()], "(a)")
        Label(fig[1, 2, TopLeft()], "(b)")

        sl_fn = get_sea_level(time_vector, insolation_vector)
        times = collect(time_axis(TIME_PROPERTIES))
        sl_values = [ustrip(u"m", sl_fn(t)) for t in times]
        times_myr = ustrip.(u"Myr", times)

        lines!(ax_left, sl_values, times_myr)
        ax_left.xlabel = "Sea level [m]"
        ax_left.ylabel = "Time [Myr]"
        ax_left.xticks = Makie.LinearTicks(4)

        sediment_profile!(ax_right, result.header, result.data_slices[:profile], show_unconformities=10)

        save("md/fig/variable-insolation.png", fig, px_per_unit=300/inch)
    end

end

result = Insolation.main()
Insolation.plot(result)
variable-insolation

Figure 13: Platform generated using the daily mean insolation during June solstice at the 25° N latitude for a period of 1 Myr starting in 1950 and using a sea level curve obtained by amplifying the insolation values. Panel (a) shows the insolation-based amplified sea level curve. Panel (b) shows the corresponding output stratigraphy.

6.3 Wave induced transport

We model the transport by waves by setting the velocity \(v_f\) and shear \(s_f\) components in Equation 2. In this analysis we forego claims on any level of realism with respect to the true long term effects of wave transport, rather we study the behaviour of the model under an imposed additional velocity component. Our use of the term wave transport should also be understood as such.

Considering the long timescales we are working with, we limit ourselves to a highly simplified model, with the goal of achieving an effect comparable with that of wave-induced transport. Given the timescales for which the model is developed, with time steps of the order of \(100\) years, a more physical representation of wave-induced transport is not possible. By necessity, the result imitates the time-averaged effect of tranport.

Our approach is illustrated with an example of an atoll, starting with a conical topography and periodic boundaries.

Here we try three different velocity profiles: first no onshore component, second a constant vector that does not depend on water depth, and third an attempt at a more realistic scenario.

The following equation is the well known phase velocity of waves as a function of depth from linear wave theory,

\[v(w) = \sqrt{\frac{\lambda g}k} {\rm tanh} (k w),\]

where \(w\) is the water depth, \(k\) the wave number (\(k = 2\pi / \lambda\)), and \(g\) is the gravitational acceleration. This velocity is the phase-velocity of surface waves, given the total depth of the water. To evaluate the transport velocity at deeper levels, we multiply the phase velocity with a factor \(\exp(-kw)\) to account for Stokes drift:

\[v_f = A_f \exp (- k w) \tanh (k w),\](3)

where \(A_f\) is the facies-dependent maximum transport velocity. The \(k\) parameter can be tweaked to set the depth at which the maximum transport velocity is attained. We assume most of the sediment transport happens close to the sea floor. This profile is chosen for its assymptotic properties: at high water depth the transport velocity converges to zero, while the decrease in wave velocity towards shallow depths ensures that there is a net influx of material close to the shore. An example of this profile is shown in Figure 14.

Depth profile

Figure 14: Depth profile of wave velocity and shear. The velocity profile was taylored to have a maximum of \(10\ \unit{m/yr}\) at a depth of \(20\ \unit{m}\). Where the shear is negative (assuming transport is directed onshore), there is a net accumulation of sediment.

velocity-profile 
#| id: velocity-profile
v_prof(v_max, max_depth, w) =
    let k = sqrt(0.5) / max_depth,
        A = 3.331 * v_max,
        α = tanh(k * w),
        β = exp(-k * w)
        (A * α * β, -A * k * β * (1 - α - α^2))
    end
#| classes: ["task"]
#| creates: md/fig/wave-transport-magnitude.pdf
#| collect: figures

module Script

using CairoMakie
using Unitful

<<velocity-profile>>

inch = 96
pt = 4/3
cm = inch / 2.54

function main()
    w = LinRange(0, 100.0, 1000)u"m"
    f = v_prof.(10.0u"m/yr", 20.0u"m", w)

    v = first.(f)
    s = last.(f)

    fig = Figure(size=(8.3cm, 6cm), fontsize=8pt)
    ax1 = Axis(fig[1, 1], title="transport velocity", yreversed=true, xlabel="velocity [m/yr]", ylabel="depth [m]")
    ax2 = Axis(fig[1, 2], title="transport shear", yreversed=true, xlabel="shear [1/yr]", ylabel="depth [m]")
    lines!(ax1, v / u"m/yr", w / u"m")
    lines!(ax2, s / u"1/yr", w / u"m")
    save("md/fig/wave-transport-magnitude.pdf", fig)
end

end

Script.main()
runs/atoll.jl 
#| file: runs/atoll.jl
#| classes: ["task"]
#| creates: data/atoll.h5
#| collect: atoll

module Atoll

using CarboKitten
using GeometryBasics
using .Threads: @threads

<<velocity-profile>>

v_none(_) = w -> (Vec2(0.0, 0.0) * u"m/yr", Vec2(0.0, 0.0) * u"1/yr")

v_prof(v_max) = w -> let (v, s) = v_prof(v_max, 10.0u"m", w)
        (Vec2(v, 0.0u"m/Myr"), Vec2(s, 0.0u"1/Myr"))
    end

v_flat(v_max) = _ -> (Vec2(v_max, 0.0u"m/yr"), Vec2(0.0u"1/yr", 0.0u"1/yr"))

initial_topography(x, y) =
    - sqrt((x - 7.5u"km")^2 + (y - 7.5u"km")^2) / 100.0

facies(v) = [
    ALCAP.Facies(
        production=BenthicProduction(
            maximum_growth_rate=500u"m/Myr",
            extinction_coefficient=0.8u"m^-1",
            saturation_intensity=60u"W/m^2"),
        transport_coefficient=20.0u"m/yr",
        wave_velocity=v(-2.0u"m/yr")),
    ALCAP.Facies(
        production=BenthicProduction(
            maximum_growth_rate=400u"m/Myr",
            extinction_coefficient=0.1u"m^-1",
            saturation_intensity=60u"W/m^2"),
        transport_coefficient=10.0u"m/yr",
        wave_velocity=v(-0.5u"m/yr")),
    ALCAP.Facies(
        production=BenthicProduction(
            maximum_growth_rate=100u"m/Myr",
            extinction_coefficient=0.005u"m^-1",
            saturation_intensity=60u"W/m^2"),
        transport_coefficient=50.0u"m/yr",
        wave_velocity=v(-2.0u"m/yr"))
]

box(res) = Box{Periodic{2}}(
    grid_size=(res, res),
    phys_scale=15.0u"km" / res)

time(steps) = TimeProperties(
    Δt=1.0u"Myr" / steps,
    steps=steps)

output(res, steps) = Dict(
    :topography => OutputSpec(write_interval = max(1, div(steps, 50))),
    :profile    => OutputSpec(slice = (:, div(res, 2)+1)))

ca_interval(steps) = max(div(steps, 5000), 1)

sea_level(t) =
    10.0u"m" * sin(2π * t / 123456.0u"yr") +
     5.0u"m" * sin(2π * t /  23456.0u"yr")

input(name, res, steps, v) = ALCAP.Input(
    tag="atoll_$(res)_$(name)",
    box=box(res),
    time=time(steps),
    output=output(res, steps),
    ca_interval=ca_interval(steps),
    initial_topography=initial_topography,
    sea_level=sea_level,
    subsidence_rate=50.0u"m/Myr",
    disintegration_rate=50.0u"m/Myr",
    lithification_time=100.0u"yr",
    insolation=400.0u"W/m^2",
    sediment_buffer_size=50,
    depositional_resolution=0.5u"m",
    transport_solver=Val{:forward_euler},
    ca_random_seed=2,
    facies=facies(v))

function main()
    CarboKitten.init()
    vs = [
        :none => v_none,
        :flat => v_flat,
        :prof => v_prof
    ]

    @threads for (n, v) in vs
        inp = input(n, 150, 5000, v)
        run_model(Model{ALCAP}, inp, "data/$(inp.tag).h5")
    end
end

end

Atoll.main()
#| file: runs/atoll-profile-plot.jl
# #| classes: ["task"]
#| requires: data/atoll.h5
#| creates: md/fig/atoll-profile.png
#| collect: figures

using CairoMakie
using CarboKitten.Visualization: sediment_profile
using HDF5
using CarboKitten.Export: read_slice, read_header

h5open("data/atoll.h5") do fid
    header = read_header(fid)
    data = read_slice(fid["profile"])
    fig = sediment_profile(header, data)
    save("md/fig/atoll-profile.png", fig)
end
#| file: runs/atoll-map-plot.jl
#|# classes: ["task"]
#|# requires: data/atoll_150_prof.h5
#|# creates: md/fig/atoll-map.png
#|# collect: figures
module AtollMapPlot

using HDF5
using Unitful
using CairoMakie
using CarboKitten: in_units_of
using CarboKitten.Visualization: sediment_profile!
using CarboKitten.Export: read_volume, read_slice

inch = 96
pt = 4/3
cm = inch / 2.54

function topography!(ax, header, data; colormap=:nuuk, levels=[-10, -5, 0, 5, 10])
    ax.aspect = DataAspect()

    h = header
    s = data.sediment_thickness[:, :, end]
    t = h.initial_topography .+ s .- (h.axes.t[end] * h.subsidence_rate)
    hm = contourf!(
        ax,
        h.axes.x |> in_units_of(u"km"),
        h.axes.y |> in_units_of(u"km"),
        t / u"m", levels=levels, colormap=colormap, extendlow=:auto, extendhigh=:auto)

    contour!(
        ax,
        h.axes.x |> in_units_of(u"km"),
        h.axes.y |> in_units_of(u"km"),
        t / u"m", levels=levels, color=:black, labels=true)
    return hm
end

function main()
    fig = Figure(size=(20cm, 12cm), fontsize=8pt)

    res= 150
    levels = -5:2
    limits = (5.5, 9.5)
    tics = 5:9

    kwargs = (
        limits=(limits, limits), aspect=DataAspect(),
        xticks=tics, yticks=tics, xlabel="x [km]", ylabel="y [km]")
    ax1 = Axis(fig[1, 1], title="(a) no onshore transport"; kwargs...)
    h1, d1 = read_volume("data/atoll_$(res)_none.h5", :topography)
    ax2 = Axis(fig[1, 2], title="(b) constant velocity"; kwargs...)
    h2, d2 = read_volume("data/atoll_$(res)_flat.h5", :topography)
    ax3 = Axis(fig[1, 3], title="(c) depth dependent velocity"; kwargs...)
    h3, d3 = read_volume("data/atoll_$(res)_prof.h5", :topography)

    _, p1 = read_slice("data/atoll_$(res)_none.h5", :profile)
    ax4 = Axis(fig[2,1], title="(d)")
    _, p2 = read_slice("data/atoll_$(res)_flat.h5", :profile)
    ax5 = Axis(fig[2,2], title="(e)")
    _, p3 = read_slice("data/atoll_$(res)_prof.h5", :profile)
    ax6 = Axis(fig[2,3], title="(f)")

    hm1 = topography!(ax1, h1, d1, levels=levels)
    topography!(ax2, h2, d2, levels=levels)
    topography!(ax3, h3, d3, levels=levels)
    Colorbar(fig[1, 4], hm1, label="height [m]")
    sediment_profile!(ax4, h1, p1)
    sediment_profile!(ax5, h2, p2)
    sediment_profile!(ax6, h3, p3)

    ax4.title = "(d)"
    ax4.limits = ((3, 12), (-100, 5))
    ax5.title = "(e)"
    ax5.limits = ((3, 12), (-100, 5))
    ax6.title = "(f)"
    ax6.limits = ((3, 12), (-100, 5))
    save("md/fig/atoll-map.png", fig)
end

end

AtollMapPlot.main()
Atoll topography and sediment profile

Figure 15: Topography and sediment profiles of an atoll. We ran the same model three times with different on-shore velocity profiles: no on-shore transport, a constant velocity and lastly the profile given in Equation 3. The top row (panels a, b, and c) show the topography of the generated island, while the bottom row (panels d, e, and f) show the corresponding sediment profiles. Small differences in water depth may get amplified exponentially by the production model, so we see some stark differences in the outcomes for the different velocity profiles. Comparing the first (without additional transport vector) and second case (flat profile), we see little change in the overall shape of the atoll, but there is a clear difference in the facies composition at the transition between oligophotic (yellow) and aphotic (green) dominated areas. In the third case we see the topography changed significantly between the leeward and windward sides of the atoll, where the slope is much steeper. Also the facies composition changed further: most notably we see a relative prominence of euphotic (blue) facies on the windward side of the island.

We model the formation of an atoll for three cases: no wave transport, constant transport directed west-ward (along the x-axis), and a depth dependent velocity profile. The results of this experiment are shown in Figure 15. Velocity functions are configured for each facies separately. We found that it was quite easy to create an unstable model by choosing on-shore velocities too high, particularly in the case where the velocity shear is non-zero. Build-up of material due to high on-shore velocity can be compensated by setting a higher facies transport coefficient.

If we assume that both the facies transport coefficient and the wave velocity are some constant times a hypothetical carrying water velocity, it would be fitting to make sure that for each facies the transport coefficient and velocity have a similar proportion. In our experiment we took the values listed in Table 4. It is very hard to find proper motivation for any of these values, but by changing them we can learn more about the mechanisms and systematic behaviours of the model and by extension possibly learn more about the formation of carbonate platforms and study sensitivities in their observed stratigraphic patterns.

Table 4: chosen parameters to generate results in Figure 15. We used a lithification time of \(100 \unit{yr}\) and a disintegration rate of \(50 \unit{m/Myr}\).
facies transport coefficient \([\unit{m/yr}]\) velocity \([\unit{m/yr}]\)
euphotic 20.0 2.0
oligophotic 10.0 0.5
aphotic 50.0 2.0

7 Conclusions

CarboKitten is a new Open Source stratigraphic forward model dedicated for carbonate depositional environments and modeling of timescales between centuries and millions of years. It integrates previous, well-tested approaches used by the community, i.e. the production model by Bosscher and Schlager (1992) and the generation of spatial heterogeneity proposed by Burgess (2013) with a new approach to sediment transport, based on the concept of the active layer by Paola et al. (1992). The software allows modeling and visualization accessible to laptop users, including attractive plotting functions for common use-cases in stratigraphy and sedimentology, such as Wheeler diagrams, age-depth models and stratigraphic columns. CarboKitten uses heuristics to approximate the dynamics of carbonate production, wave transport and biologically driven spatial heterogeneity. The algorithms do not replicate the physical and biological processes behind these phenomena, but allow obtaining results imitating them at timescales at which they cannot be observed directly.

CarboKitten fills a gap identified at the outset of this manuscript: an Open Source, modular, performant stratiraphic forward model dedicated to carbonate platforms and their idiosyncracies. It complements the landscape of Earth-surface models such as Landlab (Hobley et al. 2017; Barnhart et al. 2020) and the broader suite of tools developed and maintained by the Community Surface Dynamics Modeling System (CSDMS). Future opportunities for CarboKitten include integration with Landlab and other CSDMS components.

At this stage, CarboKitten’s primary value lies not in a realistic replication of empirical stratigraphic architectures, but in its utility for testing hypotheses on the formation of the carbonate geological record and on our own understanding on its governing processes. Further work is needed to allow more detailed reconstructions of known geological situations. Among future refinements are empirical validation of transport and production values or storing the history of sediment transport to track autochthonous and allochthonous sediment.

CarboKitten offers a powerful tool to ground-truth concepts of how time is represented in the physical rock record (e.g., Burgess 2008; Burgess et al. 2019; Sultana et al. 2022) and constrain the limits of reconstruction of processes such as evolution (Holland and Patzkowsky 1999; Hannisdal 2006; Hohmann et al. 2024), climate change, or other aspects of the changing Earth’s environment (e.g., Kemp et al. 2016; Kemp and Van Manen 2019; Myrow and Grotzinger 2000; Geyman et al. 2021; Husinec et al. 2023; Curtis et al. 2025). We hope the accessibility and reproducibility of CarboKitten simulations will encourage wider use of stratigraphic forward models towards a hypothetico-deductive research in stratigraphy.

Code and Data Availability

CarboKitten is available under the GNU Public Licencse 3.0 and is hosted on Github. Releases are also made available on Zenodo (doi:10.5281/zenodo.14051612), see Hidding et al. (2025). The full source code for this paper including all of its figures is also made available on Zenodo (doi:10.5281/zenodo.17135660).

Appendix

A Derivation of transport equations

Our basic assumption is that the sediment flux scales linearly with the local bathymetric gradient and the concentration of sediment in the active layer,

\[\vec{q}_f(x) = -d_f C_f(x) \vec{\nabla}\eta(x) + \vec{v}_f(w(x))) C_f,\]

where \(C_f\) is the active sediment amount per facies (all \(f\) suffixes indicate facies dependent quantities), \(d_f\) the transport (diffusivity) coefficient, \(\eta(x)\) the bathymetry, and \(\vec{v}_f\) the wave velocity as a function of water depth \(w(x)\). Since the water depth and bathymetry differ at any time by a constant an a minus sign, \(\vec{\nabla}w = -\vec{\nabla}\eta\), it is advantageous to write the equation completely in terms of \(w\),

\[\vec{q}_f(x) = C_f(x)\ (d_f\ \vec{\nabla}w(x) + \vec{v}_f(w(x))).\]

We can transform this assumption in to an advection equation by considering the continuity equation,

\[\frac{\partial C_f(x)}{\partial t} = -\vec{\nabla} \cdot \vec{q}_f(x) + P(x),\]

where \(P(x)\) is the sediment production rate (including disintegration). We can leave the production out of consideration for the moment. Also, for readability, we will be dropping the function dependencies on \(x\) and in our notation, as well as the \(f\) suffix, considering each facies independently. Using the product rule,

\[\frac{\partial C}{\partial t} = -\vec{\nabla} C\cdot(d\vec{\nabla}w + \vec{v}(w)) - C\ (d\nabla^2 w + \vec{\nabla} \cdot \vec{v}(w))).\]

Now, we demand that the user provide the derivative \(\vec{s}(w) = \vec{v}'(w)\), so applying the chain rule we can write,

\[\vec{\nabla}\cdot\vec{v}(w) = \vec{\nabla}w \cdot \vec{s}(w).\]

Substituting that into the previous equation and collecting terms by \(C\) and \(\nabla C\) brings us to Equation 2,

\[\frac{\partial C}{\partial t} = -\big(d\vec{\nabla}w + \vec{v}(w)\big)\cdot\vec{\nabla} C - \big(d\nabla^2 w + \vec{\nabla}w \cdot \vec{s}(w))\big)\ C.\]

B Model parameters

For completeness we give an overview of the paramaters that enter CarboKitten’s main model called ALCAP (Active Layer, Cellular Automaton and Production), at the time of writing the most extensive model in CarboKitten. We split these parameters in three categories. The global parameters deal with the grid geometry, time integration and external forcings (see Table 5). A second set of parameters are specified for each facies individually, e.g. production rates and CA behaviour (see Table 6). The third category concerns everything related to the transport model (see Table 7). Note that the transport coefficients are facies dependent, all others are global parameters.

Author Contributions

EJ, JH, PB and NH conceptualized the study. Software and visualizations were developed by JH with contributions from EJ, XL, NH, and HS. All authors contributed to the methodology. EJ was responsible for funding acquisition and project administration. JH and EJ drafted the manuscript and PB and NH contributed to the final version.

Competing Interests

The authors declare that they have no conflict of interest.

Acknowledgements

We thank Joris Eggenhuisen for discussions on the transport model and Charlotte Summers for programming support. Niels Drost provided administrative and management support during the project. We thank Ton Markus for the visualization of Figure 7.

Funded by the European Union (ERC, MindTheGap, StG project no 101041077). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

Table 5: Grid and time parameters. Parameters marked with \(\dagger\) can be specified as a constant, table or Julia function. Parameters that have \(-\) as default value are required.
Parameter Description Default Unit Tested range
Time parameters
time.t0 Simulation start time 0.0 Myr 0.0
time.\(\Delta\)t Time step — Myr 10 \(\unit{yr}\) – 0.001 \(\unit{Myr}\)
time.steps Number of time steps — — 100 – 50000
Box parameters
box.grid_size Number of grid cells (x, y) — — (100,1) – (512,512)
box.phys_scale Physical size of each cell — \(\unit{m}\) 50 – 600 \(\unit{m}\)
Water depth
sea_level Relative sea-level function \(\dagger\) \(t \to 0.0\) \(\unit{m}\) —
initial_topography Initial seafloor elevation \(\dagger\) \((x,y) \to 0\) m example: \(x \to -x/300.0\)
subsidence_rate Rate of tectonic subsidence 0.0 \(\unit{m\ Myr^{-1}}\) \(0.0\) – \(50.0\ \unit{m\ Myr^{-1}}\)
Insolation
insolation Surface insolation \(\dagger\) — \(\unit{Wm^{-2}}\) 400 \(\unit{Wm^{-2}}\)
Table 6: Facies parameters.
Parameter Description Default Unit Tested range
Benthic Production
maximum_growth_rate Maximum growth rate 0.0 m/Myr 100 – 500 m/Myr
extinction_coefficient Light attenuation coefficient 0.0 \(\unit{m^{-1}}\) 0.005 – 0.8 \(\unit{m^{-1}}\)
saturation_intensity Half-saturation light intensity 1.0 \(\unit{Wm^{-2}}\) 50 – 60 \(\unit{Wm^{-2}}\)
Pelagic Production
maximum_growth_rate Maximum pelagic growth rate 0.0 \(\unit{Myr^{-1}}\) 7.0 \(\unit{Myr^{-1}}\)
extinction_coefficient Light attenuation coefficient 0.0 \(\unit{m^{-1}}\) 0.1 \(\unit{m^{-1}}\)
saturation_intensity Half-saturation light intensity 1.0 \(\unit{Wm^{-2}}\) 60 \(\unit{Wm^{-2}}\)
Cellular Automaton
ca_interval Update CA every N time steps 1 — 1
ca_random_seed Random seed for initial CA state 0 — 0
viability_range Min–max neighbour count to stay alive (4, 10) — (4, 10)
activation_range Min–max neighbour count to be born (6, 10) — (6, 10)
active Whether facies participates in CA true — true, false
Table 7: Transport parameters. \(\star\): The user needs to specify both velocity and its derivative as function of water depth.
Parameter Description Default Unit Tested range
Transport coefficients facies dependent
transport_coefficient Facies-specific transport coefficient 0.0 m/yr 1.0 – 50.0 m/yr
wave_velocity Facies advection velocity and shear \(\star\) zero m/Myr custom function
Control parameters
intertidal_zone Upward extension of the intertidal zone 0.0 m 0.0
disintegration_rate Rate at which sediment is mobilised 50.0 m/Myr 5 – 500 m/Myr
disintegration_transfer Remapping disintegrated sediment identity — custom redistribution
lithification_time Half-life time for active-layer settling Myr 100 yr – 5000 yr
Transport solver
transport_solver ODE solver forward Euler — :forward_euler, :RK4
transport_substeps Number of sub-steps per iteration adaptive — :adaptive or integer
Sediment buffer
sediment_buffer_size Depth of sediment buffer 50 — 2 – 150
depositional_resolution Resolution of sediment buffer 0.5 m 0.5 – 1.0 m
Adams, Erwin W., and Wolfgang Schlager. 2000. “Basic Types of Submarine Slope Curvature.” Journal of Sedimentary Research 70 (4): 814–28. https://doi.org/10.1306/2DC4093A-0E47-11D7-8643000102C1865D.
Barnhart, K. R., E. W. H. Hutton, G. E. Tucker, et al. 2020. “Short Communication: Landlab V2.0: A Software Package for Earth Surface Dynamics.” Earth Surface Dynamics 8 (2): 379–97. https://doi.org/10.5194/esurf-8-379-2020.
Barrett, Samuel J., and Jody M. Webster. 2017. “Reef Sedimentary Accretion Model (ReefSAM): Understanding Coral Reef Evolution on Holocene Time Scales Using 3D Stratigraphic Forward Modelling.” Marine Geology 391: 108–26. https://doi.org/10.1016/j.margeo.2017.07.007.
Bosence, Dan W. J., Luis Pomar, David A. Waltham, and Thomas H. G. Lankester. 1994. “Computer Modeling a Miocene Carbonate Platform, Mallorca, Spain1.” AAPG Bulletin 78 (2): 247–66. https://doi.org/10.1306/BDFF9078-1718-11D7-8645000102C1865D.
Bosscher, Hemmo, and Wolfgang Schlager. 1992. “Computer Simulation of Reef Growth.” Sedimentology 39 (3): 503–12.
Bosscher, Hemmo, and John Southam. 1992. “CARBPLAT—a Computer Model to Simulate the Development of Carbonate Platforms.” Geology 20 (3): 235–38. https://doi.org/10.1130/0091-7613(1992)020<0235:CACMTS>2.3.CO;2.
Burgess, P. M., V. P. Wright, and D. Emery. 2001. “Numerical Forward Modelling of Peritidal Carbonate Parasequence Development: Implications for Outcrop Interpretation.” Basin Research 13 (1): 1–16. https://doi.org/10.1046/j.1365-2117.2001.00130.x.
Burgess, Peter M. 2013. “CarboCAT: A Cellular Automata Model of Heterogeneous Carbonate Strata.” Computers & Geosciences 53: 129–40.
Burgess, Peter M. 2008. “The Nature of Shallow-Water Carbonate Lithofacies Thickness Distributions.” Geology 36 (3): 235–38. https://doi.org/10.1130/G243326A.1.
Burgess, Peter M. 2016. “Identifying Ordered Strata: Evidence, Methods, and Meaning.” Journal of Sedimentary Research 86 (3): 148–67. https://doi.org/10.2110/jsr.2016.10.
Burgess, Peter M., and David J. Emery. 2004. “Sensitive Dependence, Divergence and Unpredictable Behaviour in a Stratigraphic Forward Model of a Carbonate System.” In Geological Prior Information: Informing Science and Engineering, edited by Andrew Curtis and Rachel Wood, vol. 239. Geological Society of London. https://doi.org/10.1144/GSL.SP.2004.239.01.06.
Burgess, Peter M., Isabella Masiero, Stephan C. Toby, and Robert A. Duller. 2019. “A Big Fan of Signals? Exploring Autogenic and Allogenic Process and Product in a Numerical Stratigraphic Forward Model of Submarine-Fan Development.” Journal of Sedimentary Research 89 (1): 1–12. https://doi.org/10.2110/jsr.2019.3.
Crucifix, Michel. 2023. Palinsol: Insolation for Palaeoclimate Studies. https://doi.org/10.32614/CRAN.package.palinsol.
Curtis, Andrew, Hugo Bloem, Rachel Wood, et al. 2025. “Natural Sampling and Aliasing of Marine Geochemical Signals.” Scientific Reports 15 (1): 760. https://doi.org/10.1038/s41598-024-84871-6.
Danisch, Simon, and Julius Krumbiegel. 2021. “Makie.jl: Flexible High-Performance Data Visualization for Julia.” Journal of Open Source Software 6 (65): 3349. https://doi.org/10.21105/joss.03349.
Demicco, Robert V. 1998. “CYCOPATH 2D—a Two-Dimensional, Forward Model of Cyclic Sedimentation on Carbonate Platforms.” Computers & Geosciences 24 (5): 405–23. https://doi.org/10.1016/S0098-3004(98)00024-7.
Ding, Xuesong, Tristan Salles, Nicolas Flament, and Patrice Rey. 2019. “Quantitative Stratigraphic Analysis in a Source-to-Sink Numerical Framework.” Geoscientific Model Development 12 (6): 2571–85. https://doi.org/10.5194/gmd-12-2571-2019.
Dormann, Sabine, Andreas Deutsch, and Anna T. Lawniczak. 2001. “Fourier Analysis of Turing-Like Pattern Formation in Cellular Automaton Models.” Future Generation Computer Systems 17 (7): 901–9. https://doi.org/https://doi.org/10.1016/S0167-739X(00)00068-6.
Drummond, Carl N., and Pamela J. Dugan. 1999. “Self-Organizing Models of Shallow-Water Carbonate Accumulation.” Journal of Sedimentary Research 69 (4): 939–46. https://doi.org/10.2110/jsr.69.939.
Dyer, Blake, Adam C. Maloof, Sam J. Purkis, and Paul M. (Mitch) Harris. 2018. “Quantifying the Relationship Between Water Depth and Carbonate Facies.” Sedimentary Geology 373 (October): 1–10. https://doi.org/10.1016/j.sedgeo.2018.05.011.
Evans, Kellie M. 1996. “Larger Than Life: It’s so Nonlinear.” {PhD}, Uniersity of Wisconsin-Madison. https://www.csun.edu/~kme52026/thesis.html.
Falivene, Oriol, Alessandro Frascati, Michele Bolla Pittaluga, and John Martin. 2019. “Three-Dimensional Reduced-Complexity Simulation of Fluvio-Deltaic Clastic Stratigraphy.” Journal of Sedimentary Research 89 (1): 46–65. https://doi.org/10.2110/jsr.2018.73.
Fathiyah Jamaludin, Siti Nur. 2025. “Quantitative Geo-History Analysis of the Luconia-Balingian Provinces, Malaysia with Emphasis on the Tectonic Subsidence Signatures.” Marine and Petroleum Geology 173: 107224. https://doi.org/https://doi.org/10.1016/j.marpetgeo.2024.107224.
Geyman, Emily C., Adam C. Maloof, and Blake Dyer. 2021. “How Is Sea Level Change Encoded in Carbonate Stratigraphy?” Earth and Planetary Science Letters 560: 116790. https://doi.org/10.1016/j.epsl.2021.116790.
Granjeon, D., and Phillip Joseph. 1999. “Concepts and Applications of A 3-D Multiple Lithology, Diffusive Model in Stratigraphic Modeling.” In Numerical Experiments in Stratigraphy: Recent Advances in Stratigraphic and Sedimentologic Computer Simulations, edited by J. W. Harbaugh, W. L. Watney, E. C. Rankey, R. Slingerland, R. H Goldstein, and E. K. Franseen, vol. 62. SEPM (Society for Sedimentary Geology). https://doi.org/10.2110/pec.99.62.0197.
Hannisdal, Bjarte. 2006. “Phenotypic Evolution in the Fossil Record: Numerical Experiments.” The Journal of Geology 114 (2): 133–53. https://doi.org/10.1086/499569.
Henglai, Puntira, Rattanaporn Fongngern, and Arthur Saller. 2024. “The Growth and Demise of a Middle Miocene Carbonate Platform in Central Luconia, Offshore Malaysia.” Marine and Petroleum Geology 163: 106763. https://doi.org/https://doi.org/10.1016/j.marpetgeo.2024.106763.
Hidding, Johan. 2023. “Entangled, a Bidirectional System for Sustainable Literate Programming.” 2023 IEEE 19th International Conference on e-Science (e-Science), 1–9. https://doi.org/10.1109/e-Science58273.2023.10254816.
Hidding, Johan, Emilia Jarochowska, Xianyi Liu, Peter Burgess, Niklas Hohmann, and Hanno Spreeuw. 2025. CarboKitten.jl. V. v0.5.2. Zenodo, released June. https://doi.org/10.5281/zenodo.15742533.
Hidding, Johannes. 2026. MindTheGap-ERC/CAExplorer.jl: CAExplorer V0.1.0. V. v0.1.0. Zenodo, released March. https://doi.org/10.5281/zenodo.18925531.
Hill, Jon, Daniel Tetzlaff, Andrew Curtis, and Rachel Wood. 2009. “Modeling Shallow Marine Carbonate Depositional Systems.” Computers & Geosciences 35 (9): 1862–74. https://doi.org/10.1016/j.cageo.2008.12.006.
Hobley, D. E. J., J. M. Adams, S. S. Nudurupati, et al. 2017. “Creative Computing with Landlab: An Open-Source Toolkit for Building, Coupling, and Exploring Two-Dimensional Numerical Models of Earth-Surface Dynamics.” Earth Surface Dynamics 5 (1): 21–46. https://doi.org/10.5194/esurf-5-21-2017.
Hohmann, Niklas, and Emilia Jarochowska. 2025. “StratPal: An r Package for Creating Stratigraphic Paleobiology Modelling Pipelines.” Methods in Ecology and Evolution 16 (4): 678–86. https://doi.org/https://doi.org/10.1111/2041-210X.14507.
Hohmann, Niklas, Joël R. Koelewijn, Peter Burgess, and Emilia Jarochowska. 2024. “Identification of the Mode of Evolution in Incomplete Carbonate Successions.” BMC Ecology and Evolution 24 (1): 113. https://doi.org/10.1186/s12862-024-02287-2.
Holland, Steven M. 2000. “The Quality of the Fossil Record: A Sequence Stratigraphic Perspective.” Paleobiology 26 (S4): 148–68. https://doi.org/10.1017/S0094837300026919.
Holland, Steven M., and Mark E. Patzkowsky. 1999. “Models for Simulating the Fossil Record.” Geology 27 (6): 491–94. https://doi.org/10.1130/0091-7613(1999)027<0491:MFSTFR>2.3.CO;2.
Holland, Steven M., and Mark E. Patzkowsky. 2002. “Stratigraphic Variation in the Timing of First and Last Occurrences.” PALAIOS 17 (2): 134–46. https://doi.org/10.1669/0883-1351(2002)017<0134:SVITTO>2.0.CO;2.
Husinec, Antun, J. Fred Read, and David B. Kemp. 2023. “Orbital Forcing of Upper Jurassic (Tithonian) Shallow-Water Carbonates, Tethyan Adriatic Platform, Croatia Evaluated Using Synthetic Vs. Real Data Sets.” Palaeogeography, Palaeoclimatology, Palaeoecology 622: 111617. https://doi.org/10.1016/j.palaeo.2023.111617.
Hutton, Eric W. H., and James P. M. Syvitski. 2008. “Sedflux 2.0: An Advanced Process-Response Model That Generates Three-Dimensional Stratigraphy.” Computers & Geosciences, Predictive Modeling in Sediment Transport and Stratigraphy, vol. 34 (10): 1319–37. https://doi.org/10.1016/j.cageo.2008.02.013.
James, Scott C, Craig A Jones, Matthew D Grace, and Jesse D Roberts. 2010. “Advances in Sediment Transport Modelling.” Journal of Hydraulic Research 48 (6): 754–63.
Jean Borgomano, Cyprien Lanteaume, Cyprien Lanteaume, Philippe Léonide, François Fournier, Lucien F. Montaggioni, and Jean-Pierre Masse. 2020. “Quantitative Carbonate Sequence Stratigraphy: Insights from Stratigraphic Forward Models.” AAPG Bulletin 104 (5): 1115–42. https://doi.org/10.1306/11111917396.
Kaufman, P., J. P. Grotzinger, and D. S. McCormick. 1991. “Depth-Dependent Diffusion Algorithm for Simulation of Sedimentation in Shallow Marine Depositional Systems.” Bulletin (Kansas Geological Survey) 233: 489–508. https://doi.org/10.17161/kgsbulletin.no.233.20474.
Kemp, David B., Wesley T. Fraser, and Kentaro Izumi. 2018. “Stratigraphic Completeness and Resolution in an Ancient Mudrock Succession.” Sedimentology 65 (6): 1875–90. https://doi.org/10.1111/sed.12450.
Kemp, David B., and Saskia M. Van Manen. 2019. “Metre-Scale Cycles in Shallow Water Carbonate Successions: Milankovitch and Stochastic Origins.” Sedimentology 66 (7): 2590–604. https://doi.org/10.1111/sed.12609.
Kemp, David B., Saskia M. Van Manen, David A. Pollitt, and Peter M. Burgess. 2016. “Investigating the Preservation of Orbital Forcing in Peritidal Carbonates.” Sedimentology 63 (6): 1701–18. https://doi.org/10.1111/sed.12282.
Kenter, Jeroen A. M. 1990. “Carbonate Platform Flanks: Slope Angle and Sediment Fabric.” Sedimentology 37 (5): 777–94. https://doi.org/https://doi.org/10.1111/j.1365-3091.1990.tb01825.x.
Knuth, D. E. 1984. “Literate Programming.” The Computer Journal 27 (2): 97–111. https://doi.org/10.1093/comjnl/27.2.97.
Laskar, Jacques. 2004. “A Long-Term Numerical Solution for the Insolation Quantities of the Earth.” Astronomy & Astrophysics 428: 261–85. https://doi.org/10.1051/0004-6361:20041335.
Lisiecki, Lorraine E., and Maureen E. Raymo. 2005. “A Pliocene-Pleistocene Stack of 57 Globally Distributed Benthic \(\delta\)18O Records.” Paleoceanography 20 (1). https://doi.org/10.1029/2004PA001071.
Liu, Jianliang, and Keyu Liu. 2021. “Estimating Stratal Completeness of Carbonate Deposition via Process-Based Stratigraphic Forward Modeling.” Science China Earth Sciences 64 (2): 253–59. https://doi.org/10.1007/s11430-020-9660-8.
Liu, Jinlong, Jody M. Webster, Tristan Salles, et al. 2022. “The Formation of Atolls: New Insights From Numerical Simulations.” Journal of Geophysical Research: Earth Surface 127 (8): e2022JF006812. https://doi.org/10.1029/2022JF006812.
Lopez-Gamundi, Cecilia, Brian B. Barnes, Christian Betzler, et al. 2025. “The Sediment Budget of Great Bahama Bank—Earth’s Largest Modern Carbonate Platform.” Geology 53 (9): 748–52. https://doi.org/10.1130/G52850.1.
Masiero, Isabella, Peter Burgess, Cathy Hollis, et al. 2021. “Syn-Rift Carbonate Platforms in Space and Time: Testing and Refining Conceptual Models Using Stratigraphic and Seismic Numerical Forward Modelling.” In Seismic Characterization of Carbonate Platforms and Reservoirs, edited by J. Hendry, P. Burgess, D. Hunt, X. Janson, and V. Zampetti. Geological Society of London. https://doi.org/10.1144/SP509-2019-217.
Masiero, Isabella, Estanislao Kozlowski, Georgios Antonatos, Haiwei Xi, and Peter Burgess. 2020. “Numerical Stratigraphic Forward Models as Conceptual Knowledge Repositories and Experimental Tools: An Example Using a New Enhanced Version of CarboCAT.” Computers & Geosciences 138 (May): 104453. https://doi.org/10.1016/j.cageo.2020.104453.
Miller, Kenneth G., Michelle A. Kominz, James V. Browning, et al. 2005. “The Phanerozoic Record of Global Sea-Level Change.” Science 310 (5752): 1293–98. https://doi.org/10.1126/science.1116412.
Mitchell, S. F., C. R. C. Paul, and A. S. Gale. 1996. “Carbon Isotopes and Sequence Stratigraphy.” Geological Society, London, Special Publications 104 (1): 11–24. https://doi.org/10.1144/GSL.SP.1996.104.01.02.
Myrow, Paul M., and John P. Grotzinger. 2000. “Chemostratigraphic Proxy Records: Forward Modeling the Effects of Unconformities, Variable Sediment Accumulation Rates, and Sampling-Interval Bias.” In Carbonate Sedimentation and Diagenesis in the Evolving Precambrian World, edited by John P. Grotzinger and Noel P. James, vol. 67. SEPM Society for Sedimentary Geology. https://doi.org/10-2110/pec.00.67.
Paola, Chris, Paul L Heller, and Charles L Angevine. 1992. “The Large-Scale Dynamics of Grain-Size Variation in Alluvial Basins, 1: Theory.” Basin Research 4 (2): 73–90.
Pastier, A.-M., L. Husson, K. Pedoja, et al. 2019. “Genesis and Architecture of Sequences of Quaternary Coral Reef Terraces: Insights From Numerical Models.” Geochemistry, Geophysics, Geosystems 20 (8): 4248–72. https://doi.org/10.1029/2019GC008239.
Paterson, Richard J., Fiona F. Whitaker, Gareth D. Jones, Peter L. Smart, David Waltham, and Graham Felce. 2006. “Accommodation and Sedimentary Architecture of Isolated Icehouse Carbonate Platforms: Insights from Forward Modeling with CARB3D+.” Journal of Sedimentary Research 76 (10): 1162–82. https://doi.org/10.2110/jsr.2006.113.
Purkis, Sam J., Johan van de Koppel, and Peter M. Burgess. 2016. “Spatial Self-Organization in Carbonate Depositional Environments.” SEPM Special Publications 106: 53–66. https://doi.org/10.2110/sepmsp.106.02.
Salles, T. 2016. “Badlands: A Parallel Basin and Landscape Dynamics Model.” SoftwareX 5 (January): 195–202. https://doi.org/10.1016/j.softx.2016.08.005.
Salles, Tristan, Xuesong Ding, and Gilles Brocard. 2018. “pyBadlands: A Framework to Simulate Sediment Transport, Landscape Dynamics and Basin Stratigraphic Evolution Through Space and Time.” PLOS ONE 13 (4): e0195557. https://doi.org/10.1371/journal.pone.0195557.
Salles, Tristan, Jodie Pall, Jody M. Webster, and Belinda Dechnik. 2018. “Exploring Coral Reef Responses to Millennial-Scale Climatic Forcings: Insights from the 1-D Numerical Tool pyReef-Core V1.0.” Geoscientific Model Development 11 (6): 2093–110. https://doi.org/10.5194/gmd-11-2093-2018.
Schlager, Wolfgang, and Oren Camber. 1986. “Submarine Slope Angles, Drowning Unconformities, and Self-Erosion of Limestone Escarpments.” Geology 14 (9): 762–65. https://doi.org/10.1130/0091-7613(1986)14<762:SSADUA>2.0.CO;2.
Schlager, Wolfgang, and Georg Warrlich. 2009. “Record of Sea‐level Fall in Tropical Carbonates.” Basin Research 21 (2): 209–24. https://doi.org/10.1111/j.1365-2117.2008.00383.x.
Strobel, John, Robert Cannon, G. St. Christopher, C. St. Kendall, Gautam Biswas, and James Bezdek. 1989. “Interactive (SEDPAK) Simulation of Clastic and Carbonate Sediments in Shelf to Basin Settings.” Computers & Geosciences 15 (8): 1279–90. https://doi.org/10.1016/0098-3004(89)90092-7.
Sultana, Daniel, Peter Burgess, and Dan Bosence. 2022. “How Do Carbonate Factories Influence Carbonate Platform Morphology? Exploring Production-Transport Interactions with Numerical Forward Modelling.” Sedimentology 69 (1): 372–93. https://doi.org/10.1111/sed.12943.
Sylvester, Zoltán, Kyle M. Straub, and Jacob A. Covault. 2024. “Stratigraphy in Space and Time: A Reproducible Approach to Analysis and Visualization.” Earth-Science Reviews 250 (March): 104706. https://doi.org/10.1016/j.earscirev.2024.104706.
Tetzlaff, Daniel. 2023. Stratigraphic Forward Modeling Software Package for Research and Education. arXiv. https://doi.org/10.48550/arXiv.2302.05272.
Warrlich, G. M. D., D. A. Waltham, and D. W. J. Bosence. 2002. “Quantifying the Sequence Stratigraphy and Drowning Mechanisms of Atolls Using a New 3‐D Forward Stratigraphic Modelling Program (CARBONATE 3D).” Basin Research 14 (3): 379–400. https://doi.org/10.1046/j.1365-2117.2002.00181.x.
Warrlich, Georg, Dan Bosence, Dave Waltham, Charlie Wood, Amy Boylan, and Beatriz Badenas. 2008. “3D Stratigraphic Forward Modelling for Analysis and Prediction of Carbonate Platform Stratigraphies in Exploration and Production.” Marine and Petroleum Geology 25 (1): 35–58. https://doi.org/https://doi.org/10.1016/j.marpetgeo.2007.04.005.
Warrlich, Georg-Mathis Dieter. 2000. “3D Computer Forward Modelling of Carbonate Platform Evolution.” PhD thesis, Royal Holloway, University of London.
Weij, Rieneke, John J. G. Reijmer, Gregor P. Eberli, and Peter K. Swart. 2019. “The Limited Link Between Accommodation Space, Sediment Thickness, and Inner Platform Facies Distribution (Holocene–Pleistocene, Bahamas).” The Depositional Record 5 (3): 400–420. https://doi.org/10.1002/dep2.50.
Wild, Thomas B., Daniel P. Loucks, and George W. Annandale. 2019. “SedSim: A River Basin Simulation Screening Model for Reservoir Management of Sediment, Water, and Hydropower.” Journal of Open Research Software 7 (1). https://doi.org/10.5334/jors.261.
Xi, Haiwei, and Peter M. Burgess. 2022. “The Stratigraphic Significance of Self-Organization: Exploring How Autogenic Processes Can Generate Cyclical Carbonate Platform Strata.” Sedimentology 69 (4): 1769–88. https://doi.org/10.1111/sed.12974.
Zimmt, Joshua B., Steven M. Holland, Seth Finnegan, and Charles R. Marshall. 2021. “Recognizing Pulses of Extinction from Clusters of Last Occurrences.” Palaeontology 64 (1): 1–20. https://doi.org/10.1111/pala.12505.