k-HLR relationships

The purpose of this vignette is to explain a way to derive “equivalent”k-HLR relationships that let you translate between Dynamic DMSTA-style P model parameters (ie. \(K_1\), \(K_2\), \(K_3\),\(C_0\) and \(C_1\)), and the familiar \(k-C^*\)/N-tank-in-series design relationships used for screening and sizing of treatment wetlands.

DMSTA was conceptualized as a dynamic extension of a steady state settling rate model. In practice, DMSTA maintains dynamic water balance, hydraulic routing and P cycling routines. Generally, P removal is proportional to the product of water column P and available P storage in vegetation and sediment. Meanwhile, P recycling is described mathematically as a quadratic to the liable pool. Finally P removal to permanent burial is 1st order from labile storage.

DMSTA has three key rate constants, \(K_1\), \(K_2\) and \(K_3\) that scale P biogeochemistry represented by the model (removal, recycle and burial). These are the principal fitting parameters for calibration to observed data. These parameters were linearly transformed to a different parameter space from \(K\), \(C_0\) and \(C_1\).

Typically, \(K\) is estimated by fitting a model to a concentration time-series, \(C_0\) and \(C_1\) are estimated from the storage and concentration relationship for a given wetlands and/or vegetative community.


The Models:

  1. The \(NKC^*\) Model: A generalized tank-in-series outflow model represented as
\[C_{out} = \frac{C_{in} - C^{*}}{(1 + \frac{k}{N \cdot HLR})^{N}} + C^*\]
Model variables
Variable Definition Units
\(Q\) Flow \(\mathrm{m^{3}\ y^{-1}}\)
\(A\) Area \(\mathrm{m^{2}}\)
\(HLR\) Hydraulic Looading Rate \(HLR = Q/A\) \(\mathrm{m\ y^{-1}}\)
\(C_{in}\),\(C_{out}\) Inflow/Outflow Concentration \(\mathrm{mg\ m^{-3}\ or \ \mu g\ L^{-1}}\)
\(C^*\) Background concentration \(\mathrm{mg\ m^{-3}\ or \ \mu g\ L^{-1}}\)
\(k\) areal net settling (removal velocity) \(\mathrm{m\ y^{-1}}\)
\(N\) Number of tanks in series unitless

When \(N = 1\) it indicates a fully mixed tank (conservative), meanwhile \(N = 3 - 5\) is the common STA design range.

Using equation above, here are some useful limited cases and/or simplifications of the model.

\[C_{out} = C^{*} + \frac{C_{in} - C^{*}}{1 + \frac{k}{N \cdot HLR}}\]

\[C_{out} = C^{*} + (C_{in} - C^{*})e^{-k / HLR}\]


Using the tank in series equation, here is the inverse relationship


Solve for required HLR (given target \(C_{out}\))

Rearrange the \(NKC^*\) formula:

\[\bigg( 1 + \frac{k}{N \cdot HLR} \bigg)^N = \frac{C_{in} - C^*}{C_{out} - C^{*}} \] \[HLR = \frac{k/N}{\big( \frac{C_{in} - C^*}{C_{out} - C^{*}} \big)^{1/N} -1}\] Then area for a specified flow \(Q\) is:

\[A = \frac{Q}{HLR}\] Domain constraints: you need \(C_{out} > C^{*}\) and \(C_{in} > C^{*}\). If \(C_{out} \le C^{*}\), the outflow is exceeding the assumed background concentration.


Solve for k (given observed \(C_{in}\),\(C_{out}\) and \(HLR\))

\[k = N \cdot HLR \bigg[ \bigg( \frac{C_{in} - C^{*}}{C_{out} - C^{*}} \bigg)^{1/N} -1 \bigg]\] This is the equivalent settling velocity (\(m \ yr^{-1}\)) implied by observed performance (at the chosen \(N\) and $C^{8}).

\[k = HLR \cdot \bigg( \frac{C_{in}-C^{*}}{C_{out}-C^{*}} -1 \bigg)\]

\[k = HLR \cdot ln\bigg( \frac{C_{in}-C^{*}}{C_{out}-C^{*}} \bigg)\]

This is often the method to estimate the upper bounds of \(k\) from observed annual data.


  1. DMSTA Model: Conceptually represented as

DMSTA State Variables and Parameters
Symbol Description Units
State Variables
M Water column P storage \(\mathrm{mg\ m^{-2}}\)
S Biomass P storage \(\mathrm{mg\ m^{-2}}\)
Z Mean water depth \(m\)
Driving Variables
L External P load \(\mathrm{mg\ m^{-2}\ yr^{-1}}\)
Q Outflow rate \(\mathrm{m\ yr^{-1}}\)
Parameter Variables
K1 Maximum uptake rate \(\mathrm{m^{3}\ mg^{-1}\ yr^{-1}}\)
K2 Recycle rate \(\mathrm{m^{2}\ mg^{-1}\ yr^{-1}}\)
K3 Burial rate \(\mathrm{yr^{-1}}\)
Steady-State Mass Balances 1 CSTR & Fz = 1
Compartment Equation
Storage: \(K1 \cdot C = K2 \cdot S + K3\)
Overall: \(L - Q \cdot C = K3 \cdot S\)
Solution for C: \(C = \frac{(K2 \cdot L + K3^{2})}{K3 \cdot K1 + Q \cdot K2}\)
Solution for S: \(S = \frac{(K1 \cdot C - K3)}{K2}\)

