Finite volume model for the simulation of 1D unsteady river flow and water quality based on the WASP

How do the kinetics process and hydraulics impact on water quality in a real environment? The answer to this question began in 1912 when the Royal Commission on Sewage Disposal put into practice the biochemical oxygen demand (BOD) and dissolved oxygen (DO) measurement and interpretation in natural environments. Since then, many works have been carried out in this field. As an example, Streeter & Phelps (1958) developed a seminal work on how to assess the impact of the discharge BOD on DO in a river system. After that, the use of mathematical models has become more and more important in DO modelling for streams and estuaries. Other works like O'Connor (1961) and Thomann (1963) observed that the BOD parameter can be divided into carbonized and nitrified or that the BOD reduction can be carried out without oxygen consumption.

Another advance in terms of water quality came from the hand of the continuous progress of computers in order to improve the prediction of pollutant transport. This prediction is still a subject of interest at present, due to the conversion processes and coupling with hydrodynamics.

Based on the conversion processes, many models have emerged, such as QUAL (Brown & Barnwell 1987), Water Quality Analysis Simulation Program (WASP) (Ambrose et al. 1993) or RWQM1 (Reichert et al. 2001), that have been used to simulate and predict the detrimental effect of the discharge of organic waste on the oxygen concentration.

With respect to the coupling, some of the aforementioned quality models have been combined with 1D and two-dimensional (2D) hydrodynamic solvers such as MIKE11 (DHI 2009), HSPF (Ambrose et al. 1989) and HEC-RAS (Brunner 2010) to simulate the kinetic processes involving DO in lowland rivers. From this combination, numerous studies on modelling have been reported, but the WASP quality model has been the most widely used in several countries with different hydrodynamic models.

The WASP model was used to evaluate the parameter sensitivity in the lower course of the Saale river (Lindenschmidt 2006). It was also reported to obtain a longitudinal profile of DO and to quantify the shift of the water quality under low flow conditions (Huang et al. 2017). The WASP model has been considered to assess the effects of installing a new wastewater treatment plant for treating wastewater from a residential area of Malaysia (Lee et al. 2017). In Jian et al. (2014), a three-dimensional unstructured mesh ecological model was developed to investigate the algal bloom in the Xiangxi River. The pollution source with the WASP was also analysed by Yao et al. (2015) to evaluate the risk of point-wise pollution on a regional scale in the Taipu River.

These mathematical models have proved to be useful for analysing or predicting the fate and transport of contaminants in rivers in the context of steady flow. However, sometimes it is necessary to evaluate the water quality evolution in the presence of discharge variation. Therefore, the water quality model and the hydrodynamic model must be applied in a coupled way, i.e., solving both at the same time rather than assuming a steady-state flow model. In the context of unsteady flow, Cea et al. (2016) evaluated their freeware simulation tool (IberWQ), which computes the spatial and temporal evolution of several species, showing some of the 2D model's capabilities and limitations. This evaluation requires processes of the calibration of coefficients, for a later numerical validation with experimental data in more complex transport environments. Despite many contributions, no combined model provides all the required functionality. In other words, all models have different limitations and strengths that still need to be studied in more detail (Cox 2003; Kannel et al. 2011; Liangliang & Daoliang 2015). Some of these limitations refer to (a) that it is not possible to correlate the temporal variability of the flow with the quality parameters in different sites, (b) it is not possible to correlate the flow variability between the main river, its tributaries and discharges or (c) that it does not properly consider the seasonal or diurnal effects on quality.

In this framework, the main novelty of the present work is the development of a simulation model in 1D to predict the fate and transport of pollutants and hydrodynamic variables in unsteady river flow. The hydraulic modelling is based on the 1D shallow water equations written to consider realistic and irregular cross-sections.

For the quality module, the equations of the WASP, including the temperature, are used. In this water quality model, the variables considered belong to three groups: suspended solids removed from the river bottom, chemical compounds by anthropogenic activities (bacteria and virus due to human activities) and chemical compounds released by weathering processes. The interaction of these biochemical species generates quality indicators (such as BOD and DO) frequently measured in the natural water bodies. The most frequent state variables in a river water quality model are algae, nutrients (phosphorus cycle and nitrogen cycle) and DO. They have been included in the present work together with their specific biochemical processes. The formulation follows the model EUTRO, as included in the WASP5 simulation package, because it allows us to predict the kinetic processes more real (Ambrose et al. 1993).

