| 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 |
| 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):
unitsone row per stratigraphic unit, oldest first, with the deposition interval (
startandendin Ma,startbeing 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).erosionsone 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_depthwater depth (m) as (age, value) knots.
heat_flowheat flow into the base of the column (mW/m2) as (age, value) knots.
surface_temptemperature 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_modelruns 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_modelcomposite 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_chartpetroleum 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
namecharacter, the label that is used by
printand by the figures.unitsdata frame, one row per stratigraphic unit, oldest first, with the columns
name(character),startandend(deposition interval in Ma,startbeing 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) andrho(grain density, kg/m3, carried through but not used in the calculation yet).erosionsdata frame with
start,end(Ma) andamount(km of section that is removed from the top of the column).water_depthtwo column matrix with the age (Ma) and the water depth (m).
heat_flowtwo column matrix with the age (Ma) and the heat flow into the base of the column (mW/m2).
surface_tempdata frame with the age (Ma) and the temperature at the sediment-water interface (degrees Celsius), at 10 Myr resolution.
source_unitsinteger, the rows of
unitsthat get a generation panel inplot_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 |
... |
For |
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 |
window |
Numeric vector of length two, the temperature window in degrees
Celsius that is used for the generation panels. Default is |
window_units |
Which units get a generation panel, given as row numbers
of |
n_levels |
Number of colour levels in the temperature raster. Default is 250. |
scale |
How the colours are spread over the temperature range, either
|
contour_n |
Approximate number of isotherms. Default is 20. |
xmax |
Oldest age shown in Ma. |
tick |
Spacing of the age axis labels in Myr. Default is 25. |
new_device |
Open a new graphics window before plotting. Leave this
|
use |
The temperature that is used for the generation window test, one of
|
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 |
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 |
source_units |
Source-rock units, given as row numbers of
|
reservoir_units |
Reservoir-rock units, given as row numbers or unit
names. Default is |
seal_units |
Seal-rock units, given as row numbers or unit names.
Default is |
trap_formation |
Optional data frame with the columns |
generation_window |
Numeric vector of length two with the temperature
window (degrees Celsius) that is used to identify generation, migration and
accumulation. Default is |
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 |
temperature_column |
Name of the temperature column of
|
xmax |
Oldest age shown on the x-axis in Ma. |
tick |
Spacing between the age labels of the x-axis in Myr. Default is 25. |
col |
Default colour of the event bars. Default is |
critical_col |
Colour of the critical moment marker. Default is
|
main |
Plot title. Default is an empty string. |
xlab |
Label of the x-axis. Default is |
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 |
Data frame of stratigraphic units, one row per unit, oldest
first. It needs the following columns: |
erosions |
Data frame of erosion events with the columns |
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
|
heat_flow |
Heat flow into the base of the column in mW/m2, given in the
same formats as |
surface_temp |
Temperature at the sediment-water interface in degrees
Celsius, given in the same formats as |
dz |
Depth step (m) of the temperature raster that is returned in
|
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
|
verbose |
Print a one line summary of the run. Default is
|
genplot |
Plot the result with |
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:
unitsthe unit table, sorted from old to young.
phaseslist of all deposition, erosion and hiatus phases. Every phase has a
type, astartandendage (Ma), theunit(row ofunits,NAwhen not a deposition phase) and the erodedamount(km).agesthe ages (Ma) of the time steps, from old to young.
water_depthwater depth (m) at each age.
heat_flowbasal heat flow (mW/m2) at each age.
surface_tempsurface temperature (degrees Celsius) at each age.
total_thicknessthickness (m) of the preserved column at each age.
profilea data frame with one row for every layer at every age:
age,unit,name,z_top,z_baseandz_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_baseandT_mid(degrees Celsius) andwater_depth(m).depth_griddepths (m below sea level) of the temperature raster.
temp_gridtemperature raster (degrees Celsius) with one row per age and one column per element of
depth_grid.radiogenicthe radiogenic heat scheme that was used.
site_namethe name of the basin.
source_unitsthe
source_unitsofsite, 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