Evaluating the performance of CE-QUAL-W2 version 4.5 sediment diagenesis model

This research evaluates the performance of the CE-QUAL-W2 v4.5 sediment diagenesis model in simulating water temperature, dissolved oxygen, total phosphorus, total nitrogen, chlorophyll a, and biochemical oxygen demand in a Portuguese reservoir over a six-year period (2016–2021). The model was calibrated using 35 observed profiles of temperature and dissolved oxygen, as well as six annual measurements of total nitrogen, total phosphorus, chlorophyll a, and biochemical oxygen demand at multiple depths. To benchmark performance to the sediment diagenesis model, three alternative sediment oxygen demand formulations – a Zero-order, First-order, and a Hybrid model combining both approaches – were also implemented and compared. All models achieved NSE and RMSE values within or near the ranges reported in the literature, effectively capturing the system's water quality dynamics. Among them, the Hybrid model yielded the best overall performance while maintaining a simpler structure (Water temperature – NSE: 0.96±0.18; RMSE: 1.09±0.23 °C; Dissolved oxygen – NSE: 0.76±0.30; RMSE: 1.87±0.72 mg L⁻¹). The sediment diagenesis model exhibited similar performance metrics (Water temperature – NSE: 0.95±0.18; RMSE: 1.13±0.28 °C; Dissolved oxygen – NSE: 0.71±0.14; RMSE: 2.01±0.59 mg L⁻¹). Overall, the results suggest that the diagenesis model may be better suited for capturing detailed process-based dynamics over extended timeframes, whereas simpler models, such as the Hybrid model, are more appropriate for short- to medium-term applications or situations with limited data availability. Hopefully, the results of this study will help improve water management strategies by supporting more informed model selection tailored to the temporal scope and data constraints of reservoir monitoring programs.

1 Introduction

Modeling water quality plays a crucial role in managing lakes and reservoirs, providing essential insights into the dynamics of nutrients, organic matter, and phytoplankton within aquatic systems (Abbaspour et al., 2015; Whitehead et al., 2009). These models simulate the physical, chemical, and biological processes that influence water quality, with examples including widely-used tools like CE-QUAL-W2 (Wells, 2021), MIKE21 (Chapman, 1996), and DYRESM (Hamilton and Schladow, 1997). The value of such modeling lies in its capacity to aid researchers and policymakers in understanding the complex interactions between various factors that impact the ecological health of water bodies (Varis et al., 1994; Loucks and Beek, 2017). However, the intricacy of these systems, combined with the substantial data requirements, often presents significant challenges for those developing and applying water quality models. Effective inflow data characterization (quantity and quality) is hard to obtain, both for major river branches and small tributaries, as is waterbody sediment characterization related to carbon and nutrients due to the significant cost associated with the sampling and laboratorial analysis process and the fact that water management stakeholders are still more focused on the classification of waterbody water quality rather than the collection of water quality forcing data. The absence of sediment initial particulate organic carbon (POC), particulate organic nitrogen (PON) and particulate organic phosphorus (POP) data can be decisive to the overall performance of a water quality model, in essence generating an imbalance right from the start of the simulation with regard to the sediment concentration of POC, PON and POP, which then has a considerable impact on the SOD and, consequently, the waterbody dissolved oxygen (DO). When calibrating the model, water quality modelers therefore need to plug this gap by evaluating the model performance considering: (i) different initial sediment oxygen demand (SOD) where a zero-order model is applied, (ii) different POC, PON and POP values where a predictive diagenesis model is considered.

The main challenge with these modeling approaches is that the sources of DO depletion – such as the inflow of organic matter and net algae growth and settling – can significantly influence DO dynamics, and these sources must be well characterized to ensure accurate predictions. While the baseline model can reproduce observed DO profiles with reasonable accuracy, its predictive reliability may be compromised if key DO sinks and sources are not well defined.

For example, the model's response to a reduction in external phosphorus loading is influenced by internal phosphorus release from sediments during anoxic periods. In CE-QUAL-W2, when a zero-order SOD model is used, the anoxic release of orthophosphate (PO₄) is modeled as a linear function of SOD: SOD [g O₂ m⁻² d⁻¹] × PO₄ release rate [g P g⁻¹ O₂]. Thus, any error in the estimation of SOD will directly affect the predicted internal phosphorus loading, and by extension, the overall phosphorus balance in the waterbody. In contrast, when using the predictive sediment diagenesis model, internal phosphorus loading depends on the organic and nutrient inputs from particulate matter in the water column and the sediment's biogeochemical response, which is highly influenced by the initial value of particulate organic carbon (POC). As a result, this approach introduces additional uncertainty when key particulate components are not adequately measured or constrained in both the water column and sediments. Calibrating other constituents, such as P-PO₄, can help reduce uncertainty. P-PO₄ is released from sediments under anaerobic conditions, and its calibration can enhance the accuracy of DO modeling. Still, this release is influenced by multiple factors, including the initial sediment P-PO₄ concentration and the release rate (in the zero-order model), or the mineralization of POP (in the diagenesis model). In both cases, significant uncertainty remains without observed data for POC, PON, and POP in both the water column and sediments. Of these, POC has the most significant influence on SOD, making access to sediment POC data essential for improving model accuracy, even when PON and POP measurements are lacking. The CE-QUAL-W2 model has been widely used to simulate various water bodies and water quality scenarios, including reservoir physical and biochemical dynamics in response to warming projections (Mi et al., 2020, 2023). This model has also been used to predict DO in a number of water bodies worldwide, although the SOD has always been modeled with a zero-order and/or 1-order model (e.g. Park et al., 2014; Zouabi-Aloui et al., 2015; Terry et al., 2017; Sadeghian et al., 2018; Lindenschmidt et al., 2019). The bibliographic research conducted before and during this study suggests that the CE-QUAL-W2 sediment diagenesis model has not been applied to any waterbodies other than the Wahiawa Reservoir in central Oahu (Berger and Wells, 2014). Moreover, no scientific publications on the evaluation of this model in other contexts have been identified, further highlighting the importance of the primary motivation for this study, namely, to evaluate the performance of the CE-QUAL-W2 model with its new sediment diagenesis component. This study benefited from having access to observed reservoir sediment total organic carbon (TOC) values, which are rare. Although, in theory, these values are typically higher than particulate organic carbon (POC) values, they provided an excellent starting point for this study. The methodological approach was, therefore, defined to evaluate the performance of the CE-QUAL-W2 model considering the new state-of the art sediment diagenesis model in modeling a reservoir, DO, Total Phosphorus (Total P); Total Nitrogen (Total N), Biochemical Oxygen Demand (BOD₅), Chlorophyll a (Chl a) and SOD.