A finite volume approximation has been used to solve the set of equations leading to the calculation of the evolution of flows, depths and all the water quality variables both in space and time. This approach is also a novelty since it allows us to predict the flow evolution under all flow regimes and the water quality variables in complex scenarios that consider both the temporal variability of the flow and the superficial changes. Previous works have reported the convenience of a careful formulation and discretization of these equations to deal with unsteady flow over a dry irregular bed (Burguete & Garcia-Navarro 2001; Murillo & Garcia-Navarro 2014). Furthermore, previous authors have indicated the necessity to implement a coupled discretization of the shallow water equations and the transport equations in order to avoid numerical difficulties. This idea was included in Burguete & Garcia-Navarro (2001); Murillo et al. (2009) for the 1D case and in Cea & Vázquez-Cendón (2012) for the 2D case using the pure advection of a passive solute. In the present work, this has been extended to include a water quality model with the ambition to be able to deal with unsteady flow in complex and realistic river cross-sections.

The numerical method applied for the advective part is an explicit first-order Roe's scheme which has been previously proved robust and well behaved for this kind of problem (Burguete et al. 2008). The novelty is the dynamical coupling of the models and their simultaneous numerical resolution using a single finite volume scheme based on a global time stepping that dynamically controls the numerical stability of the solution either under steady or unsteady conditions. This is achieved as well as offering the possibility to involve all sorts of source terms in the equations (bed slope, cross-section changes, friction, heat flux and all the water quality kinetic terms). Furthermore, the numerical formulation ensures conservative solutions which are crucial when dealing with the transport processes in water quality simulations. The model is validated with theoretical cases, with well-known analytical solutions and with a real case of the Ebro river (Spain).

The fate and transport processes of contaminants are controlled by two factors: their reactivity and their hydrodynamic transport. The formulation of both models is next presented.

The above formulation is widely accepted as a good approximation for long reaches when considering sections of real rivers.

The advective transport is a result of passive movement along with the flow. The diffusive transport caused by Brownian motion as well as turbulence is influenced by the fluid temperature and the size of the molecules. Also, this term will be affected by the slope, morphology and roughness. This means that it varies according to the direction of the flow. For this case, we assume that the dispersion coefficient is independent of any direction, so that the surface water body is homogeneous and isotropic. In order to compute this term, many empirical and analytical expressions exist in the literature (Chen et al. 2012; Wang et al. 2017).

In Equation (5), the term governs the permanence of pollutants in a body of water depending on the nature of the compound. Most chemicals undergo biological or chemical decay. For this study, the solutes are non-conservative, that is, the pollutants can chemically or biologically decay. These fate and decay processes include volatilization, hydrolysis, photolysis and biodegradation and have been addressed and characterized by Thomann & Mueller (1987); Chapra (2008); Ji (2017), among others.

The list of the kinetic processes considered is summarized in Table 1. A detailed description of each decay rate for each process of Table 1 is given in Appendix Notation in Supplementary Materials.

The formation or disappearance of each constituent in the water column, due to its inherent process, will be discussed in detail in the following subsections.

The DO (⁠⁠) is one of the most important parameters of water quality, because it is a basic requirement for a healthy aquatic ecosystem. The DO concentration in a stream can change through an exchange with the atmosphere (coefficient ⁠) and the growth of algae (photosynthesis). However, its depletion is due to the oxidation of organic carbon (affected by coefficient ⁠), nitrification (⁠⁠), the death of algae (respiration ⁠) and the sediment oxygen demand (⁠⁠), defined as the rate of DO required for the oxidation of organic matter in benthic sediments. The DO concentration frequently oscillates in the water column, but its oscillation is higher when it arises from the human activity. The complete process that includes all these gains and losses of DO in the water column can be expressed as ⁠.

The carbonaceous BOD (CBOD; ⁠) is the concentration of organic material present in the water body. This process includes the effects of sedimentation, oxidation (⁠⁠) and denitrification (⁠⁠). The principal sources of CBOD are man-made sources, algal death (⁠⁠) and natural runoff. There is a mutual interaction between the CBOD and the DO components, so that, in particular, the DO level will decrease when the injection of CBOD is continuous in time, causing an oxygen deficit.