Whereas the steady-state DMSTA solution is equivalent to a simple \(KC^{*}\) model with \(C^{*} = C_0\) and \(K\) is interpreted as a net settling rate (\(m \ yr^{-1}\)). Therefore, for “equivalent design model” purposes \(C^{*} \equiv C_{0}\) and \(k \equiv K\). This means that the DMSTAs calibrated (\(K\) and \(C_{0}\)) can be subsituted directly into the \(NKC^{*}\) model.

Where \[C = \frac{(L + K \cdot C_{0})}{(K + Q)}\]

\[S = 1000 \cdot \frac{C-C_0}{C1 - C_0}\]

Where,

\[K = \frac{K2 \cdot K1}{K2}\]

Therefore, given (\(K\),\(C_0\),\(C_1\)), (\(K1\), \(K2\), \(K3\)) can be estimated:

\[K3 = K \cdot \frac{C_{1} - C_{0}}{1000}\]

\[K1 = \frac{K3}{C_0}\]

\[K_2 = \frac{K3 \cdot K1}{K}\]

Given (\(K1\), \(K2\), \(K3\)), (\(K\),\(C_0\),\(C_1\)) can be estimated:

\[C_{0} = \frac{K_3}{K_1}\]

\[K = \frac{K_1 \cdot K_3}{K_2}\]

\[ C_1 = \frac{1000 \cdot K_{2} + K_{3}}{K_{1}}\]

Helpful R Code

To see these equation in practice here are some simple R functions and examples.

# Forward NKC* (tanks-in-series) prediction
nkc_out <- function(Cin, k, HLR, N = 1, Cstar) {
  Cstar + (Cin - Cstar) / (1 + k / (N * HLR))^N
}

# Solve for HLR given target Cout
nkc_hlr <- function(Cin, Cout, k, N = 1, Cstar) {
  r <- (Cin - Cstar) / (Cout - Cstar)
  (k / N) / (r^(1/N) - 1)
}

# Solve for k given observed Cin/Cout and HLR
nkc_k <- function(Cin, Cout, HLR, N = 1, Cstar) {
  if (any(Cout <= Cstar))
    stop("Cout must be greater than Cstar")
  
  r <- (Cin - Cstar) / (Cout - Cstar)
  N * HLR * (r^(1/N) - 1)
}

# DMSTA parameter translations 
dmsta_KC_to_K123 <- function(K, C0, C1) {
  K3 <- K * (C1 - C0) / 1000
  K1 <- K3 / C0
  K2 <- K3 * K1 / K
  list(K1 = K1, K2 = K2, K3 = K3)
}

dmsta_K123_to_KC <- function(K1, K2, K3) {
  C0 <- K3 / K1
  K  <- (K1 * K3) / K2
  C1 <- (1000 * K2 + K3) / K1
  list(K = K, C0 = C0, C1 = C1)
}

# Estimate k from plug and flow
estimate_k_plug <- function(
  Cin, Cout,
  HLR,
  Cstar
) {
  if (any(Cout <= Cstar))
    stop("Cout must be greater than Cstar")

  HLR * log( (Cin - Cstar) / (Cout - Cstar) )
}
## Example Data
Cin <-  122
Cout <- 17
k  <-  66.5
HLR <- 9.5
Cstar <- 2

K <- 16.8
C0 <- 2
C1 <- 22

K1 <- 0.168
K2 <- 0.00336
K3 <- 0.336

# Estimated Cout given Cin, k, & HLR
nkc_out(Cin, k, HLR, N = 1, Cstar)
#> [1] 17

# Estimated HLR given Cin, Cout & k
nkc_hlr(Cin,Cout,k,N = 1,Cstar)
#> [1] 9.5

# vs more tanks
nkc_hlr(Cin,Cout,k,N = 3,Cstar)
#> [1] 22.16667

## Estimate k 
nkc_k(Cin,Cout,HLR,1,Cstar)
#> [1] 66.5

# estimate model parameters
dmsta_KC_to_K123(K,C0,C1)
#> $K1
#> [1] 0.168
#> 
#> $K2
#> [1] 0.00336
#> 
#> $K3
#> [1] 0.336

dmsta_K123_to_KC(K1,K2,K3)
#> $K
#> [1] 16.8
#> 
#> $C0
#> [1] 2
#> 
#> $C1
#> [1] 22