To achieve this, the water quality of a highly productive reservoir was simulated over a six-year period (2016–2021) using the CE-QUAL-W2 v4.5 model. The simulation incorporated a Zero-order sediment model, a First-order model, a Hybrid model combining both approaches, and a sediment diagenesis model. The Zero-order, First-order, and Hybrid models were included to provide alternative representations of sediment oxygen demand, enabling comparative analysis and supporting the calibration and evaluation of the more complex sediment diagenesis model. In the case of water temperature and DO, the modeling results were compared with 35 water column profiles observed near the dam. The remaining parameters were calibrated using time series datasets collected at multiple depths, with six annual values available for each parameter. A sensitivity analysis was performed to evaluate the reservoir water quality response, namely DO, to the variation of POC, PON and POP concentration in the reservoir sediments. The results of this study will hopefully prove useful by helping to improve lake and reservoir water quality modeling and, therefore, the water management process from a practical perspective.

2 Methods

2.1 Site location and main characteristics

Portugal experiences a temperate maritime climate characterized by a wet, cool season and a dry summer. Despite most of the precipitation occurring during the winter months, there is significant inter-annual variability. Precipitation patterns are spatially and temporally heterogeneous, with annual maxima exceeding 2500 mm in the rugged highlands of the northwest, while the low-lying plains of the southeast receive around 400 mm yr⁻¹ (Cardoso et al., 2013; Soares et al., 2015) (Fig. 1). The Torrão dam, located in the northern region of mainland Portugal in the Tâmega River, is a significant hydraulic structure designed for multiple purposes, including water supply, irrigation, and hydroelectric power generation. The reservoir has a substantial storage capacity, contributing to regional water management and flood control (Table 1). This infrastructure plays a crucial role in the socio-economic development of the region, balancing resource management and environmental preservation. However, it is also important to note that the reservoir was classified as eutrophic for all the simulated years, a condition that can lead to persistent water quality issues.

Figure 1 Torrão reservoir watershed. Thiessen polygons. Water quality stations.

Table 1 Main features of Torrão dam and reservoir.

^* Classification according to OECD Trophic State limits (OECD, 1982).

2.2 CE-QUAL-W2 v4.5 model

This study employed the latest version of CE-QUAL-W2 (Version 4.5), a model originally developed in 1975 by the US Army Corps of Engineers and written in Fortran. Since its inception, the model has undergone regular updates and enhancements, primarily by researchers at Portland State University (Cole and Wells, 2006). CE-QUAL-W2 is a two-dimensional, laterally averaged hydrodynamic and water quality model capable of simulating free surface elevation, hydrostatic pressure, density, horizontal and vertical velocities, as well as constituent concentrations. The model uses the finite difference method to solve key equations, including mean transverse momentum in the x- and z-directions, the continuity equation, state equations, and water surface elevation equations (Tavera-Quiroz et al., 2023; Wells, 2021). This model represents SOD through four distinct approaches: (i) a user-defined zero-order formulation that is decoupled from the water column, (ii) a simple predictive first-order model, (iii) a hybrid approach combining the zero- and first-order methods, and (iv) a comprehensive sediment diagenesis model. The zero-order model is not a predictive approach, as, other than variations resulting from the temperature dependence of the decay rate, the rates remain constant over time (Wells, 2021). Additionally, under anoxic conditions in the water column, SOD is disabled in the model. The first-order sediment model does not function as a full sediment diagenesis model, as it lacks the capability to track the fate of organic nutrients delivered to the sediments, their breakdown, and the release of byproducts into the water column under low-oxygen conditions. However, it does represent the deposition of particulate organic matter and dead algal biomass, along with the resulting oxygen demand imposed on the water column. By including this first-order sediment process, the model becomes sensitive to increased organic loading to the sediment, which in turn influences sediment oxygen demand. A combination of the zero and first order model can be considered where organic materials accumulate and decay in the sediments under aerobic conditions and are released based on the SOD zero-order decay rate under anaerobic conditions. In contrast, the sediment diagenesis model simulates kinetic processes occurring within the sediment and at the sediment–water interface. This module originally developed for application in oil sand pit lakes, has been adapted for application in other aquatic environments and integrated into version 4.0 (Vandenberg et al., 2015). The conceptual framework of the model has been elaborated in works by Prakash et al., (2015), Berger and Wells (2014), and Vandenberg et al. (2015). It is important to note that significant enhancements to the sediment diagenesis module were introduced in version 4.5 of the model, as detailed in the User Manual (Wells, 2021). These improvements mark a substantial advancement over the initial version 4, which was more limited in its capabilities. The CE-QUAL-W2 model has demonstrated its utility in simulating hydrodynamic and ecological processes – such as stratification, internal waves, DO dynamics, and phytoplankton blooms – in lakes and reservoirs worldwide (Zhang et al., 2015; Chuo et al., 2019; Kobler et al., 2018; Uhlmann, 2017; Terry et al., 2017; Mi et al., 2020). Additional details about the model's structure, algorithms, and historical applications can be found in the user manual (Wells, 2021).

Model setup

The bathymetry of the Torrão reservoir was initially defined using a Digital Elevation Model (DEM) provided by Energies of Portugal, SA (EDP) and structured according to the methodology outlined in Wells (2021). The reservoir comprises one main branch (the Tâmega River), three tributaries and one distributed tributary (Fig. 1). Tributaries 1 and 2 are depicted in Fig. 1. Tributary 3 represents the inflow from the Douro River into the pump-back system of the Torrão Reservoir. The bathymetric map includes 27 segments, each measuring 1000 m in length, and a maximum number of 58 layers, each with a depth of 1 m. Following this preliminary step, the reservoir boundary conditions (including water quality, hydrology, meteorology, and sediment characterization) were defined according to the methods described in Sect. 2.4. Due to the lack of available information, the model structure only includes a single algae group (Diatoms).

2.3 Modeling approach

To thoroughly evaluate the capability of CE-QUAL-W2 in modeling dissolved oxygen using the sediment diagenesis module, the four available SOD modeling approaches were considered: Zero-order model; First-order model; Zero/First-order model (Hybrid model) and the sediment diagenesis model (SG model). The models were calibrated for the 2016–2021 period (see Sect. 2.5). During the results analysis, the performance metrics obtained during each model's calibration process were compared, along with the SOD values across the bottom layers of each model. A sensitivity analysis was conducted following calibration to evaluate each model's response: (a) to varying POC, PON, and POP values in the case of the SG model; (b) to different SOD values in the Zero-order and Hybrid models and (c) to varying the initial first order sediment concentration in the case of the First-order model. Section 2.6 details the methodological approach used for the sensitivity analysis. To assess the sensitivity of each model to reductions in external organic matter (OM) and P-PO₄ inputs, two separate scenario analyses were conducted. The first scenario involved an 80 % reduction in OM inflow load, while the second applied an 80 % reduction in both OM and P-PO₄ inflow loads. These reductions were implemented specifically in the main reservoir branch (Branch 1 – Tâmega River), where the majority of nutrient and organic inputs occur. Each sediment model – SD, Zero-order, First-order, and Hybrid – was run under baseline conditions and under both reduction scenarios. The impact on DO dynamics was evaluated using time series of depth- and segment-averaged DO concentrations. The evaluation of model performance, along with the results of the sensitivity analysis, provided deeper insights into simulating SOD dynamics using the sediment diagenesis approach in comparison to the other SOD formulations.

2.4 Model forcing datasets

The meteorological data used to drive the model, including hourly air temperature, dew point, solar radiation, cloud cover, and wind characteristics, were sourced from ERA5-Land, a high-resolution reanalysis dataset optimized for land applications (Muñoz-Sabater, 2019). Although no on-site meteorological stations are available in the study area for direct validation, studies by Almeida and Coelho (2023b) and Barbosa and Scotto (2022) have demonstrated a strong correlation between ERA5-Land air temperature data and observed measurements at regional scales, supporting the reliability of this dataset for our modeling purposes. Furthermore, the accuracy of water temperature predictions in our simulations indicates that the meteorological forcing was well represented, confirming the suitability of ERA5-Land data for driving the model. Reservoir data, such as daily inflow/outflow, water levels, and water quality, covering the years 2016–2021, were provided by EDP. Water quality data specific to Branch 1 originated from the Praia Aurora Station, accessed via the Portuguese National Water Resources Information System (SNIRH, 2024). With only 21 recorded measurements for Branch 1 during this period, three modeling methods were employed to address the 99.04 % of missing data. The variables include: water temperature; DO; Total P; Ammonium (N-NH₄); Nitrate+Nitrite (N-NO_x); BOD₅; Chl a; Alkalinity; Conductivity and Total Suspended Solids (SST). The first method employed regression models implemented through the LOADEST package (Runkel et al., 2004) developed by the US Geological Survey. The second method utilized the Extreme Gradient Boosting (XGBoost) machine learning algorithm, implemented using the Chen and Guestrin (2016) open-source library, a method proven effective in various environmental studies (Feigl et al., 2021; Adedeji et al., 2022; Xu et al., 2022). For additional details on the algorithm, refer to Almeida and Coelho (2023a). The third approach relied on Support Vector Regression (SVR), implemented via the scikit-learn library (Pedregosa et al., 2011), which has also demonstrated strong performance in environmental modeling applications (Adedeji et al., 2022; Ji and Lu, 2018). For machine learning approaches, datasets were split into training (80 %) and testing (20 %) sets. Hyperparameters for these models were optimized using the Tree-structured Parzen Estimators (TPE) algorithm, executed with the Hyperopt library (Bergstra et al., 2013) and 100 iterations. The Nash–Sutcliffe Efficiency (NSE) was used to determine the best model. Table A1 describes the input features of each model. Correlations derived from Branch 1 informed data extrapolation to other tributaries using flow as the predictor. Observed data for Tributary 3 was retrieved from the Crestuma–Lever reservoir monitoring station.

Water quality variables used for model inputs included water temperature, DO, P-PO₄, N-NH₄, N-NO_x, labile and refractory dissolved and particulate organic matter (LDOM, RDOM, LPOM, RPOM), alkalinity, inorganic suspended solids (ISS), total dissolved solids (TDS), total inorganic carbon (TIC), and algal biomass (diatoms). For non-monitored variables, estimations were made based on available data: (i) P-PO₄: Derived from total phosphorus, assuming inorganic phosphorus represents 70 % of the total; (ii) Organic matter: BOD₅ was converted to organic matter using a stoichiometric ratio of 1.4 g O₂ per 1.0 g organic matter, with 60 % assumed refractory and 40 % labile; (iii) ISS: Estimated as 97.4 % of TSS; (iv) TDS: Calculated from electrical conductivity (Eq. 1); (v) TIC: Estimated from alkalinity (Eq. 2); (vi) Algae biomass: Chl a was converted to biomass using the following ratio: Algal Biomass (mg L⁻¹)/Chl a(µg L⁻¹) = 5.0 $\times {\mathrm{10}}^{-\mathrm{2}}$

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{TDS}\left(\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}\right)=\mathrm{0.65}\times \mathrm{Electrical}\mathrm{Conductivity}\left(\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}\right)\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{TIC}\left(\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}\right)=\mathrm{0.2782}\times \mathrm{Alkalinity}\left(\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}\right)\mathrm{0.9706}.\end{array}$$

This equation was derived from the relationship between TIC and alkalinity values observed in four reservoirs within the United States, utilizing a dataset comprising 55232 value pairs available in the CE-QUAL-W2 v4.5 model examples (Wells, 2021). The analysis achieved an R² value of 0.99. Figure 2 illustrates the locations of the five sediment sampling sites used to define the SD model baseline run. The spatial distribution of the sediment samples depicted in the figure were linked to specific reservoir segments to characterize the initial sediment content of POC, PON, and POP, as detailed in Table 2. Sediment values were assigned as follows: site A to segments 25–28, site B to segments 20–24, site C to segments 16–19, site D to segments 11–15, and site E to segments 2–10. Several assumptions were made to establish the sediment characterization: (i) a sediment density of 960 kg m⁻³(density of dried sediment with air in the pore space) (Minear and Kondolf, 2009) was applied to convert sample values from mg kg⁻¹ to mg L⁻¹; (ii) POP values were set at 25 mg L⁻¹, based on established literature benchmarks (Wells, 2021); (iii) the TN value observed at site B was used to characterize sites C, D, and E; (iv) TOC and TN were assumed to exist entirely in particulate form, represented as particulate organic carbon and nitrogen. This approach ensured a consistent and representative characterization of sediment properties across the reservoir segments.

Figure 2 Sediment sampling sites. CE-QUAL-W2 model segments.

Table 2 Torrão Reservoir sediment chemical characterization obtained for each sampling site.

2.5 Water quality model (CE-QUAL-W2) calibration

The simulation period considered for this study spanned 2016 to 2021. This period was selected due to the availability of flow and water quality data. The trial-and-error technique was applied to calibrate the model for the simulation period, considering the default calibration parameters described in Wells, 2021. The error between observed and predicted values of six state variables was evaluated with five different metrics (see Sect. 2.7). The observed data included six annual values for water temperature, DO, TP, TN, BOD₅, and Chl a. These time series were obtained from: (a) an integrated sample between the reservoir surface and a depth of 5.8 m, (b) a depth of 23 m, and (c) a depth of 43.7 m. In addition, 35 water temperature and DO profiles – six per year from 2016 to 2021 – were also included. These profiles were observed 300 meters upstream from the Torrão Dam. Details on the models' initial conditions, parameters, constants, and forcing datasets can be found in Almeida and Coelho (2025) and in Tables A2 to A8. The models were calibrated by adjusting their parameters to improve the fit between the model output and observed data. Please refer to Wells (2021) for a detailed account of the model calibration parameters and default values. Water temperature was the first constituent to be calibrated. The wind sheltering coefficient (WSC) was manually adjusted to achieve the best fit between the modeled and observed water temperature profiles, resulting in a final value of 1. A value of 1 implies that the WSC has no effect over the wind velocity forcing the model. The zero-order model for SOD was then manually adjusted to improve DO predictions based on 35 DO profiles. The optimal result was achieved with a zero-order SOD value of 2.5 g O₂ m⁻² d⁻¹. Following this calibration, the phosphorus sediment release rate (PO₄R) in the zero-order model was modified from its default value of 0.015 to 0.001. The same process was applied to the Hybrid model, where the best results were achieved using a PO₄R of 0.001 and a zero-order SOD value of 1.0 g O₂ m⁻² d⁻¹. In the first-order model, the PO₄R parameter was adjusted to 0.001, and the initial concentration of first-order sediment was set to 0.5 g m⁻². All other parameters were kept at their relevant default values and the default settings for the sediment diagenesis model were also maintained. The observed data included water temperature, DO, TP, TN, BOD₅, and Chl a.