The quality of a body of water can be affected by the presence of phytoplankton (⁠⁠). This population can be accelerated by the addition of nutrients (nitrogen and phosphorus), either by human activities or natural processes. The excess of nutrients provides more population growth. This uncontrolled growth is commonly referred to as eutrophication. When this population becomes large, it may cause diurnal variations in DO that can be fatal to fish life. Also, the presence of phytoplankton can cause water taste and odour problems.

The model considers two of the three primary dependent systems: the phytoplankton population and the nutrient system. The external environment variables that affect those systems are temperature, advective flow and solar radiation. The classical approach is to assume that these effects are multiplicative. Thereby, the increase of phytoplankton in rivers and streams is due to the availability of nutrients and solar energy, while its reduction occurs primarily through respiration. A simplified representation of this process is given by ⁠.

Nitrogen can be found in five major forms in aquatic environments: organic nitrogen (ON), ammonia ⁠, nitrite ⁠, nitrate and dissolved nitrogen gas ⁠. The sequential processes of nitrogen compounds transforming ON to ⁠, then to and finally to ⁠.

The ON (⁠⁠) originally present in water is partially transformed into and partially settles on the bottom; meanwhile, the increase of ON is due to phytoplankton death. The ON kinetic equation describing these processes is expressed as ⁠, involving (⁠⁠).

(⁠⁠) is one of the intermediate compounds formed during biological metabolism and, together with ON, is considered as an indicator of recent pollution. The ammonia concentration is increased by the change of ON to due to mineralization and the production of nitrogen due to phytoplankton death and respiration. At the same time, its concentration is reduced by the uptake for phytoplankton growth and the change of to due to nitrification. The corresponding kinetic process for this state variable is ⁠.

(⁠⁠), the end product of nitrification, represents the sum of and (the amount of present in natural waters is usually very small). Its reduction is due to the growth of phytoplankton and the denitrification process (⁠⁠); meanwhile, its rise is generated by the nitrification process (⁠⁠). Therefore, the total nitrate concentration can be expressed as ⁠.

To model the phosphorus cycle, two single kinetic equations are taken into account: organic phosphorus (OP) (⁠⁠) and inorganic phosphorus ⁠.

The production of OP is caused by phytoplankton death and respiration. Conversely, the loss of OP is due to mineralization and the settling process (⁠⁠). The kinetics of is affected by the uptake for phytoplankton growth. However, the production of phosphorus is caused by phytoplankton death, respiration and mineralization. Hence, the mathematical equation for can be expressed as ⁠.

All the kinetic equations that describe the increase or decrease of state variables in time as well as all the coefficients controlling the rate of the processes are summarized in Tables 1 and 2.

In this study, following Edinger et al. (1968), a temperature transport model has been included assuming it as a scalar variable whose evolution can be formulated following also an advection–reaction equation. The total heat budget for a water body includes the effects of water depth, velocity and atmospheric conditions. For that purpose, it is useful to estimate the daily average stream temperature based on climate conditions (Gu & Li 2002; Herb & Stefan 2011).

The heat flux on the interface of the water column and the sediment bed also changes the water temperature, but this exchange with the sediment is much smaller than the surface exchange; hence it is neglected.

The temporal and spatial distribution of all the water quality constituents, temperature, water depth and discharge can be obtained by solving the problem from the initial conditions, assuming that boundary conditions and kinetic parameters are known.

An explicit first-order upwind scheme has been used to solve the flow Equations (1) and the advective part of (5). It has been used, because it allows us to have greater control of the numerical stability when dealing with equations where source terms are as relevant as flux terms (Murillo & García-Navarro 2010). The balance between flows and source terms is accurate, therefore ensuring conservative solutions both in unsteady state and when approaching the steady state. In this way, no extra corrections are made, such as the hydrostatic reconstruction (Audusse et al. 2004). The formulation is written in a finite volume form, so that it can easily be extended to two dimensions and/or high-order schemes (Navas-Montilla & Murillo 2018).

