The community-centered freshwater biogeochemistry model unified RIVE v1.0: a unified version for water column

Abstract

Research on mechanisms of organic matter degradation, bacterial activities, phytoplankton dynamics, and other processes has led to the development of numerous sophisticated water quality models. The earliest model, dating back to 1925, was based on first-order kinetics for organic matter degradation. The community-centered freshwater biogeochemistry model RIVE was initially developed in 1994 and has subsequently been integrated into several software programs such as Seneque-Riverstrahler, pyNuts-Riverstrahler, ProSe/ProSe-PA, and Barman. After 30 years of research, the use of different programming languages including QBasic, Visual Basic, Fortran, ANSI C, and Python, as well as parallel evolution and the addition of new formalisms, raises questions about their comparability.

This paper presents a unified version of the RIVE model for the water column, including formalisms for bacterial communities (heterotrophic and nitrifying), primary producers, zooplankton, nutrients, inorganic carbon, and dissolved oxygen cycles. The unified RIVE model is open-source and implemented in Python 3 to create pyRIVE 1.0 and in ANSI C to create C-RIVE 0.32. The organic matter degradation module is validated by simulating batch experiments. The comparability of the pyRIVE 1.0 and C-RIVE 0.32 software is verified by modeling a river stretch case study. The case study considers the full biogeochemical cycles (microorganisms, nutrients, carbon, and oxygen) in the water column, as well as the effects of light and water temperature. The results show that the simulated concentrations of all state variables, including microorganisms and chemical species, are very similar for pyRIVE 1.0 and C-RIVE 0.32. This open-source project highly encourages contributions from the freshwater biogeochemistry community to further advance the project and achieve common objectives.

1 Introduction

Modeling the water quality of a freshwater system (river, lake, or reservoir) is critical to understanding and managing its functioning. The functioning of a freshwater system is the result of complex interrelated biogeochemical processes. The first water quality model developed by describes the degradation of organic matter (OM) in a river. The organic matter, measured globally by biochemical oxygen demand in 5 d (BOD5), is considered to be degraded according to first-order kinetics. Although dating back more than a century (the study was completed in 1915, but publication was delayed to 1925 due to World War I; ), this model is still widely used to represent the dynamics of organic matter in water quality modeling .

While the role of microorganisms in the degradation of organic matter has been acknowledged since the end of the 19th century, there is an important limitation of this type of representation. The microbiological nature of the organic matter degradation process and the bacterial population dynamics intrinsically involved are completely obscured. They are implicitly taken into account only through a biodegradability constant of organic matter and its dependence on temperature. Microbial biogeochemical work in the 1980s–1990s led to the elucidation of the detailed mechanisms of the organic matter degradation process and the associated heterotrophic bacterial activities . This new corpus of knowledge led to the development and the formulation of the biogeochemical model RIVE . It is capable of simulating the degradation of OM in freshwater systems and the associated oxygen consumption by bacterial activities, which is more realistic than the model of . In the RIVE model, the HSB model is used to represent the degradation of organic matter and heterotrophic bacterial activities. This model simulates the exoenzymatic hydrolysis of particulate and dissolved organic matter (split into biodegradable and refractory pools), including high-weight polymers, into small monomeric substrates. These substrates are subsequently assimilated by bacteria for their growth and respiration.

Apart from the degradation of organic matter, the AQUAPHY model is used to simulate the dynamics of phytoplankton in the RIVE model . The model explicitly simulates photosynthesis of phytoplankton, growth, mortality, and respiration processes. In addition to water temperature, the photosynthesis depends on the light intensity, while growth is controlled by nutrient availability and small organic metabolites. The small organic metabolites are formed either directly by photosynthesis or by catabolysis of reserve products. This conceptualization allows for the growth of phytoplankton during dark periods. In addition, the model also introduces a limiting factor of nutrients in the growth of phytoplankton and considers the cycling of nutrients during the life cycle of phytoplankton.

Since its initial development by , the RIVE model has co-existed within several software packages (Table ) developed for different aquatic compartments and supported by the PIREN-Seine program (https://www.piren-seine.fr/, last access: 20 December 2023). The RIVE model was firstly applied in river systems using the Riverstrahler drainage network approach . It was initially coded in QBasic and later on piloted by a GIS graphical interface, Seneque-Riverstrahler (Visual Basic; ). And it is now fully integrated within the pyNuts-Riverstrahler (https://gitlab.in2p3.fr/rive/pynuts/, last access: 20 December 2023) Python framework . It can model the biogeochemical functioning of hydrographic networks at scales ranging from local to continental. The RIVE model was also applied to lentic freshwater systems like regulated reservoirs (Barman software, – Table ) or to simulate hydro-biodynamic functioning of highly human-impacted river systems (ProSe software – ; ProSe-PA software – , https://gitlab.com/prose-pa/prose-pa last access: 20 December 2023; developed in ANSI C coupled with a self-developed lex and yacc parser; Table ). The RIVE model is also coupled with the Soil & Water Assessment Tool (SWAT) to simulate the water quality of the Vienne basin, France , and incorporated into the QUAL-NET model to simulate river eutrophication in the drainage network of the middle Loire River corridor, France. Moreover, the RIVE model is implemented into the VEMALA V3 model to simulate phosphorus and nitrogen loading in Finnish watersheds .

Based on the above implementations, different versions of the RIVE model code have successfully simulated a large variety of freshwater systems (lake or reservoirs, river systems) across the world. The parameter values were determined through laboratory experiments or calibrated with observation data . These applications (Table ) were carried out for different networks and scales as well as various degrees of anthropogenic impacts in a wide climatic gradient using either Riverstrahler (possibly with its Seneque or pyNuts environments) or ProSe/ProSe-PA, such as the Seine River (France) , the Danube River (Romania and Bulgaria) , the Red River (China and Vietnam) and its distributary the Day–Nhue River , the Lule and Kalix rivers (Sweden) , the Scheldt river (Belgium and Netherlands) , the Zenne River (Belgium) , the Mosel River (Germany) , the Somme River (France) , the Loire River (France) , the Lot River (France) , and the Orgeval watershed (France) . Moreover, the RIVE model has also been applied to stagnant systems (e.g., sand-pit lake – Lake Crétail, France, ; reservoirs – Marne, Aube, Seine, France, ).

After 30 years of research, the parallel evolutions of these codes, the numerical adaptations inherent in programming languages (QBasic, Visual Basic, Fortran, Python, and ANSI C), and the addition of new formalisms raise the question of their comparability. The identification of a unified version of the RIVE model is then necessary. A project aiming at unifying these RIVE implementations was undertaken. The unified version brings together all recent developments, especially the ones achieved with Python and ANSI C programming languages. This action will strengthen the collaboration of the research teams involved in the development of the model. This paper presents a unified version of RIVE for the water column (called unified RIVE v1.0) with a presentation of the formalisms for the biogeochemical cycles. That integrates the bacterial communities (heterotrophic and nitrifying), primary producers, zooplankton, and fate of detritic organic matter either particulate or dissolved and as biodegradable and refractory, as well as the associated nutrients and dissolved oxygen cycles. The most recent developments in modeling inorganic forms of carbon are also presented. The unified RIVE v1.0 included in pyRIVE 1.0 (tested with Python 3 versions up to 3.10 release) and C-RIVE 0.32 is open-source and therefore available to the scientific community. A numerical experiment is then introduced to evaluate the comparability of the pyRIVE 1.0 and C-RIVE 0.32 through a systematic comparison of simulations produced under controlled conditions. We thus establish a reference framework to evaluate different implementations (programming languages, performance – comparability) of the unified RIVE v1.0 formulation that continues to evolve in several water quality models.

2 Model description

The unified RIVE v1.0 model simulates the cycling of carbon, nutrients, and oxygen within a freshwater system (river, lake, reservoir). Biogeochemical cycles are simulated with a community-centered or agent-based model. That means that the freshwater system functioning is explicitly modeled, taking into consideration the activities of microorganisms such as phytoplankton, zooplankton, heterotrophic bacteria, and nitrifying bacteria. Additionally, it accounts for physical processes like oxygen re-aeration and dilution. This modeling approach is developed in relation to water temperature, macronutrients, and organic matter, (particulate, dissolved, and biodegradable fractions). The organic matter degradation, nitrifying bacteria dynamics, primary producer dynamics, zooplankton dynamics, nutrients, and inorganic carbon cycling are described subsequently. A high number of model parameters are used to characterize the microorganisms' properties and most of them have been determined through field or laboratory experiments under controlled conditions. This paper presents a focus on the conceptualization of the unified RIVE v1.0 model in the water column exclusively. While the RIVE model does have applications for sediment dynamics and its interaction with the water column , relevant community-centered efforts need to be made in future work, which is not the focus of this study.

2.1 Organic matter degradation

The mechanisms of organic matter degradation by the activity of heterotrophic bacteria are represented using the HSB model . It contains three variables: H represents high-weight polymers (large molecules), which form the majority of dissolved and particulate organic matter but must be exoenzymatically hydrolyzed to be accessible to heterotrophic bacteria; S represents small monomeric substrates (SMSs), directly accessible to microbial uptake; and B represents heterotrophic bacteria that absorb the substrates for their growth and respiration (Fig. ). In diagrams of this paper (for instance the HSB model, Fig. ), the state variables are shown as circles and represent either concentrations or stocks entering and leaving the (biogeochemical) processes. The biogeochemical processes are represented by squares.

The high-weight polymer (total organic carbon) is conceptually divided for each phase (dissolved – HD; particulate – HP) into three pools. Each pool is characterized by a specific biodegradability: (1) rapidly biodegradable in 5 d (HD${}_{\mathrm{1}}$ and HP${}_{\mathrm{1}}$), (2) slowly biodegradable in 45 d (HD${}_{\mathrm{2}}$ and HP${}_{\mathrm{2}}$), and (3) refractory (HD${}_{\mathrm{3}}$ and HP${}_{\mathrm{3}}$).

Figure 1

Flowchart of the HSB model. HD: dissolved high-weight polymer; HP: particulate high-weight polymer; SMS: small monomeric substrate; HB: heterotrophic bacteria; PHY: phytoplankton; ZOO: zooplankton; nitr. bact.: nitrifying bacteria; extr. excretion of phytoplankton; sink.: sinking. respi.: respiration; ${\mathit{ϵ}}_{\mathrm{hp}\mathrm{1},\mathrm{2},\mathrm{3}}$ and ${\mathit{ϵ}}_{\mathrm{hd}\mathrm{1},\mathrm{2},\mathrm{3}}$: proportion to convert dead biomass to HP and HD.

[Figure omitted. See PDF]

2.1.1 Heterotrophic bacteria dynamics

The dynamic of heterotrophic bacteria is explicitly simulated, including growth, mortality, and respiration. The growth of heterotrophic bacteria depends on water temperature and the availability of small monomeric substrate (SMS). The dependence is represented by the Monod equation . A maximal substrate uptake rate at 20 ${}^{\circ }$C ($b_{\max \mathrm{20},\mathrm{hb}}$) and a bacterial growth yield ($Y_{\mathrm{hb}}$) are used to calculate the growth rate of heterotrophic bacteria (${\mathit{\mu }}_{{\mathrm{hb}}_i}$) (Eq. ). The fraction of uptake not used for growth ($\mathrm{1}-Y_{\mathrm{hb}}$) is respired. $$\begin{array}{}\text{1}& {\displaystyle }{b}_{{\mathrm{hb}}_i}=b_{\max \mathrm{20},{\mathrm{hb}}_i}f(T{)}_{{\mathrm{hb}}_i}{\displaystyle \frac{\left[\mathrm{SMS}\right]}{\left[\mathrm{SMS}\right]+K_{\mathrm{sms},{\mathrm{hb}}_i}}}\text{2}& {\displaystyle }f(T{)}_{{\mathrm{hb}}_i}={\displaystyle \frac{e^{-\frac{{\left(T-T_{\mathrm{opt},{\mathrm{hb}}_i}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{{\mathrm{hb}}_i}^{\mathrm{2}}}}}{e^{-\frac{{\left(\mathrm{20}-T_{\mathrm{opt},{\mathrm{hb}}_i}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{{\mathrm{hb}}_i}^{\mathrm{2}}}}}}\text{3}& {\displaystyle }{\mathit{\mu }}_{{\mathrm{hb}}_i}=Y_{{\mathrm{hb}}_i}{b}_{{\mathrm{hb}}_i}\end{array}$$ Here, $b_{{\mathrm{hb}}_i}$ is the effective substrate uptake rate of the $i$th species of heterotrophic bacteria [h${}^{-\mathrm{1}}$], $b_{\max \mathrm{20},{\mathrm{hb}}_i}$ is the maximal substrate uptake rate of the $i$th species of heterotrophic bacteria at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], [SMS] is the small monomeric substrate concentration [mgC L${}^{-\mathrm{1}}$], $K_{\mathrm{sms},{\mathrm{hb}}_i}$ is the half-saturation constant for small monomeric substrate of the $i$th species of heterotrophic bacteria [mgC L${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{hb}}_i}$ is the water temperature weight of the $i$th species of heterotrophic bacteria at $T$ ${}^{\circ }$C [–], $T_{\mathrm{opt},{\mathrm{hb}}_i}$ is the optimal temperature of the $i$th species of heterotrophic bacteria for its growth [${}^{\circ }$C], ${\mathit{\sigma }}_{{\mathrm{hb}}_i}$ is the range of temperature for the $i$th species of heterotrophic bacteria [${}^{\circ }$C], $Y_{{\mathrm{hb}}_i}$ is the bacterial growth yield of the $i$th species of heterotrophic bacteria [–], and ${\mathit{\mu }}_{{\mathrm{hb}}_i}$ is the effective growth rate of the $i$th species of heterotrophic bacteria [h${}^{-\mathrm{1}}$]

A sinking velocity (${\mathrm{vs}}_{\mathrm{hb}}$) is associated with each particulate species to represent particulate sinking by gravity. The mortality of heterotrophic bacteria is simulated by first-order kinetics (Eq. ). The dead biomass of living species is converted into varying types of organic matter content, including both dissolved and particulate forms, based on specified proportions (${\mathit{ϵ}}_{\mathrm{hd}}$ and ${\mathit{ϵ}}_{\mathrm{hp}}$, Fig. ).

4$$\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{HB}}_i\right]}{\mathrm{d}t}}=({\mathit{\mu }}_{{\mathrm{hb}}_i}-k_{\mathrm{d}\mathrm{20},{\mathrm{hb}}_i}f(T{)}_{{\mathrm{hb}}_i}-k_{\mathrm{sink},{\mathrm{hb}}_i}\left)\right[{\mathrm{HB}}_i]\\ & k_{\mathrm{sink},{\mathrm{hb}}_i}={\displaystyle \frac{{\mathrm{vs}}_{{\mathrm{hb}}_i}}{\mathrm{depth}}}\end{array}$$ Here, ${\mathit{\mu }}_{{\mathrm{hb}}_i}$ is the effective growth rate of the $i$th species of heterotrophic bacteria [h${}^{-\mathrm{1}}$], $k_{\mathrm{d}\mathrm{20},{\mathrm{hb}}_i}$ is the mortality rate of the $i$th species of heterotrophic bacteria at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $k_{\mathrm{sink},{\mathrm{hb}}_i}$ is the sinking rate of the $i$th species of heterotrophic bacteria [h${}^{-\mathrm{1}}$], $\left[{\mathrm{HB}}_i\right]$ is the biomass concentration of the $i$th species of heterotrophic bacteria [mgC L${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{hb}}_i}$ is the water temperature weight at $T$ ${}^{\circ }$C defined by Eq. () [–], ${\mathrm{vs}}_{{\mathrm{hb}}_i}$ is the sinking velocity of the $i$th species of heterotrophic bacteria [m h${}^{-\mathrm{1}}$], and “depth” represents the water depth [m].

2.1.2 Hydrolysis of high-weight polymer

The particulate biodegradable high-weight polymer (HP${}_{\mathrm{1}}$ and HP${}_{\mathrm{2}}$) is firstly hydrolyzed to the dissolved biodegradable high-weight polymer (HD${}_{\mathrm{1}}$ and HD${}_{\mathrm{2}}$). The dissolved biodegradable high-weight polymer is then hydrolyzed exoenzymatically to small monomeric substrate (Fig. ). The hydrolysis of HP is represented by first-order kinetics (Eq. ), while a Michaelis–Menten function is used to express the exoenzymatic hydrolysis of HD (Eq. ).

