Package {BasinmodlR}


Type: Package
Title: Burial and Thermal History Modelling of Sedimentary Basins
Version: 0.1.0
Maintainer: Michiel Arts <michiel.arts@stratigraphy.eu>
Depends: R (≥ 3.5.0)
Imports: graphics, grDevices, stats, utils
Suggests: testthat (≥ 3.1.9)
Description: Simulates the burial, compaction and thermal history of a sedimentary basin in one dimension, for geologists and stratigraphers. A basin is described by one row per stratigraphic unit and one row per erosion event, from which deposition phases, hiatuses, layer stacks, maximum-burial windows and grain (solid) thicknesses are derived. Compaction follows an exponential porosity-depth law and is irreversible, so uplifted layers keep the porosity acquired at their deepest burial. Temperatures follow the steady-state conductive geotherm of a layered column with radiogenic heat production. The methods follow the textbook treatment of Allen and Allen (2013, ISBN:978-0-470-67377-5). Functions are included to draw a composite figure of burial, temperature and hydrocarbon generation windows and a petroleum system event chart.
License: GPL-3
URL: https://github.com/stratigraphy/BasinmodlR
BugReports: https://github.com/stratigraphy/BasinmodlR/issues
Encoding: UTF-8
Language: en-GB
LazyData: true
RoxygenNote: 7.3.3
NeedsCompilation: no
Config/testthat/edition: 3
Packaged: 2026-10-28 13:34:18 UTC; Administrator
Author: Michiel Arts ORCID iD [aut, cre, cph]
Repository: CRAN
Date/Publication: 2026-11-13 17:10:33 UTC

Burial and Thermal History Modelling of Sedimentary Basins

Description

BasinmodlR is an R package for one-dimensional (1D) burial, compaction and thermal history modelling of sedimentary basins. It is aimed at geologists and stratigraphers who want to reconstruct how deep a stratigraphic unit was buried, how much it was compacted and how hot it became through geological time, starting from nothing more than a stratigraphic column. The whole basin history is described by two small tables: one row per stratigraphic unit and one row per erosion event. From these tables the package derives the deposition phases, hiatuses, layer stacks, the maximum-burial depth of each unit and its grain (solid) thickness. Compaction follows an exponential porosity-depth law and is treated as irreversible, so that uplifted layers keep the porosity they acquired at their deepest burial. Temperatures are the steady-state conductive geotherm of the layered column including radiogenic heat production. The package includes functions to draw a composite figure of the burial path, temperature and hydrocarbon generation windows and to draw a petroleum system event chart. The methods follow the textbook treatment of Allen and Allen (2013).

Details

Package: 'BasinmodlR'

Type: R package

Version: 0.1.0 (4th quarter 2026)

License: GPL-3

Describing a basin

A basin is a list with five elements (see brabant, condroz and ardennes for worked examples):

units

one row per stratigraphic unit, oldest first, with the deposition interval (start and end in Ma, start being the older age), the thickness at maximum burial (thickness, km), the depositional porosity (phi0), the porosity-depth coefficient (ck, 1/km), the thermal conductivity of the grains (K, W/m/K), the radiogenic heat production (A, uW/m3) and optionally the grain density (rho, kg/m3).

erosions

one row per erosion event, with start, end (Ma) and the amount of section removed (amount, km). Erosion is stripped off the top of the column and may cut through several units.

water_depth

water depth (m) as (age, value) knots.

heat_flow

heat flow into the base of the column (mW/m2) as (age, value) knots.

surface_temp

temperature at the sediment-water interface (degrees Celsius) as (age, value) knots.

Deposition phases follow from the unit intervals, erosion phases from the erosion table, and every remaining gap in the timeline automatically becomes a hiatus. Nothing has to be repeated for each phase: layer stacks, layer identifiers, maximum-burial windows and the amount of section previously removed are all worked out by run_basin_model.

How the model works

Porosity decreases exponentially with depth below the sediment-water interface: \phi(z) = \phi_0 e^{-c z}. For every unit the model first determines the deepest position that the unit ever occupies and from the thickness at that position it computes the volume of grains (the solid thickness), which does not change during burial. At every time step the thickness of a layer whose top is at depth z_1 is then found by solving

