Theoretical Foundation
Conservation Equations
There are two options available for computing the conservative transport of solutes (dissolved in water) or fully suspended sediments.
Since OpenWQ operates as a plug-in water quality module coupled to external host hydrological models, it does not have direct access to the spatial grid topology or neighbouring cells. Instead, the host model provides pairwise water flux exchanges between source and recipient cells. The transport equations are therefore formulated in terms of discrete mass fluxes between cell pairs rather than continuous PDEs.
Advection
The advective mass flux of a dissolved species from a source cell to a recipient cell is:
where \(Q_{s \rightarrow r}\) is the water flux from source to recipient \([MT^{-1}]\), \(W_s\) is the water mass in the source cell \([M]\), and \(m_s\) is the chemical mass in the source cell \([M]\). This is equivalent to transporting mass at the source cell concentration: \(F_{adv} = Q_{s \rightarrow r} \cdot C_s\).
Dispersion (Fickian approximation)
Since only pairwise source-recipient cell connections are available (without access to neighbouring cells for computing second-order spatial derivatives), the dispersion term is approximated using a Fickian dispersive flux between cell pairs:
where:
\(F_{disp}\) is the dispersive mass flux \([MT^{-1}]\)
\(C_s = m_s / W_s\) is the concentration in the source cell
\(C_r = m_r / W_r\) is the concentration in the recipient cell
\(W_s\) is the water mass in the source cell \([M]\)
\(D_{eff}\) is the effective dispersion rate \([T^{-1}]\)
The effective dispersion rate is pre-computed from user-provided parameters:
where \(D_{avg} = (D_x + D_y + D_z) / 3\) is the average of the directional dispersion coefficients \([L^2T^{-1}]\) (averaged because the orientation of each source-recipient connection is not known), and \(L\) is the user-provided characteristic length \([L]\) representing the distance between cell centers.
This formulation is bidirectional: the dispersive flux can move mass in either direction depending on the concentration gradient, smoothing concentration differences regardless of flow direction.
Total transport flux
The first option (OPENWQ_NATIVE_TD_ADVDISP) computes the total transport as the sum of advective and dispersive fluxes:
The second option (NATIVE_TD_ADV) accounts for transport through advection only:
Biogeochemistry: Reaction networks
These expressions are flexible and defined by the user.
They are powered by the comprehensive C++ Mathematical Expression Toolkit Library (ExprTk) developed by Arash Partow (1999-2020).
The implementation of ExprTk in OpenWQ is simple to use and provides an extremely efficient run-time mathematical expression parser and evaluation engine.
ExprTk supports numerous forms of functional, logical and vector processing semantics and is very easily extendible. In computational biogeochemistry, reaction kinetics usually take the general form of \(1^{st}\), \(2^{nd}\) or \(3^{rd}\)-multispecies reaction kinetics.
where \(c_a\) and \(c_b\) \([ML^{-3}]\) are the concentrations of chemical species \(a\) and \(b\), parameter/variable \(B\) represents weather/hydrological dependencies (such as soil moisture and temperature) and \(k\) is the reaction rate \([ML^{-3}T^{-1}]\). The reaction rate \(k\) can be provided as the reaction rate (usually using standard maximum at a reference temperature, often \(20^oC\)) or using expressions that can include relationships with the hydrological/weather dependency variables/parameters.
PHREEQC Geochemistry
As an alternative to the flexible BGC reaction networks, OpenWQ integrates the PHREEQC geochemical engine via PhreeqcRM. PHREEQC computes equilibrium speciation and kinetic reactions based on thermodynamic databases, solving the full set of mass-action equations for aqueous species, minerals, surfaces, and exchangers.
For a given solution, PHREEQC solves:
where \(K_i\) is the equilibrium constant for reaction \(i\), \(a_j\) is the activity of species \(j\), and \(\nu_{ij}\) is the stoichiometric coefficient.
This approach is more computationally expensive than the BGC-Flex approach, but provides thermodynamically consistent speciation calculations, which are essential for problems involving pH-dependent reactions, mineral interactions, or complex aqueous chemistry. See PHREEQC configuration for setup details.
Sorption Isotherms
OpenWQ includes two sorption isotherm models for simulating partitioning between dissolved and sorbed phases. Both use a kinetic approach: the equilibrium sorbed concentration is computed from the isotherm equation, and mass is transferred at a user-defined kinetic rate.
Freundlich isotherm – Models nonlinear sorption:
where \(q\) is the sorbed concentration \([ML^{-3}_{soil}]\), \(C\) is the dissolved concentration \([ML^{-3}]\), \(K_{fr}\) is the Freundlich coefficient, and \(N_{fr}\) is the Freundlich exponent. The equilibrium dissolved concentration is found via Newton-Raphson iteration on the mass conservation equation \(C_{eq} \cdot V + K_{fr} \cdot C_{eq}^{N_{fr}} \cdot \rho \cdot L = m_{total}\).
Langmuir isotherm – Models sorption with a finite number of binding sites:
where \(q_{max}\) is the maximum adsorption capacity \([MM^{-1}_{soil}]\) and \(K_L\) is the Langmuir equilibrium constant \([L^3M^{-1}]\). The equilibrium concentration is solved analytically via the quadratic formula.
Kinetic adsorption/desorption – For both isotherms, the mass flux from dissolved to sorbed phase is:
where \(K_{adsdes}\) is the kinetic adsorption/desorption rate \([T^{-1}]\), \(\rho\) is bulk density \([ML^{-3}]\), and \(L\) is layer thickness \([L]\). This flux is applied as a sink/source on the dissolved mass. See Sorption Isotherm configuration for setup details.
Note
The Freundlich/Langmuir SI module is intended for use with NATIVE_BGC_FLEX. When PHREEQC is the active biogeochemistry engine, sorption should be handled through PHREEQC’s SURFACE and EXCHANGE blocks, which provide mechanistically richer surface complexation and ion exchange models. See Process overlap with other OpenWQ modules.
Sediment Transport
OpenWQ includes two sediment transport erosion models based on the HYPE hydrological model framework. Both models compute erosion rates and transport sediments between compartments along a user-defined direction.
HYPE_HBVSED – HBV-based sediment erosion model
This model computes soil erosion as a function of precipitation, slope, soil properties, and land use. The erosion rate is computed as:
where:
\(\text{eroindex}\) – erosion index (spatially varying)
\(\text{eroindexpar}\) – scaling parameter for the erosion index
\(\text{lusepar}\) – land-use dependent soil erosion factor
\(\text{soilpar}\) – soil-type dependent erosion factor
\(\text{slope}\) – basin slope
\(\text{slopepar}\) – slope erosion exponent
\(\text{precip}\) – precipitation intensity
\(\text{precexppar}\) – precipitation erosion exponent
\(\text{erodmonth} = 1 + \text{monthpar}(m)\) – monthly erosion factor, where \(\text{monthpar}\) is a 12-value array (Jan–Dec) allowing seasonal variation
Erosion only occurs when precipitation exceeds zero and the inhibiting compartment (e.g., snow layer) has no water. Sediments are removed from the source compartment and added to the recipient compartment.
HYPE_MMF – Morgan-Morgan-Finney erosion model
This model uses the MMF approach to compute erosion based on kinetic energy from rainfall, soil cohesion, erodibility, and ground cover fractions. The erosion rate depends on:
where:
\(F\) – soil particle detachment by raindrop impact
\(H\) – resistance to erosion (soil cohesion and crop cover)
\(TC\) – transport capacity of overland flow
\(KE\) – total kinetic energy of rainfall (direct throughfall + leaf drainage)
\(SR\) – surface runoff depth
\(\text{cohesion}\) – soil cohesion (kPa)
\(\text{erodibility}\) – soil erodibility (g/J)
\(\text{sreroexp}\) – surface runoff erosion exponent
\(\text{cropcover}\) + \(\text{groundcover}\) – must sum to 1.0
\(\text{trans}_1\), \(\text{trans}_2\) – transport capacity factors
\(\text{slope}\) – basin slope
The pseudo day-of-year is dynamically computed from the simulation time and used to estimate kinetic energy from rainfall. Erosion only occurs when the inhibiting compartment has no water and crop/ground cover fractions are valid.
Both models support spatially-varying parameters via JSON or ASCII input formats, allowing cell-by-cell parameter overrides over the default values.
Numerical Solvers
OpenWQ provides two numerical solvers for integrating the water quality ODEs: