Note

This page was generated from a jupyter notebook.

The Carbonate Producer component

Overview

This notebook demonstrates the CarbonateProducer Landlab component. The component computes the creation of carbonate rocks, such as carbonate platforms and coral reefs, given a particular bathymetry. The component can be used either to calculate the rate of carbonate production (in terms of a vertical deposition rate), or to progressively add to a field representing carbonate thickness. The model does not distinguish among different types of carbonate material, or species of carbonate-producing organism.

Theory

Carbonate production rate

The CarbonateProducer uses the mathematical model of Bosscher and Schlager (1992), which represents the local carbonate growth rate, \(G(x,y,t)\) (thickness per time) as a function of local water depth. The carbonate production rate is calculated as

\[G(x,y,t) = G_m \tanh \left( I_0 e^{-kd} / I_k \right)\]

where:

• \(G_m\) is the maximum production rate

• \(I_0\) is the surface light intensity

• \(I_k\) is the saturating light intensity

• \(d\) is water depth

• \(k\) is the light extinction coefficient

By default the production rate is zero where \(d<0\), but as noted below the user can invoke an option that allows for carbonate production within the tidal range.

Bosscher and Schlager (1992) suggest plausible values or ranges for these parameters as follows: \(G_m\) = 0.010 to 0.015 m/y, \(I_0\) = 2,000 to 2,250 micro Einsteins per square meter per second in the tropics, \(I_k\) = 50 to 450 micro Einsteins per square meter per second, and \(k\) = 0.04 to 0.16 m\(^{-1}\) (corresponding to a decay depth of 6.25 to 16 m).

Smoothing near sea level using tidal range

The default form of the model involves a mathematical discontinuity at zero water depth. The user has the option of smoothing out this discontinuity by passing a positive value for the tidal_range parameter. If this parameter is given, the growth function is modified as follows:

\[G' = H(d) G\]

where \(G\) is the growth rate calculated by the growth function shown above, and \(H(d)\) is a smoothed Heaviside step function of local water depth that is defined as:

\[\begin{split}H(d) = \begin{cases} \frac{1}{1 + \exp(-2 / R_t d)} & \text{if } d > -2R_t \\ 0 & \text{otherwise.} \end{cases}\end{split}\]

Essentially, the \(H(d)\) function allows a limited amount of growth above mean sea level, while reducing the growth rate somewhat within the tidal zone.

Numerical methods

The component’s calc_carbonate_production_rate method can be used to return the current rate of production given the water depths (calculated by subtracting the topographic__elevation node field from the sea_level__elevation grid field). In this case, no numerical methods are needed. This approach might be useful, for example, when modeling simultaneous carbonate and siliciclastic sedimentation, and the user wishes to generate depositional layers that contain some fractional amount of both.

Alternatively, the user can calculate the accumulation of carbonate thickness (in node field carbonate_thickness) by calling either produce_carbonate or run_one_step (the latter simply calls the former). In this case, simple forward Euler differencing is used to add to carbonate thickness, \(C\), via

\[C_i^{k+1} = C_i^k + G_i^k \Delta t\]

where \(i\) refers to node number and \(k\) to time-step number, and \(\Delta t\) is the duration of the time step (passed as an argument to produce_carbonate or run_one_step).

import numpy as np

from landlab import RasterModelGrid
from landlab.components import CarbonateProducer
from landlab.plot import plot_layers

Information about the component

Passing the class name to the help function returns descriptions of the various methods and parameters:

help(CarbonateProducer)

Examples

Example 1: carbonate growth on a rising continental margin under sinusoidal sea-level variation

In this example, we consider a sloping ramp that rises tectonically while sea level oscillates.

# Parameters and settings
nrows = 3  # number of grid rows
ncols = 101  # number of grid columns
dx = 1000.0  # node spacing, m
sl_range = 120.0  # range of sea-level variation, m
sl_period = 40000.0  # period of sea-level variation, y
run_duration = 200000.0  # duration of run, y
dt = 100.0  # time-step size, y
initial_shoreline_pos = 25000.0  # starting position of the shoreline, m
topo_slope = 0.01  # initial slope of the topography/bathymetry, m/m
uplift_rate = 0.001  # rate of tectonic uplift, m/y
plot_ylims = [-1000.0, 1000.0]

# Derived parameters
num_steps = int(run_duration / dt)
sin_fac = 2.0 * np.pi / sl_period  # factor for sine calculation of sea-level
middle_row = np.arange(ncols, 2 * ncols, dtype=int)

# Grid and fields
#
# Create a grid object
grid = RasterModelGrid((nrows, ncols), xy_spacing=dx)

# Create sea level field (a scalar, a.k.a. a "grid field")
sea_level = grid.add_field("sea_level__elevation", 0.0, at="grid")

# Create elevation field and initialize as a sloping ramp
bedrock_elevation = topo_slope * (initial_shoreline_pos - grid.x_of_node)
elev = grid.add_field("topographic__elevation", bedrock_elevation.copy(), at="node")
# elev[:] = topo_slope * (initial_shoreline_pos - grid.x_of_node)

# Remember IDs of middle row of nodes, for plotting
middle_row = np.arange(ncols, 2 * ncols, dtype=int)

plot_layers(
    bedrock_elevation[middle_row],
    x=grid.x_of_node[middle_row],
    sea_level=grid.at_grid["sea_level__elevation"],
    x_label="Distance (km)",
    y_label="Elevation (m)",
    title="Starting condition",
    legend_location="upper right",
)

# Instantiate component
cp = CarbonateProducer(grid)

# RUN
for i in range(num_steps):
    cp.sea_level = sl_range * np.sin(sin_fac * i * dt)
    cp.produce_carbonate(dt)
    elev[:] += uplift_rate * dt

plot_layers(
    [
        elev[middle_row] - grid.at_node["carbonate_thickness"][middle_row],
        elev[middle_row],
    ],
    x=grid.x_of_node[middle_row],
    sea_level=grid.at_grid["sea_level__elevation"],
    color_layer="Blues",
    x_label="Distance (km)",
    y_label="Elevation (m)",
    title="Carbonate production",
    legend_location="upper right",
)

Example 2: tracking stratigraphy

Here we repeat the same example, except this time we use Landlab’s MaterialLayers class to track stratigraphy.

# Track carbonate strata in time bundles of the below duration:
layer_time_interval = 20000.0

# Derived parameters and miscellaneous
next_layer_time = layer_time_interval
time_period_index = 0
time_period = "0 to " + str(int(layer_time_interval) // 1000) + " ky"

# Grid and fields
grid = RasterModelGrid((nrows, ncols), xy_spacing=dx)
sea_level = grid.add_field("sea_level__elevation", 0.0, at="grid")
base_elev = grid.add_zeros("basement__elevation", at="node")
base_elev[:] = topo_slope * (initial_shoreline_pos - grid.x_of_node)
elev = grid.add_zeros("topographic__elevation", at="node")
middle_row = np.arange(ncols, 2 * ncols, dtype=int)
middle_row_cells = np.arange(0, ncols - 2, dtype=int)
carbo_thickness = grid.add_zeros("carbonate_thickness", at="node")
prior_carbo_thickness = np.zeros(grid.number_of_nodes)

# Instantiate component
cp = CarbonateProducer(grid)

# RUN
for i in range(num_steps):
    cp.sea_level = sl_range * np.sin(sin_fac * i * dt)
    cp.produce_carbonate(dt)
    base_elev[:] += uplift_rate * dt
    elev[:] = base_elev + carbo_thickness
    if (i + 1) * dt >= next_layer_time:
        time_period_index += 1
        next_layer_time += layer_time_interval
    added_thickness = np.maximum(
        carbo_thickness - prior_carbo_thickness, 0.00001
    )  # force a tiny bit of depo to keep layers consistent
    prior_carbo_thickness[:] = carbo_thickness
    # grid.material_layers.add(added_thickness[grid.node_at_cell], age=time_period_index)
    grid.event_layers.add(added_thickness[grid.node_at_cell], age=time_period_index)

First get the layers we want to plot. In this case, plot the bottom and top layers as well as layers that correspond to sea level high stands. For the sinusoidal sea level curve we used, high stands occur every 400 time steps, with the first one being at time step 100.

layers = (
    np.vstack(
        [
            grid.event_layers.z[0],
            grid.event_layers.z[100::400],
            grid.event_layers.z[-1],
        ]
    )
    + grid.at_node["basement__elevation"][grid.node_at_cell]
)

plot_layers(
    layers,
    x=grid.x_of_node[grid.node_at_cell],
    sea_level=grid.at_grid["sea_level__elevation"],
    color_layer="Oranges_r",
    legend_location="upper right",
    x_label="Distance (km)",
    y_label="Elevation (m)",
    title="Carbonates colored by age of deposition (darkest = oldest)",
)

References

Bosscher, H., & Schlager, W. (1992). Computer simulation of reef growth. Sedimentology, 39(3), 503-512.

Galewsky, J. (1998). The dynamics of foreland basin carbonate platforms: tectonic and eustatic controls. Basin Research, 10(4), 409-416.

---

Generated by nbsphinx from a Jupyter notebook.