5$$\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{HP}}_i\right]}{\mathrm{d}t}}& =-k_{{\mathrm{hp}}_i}\times \left[{\mathrm{HP}}_i\right]+\left(\sum _j{k}_{\mathrm{d}\mathrm{20},j}f\left(T{)}_j\right[\mathrm{LS}{]}_j\right)\\ & {\mathit{ϵ}}_{{\mathrm{hp}}_i}-k_{\mathrm{sink},{\mathrm{hp}}_i}\left[{\mathrm{HP}}_i\right]\end{array}$$ Here, $\left[{\mathrm{HP}}_i\right]$ is the concentration of particulate high-weight polymer $i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$ [mgC L${}^{-\mathrm{1}}$], $k_{{\mathrm{hp}}_i}$ is the hydrolysis rate of ${\mathrm{HP}}_i$ ($i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$) [h${}^{-\mathrm{1}}$], $f(T{)}_j$ is the water temperature weight of the $j$th living species at $T$ ${}^{\circ }$C defined as in Eq.  [–], $k_{\mathrm{d}\mathrm{20},j}$ is the mortality rate of the $j$th living species (such as phytoplankton, zooplankton, bacteria) at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $[\mathrm{LS}{]}_j$ is the concentration of the $j$th living species [mgC L${}^{-\mathrm{1}}$], ${\mathit{ϵ}}_{{\mathrm{hp}}_i}$ is the proportion to convert the dead biomass to ${\mathrm{HP}}_i$ ($i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$) [–], and $k_{\mathrm{sink},{\mathrm{hp}}_i}$ is the sinking rate for ${\mathrm{HP}}_i$ ($i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$) [h${}^{-\mathrm{1}}$]. 6$$\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{HD}}_i\right]}{\mathrm{d}t}}& =-\sum _k\left(e_{\max \mathrm{20},{\mathrm{hd}}_i,h{b}_k}f(T{)}_{{\mathrm{hb}}_k}{\displaystyle \frac{\left[{\mathrm{HD}}_i\right]}{\left[{\mathrm{HD}}_i\right]+K_{{\mathrm{hd}}_i,{\mathrm{hb}}_k}}}\right.\\ & \left.\left[{\mathrm{HB}}_k\right]\right)+k_{{\mathrm{hp}}_i}\times \left[{\mathrm{HP}}_i\right]\\ & +\left(\sum _j{k}_{\mathrm{d}\mathrm{20},j}f\left(T{)}_j\right[\mathrm{LS}{]}_j\right){\mathit{ϵ}}_{{\mathrm{hd}}_i}\end{array}$$ Here, $\left[{\mathrm{HD}}_i\right]$ is the concentration of dissolved high-weight polymer ($i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$) [mgC L${}^{-\mathrm{1}}$], $e_{\max \mathrm{20},{\mathrm{hd}}_i,{\mathrm{hb}}_k}$ is the maximum hydrolysis rate of ${\mathrm{HD}}_i$ at 20 ${}^{\circ }$C related to ${\mathrm{HB}}_k$ and $i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$ [h${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{hb}}_k}$ is the water temperature weight of the $k$th species of heterotrophic bacteria at $T$ ${}^{\circ }$C (Eq. ) [–], $K_{{\mathrm{hd}}_i,{\mathrm{hb}}_k}$ is the half-saturation constant for ${\mathrm{HD}}_i$ related to ${\mathrm{HB}}_k$ and $i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$ [mgC L${}^{-\mathrm{1}}$], $\left[{\mathrm{HB}}_k\right]$ is the concentration of the $k$th species of heterotrophic bacteria [mgC L${}^{-\mathrm{1}}$], $k_{\mathrm{d}\mathrm{20},j}$ is the mortality rate of the $j$th living species (such as phytoplankton, zooplankton, bacteria) at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $f(T{)}_j$ is the water temperature weight of the $j$th living species at $T$ ${}^{\circ }$C defined as in Eq. () [–], $[\mathrm{LS}{]}_j$ is the concentration of the $j$th living species [mgC L${}^{-\mathrm{1}}$], and ${\mathit{ϵ}}_{{\mathrm{hd}}_i}$ is the proportion to convert the dead biomass to ${\mathrm{HD}}_i$ and $i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$ [–].

Figure 2

Nitrifying bacteria dynamics. AOB: ammonia-oxidizing bacteria; NOB: nitrite-oxidizing bacteria; mort.: mortality; sink.: sinking.

[Figure omitted. See PDF]

2.2 Nitrifying bacteria dynamics

The unified RIVE v1.0 model includes the description of the nitrification microbial process, mediated by two types of nitrifying bacteria. They are respectively responsible for the production of nitrite (${\mathrm{NH}}_{\mathrm{4}}^{+}+\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}⟶{\mathrm{NO}}_{\mathrm{2}}^{-}+{\mathrm{2}\mathrm{H}}^{+}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$) and nitrate (${{\mathrm{NO}}_{\mathrm{2}}}^{-}+\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}⟶{\mathrm{NO}}_{\mathrm{3}}^{-}$). The nitrifying bacteria get energy by oxidizing NH${}_{\mathrm{4}}^{+}$ (ammonium) and NO${}_{\mathrm{2}}^{-}$ (nitrite) for their growth. These two bacteria are named AOB (ammonia-oxidizing bacteria) and NOB (nitrite-oxidizing bacteria) . The growth of nitrifying bacteria is limited by the availability of ammonium, nitrite, and oxygen, which is represented with Monod functions (Eq. ). The effect of water temperature is also taken into account. $$\begin{array}{}\text{7}& {\displaystyle }\begin{array}{rl}{\mathit{\mu }}_{\mathrm{aob}}& ={\mathit{\mu }}_{\max \mathrm{20},\mathrm{aob}}f(T{)}_{\mathrm{aob}}\left({\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]+K_{{\mathrm{nh}}_{\mathrm{4}},\mathrm{aob}}}}\right)\\ & \left({\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+K_{{\mathrm{o}}_{\mathrm{2}},\mathrm{aob}}}}\right)\end{array}\text{8}& {\displaystyle }\begin{array}{rl}{\mathit{\mu }}_{\mathrm{nob}}& ={\mathit{\mu }}_{\max \mathrm{20},\mathrm{nob}}f(T{)}_{\mathrm{nob}}\left({\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]}{\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]+K_{{\mathrm{no}}_{\mathrm{2}},\mathrm{nob}}}}\right)\\ & \left({\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+K_{{\mathrm{o}}_{\mathrm{2}},\mathrm{nob}}}}\right)\end{array}\end{array}$$ Here, ${\mathit{\mu }}_{\mathrm{aob}}$ and ${\mathit{\mu }}_{\mathrm{nob}}$ are the effective growth rates of AOB and NOB [h${}^{-\mathrm{1}}$], ${\mathit{\mu }}_{\max \mathrm{20},\mathrm{aob}}$ and ${\mathit{\mu }}_{\max \mathrm{20},\mathrm{nob}}$ are the maximal growth rates of AOB and NOB at 20 ${}^{\circ }$C, respectively [h${}^{-\mathrm{1}}$], $f(T{)}_{\mathrm{aob}}$ and $f(T{)}_{\mathrm{nob}}$ are the water temperature weight at $T$ ${}^{\circ }$C defined as in Eq. () [–], $K_{{\mathrm{nh}}_{\mathrm{4}},\mathrm{aob}}$ and $K_{{\mathrm{no}}_{\mathrm{2}},\mathrm{nob}}$ are the half-saturation constants for NH${}_{\mathrm{4}}^{+}$ (AOB) and for NO${}_{\mathrm{2}}^{-}$ (NOB) [mgN L${}^{-\mathrm{1}}$], and $K_{{\mathrm{o}}_{\mathrm{2}},\mathrm{aob}}$ and $K_{{\mathrm{o}}_{\mathrm{2}},\mathrm{nob}}$ are the half-saturation constants for oxygen (AOB and NOB) [mgO${}_{\mathrm{2}}$ L${}^{-\mathrm{1}}$].

The mortality and sinking of nitrifying bacteria are simulated the same way as for other living species. $$\begin{array}{}\text{9}& {\displaystyle \frac{\mathrm{d}\left[\mathrm{AOB}\right]}{\mathrm{d}t}}=({\mathit{\mu }}_{\mathrm{aob}}-k_{\mathrm{d}\mathrm{20},\mathrm{aob}}f(T{)}_{\mathrm{aob}}-k_{\mathrm{sink},\mathrm{aob}}\left)\right[\mathrm{AOB}]\text{10}& {\displaystyle \frac{\mathrm{d}\left[\mathrm{NOB}\right]}{\mathrm{d}t}}=({\mathit{\mu }}_{\mathrm{nob}}-k_{\mathrm{d}\mathrm{20},\mathrm{nob}}f(T{)}_{\mathrm{nob}}-k_{\mathrm{sink},\mathrm{nob}}\left)\right[\mathrm{NOB}]\end{array}$$ Here, ${\mathit{\mu }}_{\mathrm{aob}}$ and ${\mathit{\mu }}_{\mathrm{nob}}$ are the effective growth rates of AOB and NOB defined by Eqs. () and () [h${}^{-\mathrm{1}}$], $k_{\mathrm{d}\mathrm{20},\mathrm{aob}}$ and $k_{\mathrm{d}\mathrm{20},\mathrm{nob}}$ are the mortality rates of AOB and NOB at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $k_{\mathrm{sink},\mathrm{aob}}$ and $k_{\mathrm{sink},\mathrm{nob}}$ are the sinking rates of AOB and NOB [h${}^{-\mathrm{1}}$], $f(T{)}_{\mathrm{aob}}$ and $f(T{)}_{\mathrm{nob}}$ are the water temperature weights at $T$ ${}^{\circ }$C defined as in Eq. () [–], and [AOB] and [NOB] are the concentrations of AOB and NOB [mgC L${}^{-\mathrm{1}}$].

2.3 Primary producer dynamics

The behavior of primary producers is represented using the AQUAPHY model . Biomass of a phytoplankton species is composed of three different cellular constituents (Fig. ):

• a.

  the structural and functional macromolecules of the cell, F, mainly proteins, chlorophyll, and structural lipids (such as membranes);

• b.

  polysaccharides playing the role of reserve products, R;

• c.

  monomeric (amino acids) and oligomeric precursors for macromolecular synthesis, S.

Figure 3

Description of the AQUAPHY model. F: functional macromolecules of the cell; R: reserve products; S: monomeric (amino acids) and oligomeric precursors for macromolecular synthesis. Phytoplankton biomass equals the sum of the three cellular constituents (F, R, S). SMS: small monomeric substrate. photos.: photosynthesis; respi.: respiration; excr.: excretion; synth.: synthesis; catab.: catabolysis; sink.: sinking

[Figure omitted. See PDF]

At any time, the biomass of the $j$th phytoplankton species (mgC L${}^{-\mathrm{1}}$), $[\mathrm{PHY}{]}_j$, is equal to the sum of the three internal constituents (Eq. ), [F]${}_j$, [R]${}_j$, and [S]${}_j$ : 11$$[\mathrm{PHY}{]}_j=[\mathrm{F}{]}_j+[\mathrm{R}{]}_j+[\mathrm{S}{]}_j.$$

The most common way of measuring phytoplankton biomass is using the chlorophyll $a$ concentration ($\mathrm{\mu }$g chl $a$ L${}^{-\mathrm{1}}$). A carbon to chlorophyll $a$ ratio of 35 mgC $\mathrm{\mu }$g chl $a^{-\mathrm{1}}$ is therefore considered to convert experimental data into phytoplankton biomass. The initial proportions of different constituents (F, R, S) are fixed . They are only used to determine the initial concentrations of the three cellular constituents and their concentrations in incoming water fluxes for each phytoplankton species. According to , the structural and functional macromolecules of the cell ([F]${}_j$) account for about 85 % of the phytoplankton biomass ($[\mathrm{PHY}{]}_j$), while the reserve products ([R]${}_j$) account for about 10 % of the biomass. The remainder (5 %) of the biomass constitutes the small precursors for macromolecules synthesis ([S]${}_j$). These proportions of F, S, and R are updated at each time step for each phytoplankton species.

2.3.1 Photosynthesis

The photosynthesis process forms small precursors (S) by fixing carbon dioxide. Its rate is determined by the photosynthesis–irradiance relationship including three parameters (Eq. ) and the active irradiance ($I\left(z\right)$, $\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$).

12$$\begin{array}{rl}P(z{)}_{{\mathrm{phy}}_j}& =P_{\max \mathrm{20},{\mathrm{phy}}_j}f(T{)}_{{\mathrm{phy}}_j}\\ & \left(\mathrm{1}-e^{-\frac{{\mathit{\alpha }}_{{\mathrm{phy}}_j}I\left(z\right)}{P_{\max \mathrm{20},{\mathrm{phy}}_j}f\left(T\right)}}\right)e^{-\frac{{\mathit{\beta }}_{{\mathrm{phy}}_j}I\left(z\right)}{P_{\max \mathrm{20},{\mathrm{phy}}_j}f\left(T\right)}}\end{array}$$ Here, $P(z{)}_{{\mathrm{phy}}_j}$ is the photosynthesis rate of the $j$th phytoplankton species at water depth $z$ m [h${}^{-\mathrm{1}}$], $P_{\max \mathrm{20},{\mathrm{phy}}_j}$ is the maximal photosynthesis rate of the $j$th phytoplankton species at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{phy}}_j}$ is the water temperature weight of the $j$th phytoplankton species at $T$ ${}^{\circ }$C defined as in Eq. () [–], ${\mathit{\alpha }}_{{\mathrm{phy}}_j}$ is the photosynthetic efficiency of the $j$th phytoplankton species [h${}^{-\mathrm{1}}$ ($\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$], ${\mathit{\beta }}_{{\mathrm{phy}}_j}$ is the photoinhibition capacity of the $j$th phytoplankton species [h${}^{-\mathrm{1}}$ ($\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$], and $I\left(z\right)$ is photosynthetically active radiation (PAR) or active irradiance in the water column at depth $z$ m [$\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$] or [W m${}^{-\mathrm{2}}$].

The averaged photosynthesis rate of the $j$th phytoplankton species over water column is obtained by integrating $P\left(z\right)$: 13$$p_{{\mathrm{phy}}_j}={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{depth}}P(z{)}_{{\mathrm{phy}}_j}\mathrm{d}z}{\mathrm{depth}}},$$ where depth is the water height [m] and $p_{{\mathrm{phy}}_j}$ is the averaged photosynthesis rate over the water column [h${}^{-\mathrm{1}}$].

The active irradiance at water depth $z$ m ($I\left(z\right)$) follows the Beer–Lambert law (Eq. ). The decrease in active irradiance from the water surface to water bottom is represented by the light extinction coefficient ($\mathit{\eta }$). The extinction coefficient is composed of three parts: pure water (${\mathit{\eta }}_{\mathrm{base}}$), suspended solid (${\mathit{\eta }}_{\mathrm{ss}}$), and algal self-shading (${\mathit{\eta }}_{\mathrm{chl}a}$). 14$$\begin{array}{rl}& I\left(z\right)=I_{\mathrm{0}}{e}^{-\mathit{\eta }z}\\ & \mathit{\eta }={\mathit{\eta }}_{\mathrm{base}}+{\mathit{\eta }}_{\mathrm{chl}a}\left[\mathrm{chl}a\right]+{\mathit{\eta }}_{\mathrm{ss}}\left[\mathrm{SS}\right]\end{array}$$ Here, $I\left(z\right)$ is the active irradiance at water depth $z$ [$\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ or W m${}^{-\mathrm{2}}$], $I_{\mathrm{0}}$ is photosynthetically active radiation (PAR) or active irradiance at the water surface measured by the photosynthetic photon

Photons within the range of visible light between 400 and 700 nm.