h = h_{solid} + (\phi_0/c)(e^{-c z_1} - e^{-c (z_1 + h)}),

which is iterated to convergence. Compaction is irreversible, so the thickness that is used is the smaller of this solution and the thickness the layer had at its deepest burial; an exhumed layer keeps its porosity. During deposition grain material is added at a constant rate, so a unit reaches exactly its maximum-burial thickness at the end of its own deposition phase and hands over to the next phase without a step.

Temperatures are the steady-state conductive solution for the layered column. The bulk conductivity of a layer is the porosity-weighted mean of the conductivity of the grains and of the pore water. The surface heat flow is the basal heat flow plus the radiogenic production of the whole column, and the heat flow entering each layer is reduced by the production of the layers above it, so that within a layer of thickness h

T_{base} = T_{top} + q h / K - A h^2 / (2 K).

Main functions

run_basin_model

runs the model and returns a table with one row for every layer at every age plus a temperature raster on a common age-depth grid.

plot_basin_model

composite figure: geological periods, event bar, temperature raster with burial paths and isotherms, colour key and hydrocarbon generation-window panels. Also available as plot(<basin_model>).

plot_event_chart

petroleum system event chart.

Units

Input thicknesses and erosion amounts are in km, porosity coefficients in 1/km, radiogenic heat production in uW/m3, heat flow in mW/m2 and water depth in m. The results are returned in metres, degrees Celsius, W/m/K and Ma.

Note

If you want to use this package for publication or research purposes, please cite:

Arts, M.C.M. (2026). BasinmodlR: Burial and Thermal History Modelling of Sedimentary Basins. https://CRAN.R-project.org/package=BasinmodlR

Author(s)

Michiel Arts

Maintainer: Michiel Arts michiel.arts@stratigraphy.eu

References

The 'BasinmodlR' package builds upon the textbook treatment of basin analysis by Allen and Allen (2013). The following list of articles and books is relevant for the 'BasinmodlR' R package and its functions. Individual references are also cited in the descriptions of the functions when relevant. The list can be grouped in three subjects: (1) basin analysis, compaction and backstripping, (2) heat flow and thermal modelling and petroleum systems and (3) the example basins.

#1. Basin analysis, compaction and backstripping

Allen, P.A., and Allen, J.R. (2013). Basin Analysis: Principles and Application to Petroleum Play Assessment, 3rd ed. Wiley-Blackwell, Chichester, 632 p. ISBN 978-0-470-67377-5.

Sclater, J.G., and Christie, P.A.F. (1980). Continental stretching: an explanation of the post-mid-Cretaceous subsidence of the central North Sea Basin. Journal of Geophysical Research, 85, 3711–3739. doi:10.1029/JB085iB07p03711

Steckler, M.S., and Watts, A.B. (1978). Subsidence of the Atlantic-type continental margin off New York. Earth and Planetary Science Letters, 41, 1–13. doi:10.1016/0012-821X(78)90036-5

#2. Thermal modelling and petroleum systems

Hantschel, T., and Kauerauf, A.I. (2009). Fundamentals of Basin and Petroleum Systems Modeling. Springer, Berlin, 476 p. doi:10.1007/978-3-540-72318-9

Magoon, L.B., and Dow, W.G. (1994). The petroleum system. In: Magoon, L.B. and Dow, W.G. (eds), The Petroleum System - From Source to Trap. American Association of Petroleum Geologists Memoir 60, 3–24. doi:10.1306/M60585C1

Tissot, B.P., and Welte, D.H. (1984). Petroleum Formation and Occurrence, 2nd ed. Springer, Berlin, 699 p. doi:10.1007/978-3-642-87813-8

Turcotte, D.L., and Schubert, G. (2014). Geodynamics, 3rd ed. Cambridge University Press, Cambridge.

#3. Example basins

Arts, M.C.M., Djouder, H., Devleeschouwer, X., and Da Silva, A.-C. (2026). Solid bitumen in the Silurian Bonne Esperance Formation at Huy, Belgium: petroleum system elements reconstructed from integrated maturity indicators and 1D petroleum system modelling. Manuscript submitted to Geologica Belgica.

Gradstein, F.M., Ogg, J.G., Schmitz, M.D., and Ogg, G.M. (eds) (2020). Geologic Time Scale 2020. Elsevier, Amsterdam, 1357 p. doi:10.1016/C2020-1-02369-3

Li, X., Hu, Y., Guo, J., Lan, J., Lin, Q., Bao, X., Yuan, S., Wei, M., Li, Z., Man, K., Yin, Z., Han, J., Zhang, J., Zhu, C., Zhao, Z., Liu, Y., Yang, J., and Nie, J. (2022). A high-resolution climate simulation dataset for the past 540 million years. Scientific Data, 9, 371. doi:10.1038/s41597-022-01490-4

Examples

# Burial and thermal history of the Huy section (Condroz inlier, Belgium)
res <- run_basin_model(condroz)
res
head(res$profile)


Example basins for the 'BasinmodlR' package

Description

Three worked basin descriptions for the Lower Palaeozoic of Belgium, each ready to be handed to run_basin_model:

The brabant data set is a highly generalised 1D model of the Brabant Massif (525-0 Ma).

The condroz data set is the 1D model of the Rue Bonne Espérance section at Huy in the Condroz Inlier (520-0 Ma).

The ardennes data set is a highly generalised 1D model of the Cambro-Ordovician inliers of the Ardennes, such as the Stavelot Massif (525-0 Ma).

Usage

brabant

condroz

ardennes

Format

Each data set is a list of seven elements, described under Details.

Details

All three data sets have the same structure: a list with the elements

name

character, the label that is used by print and by the figures.

units

data frame, one row per stratigraphic unit, oldest first, with the columns name (character), start and end (deposition interval in Ma, start being the older age), thickness (km at maximum burial, that is as deposited and before any later erosion of that unit), phi0 (depositional porosity), ck (porosity-depth coefficient, 1/km), K (thermal conductivity of the grains, W/m/K), A (radiogenic heat production, uW/m3) and rho (grain density, kg/m3, carried through but not used in the calculation yet).

erosions

data frame with start, end (Ma) and amount (km of section that is removed from the top of the column).

water_depth

two column matrix with the age (Ma) and the water depth (m).

heat_flow

two column matrix with the age (Ma) and the heat flow into the base of the column (mW/m2).

surface_temp

data frame with the age (Ma) and the temperature at the sediment-water interface (degrees Celsius), at 10 Myr resolution.

source_units

integer, the rows of units that get a generation panel in plot_basin_model.

Ages between the knots of water_depth, heat_flow and surface_temp are filled in by linear interpolation, and the curves are held constant beyond their first and last knot. Gaps in the timeline, such as 480-467 Ma in the Brabant Massif and at Huy, become hiatuses automatically.

Stratigraphic units are grouped where possible, so that one modelled unit can contain several formations or groups. Unit names such as "layer 1" are generalised units without a formal lithostratigraphic name. Thermal conductivity, heat production, depositional porosity and the porosity-depth coefficient are reference values for generalised lithologies (Allen and Allen, 2013), and the basal heat flow was assigned to each time interval according to the prevailing tectonic regime. The heat flow and the eroded thicknesses are therefore scenario assumptions and not independently constrained quantities, and the models should be regarded as scenario-based reconstructions (Hantschel and Kauerauf, 2009). The palaeo-surface temperatures were extracted for the palaeo-position of Belgium from the climate simulations of Li et al. (2022) and the same curve is used for all three basins.

The maximum-burial depth window of every unit, which has to be listed by hand for every phase in layer-by-layer scripts, follows from the units and the erosions and is not part of the input. For the Brabant Massif this gives 9600 m of section at 420 Ma, 8600 m at 300 Ma and 6100 m preserved today.

Source

The stratigraphic framework (unit thickness, depositional environment and the position of hiatuses) was compiled from the published stratigraphy of the Brabant Massif, the Condroz Inlier and the Ardennes and is described in Arts et al. (2026). The surface temperature curve is derived from Li et al. (2022).

References

Allen, P.A., and Allen, J.R. (2013). Basin Analysis: Principles and Application to Petroleum Play Assessment, 3rd ed. Wiley-Blackwell, Chichester, 632 p. ISBN 978-0-470-67377-5.