2.6 Sensitivity analysis

A sensitivity analysis was conducted after the calibration process to evaluate the model's response:

• i.

  Different initial sediment values for POC, PON, and POP were used in the SD model (Table 3). It is important to note that for each of the 24 runs described in Table 2, only the corresponding parameter was modified, while the other two parameters retained their default values shown in Table 3. The number of runs varying the PON and POP values is higher than the number of runs considered for POC, with 6 versus 9 runs, respectively. This adjustment was necessary to achieve a minimal RMSE in the predictions of dissolved oxygen in the water column;

• ii.

  Different zero-order SOD values for the Zero-order model (0.5, 1.0, 1.5, 2.0, 2.5 and 3.0 g O₂ m⁻² d⁻¹);

• iii.

  Different initial first order sediment concentration (ISC) for the First-order model (0.0, 0.5, 1.0, 1.5, 2.0, 2.5 and 3.0 g m⁻²);

• iv.

  Different zero-order SOD values for the Hybrid model (0.5, 1.0, 1.5, 2.0, 2.5 and 3.0 g O₂ m⁻² d⁻¹.

In the results analysis for each run and for both scenarios (i) and (ii), the prediction error for DO was compared with the SOD values derived from each model. Specifically, runs 5, 8, and 20 were forced with the POC, PON, and POP values defined in the SD model baseline run.

Table 3 TOC, PON, and POP initial sediment values for the SD model sensitivity analyses.

2.7 Metrics

The evaluation of model calibration and the analysis of quantitative differences across simulation scenarios utilized various performance metrics. These included the root mean square error (RMSE), mean absolute error (MAE), Nash–Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970), percent bias (PBIAS), and the coefficient of determination (R²). The calculations were carried out using equations where m_i and o_i represent the simulated and observed values, respectively, and $\stackrel{\mathrm{‾}}{o_i}$ the observed values mean.

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\left(m_i-o_i\right)}^{\mathrm{2}}}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{MAE}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left|m_i-o_i\right|\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{NSE}=\mathrm{1}-\left[{\displaystyle \frac{\sum _{i=\mathrm{1}}^N{\left(o_i-m_i\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^N{\left(o_i-\stackrel{\mathrm{‾}}{o_i}\right)}^{\mathrm{2}}}}\right]\end{array}$$

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }\mathrm{PBIAS}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N\left(o_i-m_i\right)}{\sum _{i=\mathrm{1}}^N\left(o_i\right)}}\times \mathrm{100}\text{(7)}& {\displaystyle }& {\displaystyle }{R}^{\mathrm{2}}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N{\left(m_i-\stackrel{\mathrm{‾}}{o}\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^N{\left(o_i-\stackrel{\mathrm{‾}}{o}\right)}^{\mathrm{2}}}}\times \mathrm{100}\end{array}$$

3 Results

3.1 Observed inflow water quality characterization

The SVR algorithm was more effective at predicting the inflow water temperature compared to the other models. The R² and PBIAS values achieved with the SVR were 0.87, and 3.77 %, respectively, indicating that the water temperature trends and average magnitudes are well described (Table A1). Additionally, the RMSE and MAE values of 2.1 and 1.6 °C, respectively, demonstrate an accurate approximation of the observed datasets. The SVR algorithm was also the best model in predicting DO. The R², PBIAS, RMSE, and MAE values reached, 0.91, 0.92 %, 0.40 and 0.26 mg L⁻¹, respectively, indicating that the model performed well. This was not the case for the remaining parameters. In fact, the Loadest regression outperformed the other models for the remaining water quality variables. This was primarily due to the limited number of training samples. Simpler models like regressions can have lower variance (i.e., be less susceptible to overfitting) compared to SVR and XBOOST algorithms. Overall, the PBIAS obtained for NH₄, N-NO_x, and Chl a (10.88 %, 43.64 %, and 30.00 %) suggests that the average magnitude was reasonably well represented.

Figure 3 Constituents observed values at three different depths: (a) an integrated sample between the reservoir surface and an average depth of 5.8 m, (b) an average depth of 23 m, and (c) an average depth of 43.7 m. These observed values were compared with the predicted time series from the SD model (run 5 – baseline) (A to F) for the same depths.

3.2 CE-QUAL-W2 calibration

Tables A2 through A8 display the most significant CE-QUAL-W2 coefficients obtained after the calibration process. The results of the calibration process for all models, are presented in Tables 4 and A9 and illustrated in Figs. 3 to 6 and Figs. 8 and 9. The performance metrics for water temperature across the different sediment models show consistent accuracy, with NSE and R² values ranging from 0.95 to 0.96 and minimal variation across models. The RMSE and MAE for temperature also remain low, indicating reliable thermal performance regardless of the sediment model applied. In contrast, DO predictions show more variability. The Hybrid model achieved the best overall DO performance, with the highest NSE (0.76±0.30) and R² (0.76±0.31), as well as the lowest RMSE (1.87±0.72) and MAE (1.22±0.55), while maintaining a near-zero PBIAS ($-\mathrm{0.55}\pm \mathrm{11.14}$), indicating minimal systemic bias. The Zero-order model also performed reasonably well, with slightly lower error metrics than the SD model. The First-order model, however, showed the weakest DO performance, with a lower NSE (0.68±0.22), higher RMSE (2.15±0.82), and a significant negative PBIAS ($-\mathrm{12.17}\pm \mathrm{15.44}$), suggesting an underestimation of oxygen concentrations. Overall, the results suggest that while temperature simulation is robust across all models, DO dynamics are better captured using the Hybrid or Zero-order models, with the Hybrid model offering the most balanced and accurate representation under the tested conditions. However, the differences in performance metrics for DO among the models are relatively small and often fall within overlapping standard deviations, with the exception of the First-order model, which consistently shows lower accuracy and higher bias, suggesting that while the Hybrid model offers slightly better overall performance, the improvements over the SD and Zero-order models are modest and should be interpreted with caution. In terms of nutrient dynamics, the Hybrid and Zero-order models improve TN and TP predictions relative to the SD and First-order models. The Hybrid model, for example, improves TN R² to 0.31 and TP to 0.27, although the associated biases remain significant (e.g., −18.75 % for TN and +36.49 % for TP). BOD5 and Chl a remain poorly simulated across all models, with R² values consistently low (≤0.06 for Chl a and ≤0.03 for BOD5), and large PBIAS values, particularly in the SD and First-order configurations. The Zero-order model slightly reduces bias in Chl a and Total N compared to the SD model but performs poorly for TP due to a large overestimation (PBIAS = 103.43 %) (Fig. 4d). Notably, the SD and First-order models failed to reproduce observed phosphorus release events from sediments on 18 September 2018, 8 September 2020, and 31 August 2021 (Figs. 3d and 5d). In contrast, the Hybrid model successfully captured these events by modeling phosphorus release as a linear function of SOD, providing a more realistic representation of sediment–water nutrient interactions (Fig. 6d). Overall, while no model fully captures the complexity of all constituents, the Hybrid model consistently provides the most balanced and improved representation, particularly for DO and nutrient parameters.