For the numerical solution of the water quality transport (Equation (5)), three steps have been considered: advective part, diffusive part and reactive part. For the sake of clarity, the transport of only one of the solutes is presented; therefore, the sub-index i labelling the different water quality components is dropped in this section.

It is worth noting that this global step is dynamically computed and controls the stability of the coupled system.

This matrix incorporates the processes (rows) and variables of the state (columns) of the eutrophication model, and the elements within the matrix include the coefficients that establish the mass relationships between the components in the individual processes (Reichert et al. 2001). The configuration of this matrix allows an easy and fast recognition of the fate of each state variable. The term R is obtained by adding the products of the stoichiometric coefficients and the process rate for each component i that is considered.

Once the theoretical development of the equations has been presented, it is crucial to pass to the stage of testing or tuning the model. Such testing should include a consistent set of contaminant transport scenarios.

On the one hand, there exist analytical solutions to 1D advective–dispersive transport that represent basic benchmark tests for the numerical schemes, providing sensitivity analyses to the numerical parameters. They have already been validated with the proposed scheme (Gordillo et al. 2018). Conversely, it is important to evaluate the well-balanced property of the method for the advective transport in complex situations. Finally, the evaluation of the model results with field data, when available, is also indispensable to achieve a calibrated, verified and validated model.

A dambreak flows with advective solute transport in a channel of variable geometry. The formulation here presented is able to provide stable and accurate results for both steady and unsteady flows even in the presence of complex channel geometry, including discontinuous initial conditions. The proposed case presents the results obtained with the coupled scheme for the advective transport of a passive solute. The test used is a dambreak in a 100 m long frictionless channel with the variable width. The width of the channel is 5 m from = 0 to = 50 m and 30 m from = 55 to = 100 m. There is also a bottom discontinuity located at the same point of dam, the ratio is 2:0. Discontinuous initial conditions are assumed both in the water depth, with a ratio of initial heights of 4:1, and in the initial solute concentration, with a ratio of 1:0. The dam or gate is located at = 30 m. The water is assumed to be initially at rest. The initial discontinuity is equivalent to a sudden and total gate opening. No diffusion or reaction terms are included in order to evaluate the ability of the scheme to simulate the advective transport free from oscillations in this complex test case. The case used a = 0.5 m and CFL = 1 for the entire simulation.

Figure 2(a) represents the initial condition, and Figure 2(b) shows the numerical solution for the longitudinal water depth profile at several times after the dam removal. Figure 3 shows the longitudinal profile of the concentration at the same times. The discontinuity of concentration advances from 30 m from 1 to 0 without oscillations. The advection speed is altered when the concentration front enters the channel widening zone becoming slower compared to the upstream advection speed. The proposed scheme is robust and reliable giving good approaches to sharply unsteady scenarios.

An additional case is proposed based on the geometric and hydraulic conditions of Case 1. For this scenario, the initial solute concentration is assumed uniform over the full domain with a value of 1 ⁠. The objective of this case is to show that the concentration profile remains unchanged during the simulation despite the severe flow conditions. Figure 4 shows the longitudinal profile of concentration without discontinuity. As the unsteady flow develops, the passive solute maintains the initial concentration in the whole domain, verifying the proposed numerical properties again.

The model is now evaluated with two state variables (CBOD and DO) and two-point sources. The example consists of observing the influence of point discharges in a mainstream over the DO.

A trapezoidal channel 110 km long with bottom width 10 m is assumed. The bottom slope is = 0.0002 from 0 to 60 km and = 0.00018 from 60 km onwards. The roughness coefficient = 0.035 is considered in this case. A sewage treatment plant located at = 10 km releases a flow = 0.463 with 200 of CBOD and 2 of DO. The mainstream flow in the channel is = 5.785 ⁠, with an initial CBOD concentration of 2 and an initial DO concentration of 7.5 ⁠. There is a second point source (tributary) at 50 km downstream where a treatment plant discharges a flow = 1.157 conveying 5 of CBOD and 9 of DO.

The decomposition process of the CBOD is ⁠.

Even though the temperature evolution is not simulated in this steady flow case, all the coefficients have been calculated as a function of the temperature using Equation (6). Both point sources have different entry temperatures: the flow of the main stream has a temperature of 20 °C; however, that of the plant and tributary are 28 and 15 °C, respectively. The SOD often plays a significant role in affecting the DO concentration. Therefore, it is considered as something crucial in the process. Figure 5 shows this effect in the reach from = 10 to = 30 km where SOD = 5 ⁠. To see the case in detail, the reader can consult (Chapra 2008, p. 493).