Arts, M.C.M., Djouder, H., Devleeschouwer, X., and Da Silva, A.-C. (2026). Solid bitumen in the Silurian Bonne Espérance Formation at Huy, Belgium: petroleum system elements reconstructed from integrated maturity indicators and 1D petroleum system modelling. Manuscript submitted to Geologica Belgica.

Hantschel, T., and Kauerauf, A.I. (2009). Fundamentals of Basin and Petroleum Systems Modeling. Springer, Berlin, 476 p. doi:10.1007/978-3-540-72318-9

Li, X., Hu, Y., Guo, J., Lan, J., Lin, Q., Bao, X., Yuan, S., Wei, M., Li, Z., Man, K., Yin, Z., Han, J., Zhang, J., Zhu, C., Zhao, Z., Liu, Y., Yang, J., and Nie, J. (2022). A high-resolution climate simulation dataset for the past 540 million years. Scientific Data, 9, 371. doi:10.1038/s41597-022-01490-4

Examples

# Structure of an example basin
str(brabant, max.level = 1)

# The stratigraphic units of the Brabant Massif
brabant$units

# The erosion events at Huy
condroz$erosions

# The water depth curve of the Ardennes (age in Ma, water depth in m)
ardennes$water_depth


Print and plot methods for basin model results

Description

Methods for objects of class "basin_model", the result of run_basin_model. print gives a short summary of the modelled burial and temperature history and plot draws the composite figure of plot_basin_model.

Usage

## S3 method for class 'basin_model'
print(x, ...)

## S3 method for class 'basin_model'
plot(x, ...)

Arguments

x

An object of class "basin_model".

...

For plot, further arguments passed on to plot_basin_model. Not used by print.

Value

x, invisibly.

Author(s)

Michiel Arts

See Also

run_basin_model, plot_basin_model

Examples

res <- run_basin_model(brabant)
print(res)
plot(res)

Composite burial and temperature figure

Description

Draws the standard figure of a basin model run: a geological period bar, an event bar (deposition, erosion and hiatus phases), a temperature raster with the water column, the burial path of every stratigraphic unit (horizon lines) and labelled isotherms, a colour key, and one generation-window panel for every selected unit.

Usage

plot_basin_model(
  res,
  window = c(95, 135),
  window_units = NULL,
  n_levels = 250,
  scale = c("quantile", "linear"),
  contour_n = 20,
  xmax = NULL,
  tick = 25,
  new_device = FALSE,
  use = c("T_mid", "T_top", "T_base"),
  span = 5,
  verbose = FALSE
)

Arguments

res

An object of class "basin_model", the result of run_basin_model.

window

Numeric vector of length two, the temperature window in degrees Celsius that is used for the generation panels. Default is c(95, 135).

window_units

Which units get a generation panel, given as row numbers of res$units or as unit names. NULL uses res$source_units when the basin description provides it and draws no panel otherwise.

n_levels

Number of colour levels in the temperature raster. Default is 250.

scale

How the colours are spread over the temperature range, either "quantile" (equal areas get equal numbers of colours, the default) or "linear" (equal temperature steps get equal colours).

contour_n

Approximate number of isotherms. Default is 20.

xmax

Oldest age shown in Ma. NULL (default) uses the start of the model.

tick

Spacing of the age axis labels in Myr. Default is 25.

new_device

Open a new graphics window before plotting. Leave this FALSE (default) when writing to a file device or to the plot pane.

use

The temperature that is used for the generation window test, one of "T_mid" (the middle of the layer, the default), "T_top" or "T_base".

span

Half width, in Myr, of the interval on either side of a phase boundary that is used to measure the background rate for the continuity check. Default is 5.

verbose

Print the continuity check table. Default is FALSE.

Details

Burial and temperature panel. The horizontal axis is the age (Ma) running from the oldest age on the left to the present on the right. The vertical axis is the depth below sea level (m). The coloured raster is the temperature of the column, the light blue area at the top is the water column, and every black line is the base of one stratigraphic unit. Where a line moves up the unit has been uplifted and eroded, and where it moves down it has been buried deeper. The thin lines are isotherms.