Table 4 Metrics between observed and predicted values for all models. Water temperature and DO metrics were obtained from 36 observed and predicted profiles.

Figure 4 Constituents observed values at three different depths: (a) an integrated sample between the reservoir surface and an average depth of 5.8 m, (b) an average depth of 23 m, and (c) an average depth of 43.7 m. These observed values were compared with the predicted time series from the Zero-order model (zero order SOD = 2.5 g O₂ m⁻² d⁻¹ – baseline) (A to F) for the same depths.

Figure 5 Constituents observed values at three different depths: (a) an integrated sample between the reservoir surface and an average depth of 5.8 m, (b) an average depth of 23 m, and (c) an average depth of 43.7 m. These observed values were compared with the predicted time series from the First-order model (ISC = 0.5 g m⁻² – baseline) (A to F) for the same depths.

Figure 6 Constituents observed values at three different depths: (a) an integrated sample between the reservoir surface and an average depth of 5.8 m, (b) an average depth of 23 m, and (c) an average depth of 43.7 m. These observed values were compared with the predicted time series from the Hybrid model (zero order SOD = 1.0 g O₂ m⁻² d⁻¹ – baseline) (A to F) for the same depths.

3.3 Sensitivity analysis

The SOD values strongly influence the water column DO; therefore, this parameter was considered to support this analysis. Figure 7 shows the SOD values from the reservoir bottom layer, predicted by the SD model for Runs 1 to 6, compared with the RMSE (Fig. 7a) and the NSE (Fig. 7b) values obtained between the predicted water column DO profiles and the mean initial POC values (across all sites values) for each run. These results suggest that Run 4 was the best modeling solution. Considering the results obtained for Run 5 (baseline), Run 4 reduced the RMSE from 2.015 mg L⁻¹ (Run 5) to 2.011 mg L⁻¹ (Run 4) and increased the NSE from 0.714 (Run 5) to 0.716 (Run 4). The average SOD value in the bottom layer of the reservoir (across all model segments) decreased from 1.162ġ O₂ m⁻³ d⁻¹ (Run 5) to 1.071 g O₂ m⁻³ d⁻¹ (Run 4). Although the reduction is modest and had only a minor effect on the DO profile predictions (Fig. 9), it suggests that the initial POC values used in Run 5 were likely overestimated. This outcome aligns with the assumption made in Run 5, where all observed TOC was considered to exist entirely as POC. In contrast, Run 4 was characterized using a lower average sediment concentration. Specifically, the mean value used in Run 4 (14 170 mg L⁻¹) represents approximately 80 % of the TOC value used in Run 5 (17 712 mg L⁻¹), which was derived from observed TOC measurements (see Table 3). This comparison suggests that a more realistic estimate is that about 80 % of the total organic carbon exists in particulate form, with the remainder composed of dissolved organic carbon. Run 4 and Run 5 show negligible differences in the predicted water temperature and DO profiles (Figs. 8 and 9). Table A10 presents the performance metrics for water temperature, DO, TN, TP, BOD5, and Chl a obtained for Run 4. While this run improved the DO simulation in the reservoir, results for the other constituents remained very similar to those of Run 5 (baseline). Overall, the water temperature profiles are very well captured by all models (Fig. 8), reflecting their robustness in simulating thermal dynamics. In contrast, DO profiles are more complex and challenging to model due to their sensitivity to multiple interacting processes. Nevertheless, the models were able to capture the main seasonal and vertical trends in DO concentrations, including stratification patterns and general oxygen depletion in bottom layers during warmer months (Fig. 9).

Figure 7 (a) SOD values from the reservoir bottom layer, predicted by the SD_model for Runs 1 to 6, compared with the RMSE obtained between the predicted water column DO profiles and the mean initial POC values (across all sites values) for each run of the SD_model. (b) Similar to (a) but considering the NSE metric.

Figure 8 Observed water temperature profiles (300 m from the dam) compared to predicted profiles using the SD model (Run 4) and (Run 5 – baseline), Zero-order model (zero-order SOD = 2.5 g O₂ m⁻² d⁻¹ – baseline); First-order model (ISC = 0.5 g m⁻² – baseline) and the Hybrid model (zero order SOD = 1.0 g O₂ m⁻² d⁻¹ – baseline).

Figure 9 Observed DO profiles (300 m from the dam) compared to predicted profiles using the SD model (Run 4) and (Run 5 – baseline), Zero-order model (zero-order SOD = 2.5 g O₂ m⁻² d⁻¹ – baseline); First-order model (ISC = 0.5 g m⁻² – baseline) and the Hybrid model (zero order SOD = 1.0 g O₂ m⁻² d⁻¹ – baseline).

The sensitivity analysis also involved varying the initial values of PON and POP for each run. The results indicate that mean reservoir SOD values remained nearly constant, as depicted in Fig. 10, suggesting that the SD model was not significantly affected by variations in the initial PON and POP values in the sediments. However, in Runs 7, 8, and 9, where PON values were higher, there was a significant increase in the release of N-NH₄ and N-NO_x from the reservoir sediments, leading to an impact on water column DO. This is evidenced by the notable increase in RMSE and the reduction of NSE values, as shown in Fig. 10a and b.

Figure 10 (a) SOD values from the reservoir bottom layer, predicted by the SD model for Runs 7 to 15, compared with the RMSE obtained between the predicted water column DO profiles and the mean initial PON values (across all sites) for each run. (b) Similar to (a) but considering the NSE metric. (c) SOD values from the reservoir bottom layer, predicted by the SD model for Runs 16 to 24, compared with the RMSE obtained between the predicted water column DO profiles and the mean initial POP values (across all sites) for each run. (d) Similar to (c) but considering the NSE metric.