flux density (PPFD) [$\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ or W m${}^{-\mathrm{2}}$], $\mathit{\eta }$ is the light extinction coefficient [m${}^{-\mathrm{1}}$], ${\mathit{\eta }}_{\mathrm{base}}$ is the light extinction coefficient related to pure water [m${}^{-\mathrm{1}}$], ${\mathit{\eta }}_{\mathrm{chl}a}$ is the linear algal self-shading light extinction coefficient [m${}^{-\mathrm{1}}$ ($\mathrm{\mu }$g chl $a$ L${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$], ${\mathit{\eta }}_{\mathrm{ss}}$ is the light extinction coefficient related to suspended solid [m${}^{-\mathrm{1}}$ (mg L${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$], [chl $a$] is the total chlorophyll $a$ concentration [$\mathrm{\mu }$g chl $a$ L${}^{-\mathrm{1}}$], and [SS] is the suspended solid concentration [mg L${}^{-\mathrm{1}}$].

2.3.2 Growth

The growth of phytoplankton involves the transformation of small precursors (S) into structural and functional macromolecules (F), which also requires the uptake of dissolved inorganic nutrients (nitrogen – N, phosphorus – P, silicon – Si) from the environment (Fig. ). The nutrients can potentially limit phytoplankton growth if their quantities are insufficient. The limitation of nutrients is represented using multiple Monod functions (Eq. ). The maximum growth rate (${\mathit{\mu }}_{\max ,{\mathrm{F}}_{\mathrm{j}}}$) is itself weighted by a limitation based on the availability of small precursors (S). The limitation by dissolved silica (DSi) is applied only for diatoms (DIA). $$\begin{array}{}\text{15}& {\displaystyle }{\mathit{\mu }}_{{\mathrm{F}}_j}={\mathit{\mu }}_{\max \mathrm{20},{\mathrm{F}}_{\mathrm{j}}}f(T{)}_{{\mathrm{phy}}_j}\left({\displaystyle \frac{\frac{\left[{\mathrm{S}}_j\right]}{\left[{\mathrm{F}}_j\right]}}{\frac{\left[{\mathrm{S}}_j\right]}{\left[{\mathrm{F}}_j\right]}+K_{\mathrm{S},{\mathrm{phy}}_j}}}\right)\mathrm{Nut}\mathit{\_}\lim \text{16}& {\displaystyle }\begin{array}{rl}& \mathrm{Nut}\mathit{\_}\lim =\min \left({\displaystyle \frac{\left[\mathrm{DIN}\right]}{\left[\mathrm{DIN}\right]+K_{\mathrm{N},{\mathrm{phy}}_j}}},{\displaystyle \frac{\left[\mathrm{DIP}\right]}{\left[\mathrm{DIP}\right]+K_{\mathrm{P},{\mathrm{phy}}_j}}},\right.\\ & \left.{\displaystyle \frac{\left[\mathrm{DSi}\right]}{\left[\mathrm{DSi}\right]+K_{\mathrm{Si},{\mathrm{phy}}_j}}}\right)\\ & \mathrm{or}\mathrm{Nut}\mathit{\_}\lim =\min \left({\displaystyle \frac{\left[\mathrm{DIN}\right]}{\left[\mathrm{DIN}\right]+K_{\mathrm{N},{\mathrm{phy}}_j}}},{\displaystyle \frac{\left[\mathrm{DIP}\right]}{\left[\mathrm{DIP}\right]+K_{\mathrm{P},{\mathrm{phy}}_j}}}\right)\end{array}\end{array}$$ Here, ${\mathit{\mu }}_{{\mathrm{F}}_j}$ is the effective growth rate of functional macromolecules for the $j$th phytoplankton species [h${}^{-\mathrm{1}}$], ${\mathit{\mu }}_{\max \mathrm{20},{\mathrm{F}}_j}$ is the maximal growth rate of functional macromolecules for the $j$th phytoplankton species at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{phy}}_j}$ is the water temperature weight of the $j$th phytoplankton species at $T$ ${}^{\circ }$C defined as in Eq. () [–], $\left[{\mathrm{S}}_j\right]$ and $\left[{\mathrm{F}}_j\right]$ are the concentrations of small precursors and functional macromolecules for the $j$th phytoplankton species [mgC L${}^{-\mathrm{1}}$], $K_{\mathrm{S},{\mathrm{phy}}_j}$ is the half-saturation constant for small precursors of the $j$th phytoplankton species [–], $\mathrm{Nut}\mathit{\_}\lim$ is the nutrient-limiting factor [–], [DIN], [DIP], and [DSi] are the concentrations of dissolved inorganic nitrogen ($\left[\mathrm{DIN}\right]=\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]$, mgN L${}^{-\mathrm{1}}$), dissolved inorganic phosphorus ($\left[\mathrm{DIP}\right]=\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]$, mgP L${}^{-\mathrm{1}}$), and dissolved silica (DSi, mgSi L${}^{-\mathrm{1}}$), $K_{\mathrm{N},{\mathrm{phy}}_j}$ and $K_{\mathrm{P},{\mathrm{phy}}_j}$ are the half-saturation constants for dissolved inorganic nitrogen and phosphorus of the $j$th phytoplankton species [mgN L${}^{-\mathrm{1}}$ and mgP L${}^{-\mathrm{1}}$], and $K_{\mathrm{Si},{\mathrm{phy}}_j}$ is the half-saturation constant for dissolved silica in the case of diatoms [mgSi L${}^{-\mathrm{1}}$].

2.3.3 Respiration

The respiration rate of phytoplankton ($r_{\mathrm{phy}}$) is divided into two components (Eq. ): one ($R_{\mathrm{m},\mathrm{phy}}$) ensuring the survival of the cell (maintenance process) and the other ($R_{\mathit{\mu },\mathrm{phy}}$) corresponding to the energetic cost of growth.

17$$r_{{\mathrm{phy}}_j}=R_{\mathrm{m}\mathrm{20},{\mathrm{phy}}_j}f(T{)}_{{\mathrm{phy}}_j}+{\mathit{\mu }}_{{\mathrm{F}}_j}{R}_{\mathit{\mu },{\mathrm{phy}}_j}$$ Here, $r_{{\mathrm{phy}}_j}$ is the respiration rate of the $j$th phytoplankton species [h${}^{-\mathrm{1}}$], $R_{\mathrm{m}\mathrm{20},{\mathrm{phy}}_j}$ is the maintenance respiration rate of the $j$th phytoplankton species at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{phy}}_j}$ is the water temperature weight of the $j$th phytoplankton species at $T$ ${}^{\circ }$C defined as in Eq. () [–], $R_{\mathit{\mu },{\mathrm{phy}}_j}$ is the respiration for the energetic cost of the $j$th phytoplankton species [–], and ${\mathit{\mu }}_{{\mathrm{F}}_j}$ is the effective growth rate of the $j$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$].

2.3.4 Excretion

Included later by , the phytoplankton excretion ($e_{\mathrm{phy}}$) includes two terms: a constant excretion rate ($E_{\mathrm{cst},\mathrm{phy}}$) and another that depends on the photosynthesis rate ($E_{\mathrm{phot},\mathrm{phy}}$). The product of excretion is the small monomeric substrate (SMS), assimilated directly by heterotrophic bacteria for their growth and respiration (Fig. ).

18$$e_{{\mathrm{phy}}_j}=E_{\mathrm{cst},{\mathrm{phy}}_j}+p_{{\mathrm{phy}}_j}{E}_{\mathrm{phot},{\mathrm{phy}}_j}$$ Here, $e_{{\mathrm{phy}}_j}$ is the excretion rate of the $j$th phytoplankton species [h${}^{-\mathrm{1}}$], $E_{\mathrm{cst},{\mathrm{phy}}_j}$ is the basic excretion rate of the $j$th phytoplankton species [h${}^{-\mathrm{1}}$], $E_{\mathrm{phot},{\mathrm{phy}}_j}$ is the excretion of the $j$th phytoplankton species related to photosynthesis [–], and $p_{{\mathrm{phy}}_j}$ is the photosynthesis rate of the $j$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$].

The variation of small monomeric substrate (SMS) can then be established (Eq. ). 19$${\displaystyle \frac{\mathrm{d}\left[\mathrm{SMS}\right]}{\mathrm{d}t}}=\mathrm{hydr}-\sum _i{b}_{{\mathrm{hb}}_i}\left[{\mathrm{HB}}_i\right]+\sum _j{e}_{{\mathrm{phy}}_j}\left[{\mathrm{F}}_j\right]$$ Here, hydr is the hydrolysis of the dissolved high-weight polymers HD${}_{\mathrm{1}}$ and HD${}_{\mathrm{2}}$ (Eq. ) [mgC L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], $b_{{\mathrm{hb}}_i}$ is the effective rate of substrate uptake by the $i$th heterotrophic bacteria species (Eq. ) [h${}^{-\mathrm{1}}$], $\left[{\mathrm{HB}}_i\right]$ is the biomass concentration of the $i$th species of heterotrophic bacteria [mgC L${}^{-\mathrm{1}}$], $e_{{\mathrm{phy}}_j}$ is the effective excretion rate of the $j$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], and $\left[{\mathrm{F}}_j\right]$ is the concentration of functional macromolecules for the $j$th phytoplankton species [mgC L${}^{-\mathrm{1}}$].

2.3.5 Synthesis and catabolysis of reserve products

The carbon fixed in the cell by photosynthesis forms small precursors (S) that can be transformed, either into functional macromolecules (F) or into reserve products (R). The synthesis of reserve products is limited by the $\frac{\left[\mathrm{S}\right]}{\left[\mathrm{F}\right]}$ ratio based on a Michaelis–Menten-like function (Eq. ).

20$$s_{\mathrm{R},{\mathrm{phy}}_j}=s_{\mathrm{R},\max \mathrm{20},{\mathrm{phy}}_j}f(T{)}_{{\mathrm{phy}}_j}{\displaystyle \frac{\frac{\left[{\mathrm{S}}_j\right]}{\left[{\mathrm{F}}_j\right]}}{\frac{\left[{\mathrm{S}}_j\right]}{\left[{\mathrm{F}}_j\right]}+K_{\mathrm{S},{\mathrm{phy}}_j}}}$$ Here, $s_{\mathrm{R},{\mathrm{phy}}_j}$ is the synthesis rate of reserve products of $j$th phytoplankton species [h${}^{-\mathrm{1}}$], $s_{\mathrm{R},\max \mathrm{20},{\mathrm{phy}}_j}$ is the maximal synthesis rate of reserve products of $j$th phytoplankton species at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{phy}}_j}$ is the water temperature weight of $j$th phytoplankton species at $T$ ${}^{\circ }$C defined as in Eq. () [–], $\left[{\mathrm{S}}_j\right]$ and $\left[{\mathrm{F}}_j\right]$ are the concentrations of small precursors and functional macromolecules for the $j$th phytoplankton species [mgC L${}^{-\mathrm{1}}$], and $K_{\mathrm{S},{\mathrm{phy}}_j}$ is the half-saturation constant for small precursors of the $j$th phytoplankton species [–].

Reserve products (R) are likely to be catabolized to produce small precursors (S). A first-order kinetic ($c_{\mathrm{R},\mathrm{phy}}$, h${}^{-\mathrm{1}}$) is used to represent catabolysis of a reserve product.

2.3.6 Extinction of phytoplankton

Three ways of phytoplankton extinction are implemented in the unified RIVE v1.0: lysis, sinking, and grazing by zooplankton (Sect. ). The phytoplankton lysis is represented by first-order kinetics using a mortality rate ($k_{\mathrm{d}\mathrm{20},\mathrm{phy}}$, h${}^{-\mathrm{1}}$). For ease of presentation, all three processes are assumed in an overall extinction rate $d_{\mathrm{phy}}$ (h${}^{-\mathrm{1}}$).

21$$\begin{array}{rl}{d}_{{\mathrm{phy}}_j}& =k_{\mathrm{d}\mathrm{20},{\mathrm{phy}}_j}f(T{)}_{{\mathrm{phy}}_j}+k_{\mathrm{sink},{\mathrm{phy}}_j}\\ & +\sum _i{\displaystyle \frac{b_{{\mathrm{zoo}}_i}\left[{\mathrm{ZOO}}_i\right]}{{\sum }_{k=\mathrm{1}}^{\mathrm{NS}}\left[{\mathrm{PHY}}_k\right]}}\end{array}$$ Here, $d_{{\mathrm{phy}}_j}$ is the extinction rate of the $j$th phytoplankton species [h${}^{-\mathrm{1}}$], $k_{\mathrm{d}\mathrm{20},{\mathrm{phy}}_j}$ is the mortality rate of the $j$th phytoplankton species at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{phy}}_j}$ is the water temperature weight of $j$th phytoplankton species at $T$ ${}^{\circ }$C defined as in Eq. () [–], $k_{\mathrm{sink},{\mathrm{phy}}_j}$ is the sinking rate of the $j$th phytoplankton species [h${}^{-\mathrm{1}}$], $b_{{\mathrm{zoo}}_i}$ is the grazing rate of the $i$th zooplankton species (Eq. , Sect. ) [h${}^{-\mathrm{1}}$], $\left[{\mathrm{ZOO}}_i\right]$ is the zooplankton concentration of the $i$th zooplankton species [mgC L${}^{-\mathrm{1}}$], and ${\sum }_{k=\mathrm{1}}^{\mathrm{NS}}\left[{\mathrm{PHY}}_k\right]$ is the total phytoplankton concentration (with NS the number of phytoplankton species grazed by zooplankton) [mgC L${}^{-\mathrm{1}}$].

2.3.7 Phytoplankton budgets

According to the processes related to phytoplankton (photosynthesis, growth, mortality, etc.), different budgets can be established for the $j$th phytoplankton species as follows. $$\begin{array}{}\text{22}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{S}}_j\right]}{\mathrm{d}t}}& =(p_{{\mathrm{phy}}_j}-r_{{\mathrm{phy}}_j}-{\mathit{\mu }}_{{\mathrm{F}}_j}-s_{\mathrm{R},{\mathrm{phy}}_j})\left[{\mathrm{F}}_j\right]\\ & +c_{\mathrm{R},{\mathrm{phy}}_j}\left[{\mathrm{R}}_j\right]-e_{{\mathrm{phy}}_j}\left[{\mathrm{F}}_j\right]-d_{{\mathrm{phy}}_j}\left[{\mathrm{S}}_j\right]\end{array}\text{23}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{R}}_j\right]}{\mathrm{d}t}}=s_{\mathrm{R},{\mathrm{phy}}_j}\left[{\mathrm{F}}_j\right]-c_{\mathrm{R},{\mathrm{phy}}_j}\left[{\mathrm{R}}_j\right]-d_{{\mathrm{phy}}_j}\left[{\mathrm{R}}_j\right]\text{24}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{F}}_j\right]}{\mathrm{d}t}}=({\mathit{\mu }}_{{\mathrm{F}}_j}-d_{{\mathrm{phy}}_j})\left[{\mathrm{F}}_j\right]\text{25}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{PHY}}_j\right]}{\mathrm{d}t}}& =(p_{{\mathrm{phy}}_j}-r_{{\mathrm{phy}}_j}-e_{{\mathrm{phy}}_j})\left[{\mathrm{F}}_j\right]\\ & -d_{{\mathrm{phy}}_j}\left[{\mathrm{PHY}}_j\right]\end{array}\end{array}$$ Here, $p_{{\mathrm{phy}}_j}$ is the photosynthesis rate of the $j$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $r_{{\mathrm{phy}}_j}$ is the respiration rate of the $j$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], ${\mathit{\mu }}_{{\mathrm{F}}_j}$ is the growth rate of the $j$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $s_{\mathrm{R},{\mathrm{phy}}_j}$ is the synthesis rate of reserve products of the $j$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $c_{\mathrm{R},{\mathrm{phy}}_j}$ is the catabolysis rate of reserve products of the $j$th phytoplankton species [h${}^{-\mathrm{1}}$], $e_{{\mathrm{phy}}_j}$ is the excretion rate of the $j$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $d_{{\mathrm{phy}}_j}$ is the extinction rate of the $j$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $\left[{\mathrm{S}}_j\right]$, $\left[{\mathrm{F}}_j\right]$, and $\left[{\mathrm{R}}_j\right]$ are the concentrations of $\left[{\mathrm{S}}_j\right]$, $\left[{\mathrm{F}}_j\right]$, and $\left[{\mathrm{R}}_j\right]$ for the $j$th phytoplankton species [mgC L${}^{-\mathrm{1}}$], and $\left[{\mathrm{PHY}}_j\right]$ is the biomass concentration of the $j$th phytoplankton species [mgC L${}^{-\mathrm{1}}$].

2.4 Zooplankton dynamics

The zooplankton dynamics include grazing on phytoplankton as well as growth, respiration, mortality, and sinking (Fig. ).

Figure 4

Dynamics of zooplankton. PHY: phytoplankton species; ZOO: zooplankton species; respi.: respiration; sink.: sinking.

[Figure omitted. See PDF]

2.4.1 Grazing and growth

Grazing on phytoplankton by zooplankton and growth of zooplankton are expressed based on a maximal grazing rate at 20 ${}^{\circ }$C ($b_{\max \mathrm{20},\mathrm{zoo}}$) limited by the phytoplankton biomass based on a Monod function (Eq. ). The grazing of zooplankton takes place only when the total phytoplankton biomass exceeds a certain threshold ($\left[{\mathrm{PHY}}_{\mathrm{0}}\right]$). No specific preference for grazing on particular phytoplankton species is considered among zooplankton species. Instead, the phytoplankton biomass grazed by the $i$th species of zooplankton is divided proportionally among each species of phytoplankton (Eq. ). The growth rate of zooplankton is considered proportional to the grazing rate using a growth yield factor (Eq. ). $$\begin{array}{}\text{26}& {\displaystyle }\begin{array}{rl}{b}_{{\mathrm{zoo}}_i}& =b_{\max \mathrm{20},{\mathrm{zoo}}_i}f(T{)}_{{\mathrm{zoo}}_i}\\ & {\displaystyle \frac{\left({\sum }_j^{\mathrm{NS}}\left[{\mathrm{PHY}}_j\right]-[{\mathrm{PHY}}_{\mathrm{0}}{]}_{{\mathrm{zoo}}_i}\right)}{\left({\sum }_j^{\mathrm{NS}}\left[{\mathrm{PHY}}_j\right]-[{\mathrm{PHY}}_{\mathrm{0}}{]}_{{\mathrm{zoo}}_i}\right)+K_{\mathrm{phy},{\mathrm{zoo}}_i}}}\end{array}\text{27}& {\displaystyle }{\mathit{\mu }}_{{\mathrm{zoo}}_i}=Y_{{\mathrm{zoo}}_i}{b}_{{\mathrm{zoo}}_i}\end{array}$$ Here, $b_{{\mathrm{zoo}}_i}$ is the effective grazing rate of the $i$th zooplankton species [h${}^{-\mathrm{1}}$], $b_{\max \mathrm{20},{\mathrm{zoo}}_i}$ is the maximal grazing rate of the $i$th zooplankton species at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{zoo}}_i}$ is the water temperature weight of the $i$th zooplankton species at $T$ ${}^{\circ }$C defined as in Eq. () [–], ${\sum }_j^{\mathrm{NS}}\left[{\mathrm{PHY}}_j\right]$ is the total phytoplankton biomass with NS the number of phytoplankton species grazed by zooplankton [mgC L${}^{-\mathrm{1}}$], $[{\mathrm{PHY}}_{\mathrm{0}}{]}_{{\mathrm{zoo}}_i}$ is the phytoplankton biomass threshold above which grazing takes place for the $i$th zooplankton species [mgC L${}^{-\mathrm{1}}$], $K_{\mathrm{phy},{\mathrm{zoo}}_i}$ is the half-saturation constant for phytoplankton biomass of the $i$th zooplankton species [mgC L${}^{-\mathrm{1}}$], ${\mathit{\mu }}_{{\mathrm{zoo}}_i}$ is the growth rate of the $i$th zooplankton species [h${}^{-\mathrm{1}}$], and $Y_{{\mathrm{zoo}}_i}$ is the growth yield of the $i$th zooplankton species [–].

2.4.2 Respiration and mortality

Grazed phytoplankton not used for zooplankton growth is respired (Fig. ). The rate of respiration is then obtained by $(\mathrm{1}-Y_{\mathrm{zoo}})\times b_{\mathrm{zoo}}$. The mortality of zooplankton is simulated by first-order kinetics ($k_{\mathrm{d}\mathrm{20},\mathrm{zoo}}$). $$\begin{array}{}\text{28}& {\displaystyle }{r}_{{\mathrm{zoo}}_i}=(\mathrm{1}-Y_{{\mathrm{zoo}}_i})\times b_{{\mathrm{zoo}}_i}\text{29}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{ZOO}}_i\right]}{\mathrm{d}t}}& =\left(Y_{{\mathrm{zoo}}_i}{b}_{{\mathrm{zoo}}_i}-k_{\mathrm{d}\mathrm{20},{\mathrm{zoo}}_i}f(T{)}_{{\mathrm{zoo}}_i}\right.\\ & \left.-k_{\mathrm{sink},{\mathrm{zoo}}_i}\right)\left[{\mathrm{ZOO}}_i\right]\end{array}\end{array}$$ Here, $r_{{\mathrm{zoo}}_i}$ is the respiration rate of the $i$th zooplankton species [h${}^{-\mathrm{1}}$], $Y_{{\mathrm{zoo}}_i}$ is the growth yield of the $i$th zooplankton species [–], $b_{{\mathrm{zoo}}_i}$ is the effective grazing rate of the $i$th zooplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $k_{\mathrm{d}\mathrm{20},{\mathrm{zoo}}_i}$ is the mortality rate of the $i$th zooplankton species at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], $f(T{)}_{{\mathrm{zoo}}_i}$ is the water temperature weight of the $i$th zooplankton species at $T$ ${}^{\circ }$C defined as in Eq. () [–], $k_{\mathrm{sink},{\mathrm{zoo}}_i}$ is the sinking rate of the $i$th zooplankton species [h${}^{-\mathrm{1}}$], and $\left[{\mathrm{ZOO}}_i\right]$ is the biomass concentration of the $i$th zooplankton species [mgC L${}^{-\mathrm{1}}$].

2.5 Nutrient cycling

As shown above, several processes related to nutrients are taken into account: uptake by phytoplankton, mineralization, nitrification, and denitrification (Figs.  and ).

2.5.1 Uptake of nutrients (N, P, Si) by phytoplankton

The Redfield–Conley stoichiometry ($\mathrm{C}:\mathrm{N}:\mathrm{P}:\mathrm{Si}$ $=\mathrm{106}:\mathrm{16}:\mathrm{1}:\mathrm{42}$; – ) is used to determine the composition of carbon, nitrogen, and phosphorus in organic matter. Constant $\mathrm{C}/\mathrm{N}$, $\mathrm{C}/\mathrm{P}$, and $\mathrm{C}/\mathrm{Si}$ mass ratios are considered to calculate the uptake of nutrients associated with phytoplankton growth. $$\begin{array}{}\text{30}& {\displaystyle \frac{\mathrm{d}\left[\mathrm{uptN}\right]}{\mathrm{d}t}}=\sum _i{\displaystyle \frac{({\mathit{\mu }}_{{\mathrm{F}}_i}+e_{{\mathrm{phy}}_i})\left[{\mathrm{F}}_i\right]}{\mathrm{C}/\mathrm{N}}}\text{31}& {\displaystyle }{\mathrm{uptNH}}_{\mathrm{4}}^{+}=\min \left(\right[{\mathrm{NH}}_{\mathrm{4}}^{+}],\mathrm{uptN}{\left({\displaystyle \frac{\left[{{\mathrm{NH}}_{\mathrm{4}}}^{+}\right]}{\left[{{\mathrm{NH}}_{\mathrm{4}}}^{+}\right]+\left[{{\mathrm{NO}}_{\mathrm{3}}}^{-}\right]}}\right)}^{\mathrm{0.025}})\text{32}& {\displaystyle }{\mathrm{uptNO}}_{\mathrm{3}}^{-}=\mathrm{uptN}-{\mathrm{uptNH}}_{\mathrm{4}}^{+}\text{33}& {\displaystyle \frac{\mathrm{d}\left[\mathrm{uptP}\right]}{\mathrm{d}t}}=\sum _i{\displaystyle \frac{({\mathit{\mu }}_{{\mathrm{F}}_i}+e_{{\mathrm{phy}}_i})\left[{\mathrm{F}}_i\right]}{\mathrm{C}/\mathrm{P}}}\text{34}& {\displaystyle \frac{\mathrm{d}\left[\mathrm{uptSi}\right]}{\mathrm{d}t}}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{F},dia}\left[F_{\mathrm{dia}}\right]}{\mathrm{C}/\mathrm{Si}}}\end{array}$$ Here, ${\mathit{\mu }}_{{\mathrm{F}}_i}$ is the effective growth rate of the $i$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $e_{{\mathrm{phy}}_i}$ is the excretion rate of the $i$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $\left[F_i\right]$ represents the functional macromolecule concentration of the $i$th phytoplankton species [mgC L${}^{-\mathrm{1}}$], $\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]$ and $\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$ are the concentrations of ammonium and nitrate [mgN L${}^{-\mathrm{1}}$], uptN is the uptake of nitrogen for phytoplankton growth [mgN L${}^{-\mathrm{1}}$], ${\mathrm{uptNH}}_{\mathrm{4}}^{+}$ is the uptake of NH${}_{\mathrm{4}}^{+}$ for phytoplankton growth [mgN L${}^{-\mathrm{1}}$], ${\mathrm{uptNO}}_{\mathrm{3}}^{-}$ is the uptake of NO${}_{\mathrm{3}}^{-}$ for phytoplankton growth [mgN L${}^{-\mathrm{1}}$], uptP is the uptake of phosphorus for phytoplankton growth [mgP L${}^{-\mathrm{1}}$], $\mathrm{C}/\mathrm{N}$ is the carbon to nitrogen mass ratio [mgC mgN${}^{-\mathrm{1}}$], $\mathrm{C}/\mathrm{P}$ is the carbon to phosphorus mass ratio [mgC mgP${}^{-\mathrm{1}}$], $\mathrm{C}/\mathrm{Si}$ is the carbon to silica mass ratio [mgC mgSi${}^{-\mathrm{1}}$], uptSi is the uptake of silica for diatom growth [mgSi L${}^{-\mathrm{1}}$], ${\mathit{\mu }}_{\mathrm{F},\mathrm{dia}}$ is the effective growth rate of diatoms [h${}^{-\mathrm{1}}$], and $\left[F_{\mathrm{dia}}\right]$ is the functional macromolecule (F) concentration of diatoms [mgC L${}^{-\mathrm{1}}$]

Figure 5

Cycling of nitrogen. PHY: phytoplankton species; HB: heterotrophic bacteria; ZOO: zooplankton species; AOB: ammonia-oxidizing bacteria; NOB: nitrite-oxidizing bacteria; respi.: respiration; excr.: excretion; denit: denitrification.

[Figure omitted. See PDF]

2.5.2 Release of nutrients by mineralization

The mineralization of organic matter by heterotrophic bacteria and zooplankton is achieved by its oxidation through respiration (Fig. ). The process consumes organic matter and releases nitrogen and phosphorus from the fraction that is not assimilated for growth of heterotrophic bacteria and zooplankton. $$\begin{array}{}\text{35}& {\displaystyle }\mathrm{respHB}=\sum _i(\mathrm{1}-Y_{{\mathrm{hb}}_i})b_{{\mathrm{hb}}_i}\left[{\mathrm{HB}}_i\right]\text{36}& {\displaystyle }\mathrm{respZOO}=\sum _j(\mathrm{1}-Y_{{\mathrm{zoo}}_j})b_{{\mathrm{zoo}}_j}\left[{\mathrm{ZOO}}_j\right]\text{37}& {\displaystyle }\mathrm{relN}={\displaystyle \frac{\mathrm{respHB}}{\mathrm{C}/\mathrm{N}}}+{\displaystyle \frac{\mathrm{respZOO}}{\mathrm{C}/\mathrm{N}}}\text{38}& {\displaystyle }\mathrm{relP}={\displaystyle \frac{\mathrm{respHB}}{\mathrm{C}/\mathrm{P}}}+{\displaystyle \frac{\mathrm{respZOO}}{\mathrm{C}/\mathrm{P}}}\end{array}$$ Here, respHB is the respiration of heterotrophic bacteria species [mgC L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], $Y_{{\mathrm{hb}}_i}$ is the growth yield of the $i$th heterotrophic bacteria species [–], $b_{{\mathrm{hb}}_i}$ is the effective rate of substrate uptake by the $i$th heterotrophic bacteria species (Eq. ) [h${}^{-\mathrm{1}}$], $\left[{\mathrm{HB}}_i\right]$ is the concentration of the $i$th heterotrophic bacteria species [mgC L${}^{-\mathrm{1}}$], respZOO is the respiration of zooplankton species [mgC L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], $Y_{{\mathrm{zoo}}_j}$ is the growth yield of the $j$th zooplankton species [–], $b_{{\mathrm{zoo}}_j}$ is the effective grazing rate of the $j$th zooplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $\left[{\mathrm{ZOO}}_j\right]$ is the concentration of the $j$th zooplankton species [mgC L${}^{-\mathrm{1}}$], relN is the release of nitrogen [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], $\mathrm{C}/\mathrm{N}$ is the carbon to nitrogen mass ratio [mgC mgN${}^{-\mathrm{1}}$], relP is the release of phosphorus [mgP L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], and $\mathrm{C}/\mathrm{P}$ is the carbon to phosphorus mass ratio [mgC mgP${}^{-\mathrm{1}}$].

2.5.3 Nitrification and denitrification

As mentioned in Sect. , the nitrification process (Figs.  and ) is related to the growth of AOB (ammonia-oxidizing bacteria) and NOB (nitrite-oxidizing bacteria). Growth yields ($Y_{{\mathrm{aob}}_i}$ and $Y_{{\mathrm{nob}}_j}$) are used to describe the amount of nitrogen consumed by nitrifying bacteria (Eqs.  and ). Denitrification occurs when dissolved oxygen is not present in sufficient quantity (Fig. ). $$\begin{array}{}\text{39}& {\displaystyle }{\mathrm{nitr}}_{\mathrm{aob}}=\sum _i{\displaystyle \frac{{\mathit{\mu }}_{{\mathrm{aob}}_i}}{Y_{{\mathrm{aob}}_i}}}\left[{\mathrm{AOB}}_i\right]\text{40}& {\displaystyle }{\mathrm{nitr}}_{\mathrm{nob}}=\sum _j{\displaystyle \frac{{\mathit{\mu }}_{{\mathrm{nob}}_j}}{Y_{{\mathrm{nob}}_j}}}\left[{\mathrm{NOB}}_j\right]\end{array}$$ Here, ${\mathit{\mu }}_{{\mathrm{aob}}_i}$ and ${\mathit{\mu }}_{{\mathrm{nob}}_j}$ are the effective growth rates of the $i$th AOB species (Eq. ) and the $j$th NOB species (Eq. ) [h${}^{-\mathrm{1}}$], $Y_{{\mathrm{aob}}_i}$ and $Y_{{\mathrm{nob}}_j}$ are the growth yields of the $i$th AOB species and the $j$th NOB species [mgC mgN${}^{-\mathrm{1}}$], ${\mathrm{nitr}}_{\mathrm{aob}}$ is nitrification calculated as ${\mathrm{NH}}_{\mathrm{4}}^{+}+\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}⟶{\mathrm{NO}}_{\mathrm{2}}^{-}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}+\mathrm{2}{\mathrm{H}}^{+}$ [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], ${\mathrm{nitr}}_{\mathrm{nob}}$ is nitrification calculated as ${\mathrm{NO}}_{\mathrm{2}}^{-}+\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}⟶{\mathrm{NO}}_{\mathrm{3}}^{-}$ [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], and $\left[{\mathrm{AOB}}_i\right]$ and $\left[{\mathrm{NOB}}_j\right]$ are the biomass concentrations of the $i$th AOB species and the $j$th NOB species [mgC L${}^{-\mathrm{1}}$].

The budgets of ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ can then be established. $$\begin{array}{}\text{41}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]}{\mathrm{d}t}}=-\mathrm{denit}+{\mathrm{nitr}}_{\mathrm{nob}}-{\displaystyle \frac{{\mathrm{uptNO}}_{\mathrm{3}}^{-}}{\mathrm{d}t}}\text{42}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}{\mathrm{d}t}}=\mathrm{relN}-{\mathrm{nitr}}_{\mathrm{aob}}-{\displaystyle \frac{{\mathrm{uptNH}}_{\mathrm{4}}^{+}}{\mathrm{d}t}}\text{43}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]}{\mathrm{d}t}}={\mathrm{nitr}}_{\mathrm{aob}}-{\mathrm{nitr}}_{\mathrm{nob}}\end{array}$$ Here, denit is denitrification [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], ${\mathrm{nitr}}_{\mathrm{nob}}$ is nitrification by NOB (Eq. ) [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], ${\mathrm{nitr}}_{\mathrm{aob}}$ is nitrification by AOB (Eq. ) [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], $\frac{{\mathrm{uptNO}}_{\mathrm{3}}^{-}}{\mathrm{d}t}$ is the uptake of NO${}_{\mathrm{3}}^{-}$ by phytoplankton growth (Eq. ) [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], relN is the release of nitrogen by respiration of heterotrophic bacteria and zooplankton (Eq. ) [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], and $\frac{{\mathrm{uptNH}}_{\mathrm{4}}^{+}}{\mathrm{d}t}$ is the uptake of NH${}_{\mathrm{4}}^{+}$ by phytoplankton growth (Eq. ) [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$].

2.5.4 Phosphate adsorption desorption

Orthophosphate (PO${}_{\mathrm{4}}^{\mathrm{3}-}$) is released by mineralization and taken up by phytoplankton exactly as inorganic nitrogen (Fig. ). Once released in the water column, however, orthophosphates are subject to a process of adsorption–desorption on mineral suspended solids (MSSs) to form PIP (particulate inorganic phosphorus). In addition, the impact of sediment dynamics on P fluxes should be considered in future work with unified RIVE. showed that P fluxes are mainly driven by hydrological conditions and sediment-related processes in the Seine River system.

Figure 6

Phosphorus and silica dynamics. HD: dissolved high-weight polymer; PHY: phytoplankton; PIP: particulate inorganic phosphorus; MSS: mineral suspended solids; BSi: biogenic silica; DSi: dissolved silica; adsorp.: adsorption; desorp.: desorption; sink.: sinking; disso. dissolution.

[Figure omitted. See PDF]

The process is represented according to an instantaneous hyperbolic equilibrium relationship of the form 44$${\displaystyle \frac{\left[\mathrm{PIP}\right]}{\left[\mathrm{MSS}\right]}}=P_{\mathrm{ac}}\times {\displaystyle \frac{\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}{\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]+K_{\mathrm{ps}}}},$$ with $\frac{\left[\mathrm{PIP}\right]}{\left[\mathrm{MSS}\right]}$ the inorganic P content of MSS [mgP mgMSS${}^{-\mathrm{1}}$], [PIP] and [MSS] the concentrations of PIP and MSS [mgP L${}^{-\mathrm{1}}$ and mgMSS L${}^{-\mathrm{1}}$], $P_{\mathrm{ac}}$ the maximum adsorption capacity of MSS [mgP mgMSS${}^{-\mathrm{1}}$], $\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]$ the concentration of orthophosphate [mgP L${}^{-\mathrm{1}}$], and $K_{\mathrm{ps}}$ the half-saturation adsorption constant [mgP L${}^{-\mathrm{1}}$].

Considering this equilibrium to be instantaneously reached implies that a relationship exists between the variables PIP, MSS, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, and TIP (total inorganic phosphorus): 45$$\left[\mathrm{TIP}\right]=\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]+\left[\mathrm{PIP}\right].$$ This equilibrium relationship can be written as 46$$\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]={\displaystyle \frac{\begin{array}{c}\left(\left[\mathrm{TIP}\right]-P_{\mathrm{ac}}\times \left[\mathrm{MSS}\right]-K_{\mathrm{ps}}\right)+\\ \sqrt{{\left(-\left[\mathrm{TIP}\right]+P_{\mathrm{ac}}\times \left[\mathrm{MSS}\right]+K_{\mathrm{ps}}\right)}^{\mathrm{2}}+\mathrm{4}\times \left[\mathrm{TIP}\right]\times K_{\mathrm{ps}}}\end{array}}{\mathrm{2}}},$$ with $\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]$ the concentration of orthophosphate [mgP L${}^{-\mathrm{1}}$], [TIP] and [MSS] the concentrations of TIP and MSS [mgP L${}^{-\mathrm{1}}$ and mgMSS L${}^{-\mathrm{1}}$], $P_{\mathrm{ac}}$ the maximal adsorption capacity of MSS [mgP mgMSS${}^{-\mathrm{1}}$], and $K_{\mathrm{ps}}$ the half-saturation adsorption constant [mgP L${}^{-\mathrm{1}}$].

2.5.5 Silica dynamics

Dissolved silica (DSi) is produced by the dissolution of dead frustules of diatoms (designated as biogenic silica, BSi). Rock weathering contributes also dissolved silica, while it is considered null in unified RIVE v1.0. DSi is taken up by the growth of diatoms (Fig. ). Biogenic silica is produced by the lysis and grazing of diatoms; it settles down and dissolves according to first-order kinetics, dependent on water temperature : $$\begin{array}{}\text{47}& {\displaystyle }{\mathrm{BSi}}_{\mathrm{disso}.}={\mathrm{Kb}}_{\mathrm{Si}}\times \left[\mathrm{BSi}\right],\text{48}& {\displaystyle }{\mathrm{Kb}}_{\mathrm{Si}}={\mathrm{Kb}}_{\mathrm{Si}\mathrm{20}}\times {\mathrm{ftp}}_{\mathrm{Si}}\left(T\right),\text{49}& {\displaystyle }{\mathrm{ftp}}_{\mathrm{Si}}\left(T\right)=\exp \left({\displaystyle \frac{\mathrm{60}\mathrm{000}}{\mathrm{8.314}}}\times \left({\displaystyle \frac{\mathrm{1}}{\mathrm{275}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{273}+T}}\right)\right),\end{array}$$ with Kb${}_{\mathrm{Si}\mathrm{20}}$ the dissolution rate of biogenic silica at 20 ${}^{\circ }$C [h${}^{-\mathrm{1}}$], [BSi] the concentration of biogenic silica [mgSi L${}^{-\mathrm{1}}$], and $T$ the water temperature [${}^{\circ }$C].

2.6 Dissolved oxygen

Dissolved oxygen is especially influenced by photosynthesis and respiration. Re-aeration at the water–air interface is also included in the unified RIVE v1.0 model. Sediment dynamics are important for sediment oxygen demand (not shown here). showed that benthic respiration accounts for one-third of the total Seine River respiration. Relevant efforts toward sediment dynamics need to be made in future work, which is not the focus of this study. An oxygen budget can then be established (Eq. ). $$\begin{array}{}\text{50}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\mathrm{d}t}}& =\mathrm{rea}+{\displaystyle \frac{\mathrm{32}}{\mathrm{12}}}\left(\sum _i\left(p_{{\mathrm{phy}}_i}-r_{{\mathrm{phy}}_i}\right)\left[{\mathrm{F}}_i\right]-\mathrm{respHB}\right.\\ & \left.-\sum _j{r}_{{\mathrm{zoo}}_j}\left[{\mathrm{ZOO}}_j\right]\right)-{\displaystyle \frac{\mathrm{32}}{\mathrm{14}}}\left({\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\mathrm{nitr}}_{\mathrm{aob}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{nitr}}_{\mathrm{nob}}\right)\end{array}\text{51}& {\displaystyle }\mathrm{rea}={\displaystyle \frac{k_{\mathrm{rea}}}{\mathrm{depth}}}\left(\right[{\mathrm{O}}_{\mathrm{2}}{]}_{\mathrm{sat}}-\left[{\mathrm{O}}_{\mathrm{2}}\right])\end{array}$$ Here, $k_{\mathrm{rea}}$ is the re-aeration coefficient [m h${}^{-\mathrm{1}}$], “depth” indicates water height [m], $[{\mathrm{O}}_{\mathrm{2}}{]}_{\mathrm{sat}}$ is the saturated concentration of dissolved oxygen in water [mgO${}_{\mathrm{2}}$ L${}^{-\mathrm{1}}$], $\left[{\mathrm{O}}_{\mathrm{2}}\right]$ is the concentration of dissolved oxygen in water [mgO${}_{\mathrm{2}}$ L${}^{-\mathrm{1}}$], $\frac{\mathrm{32}}{\mathrm{12}}$ is the molar mass ratio between dissolved oxygen and carbon [mgO${}_{\mathrm{2}}$ mgC${}^{-\mathrm{1}}$], $p_{{\mathrm{phy}}_i}$ is the photosynthesis rate of the $i$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $r_{{\mathrm{phy}}_i}$ is the respiration rate of the $i$th phytoplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $\left[{\mathrm{F}}_i\right]$ is the functional biomass concentration of the $i$th phytoplankton species [mgC L${}^{-\mathrm{1}}$], respHB is the respiration of all heterotrophic bacteria species (Eq. ) [mgC L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], $r_{{\mathrm{zoo}}_j}$ is the respiration rate of the $j$th zooplankton species (Eq. ) [h${}^{-\mathrm{1}}$], $\left[{\mathrm{ZOO}}_j\right]$ is the biomass concentration of the $j$th zooplankton species [mgC L${}^{-\mathrm{1}}$], $\frac{\mathrm{32}}{\mathrm{14}}$ is the molar mass ratio between dissolved oxygen and nitrogen [mgO${}_{\mathrm{2}}$ mgN${}^{-\mathrm{1}}$], nitr${}_{\mathrm{aob}}$ is nitrification to produce nitrite by oxidizing ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (Eq. , with $\frac{\mathrm{3}}{\mathrm{2}}$ the stoichiometric coefficient) [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$], and nitr${}_{\mathrm{nob}}$ is nitrification to produce nitrate by oxidizing ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (Eq. , with $\frac{\mathrm{1}}{\mathrm{2}}$ the stoichiometric coefficient) [mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$].

2.7 Inorganic carbon

An inorganic carbon module is implemented in unified RIVE v1.0. The carbonate system is described by a set of equations (named the CO${}_{\mathrm{2}}$ module) based on a previous representation provided by and adapted for freshwater environments . In this module, four state variables are defined: dissolved inorganic carbon – DIC, total alkalinity – TA, acidity – pH, and aqueous carbon dioxide – CO${}_{\mathrm{2}}$((aq)).

2.7.1

CO${}_{\mathrm{2}}$ flux at the air–water interface

The DIC is defined as the sum of three dissolved carbonate species: 52$$\left[\mathrm{DIC}\right]=\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}\right]+\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right].$$

The calculation of pH is derived from using TA and DIC. Then the aqueous carbon dioxide (CO${}_{\mathrm{2}}$ (aq)) is derived from the carbonate chemical equilibrium using DIC and pH . 53$$\left[{\mathrm{CO}}_{\mathrm{2}}\right(\mathrm{aq}\left)\right]={\displaystyle \frac{\left[\mathrm{DIC}\right]\frac{\left[{\mathrm{H}}^{+}\right]}{K_{\mathrm{1}}}}{\left(\mathrm{1}+\frac{\left[{\mathrm{H}}^{+}\right]}{K_{\mathrm{1}}}+\frac{K_{\mathrm{2}}}{\left[{\mathrm{H}}^{+}\right]}\right)}}$$ Here, $K_{\mathrm{1}}$ and $K_{\mathrm{2}}$ are the equilibrium constants of carbonate equilibrium reactions [mol L${}^{-\mathrm{1}}$], $\left[{\mathrm{H}}^{+}\right]$ is the concentration of hydrogen ions with pH $=-\log \left(\right[H^{+}\left]\right)$ [mol L${}^{-\mathrm{1}}$], and [DIC] is the concentration of dissolved inorganic carbon [mgC L${}^{-\mathrm{1}}$].

The flux of CO${}_{\mathrm{2}}$ at the water–air interface ($F_{{\mathrm{CO}}_{\mathrm{2}}}$, gC m${}^{-\mathrm{2}}$ h${}^{-\mathrm{1}}$) is calculated based on Fick's first law with a gas transfer velocity of CO${}_{\mathrm{2}}$ ($k_{{\mathrm{CO}}_{\mathrm{2}}}$). 54$$F_{{\mathrm{CO}}_{\mathrm{2}}}=k_{{\mathrm{CO}}_{\mathrm{2}}}\left(\right[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sat}\right)]-[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\left]\right)$$ Here, $k_{\mathrm{co}\mathrm{2}}$ is the gas transfer velocity of CO${}_{\mathrm{2}}$ [m h${}^{-\mathrm{1}}$], $\left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sat}\right)\right]$ is the solubility of CO${}_{\mathrm{2}}$ in water calculated based on Henry's law [mgC L${}^{-\mathrm{1}}$], and $\left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]$ is the aqueous carbon dioxide concentration [mgC L${}^{-\mathrm{1}}$].

The gas transfer velocity of CO${}_{\mathrm{2}}$ ($k_{\mathrm{co}\mathrm{2}}$) depends on water temperature and $k_{\mathrm{600}}$ (gas transfer velocity of CO${}_{\mathrm{2}}$ for a Schmidt number of 600, corresponding to a temperature of 20 ${}^{\circ }$C in freshwater). According to , , and , the gas transfer velocity of CO${}_{\mathrm{2}}$ ($k_{\mathrm{co}\mathrm{2}}$) at water temperature $T$ (${}^{\circ }$C) can be calculated as 55$${\displaystyle }{k}_{\mathrm{co}\mathrm{2}}=k_{\mathrm{600}}\sqrt{{\displaystyle \frac{\mathrm{600}}{S{c}_{{\mathrm{CO}}_{\mathrm{2}}}\left(T\right)}}},$$ where $k_{\mathrm{600}}$ (m h${}^{-\mathrm{1}}$) is the gas transfer velocity of CO${}_{\mathrm{2}}$ for a Schmidt number of 600, and $S{c}_{{\mathrm{CO}}_{\mathrm{2}}}\left(T\right)$ is the Schmidt number (dimensionless) calculated with the water temperature in degrees Celsius (${}^{\circ }$C). The $S{c}_{{\mathrm{CO}}_{\mathrm{2}}}\left(T\right)$ can be determined as,= 56$$\begin{array}{rl}S{c}_{{\mathrm{CO}}_{\mathrm{2}}}\left(T\right)& =\mathrm{1911.1}-\mathrm{118.11}T+\mathrm{3.4527}{T}^{\mathrm{2}}\\ & -\mathrm{0.04132}{T}^{\mathrm{3}}.\end{array}$$

2.7.2 Budgets of TA and DIC

Processes such as respiration, photosynthesis, nitrification, denitrification, and input flows affect TA and DIC. The unified RIVE v1.0 considers these processes explicitly. $$\begin{array}{}\text{57}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dTA}}{\mathrm{d}t}}& =\left[\left({\displaystyle \frac{\mathrm{14}}{\mathrm{106}}}\times {\displaystyle \frac{(\mathrm{respPHY}+\mathrm{respHB}+\mathrm{respZOO})}{\mathrm{12}}}\right)\right.\\ & +{\displaystyle \frac{(\mathrm{denit}-\mathrm{2}\cdot {\mathrm{nitr}}_{\mathrm{aob}})}{\mathrm{14}}}+\left({\displaystyle \frac{\mathrm{17}}{\mathrm{106}}}\times {\displaystyle \frac{{\mathrm{uptNO}}_{\mathrm{3}}^{-}}{\mathrm{uptN}}}-{\displaystyle \frac{\mathrm{15}}{\mathrm{106}}}\right.\\ & \left.\left.\times {\displaystyle \frac{{\mathrm{uptNH}}_{\mathrm{4}}^{+}}{\mathrm{uptN}}}\right)\times {\displaystyle \frac{\sum ({\mathit{\mu }}_{{\mathrm{F}}_i}+e_{{\mathrm{phy}}_i})\left[{\mathrm{F}}_i\right]}{\mathrm{12}}}\right]\\ & \times \mathrm{1000}+{\mathrm{TA}}_{\mathrm{Net}\mathit{\_}\mathrm{Input}}\end{array}\text{58}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dDIC}}{\mathrm{d}t}}& =(\mathrm{respPHY}+\mathrm{respHB}+\mathrm{respZoo})+\mathrm{denit}\\ & \times {\displaystyle \frac{\mathrm{12}}{\mathrm{14}}}\times {\displaystyle \frac{\mathrm{5}}{\mathrm{4}}}-\sum p_{{\mathrm{phy}}_i}\left[{\mathrm{F}}_i\right]+{\displaystyle \frac{F_{{\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{depth}}}\\ & +{\mathrm{DIC}}_{\mathrm{Net}\mathit{\_}\mathrm{Input}}\end{array}\end{array}$$ Here, TA${}_{\mathrm{Net}\mathit{\_}\mathrm{input}}$ ($\mathrm{\mu }$mol L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$) and DIC${}_{\mathrm{Net}\mathit{\_}\mathrm{input}}$ (mgC L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$) are the net input fluxes. The respiration of all phytoplankton, bacteria, and zooplankton species (respPHY, respHB, respZOO; mgC L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$) transforms organic carbon to CO${}_{\mathrm{2}}$ by full oxidization. The denitrification (denit, mgN L${}^{-\mathrm{1}}$ h${}^{-\mathrm{1}}$) is also considered in the calculation of TA and DIC. $F_{{\mathrm{CO}}_{\mathrm{2}}}$ (gC m${}^{-\mathrm{2}}$ h${}^{-\mathrm{1}}$) is the CO${}_{\mathrm{2}}$ flux at the air–water interface, and “depth” is the water depth (m). $\frac{\mathrm{14}}{\mathrm{106}},\frac{\mathrm{17}}{\mathrm{106}},\frac{\mathrm{15}}{\mathrm{106}}$, and $\frac{\mathrm{5}}{\mathrm{4}}$ are the stoichiometry coefficients of biogeochemical processes .

2.8 Kinetic parameters in the unified RIVE model

A total of 120 parameters are used to describe the aforementioned processes considering three phytoplankton species and two heterotrophic bacteria species. Some of them depend on water temperature and are calculated with a water temperature function (Eq. ). Their definitions and reference values are provided in Appendix .

3 Results

3.1 Digital implementation with Python 3 (pyRIVE 1.0) or ANSI C (C-RIVE 0.32)

The above unified governing equations are implemented in Python 3 to create pyRIVE 1.0 (10.48579/PRO/Z9ACP1; ) and in ANSI C to create C-RIVE 0.32 (10.5281/zenodo.7849609; ), respectively. A Jupyter Notebook is used for pedagogical exercises with pyRIVE 1.0, while C-RIVE 0.32 needs to be compiled with gcc under a Linux or a MAC OS operating system. In addition, the user interface of C-RIVE 0.32 uses its own parser based on flex and bison, which allows the software to read ASCII files.

In practice, the number of living species is predefined in pyRIVE 1.0, while we have the ability to define as many species as desired in C-RIVE 0.32 (Table ). For instance, three communities of phytoplankton (DIA: diatoms; GRA: green algae; CYA: cyanobacteria), two populations of heterotrophic bacteria distinct in their growth rate and size (small one SHB and large one LHB; ), and two zooplankton communities (ZOR: rotifer and ZOC, micro-crustaceans; ) are predefined in pyRIVE 1.0.

In addition, the TIP (total inorganic phosphorus) is considered to be a state variable in pyRIVE 1.0. PO${}_{\mathrm{4}}^{\mathrm{3}-}$ and PIP are derived from it according to Eqs. () and (). TIP is subject to release by heterotrophic bacteria and zooplankton respiration (Eq. ), uptake by phytoplankton, and settling of PIP together with MSS. However, the PO${}_{\mathrm{4}}^{\mathrm{3}-}$ is treated as a state variable and released by respiration (Eq. ) in C-RIVE 0.32 and only PIP (particulate inorganic phosphorus) is derived from the Eq. ().

Table 1

Number of living species defined in pyRIVE 1.0 and C-RIVE 0.32 which implement the unified RIVE v1.0.

Species	PHY	HB	AOB	NOB	ZOO
pyRIVE 1.0	3	2	1	1	2
C-RIVE 0.32	User-defined	User-defined	User-defined	User-defined	User-defined

Table 2

Initial concentrations and boundary conditions.

Species	Description	$C_{\mathrm{init}}$	$C_{\mathrm{boundary}}$	Unit
SHB	Small heterotrophic bacteria	0.005	0.005	[mgC L${}^{-\mathrm{1}}$]
LHB	Large heterotrophic bacteria	0.004	0.004	[mgC L${}^{-\mathrm{1}}$]
AOB	Ammonia-oxidizing bacteria	0.001	0.001	[mgC L${}^{-\mathrm{1}}$]
NOB	Nitrite-oxidizing bacteria	0.0002	0.0002	[mgC L${}^{-\mathrm{1}}$]
DIA	Diatoms	0.447	0.447	[mgC L${}^{-\mathrm{1}}$]
GRA	Green algae	0.539	0.539	[mgC L${}^{-\mathrm{1}}$]
CYA	Cyanobacteria	0.662	0.662	[mgC L${}^{-\mathrm{1}}$]
ZOR	Rotifer	$\mathrm{9.33}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{9.33}\times {\mathrm{10}}^{-\mathrm{5}}$	[mgC L${}^{-\mathrm{1}}$]
ZOC	Micro-crustaceans	$\mathrm{9.33}\times {\mathrm{10}}^{-\mathrm{6}}$	$\mathrm{9.33}\times {\mathrm{10}}^{-\mathrm{6}}$	[mgC L${}^{-\mathrm{1}}$]
SMS	Small monomeric substrate	0.036	0.036	[mgC L${}^{-\mathrm{1}}$]
DOM${}_{\mathrm{1}}$	Rapidly biodegradable dissolved organic matter	0.022	0.022	[mgC L${}^{-\mathrm{1}}$]
DOM${}_{\mathrm{2}}$	Slowly biodegradable dissolved organic matter	0.174	0.174	[mgC L${}^{-\mathrm{1}}$]
DOM${}_{\mathrm{3}}$	Dissolved refractory organic matter	1.625	1.625	[mgC L${}^{-\mathrm{1}}$]
POM${}_{\mathrm{1}}$	Rapidly biodegradable particulate organic matter	0.005	0.005	[mgC L${}^{-\mathrm{1}}$]
POM${}_{\mathrm{2}}$	Slowly biodegradable particulate organic matter	0.021	0.021	[mgC L${}^{-\mathrm{1}}$]
POM${}_{\mathrm{3}}$	Particulate refractory organic matter	0.107	0.107	[mgC L${}^{-\mathrm{1}}$]
NH4	Ammonium	1.5	1.5	[mgN L${}^{-\mathrm{1}}$]
NO2	Nitrite	0.016	0.016	[mgN L${}^{-\mathrm{1}}$]
NO3	Nitrate	0.941	0.941	[mgN L${}^{-\mathrm{1}}$]
TIP	Total inorganic phosphorus	0.2	0.2	[mgP L${}^{-\mathrm{1}}$]
DSi	Dissolved silica	3.090	3.090	[mgSi L${}^{-\mathrm{1}}$]
MSS	Mineral suspended solids	2.611	2.611	[mgMSS L${}^{-\mathrm{1}}$]
OXY	Dissolved oxygen	9.446	9.446	[mgO${}_{\mathrm{2}}$ L${}^{-\mathrm{1}}$]
TA	Total alkalinity	5291	5291	[$\mathrm{\mu }$mol L${}^{-\mathrm{1}}$]
DIC	Dissolved inorganic carbon	62.728	62.728	[mgC L${}^{-\mathrm{1}}$]
CO${}_{\mathrm{2}}$(aq)	Aqueous carbon dioxide	0.343	0.343	[mgC L${}^{-\mathrm{1}}$]
pH	Acidity	8.659	8.695	[–]

Figure 7

Simulation by the HSB model (unified RIVE v1.0) of the dynamics of heterotrophic bacteria in a filtered and reinoculated sample from drainage pond water (Seine basin, France) taken in February 2021 . DOM: dissolved organic matter; SHB: small heterotrophic bacteria.

[Figure omitted. See PDF]

3.2 Modeling the organic matter degradation by unified RIVE v1.0 (HSB model)

The ability of the HSB model (Fig. ) to simulate organic matter degradation has been verified by modeling two batch experiments conducted by . Two water samples were used in the study. One sample was obtained from a drainage pond in the Seine basin, France, in February 2021 (Fig. ). The other sample was taken from an urban sewage collector in Rosny-sur-Seine, France, in February 2021 (Fig. ). These samples were incubated in the dark at a temperature of 21 ${}^{\circ }$C for a period of 45 d, during which aerobic bacteria consumed organic matter . Only DOM and bacterial biomass are measured during batch experiments and then used to show validation. The HSB model is able to effectively reproduce the concentrations of dissolved organic matter and bacterial biomass with a trial–error adjustment of its parameter values (Figs.  and ). The parameter values are kept the same for both water samples.

Figure 8

Simulation by the HSB model (unified RIVE v1.0) of the dynamics of heterotrophic bacteria in a filtered and reinoculated sample from urban sewage water (Rosny-sur-Seine, France) taken in February 2021 . DOM: dissolved organic matter; SHB: small heterotrophic bacteria; LHB: large heterotrophic bacteria.

[Figure omitted. See PDF]

Figure 9

Geometric and hydraulic description of a river stretch.

[Figure omitted. See PDF]

3.3 A river stretch simulated with unified RIVE v1.0: pyRIVE 1.0 vs. C-RIVE 0.32

A river stretch with a Strahler order of 8 (Fig. ) is designed to compare the results simulated by two versions of unified RIVE v1.0 implemented in pyRIVE 1.0 and C-RIVE 0.32. The case study allows us to compare the two versions of unified RIVE v1.0 under transient contrasting conditions (i) between species communities and (ii) temporally for each species community.

3.3.1 River stretch morphology and hydraulic conditions

The stretch measures 10 000 m long and 300 m in width. To simplify the boundary conditions, the upstream inflow and downstream outflow are fixed at 25 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, which corresponds to a residence time of 7 d. The water height is fixed at 5 m.

3.3.2 Simulation settings and evaluation strategy

The concentrations of all water quality variables of inflow are defined as their initial concentrations in the stretch and remain constant during the simulation (Table ). Since this paper focuses on the conceptualization of the unified RIVE v1.0 in the water column, no exchange between the benthic layer and water column is considered. The time step of the simulation is 6 min and a simulation period of 365 d is considered. To compare the results of the two digital implementations of unified RIVE v1.0, daily concentrations at 00:00 are plotted. Three statistical criteria (PBIAS: percent bias, %; MAD: mean absolute difference; MaAD: maximum absolute difference) are calculated to evaluate the similarity of the two sets of results. The closer the criteria are to 0, the more similar the concentrations simulated by the two software packages are (pyRIVE 1.0 and C-RIVE 0.32). $$\begin{array}{}\text{59}& {\displaystyle }\mathrm{PBIAS}=\mathrm{100}{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(C_i-{\mathrm{Py}}_i)}{{\sum }_{i=\mathrm{1}}^N{\mathrm{Py}}_i}}\text{60}& {\displaystyle }\mathrm{MAD}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N|C_i-{\mathrm{Py}}_i|}{N}}\text{61}& {\displaystyle }\mathrm{MaAD}=\max \left(\right|C_i-{\mathrm{Py}}_i\left|\right)\end{array}$$ Here, $C_i$ represents the concentrations simulated by C-RIVE 0.32 (in ANSI C) and ${\mathrm{Py}}_i$ those simulated by pyRIVE 1.0 (in Python 3). $N$ is the number of values.

Figure 10

Simulated concentrations of main species by pyRIVE 1.0 and C-RIVE 0.32. See Table  for their definitions.

[Figure omitted. See PDF]

Table 3

Statistical criteria for comparing the simulated variables by pyRIVE 1.0 and C-RIVE 0.32 which implement the unified RIVE v1.0. PBIAS: percent bias [%]; MAD: mean absolute difference; MaAD: maximum absolute The unit of MAD and MaAD depends on the species. It can be [mgC L${}^{-\mathrm{1}}$], [mgN L${}^{-\mathrm{1}}$], [mgP L${}^{-\mathrm{1}}$], [mgSi L${}^{-\mathrm{1}}$], or [$\mathrm{\mu }$mol L${}^{-\mathrm{1}}$].

Species	PBIAS	MAD	MaAD	Species	PBIAS	MAD	MaAD
SHB	$-\mathrm{0.4}$	$\mathrm{2.07}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{1.40}\times {\mathrm{10}}^{-\mathrm{4}}$	LHB	$-\mathrm{0.4}$	$\mathrm{1.10}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{4.27}\times {\mathrm{10}}^{-\mathrm{4}}$
AOB	0	$\mathrm{1.78}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{6.95}\times {\mathrm{10}}^{-\mathrm{5}}$	NOB	0	$\mathrm{4.85}\times {\mathrm{10}}^{-\mathrm{6}}$	$\mathrm{1.98}\times {\mathrm{10}}^{-\mathrm{5}}$
DIA	$-\mathrm{0.1}$	$\mathrm{1.07}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{3.63}\times {\mathrm{10}}^{-\mathrm{3}}$	GRA	$-\mathrm{0.1}$	$\mathrm{4.98}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{1.87}\times {\mathrm{10}}^{-\mathrm{3}}$
CYA	$-\mathrm{1.2}$	$\mathrm{7.26}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{3.15}\times {\mathrm{10}}^{-\mathrm{2}}$	ZOR	$-\mathrm{0.2}$	$\mathrm{2.26}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{2.89}\times {\mathrm{10}}^{-\mathrm{3}}$
ZOC	0	$\mathrm{5.25}\times {\mathrm{10}}^{-\mathrm{9}}$	$\mathrm{1.69}\times {\mathrm{10}}^{-\mathrm{7}}$	SMS	$-\mathrm{0.6}$	$\mathrm{6.65}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{4.31}\times {\mathrm{10}}^{-\mathrm{3}}$
DOM${}_{\mathrm{1}}$	0.1	$\mathrm{6.54}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{3.61}\times {\mathrm{10}}^{-\mathrm{4}}$	DOM${}_{\mathrm{2}}$	$-\mathrm{0.1}$	$\mathrm{3.94}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{1.34}\times {\mathrm{10}}^{-\mathrm{3}}$
DOM${}_{\mathrm{3}}$	0	$\mathrm{4.06}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{2.06}\times {\mathrm{10}}^{-\mathrm{3}}$	OXY	0	$\mathrm{6.78}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{3.07}\times {\mathrm{10}}^{-\mathrm{2}}$
POM${}_{\mathrm{1}}$	$-\mathrm{0.3}$	$\mathrm{4.68}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{2.27}\times {\mathrm{10}}^{-\mathrm{3}}$	POM${}_{\mathrm{2}}$	$-\mathrm{0.3}$	$\mathrm{7.80}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{3.97}\times {\mathrm{10}}^{-\mathrm{3}}$
POM${}_{\mathrm{3}}$	$-\mathrm{0.2}$	4$.06\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{2.06}\times {\mathrm{10}}^{-\mathrm{3}}$	NH4	0	$\mathrm{4.66}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{5.28}\times {\mathrm{10}}^{-\mathrm{4}}$
NO2	0	$\mathrm{2.13}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{3.42}\times {\mathrm{10}}^{-\mathrm{4}}$	NO3	0	$\mathrm{4.04}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{1.90}\times {\mathrm{10}}^{-\mathrm{3}}$
TIP	0	$\mathrm{1.45}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{3.62}\times {\mathrm{10}}^{-\mathrm{4}}$	DSi	0	$\mathrm{1.81}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{8.46}\times {\mathrm{10}}^{-\mathrm{3}}$
TA	0	$\mathrm{4.39}\times {\mathrm{10}}^{-\mathrm{2}}$	$\mathrm{2.65}\times {\mathrm{10}}^{-\mathrm{1}}$	pH	0	$\mathrm{1.14}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{7.22}\times {\mathrm{10}}^{-\mathrm{3}}$
DIC	0	$\mathrm{8.49}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{6.54}\times {\mathrm{10}}^{-\mathrm{2}}$	CO${}_{\mathrm{2}}$(aq)	0	$\mathrm{2.24}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{1.65}\times {\mathrm{10}}^{-\mathrm{2}}$

3.3.3 Simulated concentrations of water quality variables

The concentrations simulated by pyRIVE 1.0 and C-RIVE 0.32 are very similar (and superimposed) for all water quality variables (Fig. ). A maximum absolute difference (MaAD) of 0.0307 mgO${}_{\mathrm{2}}$ L${}^{-\mathrm{1}}$, which is relatively low, is obtained for dissolved oxygen concentration. The mean absolute difference (MAD) for dissolved oxygen concentration is 0.00678 mgO${}_{\mathrm{2}}$ L${}^{-\mathrm{1}}$ (Table ) and the corresponding percent bias (PBIAS) is 0 %. The MaAD of 0.0307 mgO${}_{\mathrm{2}}$ L${}^{-\mathrm{1}}$ for dissolved oxygen is the cause of the depletion of CYA S (small precursors S of cyanobacteria, Fig. ) at the beginning of the simulation (not shown here). To correct this depletion of CYA S, the growth of functional macromolecules (CYA F) is reduced according to the availability of CYA S in C-RIVE 0.32. This is why the simulated concentrations of CYA (cyanobacteria) depict a MaAD of 0.0321 mgC L${}^{-\mathrm{1}}$ between pyRIVE 1.0 and C-RIVE 0.32. Due to this auto-correction in C-RIVE 0.32, the simulated concentrations of CYA by C-RIVE 0.32 are slightly smaller than those simulated by pyRIVE 1.0 (PBIAS $=-\mathrm{1.2}\mathit{\%}$). The values of PBIAS also indicate the similarity between the concentrations simulated by pyRIVE 1.0 and C-RIVE 0.32. Except for CYA, the discrepancies of other variables are extremely low compared to their concentrations (PBIAS $\le \mathrm{0.6}\mathit{\%}$). More than half of simulated variables have a PBIAS of 0 %.

4 Discussion

The results show the ability of the unified RIVE v1.0 to correctly simulate the organic matter degradation and the similarity of its two digital implementations (pyRIVE 1.0 and C-RIVE 0.32). Here, we discuss the biogeochemical cycling simulated by unified RIVE v1.0 in the water column (Sect. ), the model limitations, and the future developments (Sect. ), as well as its benefits for the scientific community (Sect. ).

Figure 11

Budget fluxes of DIA (ton C d${}^{-\mathrm{1}}$). Grazing: zooplankton grazing fluxes; net production: Eq. (); net input: input flux–output flux; net change: daily variation of DIA in the river stretch.

[Figure omitted. See PDF]

Figure 12

Simulated water temperature, active irradiance, and zoom of active irradiance for days 111–114.

[Figure omitted. See PDF]

4.1 Biogeochemical cycling in water column simulated by unified RIVE v1.0

The unified RIVE v1.0 simulates the dynamics of microorganisms involving biogeochemical cycling, although the boundary conditions are defined as constant for modeling a river stretch (Fig. ). Here we interpret the dynamics of diatoms (DIA) and large heterotrophic bacteria (LHB). For this purpose, the budget fluxes of DIA and LHB are calculated.

Figure 13

Large heterotrophic bacteria (LHB) dynamics and simulated concentrations of small monomeric substrate (SMS). Net change: daily variation of LHB in river stretch; net input: input flux–output flux.

[Figure omitted. See PDF]