Generation windows. The lower panels show for a selected unit when it is inside a temperature window (for example the oil window) for the first time. Because organic matter that has been heated cannot go back to an immature state, generation can only continue when the unit reaches a temperature that is higher than any temperature it has had before. A unit is therefore drawn as "gen" (generating) at an age when its temperature lies inside window and is equal to the highest temperature that the unit has reached up to that age (Tissot and Welte, 1984; Allen and Allen, 2013). The default window of 95 to 135 degrees Celsius is a nominal choice and should be adjusted to the kerogen type and the maturity model of the study at hand.

Continuity check. With verbose = TRUE the function prints a table that compares the change in layer depth and temperature over the one million year step that straddles every phase boundary with the background rate of change on either side of it. Ratios close to 1 show that no artificial step is introduced at the boundary between two phases.

The geological periods and their colours follow the International Chronostratigraphic Chart (Gradstein et al., 2020). The colour scale of the temperature raster is a rainbow-like palette (the 'turbo' colour scheme).

Value

res, invisibly. The function is called for its side effect, the figure.

Author(s)

Michiel Arts

References

Allen, P.A., and Allen, J.R. (2013). Basin Analysis: Principles and Application to Petroleum Play Assessment, 3rd ed. Wiley-Blackwell, Chichester, 632 p. ISBN 978-0-470-67377-5.

Gradstein, F.M., Ogg, J.G., Schmitz, M.D., and Ogg, G.M. (eds) (2020). Geologic Time Scale 2020. Elsevier, Amsterdam, 1357 p. doi:10.1016/C2020-1-02369-3

Tissot, B.P., and Welte, D.H. (1984). Petroleum Formation and Occurrence, 2nd ed. Springer, Berlin, 699 p. doi:10.1007/978-3-642-87813-8

See Also

run_basin_model, plot_event_chart

Examples

# Burial and temperature history of the Huy section (Condroz inlier).
# The generation panels are drawn for the units listed in condroz$source_units.
res <- run_basin_model(condroz)
plot_basin_model(res)

# Show the generation window of the Mousty equivalent only, using a
# narrower window, and print the continuity check.
plot_basin_model(res,
                 window = c(100, 130),
                 window_units = "Mousty eq. (unit 2)",
                 verbose = TRUE)

Plot a petroleum system event chart from a basin model result

Description

Draws a petroleum system event chart (Magoon and Dow, 1994) from the output of run_basin_model. The chart shows, against geological time, when the source rocks, reservoir rocks and seal rocks were deposited, when the overburden was deposited, when the section was uplifted and eroded, when traps formed, when hydrocarbons were generated, migrated and accumulated, how long the accumulation has been preserved, and the critical moment of the system.

Usage

plot_event_chart(
  res,
  source_units = NULL,
  reservoir_units = NULL,
  seal_units = NULL,
  trap_formation = NULL,
  generation_window = c(95, 135),
  critical_moment = NULL,
  temperature_column = "T_mid",
  xmax = NULL,
  tick = 25,
  col = "darkred",
  critical_col = "black",
  main = "",
  xlab = "Age (Ma)"
)

Arguments

res

An object returned by run_basin_model. It supplies the model units, the deposition, erosion and hiatus phases, the ages and the temperature of every layer.

source_units

Source-rock units, given as row numbers of res$units or as unit names. Default is NULL (no source rock).

reservoir_units

Reservoir-rock units, given as row numbers or unit names. Default is NULL.

seal_units

Seal-rock units, given as row numbers or unit names. Default is NULL.

trap_formation

Optional data frame with the columns start and end (Ma, start being the older age) that gives the interval(s) of trap formation. Default is NULL.

generation_window

Numeric vector of length two with the temperature window (degrees Celsius) that is used to identify generation, migration and accumulation. Default is c(95, 135).

critical_moment

Age (Ma) of the critical moment of the petroleum system, the time that best represents the generation, migration and accumulation of most of the hydrocarbons (Magoon and Dow, 1994). Default is NULL, in which case neither the critical moment nor the preservation time is drawn.

temperature_column

Name of the temperature column of res$profile that is used, one of "T_mid" (default), "T_top" or "T_base".

xmax

Oldest age shown on the x-axis in Ma. NULL (default) uses the oldest model age.

tick

Spacing between the age labels of the x-axis in Myr. Default is 25.

col

Default colour of the event bars. Default is "darkred".

critical_col

Colour of the critical moment marker. Default is "black".

main

Plot title. Default is an empty string.

xlab

Label of the x-axis. Default is "Age (Ma)".

Details

Source rocks. Source rocks are displayed as separate bars and may be assigned different shades of red. Source-rock bars are kept separate even when their depositional intervals overlap; overlapping bars are stacked within the Source Rock row.

Overburden. Overburden begins after the first specified source rock has finished being deposited. From that point onward, every deposition phase is represented as overburden, regardless of whether the deposited material is source rock, reservoir rock, seal rock or another unit. Hiatuses, uplift and erosion interrupt the overburden record and therefore appear as gaps.

Generation, migration and accumulation. The interval is calculated from the modelled temperature history of each specified source rock. Only periods within generation_window during which the temperature reaches or exceeds all previous temperatures in the history of the source rock are plotted (organic matter that has been heated does not generate again until it gets hotter than before, see plot_basin_model). Generation events are handled per source rock: all events that belong to the same source rock are drawn in the same sub-row, so a source rock with several separate generation events does not create additional rows.

Preservation. Preservation is plotted from the critical moment to the present only when source rock, reservoir rock, seal rock and trap formation have all been specified.

Time axis. Ages are plotted from old to young, with the present at the right-hand side of the chart. Internally the x-axis is transformed so that an increasing plotted position corresponds to a decreasing geological age.

Value

Invisibly returns a named list that contains, for every row of the chart, the list of data frames with the event intervals (start and end in Ma) that were used to draw the chart.

Author(s)

Michiel Arts

References

Magoon, L.B., and Dow, W.G. (1994). The petroleum system. In: Magoon, L.B. and Dow, W.G. (eds), The Petroleum System - From Source to Trap. American Association of Petroleum Geologists Memoir 60, 3–24. doi:10.1306/M60585C1

Allen, P.A., and Allen, J.R. (2013). Basin Analysis: Principles and Application to Petroleum Play Assessment, 3rd ed. Wiley-Blackwell, Chichester, 632 p. ISBN 978-0-470-67377-5.

See Also

run_basin_model, plot_basin_model

Examples

# Event chart of the Huy section (Condroz inlier). The Mousty equivalent
# (unit 2) and the Huy-Vitrival-Bruyere Fms. (unit 4) are the candidate
# source rocks.
res <- run_basin_model(condroz)
plot_event_chart(res,
                 source_units = c(2, 4),
                 reservoir_units = 5,
                 seal_units = 5,
                 trap_formation = data.frame(start = 320, end = 295),
                 critical_moment = 310,
                 main = "Condroz Inlier (Huy)")

Run a one-dimensional burial and thermal history model

Description

Forward models the burial, compaction and temperature history of a one-dimensional sedimentary column (a single "well" or stratigraphic section) through geological time. The basin is described by one row per stratigraphic unit and one row per erosion event. Deposition phases, hiatuses, layer stacks, maximum-burial windows and grain (solid) thicknesses are all derived from that description, so nothing has to be re-specified for each time interval.

Usage

run_basin_model(
  site = NULL,
  units = site$units,
  erosions = site$erosions,
  water_depth = site$water_depth,
  heat_flow = site$heat_flow,
  surface_temp = site$surface_temp,
  dz = 10,
  k_water = 0.6,
  radiogenic = c("cumulative", "legacy"),
  verbose = FALSE,
  genplot = FALSE
)

Arguments

site

A basin description: a list with the elements units, erosions, water_depth, heat_flow and surface_temp, and optionally name (used by the print method and the figures) and source_units (used by plot_basin_model). See brabant, condroz and ardennes for worked examples. Can be left NULL when units, erosions, water_depth, heat_flow and surface_temp are supplied separately.

units

Data frame of stratigraphic units, one row per unit, oldest first. It needs the following columns: name (character label), start and end (deposition interval in Ma, start being the older age), thickness (km, the thickness at maximum burial, that is as deposited and before any later erosion of that unit), phi0 (depositional porosity, between 0 and 1), ck (porosity-depth coefficient in 1/km), K (thermal conductivity of the grains in W/m/K) and A (radiogenic heat production in uW/m3). The optional column rho (grain density in kg/m3) is carried through but is not used in the calculation yet. Typical values of phi0 and ck for common lithologies are given by Sclater and Christie (1980) and Allen and Allen (2013). Defaults to site$units.

erosions

Data frame of erosion events with the columns start and end (Ma, start being the older age) and amount (km of section removed from the top of the column), or NULL when no erosion occurred. Defaults to site$erosions.

water_depth

Water depth in metres. Either a two column matrix or data frame with age (Ma) and water depth (m), the path of a csv file in the same layout, or a single number for a constant water depth. Defaults to site$water_depth.

heat_flow

Heat flow into the base of the column in mW/m2, given in the same formats as water_depth. Defaults to site$heat_flow.

surface_temp

Temperature at the sediment-water interface in degrees Celsius, given in the same formats as water_depth. Defaults to site$surface_temp.

dz

Depth step (m) of the temperature raster that is returned in temp_grid. Default is 10.

k_water

Thermal conductivity of the pore water in W/m/K. Default is 0.6.

radiogenic

How radiogenic heat is carried through the column, either "cumulative" (default) or "legacy". See Details.

verbose

Print a one line summary of the run. Default is FALSE.

genplot

Plot the result with plot_basin_model. Default is FALSE.

Details

Time and depth. Time runs from the oldest age of the oldest unit to the youngest age of the last phase in steps of 1 Myr. Depths are measured from the sediment-water interface, so the water column is not part of the burial depth. A changing water depth therefore does not compact or decompact the column; it only shifts the column relative to sea level.

Phases. The timeline is divided into deposition phases (one for each unit in units), erosion phases (one for each row of erosions) and hiatuses. Every part of the timeline that is covered by neither a unit nor an erosion event is automatically treated as a hiatus (no deposition, no erosion). Phases may not overlap.

Compaction. Porosity decreases exponentially with burial depth (Sclater and Christie, 1980; Allen and Allen, 2013):

\phi(z) = \phi_0 e^{-c z}

where \phi_0 is the depositional porosity (phi0), c is the porosity-depth coefficient (ck) and z the burial depth. The volume of solid grains in a layer does not change during burial, so the present thickness of a layer is found by "decompaction" (the reverse of backstripping, Steckler and Watts, 1978; Allen and Allen, 2013). The decompaction equation is solved by fixed-point iteration for every layer at every time step, rather than on a search grid, so layer boundaries are not quantised.

Irreversible compaction. Compaction is treated as irreversible: the thickness of a layer is never larger than the thickness it had at its deepest burial. An uplifted (exhumed) layer therefore keeps the porosity it acquired at maximum burial instead of re-expanding.

Deposition. During deposition grain material is added at a constant rate. A unit therefore reaches exactly its maximum-burial thickness at the end of its own deposition phase and hands over to the next phase without a step.

Erosion. Erosion removes the stated amount (km) from the top of the column. It may cut through several units, and the layer that is partly eroded keeps the porosity of its remaining part.

Temperature. Temperatures are the steady-state conductive solution for a layered column with radiogenic heat production (Allen and Allen, 2013; Turcotte and Schubert, 2014). The bulk thermal conductivity of a layer is the porosity-weighted arithmetic mean of the conductivity of the grains (K) and of the pore water (k_water). Within a layer of thickness h, the temperature at depth d below the top of the layer is

T(d) = T_{top} + q_{top} d / K - A d^2 / (2 K)

where q_{top} is the heat flow entering the top of the layer and A is its radiogenic heat production. With radiogenic = "cumulative" the surface heat flow is the basal heat flow plus the radiogenic production of the whole column, and the heat flow entering each layer is reduced by the production of the layers above it. With radiogenic = "legacy" every layer instead receives the basal heat flow plus its own production over the full column thickness; use this option only to reproduce results that were calibrated against that older scheme.

Time-dependent boundary conditions. water_depth, heat_flow and surface_temp can be given as a table of (age, value) knots. Values between the knots are linearly interpolated and the curves are held constant beyond the first and last knot. A single number is used as a constant for the whole run. A path to a csv file with the age in the first and the value in the second column is also accepted.

