FESDIAsolve {FESDIA}R Documentation

Steady-state solution of the FESDIA model, calculating C, N, P, O2, Fe and S dynamics in the sediment.

Description

FESDIAsolve finds the steady-state solution of the FESDIA model; pH can be calculated afterwards, ie ignoring carbonate dynamics.

Usage

  FESDIAsolve (parms = list(), yini = NULL, gridtype = 1, Grid = NULL, 
      porosity = NULL, bioturbation = NULL, irrigation = NULL,   
      surface = NULL, diffusionfactor = NULL, 
      dynamicbottomwater = FALSE, ratefactor = NULL, 
      calcpH = FALSE, verbose = FALSE, 
      method = NULL, times = c(0, 1e+06), ...)

Arguments

parms

A list with parameter values.Available parameters can be listed using function FESDIAparms.

See details of FESDIAparms.

gridtype

Type of grid: 1 for cartesian, 2 for cylindrical, 3 for spherical.

Grid

If specified: either an object, as returned by setup.grid.1D from the package ReacTran, a vector of length 101 with the transport distances (from mid to mid of layers, upper interface = diffusive boundary layer), or one number with the constant layer thickness. If NULL, it is defined as setup.grid.1D(x.up = 0, dx.1 = 0.01, N = 100, L = 100), i.e. the total length is 100 cm, the first layer is 0.01 cm thick and layers are increasing with depth for 100 layers.

porosity

If specified, either an object with porosities ([-]) as returned by setup.prop.1D from the package ReacTran, a vector of length 101 with the porosities defined at the layer interfaces, or one number with the constant porosity. If NULL, it is defined by the parameters por0, pordeep and porcoeff as: (pordeep + (por0 - pordeep) * exp(-pmax(0, x.int)/porcoeff)), where x.int is the distance, from the surface of the layer interface. Note that the porosity values should be consistent withe the Grid - and should be inbetween 0 and 1.

bioturbation

If specified, either an object with bioturbation rates (units [cm2/d]) as returned by setup.prop.1D from the package ReacTran, a vector of length 101 with the bioturbation defined at the layer interfaces, or one number with the constant bioturbation. If NULL, it is defined by the parameters biot, biotdepth and biotatt as: biot * exp(-pmax(0, (x.int-biotdpeth))/biotatt), where x.int is the distance, from the surface of the layer interface. Note that the bioturbation values should be consistent with the Grid.

irrigation

If specified, either an object with irrigation rates (units [/d]) as returned by setup.prop.1D from the package ReacTran, a vector of length 100 with the irrigation rates defined at the layer centres, or one number with the constant rates. If NULL, it is defined by the parameters irr, irrdepth and irratt as: irr * exp(-pmax(0, (x-irrdepth)/irratt)), where x is the distance, from the surface, of the layer centres. Note that the irrigation values should be consistent with the Grid.

surface

If specified, either an object with surface areas (units [cm2]) as returned by setup.prop.1D from the package ReacTran, a vector of length 101 with the surface areas defined at the layer interfaces, or one number with the constant surface area. If NULL, it is defined by the parameter gridtype, and the Grid as: surface = 1 for gridtype == 1, surface = rev(2*pi*Grid$x.int) for gridtype == 2 and surface = rev(pi*(Grid$x.int)^2) for gridtype == 3.

Note that the surface values should be consistent with the Grid.

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 setup.prop.1D from the package ReacTran, a vector of length 101 with these factors defined at the layer interfaces, or one number with the constant factor. If NULL, it is set equal to the porosity. Note that the factors should be consistent with the Grid.

yini

Initial guess of the steady-state solution.

dynamicbottomwater

If TRUE, then the concentrations in the water overlying the sediment will also be dynamically described, and with water height equal to Hwater. Note that this will slow down the simulation.

ratefactor

NULL or a list, detailing the forcing function for the biogeochemical rate multiplication factor. If not specified (or NULL), then it is assumed to be 1, if a data series, then the mean of the data will be used. If a list, it should contain either a data time series (list(data = )) or parameters determining the periodicity of the seasonal signal (defined as list(data = NULL, amp = 0, period = 365, phase = 0, pow = 1, min = 0). see details.

method

The method to be used for estimating steady-state, will be passed as argument to the solver steady.1D. The default is to use the "stode" method. Other options are to use "stodes", "runsteady" or "mix". The last option combines runsteady with stode, i.e. the model is first dynamically run for the times specified in argument times, and then the final value from this run used as initial guess for the steady-state estimated using stode

times

start and end time of the dynamic run - only if method is one of "runsteady" or "mix".

verbose

If TRUE, will write progession to the screen .

calcpH

if TRUE, the pH will be calculated - note that this will not include the effects of calcium carbonate precipitation or dissolution on pH. Note also that the pH can be estimated afterwards by running FESDIApH. See examples.

...

Any argument passed to the steady-state solver.

Details

To solve the model, a steady-state solver from package rootSolve (here we used steady.1D) is used.

Value

FESDIAsolve returns an object of class FESDIAstd or PHDIAstd, and of class steady1D, as generated by the solvers from R-package rootSolve (steady.1D[rootSolve]).

It contains a.o. the elements:

The meaning and units of these columns can be assessed via the R-functions:

FESDIAsvar(), FESDIA1D(), FESDIA0D(). See FESDIA0D.

Author(s)

Karline Soetaert

References

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.

Examples


#===========================================
# Show parameter values
#===========================================

  FESDIAparms()

#===========================================
# Runs with different carbon fluxes
#===========================================

  out  <- FESDIAsolve(parms = list(Cflux = 2*1e5/12/365), calcpH = TRUE)
  out2 <- FESDIAsolve(parms = list(Cflux = 20*1e5/12/365), calcpH = TRUE)
  out3 <- FESDIAsolve(parms = list(Cflux = 200*1e5/12/365), calcpH = TRUE)

  par(mar = c(3,3,3,2))
  plot(out, out2, out3, ylim = c(20, 0), mfrow = c(4, 4), 
    which = c(1:14,16), lwd = 2, lty = 1)

  plot(out, out2, out3, which = "pH", ylim = c(20,0))

#===========================================
# Runs that are difficult to solve
#===========================================

# simple run
  out1 <- FESDIAsolve(parms = c(por0 = 0.5, MPBprod = 0))

# simple run used as initial condition for second run
  out2 <- FESDIAsolve(parms = c(por0 = 0.5, MPBprod = 1e3), yini = out1$y)

# second run used as initial condition for third run
# use mixed method: first dynamic run, then steady-state solver
  out3 <- FESDIAsolve(parms = c(por0 = 0.5, MPBprod = 1e4), yini = out2$y, method = "mixed")

  out4 <- FESDIAsolve(parms = c(por0 = 0.5, MPBprod = 5e4), yini = out3$y, 
     method = "runsteady", times = c(0, 1e6))
  
  FESDIAbudgetO2(out1, out2, out3, out4, which = "Rates") 


[Package FESDIA version 1.0 Index]