Phosphorus Kinetics in DMSTAr

Summary

This vignette documents the DMSTA-style phosphorus kinetics available in DMSTAr: how kinetic parameters are constructed, what the coefficients mean, how tank-in-series settings interact with kinetics, and how to customize kinetics in a reproducible R workflow.

This vignette focuses on within-case phosphorus dynamics. Network routing and hydrology setup are covered in vignette("Getting Started with DMSTAr").

Introduction

DMSTAr implements a modular version of DMSTA-style phosphorus dynamics. The goal is to keep the scientific intent clear (what the parameters mean) while supporting extensible implementation (how users provide parameters).

Background

Conceptual Structure

DMSTA represents phosphorus cycling with two storages:

• water-column concentration/state (C)

• labile biomass/sediment storage (S)

Three internal fluxes are represented with rate coefficients:

• removal (K1), recycle (K2), burial (K3)

Core internal flux forms (conceptually) are:

• \(Removal \propto K1 \cdot F_{c}(C) \cdot F_{z}(z) \cdot C \cdot S\)

• \(Recycle \propto K2 \cdot S^{2}\)

• \(Burial \propto K3 \cdot S\)

Where

• \(K1\), \(K2\), \(K3\) = removal, recycling and burial rate coefficients

• C = water-column P concentration

• S = liable P Storage

• \(F_{c}\) and \(F_{z}\) = efficiency modifiers

Concentration-dependent scale factor (\(F_{c}\))

\(F_{c}\) limits removal efficiency at low water column P concentration. DMSTA, at its core, assumes that biological and physicochemical uptake processes become less efficient as C approaches background levels. Therefore, \(F_{c}\) represents declining phosphorus uptake efficiency as water‑column concentrations approach low‑P conditions. Ultimately, this factor prevents the model from unrealistically removing phosphorus at full strength when concentrations are already very low.

In practice, \(F_{c}\) is governed by the parameter “C2” also called Chalf. It is defined as the concentration at which removal efficiency is ~50% of the maximum. Conceptually, at high concentrations, \(F_{c} \approx 1\), whereas at low concentrations, \(F_{c}\) is effectively 0.

Graphical representation of \(F_c\).

Depth Dependent scale factor (\(F_{z}\))

\(F_{z}\) adjusts removal efficiency based on water depth. Phosphorus removal processes in STAs are most effective over a preferred depth range. Very shallow or very deep water reduces contact with vegetation, periphyton, and sediments. \(F_z\) encodes that assumption without adding new state variables accounting for the reducted P removal efficiency below optimal water depths.

In DMSTA (and DMSTAr), \(F_z\) is controlled by three depth parameters, depth of maximum or saturated uptake (\(Z1\)), lower penalty depth (\(Z2\)) and upper penalty depth (\(Z3\)). Conceptually, depth near \(Z1\), \(F_z \approx 1\). However, when depth is \(<Z2\) or \(>Z3\) \(F_z <1\).

Graphical representation of \(F_z\).

Parameterization

In DMSTAr, phosphorus kinetics are created using a model-type builder (e.g. “STA”, “PSTA”, “RES”). Builders return a standardized set of raw parameters (including at least C1000, Cstar, Ks), and DMSTAr derives \(K1\), \(K2\) and \(K3\) plus a model indicator PModel.

DMSTAr, consistent with DMSTA, supports two internal formulations

• PModel 1: standard model

• PModel 2: special case, triggered when Cstar < 0

Phosphorus kinetics are calculated from a series of functions, some internal to the package but can be called using DMSTAr:::. Here is an example of a function to calculate the various K values and identifies which PModel is used in downstream analytical processes

DMSTAr:::compute_DMSTA_kvals(C1000 = 22, Cstar = 3, Ks = 16)
#> $K1
#> [1] 0.1013333
#> 
#> $K2
#> [1] 0.001925333
#> 
#> $K3
#> [1] 0.304
#> 
#> $PModel
#> [1] 1

When Cstar < 0, PModel 2 methodology is used to estimate K value

DMSTAr:::compute_DMSTA_kvals(C1000 = 22, Cstar = -3, Ks = 16)
#> $K1
#> [1] 5.333333
#> 
#> $K2
#> [1] -0.2346667
#> 
#> $K3
#> [1] 0.352
#> 
#> $PModel
#> [1] 2

Similar to DMSTA, DMSTAr calculates kinetic variables into three slots to perform P simulations depending on what models are applied. To build the necessary kinetics, build_P_kin_slots is used where parameters and model types are supplied. Currently the default models are STA, PSTA and RES but there is functionality to change the model types in the three slot configuration. Currently the function includes a Dpy argument which is a factor to convert per day to pre year. It is possible in the this argument will be moved into the function as a default.

pparams <- list(
  # shared / STA
  C1000 = 22, Cstar = 3, Ks_per_yr = 16,
  Z1 = 40, Z2 = 100, Z3 = 200,
  K2Coef1 = 0.1, Chalf = 50, SeasonalFactor = 1,
  # PSTA
  Ytrans = 1, Ysigma = 1, C1000_2 = 50, ks_2 = 20, zh_2 = 10,
  # RES
  k_depth_penalty = 0.5,
  DutyCycle = 0.95
)

kin_out <- build_P_kin_slots(
  mods = c("STA", "PSTA", "RES"),
  pparams = pparams,
  Dpy = 365.25
)
kin_out
#> $K1
#> [1] 2.635638e-04 2.548939e-03 8.236368e-06
#> 
#> $K2
#> [1] 5.007712e-06 1.248980e-04 1.564910e-07
#> 
#> $K3
#> [1] 0.0007906913 0.0025489391 0.0000247091
#> 
#> $Chalf
#> [1] 50  0 50
#> 
#> $Z_1
#> [1] 0.4 0.1 0.4
#> 
#> $Z_2
#> [1] 1 0 0
#> 
#> $Z_3
#> [1] 2 0 0
#> 
#> $K2Coef
#> [1] 0.0002737851 0.0000000000 0.0000000000
#> 
#> $Kslots
#> [1] 3
#> 
#> $SeasonalFactor
#> [1] 1
#> 
#> $Ytrans
#> [1] 1
#> 
#> $Ysigma
#> [1] 1
#> 
#> $CZero
#> [1] 0
#> 
#> $PModel
#> [1] 1
#> 
#> $mods
#> [1] "STA"  "PSTA" "RES"