tran.1D: General One-Dimensional Advective-Diffusive Transport
In ReacTran: Reactive Transport Modelling in 1d, 2d and 3d
Description
Usage
Arguments
Details
Value
Note
Author(s)
References
See Also
Examples
Description
Estimates the transport term (i.e. the rate of change of a concentration
due to diffusion and advection) in a one-dimensional model of a
liquid (volume fraction constant and equal to one) or in a porous medium
(volume fraction variable and lower than one).
The interfaces between grid cells can have a variable cross-sectional area,
e.g. when modelling spherical or cylindrical geometries (see example).
Usage
1
2
3
4
5
Arguments
C
concentration, expressed per unit of phase volume, defined at the
centre of each grid cell. A vector of length N [M/L3]
C.up
concentration at upstream boundary. One value [M/L3]
C.down
concentration at downstream boundary. One value [M/L3]
flux.up
flux across the upstream boundary, positive = INTO model
domain. One value, expressed per unit of total surface [M/L2/T].
If NULL, the boundary is prescribed as
a concentration or a convective transfer boundary.
flux.down
flux across the downstream boundary, positive = OUT
of model domain. One value, expressed per unit of total surface [M/L2/T].
If NULL, the boundary is prescribed as
a concentration or a convective transfer boundary.
a.bl.up
convective transfer coefficient across the upstream
boundary layer. Flux = a.bl.up*(C.up-C0). One value [L/T]
a.bl.down
convective transfer coefficient across the downstream
boundary layer (L). Flux = a.bl.down*(CL-C.down).
One value [L/T]
D
diffusion coefficient, defined on grid cell interfaces.
One value, a vector of length N+1 [L2/T], or a 1D property list; the list
contains at least the element int (see setup.prop.1D)
[L2/T]
v
advective velocity, defined on the grid cell
interfaces. Can be positive (downstream flow) or negative (upstream flow).
One value, a vector of length N+1 [L/T], or a 1D property list; the list
contains at least the element int (see setup.prop.1D)
[L/T]
AFDW
weight used in the finite difference scheme for advection,
defined on grid cell interfaces; backward = 1, centred = 0.5, forward = 0;
default is backward. One value, a vector of length N+1, or a
1D property list; the list contains at least the element int
(see setup.prop.1D) [-]
VF
Volume fraction defined at the grid cell interfaces. One value,
a vector of length N+1, or a 1D property list; the list
contains at least the elements int and mid
(see setup.prop.1D) [-]
A
Interface area defined at the grid cell interfaces. One value,
a vector of length N+1, or a 1D grid property list; the list
contains at least the elements int and mid
(see setup.prop.1D) [L2]
dx
distance between adjacent cell interfaces (thickness of grid
cells). One value, a vector of length N, or a 1D grid list containing
at least the elements
dx and dx.aux (see setup.grid.1D) [L]
full.check
logical flag enabling a full check of the consistency
of the arguments (default = FALSE; TRUE slows down execution
by 50 percent)
full.output
logical flag enabling a full return of the output
(default = FALSE; TRUE slows down execution by 20 percent)
Details
The boundary conditions are either
(1) zero-gradient.
(2) fixed concentration.
(3) convective boundary layer.
(4) fixed flux.
The above order also shows the priority. The default condition is the
zero gradient. The fixed concentration condition overrules the zero gradient.
The convective boundary layer condition overrules the fixed concentration
and zero gradient. The fixed flux overrules all other specifications.
Ensure that the boundary conditions are well defined: for instance, it
does not make sense to specify an influx in a boundary cell with the advection
velocity pointing outward.
Transport properties:
The diffusion coefficient (D),
the advective velocity (v),
the volume fraction (VF), the interface surface (A),
and the advective finite difference weight (AFDW)
can either be specified as one value, a vector or a 1D property list
as generated by setup.prop.1D.
When a vector, this vector must be of length N+1, defined at all grid
cell interfaces, including the upper and lower boundary.
The finite difference grid (dx) is specified either as
one value, a vector or a 1D grid list, as generated by setup.grid.1D.
Value
dC
the rate of change of the concentration C due to transport,
defined in the centre of each grid cell. The rate of change is expressed
per unit of phase volume [M/L3/T]
C.up
concentration at the upstream interface. One value [M/L3]
only when (full.output = TRUE)
C.down
concentration at the downstream interface. One value [M/L3]
only when (full.output = TRUE)
dif.flux
diffusive flux across at the interface of each grid cell.
A vector of length N+1 [M/L2/T]
only when (full.output = TRUE)
adv.flux
advective flux across at the interface of each grid cell.
A vector of length N+1 [M/L2/T]
only when (full.output = TRUE)
flux
total flux across at the interface of each grid cell. A vector
of length N+1 [M/L2/T].
only when (full.output = TRUE)
flux.up
flux across the upstream boundary, positive = INTO model
domain. One value [M/L2/T]
flux.down
flux across the downstream boundary, positive = OUT of
model domain. One value [M/L2/T]
Note
The advective equation is not checked for mass conservation. Sometimes, this is
not an issue, for instance when v represents a sinking velocity of
particles or a swimming velocity of organisms.
In others cases however, mass conservation needs to be accounted for.
To ensure mass conservation, the advective velocity must obey certain
continuity constraints: in essence the product of the volume fraction (VF),
interface surface area (A) and advective velocity (v) should be constant.
In sediments, one can use setup.compaction.1D to ensure that
the advective velocities for the pore water and solid phase meet these
constraints.
In terms of the units of concentrations and fluxes we follow the convention
in the geosciences.
The concentration C, C.up, C.down as well at the rate of
change of the concentration dC are always expressed per unit of
phase volume (i.e. per unit volume of solid or liquid).
Total concentrations (e.g. per unit volume of bulk sediment) can be obtained by
multiplication with the appropriate volume fraction. In contrast, fluxes are
always expressed per unit of total interface area (so here the volume fraction
is accounted for).
Author(s)
Filip Meysman <filip.meysman@nioz.nl>,
Karline Soetaert <karline.soetaert@nioz.nl>
References
Soetaert and Herman (2009). A practical guide to ecological modelling -
using R as a simulation platform. Springer
See Also
tran.volume.1D for a discretisation the transport equation using finite volumes.
advection.1D, for more sophisticated advection schemes
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
## =============================================================================
## EXAMPLE 1: O2 and OC consumption in sediments
## =============================================================================
# this example uses only the volume fractions
# in the reactive transport term
#====================#
# Model formulation #
#====================#
# Monod consumption of oxygen (O2)
O2.model <- function (t = 0, O2, pars = NULL) {
tran <- tran.1D(C = O2, C.up = C.ow.O2, D = D.grid,
v = v.grid, VF = por.grid, dx = grid)$dC
reac <- - R.O2*(O2/(Ks+O2))
return(list(dCdt = tran + reac))
}
# First order consumption of organic carbon (OC)
OC.model <- function (t = 0, OC, pars = NULL) {
tran <- tran.1D(C = OC, flux.up = F.OC, D = Db.grid,
v = v.grid, VF = svf.grid, dx = grid)$dC
reac <- - k*OC
return(list(dCdt = tran + reac))
}
#======================#
# Parameter definition #
#======================#
# Parameter values
F.OC <- 25 # input flux organic carbon [micromol cm-2 yr-1]
C.ow.O2 <- 0.25 # concentration O2 in overlying water [micromol cm-3]
por <- 0.8 # porosity
D <- 400 # diffusion coefficient O2 [cm2 yr-1]
Db <- 10 # mixing coefficient sediment [cm2 yr-1]
v <- 1 # advective velocity [cm yr-1]
k <- 1 # decay constant organic carbon [yr-1]
R.O2 <- 10 # O2 consumption rate [micromol cm-3 yr-1]
Ks <- 0.005 # O2 consumption saturation constant
# Grid definition
L <- 10 # depth of sediment domain [cm]
N <- 100 # number of grid layers
grid <- setup.grid.1D(x.up = 0, L = L, N = N)
# Volume fractions
por.grid <- setup.prop.1D(value = por, grid = grid)
svf.grid <- setup.prop.1D(value = (1-por), grid = grid)
D.grid <- setup.prop.1D(value = D, grid = grid)
Db.grid <- setup.prop.1D(value = Db, grid = grid)
v.grid <- setup.prop.1D(value = v, grid = grid)
#====================#
# Model solution #
#====================#
# Initial conditions + simulation O2
yini <- rep(0, length.out = N)
O2 <- steady.1D(y = yini, func = O2.model, nspec = 1)
# Initial conditions + simulation OC
yini <- rep(0, length.out = N)
OC <- steady.1D(y = yini, func = OC.model, nspec = 1)
# Plotting output, using S3 plot method of package rootSolve"
plot(O2, grid = grid$x.mid, xyswap = TRUE, main = "O2 concentration",
ylab = "depth [cm]", xlab = "", mfrow = c(1,2), type = "p", pch = 16)
plot(OC, grid = grid$x.mid, xyswap = TRUE, main = "C concentration",
ylab = "depth [cm]", xlab = "", mfrow = NULL)
## =============================================================================
## EXAMPLE 2: O2 in a cylindrical and spherical organism
## =============================================================================
# This example uses only the surface areas
# in the reactive transport term
#====================#
# Model formulation #
#====================#
# the numerical model - rate of change = transport-consumption
Cylinder.Model <- function(time, O2, pars)
return (list(
tran.1D(C = O2, C.down = BW, D = Da, A = A.cyl, dx = dx)$dC - Q
))
Sphere.Model <- function(time, O2, pars)
return (list(
tran.1D(C = O2, C.down = BW, D = Da, A = A.sphere, dx = dx)$dC - Q
))
#======================#
# Parameter definition #
#======================#
# parameter values
BW <- 2 # mmol/m3, oxygen conc in surrounding water
Da <- 0.5 # cm2/d effective diffusion coeff in organism
R <- 0.0025 # cm radius of organism
Q <- 250000 # nM/cm3/d oxygen consumption rate/ volume / day
L <- 0.05 # cm length of organism (if a cylinder)
# the numerical model
N <- 40 # layers in the body
dx <- R/N # thickness of each layer
x.mid <- seq(dx/2, by = dx, length.out = N) # distance of center to mid-layer
x.int <- seq(0, by = dx, length.out = N+1) # distance to layer interface
# Cylindrical surfaces
A.cyl <- 2*pi*x.int*L # surface at mid-layer depth
# Spherical surfaces
A.sphere <- 4*pi*x.int^2 # surface of sphere, at each mid-layer
#====================#
# Model solution #
#====================#
# the analytical solution of cylindrical and spherical model
cylinder <- function(Da, Q, BW, R, r) BW + Q/(4*Da)*(r^2-R^2)
sphere <- function(Da, Q, BW, R, r) BW + Q/(6*Da)*(r^2-R^2)
# solve the model numerically for a cylinder
O2.cyl <- steady.1D (y = runif(N), name = "O2",
func = Cylinder.Model, nspec = 1, atol = 1e-10)
# solve the model numerically for a sphere
O2.sphere <- steady.1D (y = runif(N), name = "O2",
func = Sphere.Model, nspec = 1, atol = 1e-10)
#====================#
# Plotting output #
#====================#
# Analytical solution - "observations"
Ana.cyl <- cbind(x.mid, O2 = cylinder(Da, Q, BW, R, x.mid))
Ana.spher <- cbind(x.mid, O2 = sphere(Da, Q, BW, R, x.mid))
plot(O2.cyl, O2.sphere, grid = x.mid, lwd = 2, lty = 1, col = 1:2,
xlab = "distance from centre, cm",
ylab = "mmol/m3", main = "tran.1D",
sub = "diffusion-reaction in a cylinder and sphere",
obs = list(Ana.cyl, Ana.spher), obspar = list(pch = 16, col =1:2))
legend ("topleft", lty = c(1, NA), pch = c(NA, 18),
c("numerical approximation", "analytical solution"))
legend ("bottomright", pch = 16, lty = 1, col = 1:2,
c("cylinder", "sphere"))
## =============================================================================
## EXAMPLE 3: O2 consumption in a spherical aggregate
## =============================================================================
# this example uses both the surface areas and the volume fractions
# in the reactive transport term
#====================#
# Model formulation #
#====================#
Aggregate.Model <- function(time, O2, pars) {
tran <- tran.1D(C = O2, C.down = C.ow.O2,
D = D.grid, A = A.grid,
VF = por.grid, dx = grid )$dC
reac <- - R.O2*(O2/(Ks+O2))*(O2>0)
return(list(dCdt = tran + reac, consumption = -reac))
}
#======================#
# Parameter definition #
#======================#
# Parameters
C.ow.O2 <- 0.25 # concentration O2 water [micromol cm-3]
por <- 0.8 # porosity
D <- 400 # diffusion coefficient O2 [cm2 yr-1]
v <- 0 # advective velocity [cm yr-1]
R.O2 <- 1000000 # O2 consumption rate [micromol cm-3 yr-1]
Ks <- 0.005 # O2 saturation constant [micromol cm-3]
# Grid definition
R <- 0.025 # radius of the agggregate [cm]
N <- 100 # number of grid layers
grid <- setup.grid.1D(x.up = 0, L = R, N = N)
# Volume fractions
por.grid <- setup.prop.1D(value = por, grid = grid)
D.grid <- setup.prop.1D(value = D, grid = grid)
# Surfaces
A.mid <- 4*pi*grid$x.mid^2 # surface of sphere at middle of grid cells
A.int <- 4*pi*grid$x.int^2 # surface of sphere at interface
A.grid <- list(int = A.int, mid = A.mid)
#====================#
# Model solution #
#====================#
# Numerical solution: staedy state
O2.agg <- steady.1D (runif(N), func = Aggregate.Model, nspec = 1,
atol = 1e-10, names = "O2")
#====================#
# Plotting output #
#====================#
par(mfrow = c(1,1))
plot(grid$x.mid, O2.agg$y, xlab = "distance from centre, cm",
ylab = "mmol/m3",
main = "Diffusion-reaction of O2 in a spherical aggregate")
legend ("bottomright", pch = c(1, 18), lty = 1, col = "black",
c("O2 concentration"))
# Similar, using S3 plot method of package rootSolve"
plot(O2.agg, grid = grid$x.mid, which = c("O2", "consumption"),
xlab = "distance from centre, cm", ylab = c("mmol/m3","mmol/m3/d"))
ReacTran documentation built on May 2, 2019, 9:38 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Note Author(s) References See Also Examples
Description
Estimates the transport term (i.e. the rate of change of a concentration due to diffusion and advection) in a one-dimensional model of a liquid (volume fraction constant and equal to one) or in a porous medium (volume fraction variable and lower than one).
The interfaces between grid cells can have a variable cross-sectional area, e.g. when modelling spherical or cylindrical geometries (see example).
Usage
1 2 3 4 5 |
Arguments
C |
concentration, expressed per unit of phase volume, defined at the centre of each grid cell. A vector of length N [M/L3] |
C.up |
concentration at upstream boundary. One value [M/L3] |
C.down |
concentration at downstream boundary. One value [M/L3] |
flux.up |
flux across the upstream boundary, positive = INTO model
domain. One value, expressed per unit of total surface [M/L2/T].
If |
flux.down |
flux across the downstream boundary, positive = OUT
of model domain. One value, expressed per unit of total surface [M/L2/T].
If |
a.bl.up |
convective transfer coefficient across the upstream
boundary layer. |
a.bl.down |
convective transfer coefficient across the downstream
boundary layer (L). |
D |
diffusion coefficient, defined on grid cell interfaces.
One value, a vector of length N+1 [L2/T], or a |
v |
advective velocity, defined on the grid cell
interfaces. Can be positive (downstream flow) or negative (upstream flow).
One value, a vector of length N+1 [L/T], or a |
AFDW |
weight used in the finite difference scheme for advection,
defined on grid cell interfaces; backward = 1, centred = 0.5, forward = 0;
default is backward. One value, a vector of length N+1, or a
|
VF |
Volume fraction defined at the grid cell interfaces. One value,
a vector of length N+1, or a |
A |
Interface area defined at the grid cell interfaces. One value,
a vector of length N+1, or a |
dx |
distance between adjacent cell interfaces (thickness of grid
cells). One value, a vector of length N, or a |
full.check |
logical flag enabling a full check of the consistency
of the arguments (default = |
full.output |
logical flag enabling a full return of the output
(default = |
Details
The boundary conditions are either
(1) zero-gradient.
(2) fixed concentration.
(3) convective boundary layer.
(4) fixed flux.
The above order also shows the priority. The default condition is the zero gradient. The fixed concentration condition overrules the zero gradient. The convective boundary layer condition overrules the fixed concentration and zero gradient. The fixed flux overrules all other specifications.
Ensure that the boundary conditions are well defined: for instance, it does not make sense to specify an influx in a boundary cell with the advection velocity pointing outward.
Transport properties:
The diffusion coefficient (D),
the advective velocity (v),
the volume fraction (VF), the interface surface (A),
and the advective finite difference weight (AFDW)
can either be specified as one value, a vector or a 1D property list
as generated by setup.prop.1D.
When a vector, this vector must be of length N+1, defined at all grid cell interfaces, including the upper and lower boundary.
The finite difference grid (dx) is specified either as
one value, a vector or a 1D grid list, as generated by setup.grid.1D.
Value
dC |
the rate of change of the concentration C due to transport, defined in the centre of each grid cell. The rate of change is expressed per unit of phase volume [M/L3/T] |
C.up |
concentration at the upstream interface. One value [M/L3]
only when ( |
C.down |
concentration at the downstream interface. One value [M/L3]
only when ( |
dif.flux |
diffusive flux across at the interface of each grid cell.
A vector of length N+1 [M/L2/T]
only when ( |
adv.flux |
advective flux across at the interface of each grid cell.
A vector of length N+1 [M/L2/T]
only when ( |
flux |
total flux across at the interface of each grid cell. A vector
of length N+1 [M/L2/T].
only when ( |
flux.up |
flux across the upstream boundary, positive = INTO model domain. One value [M/L2/T] |
flux.down |
flux across the downstream boundary, positive = OUT of model domain. One value [M/L2/T] |
Note
The advective equation is not checked for mass conservation. Sometimes, this is
not an issue, for instance when v represents a sinking velocity of
particles or a swimming velocity of organisms.
In others cases however, mass conservation needs to be accounted for.
To ensure mass conservation, the advective velocity must obey certain
continuity constraints: in essence the product of the volume fraction (VF),
interface surface area (A) and advective velocity (v) should be constant.
In sediments, one can use setup.compaction.1D to ensure that
the advective velocities for the pore water and solid phase meet these
constraints.
In terms of the units of concentrations and fluxes we follow the convention
in the geosciences.
The concentration C, C.up, C.down as well at the rate of
change of the concentration dC are always expressed per unit of
phase volume (i.e. per unit volume of solid or liquid).
Total concentrations (e.g. per unit volume of bulk sediment) can be obtained by multiplication with the appropriate volume fraction. In contrast, fluxes are always expressed per unit of total interface area (so here the volume fraction is accounted for).
Author(s)
Filip Meysman <filip.meysman@nioz.nl>, Karline Soetaert <karline.soetaert@nioz.nl>
References
Soetaert and Herman (2009). A practical guide to ecological modelling - using R as a simulation platform. Springer
See Also
tran.volume.1D for a discretisation the transport equation using finite volumes.
advection.1D, for more sophisticated advection schemes
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 | ## =============================================================================
## EXAMPLE 1: O2 and OC consumption in sediments
## =============================================================================
# this example uses only the volume fractions
# in the reactive transport term
#====================#
# Model formulation #
#====================#
# Monod consumption of oxygen (O2)
O2.model <- function (t = 0, O2, pars = NULL) {
tran <- tran.1D(C = O2, C.up = C.ow.O2, D = D.grid,
v = v.grid, VF = por.grid, dx = grid)$dC
reac <- - R.O2*(O2/(Ks+O2))
return(list(dCdt = tran + reac))
}
# First order consumption of organic carbon (OC)
OC.model <- function (t = 0, OC, pars = NULL) {
tran <- tran.1D(C = OC, flux.up = F.OC, D = Db.grid,
v = v.grid, VF = svf.grid, dx = grid)$dC
reac <- - k*OC
return(list(dCdt = tran + reac))
}
#======================#
# Parameter definition #
#======================#
# Parameter values
F.OC <- 25 # input flux organic carbon [micromol cm-2 yr-1]
C.ow.O2 <- 0.25 # concentration O2 in overlying water [micromol cm-3]
por <- 0.8 # porosity
D <- 400 # diffusion coefficient O2 [cm2 yr-1]
Db <- 10 # mixing coefficient sediment [cm2 yr-1]
v <- 1 # advective velocity [cm yr-1]
k <- 1 # decay constant organic carbon [yr-1]
R.O2 <- 10 # O2 consumption rate [micromol cm-3 yr-1]
Ks <- 0.005 # O2 consumption saturation constant
# Grid definition
L <- 10 # depth of sediment domain [cm]
N <- 100 # number of grid layers
grid <- setup.grid.1D(x.up = 0, L = L, N = N)
# Volume fractions
por.grid <- setup.prop.1D(value = por, grid = grid)
svf.grid <- setup.prop.1D(value = (1-por), grid = grid)
D.grid <- setup.prop.1D(value = D, grid = grid)
Db.grid <- setup.prop.1D(value = Db, grid = grid)
v.grid <- setup.prop.1D(value = v, grid = grid)
#====================#
# Model solution #
#====================#
# Initial conditions + simulation O2
yini <- rep(0, length.out = N)
O2 <- steady.1D(y = yini, func = O2.model, nspec = 1)
# Initial conditions + simulation OC
yini <- rep(0, length.out = N)
OC <- steady.1D(y = yini, func = OC.model, nspec = 1)
# Plotting output, using S3 plot method of package rootSolve"
plot(O2, grid = grid$x.mid, xyswap = TRUE, main = "O2 concentration",
ylab = "depth [cm]", xlab = "", mfrow = c(1,2), type = "p", pch = 16)
plot(OC, grid = grid$x.mid, xyswap = TRUE, main = "C concentration",
ylab = "depth [cm]", xlab = "", mfrow = NULL)
## =============================================================================
## EXAMPLE 2: O2 in a cylindrical and spherical organism
## =============================================================================
# This example uses only the surface areas
# in the reactive transport term
#====================#
# Model formulation #
#====================#
# the numerical model - rate of change = transport-consumption
Cylinder.Model <- function(time, O2, pars)
return (list(
tran.1D(C = O2, C.down = BW, D = Da, A = A.cyl, dx = dx)$dC - Q
))
Sphere.Model <- function(time, O2, pars)
return (list(
tran.1D(C = O2, C.down = BW, D = Da, A = A.sphere, dx = dx)$dC - Q
))
#======================#
# Parameter definition #
#======================#
# parameter values
BW <- 2 # mmol/m3, oxygen conc in surrounding water
Da <- 0.5 # cm2/d effective diffusion coeff in organism
R <- 0.0025 # cm radius of organism
Q <- 250000 # nM/cm3/d oxygen consumption rate/ volume / day
L <- 0.05 # cm length of organism (if a cylinder)
# the numerical model
N <- 40 # layers in the body
dx <- R/N # thickness of each layer
x.mid <- seq(dx/2, by = dx, length.out = N) # distance of center to mid-layer
x.int <- seq(0, by = dx, length.out = N+1) # distance to layer interface
# Cylindrical surfaces
A.cyl <- 2*pi*x.int*L # surface at mid-layer depth
# Spherical surfaces
A.sphere <- 4*pi*x.int^2 # surface of sphere, at each mid-layer
#====================#
# Model solution #
#====================#
# the analytical solution of cylindrical and spherical model
cylinder <- function(Da, Q, BW, R, r) BW + Q/(4*Da)*(r^2-R^2)
sphere <- function(Da, Q, BW, R, r) BW + Q/(6*Da)*(r^2-R^2)
# solve the model numerically for a cylinder
O2.cyl <- steady.1D (y = runif(N), name = "O2",
func = Cylinder.Model, nspec = 1, atol = 1e-10)
# solve the model numerically for a sphere
O2.sphere <- steady.1D (y = runif(N), name = "O2",
func = Sphere.Model, nspec = 1, atol = 1e-10)
#====================#
# Plotting output #
#====================#
# Analytical solution - "observations"
Ana.cyl <- cbind(x.mid, O2 = cylinder(Da, Q, BW, R, x.mid))
Ana.spher <- cbind(x.mid, O2 = sphere(Da, Q, BW, R, x.mid))
plot(O2.cyl, O2.sphere, grid = x.mid, lwd = 2, lty = 1, col = 1:2,
xlab = "distance from centre, cm",
ylab = "mmol/m3", main = "tran.1D",
sub = "diffusion-reaction in a cylinder and sphere",
obs = list(Ana.cyl, Ana.spher), obspar = list(pch = 16, col =1:2))
legend ("topleft", lty = c(1, NA), pch = c(NA, 18),
c("numerical approximation", "analytical solution"))
legend ("bottomright", pch = 16, lty = 1, col = 1:2,
c("cylinder", "sphere"))
## =============================================================================
## EXAMPLE 3: O2 consumption in a spherical aggregate
## =============================================================================
# this example uses both the surface areas and the volume fractions
# in the reactive transport term
#====================#
# Model formulation #
#====================#
Aggregate.Model <- function(time, O2, pars) {
tran <- tran.1D(C = O2, C.down = C.ow.O2,
D = D.grid, A = A.grid,
VF = por.grid, dx = grid )$dC
reac <- - R.O2*(O2/(Ks+O2))*(O2>0)
return(list(dCdt = tran + reac, consumption = -reac))
}
#======================#
# Parameter definition #
#======================#
# Parameters
C.ow.O2 <- 0.25 # concentration O2 water [micromol cm-3]
por <- 0.8 # porosity
D <- 400 # diffusion coefficient O2 [cm2 yr-1]
v <- 0 # advective velocity [cm yr-1]
R.O2 <- 1000000 # O2 consumption rate [micromol cm-3 yr-1]
Ks <- 0.005 # O2 saturation constant [micromol cm-3]
# Grid definition
R <- 0.025 # radius of the agggregate [cm]
N <- 100 # number of grid layers
grid <- setup.grid.1D(x.up = 0, L = R, N = N)
# Volume fractions
por.grid <- setup.prop.1D(value = por, grid = grid)
D.grid <- setup.prop.1D(value = D, grid = grid)
# Surfaces
A.mid <- 4*pi*grid$x.mid^2 # surface of sphere at middle of grid cells
A.int <- 4*pi*grid$x.int^2 # surface of sphere at interface
A.grid <- list(int = A.int, mid = A.mid)
#====================#
# Model solution #
#====================#
# Numerical solution: staedy state
O2.agg <- steady.1D (runif(N), func = Aggregate.Model, nspec = 1,
atol = 1e-10, names = "O2")
#====================#
# Plotting output #
#====================#
par(mfrow = c(1,1))
plot(grid$x.mid, O2.agg$y, xlab = "distance from centre, cm",
ylab = "mmol/m3",
main = "Diffusion-reaction of O2 in a spherical aggregate")
legend ("bottomright", pch = c(1, 18), lty = 1, col = "black",
c("O2 concentration"))
# Similar, using S3 plot method of package rootSolve"
plot(O2.agg, grid = grid$x.mid, which = c("O2", "consumption"),
xlab = "distance from centre, cm", ylab = c("mmol/m3","mmol/m3/d"))
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.