The longitudinal profiles in Figure 5 show the concentrations of DO and CBOD obtained with the proposed model for the steady state. The solutions presented are very similar in comparison to the QUAL2E output. According to the DO kinetic equation described previously, the concentration falls to 4 approximately, due to the amount of CBOD present and the SOD process in the 30 km reach. Also, the temperature influences the DO concentration through the dependence of the different reaction speeds due to the increase in metabolic speeds of aquatic life.

Although the example is simple, it is useful to show the potential application of the model to assess the impact of point sources. The approach for this kind of scenario is satisfactory, encouraging the extension to a real case for water utility managers.

To continue the evaluation, the model now also incorporates the nitrogen cycle. Therefore, the state variables taken into account for this case are DO, CBOD, and ⁠. In the same way, as in Case 2, the main interest is focused on the variation in the concentration of DO following a point discharge. This case is based on a problem from Thomann & Mueller (1987) where the existence of field data is reported. All the data of the hydraulics characteristics of the channel and of the source are summarized in Table 4.

The predictions of the state variables for nitrogen are compared with the reported field data. The results are shown in Figure 6. The evolution of all the modelled solutes along the reach reproduces correctly the observed in-stream concentrations in Thomann & Mueller (1987).

As in the previous cases, the kinetic coefficients for organic matter degradation and nitrification were constant during all the simulation, but including the correction of these rates with the average flow temperature.

The eutrophication in rivers and streams owing to the excessive discharge of nutrients, such as nitrogen and phosphorus, is the main objective of this test case. The focus is put on the predictive character of the model, in order to control the growth of aquatic plants in surface water body.

According to the description in Chapra (2008, p. 653), the kinetics of the nitrogen and phosphorus are incorporated in order to analyse the impact of nutrient-rich discharges on the phytoplankton. The participation of these new processes implies that other kinetics, such as CBOD, are modified to account for the effects of nitrification and plant growth/respiration.

Besides the changes in some kinetic processes, it is necessary to include the effect of solar radiation, because, in the natural environment, the light intensity to which the aquatic plants are exposed is not uniform (Di Toro et al. 1971). The empirical equation that determines this relation is given in Table 2. Another factor of this test case is the nutrient limitation at all times.

The heat balance Equation (7) was therefore included in this example in order to be able to follow the spatial and temporal evolution of the water temperature as a consequence of the changes in radiation and atmospheric conditions.

All the parameters for the nutrient/plant interactions in this test case and some others already mentioned have been extracted from the example nutrients and algae in QUAL2E (Chapra 2008). The results are displayed in Figures 7 and 8. In Figure 7, the water temperature profiles obtained from our model and the QUAL2E model are plotted, while that in Figure 8, the profiles of DO and saturation oxygen (⁠⁠) are displayed. These profiles have been obtained after convergence to the steady state of the full set of transport equations involving the state variables, such as ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, and temperature T.

The results show that the interaction of all the processes leads to a decrease in the oxygen downstream, causing hypoxic conditions at = 30 km. In this condition of DO levels (not very far from reality), it is likely that aquatic life is put under stress, so that foul odours, tastes and undesired colours might be produced if the low DO conditions are maintained for a prolonged period of time.

This test case is useful to stress the importance of considering the actual temperature of lateral streams acting as tributary sources of contaminant as they can alter not only the temperature of the mainstream but also the distribution and concentration of other water quality components.

To assess the applicability of the proposed water quality model to a real case, the unsteady flow in a reach of the Ebro river (Spain) is next considered.

The Ebro river is the largest and longest river in Spain and has a total length of 930 km. Only the reach Alagón-Zaragoza with 44 km was considered in this example (Figure 9).

The period of validation of the model was from 3 July 2006 to 21 July 2006, and the data series for the flow rate and depth were provided by the Ebro Basin Water Authority (CHE). This public body has a fundamental work in terms of the hydrological planning, the quality of the waters, the control of the discharges, the environmental improvement and the execution of the fundamental works for the development of the territory.