The decrease in DIA biomass from day 1 to day 15 is related to the low water temperature and low active irradiance, which limit its photosynthesis (Figs. , ). The optimal temperature for the growth of DIA is 21 ${}^{\circ }$C, while the water temperature is lower than 3 ${}^{\circ }$C (Fig. ). The low photosynthesis rate leads to negative net production (Fig. , green line), which is the difference between the fluxes of photosynthesis and the combined fluxes of respiration, mortality, and excretion. During this period, while the input factors play a positive role, the net change in DIA is still negative (Fig. , black line). Over the following days, as the water temperature and active irradiance increase (Fig. ), the net production shows an increase. However, it still remains negative. The net change in DIA shifts to a positive direction due to a combination of net input and net production, leading to a simulated increase in DIA biomass. This trend continues until day 130 when the maximum DIA biomass is reached (Fig. ). The decline in DIA biomass simulated from day 130 onwards is due to a combination of factors. Firstly, the input factor could be contributing to the decrease when the DIA biomass exceeds the concentration of DIA in input flow (0.447 mgC L${}^{-\mathrm{1}}$, Table , Fig. ). Additionally, the net production rate also plays a role (Fig. ). Although the photosynthesis rate increases with water temperature and active irradiance until day 179 (not shown here), it is not enough to compensate for the other processes occurring in the diatom population (days 150–170), resulting in an overall decrease in biomass. Despite the positive contributions of net input and net production to DIA biomass around day 175, a significant decrease in biomass occurred due to zooplankton grazing (Fig. , red line). Two factors impact the zooplankton dynamics: water temperature and the half-saturation constant of grazing (Eq. ). The optimal temperature for zooplankton is 25 ${}^{\circ }$C and the half-saturation constant of grazing for zooplankton is set to 0.4 mgC L${}^{-\mathrm{1}}$. Then, equilibrium of DIA biomass is simulated until day 260 (Fig. ), which means that the net production and net input in DIA biomass are balanced by the grazing of zooplankton. The input in DIA biomass primarily contributes to the increase in DIA biomass from day 260 (Fig. ). As the water temperature and active irradiance decrease during this time, the net production of DIA decreases and changes to negative by day 292.

The growth rate of large heterotrophic bacteria (LHB) increases (Fig. ) with the increase in water temperature, causing a rise in LHB biomass until day 170 (Fig. ). The fast decrease in small monomeric substrate (SMS) around day 175, synchronized with the grazing of zooplankton (Fig. ), causes a decrease in growth rate of LHB. Its mortality rate is not impacted (not shown here). Consequently, this leads to a significant reduction in LHB biomass around day 175 (Fig. ). The biomass of LHB remains stable until day 260. After that, it increases in conjunction with the rise in SMS concentration, which is synchronized with the increase in phytoplankton biomass.

4.2 Complexity and strengths of the RIVE model

Complexity can be understood in terms of the large number of variables represented and interacting with each other. The RIVE model is a multi-element, multi-form model and the kinetics it represents inevitably incorporate a large number of parameters. This is especially true as the RIVE model has opted for an explicit representation of the living communities (bacteria, phytoplankton, zooplankton, etc.) involved in the carbon and nutrient cycles. The model has thus become more complex over time and the addition of new processes (and therefore new parameters) has, as far as possible, been systematically based on experimental work in the laboratory or in the field to reduce the ranges of uncertainty around the kinetic parameters.

The RIVE model is designed as a tool for generating knowledge about the functioning of freshwater ecosystems, and therefore it documents a large number of biogeochemical processes, whether they are expressed weakly or strongly in a given freshwater ecosystem. The underlying hypothesis is that environmental factors control the intensity with which the various processes involved in the overall functioning of a hydrosystem are expressed.

Nevertheless, some work has specifically focused on analyzing the influence of RIVE parameters, particularly those controlling oxygen levels . This work identified key physical and physiological parameters. Based on the result of sensitivity analysis, a continuous oxygen data assimilation scheme has been developed (Prose-PA, ). The data assimilation allows determining the physiological properties of microorganisms by integrating the associated uncertainties over time. The recent work of has also helped to better quantify the sensitivity of oxygen to bacterial kinetics parameters as well as those relating to the composition of organic matter with the aim of parsimonious simplification of the number of parameters.

In these two examples, RIVE (C-RIVE) biogeochemical modeling is implemented in much more complex modeling platforms (particle filter, data assimilation, etc.) and the various analyses (sensitivity, uncertainties, etc.) are also supported by an overall assessment of the performance of the model applied to the Seine River.

4.3 Model limitations and future developments of unified RIVE

Currently, the unified RIVE v1.0 presented in this paper describes only the biogeochemical processes in the water column. Comparison of benthic processes and simulations have not been investigated yet. Previous studies showed that sediment plays an important role in the metabolism of river and lakes . A unified sediment module should be further elaborated based on existing modules and implemented into unified RIVE. This sediment module will have to take into account not only the dissolved exchanges between the water column and the sediment but also the resuspension of particulates.

In addition, the unified RIVE v1.0 simulates phytoplankton dynamics, but periphyton or macrophyte development is not implemented in current versions. showed that periphyton plays a major role in carbon cycling (primary productivity) in small rivers, not only in the carbon stock fixed at the bottom of the river but also in the carbon enrichment downstream of the river. These limitations should be considered in future developments.

4.4 Benefit of the unified RIVE model

The unified RIVE provides a set of governing equations of freshwater biogeochemical processes across different software platforms, such as pyNuts-Riverstrahler , ProSe-PA , SWAT-RIVE , QUAL-NET , VEMALA V3 , and Barman , while incorporating the latest developments. The unity of the kinetics is important for facilitating and reinforcing the collaboration nationally or internationally within different research teams. Thanks to the unity property formerly pointed out by the river continuum concept , the software packages based on unified RIVE can leverage the already identified parameter values, regardless of the location in the network . This feature is of great interest to the different research teams involved in freshwater quality research, for instance river metabolism or nutrient cycling .

Open science has become increasingly popular and even indispensable in the scientific community as it allows for easier accessibility and the reproduction of scientific results. The unified RIVE project, as an open-source project, allows for the dissemination and wider use of the RIVE biogeochemical model by creating a public repository with different programming languages.

5 Conclusions

This paper presents a conceptual freshwater biogeochemistry model: unified RIVE v1.0, programmed in Python 3 and ANSI C. The degradation of organic matter by heterotrophic bacteria and the dynamics of primary producers (phytoplankton) and zooplankton including carbon cycling and nutrient cycling are described exhaustively. In unified RIVE v1.0, organic matter is degraded via bacteria activity, which is simulated by an HSB model. According to the results, the HSB model is able to model the organic matter degradation and bacterial dynamics in batch experiments. A case study is designed to compare the simulations of the two digital implementations (Python 3 for pyRIVE 1.0 and ANSI C for C-RIVE 0.32). These implementations simulate similar concentrations of all state variables including microorganisms, organic carbon, nutrients, and inorganic carbon.

The river stretch case study allows us to compare the two implementations of unified RIVE V1.0 under transient contrasting conditions involving complex biogeochemical cycles. The specific dynamics of each simulated species depend on different limiting factors. The calculation of photosynthesis of phytoplankton (diatoms, Chlorophyceae, cyanobacteria) takes into account the light that naturally presents a day–night variation. The development of diatoms specifically takes into account the dissolved silica in the simulated environment. The growth of microorganisms depends on the quantity of nutrients (primary producer, nitrifying bacteria) and the small monomeric substrate (heterotrophic bacteria). In addition, the effect of water temperature is also taken into account for the physiology of the simulated microorganism communities (photosynthesis, growth, respiration, mortality).

Finally, unified RIVE being an open-source project, contributions from the freshwater biogeochemistry community are strongly encouraged to achieve a better understanding of freshwater ecosystem functioning and further investigate the future of river systems in a changing world.

Appendix A Various implementations of the RIVE model and its applications in different freshwater systems

Table A1

Various implementations of the RIVE model and its applications in different freshwater systems.

Freshwater systems	Climates	Software platforms	References
Danube River (Romania and Bulgaria)	Continental	Riverstrahler
Day–Nhue River (Vietnam)	Tropical	Seneque-Riverstrahler
Grand Morin River (France)	Temperate	ProSe
Loire River (France)	Temperate	Grafs-Seneque/Riverstrahler
Lot River (France)	Temperate	Grafs-Seneque/Riverstrahler
Lule and Kalix rivers (Sweden)	Subarctic	Riverstrahler
Mosel River (Germany)	Temperate	Riverstrahler
Orgeval watershed (France)	Temperate	Seneque-Riverstrahler
Red River (China and Vietnam)	Tropical	Seneque-Riverstrahler
Scheldt River (Belgium and Netherlands)	Temperate	Seneque-Riverstrahler
Somme River (France)	Temperate	Seneque-Riverstrahler
Seine River (France)	Temperate	Seneque-Riverstrahler
Seine River (France)	Temperate	pyNuts-Riverstrahler
Seine River (France)	Temperate	ProSe/ProSe-PA
Zenne River (Belgium)	Temperate	Seneque-Riverstrahler
Sand-pit lake, reservoirs (France)	Temperate	Barman

Appendix B Parameter values for unified RIVE v1.0

The 120 parameter values necessary for running unified RIVE v1.0 are provided hereafter.

Table B1

Parameters related to heterotrophic bacteria.

Parameter	Description	Value	Unit
$e_{\max \mathrm{20},{\mathrm{dom}}_{\mathrm{1}},\mathrm{lhb}}^{*}$	Maximal hydrolysis rate of DOM${}_{\mathrm{1}}$ at 20 ${}^{\circ }$C related to LHB	0.75	[h${}^{-\mathrm{1}}$]
$e_{\max \mathrm{20},{\mathrm{dom}}_{\mathrm{1}},\mathrm{shb}}^{*}$	Maximal hydrolysis rate of DOM${}_{\mathrm{1}}$ at 20 ${}^{\circ }$C related to SHB	0.75	[h${}^{-\mathrm{1}}$]
$e_{\max \mathrm{20},{\mathrm{dom}}_{\mathrm{2}},\mathrm{lhb}}^{*}$	Maximal hydrolysis rate of DOM${}_{\mathrm{2}}$ at 20 ${}^{\circ }$C related to LHB	0.25	[h${}^{-\mathrm{1}}$]
$e_{\max \mathrm{20},{\mathrm{dom}}_{\mathrm{2}},\mathrm{shb}}^{*}$	Maximal hydrolysis rate of DOM${}_{\mathrm{2}}$ at 20 ${}^{\circ }$C related to SHB	0.25	[h${}^{-\mathrm{1}}$]
$b_{\max \mathrm{20},\mathrm{lhb}}^{*}$	Maximal substrate (SMS) uptake rate of LHB at 20 ${}^{\circ }$C	0.6	[h${}^{-\mathrm{1}}$]
$b_{\max \mathrm{20},\mathrm{shb}}^{*}$	Maximal substrate (SMS) uptake rate of SHB at 20 ${}^{\circ }$C	0.16	[h${}^{-\mathrm{1}}$]
$Y_{\mathrm{lhb}}$	Growth yield of LHB	0.25	[–]
$Y_{\mathrm{shb}}$	Growth yield of SHB	0.25	[–]
$k_{\mathrm{d}\mathrm{20},\mathrm{lhb}}^{*}$	Mortality rate of LHB at 20 ${}^{\circ }$C	0.05	[h${}^{-\mathrm{1}}$]
$k_{\mathrm{d}\mathrm{20},\mathrm{shb}}^{*}$	Mortality rate of SHB at 20 ${}^{\circ }$C	0.02	[h${}^{-\mathrm{1}}$]
${\mathrm{vs}}_{\mathrm{shb}}$	Sinking velocity of LHB	0.0	[m h${}^{-\mathrm{1}}$]
${\mathrm{vs}}_{\mathrm{lhb}}$	Sinking velocity of LHB	0.02	[m h${}^{-\mathrm{1}}$]
$T_{\mathrm{opt},\mathrm{shb}}$	Optimal temperature of SHB	20	${}^{\circ }$C
$T_{\mathrm{opt},\mathrm{lhb}}$	Optimal temperature of LHB	22	${}^{\circ }$C
${\mathit{\sigma }}_{\mathrm{shb}}$	Range of temperature for SHB	17	${}^{\circ }$C
${\mathit{\sigma }}_{\mathrm{lhb}}$	Range of temperature for LHB	12	${}^{\circ }$C
$K_{\mathrm{sms},\mathrm{lhb}}$	Half-saturation constant of LHB for small monomeric substrate	0.1	[mgC L${}^{-\mathrm{1}}$]
$K_{\mathrm{sms},\mathrm{shb}}$	Half-saturation constant of SHB for small monomeric substrate	0.1	[mgC L${}^{-\mathrm{1}}$]
$K_{\mathrm{dom}\mathrm{1},\mathrm{shb}}$	Half-saturation constant for DOM${}_{\mathrm{1}}$ hydrolysis related to SHB	0.25	[mgC L${}^{-\mathrm{1}}$]
$K_{\mathrm{dom}\mathrm{1},\mathrm{lhb}}$	Half-saturation constant for DOM${}_{\mathrm{1}}$ hydrolysis related to LHB	0.25	[mgC L${}^{-\mathrm{1}}$]
$K_{\mathrm{dom}\mathrm{2},\mathrm{shb}}$	Half-saturation constant for DOM${}_{\mathrm{2}}$ hydrolysis related to SHB	2.5	[mgC L${}^{-\mathrm{1}}$]
$K_{\mathrm{dom}\mathrm{2},\mathrm{lhb}}$	Half-saturation constant for DOM${}_{\mathrm{2}}$ hydrolysis related to LHB	2.5	[mgC L${}^{-\mathrm{1}}$]