Figure 11 shows the RMSE (Fig. 11a) and the NSE (Fig. 11b) values between observed and predicted water column DO profiles for all models: SD model (Runs 1 to 6), Zero-order model and Hybrid model, each with six different SOD values ranging from 0.5 to 3.0 g m⁻² d⁻¹, along with the corresponding reservoir SOD values. Additionally, this figure illustrates how the First-order model varies with the initial sediment concentration. Among the four models evaluated, the Hybrid model demonstrated the best overall performance in predicting DO concentrations in the reservoir. With an average SOD of 1.49 g O₂ m⁻² d⁻¹, the hybrid model achieved the lowest RMSE (1.87 mg L⁻¹) and highest NSE (0.76), demonstrating superior predictive accuracy. The Zero-order model followed closely, reaching optimal performance at an average zero-order SOD of 1.43 g O₂ m⁻² d⁻¹, with an RMSE of 1.965 mg L⁻¹ and an NSE of 0.732. The SD model also performed well, attaining its best accuracy at an average SOD of 1.07 g O₂ m⁻² d⁻¹, where the RMSE decreased to 2.011 mg L⁻¹ and the NSE peaked at 0.716; however, further improvements plateaued beyond this point. In contrast, the First-order model consistently exhibited higher RMSE values (ranging from 2.15 to 2.22 mg L⁻¹) and lower NSE values (between 0.66 and 0.68), regardless of the initial sediment concentration. Moreover, its SOD at the bottom layer remained relatively stable, indicating limited sensitivity to input variations. Overall, these results underscore the hybrid model's robustness and accuracy, followed by the Zero-order and SD models, while the First-order model demonstrated the weakest performance in this context.

Figure 11 (a) RMSE between observed and simulated DO profiles in the water column for all models: the SD model (Runs 1–6), the Zero-order model, the Hybrid model with six SOD values ranging from 0.5 to 3.0 g O₂ m⁻² d⁻¹, and the First-order model with initial sediment organic matter concentrations from 0.0 to 3.0 g m?⁻². (b) Same as (a), but using the Nash–Sutcliffe Efficiency (NSE) as the performance metric.

Figure 12 Time series of DO, averaged across all model layers and segments, for each baseline model scenario: SD model (Run 5), Zero-order model (SOD = 2.5 g O₂ m⁻² d⁻¹), First-order model (initial sediment concentration = 0.5 g m⁻²), and Hybrid model (Zero-order SOD = 1.0 g O₂ m⁻² d⁻¹). The figure compares baseline conditions with an 80 % reduction in organic matter inflow load in the main reservoir branch (Branch 1 – Tâmega River). Performance metrics (R², RMSE, MAE, NSE, and PBIAS) are also shown for each case.

Figure 13 Time series of DO, averaged across all model layers and segments, for each baseline model scenario: SD model (Run 5), Zero-order model (SOD = 2.5 g O₂ m⁻² d⁻¹), First-order model (initial sediment concentration = 0.5 g m⁻²), and Hybrid model (Zero-order SOD = 1.0 g O₂ m⁻² d⁻¹). The figure compares baseline conditions with an 80 % reduction in organic matter and P-PO₄ inflow loads in the main reservoir branch (Branch 1 – Tâmega River). Performance metrics (R², RMSE, MAE, NSE, and PBIAS) are also shown for each case.

3.4 Inflow organic matter and phosphorus load reduction scenarios

The results reveal clear differences in model sensitivity to inflow load reductions, with the First-order and Hybrid models exhibiting a stronger response compared to the SD and Zero-order models (Figs. 12 and 13). The SD model showed minimal change, indicating limited sensitivity to external loading (Figs. 12a and 13a), likely due to strong internal loading feedback from legacy phosphorus and organic matter stored in sediments. The Zero-order model demonstrated limited utility for management scenarios because it is decoupled from the water column, reducing its responsiveness to external changes. The First-order model may overestimate sensitivity as it tends to underestimate internal loading contributions. The Hybrid model, which combines both approaches, is less reactive than the First-order model due to the influence of the Zero-order component, offering a more balanced response. However, the Zero-order SOD component in the Hybrid model depends solely on temperature and remains decoupled from water column conditions; this limitation may gradually reduce the model's accuracy in long-term simulations. These differences in model sensitivity are further reflected in the evolution of average SOD across scenarios (Table 5). While the Zero-order and SD models show virtually no change in bottom-layer SOD under reduced loading conditions, the First-order and Hybrid models register clear declines. The First-order model's SOD drops from 0.87 to 0.42 g O₂ m⁻² d⁻¹ (80 % OM reduction) and 0.29 g O₂ m⁻² d⁻¹ (80 % OM and P reduction) and the Hybrid model from 1.49 to 1.07 g O₂ m⁻² d⁻¹ (80 % OM reduction) and 0.94 g O₂ m⁻² d⁻¹ (80 % OM and P reduction).

Table 5 Average sediment oxygen demand (SOD) in the bottom layers of the reservoir, calculated across all segments, for each model under three scenarios: Reference (baseline conditions), 80 % reduction in organic matter inflow (OM 80 %), and combined 80 % reduction in organic matter and phosphorus inflow (OM and P %) in the in the main reservoir branch (Branch 1 – Tâmega River).

4 Discussion

Overall, the temperature and DO predictions for the reservoir boundary conditions (Tâmega river) were quite good: PBIAS: 0.76 % and 0.92 %, respectively. When a significant number of samples and forcing variables are available the accuracy of machine learning algorithms can be greatly enhanced. This was demonstrated in the studies by Lu and Ma (2020), Rajesh and Rehana (2021), and Feigl et al. (2021), where the RMSE for river water temperature prediction reached 1.04, 1.03, and 0.58 °C, respectively. The results obtained for alkalinity, conductivity and TSS were also good: Alkalinity-PBIAS: 17.44 %; Conductivity – PBIAS: 8.23 %; TSS – PBIAS: 11.86 %. However, as expected, the PBIAS values obtained for the remaining constituents were not as favorable (Total P – PBIAS: 7.11 %; N-NO_x – PBIAS: 3.92 %; BOD₅ – PBIAS: 6.93 %; Chl a – PBIAS: 30 %). The modeling of these constituents involves complex biological, chemical, and physical processes that are harder to model accurately. However, except for Chl a, the PBIAS values were generally less than 10 %, reflecting acceptable levels of bias. Ammonium (N-NH₄) was the only parameter for which performance was significantly lower, generating a PBIAS of 28.27 %. Moriasi et al. (2015) suggest that ±10 ≤ PBIAS ≤ ±25 is indicative of a satisfactory model performance.