The meteorological variables, such as air temperature, specific humidity and radiation, were collected from the database of the Aragonese Society of Agro-environmental Management (SARGA), which provides average daily and monthly values.

The reach Alagón-Zaragoza was considered in order to observe if the model is able to represent correctly the actual hydrodynamic flow variables of the river, to analyse the longitudinal flow behaviour in a low range of 20–50 and to verify the sensitivity of the temperature in front of these flows.

The discharges are low in this period, because it is a summer period with almost no rainfall, and there is an average discharge derivation of 20–50 approximately towards an adjacent irrigation canal. This canal runs parallel to the Ebro river, by its right margin along about 108 km. The surface irrigated by this canal is 26.500 Ha.

A contribution of flow to be considered in this section is from the Jalón river. This tributary is located 13 km downstream the Alagón station. Both the flow equations and the water temperature equation were included in the model.

The comparison of the measured and simulated time evolution of the water temperature at the downstream section located at Zaragoza in that period is shown in Figure 10. The blue line corresponds to the temperature evolution obtained using the raw daily average meteorological data available. The difference between the numerical and measured temperature can be due to an error in the measurement equipment as well as to the temporal variation in the meteorological variables throughout the day that was not properly recorded. Therefore, a calibration is necessary to obtain satisfactory results. The solar radiation and the average air temperature were adjusted over the daily average in order to generate a reasonable hourly variation (see Table 5), whereas the wind speed (3.11 m/s) and the relative humidity (0.45) were assumed constant. The calibration process consisted of performing trial and error simulations until reaching a set of parameters that adequately fit the measured and simulated values. The minimum values of solar radiation and air temperature are 50 and 20 °C, respectively. The reason is that heat exchange is also determined by wind speed and relative humidity variations that have been considered as constant in this case. The result after calibration is plotted in Figure 10 (red line).

It is important to emphasize that this data manipulation was only performed when the transported flows were in the range of 20–50 (Gu & Li 2002).

Concerning the hydraulic variables, the model was run from a dry bed initial condition to the steady state by imposing a constant discharge (⁠ = 41 ⁠) upstream boundary condition and a gauging curve downstream boundary condition. Then, a second run was performed using the actual discharge hydrograph measured by CHE at the upstream section every 15 min. The model responds reasonably well when comparing the time evolution of the flow discharge with respect to the measured data (see Figure 11).

The last case is evaluated in the same reach as Case 5, and the interest is focused on the evolution of and the temperature under high water discharges and unsteady conditions. There were no data for the rest of the components.

Monitoring stations are located in Alagón, Jalón (at Grisén) and Zaragoza, as shown in Figure 9. The station in Grisén is located approximately 11 km upstream the junction of the tributary Jalón river with the Ebro river.

All gauging stations provide the measurement of discharge, DO and water temperature every 15 min (see Figure 12). The rest of the water quality model components are characterized or directly measured with far less frequency. Therefore, they cannot be used for verification. In order to incorporate them into the full water quality model, initial values were considered zero.

The model was calibrated in the time period from 1 January 2011 to 30 June 2011 (181 days). This process consists in finding the values involved in the considered reaction mechanisms that yield the most accurate simulated results with respect to field measured data.

To begin with the procedure, initial and boundary values are required for both the hydrodynamic variables and the state variables. Conversely, values for the different reaction processes are also required. These values were initially taken from those suggested by Bowie et al. (1985) and Ambrose et al. (1989). All these values are imposed at the beginning of the simulation. Then, a trial and error process was performed to adjust the measured values with the simulated ones. This verification was carried out at the measurement station located in Zaragoza, both for hydrodynamics and quality.

The parameters obtained from this calibration are shown in Table 6. As can be seen, almost all decay rate and temperature correction factors do not vary too much from the recommended ranges. Therefore, few simulations are necessary to converge to a satisfactory solution.

The numerical solution, presented in Figure 13, depicts the simulated and measured flow-rate variation, while Figure 14 represents the water surface level variation. The flow discharge ranges between 57 and 1,087 ⁠. This period was chosen to check the efficiency and robustness of the method in the unsteady case due to the strong variations in the flow. Moreover, it included the best and most complete field data set of all the historical data. The results of this scenario provide a reasonable good agreement of the variables flow-rate and water-level surface in the reach under study with a Manning coefficient around 0.035.

Figure 15 (left) shows the comparison of the numerical temperature predictions with the observed temperature data. It showed a good agreement throughout the simulation period. Conversely, Figure 15 (right) depicts a comparison of the computed and measured DO concentration. Marked differences appeared between observed DO concentrations and those numerically predicted, especially at the concentration peaks. Nevertheless, measured and computed DO concentration distributions showed a reasonably good agreement.

Figure 15 also shows how the increase in the temperature has an effect on the DO concentration. The increase in the temperature on the one hand increases the metabolic rates of aquatic life, at the same time, the available DO is reduced. Thus, as the temperature increases, the oxygen demand increases, while the DO available concentration goes down. If this cycle is prolonged in time it could be harmful, because many organisms would be unable to adapt to the new temperature and therefore disappear.

The longitudinal profile of ammonium was not shown due to its low concentration. Therefore, the concentration of ammonia is also low because of its direct relation with the ammonium.

Finally, an analysis of the numerical model was performed with the root mean square error (RMSE) coefficients and Nash–Sutcliffle (N–S) coefficient. The evaluation presented in Table 7 shows that in both Cases 3 and 6 for and ⁠, the adjustment is satisfactory according to the RMSE. For the rest of the proposed cases, the adjustment is very good. Paying attention to the N–S coefficient, the performance of the method is acceptable in Cases 3 and 6 (for the DO variable) and optimal for other cases. The performance indexes for Case 6 can be increased by improving the control over the following issues:

• (a)

  The quality of the available data to calibrate the water quality model was limited. The data required at the boundaries for the state variables considered were not provided at a rate of 15 min as the hydraulic data. The sampling period for some variables of quality is typically characterized by a too low measuring frequency (quarterly). Therefore, the missing information was generated by linear interpolation.

• (b)

  The reaction rates for the different processes were not calibrated in situ, only adjusted with the recommended values. These rates may vary depending on the water body, so it is convenient to calibrate in situ or in the laboratory.

• (c)

  The quality of the data measured and protocols for the design of data collection at different locations along the river are incomplete. This deficiency is often due to the lack of investment on a suitable monitoring for the quality control stations.

A complete hydrodynamic and water quality model has been presented as a useful numerical tool for the numerical simulation of river water quality and temperature evolution under steady and unsteady conditions. This cross-section averaged model is able to estimate complex temporal scenarios, in each of the physical processes present in real environments. Physical processes, such as kinetic equations, tributaries, abstractions and heat fluxes, are included to model the seasonal evolution of the water depth, discharge, temperature and concentration of several components. The resulting system of equations has been solved using a finite volume formulation that ensures stable and conservative solutions under all flow regimes.

Through calibration and validation, the proposed model is able to represent the reality of the interactions between hydrodynamics and quality accurately and reliably. Therefore, the combined model serves as a predictive tool for water quality in various scenarios. This tool could be useful to the entities providing this service, because it would allow them to reduce and define correctly the location of the quality measurement stations.

A particularity of the model when evaluated during the summer period was the sensitivity of the water temperature to the flows. This sensitivity appears when the transported flow rate is less than 50 ⁠. The results showed that a careful calibration must be carried out on the daily variation of air temperature and solar radiation. It is worth mentioning that this sensitivity disappears when the transported flows are higher than 50 ⁠.

In all the proposed scenarios, the response of the model has been satisfactorily evaluated with the field data as well as the literature tests. The temporal distributions of both quality and flow variables in all cases are free of oscillations throughout the period of time under study.

Finally, another characteristic of the model is its usefulness in the academic field, because it allows the comparison of the different processes with known analytical solutions.

The suitable way to apply these alterations allows us to know the real characteristics of the ecosystems, at the same time as it helps us in the analysis of environmental problems. Being able to identify these problems is of vital importance, because they allow us to elaborate environmental management plans especially when the bodies of water are large and complex.

The author is grateful to the University of Zaragoza and Banco Santander grant (convocatoria de ayudas de movilidad para latinoamericanos 2016–2017) for the support granted. This work is also a part of the project CGL2015-66114-R granted by the Spanish Ministry of Science and Competitiveness (MINECO).

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/hydro.2019.080.