${}^{*}$ Parameters depend on water temperature and are multiplied by $f\left(T\right)=\frac{e^{-\frac{(T-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}{e^{-\frac{(\mathrm{20}-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}$, where $T$ is water temperature in ${}^{\circ }$C.

Table B2

Parameters related to nitrifying bacteria.

Parameter	Description	Value	Unit
${\mathit{\mu }}_{\max \mathrm{20},\mathrm{aob}}^{*}$	Maximal growth rate of AOB at 20 ${}^{\circ }$C	0.07	[h${}^{-\mathrm{1}}$]
${\mathit{\mu }}_{\max \mathrm{20},\mathrm{nob}}^{*}$	Maximal growth rate of NOB at 20 ${}^{\circ }$C	0.05	[h${}^{-\mathrm{1}}$]
$K_{{\mathrm{o}}_{\mathrm{2}},\mathrm{aob}}$	Half-saturation constant of AOB for O${}_{\mathrm{2}}$	0.64	[mgO${}_{\mathrm{2}}$ L${}^{-\mathrm{1}}$]
$K_{{\mathrm{o}}_{\mathrm{2}},\mathrm{nob}}$	Half-saturation constant of NOB for O${}_{\mathrm{2}}$	1.088	[mgO${}_{\mathrm{2}}$ L${}^{-\mathrm{1}}$]
$K_{{\mathrm{nh}}_{\mathrm{4}},\mathrm{aob}}$	Half-saturation constant of AOB for NH${}_{\mathrm{4}}^{+}$	0.75	[mgN L${}^{-\mathrm{1}}$]
$K_{{\mathrm{no}}_{\mathrm{2}},\mathrm{nob}}$	Half-saturation constant of NOB for NO${}_{\mathrm{2}}^{-}$	0.05	[mgN L${}^{-\mathrm{1}}$]
$Y_{\mathrm{aob}}$	Growth yield of AOB	0.07	[mgC mgN${}^{-\mathrm{1}}$]
$Y_{\mathrm{nob}}$	Growth yield of AOB	0.02	[mgC mgN${}^{-\mathrm{1}}$]
$k_{\mathrm{d}\mathrm{20},\mathrm{aob}}^{*}$	Mortality rate of AOB at 20 ${}^{\circ }$C	0.005	[h${}^{-\mathrm{1}}$]
$k_{\mathrm{d}\mathrm{20},\mathrm{nob}}^{*}$	Mortality rate of NOB at 20 ${}^{\circ }$C	0.005	[h${}^{-\mathrm{1}}$]
${\mathrm{vs}}_{\mathrm{aob}}$	Sinking velocity of AOB	0.005	[m h${}^{-\mathrm{1}}$]
${\mathrm{vs}}_{\mathrm{nob}}$	Sinking velocity of NOB	0.005	[m h${}^{-\mathrm{1}}$]
$T_{\mathrm{opt},\mathrm{aob}}$	Optimal temperature of AOB	23	[${}^{\circ }$C]
${\mathit{\sigma }}_{\mathrm{aob}}$	Range of temperature for NOB	18	[${}^{\circ }$C]
$T_{\mathrm{opt},\mathrm{nob}}$	Optimal temperature of NOB	23	[${}^{\circ }$C]
${\mathit{\sigma }}_{\mathrm{nob}}$	Range of temperature for NOB	18	[${}^{\circ }$C]

${}^{*}$ Parameters depend on water temperature and are multiplied by $f\left(T\right)=\frac{e^{-\frac{(T-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}{e^{-\frac{(\mathrm{20}-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}$, where $T$ is water temperature in ${}^{\circ }$C.

Table B3

Parameters related to primary producer dynamics.

Parameter	Description	Value	Unit
$P_{\max \mathrm{20},\mathrm{dia}}^{*}$	Maximal photosynthesis rate of diatoms at 20 ${}^{\circ }$C	0.2	[h${}^{-\mathrm{1}}$]
$P_{\max \mathrm{20},\mathrm{gra}}^{*}$	Maximal photosynthesis rate of green algae at 20 ${}^{\circ }$C	0.25	[h${}^{-\mathrm{1}}$]
$P_{\max \mathrm{20},\mathrm{cya}}^{*}$	Maximal photosynthesis rate of cyanobacteria at 20 ${}^{\circ }$C	0.1	[h${}^{-\mathrm{1}}$]
${\mathit{\alpha }}_{\mathrm{dia}}$	Photosynthetic efficiency of diatoms	0.0012	[h${}^{-\mathrm{1}}$ ($\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$]
${\mathit{\alpha }}_{\mathrm{gra}}$	Photosynthetic efficiency of green algae	0.0012	[h${}^{-\mathrm{1}}$ ($\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$]
${\mathit{\alpha }}_{\mathrm{cya}}$	Photosynthetic efficiency of cyanobacteria	0.0012	[h${}^{-\mathrm{1}}$ ($\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$]
${\mathit{\beta }}_{\mathrm{dia}}$	Photoinhibition capacity of diatoms	0.0	[h${}^{-\mathrm{1}}$ ($\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$]
${\mathit{\beta }}_{\mathrm{gra}}$	Photoinhibition capacity of green algae	0.0	[h${}^{-\mathrm{1}}$ ($\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$]
${\mathit{\beta }}_{\mathrm{cya}}$	Photoinhibition capacity of cyanobacteria	0.0	[h${}^{-\mathrm{1}}$ ($\mathrm{\mu }$E m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$]
${\mathit{\eta }}_{\mathrm{base}}$	Light-extinction-related coefficient for pure water	0.2	[m${}^{-\mathrm{1}}$]
${\mathit{\eta }}_{\mathrm{chl}a}$	Light algal self-shading light extinction coefficient	0.02	[m${}^{-\mathrm{1}}$ ($\mathrm{\mu }$g chl $a$ L${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$]
${\mathit{\eta }}_{\mathrm{ss}}$	Light extinction coefficient related to suspended solid	0.042	[m${}^{-\mathrm{1}}$ (mg L${}^{-\mathrm{1}}$)${}^{-\mathrm{1}}$]
${\mathit{\mu }}_{\max \mathrm{20},\mathrm{dia}}^{*}$	Maximal growth rate of diatoms at 20 ${}^{\circ }$C	0.05	[h${}^{-\mathrm{1}}$]
${\mathit{\mu }}_{\max \mathrm{20},\mathrm{gra}}^{*}$	Maximal growth rate of green algae at 20 ${}^{\circ }$C	0.05	[h${}^{-\mathrm{1}}$]
${\mathit{\mu }}_{\max \mathrm{20},\mathrm{cya}}^{*}$	Maximal growth rate of cyanobacteria at 20 ${}^{\circ }$C	0.025	[h${}^{-\mathrm{1}}$]
$K_{\mathrm{S},\mathrm{dia}}$	Half-saturation constant for small precursors of diatoms	0.06	[–]
$K_{\mathrm{S},\mathrm{gra}}$	Half-saturation constant for small precursors of green algae	0.06	[–]
$K_{\mathrm{S},\mathrm{cya}}$	Half-saturation constant for small precursors of cyanobacteria	0.06	[–]
$K_{\mathrm{N},\mathrm{dia}}$	Half-saturation constant for nitrogen of diatoms	0.014	[mgN L${}^{-\mathrm{1}}$]
$K_{\mathrm{N},\mathrm{gra}}$	Half-saturation constant for nitrogen of green algae	0.014	[mgN L${}^{-\mathrm{1}}$]
$K_{\mathrm{N},\mathrm{cya}}$	Half-saturation constant for nitrogen of cyanobacteria	0.014	[mgN L${}^{-\mathrm{1}}$]
$K_{\mathrm{P},\mathrm{dia}}$	Half-saturation constant for phosphorus of diatoms	0.0155	[mgP L${}^{-\mathrm{1}}$]
$K_{\mathrm{P},\mathrm{gra}}$	Half-saturation constant for phosphorus of green algae	0.062	[mgP L${}^{-\mathrm{1}}$]
$K_{\mathrm{P},\mathrm{cya}}$	Half-saturation constant for phosphorus of cyanobacteria	0.062	[mgP L${}^{-\mathrm{1}}$]
$K_{\mathrm{Si},\mathrm{dia}}$	Half-saturation constant for silica of diatoms	0.196	[mgSi L${}^{-\mathrm{1}}$]
$R_{\mathrm{m}\mathrm{20},\mathrm{dia}}^{*}$	Maintenance respiration coefficient of diatoms at 20 ${}^{\circ }$C	0.002	[h${}^{-\mathrm{1}}$]
$R_{\mathrm{m}\mathrm{20},\mathrm{gra}}^{*}$	Maintenance respiration coefficient of green algae at 20 ${}^{\circ }$C	0.002	[h${}^{-\mathrm{1}}$]
$R_{\mathrm{m}\mathrm{20},\mathrm{cya}}^{*}$	Maintenance respiration coefficient of cyanobacteria at 20 ${}^{\circ }$C	0.002	[h${}^{-\mathrm{1}}$]
$R_{\mathit{\mu },\mathrm{dia}}$	Energetic cost of growth of diatoms	0.5	[–]
$R_{\mathit{\mu },\mathrm{gra}}$	Energetic cost of growth of green algae	0.5	[–]
$R_{\mathit{\mu },\mathrm{cya}}$	Energetic cost of growth of cyanobacteria	0.5	[–]
$E_{\mathrm{cst},\mathrm{dia}}$	Basic excretion rate of diatoms	0.006	[h${}^{-\mathrm{1}}$]
$E_{\mathrm{cst},\mathrm{gra}}$	Basic excretion rate of green algae	0.006	[h${}^{-\mathrm{1}}$]
$E_{\mathrm{cst},\mathrm{cya}}$	Basic excretion rate of cyanobacteria	0.006	[h${}^{-\mathrm{1}}$]
$E_{\mathrm{phot},\mathrm{dia}}$	Excretion constant of diatoms related to photosynthesis	0.001	[–]
$E_{\mathrm{phot},\mathrm{gra}}$	Excretion constant of green algae related to photosynthesis	0.001	[–]
$E_{\mathrm{phot},\mathrm{cya}}$	Excretion constant of cyanobacteria related to photosynthesis	0.001	[–]
$S_{\mathrm{R},\max \mathrm{20},\mathrm{dia}}^{*}$	Maximal rate of reserve product synthesis for diatoms at 20 ${}^{\circ }$C	0.15	[h${}^{-\mathrm{1}}$]
$S_{\mathrm{R},\max \mathrm{20},\mathrm{gra}}^{*}$	Maximal rate of reserve product synthesis for green algae at 20 ${}^{\circ }$C	0.2	[h${}^{-\mathrm{1}}$]
$S_{\mathrm{R},\max \mathrm{20},\mathrm{cya}}^{*}$	Maximal rate of reserve product synthesis for cyanobacteria at 20 ${}^{\circ }$C	0.075	[h${}^{-\mathrm{1}}$]
$C_{\mathrm{R},\max \mathrm{20},\mathrm{dia}}^{*}$	Maximal rate of reserve product catabolism for diatoms at 20 ${}^{\circ }$C	0.2	[h${}^{-\mathrm{1}}$]
$C_{\mathrm{R},\max \mathrm{20},\mathrm{gra}}^{*}$	Maximal rate of reserve product catabolism for green algae at 20 ${}^{\circ }$C	0.2	[h${}^{-\mathrm{1}}$]
$C_{\mathrm{R},\max \mathrm{20},\mathrm{cya}}^{*}$	Maximal rate of reserve product catabolism for cyanobacteria at 20 ${}^{\circ }$C	0.2	[h${}^{-\mathrm{1}}$]
$k_{\mathrm{d}\mathrm{20},\mathrm{dia}}^{*}$	Rate of diatom mortality at 20 ${}^{\circ }$C	0.025	[h${}^{-\mathrm{1}}$]
$k_{\mathrm{d}\mathrm{20},\mathrm{gra}}^{*}$	Rate of green algae mortality at 20 ${}^{\circ }$C	0.025	[h${}^{-\mathrm{1}}$]
$k_{\mathrm{d}\mathrm{20},\mathrm{cya}}^{*}$	Rate of cyanobacteria mortality at 20 ${}^{\circ }$C	0.015	[h${}^{-\mathrm{1}}$]
${\mathrm{vs}}_{\mathrm{dia}}$	Sinking velocity of diatoms	0.006	[m h${}^{-\mathrm{1}}$]
${\mathrm{vs}}_{\mathrm{gra}}$	Sinking velocity of green algae	0.001	[m h${}^{-\mathrm{1}}$]
${\mathrm{vs}}_{\mathrm{cya}}$	Sinking velocity of cyanobacteria	0.006	[m h${}^{-\mathrm{1}}$]
$T_{\mathrm{opt},\mathrm{dia}}$	Optimal temperature of diatoms	21	[${}^{\circ }$C]
$T_{\mathrm{opt},\mathrm{gra}}$	Optimal temperature of green algae	37	[${}^{\circ }$C]
$T_{\mathrm{opt},\mathrm{cya}}$	Optimal temperature of cyanobacteria	37	[${}^{\circ }$C]
${\mathit{\sigma }}_{\mathrm{dia}}$	Range of temperature for diatoms	13	[${}^{\circ }$C]
${\mathit{\sigma }}_{\mathrm{gra}}$	Range of temperature for green algae	15	[${}^{\circ }$C]
${\mathit{\sigma }}_{\mathrm{cya}}$	Range of temperature for cyanobacteria	12	[${}^{\circ }$C]

${}^{*}$ Parameters depend on water temperature and are multiplied by $f\left(T\right)=\frac{e^{-\frac{(T-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}{e^{-\frac{(\mathrm{20}-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}$, where $T$ is water temperature in ${}^{\circ }$C.

Table B4

Parameters for organic matter dynamics.

Parameter	Description	Value	Unit
${\mathit{ϵ}}_{{\mathrm{dom}}_{\mathrm{1}}}$	DOM${}_{\mathrm{1}}$ fraction in lysis products	0.2	[–]
${\mathit{ϵ}}_{{\mathrm{dom}}_{\mathrm{2}}}$	DOM${}_{\mathrm{2}}$ fraction in lysis products	0.2	[–]
${\mathit{ϵ}}_{{\mathrm{dom}}_{\mathrm{3}}}$	DOM${}_{\mathrm{3}}$ fraction in lysis products	0.1	[–]
${\mathit{ϵ}}_{{\mathrm{pom}}_{\mathrm{1}}}$	POM${}_{\mathrm{1}}$ fraction in lysis products	0.2	[–]
${\mathit{ϵ}}_{{\mathrm{pom}}_{\mathrm{2}}}$	POM${}_{\mathrm{2}}$ fraction in lysis products	0.2	[–]
${\mathit{ϵ}}_{{\mathrm{pom}}_{\mathrm{3}}}$	POM${}_{\mathrm{3}}$ fraction in lysis products	0.1	[–]
$k_{{\mathrm{pom}}_{\mathrm{1}},\mathrm{20}}^{*}$	POM${}_{\mathrm{1}}$ hydrolysis rate constant at 20 ${}^{\circ }$C	0.005	[h${}^{-\mathrm{1}}$]
$k_{{\mathrm{pom}}_{\mathrm{2}}}$	POM${}_{\mathrm{2}}$ hydrolysis rate constant	0.00025	[h${}^{-\mathrm{1}}$]

${}^{*}$ Parameters depend on water temperature and are multiplied by $f\left(T\right)=\frac{e^{-\frac{(T-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}{e^{-\frac{(\mathrm{20}-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}$, where $T$ is water temperature in ${}^{\circ }$C.

Table B5

Zooplankton parameters.

Parameter	Description	Value	Unit
${\mathit{\mu }}_{\max \mathrm{20},\mathrm{zor}}^{*}$	Maximal growth rate of ZOR at 20 ${}^{\circ }$C	0.025	[h${}^{-\mathrm{1}}$]
${\mathit{\mu }}_{\max \mathrm{20},\mathrm{zoc}}^{*}$	Maximal growth rate of ZOC at 20 ${}^{\circ }$C	0.015	[h${}^{-\mathrm{1}}$]
$b_{\max \mathrm{20},\mathrm{zor}}^{*}$	Maximal grazing rate of ZOR at 20 ${}^{\circ }$C	0.1	[h${}^{-\mathrm{1}}$]
$b_{\max \mathrm{20},\mathrm{zoc}}^{*}$	Maximal grazing rate of ZOC at 20 ${}^{\circ }$C	0.05	[h${}^{-\mathrm{1}}$]
$K_{\mathrm{phy},\mathrm{zor}}$	Half-saturation constant for grazing phytoplankton of ZOR	0.1	[mgC L${}^{-\mathrm{1}}$]
$K_{\mathrm{phy},\mathrm{zoc}}$	Half-saturation constant for grazing phytoplankton of ZOC	0.1	[mgC L${}^{-\mathrm{1}}$]
PHY${}_{\mathrm{0},\mathrm{zor}}$	Threshold phytoplankton concentration for grazing of ZOR	0.1	[mgC L${}^{-\mathrm{1}}$]
PHY${}_{\mathrm{0},\mathrm{zoc}}$	Threshold phytoplankton concentration for grazing of ZOC	0.1	[mgC L${}^{-\mathrm{1}}$]
${}^{\ast }{k}_{\mathrm{d}\mathrm{20},\mathrm{zor}}$	Mortality rate of ZOR at 20 ${}^{\circ }$C	0.007	[h${}^{-\mathrm{1}}$]
${}^{\ast }{k}_{\mathrm{d}\mathrm{20},\mathrm{zoc}}$	Mortality rate of ZOC at 20 ${}^{\circ }$C	0.007	[h${}^{-\mathrm{1}}$]
$T_{\mathrm{opt},\mathrm{zor}}$	Optimal temperature of ZOR	25	[${}^{\circ }$C]
$T_{\mathrm{opt},\mathrm{zoc}}$	Optimal temperature of ZOC	25	[${}^{\circ }$C]
${\mathit{\sigma }}_{\mathrm{zor}}$	Range of temperature for ZOR	10	[${}^{\circ }$C]
${\mathit{\sigma }}_{\mathrm{zoc}}$	Range of temperature for ZOC	10	[${}^{\circ }$C]
${\mathrm{vs}}_{\mathrm{zor}}$	Sinking velocity of ZOR	0.02	[m h${}^{-\mathrm{1}}$]
${\mathrm{vs}}_{\mathrm{zoc}}$	Sinking velocity of ZOC	0.02	[m h${}^{-\mathrm{1}}$]

${}^{*}$ Parameters depend on water temperature and are multiplied by $f\left(T\right)=\frac{e^{-\frac{(T-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}{e^{-\frac{(\mathrm{20}-T_{\mathrm{opt}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}}$, where $T$ is water temperature in ${}^{\circ }$C.

Table B6

Phosphate- and silica-related parameters.

Parameter	Description	Value	Unit
Phosphate adsorption desorption
$P_{\mathrm{ac}}$	Maximal adsorption capacity of mineral suspended solids (MSSs)	0.00558	[mgP mgMSS${}^{-\mathrm{1}}$]
$K_{\mathrm{ps}}$	Half-saturation adsorption constant	0.682	[mgP L${}^{-\mathrm{1}}$]
Silica dynamics
Kb${}_{\mathrm{Si}\mathrm{20}}$	Biogenic silica dissolution rate at 20 ${}^{\circ }$C	0.0001	[h${}^{-\mathrm{1}}$]

Code and data availability

C-RIVE 0.32 implements the unified RIVE v1.0 in ANSI C. It is available under Eclipse Public License 2.0 in the following Zenodo repository: 10.5281/zenodo.7849609 . pyRIVE 1.0 implements the unified RIVE v1.0 in Python 3 and is available under Eclipse Public License 2.0 in the InDoRES repository: 10.48579/PRO/Z9ACP1 . The dataset used in this paper is available in the following Zenodo repository: 10.5281/zenodo.10490669 ().

Author contributions

SW: conceptualization, methodology, software, validation, formal analysis, investigation, visualization, writing (original draft). VT: conceptualization, methodology, software, formal analysis, writing (review and editing), supervision. GB: conceptualization, methodology, data curation, formal analysis, writing (review and editing). JG: conceptualization, data curation, formal analysis, writing (review and editing). MS: software, writing (review and editing). AM: conceptualization, software. XY: software. NF: conceptualization, methodology, formal analysis, software, writing (review and editing), supervision, funding acquisition.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

This work was carried out under the PIREN-Seine program (https://www.piren-seine.fr/, last access: 20 December 2023), part of French Long-Term Socio-Ecological Research (LTSER) (in French: “Zone Atelier Seine” – INEE, CNRS). We thank Lou Weidenfeld, Benjamin Mercier, and Anun Marninez for carrying out organic matter batch experiments.

Particularly, we would like to dedicate this paper to the memory of Pierre Servais: his knowledge, expertise, and passion for this field of study will continue to guide us in the pursuit of our scientific work.

Review statement

This paper was edited by Vassilios Vervatis and reviewed by two anonymous referees.