In dynamic-R/RTM: Reaction Transport Modelling: Course Tutorials and Exercises
Answers
Mass balance equations
The mass balances for ammonia ($NH_3$), nitrate ($NO_3^-$) and oxygen ($O_2$) can be written as the following differential equations:
$$\frac{\partial [NO_3^-]}{\partial t} = D\cdot \frac{\partial^2 [NO_3^-]}{\partial x^2 } - v\cdot \frac{\partial [NO_3^-]}{\partial x } + Nitrification$$
$$\frac{\partial [NH_3]}{\partial t} = D\cdot \frac{\partial^2 [NH_3]}{\partial x^2 } - v\cdot \frac{\partial [NH_3]}{\partial x } - Nitrification$$
$$\frac{\partial [O_2]}{\partial t} = D\cdot \frac{\partial^2 [O_2]}{\partial x^2 } - v\cdot \frac{\partial [O_2]}{\partial x } + Aeration - 2 \cdot Nitrification$$
Here, $D$ is the dispersion coefficient due to tidal mixing, and $v$ is the advective velocity (in $m~s^{-1}$), calculated from the discharge (Q) and cross-sectional area (A) according to $v=Q/A$. Thus, the 1st and 2nd term in the differential equations describe the rate of change due to transport by dispersion (diffusion-like) and advection, respectively. Note that the parameters characterizing the magnitude of this change are the same for each component, which is because all components are dissolved substances affected in the same way by the water flow. Also note that we assumed that $D$ and $u$ are not varying along the estuary, which is why we could take them out of the spatial derivatives.
The boundary conditions at the upper boundary (river side): $[O_2]{x=0} = 0.1$, $[NO_3^-]{x=0} = 0.3$, $[NH_3]{x=0} = 0.1$.
The boundary conditions at the lower boundary (sea side): $[O_2]{x=1e5} = 0.3$, $[NO_3^-]{x=1e5} = 0.050$, $[NH_3]{x=1e5} = 0.010$. All values are in $mol~m^{-3}$.
It is not necessary to model the proton ($H^+$) or water ($H_2O$) concentrations in this exercise.
R implementation
Note how the vector of state variables is defined in models with transport. From zero-dimensional models (models without transport), we are used to define state variables in a vector with named components, for example like this:
# state variables in a zero-dimensional models
state.ini <- c(O2=0.1, NO3=0.3, NH3=0.1)
This made it easy to refer to the state variables in the model function using the \texttt{with(as.list(...} function.
In models with transport, each state variable is a vector by itself! Thus, this "trick" is no longer possible. Instead, we concatenate the vectors for each state variable into a long vector called state, and then "unpack" the individual vectors from this variable within the model function. We will assign the initial values of the state variables into a vector called state.ini.
require(ReacTran) # package with solution methods
# Note: units are: m, days, mol/m3.
# model grid
Length <- 100000 # m, domain length
N <- 500 # -, number of boxes
Grid <- setup.grid.1D(L = Length, N = N) # grid of equally sized boxes
Area <- 20000 # m2, cross-sectional area of the estuary
day2sec <- 24*3600 # number of seconds in a day
# State variables, each is a vector of length=N - initial concentrations set to 0
Oxygen <- rep(0, times = N)
Nitrate <- rep(0, times = N)
Ammonia <- rep(0, times = N)
state.ini <- c(Oxygen, Nitrate, Ammonia)
SVnames <- c("Oxygen", "Nitrate", "Ammonia")
# model parameters
pars <- c(
riverO2 = 0.1, # river oxygen conc [mol/m3]
seaO2 = 0.3, # marine oxygen conc [mol/m3]
riverNit = 0.3, # river nitrate conc [mol/m3]
seaNit = 0.05, # marine nitrate conc
riverAmm = 0.1, # river ammonium conc
seaAmm = 0.01, # marine ammonium conc
depth = 10, # [m]
v = 100/Area*day2sec, # advection velocity [m/d], discharge/cross-sect. area
Ddisp = 350*day2sec, # dispersion coefficient [m2/d], tidal dispersion
rnitri = 0.1, # nitrification rate constant, [/d]
ksO2 = 1e-3, # Monod ct for O2 limitation of nitrification [mol/m3]
piston = 1.0, # piston velocity [m/day]
O2sat = 0.3 # saturated oxygen concentration [mol/m3], solubility
)
# Model function
Scheldt1D <- function(t, state, parms) { # state is a LONG vector of state variables at time t
with(as.list(parms),{
# here we "unpack" the vectors of individual state variables from the LONG input vector
Oxygen <- state[(0*N+1) : (1*N)] # note parentheses around (0*N+1) and (1*N)!!
Nitrate <- state[(1*N+1) : (2*N)]
Ammonia <- state[(2*N+1) : (3*N)]
# Transport - tran.1D solves the "transport terms" (spatial derivatives).
# It also implements the boundary conditions!
TranOxygen <- tran.1D(C = Oxygen,
C.up = riverO2, # upstream concentration
C.down = seaO2, # downstream concentration
D = Ddisp, v = v, dx = Grid)
TranNitrate <- tran.1D(C = Nitrate,
C.up = riverNit, C.down = seaNit,
D = Ddisp, v = v, dx = Grid)
TranAmmonia <- tran.1D(C = Ammonia,
C.up = riverAmm, C.down = seaAmm,
D = Ddisp, v = v, dx = Grid)
# Nitrification rate (calculated in every location)
Nitrification <- rnitri * Ammonia * Oxygen/(Oxygen+ksO2)
# air-water exchange rate (calculated in every location)
Aeration <- -piston/depth * (Oxygen - O2sat)
# the rates of change = transport + net reaction
dNitrate.dt <- TranNitrate$dC + Nitrification
dAmmonia.dt <- TranAmmonia$dC - Nitrification
dOxygen.dt <- TranOxygen$dC - 2*Nitrification + Aeration
return(list(c(dOxygen.dt, dNitrate.dt, dAmmonia.dt), # the rates of change
TotalNitrogen = mean(Nitrate + Ammonia),
# process rates per volume (mol m-3 d-1):
Nitrification = Nitrification,
Aeration = Aeration,
# import/export rates at the domain boundaries (mol d-1):
NH3import = TranAmmonia$flux.up * Area,
NO3import = TranNitrate$flux.up * Area,
O2import = TranOxygen$flux.up * Area,
NH3export = TranAmmonia$flux.down * Area,
NO3export = TranNitrate$flux.down * Area,
O2export = TranOxygen$flux.down * Area,
# total import/export rates at the domain boundaries (mol d-1):
TotNimport = (TranAmmonia$flux.up +
TranNitrate$flux.up) * Area,
TotNexport = (TranAmmonia$flux.down +
TranNitrate$flux.down) * Area,
# integrated rates over the volume of the modeled domain (mol d-1):
TotalNitrification = sum(Nitrification * Area * Grid$dx),
TotalAeration = sum(Aeration * Area * Grid$dx)
))
})
}
Comparison of the 1970s and the current situation
We estimate the steady-state solution for the current situation and for the 1970s.
Scheldt00 <- steady.1D(y = state.ini, parms = pars, func = Scheldt1D,
positive = TRUE, nspec = length(SVnames),
dimens = N, names = SVnames)
# The parameters for the 1970s
par70 <- pars
par70["riverNit"] <- 0.05 # river nitrate concentration
par70["seaNit"] <- 0.02 # marine nitrate concentration
par70["riverAmm"] <- 0.60 # river ammonium concentration
par70["seaAmm"] <- 0.02 # marine ammonium concentration
Scheldt70 <- steady.1D(y = state.ini, parms = par70, func = Scheldt1D,
positive = TRUE, nspec = length(SVnames),
dimens = N, names = SVnames)
plot(Scheldt00, Scheldt70, grid=Grid$x.mid/1000, lwd=2, lty=1,
mfrow=c(2,3), ylab="mol/m3", xlab="distance (km)")
plot(Scheldt00, Scheldt70, grid=Grid$x.mid/1000, lwd=2, lty=1, mfrow=NULL,
which = "Nitrification", ylab="molN/m3/d", xlab="distance (km)")
plot(Scheldt00, Scheldt70, grid=Grid$x.mid/1000, lwd=2, lty=1, mfrow=NULL,
which = "Aeration", ylab="molO2/m3/d", xlab="distance (km)")
legend("topright", legend = c("2000", "1970s"), lty=1, lwd=2, col=1:2, bty="n")
Budget in Mmol/year
The current model function returns the rates in $mol~d^{-1}$. To construct the annual budget in $Mmol~yr^{-1}$, we use the conversion factor $365/1e6$. The budget values can be added to a diagram that represents the relevant state variables for the entire estuary as boxes, similar to the one shown in Figure 1.
f1 <- 365/1e6 # from mol d-1 to 10^6 mol yr-1
toselect <- c("NH3import", "NO3import", "O2import",
"NH3export", "NO3export", "O2export",
"TotalNitrification", "TotalAeration" )
BUDGET <- data.frame(Megamol_per_yr_2000 = unlist(Scheldt00[toselect]),
Megamol_per_yr_1970s = unlist(Scheldt70[toselect]))*f1
knitr::kable(BUDGET, digits = 0)
The impact of ammonia concentration in inflowing waters
The sensitivity analysis is run using a for-loop, where the minimal value in the steady state solution is found for each run of the loop.
Sens <- function(NH3up){
parsNH3 <- pars
parsNH3["riverAmm"] <- NH3up # river ammonium conc
ScheldtNH3 <- steady.1D(y = state.ini, parms = parsNH3, func = Scheldt1D,
positive = TRUE, nspec = length(SVnames),
dimens = N, names = SVnames)
return(min(ScheldtNH3$y[,"Oxygen"]))
}
NH3_vect <- seq(from = 0, to = 0.8, by = 0.04)
O2_vect <- NULL # initialize to an "empty vector"
# use a for-loop to find the minimal O2 concentration for each value in NH3_vect
for (NH3up in NH3_vect)
O2_vect <- c(O2_vect, Sens(NH3up)) # add a value in every step of the loop
plot(NH3_vect, O2_vect, type="b", xlab = "upstream NH3 concentration, [molN/m3]",
ylab = "O2, [molO2/m3]", main = "Minimal O2 concentration in the estuary")
dynamic-R/RTM documentation built on Feb. 28, 2025, 1:23 p.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.
Answers
Mass balance equations
The mass balances for ammonia ($NH_3$), nitrate ($NO_3^-$) and oxygen ($O_2$) can be written as the following differential equations:
$$\frac{\partial [NO_3^-]}{\partial t} = D\cdot \frac{\partial^2 [NO_3^-]}{\partial x^2 } - v\cdot \frac{\partial [NO_3^-]}{\partial x } + Nitrification$$ $$\frac{\partial [NH_3]}{\partial t} = D\cdot \frac{\partial^2 [NH_3]}{\partial x^2 } - v\cdot \frac{\partial [NH_3]}{\partial x } - Nitrification$$ $$\frac{\partial [O_2]}{\partial t} = D\cdot \frac{\partial^2 [O_2]}{\partial x^2 } - v\cdot \frac{\partial [O_2]}{\partial x } + Aeration - 2 \cdot Nitrification$$
Here, $D$ is the dispersion coefficient due to tidal mixing, and $v$ is the advective velocity (in $m~s^{-1}$), calculated from the discharge (Q) and cross-sectional area (A) according to $v=Q/A$. Thus, the 1st and 2nd term in the differential equations describe the rate of change due to transport by dispersion (diffusion-like) and advection, respectively. Note that the parameters characterizing the magnitude of this change are the same for each component, which is because all components are dissolved substances affected in the same way by the water flow. Also note that we assumed that $D$ and $u$ are not varying along the estuary, which is why we could take them out of the spatial derivatives.
The boundary conditions at the upper boundary (river side): $[O_2]{x=0} = 0.1$, $[NO_3^-]{x=0} = 0.3$, $[NH_3]{x=0} = 0.1$. The boundary conditions at the lower boundary (sea side): $[O_2]{x=1e5} = 0.3$, $[NO_3^-]{x=1e5} = 0.050$, $[NH_3]{x=1e5} = 0.010$. All values are in $mol~m^{-3}$.
It is not necessary to model the proton ($H^+$) or water ($H_2O$) concentrations in this exercise.
R implementation
Note how the vector of state variables is defined in models with transport. From zero-dimensional models (models without transport), we are used to define state variables in a vector with named components, for example like this:
# state variables in a zero-dimensional models state.ini <- c(O2=0.1, NO3=0.3, NH3=0.1)
This made it easy to refer to the state variables in the model function using the \texttt{with(as.list(...} function.
In models with transport, each state variable is a vector by itself! Thus, this "trick" is no longer possible. Instead, we concatenate the vectors for each state variable into a long vector called state, and then "unpack" the individual vectors from this variable within the model function. We will assign the initial values of the state variables into a vector called state.ini.
require(ReacTran) # package with solution methods # Note: units are: m, days, mol/m3. # model grid Length <- 100000 # m, domain length N <- 500 # -, number of boxes Grid <- setup.grid.1D(L = Length, N = N) # grid of equally sized boxes Area <- 20000 # m2, cross-sectional area of the estuary day2sec <- 24*3600 # number of seconds in a day # State variables, each is a vector of length=N - initial concentrations set to 0 Oxygen <- rep(0, times = N) Nitrate <- rep(0, times = N) Ammonia <- rep(0, times = N) state.ini <- c(Oxygen, Nitrate, Ammonia) SVnames <- c("Oxygen", "Nitrate", "Ammonia") # model parameters pars <- c( riverO2 = 0.1, # river oxygen conc [mol/m3] seaO2 = 0.3, # marine oxygen conc [mol/m3] riverNit = 0.3, # river nitrate conc [mol/m3] seaNit = 0.05, # marine nitrate conc riverAmm = 0.1, # river ammonium conc seaAmm = 0.01, # marine ammonium conc depth = 10, # [m] v = 100/Area*day2sec, # advection velocity [m/d], discharge/cross-sect. area Ddisp = 350*day2sec, # dispersion coefficient [m2/d], tidal dispersion rnitri = 0.1, # nitrification rate constant, [/d] ksO2 = 1e-3, # Monod ct for O2 limitation of nitrification [mol/m3] piston = 1.0, # piston velocity [m/day] O2sat = 0.3 # saturated oxygen concentration [mol/m3], solubility ) # Model function Scheldt1D <- function(t, state, parms) { # state is a LONG vector of state variables at time t with(as.list(parms),{ # here we "unpack" the vectors of individual state variables from the LONG input vector Oxygen <- state[(0*N+1) : (1*N)] # note parentheses around (0*N+1) and (1*N)!! Nitrate <- state[(1*N+1) : (2*N)] Ammonia <- state[(2*N+1) : (3*N)] # Transport - tran.1D solves the "transport terms" (spatial derivatives). # It also implements the boundary conditions! TranOxygen <- tran.1D(C = Oxygen, C.up = riverO2, # upstream concentration C.down = seaO2, # downstream concentration D = Ddisp, v = v, dx = Grid) TranNitrate <- tran.1D(C = Nitrate, C.up = riverNit, C.down = seaNit, D = Ddisp, v = v, dx = Grid) TranAmmonia <- tran.1D(C = Ammonia, C.up = riverAmm, C.down = seaAmm, D = Ddisp, v = v, dx = Grid) # Nitrification rate (calculated in every location) Nitrification <- rnitri * Ammonia * Oxygen/(Oxygen+ksO2) # air-water exchange rate (calculated in every location) Aeration <- -piston/depth * (Oxygen - O2sat) # the rates of change = transport + net reaction dNitrate.dt <- TranNitrate$dC + Nitrification dAmmonia.dt <- TranAmmonia$dC - Nitrification dOxygen.dt <- TranOxygen$dC - 2*Nitrification + Aeration return(list(c(dOxygen.dt, dNitrate.dt, dAmmonia.dt), # the rates of change TotalNitrogen = mean(Nitrate + Ammonia), # process rates per volume (mol m-3 d-1): Nitrification = Nitrification, Aeration = Aeration, # import/export rates at the domain boundaries (mol d-1): NH3import = TranAmmonia$flux.up * Area, NO3import = TranNitrate$flux.up * Area, O2import = TranOxygen$flux.up * Area, NH3export = TranAmmonia$flux.down * Area, NO3export = TranNitrate$flux.down * Area, O2export = TranOxygen$flux.down * Area, # total import/export rates at the domain boundaries (mol d-1): TotNimport = (TranAmmonia$flux.up + TranNitrate$flux.up) * Area, TotNexport = (TranAmmonia$flux.down + TranNitrate$flux.down) * Area, # integrated rates over the volume of the modeled domain (mol d-1): TotalNitrification = sum(Nitrification * Area * Grid$dx), TotalAeration = sum(Aeration * Area * Grid$dx) )) }) }
Comparison of the 1970s and the current situation
We estimate the steady-state solution for the current situation and for the 1970s.
Scheldt00 <- steady.1D(y = state.ini, parms = pars, func = Scheldt1D, positive = TRUE, nspec = length(SVnames), dimens = N, names = SVnames) # The parameters for the 1970s par70 <- pars par70["riverNit"] <- 0.05 # river nitrate concentration par70["seaNit"] <- 0.02 # marine nitrate concentration par70["riverAmm"] <- 0.60 # river ammonium concentration par70["seaAmm"] <- 0.02 # marine ammonium concentration Scheldt70 <- steady.1D(y = state.ini, parms = par70, func = Scheldt1D, positive = TRUE, nspec = length(SVnames), dimens = N, names = SVnames)
plot(Scheldt00, Scheldt70, grid=Grid$x.mid/1000, lwd=2, lty=1, mfrow=c(2,3), ylab="mol/m3", xlab="distance (km)") plot(Scheldt00, Scheldt70, grid=Grid$x.mid/1000, lwd=2, lty=1, mfrow=NULL, which = "Nitrification", ylab="molN/m3/d", xlab="distance (km)") plot(Scheldt00, Scheldt70, grid=Grid$x.mid/1000, lwd=2, lty=1, mfrow=NULL, which = "Aeration", ylab="molO2/m3/d", xlab="distance (km)") legend("topright", legend = c("2000", "1970s"), lty=1, lwd=2, col=1:2, bty="n")
Budget in Mmol/year
The current model function returns the rates in $mol~d^{-1}$. To construct the annual budget in $Mmol~yr^{-1}$, we use the conversion factor $365/1e6$. The budget values can be added to a diagram that represents the relevant state variables for the entire estuary as boxes, similar to the one shown in Figure 1.
f1 <- 365/1e6 # from mol d-1 to 10^6 mol yr-1 toselect <- c("NH3import", "NO3import", "O2import", "NH3export", "NO3export", "O2export", "TotalNitrification", "TotalAeration" ) BUDGET <- data.frame(Megamol_per_yr_2000 = unlist(Scheldt00[toselect]), Megamol_per_yr_1970s = unlist(Scheldt70[toselect]))*f1 knitr::kable(BUDGET, digits = 0)
The impact of ammonia concentration in inflowing waters
The sensitivity analysis is run using a for-loop, where the minimal value in the steady state solution is found for each run of the loop.
Sens <- function(NH3up){ parsNH3 <- pars parsNH3["riverAmm"] <- NH3up # river ammonium conc ScheldtNH3 <- steady.1D(y = state.ini, parms = parsNH3, func = Scheldt1D, positive = TRUE, nspec = length(SVnames), dimens = N, names = SVnames) return(min(ScheldtNH3$y[,"Oxygen"])) } NH3_vect <- seq(from = 0, to = 0.8, by = 0.04) O2_vect <- NULL # initialize to an "empty vector" # use a for-loop to find the minimal O2 concentration for each value in NH3_vect for (NH3up in NH3_vect) O2_vect <- c(O2_vect, Sens(NH3up)) # add a value in every step of the loop
plot(NH3_vect, O2_vect, type="b", xlab = "upstream NH3 concentration, [molN/m3]", ylab = "O2, [molO2/m3]", main = "Minimal O2 concentration in the estuary")
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.