CNPDIAdyna {CNPDIA} | R Documentation |
CNPDIAdyna
dynamically runs the CNPDIA model.
CNPDIAdyna (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, 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, ...)
parms |
A list with parameter values. Available parameters can be listed using function CNPDIAparms. See details of CNPDIAparms. |
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 |
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 |
diffusionfactor |
The multiplication factor necessary to go from molecular diffusion to effective sediment diffusion, i.e. that takes into account tortuosity. Effective diffusion is then calculated as molecular diffusion*diffusionfactor.
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
|
verbose |
If TRUE, will write progession to the screen . |
... |
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 imposed by a data series are set with item data
, and contains a time series to replace the parameter. data
should be a two-columned 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, and in the same units as the parameters.
From this the forcing function time series is estimated, e.g. for CfluxForc
as follows: CfluxForc = 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.
A matrix of class CNPDIAdyn
and deSolve
, as generated by the solver from R-package deSolve
(ode.1D).
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 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 sinusoidal forcing of organic carbon deposition flux (truncated to be >= 0).
If only one year is run, the dynamic run is likely not at equilibrium with the transient forcing.
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.
#=========================================== # Show parameter values #=========================================== CNPDIAparms() #=========================================== # Default run #=========================================== out <- CNPDIAdyna() par(mar = c(3,3,3,3)) image2D(out, ylim = c(20, 0), mfrow = c(3, 4), which = 1:12) CNPDIAbudgetO2(out) #=========================================== # Altered parameter values #=========================================== out2 <- CNPDIAdyna(parms = list(pFast = 0.95, Cflux = 50*1e5/12/365), spinup = 0:365) par(mar = c(3,3,3,3)) image2D(out2, ylim = c(20, 0), mfrow = c(2, 2), which = c("TOC", "O2", "NO3", "PO4")) #=========================================== # A forcing function flux timeseries #=========================================== Cflux <- cbind (time = c(0, 100, 150, 175, 200, 250, 365), flux = c(1, 1, 1000, 800, 1200, 10, 1)) out3 <- CNPDIAdyna(parms = list(pFast = 0.9), CfluxForc = list(data = Cflux), times = 0:365) image2D(out3, mfrow = c(3, 3), ylim = c(20, 0), which = 1:9) plot(out3, which = c("OrgCflux", "O2flux"), mfrow = NULL) head(CNPDIA0D(out3)) CNPDIAbudgetC(out3) #=========================================== # a simple forcing function # (varying the part of deposition in fast decay C) #=========================================== out4 <- CNPDIAdyna(parms = list(pFast = 0.9), pFastForc = list(amp = 0.8), times = 0:365) image2D(out4, mfrow = c(3, 3), which = 1:9, ylim = c(20, 0)) plot(out4, which = c("OrgCflux", "O2flux"), mfrow = NULL) head(CNPDIA0D(out4)) CNPDIAbudgetC(out4) #=========================================== # long-distance reactions #=========================================== out <- CNPDIAsolve(parms = list(Cflux = 100*1e5/12/365)) dyn <- CNPDIAdyna(parms = list(Cflux = 100*1e5/12/365), yini = out$y, times = seq(0,5, length.out = 30)) image2D(dyn, which = c("O2","ODU"), ylim = list(c(1,0), c(10,0))) out2 <- CNPDIAsolve(parms = list(Cflux = 100*1e5/12/365, rSurfODUox = 5, ODUoxdepth = 5), yini = out$y) plot(out, out2, which = c("O2","ODU"), ylim = list(c(1,0), c(10,0))) dyn <- CNPDIAdyna(parms = list(Cflux = 100*1e5/12/365, rSurfODUox = 5, ODUoxdepth = 5), yini = out$y, times = seq(0,5, length.out = 30)) image2D(dyn, which = c("O2","ODU"), ylim = list(c(1,0), c(10,0)))