CNPDIAperturb {CNPDIA} | R Documentation |
CNPDIAperturb
dynamically runs the CNPDIA model with event-like perturbations.
CNPDIAperturbSettings
and CNPDIAperturbFluxes
retrieve the settings and perturbation fluxes.
CNPDIAperturb (parms = list(), times = 0:365, spinup = NULL, yini = NULL, gridtype = 1, Grid = NULL, porosity = NULL, bioturbation = NULL, irrigation = NULL, surface = NULL, diffusionfactor = NULL, dynamicbottomwater = FALSE, perturbType = "mix", perturbTimes = seq(from = 0, to = max(times), by = 365), perturbDepth = 5, concfac = 1, CfluxForc = NULL, FePfluxForc = NULL, CaPfluxForc = NULL, O2bwForc = NULL, NO3bwForc = NULL, NO2bwForc = NULL, NH3bwForc = NULL, ODUbwForc = NULL, PO4bwForc = NULL, DICbwForc = NULL, wForc = NULL, biotForc = NULL, irrForc = NULL, rFastForc = NULL, rSlowForc = NULL, pFastForc = NULL, MPBprodForc = NULL, gasfluxForc = NULL, HwaterForc = NULL, ratefactor = NULL, verbose = FALSE, ...) mixSolid (C, depth = 5, grid = setup.grid.1D(N = 100, L = 30), por = 0.5) mixLiquid (C, depth = 5, grid = setup.grid.1D(N = 100, L = 30), por = 0.5, BW) erodeSolid (C, depth = 5, grid = setup.grid.1D(N = 100, L = 30), por = 0.5) erodeLiquid (C, depth = 5, grid = setup.grid.1D(N = 100, L = 30), por = 0.5) depositSolid (C, depth = 5, grid = setup.grid.1D(N = 100, L = 30), por = 0.5, concfac = 1) depositLiquid (C, depth = 5, grid = setup.grid.1D(N = 100, L = 30), por = 0.5, BW) CNPDIAperturbFluxes(out, which = NULL) CNPDIAperturbSettings(out)
parms |
A list with parameter names and values. Possible parameters are:
See details of CNPDIAsolve. |
gridtype |
Type of grid: |
times |
Output times for the dynamic run |
spinup |
Spinput times for the dynamic run; not used for output; the outputted simulation starts from the final values of the spinup run. |
CfluxForc, FePfluxForc, CaPfluxForc |
|
O2bwForc, NO3bwForc, NO2bwForc, NH3bwForc, ODUbwForc, PO4bwForc, DICbwForc |
|
wForc, biotForc, irrForc |
|
rFastForc, rSlowForc, pFastForc, MPBprodForc |
|
gasfluxForc |
|
HwaterForc |
|
dynamicbottomwater |
If TRUE, then the concentrations in the water overlying the sediment will also be dynamically described,
and with water height equal to |
ratefactor |
|
Grid |
If specified: either an object, as returned by |
porosity |
If specified, either an object with porosities ([-]) as returned by |
bioturbation |
If specified, either an object with bioturbation rates (units [cm2/d]) as returned by |
irrigation |
If specified, either an object with irrigation rates (units [/d]) as returned by |
surface |
If specified, either an object with surface areas (units [cm2]) as returned by Note that the surface values should be consistent withe the |
diffusionfactor |
The multiplication factor necessary to go from molecular diffusion to effective sediment diffusion, i.e. that takes into account tortuosity.
If specified, either an object with these factors ([-]) as returned by |
yini |
The condition at which to inialise the dynamic simulation, i.e. a vector or matrix, with the values of
|
perturbType |
how to perturb, one of "mix", "deposit", "erode". Several options can be chosen, in which case the perturbations will be done in the order defined. |
perturbTimes |
times at which the perturbations should take place. Either one value, or a vector of the same length as |
perturbDepth, depth |
the depth of the perturbation, in cm. |
concfac |
only when perturbType = "deposit": the factor at which the available concentration should be increased or decreased. |
C |
The concentration vector to be perturbed; a vector that is compatible with |
BW |
The concentration in the overlying water, one value. |
grid |
a list as generated with setup.grid.1D[ReacTran] that contains the grid definition |
por |
the sediment porosity; either as generated with setup.prop.1D[ReacTran],
or a vector of length = |
verbose |
If TRUE, will write progession to the screen . |
out |
an output object returned by CNPDIAperturb or CNPDIAdyna. |
which |
if not |
... |
Any argument passed to the dynamic solver (ode.1D[deSolve]) |
Several parameters can also be described as forcing functions. They are:
Cflux, FePflux, CaPflux, O2bw, NO3bw, NO2bw, NH3bw, ODUbw, PO4bw, DICbw, w,
biot, irr, rFast, rSlow, pFast, MPBprod, gasflux, Hwater
.
The forcing functions are prescribed as a list
that either contains a data series or specifies a periodic signal.
The list is defined as: list(data = NULL, amp = 0, period = 365, phase = 0, pow = 1, min = 0)
Forcing functions a data series are set with item data
contains time series for the parameter - a matrix with times (first column) and values, second column. The values should be in the same units as the parameters.
The time series should embrace both arguments times
and spinup
.
if a periodic signal, the list should contain amp, period, phase, pow
and min
: parameters determining the periodicity of the seasonal signal in the same units as the parameters.
From this the forcing function time series is estimated, e.g. for CfluxForc as: max(min, Cflux*(1 + (amp*sin((times-phase)/period*2*pi))^pow)
, where Cflux
is the parameter value. Used only if data
is NULL
. If amp
= 0, or pow
= 0, then the forcing will be kept constant and equal to the parameter value.
CNPDIAperturb
returns a matrix of type deSolve
, as generated by the
solver ode.1D from R-package deSolve
.
It contains several output columns, the first is time
. The meaning and units of these columns can be assessed via the R-functions:
CNPDIAsvar()
, CNPDIA1D()
, CNPDIA0D()
. See CNPDIA0D.
The instantaneous release/gain is saved in the attribute perturbFluxes
and the settings in attribute perturbSettings
. They can be retrieved with functions
CNPDIAperturbFluxes
, and CNPDIAperturbSettings
Functions depositLiquid, depositSolid, mixLiquid, mixSolid, erodeLiquid, erodeSolid
return
a list, with the values:
C
, the altered concentration profile
Flux
, the instantaneous gain in mass, estimated as (integral(Cnew) - integral(C)).
The model application starts by estimating the steady-state condition
of the model. This steady-state condition is then used as a starting
condition for a dynamic simulation, with perturbations as provided in perturbTimes
.
Mixing will homogenise the perturbed depth of the sediment (perturbType = "mix"
).
Erosion will remove the perturbed depth of the sediment (perturbType = "erode"
.
Deposition will add a layer of sediment (perturbType = "deposit"
.
All these events can be combined; the sequence of events is as provided, i.e.
perturbType = c("mix", "erode")
will not give the same results as perturbType = c("erode", "mix")
.
Karline Soetaert
Soetaert K, PMJ Herman and JJ Middelburg, 1996a. A model of early diagenetic processes from the shelf to abyssal depths. Geochimica Cosmochimica Acta, 60(6):1019-1040.
Soetaert K, PMJ Herman and JJ Middelburg, 1996b. Dynamic response of deep-sea sediments to seasonal variation: a model. Limnol. Oceanogr. 41(8): 1651-1668.
# ======================================== # One perturbation at the start # ======================================== out <- CNPDIAperturb() par(mar = c(3,3,3,3)) image2D(out, ylim = c(20, 0), mfrow = c(4, 3)) # ======================================== # Mixing at selected times # ======================================== out2 <- CNPDIAperturb(perturbTime = c(0, 100, 200, 300), perturbType = "mix", perturbDepth = 10) image2D(out2, ylim = c(20, 0), mfrow = c(4, 3)) CNPDIAbudgetO2(out) # ======================================== # Erosion at selected times # ======================================== out3 <- CNPDIAperturb(perturbTime = c(0, 100, 200, 300), perturbType = "erode", perturbDepth = 5) image2D(out3, ylim = c(20, 0), mfrow = c(4, 3)) (PertFluxes <- CNPDIAperturbFluxes(out3)) CNPDIAperturbSettings(out3) # ======================================== # Several subsequent perturbations # ======================================== out4 <- CNPDIAperturb(perturbTime = c(0, 100, 200, 300), perturbType = c("mix", "erode"), perturbDepth = c(10, 5)) image2D(out4, ylim = c(20, 0), mfrow = c(4, 3)) # ======================================== # Create perturbation profiles # ======================================== # create reference solution pm <- par(mfrow = c(1, 2)) DIA <- CNPDIAsolve() # ------------------------ # perturb solid fraction # ------------------------ SDET <- DIA$y[,"SDET"] Grid <- DIA$Grid por <- DIA$porosity # perturb it mix <- mixSolid (SDET, grid = Grid, por = por) deposit <- depositSolid(SDET, grid = Grid, por = por) erode <- erodeSolid (SDET, grid = Grid, por = por) matplot(x = cbind(SDET, mix$C, erode$C, deposit$C), main = "Solid", y = Grid$x.mid, type = "l", lwd = 3, ylim = c(20, 0)) legend("bottomright", col = 1:4, lty = 1:4, legend = c("default", "mix", "erode", "deposit")) cat("instantaneous fluxes, nmol/cm2\n", "mixing: ", mix$Flux, "erosion: ", erode$Flux, "deposition: ", deposit$Flux, "\n") # ------------------------ # perturb liquid fraction # ------------------------ NO3 <- DIA$y[,"NO3"] Grid <- DIA$Grid por <- DIA$porosity # perturb it mix <- mixLiquid (NO3, grid = Grid, BW = 10, por = por) deposit <- depositLiquid(NO3, grid = Grid, BW = 10, por = por) erode <- erodeLiquid (NO3, grid = Grid, por = por) matplot(x = cbind(NO3, mix$C, erode$C, deposit$C), main = "Liquid", y = Grid$x.mid, type = "l", lwd = 3, ylim = c(20, 0)) legend("bottomright", col = 1:4, lty = 1:4, legend = c("default", "mix", "erode", "deposit")) cat("instantaneous fluxes, nmol/cm2\n", "mixing: ", mix$Flux, "erosion: ", erode$Flux, "deposition: ", deposit$Flux) par(mfrow = pm)