Description Usage Arguments Details Value Note Author(s) References See Also Examples
Estimates the advection term in a onedimensional 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 crosssectional area, e.g. when modelling spherical or cylindrical geometries (see example).
TVD (total variation diminishing) slope limiters ensure monotonic and positive schemes in the presence of strong gradients.
advection.1D
: uses finite differences.
This implies the use of velocity (length per time) and fluxes (mass per unit of area per unit of time).
advection.volume.1D
Estimates the volumetric advection term in a onedimensional model
of an aquatic system (river, estuary). This routine is particularly
suited for modelling channels (like rivers, estuaries) where the
crosssectional area changes, and hence the velocity changes.
Volumetric transport implies the use of flows (mass per unit of time).
When solved dynamically, the euler method should be used, unless the firstorder upstream method is used.
1 2 3 4 5 6 7 8 9  advection.1D(C, C.up = NULL, C.down = NULL,
flux.up = NULL, flux.down = NULL, v, VF = 1, A = 1, dx,
dt.default = 1, adv.method = c("muscl", "super", "quick", "p3", "up"),
full.check = FALSE)
advection.volume.1D(C, C.up = C[1], C.down = C[length(C)],
F.up = NULL, F.down = NULL, flow, V,
dt.default = 1, adv.method = c("muscl", "super", "quick", "p3", "up"),
full.check = FALSE)

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]. If 
C.down 
concentration at downstream boundary. One value [M/L3]. If 
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 
F.up 
total input across the upstream boundary, positive = INTO model
domain; used with 
F.down 
total input across the downstream boundary, positive = OUT
of model domain; used with 
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 
flow 
water flow rate, defined on grid cell interfaces.
One value, a vector of length N+1, or a list as defined by

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 
dt.default 
timestep to be used, if it cannot be estimated (e.g. when calculating steadystate conditions. 
V 
volume of cells. One value, or a vector of length N [L^3]. 
adv.method 
the advection method, slope limiter used to reduce the numerical dispersion. One of "quick","muscl","super","p3","up"  see details. 
full.check 
logical flag enabling a full check of the consistency
of the arguments (default = 
This implementation is based on the GOTM code
The boundary conditions are either
zerogradient.
fixed concentration.
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 fixed flux overrules the 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 advective velocity (v
),
the volume fraction (VF), and the interface surface (A
),
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
.
Several slope limiters are implemented to obtain monotonic and positive schemes also in the presence of strong gradients, i.e. to reduce the effect of numerical dispersion. The methods are (Pietrzak, 1989, Hundsdorfer and Verwer, 2003):
"quick": thirdorder scheme (TVD) with ULTIMATE QUICKEST limiter (quadratic upstream interpolation for convective kinematics with estimated stream terms) (Leonard, 1988)
"muscl": thirdorder scheme (TVD) with MUSCL limiter (monotonic upstream centered schemes for conservation laws) (van Leer, 1979).
"super": thirdorder scheme (TVD) with Superbee limiter (method=Superbee) (Roe, 1985)
"p3": thirdorder upstreambiased polynomial scheme (method=P3)
"up": firstorder upstream ( method=UPSTREAM)  this is the same method as implemented in tran.1D or tran.volume.1D
where "TVD" means a total variation diminishing scheme
Some schemes may produce artificial steepening. Scheme "p3" is not necessarily monotone (may produce negative concentrations!).
If during a certain time step the maximum Courant number is larger than one, a split iteration will be carried out which guarantees that the split step Courant numbers are just smaller than 1. The maximal number of such iterations is set to 100.
These limiters are supposed to work with explicit methods (euler). However, they will also work with implicit methods, although less effectively. Integrate ode.1D only if the model is stiff (see first example).
dC 
the rate of change of the concentration C due to advective transport, defined in the centre of each grid cell. The rate of change is expressed per unit of (phase) volume [M/L^3/T]. 
adv.flux 
advective flux across at the interface of each grid cell.
A vector of length N+1 [M/L2/T]  only for 
flux.up 
flux across the upstream boundary, positive = INTO model
domain. One value [M/L2/T]  only for 
flux.down 
flux across the downstream boundary, positive = OUT of
model domain. One value [M/L2/T]  only for 
adv.F 
advective mass flow across at the interface of each grid cell.
A vector of length N+1 [M/T]  only for 
F.up 
mass flow across the upstream boundary, positive = INTO model
domain. One value [M/T]  only for 
F.down 
flux across the downstream boundary, positive = OUT of
model domain. One value [M/T]  only for 
it 
number of split time iterations that were necessary. 
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).
Karline Soetaert <karline.soetaert@nioz.nl>
Pietrzak J (1998) The use of TVD limiters for forwardintime upstreambiased advection schemes in ocean modeling. Monthly Weather Review 126: 812 .. 830
Hundsdorfer W and Verwer JG (2003) Numerical Solution of TimeDependent AdvectionDiffusionReaction Equations. Springer Series in Computational Mathematics, SpringerVerlag, Berlin, 471 pages
Burchard H, Bolding K, Villarreal MR (1999) GOTM, a general ocean turbulence model. Theory, applications and test cases. Tech Rep EUR 18745 EN. European Commission
Leonard BP (1988) Simple high accuracy resolution program for convective modeling of discontinuities. Int. J. Numer. Meth.Fluids 8: 1291–1318.
Roe PL (1985) Some contributions to the modeling of discontinuous flows. Lect. Notes Appl. Math. 22: 163193.
van Leer B. (1979) Towards the ultimate conservative difference scheme V. A second order sequel to Godunov's method. J. Comput. Phys. 32: 101136
tran.1D
, for a discretisation of the general transport equations
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  ## =============================================================================
## EXAMPLE 1: Testing the various methods  moving a square pulse
## use of advection.1D
## The tests as in Pietrzak
## =============================================================================
##
# Model formulation #
##
model < function (t, y, parms,...) {
adv < advection.1D(y, v = v, dx = dx,
C.up = y[n], C.down = y[1], ...) # out on one side > in at other
return(list(adv$dC))
}
##
# Parameters #
##
n < 100
dx < 100/n
y < c(rep(1, 5), rep(2, 20), rep(1, n25))
v < 2
times < 0:300 # 3 times out and in
##
# model solution #
##
## a plotting function
plotfun < function (Out, ...) {
plot(Out[1, 1], type = "l", col = "red", ylab = "y", xlab = "x", ...)
lines(Out[nrow(Out), 2:(1+n)])
}
# courant number = 2
pm < par(mfrow = c(2, 2))
## third order TVD, quickest limiter
out < ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
method = "euler", nspec = 1, adv.method = "quick")
plotfun(out, main = "quickest, euler")
## thirdorder ustreambiased polynomial
out2 < ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
method = "euler", nspec = 1, adv.method = "p3")
plotfun(out2, main = "p3, euler")
## third order TVD, superbee limiter
out3 < ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
method = "euler", nspec = 1, adv.method = "super")
plotfun(out3, main = "superbee, euler")
## third order TVD, muscl limiter
out4 < ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
method = "euler", nspec = 1, adv.method = "muscl")
plotfun(out4, main = "muscl, euler")
## =============================================================================
## upstream, different timesteps , i.e. different courant number
## =============================================================================
out < ode.1D(y = y, times = times, func = model, parms = 0,
method = "euler", nspec = 1, adv.method = "up")
plotfun(out, main = "upstream, courant number = 2")
out2 < ode.1D(y = y, times = times, func = model, parms = 0, hini = 0.5,
method = "euler", nspec = 1, adv.method = "up")
plotfun(out2, main = "upstream, courant number = 1")
## Now muscl scheme, velocity against xaxis
y < rev(c(rep(0, 5), rep(1, 20), rep(0, n25)))
v < 2.0
out6 < ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
method = "euler", nspec = 1, adv.method = "muscl")
plotfun(out6, main = "muscl, reversed velocity, , courant number = 1")
image(out6, mfrow = NULL)
par(mfrow = pm)
## =============================================================================
## EXAMPLE 2: moving a square pulse in a widening river
## use of advection.volume.1D
## =============================================================================
##
# Model formulation #
##
river.model < function (t=0, C, pars = NULL, ...) {
tran < advection.volume.1D(C = C, C.up = 0,
flow = flow, V = Volume,...)
return(list(dCdt = tran$dC, F.down = tran$F.down, F.up = tran$F.up))
}
##
# Parameters #
##
# Initialising morphology river:
nbox < 100 # number of grid cells
lengthRiver < 100000 # [m]
BoxLength < lengthRiver / nbox # [m]
Distance < seq(BoxLength/2, by = BoxLength, len = nbox) # [m]
# Cross sectional area: sigmoid function of distance [m2]
CrossArea < 4000 + 72000 * Distance^5 /(Distance^5+50000^5)
# Volume of boxes (m3)
Volume < CrossArea*BoxLength
# Transport coefficients
flow < 1000*24*3600 # m3/d, main river upstream inflow
##
# Model solution #
##
pm < par(mfrow=c(2,2))
# a square pulse
yini < c(rep(10, 10), rep(0, nbox10))
## third order TVD, muscl limiter
Conc < ode.1D(y = yini, fun = river.model, method = "euler", hini = 1,
parms = NULL, nspec = 1, times = 0:40, adv.method = "muscl")
image(Conc, main = "muscl", mfrow = NULL)
plot(Conc[30, 2:(1+nbox)], type = "l", lwd = 2, xlab = "x", ylab = "C",
main = "muscl after 30 days")
## simple upstream differencing
Conc2< ode.1D(y = yini, fun = river.model, method = "euler", hini = 1,
parms = NULL, nspec = 1, times = 0:40, adv.method = "up")
image(Conc2, main = "upstream", mfrow = NULL)
plot(Conc2[30, 2:(1+nbox)], type = "l", lwd = 2, xlab = "x", ylab = "C",
main = "upstream after 30 days")
par(mfrow = pm)
# Note: the more sophisticated the scheme, the more mass lost/created
# increase tolerances to improve this.
CC < Conc[ , 2:(1+nbox)]
MASS < t(CC)*Volume
colSums(MASS)
## =============================================================================
## EXAMPLE 3: A steadystate solution
## use of advection.volume.1D
## =============================================================================
Sink < function (t, y, parms, ...) {
C1 < y[1:N]
C2 < y[(N+1):(2*N)]
C3 < y[(2*N+1):(3*N)]
# Rate of change= Flux gradient and firstorder consumption
# upstream can be implemented in two ways:
dC1 < advection.1D(C1, v = sink, dx = dx,
C.up = 100, adv.method = "up", ...)$dC  decay*C1
# same, using tran.1D
# dC1 < tran.1D(C1, v = sink, dx = dx,
# C.up = 100, D = 0)$dC 
# decay*C1
dC2 < advection.1D(C2, v = sink, dx = dx,
C.up = 100, adv.method = "p3", ...)$dC 
decay*C2
dC3 < advection.1D(C3, v = sink, dx = dx,
C.up = 100, adv.method = "muscl", ...)$dC 
decay*C3
list(c(dC1, dC2, dC3))
}
N < 10
L < 1000
dx < L/N # thickness of boxes
sink < 10
decay < 0.1
out < steady.1D(runif(3*N), func = Sink, names = c("C1", "C2", "C3"),
parms = NULL, nspec = 3, bandwidth = 2)
matplot(out$y, 1:N, type = "l", ylim = c(10, 0), lwd = 2,
main = "Steadystate")
legend("bottomright", col = 1:3, lty = 1:3,
c("upstream", "p3", "muscl"))

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.