Value

An object of class "basin_model". This is a list with the elements:

units

the unit table, sorted from old to young.

phases

list of all deposition, erosion and hiatus phases. Every phase has a type, a start and end age (Ma), the unit (row of units, NA when not a deposition phase) and the eroded amount (km).

ages

the ages (Ma) of the time steps, from old to young.

water_depth

water depth (m) at each age.

heat_flow

basal heat flow (mW/m2) at each age.

surface_temp

surface temperature (degrees Celsius) at each age.

total_thickness

thickness (m) of the preserved column at each age.

profile

a data frame with one row for every layer at every age: age, unit, name, z_top, z_base and z_mid (burial depth in m below the sediment-water interface), h (thickness in m), porosity, K (bulk thermal conductivity in W/m/K), A (heat production in W/m3), q_top (heat flow into the top of the layer in W/m2), T_top, T_base and T_mid (degrees Celsius) and water_depth (m).

depth_grid

depths (m below sea level) of the temperature raster.

temp_grid

temperature raster (degrees Celsius) with one row per age and one column per element of depth_grid.

radiogenic

the radiogenic heat scheme that was used.

site_name

the name of the basin.

source_units

the source_units of site, if any.

Author(s)

The porosity-depth law, the decompaction approach and the heat conduction equation follow the textbook treatment of Allen and Allen (2013).

References

Allen, P.A., and Allen, J.R. (2013). Basin Analysis: Principles and Application to Petroleum Play Assessment, 3rd ed. Wiley-Blackwell, Chichester, 632 p. ISBN 978-0-470-67377-5.

Sclater, J.G., and Christie, P.A.F. (1980). Continental stretching: an explanation of the post-mid-Cretaceous subsidence of the central North Sea Basin. Journal of Geophysical Research, 85, 3711–3739. doi:10.1029/JB085iB07p03711

Steckler, M.S., and Watts, A.B. (1978). Subsidence of the Atlantic-type continental margin off New York. Earth and Planetary Science Letters, 41, 1–13. doi:10.1016/0012-821X(78)90036-5

Turcotte, D.L., and Schubert, G. (2014). Geodynamics, 3rd ed. Cambridge University Press, Cambridge.

See Also

plot_basin_model, plot_event_chart, BasinmodlR_datasets

Examples

# Burial and thermal history of the Huy section (Condroz inlier, Belgium)
# using the example basin that is shipped with the package.
res <- run_basin_model(condroz)
res

# The layer table has one row for every layer at every age.
# Show the burial depth and temperature of the Mousty equivalent at 300 Ma.
prof <- res$profile
prof[prof$age == 300 & prof$unit == 2, c("age", "name", "z_top", "z_base", "T_mid")]

# Build your own basin: three units, one erosion event and constant boundary
# conditions. Thicknesses and erosion amounts are in km, ages in Ma.
my_basin <- list(
  name = "Synthetic foreland basin",
  units = data.frame(
    name = c("Sandstone", "Shale", "Shaly sand"),
    start = c(300, 250, 200),   # older age of the deposition interval (Ma)
    end = c(250, 200, 150),     # younger age of the deposition interval (Ma)
    thickness = c(1.2, 0.8, 0.5),
    phi0 = c(0.49, 0.63, 0.56),  # depositional porosity
    ck = c(0.27, 0.51, 0.39),    # porosity-depth coefficient (1/km)
    K = c(3.0, 1.7, 2.4),        # grain thermal conductivity (W/m/K)
    A = c(1.0, 1.8, 1.4)         # radiogenic heat production (uW/m3)
  ),
  # 0.4 km of section is removed between 100 and 0 Ma. The interval between
  # 150 and 100 Ma is not covered by a unit or an erosion event and
  # becomes a hiatus automatically.
  erosions = data.frame(start = 100, end = 0, amount = 0.4),
  water_depth = 0,     # constant water depth (m)
  heat_flow = 60,      # constant basal heat flow (mW/m2)
  surface_temp = 15    # constant surface temperature (degrees Celsius)
)
my_res <- run_basin_model(my_basin)
max(my_res$profile$T_base)   # maximum temperature reached by any layer