Based on the RMSE, the overall reservoir calibration results obtained for all constituents with all models for the 2016–2021 period were consistent with the results seen in other studies (see Table A11). The mean RMSE values for Chl a obtained with all models (SD model (run 5 – baseline): 17.72 µg L⁻¹; Zero-order model (zero-order SOD = 2.5 g O₂ m⁻² d⁻¹ – baseline): 17.78 µg L⁻¹; First-order model (ISC = 0.5 g m⁻² – baseline): 14.88 µg L⁻¹ and the Hybrid model (zero-order SOD =1.0 g O₂ m⁻² d⁻¹ – baseline): 14.88 µg L⁻¹) are aligned with the results of other modeling studies (Brito et al., 2018: 62.9 µg L⁻¹; Kim et al., 2019: 6.7 to 13.2 µg L⁻¹; Tasnim et al., 2021: 0.6 to 27.6 µg L⁻¹; Almeida et al., 2023: 19.36 to 25.57 µg L⁻¹). For TP, the mean RMSE values were 0.03 mg L⁻¹ for both the SD model (Run 5 – baseline) and the First-order model (ISC = 0.5 g m⁻² – baseline), while the Hybrid model (zero-order SOD = 1.0 g O₂ m⁻² d⁻¹ – baseline) showed a slightly higher value of 0.04 mg L⁻¹. These results fall within the range reported in previous studies, including Brett et al. (2016) at 0.012 mg L⁻¹, Kim et al. (2019) between 0.014 and 0.068 mg L⁻¹, Tasnim et al. (2021) from 0.005 to 0.036 mg L⁻¹, and Almeida et al. (2023) ranging from 0.07 to 0.09 mg L⁻¹. The only exception was the Zero-order model (SOD = 2.5 g O₂ m⁻² d⁻¹ – baseline), which overestimated phosphorus export from sediments during the summer months (July to September) of 2018 to 2021, resulting in a notably higher RMSE of 0.1 mg L⁻¹. Even with a very low phosphorus release rate from the sediments – representing a fraction of the SOD (0.001) – the Zero-order model still overestimated phosphorus concentrations, particularly during periods of elevated sediment oxygen demand. This suggests that the model may lack the sensitivity needed to accurately simulate low-level sediment-phosphorus interactions under such conditions. The mean RMSE values obtained for TN were lower than the only reference value available in the literature – 0.77 mg L⁻¹ reported by Deliman and Gerald (2002). Specifically, the SD model (Run 5 – baseline) yielded an RMSE of 0.33 mg L⁻¹, the First-order model (ISC = 0.5 g m⁻² – baseline) produced 0.36 mg L⁻¹, and the Hybrid model (zero-order SOD = 1.0 g O₂ m⁻² d⁻¹ – baseline) resulted in 0.35 mg L⁻¹. The only exception was the Zero-order model (SOD = 2.5 g O₂ m⁻² d⁻¹ – baseline), which had a significantly higher RMSE of 0.79 mg L⁻¹ – slightly exceeding the value reported by Deliman and Gerald (2002) yet still within a comparable range. The RMSE obtained with the SD model (Run 5 - baseline), Zero-order model (zero-order SOD = 2.5 g O₂ m⁻² d⁻¹ – baseline); First-order model (ISC = 0.5 g m⁻² – baseline) and the Hybrid model (zero order SOD = 1.0 g O₂ m⁻² d⁻¹ – baseline) for DO, 2.01, 1.97, 2.15 and 1.87 mg L⁻¹ respectively) are also in line with the results obtained in other studies (e.g., Deliman and Gerald (2002): 1.34 mg L⁻¹; Brett et al., 2016: 1.2 mg L⁻¹; Brito et al., 2018: 7.6 mg L⁻¹; Luo et al., 2018: 1.78 mg L⁻¹; Tasnim et al., 2021: 2.33 mg L⁻¹). In the SD model (Run 5 – baseline), bottom-layer SOD values ranged from 0.015 to 5.152 g O₂ m⁻² d⁻¹ ($\mathit{\mu }=\mathrm{1.162};$ σ=0.823), reflecting moderate variability driven by seasonal biogeochemical processes. In comparison, the Zero-order model (SOD = 2.5 g O₂ m⁻² d⁻¹ – baseline) showed a broader but more temperature-driven range, from 0.000 to 15.640 g O₂ m⁻² d⁻¹ (μ=1.432; σ=2.122). The First-order model (ISC = 0.5 g m⁻² – baseline) yielded values between 0.000 and 20.000 g O₂ m⁻² d⁻¹, with a much lower mean (μ=0.870) and relatively high variability (σ=1.920), consistent with its sensitivity to organic matter loading. The Hybrid model (zero-order SOD = 1.0 g O₂ m⁻² d⁻¹ – baseline) incorporated both zero- and first-order processes and produced the widest overall range, from 0.000 to 21.938 g O₂ m⁻² d⁻¹ (μ=1.491; σ=2.024), highlighting its enhanced responsiveness to both physical (e.g., temperature) and biogeochemical (e.g., organic matter) drivers. The monthly variation in SOD across the four models reveals distinct seasonal patterns influenced by their underlying formulations (Fig. A2). All models show notable peaks in May and October, corresponding to periods of elevated organic matter inflow, while a consistent decline is observed during the summer months (June to August), when external organic inputs are comparatively low. The Zero-order model (baseline SOD = 2.5 g O₂ m⁻² d⁻¹) exhibits a sharp rise from winter to a peak of 1.919 g O₂ m⁻² d⁻¹ in May, then gradually declines over the summer, before increasing again in October (1.910 g O₂ m⁻² d⁻¹). A similar double-peak pattern is observed in the Hybrid model (zero-order SOD = 1.0 g O₂ m⁻² d⁻¹, baseline), with SOD reaching 1.715 g O₂ m⁻² d⁻¹ in May and a more pronounced maximum of 2.338 g O₂ m⁻² d⁻¹ in October, reflecting the combined effects of temperature and organic matter availability. The SD model (Run 5 – baseline) shows more moderate seasonal variation, with values dipping to 0.679 g O₂ m⁻² d⁻¹ in August, then rising to 1.501 g O₂ m⁻² d⁻¹ in November, consistent with internal sediment dynamics. The First-order model (ISC = 0.5 g m⁻², baseline), which is most sensitive to organic matter loading, also mirrors this seasonal structure, peaking in October (1.235 g O₂ m⁻² d⁻¹) after a gradual summer decline. Collectively, these patterns underscore the importance of organic matter availability – particularly in spring and autumn – as a key driver of SOD across the different modeling approaches. This pattern indicates the model's responsiveness to both organic matter inputs and temperature, leading to a more nuanced representation of seasonal variation compared to the other models. These values are consistent with the SOD values obtained in other studies, such as those of Schnoor and Fruh (1979), which concluded that the SOD values of Lake Lydon B. Johnson (located in the US) ranged from 1.7 to 5.8 g O₂ m⁻² d⁻¹, and of Beutel (2015), which measured SOD values in different locations around Lake Hodges (located in the US) ranging from 0.6 to 2.3 g O₂ m⁻² d⁻¹. It would be useful to be able to compare these results with SOD values measured at different sites within the Torrão reservoir.

