CNPDIAsolve {CNPDIA}R Documentation

Steady-state solution of C, N, P and O2 dynamics in the sediment.

Description

CNPDIAsolve finds the steady-state solution of the CNPDIA model.

Usage

  CNPDIAsolve (parms = list(), yini = NULL, gridtype = 1, Grid = NULL, 
     porosity = NULL, bioturbation = NULL, irrigation = NULL, surface = NULL, 
     diffusionfactor = NULL, dynamicbottomwater = FALSE, 
     ratefactor = NULL, verbose = FALSE, ...)
  

Arguments

parms

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

See details.

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, mixdepth and mixatt as: biot * exp(-pmax(0, (x.int-mixL))/mixatt), where x.int is the distance, from the surface of the layer interface. Note that the bioturbation values should be consistent withe 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, mixdepth and mixatt as: irr * exp(-pmax(0, (x-mixdepth)/mixxatt)), where x is the distance, from the surface, of the layer centres. Note that the irrigation values should be consistent withe 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 = rep(1, 101) 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 withe 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 withe 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 and constant. 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.

verbose

If TRUE, will write progession to the screen .

...

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

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

It contains, a.o. the elements:

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
#===========================================

  CNPDIAparms()
  par(mar = c(3,3,3,3))

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

  out  <- CNPDIAsolve()
  out2 <- CNPDIAsolve(parms = list(Cflux = 200*1e5/12/365))
  out3 <- CNPDIAsolve(parms = list(Cflux = 2*1e5/12/365))
  plot(out, out2, out3, xyswap = TRUE, grid = CNPDIAdepth(out), 
    ylim = c(20,0), mfrow = c(3, 4))

  CNPDIAbudgetO2(out)  
  
#===========================================
# long-distance reactions
#===========================================

  out2 <- CNPDIAsolve(parms = list(Cflux = 100*1e5/12/365))
  out2b <- CNPDIAsolve(parms = list(Cflux = 100*1e5/12/365, 
               rSurfODUox = 10, ODUoxdepth = 1), yini = out2$y)
  out2c <- CNPDIAsolve(parms = list(Cflux = 100*1e5/12/365,  
               rSurfODUox = 10, ODUoxdepth = 5), yini = out2b$y)
  plot(out2, out2b, out2c, xyswap = TRUE, grid = CNPDIAdepth(out2), 
     which = c("O2","ODU"), ylim = list(c(1,0), c(10,0)), 
     log = list("","x"))

  plot(out2b, out2c, xyswap = TRUE, grid = CNPDIAdepth(out2), 
     which = c("O2","ODU"), ylim = list(c(1,0), c(10,0)))

#===========================================
# Dynamic bottom water concentrations
#===========================================

  std <- CNPDIAsolve(parms = list(Cflux = 20*1e5/12/365))
  BWs <- CNPDIAparms(std, as.vector = TRUE)[c("O2bw", "NO3bw", 
       "NO2bw", "NH3bw", "ODUbw", "DICbw", "PO4bw")]
  yini <- rbind(c(FDETbw = 0, SDETbw = 0, BWs, FePbw = 0, 
        CaPbw = 0, Padsbw=0), std$y)
  out2 <- CNPDIAdyna(parms = list(Cflux = 20*1e5/12/365, Hwater = 50), 
      dynamicbottomwater = TRUE, times = seq(0, 1, by = 0.01), 
      yini = yini)
  
  plot(out2, which = c("O2bw", "NO3bw", "O2flux", "NO3flux"))
  image(out2, which = c("O2","NO3","NH3","DIC"), legend = TRUE)
## Not run: 
# 
  
## End(Not run)

[Package CNPDIA version 1.0 Index]