| Type: | Package |
| Title: | Dynamic Model for Stormwater Treatment Areas |
| Version: | 0.1.1 |
| Date: | 2026-07-23 |
| Maintainer: | Paul Julian <pjulian@evergladesfoundation.org> |
| URL: | https://github.com/SwampThingPaul/DMSTAr |
| BugReports: | https://github.com/SwampThingPaul/DMSTAr/issues |
| Description: | Performs treatment wetland modeling consistent with the original 'Excel' spreadsheet and VBA code version of DMSTA2 developed by Dr. Bill Walker http://wwwalker.net/dmsta/. |
| License: | MIT + file LICENSE |
| Encoding: | UTF-8 |
| LazyData: | true |
| Imports: | stats, utils |
| Depends: | R (≥ 4.1.0) |
| Suggests: | knitr, rmarkdown, kableExtra, testthat (≥ 3.0.0) |
| Config/testthat/edition: | 3 |
| RoxygenNote: | 7.3.3 |
| VignetteBuilder: | knitr |
| NeedsCompilation: | no |
| Packaged: | 2026-07-23 11:45:54 UTC; PaulJulian |
| Author: | Paul Julian |
| Repository: | CRAN |
| Date/Publication: | 2026-08-02 16:30:18 UTC |
Integrate DMSTA hydrology over one day using RKF45
Description
Integrates daily hydrology using adaptive Runge–Kutta–Fehlberg stepping. Intended for numerical diagnostics and sensitivity analysis rather than DMSTA parity runs.
Usage
.dmsta_ode_one_step(
Vo,
k,
Nsteps,
inputs,
params,
atol = 1e-08,
rtol = 1e-06,
max_steps = 200000L
)
Arguments
Vo |
Numeric scalar. Volume at start of day. |
k |
Numeric. number of steps. |
Nsteps |
Integer. Nominal steps per day. |
inputs |
List of daily hydrologic forcings. |
params |
List of hydrologic parameters. |
atol |
Absolute tolerance. |
rtol |
Relative tolerance. |
Value
Hydrology result object with adaptive diagnostics.
Resolve offline diversion Qin_Frac for a given date (DMSTA 2C2B)
Description
Implements the DMSTA 2C2B "offline" scheduling logic that temporarily
overrides Qin_Frac during an offline window of length offline_dur
starting at a fixed month/day each year, repeating every offline_freq years.
Usage
.dmsta_offline_qin_frac(
date,
base_frac,
offline_trigger = FALSE,
offline_start = as.Date("1965-03-15"),
offline_freq = 3L,
offline_dur = 45L,
offline_fracs = NULL,
frac_1 = NULL,
frac_2 = NULL,
frac_3 = NULL,
frac_4 = NULL,
frac_5 = NULL,
frac_6 = NULL
)
Details
The fraction used during the offline window is selected from a vector of candidate fractions (frac_1..frac_6 in DMSTA 2C2B).
Adaptive Runge–Kutta–Fehlberg (RKF45) ODE integrator integration diagnostics.
Description
Integrates an ordinary differential equation using adaptive Runge–Kutta–Fehlberg (4/5) stepping with local error control.
Usage
.dmsta_rkf45_integrate(
rhs,
y0,
t0,
t1,
h0 = NULL,
atol = 1e-08,
rtol = 1e-06,
hmin = 1e-10,
hmax = NULL,
max_steps = 200000L
)
Arguments
rhs |
Function computing dy/dt given (t, y). |
y0 |
Numeric. Initial state. |
t0 |
Numeric. Start time. |
t1 |
Numeric. End time. |
h0 |
Optional initial step size. |
atol |
Absolute tolerance. |
rtol |
Relative tolerance. |
hmin |
Minimum step size. |
hmax |
Maximum step size. |
max_steps |
Integer. Maximum number of integration steps. |
Details
This solver is intended for diagnostic and research use and is not required for DMSTA parity.
Trapezoidal integration helper
Description
Computes the trapezoidal integral of a time series. Used internally for diagnostic aggregation of sub-step outputs.
Usage
.dmsta_trapz_integrate(t, y)
Arguments
t |
Numeric vector. Time values. |
y |
Numeric vector. Values to integrate. |
Value
Numeric scalar. Trapezoidal integral.
Build multi-slot phosphorus kinetics parameter vectors
Description
Builds phosphorus kinetics for multiple model "slots" and assembles the
results into K-length vectors (where K = length(mods)).
Usage
build_P_kin_slots(
mods,
registry = NULL,
pparams,
Dpy = 365.25,
DutyCycle = NULL,
derive_PModel = TRUE,
default_PModel = 1L,
...
)
Arguments
mods |
Character vector of model identifiers (e.g., |
registry |
Optional model registry (named list of builder functions).
If |
pparams |
Named |
Dpy |
Numeric scalar. Time steps per year (default |
DutyCycle |
Numeric scalar in |
derive_PModel |
Logical; reserved for future use. |
default_PModel |
Integer scalar; reserved for future use. |
... |
Additional arguments passed to |
Details
Each slot is constructed by calling build_P_kinetics() with the same
parameter list pparams. This is useful when a simulation uses multiple
phosphorus modules that differ by model type but share the same raw parameter set.
Value
A named list containing:
- K1, K2, K3
Numeric vectors of length
K.- Chalf, Z_1, Z_2, Z_3, K2Coef
Numeric vectors of length
K.- Kslots
Integer scalar giving
K.- SeasonalFactor, Ytrans, Ysigma, Czero, PModel
Global scalars.
- mods
Character vector of model identifiers.
Examples
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
)
out <- build_P_kin_slots(
mods = c("STA", "PSTA", "RES"),
pparams = pparams,
Dpy = 365.25
)
out$K1
out$Z_1
Build phosphorus kinetics parameters for a registered model
Description
Builds a standardized parameter list for a selected phosphorus model type
(e.g., "STA", "PSTA", "RES") using a registered
model-builder function, then derives kinetic coefficients K1, K2,
K3 and the kinetic model indicator PModel via
compute_DMSTA_kvals().
Usage
build_P_kinetics(mod_type, Dpy = 365.25, DutyCycle = NULL, pparams, ...)
Arguments
mod_type |
Character scalar. Model identifier (e.g., |
Dpy |
Numeric scalar. Time steps per year (default |
DutyCycle |
Numeric scalar in |
pparams |
Named |
... |
Additional arguments passed through to the underlying model builder. |
Details
Model builders are obtained from the package phosphorus model registry.
The selected model builder should return a list containing at minimum
C1000, Cstar, and Ks (on the target time step).
Additional fields returned by the builder (e.g., Z1, Z2, Chalf)
are preserved. The function appends K1, K2, K3, and PModel.
The result is assigned class "P_kinetics" (prepended to any existing classes).
Value
A named list of standardized parameters with elements
K1, K2, K3, and PModel. The result has class
"P_kinetics".
Examples
# Example assumes DMSTAr ships with a registered "STA" builder.
pparams <- list(
C1000 = 1, Cstar = 0.2, Ks_per_yr = 0.5,
Z1 = 10, Z2 = 30, Z3 = 60,
K2Coef1 = 0.1, Chalf = 0.2, SeasonalFactor = 1,
DutyCycle = 0.95
)
pars <- build_P_kinetics("STA", Dpy = 365.25, pparams = pparams)
pars$K1
pars$PModel
Convert Cubic Feet per Second to Cubic Hectometers per Day
Description
Converts volumetric flow from cubic feet per second (cfs) to cubic hectometers per day (hm^3/day).
Usage
cfs_to_hm3d(x)
Arguments
x |
Numeric vector of flow values in cubic feet per second. |
Value
Numeric vector of flow values in cubic hectometers per day.
centimeters to meters
Description
centimeters to meters
Usage
cm_to_m(x)
Arguments
x |
numeric value |
Value
converted numeric value
Build DMSTA tank partitioning for a cell
Description
Internal helper that partitions a single cell area into a series of
conceptual "tanks" used by DMSTA. The number of tanks is derived from
ttankS; if ttankS is fractional, the final tank receives
the fractional area share and all preceding tanks receive equal shares.
Usage
dmsta_build_tanks(A_cell, ttankS, snap_last = TRUE)
Arguments
A_cell |
Numeric scalar > 0. Cell area (units consistent with the rest of the DMSTA implementation). |
ttankS |
Numeric scalar > 0. Effective number of tanks. May be fractional. |
snap_last |
Logical; if |
Details
The function also returns per-tank area fractions and cumulative fractions (useful for mapping depth/area relationships). Optionally, the last cumulative fraction can be "snapped" to exactly 1.0 for numerical stability.
Value
A named list with elements:
- Ntanks
Integer number of tanks.
- A_Tank
Numeric vector of tank areas, length
Ntanks.- F_Tank
Numeric vector of tank area fractions (
A_Tank / A_cell).- Fcum
Numeric vector of cumulative area fractions (
cumsum(F_Tank)).
See Also
dmsta_p_init_state() for initializing state vectors
compatible with the returned Ntanks.
Parameter Sets for Stormwater Treatment Area Dynamic Model
Description
Internal lookup table of calibration sets used by the Dynamic Model for Stormwater Treatment Areas (DMSTA). Each row corresponds to a vegetation / compartment set (e.g., emergent marsh, SAV, PSTA, reservoir) and provides model coefficients, depth/flow/concentration ranges, and variability terms.
Usage
dmsta_cals
Format
A data frame with 5 rows and 25 variables:
- Set
Character. Parameter set identifier (e.g.,
"EMG_3").- Descript
Character. Human-readable description of the set.
- C0
Integer. Concentration parameter C0 = Conc at 0 g/m2 P Storage (ug/L).
- C1
Integer. Concentration parameter C1 = Conc at 1 g/m2 P storage (ug/L).
- C2
Integer. Concentration parameter C2 = Conc at Half-Max Uptake.
- Ks
Numeric. Net Settling Rate at Steady State (m/yr).
- Z1
Integer. Saturated Uptake Depth (cm).
- Z2
Integer. Lower Penalty Depth (cm).
- Z3
Integer. Upper Penalty Depth (cm).
- K1
Numeric. First Order Removal Rate (m/yr).
- C0_NEWS_SF
Numeric. Periphyton for NEWS or Seasonal Adjustment for P Uptake (ug/L).
- C1_Peri
Numeric. Periphyton system - concentration parameter (ug/L).
- Ks_Peri
Numeric. Periphyton system - settling rate (1/yr).
- Zx_Peri
Numeric. Periphyton system - saturated uptake depth (cm).
- Sm
Numeric. Transition Storage Midpoint (mg/m2).
- Sb
Numeric. Transition Storage Bandwidth (mg/m2).
- MinDepth
Integer. Minimum depth observed/allowed for the set.
- MaxDepth
Integer. Maximum depth observed/allowed for the set.
- MinQW
Integer. Minimum flow (QW) observed/allowed for the set.
- MaxQW
Integer. Maximum flow (QW) observed/allowed for the set.
- MinConc
Numeric. Minimum concentration observed/allowed for the set.
- MaxConc
Numeric. Maximum concentration observed/allowed for the set.
- MinFreqZ_LT10cm
Numeric. Minimum frequency of Z < 10 cm.
- MaxFreqZ_LT10cm
Numeric. Maximum frequency of Z < 10 cm.
- K_CV
Numeric. Coefficient of variation for Ks (or related K term).
Source
dmsta xlsm file, sheet Calibrations
Examples
data(dmsta_cals)
dmsta_cals
Case-level network simulation utilities for DMSTAr
Description
A collection of helper functions that support DMSTA-style CASE networks in DMSTAr. These utilities provide:
Converts a DMSTA workbook "Networks" table into a normalized routing
data.frame suitable for downstream case-level network simulation.
Each non-blank destination in the network table becomes one route row.
Flexible extractor for pulling either the case-level output (level = "case")
or a specific cell-level output (level = "cell") from a list returned by
dmsta_flowP_case().
Resolves a stream identifier (e.g., "bypass", "release1",
"outflow") to the corresponding discharge (Q_*) and load
(L_*) columns in a DMSTAr case output table.
Given a routing table and a set of case names, computes a topological order (upstream to downstream) under the assumption that the case network is a DAG (i.e., contains no directed cycles).
Applies a simple discrete lag/lead to a numeric vector, padding with zeros.
Positive lag_days shifts values later in time (prepends zeros);
negative values shift earlier in time (appends zeros).
Executes a set of CASE simulations in topological (upstream-to-downstream) order and routes configured outflows/loads from upstream CASE outputs into downstream CASE inflow forcing, emulating DMSTA workbook "Networks" behavior.
Usage
build_routes_from_net_table(
net_table,
outlet_count = 5L,
case_col = "CaseName"
)
extract_df(res, level = c("case", "cell"), cell_index = 1L)
stream_map_cols(stream, outflow_def = c("treated", "total"))
topo_order_cases(routes, case_names)
lag_vec(x, lag_days)
run_network_of_cases(
cases,
net_table = NULL,
routes = NULL,
outlet_count = 5L,
verbose = TRUE,
check_route = FALSE,
Nsteps_case = NULL,
...
)
Arguments
net_table |
Optional DMSTA-style network table. If provided and |
outlet_count |
Integer number of outlet bins (used when parsing numeric destinations). |
case_col |
Character name of the column in |
res |
A list-like object, typically returned by |
level |
Character, one of |
cell_index |
Integer index of the cell to extract when |
stream |
Character stream name. Supported values:
|
outflow_def |
Character definition for |
routes |
Optional normalized routing table. If provided, |
case_names |
Character vector of all CASE IDs participating in the network. |
x |
Numeric vector. |
lag_days |
Integer number of days to shift. Positive shifts forward in time;
negative shifts backward. Zero returns |
cases |
Named list of CASE definitions. Names must be DMSTA CASE IDs. Each element must contain at least:
|
verbose |
Logical; if |
check_route |
Logical; if |
... |
Additional arguments passed through to |
Details
Parsing of DMSTA workbook "Networks" tables into routing definitions
Topological ordering of CASEs for downstream execution
Robust extraction of case- and cell-level output tables
Stream-to-column mapping for routing discharge and load
Vector lagging utilities for travel-time effects
A high-level wrapper to run downstream CASE networks
Together, these functions allow multiple DMSTA CASE simulations to be linked into a strictly downstream network, emulating the behavior of the DMSTA VBA workbook network framework while preserving DMSTAr's case- and cell-level APIs.
The network execution model follows standard DMSTA semantics:
CASEs are executed in upstream-to-downstream order.
Routed discharge and load from upstream CASE outputs are added to downstream CASE inflows.
Concentrations are recomputed safely as
C = L / Q.Routes may terminate at downstream CASEs or numeric outlet bins.
These helpers are primarily used by run_network_of_cases(),
but many are also useful independently for diagnostics, testing, and
custom orchestration.
Destinations are interpreted as:
A positive integer: route to an
"OUTLET"bin.Otherwise: route to another
"CASE"identified by string ID.Blank / whitespace: ignored (no route for that stream).
The returned routes include frac (default 1) and lag_days
(default 0) to support future extensions such as flow splitting or travel-time.
Current DMSTAr documentation shows dmsta_flowP_case() returning a list
with out$results$case and out$results$cells[[i]] (or similarly named
cell outputs).
This function supports that nested style and also checks common legacy fields
such as res$case / res$case_out.
DMSTAr case outputs include standardized Q_* and L_* fields for
bypass, releases, seepage discharge, and outflows.
This uses Kahn's algorithm over CASE-to-CASE edges (routes with to_type == "CASE").
If not all nodes can be ordered, the function stops, indicating a directed cycle.
For each CASE:
Base inflow discharge and concentration (
Qi,Ci) are sanitized.Routed inflow discharge/load from upstream CASEs are added.
A safe routed inflow concentration is computed (
Cin = Lin / Qin).-
dmsta_flowP_case()is executed for that CASE. Specified stream outputs are routed to downstream CASEs or outlet bins.
This wrapper relies on the standardized case output columns produced by
dmsta_flowP_case() (e.g., Q_out_bypass, Q_out_treated,
Q_out_total, and corresponding L_* fields).
It also expects that all CASEs share the same Date sequence (same length and alignment).
The outlet summary uses fw() to compute flow-weighted concentration (FWC) from total
load and discharge; ensure fw() is available in your namespace.
Value
A data.frame with one row per configured route, containing:
- from_case
Upstream CASE ID (character).
- stream
One of
"bypass","release1","release2","outflow","seepage".- to_type
"CASE"or"OUTLET".- to_id
Downstream CASE ID (character) or outlet index (integer stored as character).
- frac
Routing fraction (numeric), currently always 1.
- lag_days
Integer travel-time lag in days, currently always 0.
A data.frame containing the requested daily time series.
A named character vector of length 2 with names Q and L,
giving the discharge and load column names.
Character vector of CASE IDs ordered such that upstream cases appear before any downstream cases they route to.
Numeric vector of the same length as x.
A named list with elements:
- order
Character vector of CASE execution order.
- case_results
Named list of raw
dmsta_flowP_case()results for each CASE.- routed_in
Named list of routed inflow time series per CASE (
Q,L).- outlets
List of outlet-bin time series (
Q,L) of lengthoutlet_count.- outlet_summary
data.framesummarizing total outlet discharge/load and FWC.- ledger
data.frameof per-route totals (total_Q,total_L) for auditing.
Main user-facing functions
-
build_routes_from_net_table() -
extract_df() -
run_network_of_cases()
Internal helpers
-
extract_case_df -
stream_map_cols -
topo_order_cases -
lag_vec
See Also
dmsta_flowP_case()
Examples
net <- data.frame(
CaseName = c("STA1_DW", "STA1W"),
Bypass_to = c("", ""),
Release1_to = c("", ""),
Release2_to = c("", ""),
Outflow_to = c("STA1W", "1"),
Seepage_to = c("", ""),
stringsAsFactors = FALSE
)
build_routes_from_net_table(net, outlet_count = 1L)
out <- dmsta_flowP_case(series, cells, return_cell_series = TRUE)
df_case <- extract_df(out, level = "case")
df_cell1 <- extract_df(out, level = "cell", cell_index = 1)
lag_vec(1:5, 2) # 0 0 1 2 3
lag_vec(1:5, -2) # 3 4 5 0 0
# 1) Build cases
cases <- list(
STA1_DW = list(
series_base = STA1_DW_input, # data.frame with Date, Qi, Ci, Rain, Et, Zcontrol
cells = STA1_DW_cell # list of dmsta_make_cell(...) objects
),
STA1W = list(
series_base = STA1W_input,
cells = STA1W_cell
)
)
# 2) Build/parse routes
net <- data.frame(
CaseName = c("STA1_DW", "STA1W"),
Bypass_to = c("", ""),
Release1_to = c("", ""),
Release2_to = c("", ""),
Outflow_to = c("STA1W", "1"),
Seepage_to = c("", ""),
stringsAsFactors = FALSE
)
routes <- build_routes_from_net_table(net_table, outlet_count = 1L)
# 3) Run network
out <- run_network_of_cases(
cases = cases,
routes = routes,
Nsteps = 4L,
return_cell_series = TRUE
)
out$outlet_summary
head(out$ledger)
Integrate DMSTA hydrology over one day using Euler method
Description
Performs single-step Euler integration of DMSTA hydrology over one day. This integrator is intended for diagnostic and debugging purposes and does not provide DMSTA numerical parity.
Usage
dmsta_euler_hydro_day(V, inputs, params)
Arguments
V |
Numeric scalar. Volume at the start of the day (hm^3). |
inputs |
List of daily hydrologic forcings. |
params |
List of hydrologic parameters. |
Value
A list containing end-of-day volume and daily flux totals.
Run a networked DMSTA hydrology–phosphorus simulation
Description
Simulates coupled hydrology and phosphorus dynamics for a network of interconnected DMSTA cells over a daily time series. Each cell is simulated sequentially within each day, with treated outflows routed downstream according to network topology, splitter rules, and recycle indices.
Usage
dmsta_flowP_case(
series,
cells,
Nsteps = 4L,
N_plant = 30L,
Qmethod = c("RK4", "Euler", "RKF45", "custom"),
Pmethod = c("RK4", "Euler"),
integrator_fun = NULL,
interp_option = 2L,
max_iter = 1L,
conv_tol = 0.01,
return_cell_series = TRUE,
keep_Q17 = TRUE,
...
)
Arguments
series |
Data frame of daily watershed inputs. Must include
|
cells |
List of cell definitions created by
|
Nsteps |
Integer. Number of hydrology sub‑steps per day. |
N_plant |
Integer. Window length (days) for rolling mean depth used in reservoir penalty blending. |
Qmethod |
Character string specifying the hydrology integrator
( |
Pmethod |
Character string specifying the phosphorus integrator
( |
integrator_fun |
Optional custom hydrology integrator function. |
interp_option |
Control‑depth interpolation option (DMSTA semantics). |
max_iter |
Integer. Maximum number of network convergence iterations. |
conv_tol |
Numeric. Relative convergence tolerance on external phosphorus loads. |
return_cell_series |
Logical. If |
keep_Q17 |
Logical. If |
... |
Additional arguments passed to daily hydrology integration. |
Details
This function manages network‑level state, routing, lagged recycle
bookkeeping, and optional convergence iteration, while delegating
per‑cell daily physics to dmsta_flowP_day().
Network routing follows DMSTA conventions. Treated outflows are routed either to a downstream cell or out of system, while bypass, release, and seepage discharge streams leave the system immediately. Lagged seepage recycle is tracked explicitly as an internal transit reservoir.
For strict DMSTA parity, use Qmethod = "RK4" and
Pmethod = "RK4" with no operational overrides.
Value
An object of class "dmsta_network_result" containing:
- results
Case‑level and optional per‑cell daily output series.
- budgets
Water and phosphorus mass budgets at case and cell level.
- meta
Convergence diagnostics and model configuration metadata.
Examples
# Read data (internal)
data(series)
series <- series[1:370,]; # for example, limit input file
# Data formatting
series$Qi <- cfs_to_hm3d(series$Flow) # cfs to hm3/d
series$Rain <- in_to_m(series$Rainfall) # inches to meters per day
series$Et <- in_to_m(series$ET)
series$Zcontrol <- 0/100 # meters; setting to zero to see what happens
# If you have release series; otherwise set to 0
series$Qr0 <- 0 # constrained outflow (forced Q) if used
series$Qr1 <- 0 # release 1
series$Qr2 <- 0 # release 2
series$Ci <- series$Conc
# input parameters
# --- 1) Base hydrology params (shared structure) ---
hydro_base <- list(
A_cell = 2.19, # km2
# depths in cm (engine converts /100 to meters internally)
Zmin = 2, # cm
Zinit = 40, # cm
Zweir = 0, # cm
Q_zmin = 38, # cm
Zrelease = 0, # cm
Bypass_elev = 121.92, # cm (approx 1.2192 m)
# hydraulics
Q_a = 1.0,
Q_b = 4.0,
Width = 1.55, # km
Qomax = 0.0, # hm3/day; 0 disables max cap in this implementation
Qimax = 0.0, # hm3/day; 0 disables inflow cap
# seepage (rates in m/day per m head; elevations in cm)
Seepout_Rate = 0.0,
Seepout_Elev = 0.0, # cm
Seepin_Rate = 0.0,
Seepin_Elev = 0.0, # cm
ShutdownET = TRUE,
force_Q_out = FALSE,
DutyCycle = 0.95,
Cmax = 2000
)
# 2) Base P params (shared structure)
P_base <- list(
# STA module
C1000 = 22,
Cstar = 3,
Ks_per_yr = 16.8,
Z1 = 40,
Z2 = 100,
Z3 = 200,
Chalf = 300,
K2Coef1 = 0,
SeasonalFactor = 0, # keep 0 for base parity
# PSTA (NEWS transition)
Ytrans = 0,
Ysigma = 0,
Czero = 0,
C1000_2 = NULL,
ks_2 = 0,
zh_2 = 0,
# RES depth penalty
k_depth_penalty = 1,
# atmos + seepage water quality
C_rain = 10, # ppb (ug/L)
DryDepo = 20, # mg/m2-yr
seepage_c = 20, # ppb cap for seep outflow
seepin_conc = 0, # ppb
# initial P state
C_init_ppb = 30,
Y_init_mgm2 = 1000
)
# 3) Cell-specific params
params_cell1 <- modifyList(hydro_base, modifyList(P_base, list(
Qin_Frac = 0.22,
Seepout_Rate = 0.00789,
Ks_per_yr = 16.8,
Y_init_mgm2 = 3387.67297548954
)))
params_cell2 <- modifyList(hydro_base, modifyList(P_base, list(
Qin_Frac = 0,
Seepout_Rate = 0.00155,
Ks_per_yr = 52.5,
Y_init_mgm2 = 768.480041186681
)))
# 4) Build cells
cells <- list(
dmsta_make_cell(
label = "CELL1",
params = params_cell1,
ttankS = 3.0,
DownCell = 2L,
Qin_Frac = params_cell1$Qin_Frac,
RecycleIndex = 1L # self; can omit if your validator maps NA/0 -> self
),
dmsta_make_cell(
label = "CELL2",
params = params_cell2,
ttankS = 3.0,
DownCell = 0L,
Qin_Frac = params_cell2$Qin_Frac,
RecycleIndex = 2L # self
)
)
cells <- dmsta_validate_cells(cells)
# Run case
out <- dmsta_flowP_case(
series = series,
cells = cells,
Nsteps = 4L,
max_iter = 1L,
return_cell_series = TRUE,
keep_Q17 = TRUE
)
head(out$results$case)
head(out$results$cells[[1]])
head(out$results$cells[[2]])
Simulate one day of coupled hydrology and phosphorus dynamics
Description
Simulates a single day of DMSTA hydrology and phosphorus behavior for either a storage cell or a node. Node behavior (no storage, no phosphorus state) is handled algebraically, while storage cells are integrated using sub-day hydrology steps coupled to phosphorus integration.
Usage
dmsta_flowP_day(
V,
P_state,
tanks,
inputs,
params,
ppar,
constants,
Qmethod = c("RK4", "Euler", "RKF45", "custom"),
Pmethod = c("RK4", "Euler"),
Nsteps = 4L,
Z_plant = 0,
integrator_fun = NULL,
interp_option = 2L,
...
)
Arguments
V |
Numeric scalar. Volume at the start of the day (hm^3). |
P_state |
List with elements |
tanks |
List describing the tanks-in-series configuration. |
inputs |
List of daily forcing inputs (Qi, Ci, Rain, Et, Zcontrol, releases, recycle flows, etc.). |
params |
List of hydrologic parameters for the cell. |
ppar |
List of phosphorus kinetic parameters. |
constants |
List of physical and chemical constants used by phosphorus derivative calculations. |
Qmethod |
Character string specifying the hydrology integrator
( |
Pmethod |
Character string specifying the phosphorus integrator
( |
Nsteps |
Integer. Number of hydrology sub-steps per day. |
Z_plant |
Numeric scalar. Rolling-mean depth (m) used for reservoir penalty blending. |
integrator_fun |
Optional custom hydrology integrator function. |
interp_option |
Integer. Control-depth interpolation option (DMSTA semantics). |
... |
Additional arguments passed to hydrology integrators. |
Details
This function acts as the top-level daily orchestrator for coupled hydrology–phosphorus dynamics. Numerical integration details are delegated to lower-level step functions.
For strict DMSTA parity, use Qmethod = "RK4" and
Pmethod = "RK4" with no operational overrides. Optional operational
features (e.g., release pauses) are applied upstream in series-level
drivers and are not part of the DMSTA core formulation.
Value
A list with elements:
- results
List of daily end-of-day states and aggregated hydrology and phosphorus outputs.
- budgets
List containing daily water and phosphorus mass budgets.
- meta
Metadata describing the simulation day, methods used, and configuration flags.
Integrate hydrology–phosphorus dynamics over one day using sub-day steps
Description
Performs coupled hydrology and phosphorus integration over a single day
for a storage cell by looping over hydrologic sub-steps and integrating
phosphorus dynamics within each sub-step. Hydrology stepping is provided
by dmsta_flow_day_steps(), and phosphorus integration is delegated
to dmsta_P_step().
Usage
dmsta_flowP_day_steps(
V,
P_state,
tanks,
inputs,
params,
ppar,
constants,
Qmethod = c("RK4", "Euler", "RKF45", "custom"),
Pmethod = c("RK4", "Euler"),
Nsteps = 4L,
Z_plant = 0,
integrator_fun = NULL,
interp_option = 2L,
...
)
Arguments
V |
Numeric scalar. Volume at the start of the day (hm^3). |
P_state |
List with elements |
tanks |
List describing the tanks-in-series configuration
(e.g., output of |
inputs |
List of daily forcing inputs (Qi, Ci, Rain, Et, Zcontrol, releases, recycle flows, etc.), already prepared for the current day. |
params |
List of hydrologic parameters for the cell. |
ppar |
List of phosphorus kinetic parameters (output of
|
constants |
List of physical and chemical constants used by
|
Qmethod |
Character string specifying the hydrology integrator
(e.g., |
Pmethod |
Character string specifying the phosphorus integrator
(currently |
Nsteps |
Integer. Number of sub-day hydrology steps per day. |
Z_plant |
Numeric scalar. Rolling-mean depth (m) used for reservoir penalty blending in phosphorus kinetics. |
integrator_fun |
Optional custom hydrology integrator function,
used when |
interp_option |
Integer. Control-depth interpolation option (DMSTA semantics; default is mid-day). |
... |
Additional arguments passed to hydrology integrators. |
Details
This function assumes a storage cell with one or more tanks in series.
Node behavior (cells without storage or P state) is handled upstream
in dmsta_flowP_day() and should not be routed through this function.
Value
A list with elements:
- hyd
Normalized hydrology object returned by
dmsta_flow_day_steps().- P_state_end
List with updated
MandSvectors at end of day.- accum
Named list of daily accumulated flows, loads, and phosphorus mechanism totals (uptake, recycle, sedimentation, burial).
- storage
List of starting and ending total phosphorus storage for the day.
Run a coupled hydrology–phosphorus simulation over a time series
Description
Simulates DMSTA hydrology and phosphorus dynamics over a multi-day
time series for a single cell. This function manages time-series
iteration, initialization, rolling diagnostics (e.g., Z_plant),
and aggregation of daily results, while delegating daily integration
to dmsta_flowP_day().
Usage
dmsta_flowP_series(
series,
params,
pparams = NULL,
ttankS = 3,
Nsteps = 4L,
N_plant = 30L,
Qmethod = c("RK4", "Euler", "RKF45", "custom"),
Pmethod = c("RK4", "Euler"),
integrator_fun = NULL,
interp_option = 2L,
ppar = NULL,
constants = NULL,
tanks = NULL,
V_init = NULL,
init_P_state = NULL,
return_steps = FALSE,
...
)
Arguments
series |
Data frame containing daily forcing inputs. Must include
at least |
params |
List of hydrologic and phosphorus parameters. |
pparams |
Optional list of phosphorus parameters to merge into
|
ttankS |
Numeric. Number of tanks in series (may be fractional). |
Nsteps |
Integer. Number of hydrology sub-steps per day. |
N_plant |
Integer. Window length (days) used to compute rolling
mean depth for |
Qmethod |
Character string specifying the hydrology integrator. |
Pmethod |
Character string specifying the phosphorus integrator. |
integrator_fun |
Optional custom hydrology integrator function. |
interp_option |
Integer. Control-depth interpolation option. |
ppar |
Optional precomputed phosphorus kinetic parameter list. |
constants |
Optional list of constants used by phosphorus derivatives. |
tanks |
Optional pre-built tanks-in-series configuration. |
V_init |
Optional initial volume (hm^3). If |
init_P_state |
Optional initial phosphorus state. If |
return_steps |
Logical. If |
... |
Additional arguments passed to daily hydrology integration. |
Details
Series-level operational semantics (e.g., HydroIndex-style presence flags, optional release warm-up periods) are applied here and passed downstream as daily inputs.
Value
An object of class "dmsta_result" with elements:
- results
Data frame of daily hydrology and phosphorus outputs.
- budgets
List containing water and phosphorus mass budget data frames.
- meta
Metadata describing model configuration, initialization, and methods used.
Simulate one day of DMSTA hydrology
Description
High-level daily hydrology orchestrator applying DMSTA control-depth semantics, release gating, and integrator selection.
Usage
dmsta_flow_day(
V,
inputs,
params,
Qmethod = c("RK4", "Euler", "RKF45", "custom"),
Nsteps = 4L,
integrator_fun = NULL,
interp_option = 2L,
...
)
Arguments
V |
Numeric scalar. Volume at start of day. |
inputs |
List of daily inputs. |
params |
List of hydrologic parameters. |
Qmethod |
Character. Hydrology integrator. |
Nsteps |
Integer. Sub-steps per day. |
integrator_fun |
Optional custom integrator. |
interp_option |
Control-depth interpolation mode. |
... |
Additional arguments. |
Value
Daily hydrology results and diagnostics.
Return daily hydrology results with normalized sub-steps
Description
Wrapper around dmsta_flow_day() that guarantees
consistent access to hydrology sub-step records for
coupling with constituent models.
Usage
dmsta_flow_day_steps(
V,
inputs,
params,
Qmethod = c("RK4", "Euler", "RKF45", "custom"),
Nsteps = 4L,
integrator_fun = NULL,
interp_option = 2L,
...
)
Arguments
V |
Numeric scalar. Volume at start of day. |
inputs |
List of daily inputs. |
params |
List of hydrologic parameters. |
Qmethod |
Character. Hydrology integrator. |
Nsteps |
Integer. Sub-steps per day. |
integrator_fun |
Optional custom integrator. |
interp_option |
Control-depth interpolation mode. |
... |
Additional arguments. |
Value
Hydrology object including steps.
Run a DMSTA hydrology simulation over a time series
Description
Simulates hydrologic behavior for a single cell over a multi-day time series, managing initialization, rolling diagnostics, and result aggregation.
Usage
dmsta_flow_series(
series,
params,
V_init = NULL,
Qmethod = c("RK4", "Euler", "RKF45", "custom"),
Nsteps = 4L,
integrator_fun = NULL,
interp_option = 2L,
...
)
Arguments
series |
Data frame of daily inputs. |
params |
List of hydrologic parameters. |
V_init |
Optional initial volume. |
Qmethod |
Character. Hydrology integrator. |
Nsteps |
Integer. Sub-steps per day. |
integrator_fun |
Optional custom integrator. |
interp_option |
Control-depth interpolation mode. |
... |
Additional options. |
Value
Object of class "dmsta_result".
Examples
# Read data (internal)
data(series)
series <- series[1:370,]; # for example, limit input file
# Data formatting
series$Qi <- cfs_to_hm3d(series$Flow) # cfs to hm3/d
series$Rain <- in_to_m(series$Rainfall) # inches to meters per day
series$Et <- in_to_m(series$ET)
series$Zcontrol <- 0/100 # meters; setting to zero to see what happens
# If you have release series; otherwise set to 0
series$Qr0 <- 0 # constrained outflow (forced Q) if used
series$Qr1 <- 0 # release 1
series$Qr2 <- 0 # release 2
series$Ci <- series$Conc
# input parameters
params <- list(
A_cell = 2.19, # km2
Zmin = 2, # cm
Vmin = 0, # hm3
Q_a = 1.0, # qcoef_a; discharge coef
Q_b = 4.0, # qcoef_b; discharge exponent
Zweir = 0, # cm; qcoef_offset; depth offset for outflow computation
Q_zmin = 38, # cm; qcoef_zmin
Qomax = 0.0, # maximum discharge hm3/day
Qimax = 0, # maximum inflow (hm3/day)
Width = 1.55, # km
Bypass_elev = 121.92, # z_byp; mean depth at which bypass begins (m) from input (cm)
Seepout_Rate = 0.00789, # outflow seepage rate per unit head (m/day)/m from input cm/d/cm
Seepout_Elev = 0.0, # elevation controling outflow seepage rate (input cm)
Seepin_Rate = 0.0,
Seepin_Elev = 0.0,
ShutdownET = TRUE,
force_Q_out = FALSE,
wrap_interp = TRUE,
Zinit = 40, # cm; initial water column depth
Qin_Frac = 0.22, # inflow_frac; fraction of basin flows going into this cell
Zrelease = 0, # cm; z_release; minimum depth for releases
RecycleQ = 0,
IsaNode = NULL,
enable_P_release = FALSE,
K_release = 0,
IsaNode = NULL
)
V_init <- (cm_to_m(params$Zinit) * params$A_cell)
out <- dmsta_flow_series(V_init, series, params, Nsteps = 4)
out$results
Gate release flows based on storage and (optionally) hydraulic availability (DMSTA 2C2B)
Description
Applies DMSTA-style release gating logic to determine which release components are active based on current storage volume relative to a specified release depth. Optionally applies the DMSTA 2C2B hydraulic availability cap (Qorel) that scales releases when specified releases exceed hydraulically available discharge.
Usage
dmsta_gate_releases(
Vo,
A_cell,
Zrelease,
Qr_0,
Qr_1,
Qr_2,
dmsta_version = c("2E", "2C2B"),
Zrelsign = NA_real_,
Z = NULL,
Q_a = 0,
Q_b = 1,
Width = 1,
Zweir = 0,
Qomax = 0
)
Arguments
Vo |
Numeric scalar. Current storage volume (hm^3). |
A_cell |
Numeric scalar. Cell area (km^2). |
Zrelease |
Numeric scalar. Release depth threshold (m). Discretionary
releases are suppressed when |
Qr_0 |
Numeric scalar. Fixed outflow / gate-controlled component (hm^3/day). |
Qr_1 |
Numeric scalar. Discretionary release 1 (hm^3/day). |
Qr_2 |
Numeric scalar. Discretionary release 2 (hm^3/day). |
dmsta_version |
Character. DMSTA version semantics. |
Zrelsign |
Numeric scalar. Raw (pre-unit-conversion) z_release value used
by DMSTA 2C2B to enable hydraulic scaling when negative. Ignored unless
|
Z |
Numeric scalar. Current depth used to compute hydraulic availability (m).
If NULL, computed as |
Q_a, Q_b, Width, Zweir, Qomax |
Hydraulic parameters for Qorel computation, as in DMSTA. |
Value
A named list with elements:
- QrU_0
Effective fixed/gated component (hm^3/day).
- QrU_1
Effective discretionary release 1 (hm^3/day).
- QrU_2
Effective discretionary release 2 (hm^3/day).
- Qrelease
Total discretionary release (
QrU_1 + QrU_2, hm^3/day).- Sspec
Total specified outflow including all components (
QrU_0 + QrU_1 + QrU_2, hm^3/day).- Qorel
Hydraulically available discharge used for scaling (hm^3/day) or NA if not used.
Dispatch daily hydrology integration
Description
Selects and executes the requested hydrology integrator (RK4, Euler, RKF45, or custom) for one simulation day.
Usage
dmsta_hydro_day(
V,
inputs,
params,
method = c("RK4", "Euler", "RKF45", "custom"),
Nsteps = 4L,
integrator_fun = NULL,
...
)
Arguments
V |
Numeric scalar. Volume at start of day. |
inputs |
List of daily hydrologic forcings. |
params |
List of hydrologic parameters. |
method |
Character. Hydrology integrator method. |
Nsteps |
Integer. Sub-steps per day. |
integrator_fun |
Optional custom integrator. |
... |
Additional arguments passed to integrator. |
Value
Hydrology result object.
Initialize per-cell state for a network/case run
Description
Builds tank geometry and initializes hydrologic and phosphorus state for each cell using the cell parameters:
-
V[ic]is initialized asA_cell * (Zinit/100)(depth cm -> m), tank geometry is built using
dmsta_build_tanks(A_cell, ttankS),phosphorus states
MandSare initialized withdmsta_p_init_state()usingC_init_ppbandY_init_mgm2.
Usage
dmsta_init_case_state(cells)
Arguments
cells |
A validated list of cell definitions. Each cell must contain
|
Value
A list with elements:
- V
Numeric vector of initial volumes, length = number of cells.
- tanks
List of per-cell tank geometry objects.
- Pstate
List of per-cell phosphorus state objects (each with
MandSvectors).
Create a DMSTA network cell definition
Description
Constructs a cell definition used by dmsta_flowP_case() network/case
simulations. The function attaches:
per-cell tank specification
ttankS,routing fields (
DownCell,Qin_Frac,RecycleIndex),optional splitter specification (
SplitterFrac),per-cell kinetics parameters (
ppar) using built-in modules (STA,PSTA,RES),per-cell constants list used by the phosphorus derivative.
Usage
dmsta_make_cell(
label,
params,
ttankS,
DownCell = 0L,
Qin_Frac = 0,
RecycleIndex = NULL,
SplitterFrac = NULL
)
Arguments
label |
Character scalar. Cell label (use |
params |
Named list of per-cell parameters (must include at least
|
ttankS |
Numeric scalar > 0. Effective number of tanks (may be fractional). |
DownCell |
Integer index of the downstream cell (0 means terminal/out of system). |
Qin_Frac |
Numeric scalar. Fraction of watershed inflow assigned to this cell in network routing (case runner may override internally for “already allocated” inflows). |
RecycleIndex |
Integer index of the cell receiving recycled seepage flow/mass.
|
SplitterFrac |
Optional. Splitter routing fractions (only used if
|
Value
A named list representing a cell definition, suitable for inclusion
in the cells argument of dmsta_flowP_case().
See Also
dmsta_flowP_case() for running a multi-cell simulation.
Initialize DMSTA phosphorus state vectors for tanks
Description
Internal helper that initializes per-tank state vectors for the DMSTA phosphorus module given tank geometry and initial conditions.
Usage
dmsta_p_init_state(tanks, Z_init_m, C_init_ppb, Y_init_mgm2)
Arguments
tanks |
A list as returned by |
Z_init_m |
Numeric scalar. Initial water column depth (meters). |
C_init_ppb |
Numeric scalar. Initial concentration (ppb). |
Y_init_mgm2 |
Numeric scalar. Initial areal mass/loading (mg/m^2). |
Details
For each tank i, the function computes:
-
M[i] = C_init_ppb * A_Tank[i] * Z_init_m -
S[i] = Y_init_mgm2 * A_Tank[i]
where A_Tank[i] is the area of tank i. Units are assumed to be
consistent with the DMSTA implementation (e.g., depth in meters).
Value
A named list with elements:
- M
Numeric vector (length
tanks$Ntanks) of initialized water-column masses.- S
Numeric vector (length
tanks$Ntanks) of initialized sediment/areal stores.
See Also
dmsta_build_tanks() to generate tank geometry.
Integrate DMSTA hydrology over one day using RK4
Description
Integrates storage volume and hydrologic fluxes over a single day using fixed-step Runge–Kutta 4th-order integration. This function loops over sub-steps and aggregates daily totals.
Usage
dmsta_rk4_hydro_day(V, inputs, params, Nsteps = 4L)
Arguments
V |
Numeric scalar. Volume at the start of the day (hm^3). |
inputs |
List of daily hydrologic forcings. |
params |
List of hydrologic parameters. |
Nsteps |
Integer. Number of RK4 sub-steps per day. |
Details
This is the primary DMSTA-parity hydrology integrator.
Value
A list containing end-of-day volume, daily flux totals, and optional per-substep diagnostics.
Perform one RK4 hydrology sub-step (DMSTA DerivFlow kernel)
Description
Executes a single Runge–Kutta 4th-order (RK4) integration step for
DMSTA hydrology over a sub-day interval. This function evaluates
dmsta_DerivFlow() at multiple within-step locations and
returns the updated volume and diagnostic fluxes.
Usage
dmsta_rk4_hydro_step(V, step_index, Nsteps, Dt, inputs, params)
Arguments
V |
Numeric scalar. Volume at the start of the sub-step (hm^3). |
step_index |
Integer. Sub-step index within the day (1..Nsteps). |
Nsteps |
Integer. Total number of sub-steps per day. |
Dt |
Numeric scalar. Sub-step duration (days), typically 1/Nsteps. |
inputs |
List of hydrologic forcings and control variables (Qi, Rain, Et, Zcontrol, Zcontrol_prev, Zcontrol_next, releases, recycle). |
params |
List of hydrologic parameters passed to
|
Details
This is a low-level integrator used internally by daily hydrology drivers. It performs no looping over days.
Value
A list containing updated volume and instantaneous flux diagnostics for the sub-step.
Validate and normalize network cell definitions (internal)
Description
Checks that cells is a non-empty list of cell definitions and that each
cell contains the minimum required fields for a coupled hydrology + phosphorus
simulation.
Usage
dmsta_validate_cells(cells)
Arguments
cells |
A list of cell definitions (as created by |
Details
The validator enforces:
required structural fields:
params,ttankS,ppar,constantsrequired
pparfields used by the phosphorus derivativerequired
constantsfields used by the phosphorus derivativeDMSTA convention for recycle indices (0 means self)
at most one splitter cell labeled
"SPLITTER"(case-insensitive)
The function also fills defaults for some routing fields if missing.
Value
The normalized cells list (with defaults filled in). Invisibly
returns cells but typically used as cells <- dmsta_validate_cells(cells).
Robust neighbor-day control-depth lookup (internal)
Description
Convenience wrapper around neighbors_zcontrol() that guarantees
finite numeric control depths for today, previous day, and next day.
Usage
dmsta_zneighbors(i, z)
Arguments
i |
Integer index into |
z |
Numeric vector of daily control depths (m). |
Details
DMSTAr hydrology and coupled hydrology–phosphorus routines require
Zcontrol (today), Zcontrol_prev (yesterday), and
Zcontrol_next (tomorrow) to be finite numerics. This helper
enforces that contract by applying defensive fallbacks:
If
todayis not finite, it is set to 0.If
prev_dayornxtare not finite, they fall back totoday.
The returned list uses the same names as neighbors_zcontrol()
(today, prev_day, nxt) so callers can assign
consistently.
Value
Named list with numeric scalars:
- today
Control depth for day
i.- prev_day
Control depth for day
i-1(ortodayat the start).- nxt
Control depth for day
i+1(ortodayat the end).
Create a default DMSTAr parameter list
Description
Constructs a named list of model parameters used by DMSTAr. Most arguments
are scalar numeric values (or logical flags). Additional named parameters
may be supplied via ...; these will be appended to the returned list
and will override any existing defaults with the same name.
Usage
dmstar_default_params(
MT = FALSE,
Qin_Frac = 0,
A_cell = 0,
Width = 0,
Ntanks = 1,
Zrelease = 0,
Q_zmin = 0,
Zweir = 0,
Q_b = 0,
Q_a = 0,
Bypass_elev = 0,
Qomax = 0,
Qimax = 0,
Seepin_Rate = 0,
Seepin_Elev = 0,
seepin_conc = 0,
Seepout_Rate = 0,
Seepout_Elev = 0,
seepage_c = 20,
C_init_ppb = 0,
Y_init_mgm2 = 0,
Zinit = 40,
Cstar = 3,
C1000 = 22,
Chalf = 300,
Ks_per_yr = 0,
Z1 = 40,
Z2 = 100,
Z3 = 200,
K2Coef1 = 0,
SeasonalFactor = 0,
Ytrans = 0,
Ysigma = 0,
Czero = 0,
ks_2 = 0,
zh_2 = 0,
k_depth_penalty = 1,
C_rain = 10,
DryDepo = 20,
ShutdownET = TRUE,
force_Q_out = FALSE,
DutyCycle = 0.95,
Zmin = 2,
Cmax = 2000,
enable_P_release = FALSE,
K_release = 0,
dmsta_version = "2E",
offline_trigger = FALSE,
offline_start = NULL,
offline_freq = NULL,
offline_dur = NULL,
frac_1 = 0,
frac_2 = 0,
frac_3 = 0,
frac_4 = 0,
frac_5 = 0,
frac_6 = 0,
IsaNode = NULL,
...
)
Arguments
MT |
Logical. If |
Qin_Frac |
Numeric. Fraction of flow entering cell from basin. |
A_cell |
Numeric. Cell effective treatment are (units: km2). |
Width |
Numeric. mean width of flow path (units: km. |
Ntanks |
Integer-ish numeric. Number of tanks in series. |
Zrelease |
Numeric. Release elevation (units: cm). |
Q_zmin |
Numeric. Zc; no outflow below this water level , added to control depth specified in input series (units: cm). |
Zweir |
Numeric. Zw; fixed weir depth, use for reservoirs with outflow hydraulics controlled by outlet structure; = 0 for shallow systems when outflow is usually controlled by vegetation resistance (units: cm). |
Q_b |
Numeric. Rating curve parameter b; q / w = a (Z - Zw) ^b for Z >= Zc; 'typicallly ~ 3 to 4 for marsh control; 1.5 for weir control (reservoirs & other deep cells). |
Q_a |
Numeric. Rating curve parameter a; flow/width at water depth of 1 m (a); typically ~ 0.5 to 2. |
Bypass_elev |
Numeric. depth at which bypass begins ( 0 = no limit ). |
Qomax |
Numeric. Inflow capacity (triggers bypass) ( 0 = no limit ) (units: hm3/d). |
Qimax |
Numeric. Outflow capacity (0 = no limit) (units: hm3/d). |
Seepin_Rate |
Numeric. centimeters per day per centimeter of head, reflects transmissivity of soils (units: cm/d/cm). |
Seepin_Elev |
Numeric. 'drives inflow seepage rate, depth relative to mean ground surface elev (units: cm). |
seepin_conc |
Numeric. Concentration associated with seepage inflow (units: mg/m3). |
Seepout_Rate |
Numeric. centimeters per day per centimeter of head, reflects transmissivity of soils (units: cm/d/cm). |
Seepout_Elev |
Numeric. drives outflow seepage rate, relative to mean ground surface elev, can be < 0 (units: cm). |
seepage_c |
Numeric. Seepage concentration term (units: mg/m3). |
C_init_ppb |
Numeric. Initial concentration (units: mg/m3). |
Y_init_mgm2 |
Numeric. Initial biomass P storage (units: mg/m^2). |
Zinit |
Numeric. Initial water column depth relative to mean ground elevation (units: cm). |
Cstar |
Numeric. Reference concentration parameter also defined as water column conc at storage = 0 mg/m2 at steady-state (units: mg/m3). |
C1000 |
Numeric. water column conc at storage = 1000 mg/m2 at steady-state (units: mg/m3). |
Chalf |
Numeric. Half-saturation / half-response concentration parameter or water column concentration at 1/2 maximum uptake (units: mg/m3). |
Ks_per_yr |
Numeric. net settling rate in steady state in K/C* model (units: m/yr). |
Z1 |
Numeric. Uptake rate decreases below this depth; flat between Z1 and Z2; =0 ignored (units: cm). |
Z2 |
Numeric. uptake rate starts to decrease above this depth (30-day average); decreases linearly between Z2 & Z3; 0 = ignored; reflects damage to vegetation at high depths (units: cm) |
Z3 |
Numeric. upper end of depth penalty range; K = 1 m/yr at depths above this value, regardless of calibration; 0 = ignored (units: cm). |
K2Coef1 |
Numeric. Coefficient for secondary rate term (model-specific). |
SeasonalFactor |
Numeric. Seasonal factor multiplier (model-specific). |
Ytrans |
Numeric. Transform parameter for Y term (model-specific). |
Ysigma |
Numeric. Sigma/spread parameter for Y term (model-specific). |
Czero |
Numeric. Baseline concentration offset (model-specific). |
ks_2 |
Numeric. Secondary coefficient (model-specific). |
zh_2 |
Numeric. Secondary depth/elevation parameter (model-specific). |
k_depth_penalty |
Numeric. Depth penalty multiplier/parameter. |
C_rain |
Numeric. Rain concentration (model-specific units). |
DryDepo |
Numeric. Dry deposition loading term (model-specific units). |
ShutdownET |
Logical. If |
force_Q_out |
Logical. If |
DutyCycle |
Numeric. Duty cycle (0–1) applied to relevant process(es). |
Zmin |
Numeric. Minimum elevation (model-specific). |
Cmax |
Numeric. Maximum concentration cap (model-specific). |
enable_P_release |
Logical. If |
K_release |
Numeric. Currently a kg/day rate. For future implementation |
dmsta_version |
Text. Version of DMSTA to implement default is set to "2E", |
offline_trigger |
Logical. If |
offline_start |
Date. Implement STA rest period start date. Default is |
offline_freq |
Numeric. frequency of rest period. Default |
offline_dur |
Numeric. Duration of rest period. Default |
frac_1 |
Numeric. value used to update Qfrac during rest period implement |
frac_2 |
Numeric. value used to update Qfrac during rest period implement |
frac_3 |
Numeric. value used to update Qfrac during rest period implement |
frac_4 |
Numeric. value used to update Qfrac during rest period implement |
frac_5 |
Numeric. value used to update Qfrac during rest period implement |
frac_6 |
Numeric. value used to update Qfrac during rest period implement |
IsaNode |
Logical or |
... |
Additional named parameters to add to the returned list, or to
override existing defaults by name. All must be named (e.g., |
Details
-
Derived
IsaNode: IfIsaNodeisNULL, it is set to(A_cell <= 0); otherwise the providedIsaNodevalue is used. -
MTbehavior: WhenMT = TRUE, all numeric values in the returned parameter list are replaced with0, providing an empty list. Logical and non-numeric entries are left unchanged. -
Extra parameters via
...: All arguments passed in...must be named. New names trigger a warning and are added to the returned list. Names that already exist in the defaults override the default value.
Value
A named list of DMSTAr parameters.
Examples
# Get defaults
p <- dmstar_default_params()
# Override a few defaults
p <- dmstar_default_params(Zmin = 5, DutyCycle = 0.9, Ntanks = 3)
# Add a new parameter via ...
p <- dmstar_default_params(NewParam = 123)
# Derive IsaNode automatically when IsaNode is NULL
p <- dmstar_default_params(A_cell = 0) # IsaNode becomes TRUE
p <- dmstar_default_params(A_cell = 10) # IsaNode becomes FALSE
# MT = TRUE zeroes numeric parameters
p0 <- dmstar_default_params(MT = TRUE)
inches to meters
Description
inches to meters
Usage
in_to_m(x)
Arguments
x |
numeric value |
Value
converted numeric value
DMSTA network/case helpers (internal)
Description
Internal helpers used for DMSTA multi-cell / network (case) simulations, including per-cell initialization, validation, per-cell kinetics preparation, daily-mean volume calculations, and output column normalization.
Computes the DMSTA-style daily average volume as the time integral of the substep mean volume:
V_{day} = \sum_k \left(\frac{V_{o,k} + V_{k}}{2}\right)\Delta t_k
where each substep record provides Vo, V, and Dt.
Populates each cell definition with:
-
ppar: per-cell phosphorus kinetic parameters (typically frombuild_P_kin_slots(mods = c("STA","PSTA","RES"), ...)), -
constants: a per-cell constants list used by the phosphorus mass-balance derivative (e.g., concentration caps, deposition terms, seepage bookkeeping fractions).
Utility to standardize output column names and structure between per-cell and case-level result data frames. The function:
maps known aliases (e.g.,
RainVol_total->RainVol),builds a consistent set of inflow/outflow component columns,
optionally keeps seep recycle (Q17/L17/C17),
optionally carries through P-budget fields (e.g., columns prefixed with
P_),optionally preserves specified list-columns (e.g.,
mass_budget).
Usage
dmsta_daily_avg_volume(step_list)
dmsta_prepare_cells_modules(cells)
dmsta_case_components(
df,
keep_Q17 = TRUE,
keep_P = TRUE,
P_prefix = "^P_",
keep_list_cols = c("mass_budget"),
keep_extra = TRUE
)
Arguments
step_list |
A list of substep records. Each element should be a list
containing numeric fields |
cells |
A list of cell definitions. Each element must contain
|
df |
A data.frame containing DMSTA-style result columns. |
keep_Q17 |
Logical; if |
keep_P |
Logical; if |
P_prefix |
Regular expression used to match P-budget fields (default |
keep_list_cols |
Character vector of list-column names to carry through if present. |
keep_extra |
Logical; if |
Details
These functions are used by dmsta_flowP_case() and related
network orchestration logic.
The function returns the updated cells object.
This function uses the phosphorus model registry exposed as
P_MODEL_BUILDERS and assumes the STA, PSTA, and RES
builders are available in the registry.
Value
Numeric scalar. Daily-average volume integral (same volume units as
Vo and V).
The updated cells list with $ppar and $constants
added to each element.
A data.frame with standardized columns, plus optional carry-through fields.
See Also
dmsta_validate_cells() for strict validation of the
resulting cell definitions.
DMSTA hydrology engine (internal)
Description
Internal helpers implementing DMSTA hydrologic integration, including within-day control depth interpolation, derivative evaluation, and RK4 stepping.
Computes instantaneous hydrologic fluxes and the rate of change of
storage volume for a single within-day sub-step, following the
DMSTA DerivFlow logic implemented in VBA (Module1).
Usage
dmsta_DerivFlow(
V,
A_cell,
Qi,
Rain,
Et,
Zcontrol,
Zcontrol1 = Zcontrol,
Zcontrol2 = Zcontrol,
Step,
StepFrac,
Nsteps,
Ddt,
params,
Q_zmin = 0,
Zweir = 0,
Zmin = 0,
Vmin = NULL,
Seepin_Rate = 0,
Seepin_Elev = 0,
Seepout_Rate = 0,
Seepout_Elev = 0,
RecycleQ = 0,
Qrelease = 0,
Bypass_elev = 0,
Qimax = 0,
Qomax = 0,
Q_a = 0,
Q_b = 1,
Width = 1,
ShutdownET = TRUE,
has_outflow_constraint = FALSE,
has_depth_constraint = FALSE,
Qr_0 = 0
)
Arguments
V |
Numeric scalar. Current storage volume (hm^3). |
A_cell |
Numeric scalar. Cell area (km^2). |
Qi |
Numeric scalar. Inflow rate (hm^3/day). |
Rain |
Numeric scalar. Rainfall rate (m/day). |
Et |
Numeric scalar. Evapotranspiration rate (m/day). |
Zcontrol |
Numeric scalar. Control depth for the current day (m). |
Zcontrol1 |
Numeric scalar. Control depth for the previous day (m).
Defaults to |
Zcontrol2 |
Numeric scalar. Control depth for the next day (m).
Defaults to |
Step |
Integer. Sub-step index within the day (1-based). |
StepFrac |
Numeric scalar. Fractional position within the current sub-step (0–1). |
Nsteps |
Integer. Total number of sub-steps per day. |
Ddt |
Numeric scalar. Sub-step duration (days). |
params |
List of hydrologic parameters. Must include
|
Q_zmin |
Numeric scalar. Minimum depth offset required for outflow (m). |
Zweir |
Numeric scalar. Weir crest elevation (m). |
Zmin |
Numeric scalar. Minimum allowable depth (m). |
Vmin |
Numeric scalar. Minimum allowable volume (hm^3).
If |
Seepin_Rate |
Numeric scalar. Seep-in rate coefficient. |
Seepin_Elev |
Numeric scalar. External groundwater elevation for seep-in (m). |
Seepout_Rate |
Numeric scalar. Seep-out rate coefficient. |
Seepout_Elev |
Numeric scalar. External groundwater elevation for seep-out (m). |
RecycleQ |
Numeric scalar. Internal recycle inflow (hm^3/day). |
Qrelease |
Numeric scalar. Explicit release flow (hm^3/day). |
Bypass_elev |
Numeric scalar. Elevation above which full bypass occurs (m). |
Qimax |
Numeric scalar. Maximum allowable inflow before bypass (hm^3/day). |
Qomax |
Numeric scalar. Maximum allowable outflow (hm^3/day). Negative values indicate a minimum bypass constraint. |
Q_a |
Numeric scalar. Hydraulic discharge coefficient. |
Q_b |
Numeric scalar. Hydraulic discharge exponent. |
Width |
Numeric scalar. Effective weir width (km). |
ShutdownET |
Logical. If |
has_outflow_constraint |
Logical. Indicates presence of an explicit outflow constraint (DMSTA HydroIndex(2) > 0). |
has_depth_constraint |
Logical. Indicates presence of a depth constraint (DMSTA HydroIndex(1) > 0). |
Qr_0 |
Numeric scalar. Fixed outflow rate used when
|
Details
This function performs no time integration. It evaluates the
hydrologic balance at a specific fractional position within the day
(StepFrac) and returns the instantaneous derivative and
associated flux diagnostics. Numerical integration over a day is
handled by higher-level integrators (e.g., RK4, Euler).
Control-depth interpolation, outflow constraints, bypass logic, seepage, evapotranspiration shutdown, and release semantics follow DMSTA conventions and are governed by the supplied parameters and HydroIndex-style flags.
This function corresponds closely to the DMSTA VBA routine
DerivFlow and is intended for internal use by hydrologic
integrators. It should not be called directly by users.
For strict DMSTA parity, control-depth interpolation behavior is governed
by params$interp_option, with 2 (mid-day values) as the
default.
Value
A named list containing instantaneous derivatives and diagnostic quantities:
- Dvdt
Rate of change of storage volume (hm^3/day).
- Qot
Total outflow including releases (hm^3/day).
- Qo
Overflow (depth-dependent) outflow only (hm^3/day).
- Qnet
Net inflow excluding overflow (hm^3/day).
- Vnext
Trial volume at end of sub-step before trimming (hm^3).
- Z
Current depth (m).
- Zcont
Interpolated and constrained control depth (m).
- Vcontrol
Control volume corresponding to
Zcont(hm^3).- Etest
Effective evapotranspiration rate after any shutdown (m/day).
- AtmoS
Net atmospheric flux (rain – ET) as volume (hm^3/day).
- Seepin
Seep-in flux (hm^3/day).
- Seepout
Seep-out flux (hm^3/day).
- Bypass
Bypass flow (hm^3/day).
- Delta
Fractional position within the day used for interpolation.
DMSTA phosphorus engine (internal)
Description
Internal helpers implementing DMSTA coupled phosphorus mass balance within a tank series, including derivative evaluation and RK4 stepping.
Computes instantaneous derivatives dMdt and dSdt for one tank at
a single RK stage. This corresponds to DMSTA "DerivMass" logic and includes
uptake, recycle, sediment/burial, direct sedimentation, and advective transport
driven by tank outflow and seepage.
Advances tank states M and S forward by one substep of length
Dt (days) using a 4th-order Runge-Kutta integrator. Stage-specific
driver fields StepFrac, Ddt, Mo_fix, and So_fix
are injected before calling dmsta_deriv_mass().
Performs a single Euler (forward) integration step for DMSTA phosphorus
mass balances over a sub-day interval. This integrator evaluates
dmsta_DerivMass() once at mid-step (StepFrac = 0.5),
matching DMSTA conventions for within-day process evaluation.
Usage
dmsta_DerivMass(state, drivers, ppar, constants)
dmsta_rk4_P_step(
M,
S,
args_base,
Dt,
ppar,
constants,
clamp = TRUE,
Mmin = 0,
Smin = 0
)
dmsta_euler_P_step(
M,
S,
args_base,
Dt,
ppar,
constants,
clamp = TRUE,
Mmin = 0,
Smin = 0
)
dmsta_P_step(
method = c("RK4", "Euler"),
M,
S,
args_base,
Dt,
ppar,
constants,
...
)
Arguments
state |
Named list containing tank state:
|
drivers |
Named list of stage-specific drivers and geometry. Expected fields include:
|
ppar |
List of phosphorus model parameters (output of
|
constants |
List of physical and chemical constants used by
|
M |
Numeric scalar. Current water-column phosphorus mass (kg). |
S |
Numeric scalar. Current labile (biota/soil) phosphorus mass (kg). |
args_base |
List containing at least |
Dt |
Numeric scalar. Sub-step duration (days). |
clamp |
Logical. If |
Mmin |
Numeric scalar. Minimum allowable water-column P mass (kg). Default is 0. |
Smin |
Numeric scalar. Minimum allowable labile P mass (kg). Default is 0. |
method |
Character string specifying the P integrator to use.
One of |
... |
Additional arguments passed through to the selected integrator (e.g., clamping options). |
Details
These functions are used by dmsta_flowP_series().
Internal diagnostics are computed from stage volume vV:
-
C = M / vV(ppb-like concentration) -
Y = S / A_tank(areal store, mg/m^2 equivalent) -
z = vV / A_tank(m)
Kinetics differ by ppar$PModel:
-
PModel == 2L: transformed kinetics using slot-1 coefficients. otherwise: blended "transition" kinetics combining slots 1–2 and optionally slot-3 reservoir blending using
Z_plant.
The Euler P integrator is intended primarily for diagnostic,
debugging, and unit-testing purposes. For production runs and
DMSTA parity comparisons, dmsta_rk4_P_step() should be used.
This function does not implement any numerical integration logic
itself. It simply dispatches to the appropriate integrator based on
method. For strict DMSTA parity and production simulations,
method = "RK4" is recommended.
Dispatches phosphorus mass integration to the requested numerical integrator. This function provides a unified interface for multiple P integrators, mirroring the design of the DMSTAr hydrology engine.
Currently supported methods are:
-
"RK4": Fourth-order Runge–Kutta integrator (default, DMSTA parity). -
"Euler": Single-step Euler integrator (diagnostic/debugging use).
Value
A named list with elements:
- dMdt
Numeric scalar. Derivative of
M(kg/day).- dSdt
Numeric scalar. Derivative of
S(kg/day).- flux
List of per-area flux terms (rates) used internally:
P_uptake,P_recycle,P_sed,P_direct.- diag
List of diagnostics:
C,Y,z.
A named list containing updated states, mean derivatives, and averaged diagnostics over the step:
- M_new, S_new
Updated tank states (kg).
- dMdt_ts, dSdt_ts
Step-mean derivatives (kg/day).
- P_uptake_ts, P_recycle_ts, P_sed_ts, P_direct_ts
Step-mean flux terms (rates).
- C_ts, Y_ts, z_ts
Step-mean diagnostics.
A list with elements:
- M_new
Updated water-column phosphorus mass (kg).
- S_new
Updated labile phosphorus mass (kg).
- dMdt_ts
Instantaneous rate of change of
M(kg/day).- dSdt_ts
Instantaneous rate of change of
S(kg/day).- flux
Named list of phosphorus flux rates (per area), including uptake, recycle, sedimentation, burial, and release.
- diag
Named list of diagnostic variables (e.g., concentration, storage density, depth) evaluated at the derivative call.
A list returned by the selected P integrator, containing updated phosphorus state variables, instantaneous rates, fluxes, and diagnostic values.
DMSTA tank geometry and state initialization (internal)
Description
Internal helpers for constructing tank partitions and initializing DMSTA phosphorus state vectors.
Internal phosphorus kinetics helpers
Description
Internal helper that derives kinetic coefficients from C1000,
Cstar, and Ks. Two formulations are supported:
Usage
compute_DMSTA_kvals(C1000, Cstar, Ks)
Arguments
C1000 |
Numeric scalar. Concentration-like parameter (units consistent with the DMSTAr kinetic formulation). |
Cstar |
Numeric scalar. Kinetic threshold parameter. If negative, the transformed kinetics (PModel 2) are used. |
Ks |
Numeric scalar. Rate parameter on the time step of interest. |
Details
-
PModel 1 (standard): used when
C1000 > 0,Ks > 0, andCstar > 0. -
PModel 2 (transformed/special): triggered when
Cstar < 0, interpreted asG = -Cstar.
If inputs do not define a valid kinetic model, the function returns
K1 = K2 = K3 = 0 and PModel = NA_integer_.
Value
A named list with elements:
- K1
Numeric scalar.
- K2
Numeric scalar.
- K3
Numeric scalar.
- PModel
Integer scalar:
1L,2L, orNA_integer_.
meters to centimeters
Description
meters to centimeters
Usage
m_to_cm(x)
Arguments
x |
numeric value |
Value
converted numeric value
Build parameter lists for supported DMSTAr models
Description
These functions convert a raw parameter list (pparams) into a
standardized parameter list used by DMSTAr model implementations.
Usage
build_STA(Dpy, DutyCycle, pparams, ...)
build_PSTA(Dpy, DutyCycle, pparams, ...)
build_RES(Dpy, DutyCycle, pparams, ...)
Arguments
Dpy |
Number of time steps per year (e.g., 365 for daily). |
DutyCycle |
Fraction of time active in each step (0–1). |
pparams |
A named |
... |
Reserved for future extensions; ignored. |
Details
Builders perform unit conversions and time-step scaling:
rate parameters per year are converted to per-timestep using
Dpydepth parameters in centimeters are converted to meters
duty cycle is applied to rate parameters where applicable
Value
A named list of standardized model parameters.
Functions
-
build_STA(): Build parameters for the STA model. -
build_PSTA(): Build parameters for the PSTA model. -
build_RES(): Build parameters for the RES model.
Example DMSTA daily forcing time series
Description
A example time series of daily hydrology and phosphorus forcing inputs suitable for demonstrating DMSTAr case and cell simulations.
Usage
series
Format
A data frame with daily records and the following columns:
- Date
Date for the record. Stored as a character vector in
"YYYY-MM-DD"format in the raw object shown; you may convert toDatewithas.Date().- Flow
Daily inflow/flow rate (cfs).
- Conc
Daily inflow concentration (e.g., ug/L).
- Rainfall
Daily rainfall rate (e.g., in/day).
- ET
Daily evapotranspiration rate (e.g., in/day).
Details
This dataset is intended as a lightweight example for testing and
documentation. For use with functions that expect columns named
Qi, Ci, and Rain, you may need to rename columns
(e.g., Flow -> Qi, Conc -> Ci, Rainfall -> Rain).
Source
Source workbook: PROJECT_SFWMD_EC_01MAR2012_NET_EAA_STA1E.xls
sheet Series_Input
Examples
data(series)
# Convert Date column (if needed)
series$Date <- as.Date(series$Date)
# If your simulation expects Qi/Ci/Rain column names:
sim_series <- within(series, {
Qi <- Flow
Ci <- Conc
Rain <- Rainfall
})
sim_series$Flow <- NULL
sim_series$Conc <- NULL
sim_series$Rainfall <- NULL
head(sim_series)
Validate alignment of K-length phosphorus parameter vectors
Description
Internal helper that checks whether the K-length parameter vectors in
params are aligned with params$Kslots.
Usage
validate_P_paramsK(params)
Arguments
params |
A named |
Value
Invisibly returns TRUE if validation passes; otherwise throws
an error.