It is important to emphasize that this study was primarily designed to evaluate the performance of the sediment diagenesis model. However, by incorporating alternative SOD modeling approaches, it inevitably allowed for a comparative ranking of model performance, highlighting the relative strengths and limitations of each formulation. The performance limitations of the Zero-order and First-order models can be attributed to their structural simplifications. Specifically, the Zero-order model's strong temperature dependence, coupled with its disregard for the dynamics of organic matter loading, reduces its ability to capture temporal variability driven by external inputs. Similarly, the lower accuracy of the First-order model likely stems from its exclusion of anaerobic decay processes and limited representation of sediment biogeochemistry, which becomes especially relevant under low-oxygen conditions. The Hybrid model outperformed all other approaches. Considering the principle of parsimony (Occam's razor) (Burnham and Anderson, 2002), the simpler Hybrid model proved more effective than the complex SD model, making it the preferred choice for simulating SOD dynamics in the reservoir. These findings underscore the importance of selecting models that align with the specific characteristics of the system being studied. Simpler models, such as the Hybrid model, may be adequate for steady-state conditions, short- to medium-term forecasts, or scenarios with limited data. The zero-order SOD component of the Hybrid model relies solely on temperature and is decoupled from the water column; therefore, in long-term simulations, this limitation can gradually undermine the model's accuracy. In contrast, the SD model may be more appropriate when the goal is to explore system-wide feedbacks and temporal dynamics over extended periods – especially those involving sediment accumulation and nutrient cycling – where it may provide valuable insight into underlying processes, provided that sufficient observational data become available to support its additional state variables. Moreover, a model's effectiveness heavily depends on the user's familiarity with its structure and their skill in calibration. Yet, it is unrealistic to expect researchers to master the implementation of every available modeling approach. As such, comparisons between models should be interpreted carefully, acknowledging the influence of user expertise on performance outcomes (Piccolroaz et al., 2024). Overall, to strengthen the analysis, it is recommended that users apply all available SOD modeling approaches in the case of the CE-QUAL-W2 model and assess the model's behavior. This comprehensive evaluation provides a solid foundation for further modeling efforts and helps ensure that the chosen approach is well-suited to the system's specific conditions and objectives.

The results also revealed that the particulate fraction of organic carbon in the reservoir sediments corresponded to 80 % of the TOC. This value is small compared to the results obtained for Taihu Lake by Yu et al. (2022), where the ratio of POP to TOC varied from 97.85 % to 89.53 %. However, this value (80 %) was obtained indirectly through the analysis of the reservoir's predicted SOD values as a function of different initial POC values and may, therefore, reflect other sources of uncertainty, such as inflow organic matter characterization. Given the fact that the magnitude of TOC in the sediment can be affected by numerous factors, including water column productivity, terrestrial inputs of organic materials, sediment properties, and microbial activity rates (Gireeshkumar et al., 2013), and that, partly due to differences in reservoir productivity and morphology, the spatial distribution and sources of organic carbon vary greatly across regions (Anderson et al., 2009), it is reasonable to assume that the only way to accurately assess the POC prediction is by monitoring the reservoir POC content. Furthermore, this study has highlighted the need to expand research to additional waterbodies across diverse regions to improve our understanding of the CE-QUAL-W2 diagenesis model's performance under varying environmental conditions. This includes evaluating its applicability in long-term scenarios, which are essential for capturing cumulative sediment dynamics and climate-driven trends. Additional SOD monitoring studies need to be conducted in lakes and reservoirs and extended to other latitudes, with particular focus on the chemical characterization of sediments and the definition of sediment burial rates.

5 Conclusions

This research evaluates the performance of the CE-QUAL-W2 v4.5 sediment diagenesis model in simulating water temperature, dissolved oxygen, total phosphorus, total nitrogen, chlorophyll a, and biochemical oxygen demand in a Portuguese reservoir over the period from 2016 to 2021. Calibration was based on 35 sets of observed temperature and dissolved oxygen profiles, supplemented by six annual measurements of total nitrogen, total phosphorus, chlorophyll a, and biochemical oxygen demand collected at various depths. To evaluate model accuracy, three alternative sediment oxygen demand formulations – a Zero-order model, a First-order model, and a Hybrid approach combining features of both – were also applied and compared. The Hybrid model consistently outperformed the other formulations, striking an effective balance between accuracy and simplicity. It therefore represents the most suitable choice for modeling the reservoir. In contrast, the Zero- and First-order models exhibited limitations related to temperature dependence and inadequate sediment process representation, respectively. Simpler models, such as the Hybrid model, may be adequate for steady-state conditions, short- to medium-term forecasts, or scenarios with limited data. In contrast, the SD model –despite its good performance – may be more appropriate when the goal is to explore system-wide feedbacks and temporal dynamics over extended periods, especially in cases involving sediment accumulation and nutrient cycling. In such contexts, it may offer valuable insights, provided that sufficient observational data are available to support its additional state variables. Overall, the study reinforces the importance of choosing models based on site characteristics, available data, and simulation goals. Future work should broaden the evaluation of these models across various waterbodies and extended timeframes, while highlighting the need for enhanced sediment monitoring to support detailed process-based modelling.

Appendix A

Figure A1 CE-QUAL-W2 bathymetry – Cross section of the Tâmega River with the average segment width.

Figure A2 Average monthly sediment oxygen demand (SOD) at the reservoir bottom layer predicted by the SD model (Run 5 – baseline), Zero-order model (baseline: SOD = 2.5 g O₂ m⁻² d⁻¹), First-order model (baseline: ISC = 0.5 g m⁻²), and Hybrid model (baseline: zero-order SOD = 1.0 g O₂ m⁻² d⁻¹). Also shown is the inflow BOD5 load from the reservoir's main branch.

Table A1 Model metrics. Characterization of Tâmega river inflow.

Table A2 Rates and constants: Hydraulic coefficients.

Table A3 Rates and constants: Extinction coefficients.

Table A4 Rates and constants: phytoplankton (Diatoms).

Table A5 Rates and constants: organic matter.

Table A6 Rates and constants: nutrients.

Table A7 Rates and constants: SOD rates.

Table A8 Rates and constants: Sediment diagenesis model (considering 5 regions).

Table A9 Metrics between observed and predicted values for both models. The predicted values were compared with observed values at three different depths: (a) an integrated sample between the reservoir surface and an average depth of 5.8 m, (b) an average depth of 23 m, and (c) an average depth of 43.  m. The values shown in this table represent the mean value of the metrics obtained for each date and the corresponding standard deviation.

Table A10 Metrics between observed and predicted values for SD model (run 4). Water temperature and DO metrics were obtained from 36 observed and predicted profiles. The predicted values for the remaining constituents were compared with observed values at three different depths: (a) an integrated sample between the reservoir surface and an average depth of 5.8 m, (b) an average depth of 23 m and (c) an average depth of 43.7 m. The values in this table represent the mean value of the metrics obtained at each date and the corresponding standard deviation or, in the case of water temperature and DO, the mean value of the metrics obtained for each profile and the standard deviation.

Table A11 Root mean square error values obtained with different models and across different time frames.