Water Quality Model for Streams: A Review

S. Ranjith¹, Anand V. Shivapur², P. Shiva Keshava Kumar³, Chandrashekarayya G. Hiremath⁴, Santhosh Dhungana⁵

¹VTU-PG Studies, Belagavi, India.
²Department of Civil Engineering, VTU-PG Studies, Belagavi, India.
³Department of Civil Engineering, PDIT Engineering College, Hosapete, India.
⁴Department of Water and Land Management, VTU, Belagavi, India.
⁵School of Environment, Resources and Development (SERD), Pathumthani, Thailand.

DOI: 10.4236/jep.2019.1012097                

Abstract

Six main public domain water quality models which are presently available for Rivers and streams are being captured in this article. These main models could produce important results if they are used in the correct manner, because they are different in terms of assumptions, strength and weaknesses, processes they represent, modeling capability and data in-put requirements. The Model review discussed includes, water quality analysis simulation program (WASP7), simulation catchment (SIMCAT), quality simulation along Rivers (QUASAR), and the temporal overall model for catchment (TOMCAT), QUAL2KW, QUAL2EU. The models are described individually according to a consistent set of criteria-conceptualization, model capability, model strengths, limitations, input data and how it utilized. The outcome showed that TOMCAT and SIMCAT are important in ASSESSING effect of point sources in a very simple way. The QUAL2KW, unlike the QUAL2EU where macrophytes play a major interaction, it can convert algal death to carbonaceous Bio-chemical Oxygen Demand (CBOD), thereby making it more suitable. In addition to the extensive requirement of data, it is expensive and time consuming to set up these complex models such as QUASAR and WASP7. Therefore, one model cannot be used for all the required functionalities. Choosing a model would depend on a specific application, financial cost and time availability. This article may be of help in choosing a suitable model for a specific water quality problem.

Keywords

Water Quality Models, Dissolved Oxygen, SIMCAT, TOMCAT, QUAL2Kw, QUAL2EU, WASP7 and QUASAR

1. Introduction

There are various ways of describing the mathematical model. A model described by the encyclopedia of life support system portrays it as an approximate description of a class of real-world objects and phenomena expressed by mathematical symbolisms [1]. The model was described by Concise Oxford dictionary (1990) as a simplified form of a system which aids calculations and predicts the Conditions of a system in a particular situation. Water quality models were described by US-EPA [2] as tools for simulating movement of pollutants and precipitation from the ground surface via pipe and network channels, storage treatment units and finally to receiving waters. Chapra [3] explained a mathematical model as an idealized formulation which shows the response of a physical system from external stimuli. Similarly, a water quality model was described by Cox [4] as anything derived from a simple empirical relationship, via a set of mass balance equations, to a complex software suite by which water quality in streams and rivers can be simulated by a user, by providing chemical and physical data.

An essential material required for the survival of aquatic life is Dissolved oxygen (DO). Reduced concentration of dissolved oxygen leads to an imbalance in the ecosystem with odors, aesthetic nuisances and fish mortality [4]. A strategy for water quality management, involving the assessment of pollutant effects on the dissolved oxygen concentration along the river systems can be used to attain a particular water quality. In 1925, quantitative methods were used to assess the effects of pollutants on the river systems, when Phelps and Streeter created a model for the simulation of DO in the river systems.

Following the advancement in computer technology in the 20^th and 21^st centuries, there has been numerous development in the aspect of water quality modelling which has brought about different models. Due to continuous studies and construction of new models for particular events around the globe, there is now a variety of water quality models [4]. A compromise between desirability and feasibility should be considered in selecting a simulation model [5]. Six most common and freely accessible water quality models have been addressed in this article, they include; QUAL2KW; QUAL2EU; Simulation catchment SIMCAT; water quality analysis program, WASP7; quality simulation along Rivers, QUASAR; and temporal overall model for catchments, TOMCAT, which are used for simulating dissolved oxygen along river systems and assessed each their potential for use in applications.

2. Model Review

The chosen water quality models QUAL2KW, QUAL2EU, QUAL2EU, TOMCAT, WASP7, SIMCAT are assessed with a consistent set of CRITERIA: model capability, processes, model strengths, limitations, conceptualization, input data and if their application.

2.1. SIMCAT

The SIMCAT is seen as a stochastic, one dimensional (1D), Steady state and deterministic model. Anglian Water [6], a leading provider of water and wastewater services in the United Kingdom (UK), developed SIMCAT. For over 20 years, it has been used commonly in the UK, and is known to be a not expensive, practical tool for management of water quality. The quality of the river water has been described by this toll, throughout a catchment, using the Monte Carlo simulation technique, like mean and range of percentiles [7]. A well detailed article on SIMCAT can be seen at Cox [4].

2.1.1. Conceptualization

According to SIMCAT, the stream system can be divided into diverse user well-defined spreads into any distance, which is occupied as the distance between branches or supplementary points of interest. This model can perform greater than one influence in a particular reach. A diffuser run off can be classified as quality and flow rate. The model depicts the river reaches as a continually stirred tank reactors in series (CSTRS) model. This model doesn’t utilize and advection-dispersion transport equation, but simulated instantaneous and perfect mixing along the reach, with the movement of water and solutes at the same velocity. At the beginning of every reaching, a mass balance is performed. The solute mass-balances and the flow for a specific reach include;

$Q_0=Q_i+Q_t+Q_c-Q_a$ (1)

And

$C_0{Q}_0=C_i{Q}_i+C_t{Q}_t+C_e{Q}_e-c_i{Q}_a\pm \Delta C$ (2)

where Q = flow and C = concentration of determinants o, i, t, e = reach flow, upstream input, tributary inputs, effluents discharges, and abstractions, respectively. Physical, biological, or chemical processes are internal transformations which are depicted by the term ΔC. An empirical velocity–flow relationship (V = _aQ^b, where V is the velocity, Q is the flow rate and a, b are constants) is used to get the water velocity in a particular reach. The residence time can be derived from the calculated velocity.

2.1.2. Processes

In order to calculate the concentration of determinants that will move to the next reach, the concentrations of the solute pass through first-order decays. The determinants which are being modeled might be treated as having first-order decay or conservatively, these models include; ammonium (first order), Biochemical oxygen demand (BOD first order), chloride (conservative). SIMCAT describes advective transport and chemical fate as

$\frac{\text{d}c}{\text{d}t}=-kc,\text{d}t=\frac{\text{d}x}{v}\text{and}{C}_x=C_0\cdot e-k\cdot \frac{x}{v}$ (3)

where C is the concentration, k is the rate of decay constant, x is the distance, and v is the velocity. In the case of DO modeling, the atmospheric re-aeration is included and the method of Elmore and Hayes [8] is used to estimate the DO saturation concentration as expressed by the equations,

$\frac{\text{d}c}{\text{d}t}=-k_{r}L+k_a\left(c_s-c\right)$ (4)

$c_s=14652-0.41022T+0.0079910T^2-0.000077774T^3$ (5)

where C is the DO concentration, C_s is the saturation concentration, L is the BOD, k_r is the rate of removal of BOD, k_a is the re-aeration rate coefficient, and T is the temperature in degrees Celsius.

2.1.3. Input Data

The top of the main river is where flow and quality data are entered. Each reach has it’s specific abstraction, tributaries and discharges given to it. The model utilizes the Monte Carlo technique and is stochastic; the inputs therefore, describe the statistical distribution for that determinants, and not single values. Distribution as annual means and standard deviations are accepted by the model, ranging from probability distributions, such as; normal, lognormal, constant (or uniform), 3-parameter lognormal, type III distributions, Pearson and LogPearson, and others.

2.1.4. Model Capability

Over 600 reaches and about 1400 features like rivers, abstractions, weirs, diffuse pollution and discharges can be modeled by the SIMCAT. There are four ways in which the model operates: first method utilizes data provided by the user and is used for manual calibration; second method utilizes auto-calibration algorithms to monitor quality and flow; the third method sets effluent quality so as to obtain specific river water quality objectives; the fourth method sets effluent standards to prevent reduction in water quality. This model can be assessed for auto-calibration. When auto-calibration is being utilized, the model feeds in extra flows which are a function of the river length, up to the point when the simulated flows are the same as this observed in the river at flow gauges. In order to balance observed and simulated water quality, it passes through various adjustments to the quality parameters. A summary of the statistics (means and 90^th or 95^th percentile) is compiled by the model for specific determinants for each reach. It also estimates the limit of confidence of the results, while assuming it is either a normal or lognormal distribution. The model simulates: ammonia, dissolved oxygen, Biochemical oxygen demand and parameters which are user defined. First order decay or conservative parameters can be assessed using this model.

2.1.5. Limitations of Model

Although fast, the modeling approach used in SIMCAT is too easy. No accounting of photosynthesis, sediment oxygen demand, respiration, and no allowance for temporal variability and changes in re-aeration rate with flow is it’s limiting factor. It is therefore not likely that a good result for productive rivers, would be produced by the DO model. Nevertheless, it is adequate for modeling determinants in freshwater which do not depend on sediment interactions. In addition to running the impacts of changes in the conditions of discharge of waste, it gives the user annual statistics.

2.1.6. Strengths and Its Application

Amongst all models, SIMCAT needs limited data most times. Its application is in a catchment scale. Trained non-specialist staff utilize this as a routine tool for quick assessment of management options [7]. This model was used on the river Tame, UK to assess its integrated water quality and environmental cost-benefit modeling [9] [10] and also for integrated catchment planning study for the four river catchments (Derwent, Ehen, Eden and Kent) located in North West of England [7].

2.2. TOMCAT

Thames Water of the UK water utility company, in 1980s, produced the TOMCAT [10] in order to help in reviewing waste water quality standards in every Thames Water site, in a bid to attain the aim of river-water quality. A well detailed review of TOMCAT can be seen at Cox [4].

2.2.1. Conceptualization

TOMCAT and SIMCAT have identical conceptualization, that is, a steady state CSTRS model and the approach it assumes is the Monte Carlo stochastic approach. Although, more temporal correlations are performed by TOMCAT.

2.2.2. Processes

The process equations that describe the concentrations of solutes are the same with that of SIMCAT, apart from the ones used to simulate DO and temperature. The temperature of the river (T) is said to move towards the temperature of air (T_air). The DO model includes, atmospheric aeration, oxidation of BOD, and nitrification, as explained by the equation below:

$\frac{\text{d}T}{\text{d}t}=-K_T\left(T-T_{\text{air}}\right)$ (6)

$\frac{\text{d}c}{\text{d}t}=\left(c_s-c\right)k_a-\frac{\text{d}L}{\text{d}t}-\frac{\text{d}\left[{\text{NH}}_4\right]}{\text{d}t}$ (7)

where K_T is the first-order rate coefficient, C is the concentration of DO, K_a is the re-aeration rate coefficient, C_s is the saturation concentration of DO, L is the BOD concentration and [NH₄] is the ammonium concentration. The re-aeration rate coefficient is determined from a “user supplied” re-aeration parameter (K_u), the river width (W) and the cross-sectional flow area of the channel (A) i.e.: Ka = K_U × W/A. Temperature dependency is included as a linear increase in K_a with increase in temperature.

2.2.3. Input Data

There are two types of input data required to operate TOMCAT.

1) Fixed values, usually physical parameters

2) Flow and quality data

The flow and quality data are presented as standard deviations and means of normal distribution or if logged data, as single values, or as percentage points on a nonparametric distribution. Flow and quality boundary conditions are given as seasonal or single distributions at instances, and each user-defined reach is given specific amount if reach parameters. These parameters include: depth, scale factor for runoff (i.e., the amount of the total runoff for each kilometer, which is received by the reach), mean cross-sectional area, oxygen exchange rate parameter, catchment number for calculating the (diffuse) catchment runoff, thermal equilibrium rate constant, ammonium decay rate parameter, ultimate ammonium concentration, BOD decay rate parameter, ultimate BOD concentration. For proper calibration, the observed data are added as seasonal distributions. Its limitation is that it cannot be utilized for predictive frameworks in terms of flow.

2.2.4. Model Capability

The present conditions of water-quality and flow in the catchment can be imitated using the TOMCAT model, and can also be used in evaluating what is necessary to improve water quality in the catchment. Annual and monthly statistics are given to the user when using this model. It can access the impacts of changes in waste discharge conditions rapidly. By “diverting” waste discharges to an alternative outlet, it is able to simulate the action of storm water to see if there is an increase above specific limits.

2.2.5. Limitations of Model

In terms of the processes involved, the functions of TOMCAT is limited, although seasonal statistics creates room for potentially greater accuracy compared to what the SIMCAT can do. It permits the user to obtain the results of each model run, so that statistical analyses may be done using techniques that are not in cooperated into the model. The manner in which this model simulates the flow velocity is not as accurate as that of SIMCAT, as it relies on the cross-sectional area of the river only. Where the easy processes are a reasonable approximation and don’t rely on interactions of sediment, TOMCAT is appropriate for modeling determinants in freshwater. Respiration and photosynthesis are exceptions.

2.2.6. Strengths and Its Application

This model requires smaller data and is easier than QUAL2EU. Unlike the SIMCAT model, it utilizes seasonal statistics which creates room for better accuracy. TOMCAT is a better model for assessing down-the-drain chemicals at catchment scale than QUAL2EU [11]. Assessment of the orthophosphate concentration in the Thames River was performed using this model [12].

2.3. QUAL2E

This enhanced stream water quality model with uncertainty analysis is a model for conventional pollutants in branching streams and well-mixed lakes and is used by the United States Environmental Protection Agency (USEPA). Produced in 1985, this model has been utilized and improved immensely by USEPA. QUAL2E stems from historical development of nitrogen (N), phosphorus (P) and oxygen (O) models [13] [14]. The first Streeter-Phelps model was its kick-off and it assessed the relationship between BOD and DO. Subsequently, it assessed the temporal and spatial differences in conservative mineral concentrations and water temperature as well as the DO and BOD concentrations, giving rise to the QUAL1 model. Lastly, in addition to the above, nitrogen cycling, phosphorus cycling, algae and several variables were used in creating the QUAL2E model family [15] [16]. Difficulties had been found with considerable use of QUAL2 which needed corrections in the interactions between algal, light and nutrient. QUAL2 was changed to QUAL2E after enhancement by a couple of modifications. With more enhancement, QUAL2E was changed to QUAL2EU. QUAL2EU has explicit features of the QUAL2E version plus the option for uncertainty analysis [13]. This model, which is a 1D steady state model (also operated as a dynamic model) is a useful tool for management of water quality. It is useful in assessing the effect of effluents on in-stream water quality. Also useful in spotting quality characteristics and magnitude of non-point effluents as part of a field-sampling program.

2.3.1. Conceptualization

A 1D conservative advection-dispersion is the basic equation used in explaining the Patten of a pollutant in the river. Its thesis includes: dispersive transport is proportional to concentration gradient, advective transport is within the mean flow and there is a complete mixture of solutes across the cross section. Algae, phosphate, nitrate, DO are assessed in detail while most determinants are simulated as first-order decays. Sediment processes as an oxygen source or as sink for substances are added. Coliforms constituents are non-conservative. A function of first-order decay, which utilizes coliform die-off only is used. Relying on the temperature, the coliforms and non-conservative constituents are assessed as first-order decay. A couple of sub-reaches (computational element) which are divisions of a stream reach equal to finite difference elements, are utilized in QUAL2EU. For every computational element, a heat balance according to temperature, a hydrologic balance according to flow and a material balance according to concentration is documented. Material balance takes into account both dispersive and advective transport. The model utilizes equations of reactions and finite-difference solution of mass transport and also utilizes unique steady-state use if a backward difference numerical scheme, rendering permanent stability to the model [17]. About 15 state variables if major interaction is assessed by each compartment of the model. The 1D conservative advection-dispersion equation is used by QUAL2EU to explain the character of a pollutant in the river (Equation (8)).

$\frac{\partial C}{\partial t}=\frac{\partial \left(A_x{D}_L\frac{\partial C}{\partial x}\right)}{A_x\partial x}+\frac{\partial \left(A_{x}UC\right)}{A_x\partial x}+\frac{\text{d}c}{\text{d}t}+\Delta S$ (8)

$\frac{\text{d}c}{\text{d}t}=Q_i\left(C_i-C\right)/V+\Delta C$ (9)

where D_L is the dispersion coefficient, A_x is the cross-sectional area, C is the concentration of the determinant, U is the mean velocity, t is the time, x is the distance along the element and ΔS is the net concentration influence of outside sources and descends. The alterations happening to separate determinants self-governing of advection, dispersal and external involvements are defined by the term dC/dt and these changes include physical-chemical and biological procedures that happen in the watercourse.

2.3.2. Processes

Figure 1 clarifies in what way a lot of determinants are imitation as first-order BOD decays, but then again phosphate, DO and nitrate are exposed in more aspect with an algal model. Chlorophyll-a is the indicator of algal biomass, utilized by algae. Biomass accumulation is derived from a balance between settling of the

Figure 1. Schematic diagram of interacting water quality state variables in QUAL2EU (ORG-N organic nitrogen, ORG-P organic phosphorous, DIS-P dissolved inorganic phosphorous, CBOD carbonaceous biochemical oxygen demand, SOD sediment oxygen demand, NH₃ ammonia, NO₃ nitrate, NO₂ nitrite, chlorophyll (a) [13].

algae, respiration and growth. The highest manner of growth is nutrient and light limited. The algal is composed of growth by respiration,, photosynthesis and deposition of algal into the river bed settlements.

$\frac{\text{d}L}{\text{d}t}=-K_{1}L-K_{3}L$ (10)

where L is the concentration of the CBODu, K₁ is the rate of oxidation of the CBOD and K₃ is the rate of CBOD loss due to settling. For 5-day CBOD (CBOD5) data, QUAL2EU uses the following conversion for ultimate CBODu (where KCBOD is the CBOD conversion coefficient)

$\text{CBODu}=\frac{\text{CBOD5}}{1-e^{-\left(e_{k\left(\text{CBOD}\right)}\right)}}$ (11)

The primary internal sink of DO in the QUAL2EU includes BOD process, sediment oxygen demand (modeled as a zero order reaction), respiration by algae, and nitrification which includes the oxidation of both ammonia and nitrite. The major source of dissolved oxygen, in addition to that supplied from algal photosynthesis, is atmospheric reaeration. Nine methods are available to calculate the reaeration coefficient in case of free water surface. The DO is modeled by,

$\frac{\text{d}C}{\text{d}t}=\left(C_s-C\right)k_2-\left({\alpha }_3\mu -{\alpha }_4\rho \right)A-K_{1}L-\frac{k_4}{D}-{\alpha }_5\beta \left[{\text{NH}}_4\right]-{\alpha }_6{\beta }_2\left[{\text{NO}}_2\right]$(12)

where C is the concentration of dissolved oxygen, C_s is the saturation concentration, K₂ is the reaeration rate, α₃ is the rate of photosynthetic oxygen production per unit of algal growth, μ is the growth rate of algal biomass (which is affected by nutrient availability, light, temperature, and self-shading), A is the algal biomass (which is directly proportional to chlorophyll-a), α₄ is the rate of respiratory oxygen uptake per unit of algal respiration, K₄ is the rate of the sediment oxygen demand, α₅ is the rate of oxygen utilization per unit of ammonium oxidized during nitrification, α₆ is the rate of oxygen uptake per unit of nitrite oxidized, [NO₂N] is the concentration of nitrite nitrogen, β₂ is the rate coefficient for the oxidation of nitrite nitrogen, [NH₄N] is the concentration of ammonium and β₁ is the rate coefficient parameter for the biological oxidation of ammonium (i.e., nitrification).

The nitrogen cycle is divided into four components: nitrate-nitrogen (NO₃N), ammonium-nitrogen (NH₄N), organic-nitrogen (NO) and nitrate-nitrogen (NO₂N). Mineralization and settling of organic nitrogen, Nitrification ( which consists of ammonia oxidation and oxidation of nitrite into nitrate), algal uptake, renewal from alagal respiration and the sediment are all CONSIDERED in the balance of nitrogen. The inhibitions which occur at low DO, are noted during Nitrification reaction rates. Algse produces organic nitrogen and it is REMOVED by hydrolysis and settling. The equations for nitrogen balance include:

$\frac{\text{d}{N}_0}{\text{d}t}={\alpha }_1\rho A-{\beta }_3{N}_0-{\sigma }_4{N}_0$ (13)

(14)

$\frac{\text{d}\left[{\text{NO}}_{\text{2}}\text{N}\right]}{\text{d}t}={\beta }_1\left[{\text{NH}}_{\text{4}}\text{N}\right]+{\beta }_2\left[{\text{NO}}_{\text{2}}\text{N}\right]$ (15)

And

$\frac{\text{d}\left[{\text{NO}}_{\text{3}}\text{N}\right]}{\text{d}t}={\beta }_2\left[{\text{NO}}_{\text{2}}\text{N}\right]-\left(1-F\right)\alpha \mu A$ (16)

where α₁ is the fraction of the algal biomass that is nitrogen, ρ is the algal respiration rate, β₃ is the rate coefficient parameter for the hydrolysis of organic nitrogen to ammonium, σ₄ is the rate coefficient parameter for the settling of organic nitrogen and F is fraction of algal nitrogen taken from ammonia pool. The phosphorus cycle is similar to, but simpler than the nitrogen cycle, having only two compartments.

Balance of phosphorus requires mineralization of organic phosphorus into inorganic phosphorus, renewal from sediment, uptake and respiration from algae. Due to the decay of organic phosphorus and SIMILAR actions in sediment, dissolved phosphorus is formed and it is used up during the process of photosynthesis. The equations for phosphorus balance are:

$\frac{\text{d}{P}_0}{\text{d}t}={\alpha }_2\rho A-\beta {\rho }_0-{\sigma }_5{p}_0$ (17)

And

$\frac{\text{d}{P}_d}{\text{d}t}={\beta }_4{P}_0+\frac{{\sigma }_2}{D}-{\alpha }_2\mu A$ (18)

where P₀ is the concentration of organic phosphorous, α₂ is the phosphorous content of algae, β₄ is the organic phosphorous decay rate and σ₅ is the organic phosphorous settling rate, P_d is the concentration of inorganic or dissolved phosphorous, and α₂ is the source rate of phosphorous from the sediments. The temperature effect for all first-order reactions used in the model is represented by:

$k\left(t\right)=k\left(20\right){\theta }^{T-20}$ (19)

where k(T) = the reaction rate [/day] at temperature T [˚C] and θ = the temperature coefficient for the reaction.

Every reaction between all expressed state variables depend on temperature, and a correction factor for all coefficients in the source/sink terms are calculated by QUAL2EU with the aid of a Streeter-Phelps type formulation. PERFORMING a heat balance on every element models the temperature. In each compartment, a complete heat balance at the interface of air and water is documented between the total incoming atmospheric radiation, the total incoming short-wave radiation, the loss of heat by conduction to the atmosphere, loss of heat by evaporation and the back radiation from the water surface.

2.3.3. Input Data

Several Windows-based or DOS-based user interface programs are available and the Windows-based version is the most recent and most utilized. Computing performed by the Windows interface user gives a couple of guidelines for choosing inputs [18] [19]. There are three groups of input data: forcing functions, stream/river system and global variables. The global variable group explains the general simulation variables like water quality constituents, units, simulation type and some of the basin’s physical characteristics. The forcing functions are user-specified inputs that DRIVE the modeled system. The input data for the stream/river system, DESCRIBES the stream system into a format that can be processed by the model.

2.3.4. Model Capability

It can simulate about 15 parameters if water quality: temperature, ammonia as N, organic phosphorus as P, nitrite as N, nitrogen as N, nitrate as N, dissolved phosphorus as P, Biochemical oxygen demand, dissolved oxygen, colifirm bacteria, algae as chlorophyll-a, one ARBITRARY constituent solute which is non-conservative, and three constituent solutes which are conservative. The QUAL2EU has in-built tools for uncertainty analysis, to aid in identification of sensitive parameters for a specific use. There are three options for uncertainty analysis: first order error analysis, Monte Carlo simulation and sensitivity analysis. The user is able to assess uncertain input data and the impact of model sensitiveness on model forecast.

2.3.5. Limitations of Model

The QUAL2EU model is unable to convert algal death to CBOD [19] [20] therefore it is inadequate where macrophytes (rooted aquatic plants) have importance [20]. The highest number if junctions, reaches and computational elements are limited and are the only version available of the QUAL2EU model, and therefore is unable to simulate the large river system correctly [21]. In this model, there are certain imposeed dimensional limitations: junction elements, number of reaches, headwater elements and computational elements should be maximum 6, 25, 7, and 250 respectively. A maximum of 25 bis required for both withdrawal and input elements.

2.3.6. Strengths and It’s Application

QUAL2EU can be described as a simple model with COMPREHENSIVE dissolved oxygen, nutrient and algal Dynamics. Its method of use is simple and it has well detailed documentation materials. Materials relating to QUAL2EU are available in textbooks, technical reports and journal papers. The in-built uncertainty analysis tools in QUAL2EU, aids in identification of parameters which are sensitive for a specific application. The QUAL2EU has been used for analysis of systems in the management of water quality, and has been used on a couple of rivers and streams worldwide [22] [23] [24] [25]. North America, Europe, and Asia are places where it has been used.

2.4. QUAL2KW

The QUAL2KW was DEVELOPED in 2002, by Park and Lee [21] after establishing some limitations of the QUAL2E/QUAL2EU. Amongst it’s limitations is the inability to convert algal death to carbonaceous Biochemical oxygen demand [19] [20]. EXPANSION of computational structure plus new constituent interactions, like algal BOD, DO change and denitrification as a result of a fixed plant, are major enhancements of QUAL2K, 2002. The QUAL2K, 2003 model developed by Chapra and Pelletier, was modified to the QUAL2KW. The purpose of QUAL2K and QUAL2KW is to represent the modernized versions of QUAL2E. Some processes which are not encompassed in QUAL2K 2002, and QUAL2K 2003 are included in QUAL2KW. HOWEVER, it is worthy of note that there is no relationship between QUAL2K and QUAL2KW with QUAL2K, 2002.

QUAL2KW is a 1D, steady state (particularly, steady flow model as flows are presumed to be steady, but the heat budget and diel water quality kinetics are used to calculate water qualities dynamically) water quality model in stream water and is utilized in the Microsoft Windows environment. It is available and well detailed in http://www.ecy.wa.gov/.

There are many new elements in QUAL2KW [26]. It utilizes reaches which are unequally spaced. Input for any reach can be derived from multiple loadings and abstractions. It makes use of two types of carbonaceous CBOD to REPRESENT organic carbon. The types include: rapidly oxidizing form (fast CBOD) and slowly oxidizing form (slow CBOD). Also, there is simulation of non-living particulate organic matter (detritus).

It decreases oxidation reactions to zero at low oxygen level, therefore accommodating anoxia. Denitrification is modeled as a firs-order reaction, which rises at low oxygen concentrations. Instead if beinng prescribed, sediment water fluxes of nutrients and dissolved oxygen are simulated. The model represents bottom algae EXPLICITLY. Inorganic acids, algae and detritus are calculated as functions of light extinction. There is simulation of river pH, alkalinity and total inorganic carbon.

QUAL2KW encompasses new options and processes, in addition to the processes in QUAL2K (2003). It assesses water exchange between hyporheic zone, surface water column and sediment pore-water quality. A generic pathogen is simulated as a function of settling, temperature and light. It consists of a genetic algorithm to calibrate parameters of kinetic rate automatically. A well detailed documentation and theory for QUAL2KW can be found at Pelletier and Chapra [26].

2.4.1. Conceptualization

The main branch of a river is seen as a group of reaches (unequal or equal lengths). Tributaries are REPRESENTED as point sources. For a constituent in the water column of a reach i (Figure 2), the equation for the general mass balance is (Equation (20)):

Figure 2. Mass balance in a reach segment i in QUAL2Kw [26].

$\begin{array}{c}\frac{\text{d}{c}_i}{\text{d}t}=\frac{Q_{i-1}}{V_i}{c}_{i-1}-\frac{Q_i}{V_i}{c}_i-\frac{Q_{ab,i}}{V_i}{c}_i+\frac{{E^{\prime }}_{i-1}}{V_i}\left(c_{i-1}-c_i\right)+\frac{{E^{\prime }}_i}{V_i}\left(c_{i+1}-c_i\right)\\ +\frac{W_i}{V_i}+S_i+\frac{{E^{\prime }}_{hyp,i}}{V_i}\left(c_{2,i}-c_i\right)\end{array}$ (20)

The general mass balance for a constituent concentration in the hyporheic sediment zone of a reach (c_2,i) (for all but the heterotrophic bacteria biofilm) is written as (Equation (21))

$\frac{\text{d}{c}_{2,i}}{\text{d}t}=S_{2,i}+\frac{{E^{\prime }}_{hyp,i}}{V_{2,i}}\left(c_i-c_{2,i}\right)$ (21)

For bottom algae, the transport, heterotrophic bacteria and loading terms are omitted. For heterotrophic bacteria, the transport and loading terms are omitted. The mass balance equations for bottom algae (Equations (22a)-(22c)) and bacteria (Equation (23), Level 2 option) are:

$\frac{\text{d}{a}_{b,i}}{\text{d}t}=S_{b,i}$ (22a)

$\frac{\text{d}I{N}_b}{\text{d}t}=S_{bN,i}$ (22b)

$\frac{\text{d}I{P}_b}{\text{d}t}=S_{bP,i}$ (22c)

$\frac{\text{d}{a}_{h,i}}{\text{d}t}=S_{ah,i}$ (23)

where, S_i is the sources and sinks of constituent i due to reactions and mass transfers (mg/L/d), a_p is the phytoplankton concentration (mgA/m³), a_b is the bottom algae concentration (gD/m²), a_h is the biofilm of attached heterotrophic bacteria in the hyporheic sediment zone, Q_i is the outflow from reach i, Q_i–1 is the inflow from the upstream reach i – 1, Q_ab,i is the total flow abstractions from the reach i, W_i is the external loading of the constituent to reach i (mg/d), ${E^{\prime }}_{i-1}$, ${E^{\prime }}_i$ are bulk dispersion coefficients between reaches i – 1 and i and i and i + 1 (m³/day), respectively, ${E^{\prime }}_{hyp,j}$ is the bulk dispersion coefficients between hyporheic zone and reach i (m³/day), V_i is the volumes of reach i (m³), t is time (day), V_2,i (=Φ_s,iA_st,iH_2,i/100) is the volume of pore water in the hyporheic sediment zone (m³), A_st,i is the surface area of the reach (m²), H_2,i = the thickness of the hyporheic zone (cm). Φ_s,i is the porosity of the hyporheic sediment zone. c_i is the concentration in the surface water in reach i (mg/L), c_i−1 is the concentration in the upstream reach i − 1 (mg/L), c_2,i is the concentration in hyporheic sediment zone (mg/L), IN_b is the intracellular nitrogen concentration in bottom algae (mgN/m²), IP_b is the intracellular phosphorus concentration in bottom algae (mgP/m²), S_b,i is the sources and sinks of bottom algae biomass due to reactions (gD/m²/day), S_bN,i is the sources and sinks of bottom algae nitrogen due to reactions (mgN/m²/day), S_bP,i is the sources and sinks of bottom algae phosphorus due to reactions (mgP/m²/day), S_ah,i is the sources and sinks of heterotrophic bacteria in the hyporheic sediment zone due to reactions (gD/m²/day) and S_2,i is the sources and sinks of the constituent in the hyporheic sediment zone due to reactions. Sources and sinks of the water quality constituent due to reactions and mass transfer mechanisms are depicted in the algal model (Figure 3). The model computes sediment oxygen demand, in which sediment-water fluxes are computed based on the downward flux of particulate organic matter from the overlying water. Organic carbon, nitrogen, and phosphorus delivered to the anaerobic sediments are transformed by mineralization into dissolved methane, ammonium and inorganic phosphorus, which are transported to the aerobic layer oxidizing some of them, causing sediment oxygen demand.

The algal model shows the sources and sinks of the water quality constituent, which are due to reactions and mass transfer mechanisms (Figure 3). The model documents sediment oxygen demand, whereby the sediment-water fluxes are

Figure 3. Schematic diagram of interacting water quality state variables in QUAL2Kw (a_b bottom algae, a_p phytoplankton, m_o detritus, c_s slow CBOD, c_f fast CBOD, c_T total inorganic carbon, O oxygen, n_O organic nitrogen, n_a ammonia nitrogen, n_n nitrite and nitrate nitrogen, p_o organic phosphorous and p_i inorganic phosphorous) [26].

documented according to the downward flux of particulate organic matter from the overlying water. Delivery of nitrogen, phosphorus and organic carbon to the anaerobic sediments are converted into inorganic phosphorus, dissolved methane and ammonium by mineralization, which are taken to the aerobic layer and oxidizing some of them, leading to sediment oxygen demand.

2.4.2. Processes

Several determinants are simulated as first-order decays, though nitrate, phosphate and DO are explained in more detail. The algal model for water quality variables interaction, is explained by Figure 3. The model comprises accumulation of detritus (m_o), due to death of plants ( bottom algae and phytoplankton), dissolution and settling (Equation (24)). There is a gradual increase in reacting CBOD (CS) resulting from dissolution of detritus and is lost through oxidation and hydrolysis (Equation (25)). Changes resulting from low oxygen F_oxcf (denoted by F), is modeled by Equations (26a)-(26c) with the oxygen inhibition parameter K_socf (denoted by K). Hydrolysis of slow reacting CBOD gives rise to the fast reacting CBOD (c_f) and is lost through oxidation and denitrification (Equation (27)). Below are the equations representing C_s & C_f increase and loss, detritus (m_o) accumulation, and attenuation reactions

$S_{mo}=r_{da}{k}_{dp}{a}_p+\frac{k_{db}{a}_b}{H}-k_{dt}{m}_0-\frac{v_{dt}{m}_0}{H}$ (24)

$S_{cs}=r_{od}{k}_{di}{M}_0-K_{h_C}{C}_s-F_{0xcf}\cdot k_{dCs}{c}_S$ (25)

$F=\frac{0}{k+0}$ (half saturation) (26a)

$F=\left[1-\exp \left(-k\cdot o\right)\right]$ (exp) (26b)

$F=\frac{0^2}{k+0^2}$ (26c)

$S_{cf}=K_{hc}{C}_s-F_{oxcf}{k}_{dC}{C}_f-r_{ondn}\left(1-F_{oxndn}\right)k_{dn}{n}_n+J_{c{H}_{4d}}\frac{A_{sti}}{v_i}$ (27)

where a_p = phytoplankton concentration (mgA/m³), a_b = bottom algae concentration (mgD/m³), A = chlorophyll a, D = detritus, r_da = ratio of dry weight of detritus to chlorophyll a (gD/mgA), r_od =ratio of CBOD generated to detritus dissolution (mgO₂/mgD), r_ondn = ratio of oxygen equivalents lost per nitrate nitrogen that is denitrified (gO₂/g N), k_dp = phytoplankton death rate (/day), k_db = bottom algae death rate (/day), k_dt = detritus dissolution rate (/day), k_hc = slow CBOD hydrolysis rate (/day), k_dc = fast CBOD oxidation rate (/d), k_dcs = slow CBOD oxidation rate (/day), k_dn = denitrification rate of nitrate nitrogen (/day), n_n = nitrate + nitrite nitrogen (mg/L), F_oxdn = effect of low oxygen on denitrification (dimensionless), F_oxcf = attenuation due to low oxygen on fast CBOD oxidation (dimensionless), O = oxygen concentration (mg/L), v_dt = detritus settling velocity (m/day), J_CH4,d =the sediment flux of dissolved methane in oxygen equivalents (gO₂/m²/day) and A_st,i =the surface area of the reach (m²).

The concentration of organic nitrogen (n_o) rises as a result of plants death and is lost through hydrolysis and settling (Equation (28)). There is an increase in ammonia nitrogen(n_a) from hydrolysis of organic nitrogen and phytoplankton respiration and lost via photosynthesis and Nitrification (Equations (29a)-(29c)). There is a rise in nitrogen (n_n) from the Nitrification of ammonia and us list via photosynthesis and denitrification (Equation (30)). Change in Nitrification resulting from reduced oxygen, F_oxna (depicted as F) is represented by Equations (26a)-(26c) with k_sona (depicted as K), which is an oxygen dependency parameter. The equations which involve nitrate nitrogen (n_n), organic nitrogen (n_o), and ammonia nitrogen (n_a) are stated below.

$S_{no}=r_{na}\left(k_{dp}{a}_p\right)+q_N\left(\frac{k_{bd}{a}_b}{H}\right)-k_{h_n}{n}_0-\frac{v_{0}n}{H}{n}_0$ (28)

$\begin{array}{c}{S}_{na}=k_{hn}{n}_0-F_{0xna}{k}_{na}{n}_a+v_{na}\left(F_0{x}_P{K}_{{\gamma }_P}\right)a_p-r_{\eta a}{F}_{ap}\left[k{g}_P{\varphi }_{Np}{\varphi }_{Lp}\right]a_p+J_{N{H}_4}\\ -\left[p_{ab}\text{botalguptakeN}\left(\text{NupWcFrac}\right)\right]\frac{A_{sti}}{v_i}\end{array}$ (29)

$P_{ap}=\frac{n_a{n}_n}{\left(k_{hnxp}+n_a\right)\left(k_{hnxp}+n_n\right)}+\frac{n_a{k}_{hnxp}}{\left(n_a+n_n\right)\left(k_{hnxp}+n_n\right)}$ (29a)

$P_{ab}=\frac{n_a{n}_n}{\left(k_{hnxb}+n_a\right)\left(k_{hnxb}+n_n\right)}+\frac{n_a{k}_{hnxb}}{\left(n_a+n_n\right)\left(k_{hnxb}+n_n\right)}$ (29b)

$\begin{array}{c}{S}_{nn}=F_{0xna}{k}_{na}{n}_a-\left(1-F_{0xna}\right)k_{dn}{n}_n-r_{\eta a}{F}_{ap}\left[k{g}_P{\varphi }_{Np}{\varphi }_{Lp}\right]a_p+J_{N{O}_3}\\ -\left[\left(1-p_{ab}\right)\text{botalguptakeN}\left(\text{NupWcFrac}\right)\right]\frac{A_{sti}}{v_i}\end{array}$ (30)

where r_na = ratio of nitrogen to chlorophyll a (mgN/mgA), khn = organic nitrogen hydrolysis rate (1/day), k_na = nitrification rate for ammonia nitrogen (1/day), k_rp = phytoplankton respiration rate (1/day), k_gp = maximum photosynthesis rate at temperature (/d), k_dp = phytoplankton death rate (/day), F_oxna = attenuation due to low oxygen on ammonia nitrification (dimensionless), F_oxp = attenuation due to low oxygen of phytoplankton respiration, v_on = organic nitrogen settling velocity (m/d), q_N = actual cell quotas of nitrogen (mgN/gD), NUpWCfrac = fraction of N uptake from the water column by bottom plants. $J_{N{H}_4}$ = sediment flux of ammonia (mgN/m²/day), BotAlgUptakeN = uptake rate for nitrogen in bottom algae (mgN/m²/day) as defined in Equation (35f), $J_{N{O}_3}$ is the sediment flux of nitrate (mgN/m²/day), A_st,i = the surface area of the reach (m²), P_ap = preferences for ammonium as nitrogen source for phytoplankton, P_ab = preferences for ammonium as nitrogen source for bottom algae, Φ_Np = phytoplankton nutrient attenuation factor (dimensionless), Φ_Lp = phytoplankton light attenuation factor (dimensionless), k_hnxp = preferences coefficient of phytoplankton for ammonium (mgN/m³), k_hnxb = preferences coefficient of bottom algae for ammonium (mgN/m³) and k_nn = temperature-dependent nitrification rate for ammonia nitrogen (1/day).

Inorganic phosphorus (p_i) rises as a result of phytoplankton respiration and the hydrolysis of organic phosphorus, and is lost via plant photosynthesis, settling and uptake of bottom algae (Equation (32)). There is an increase in organic phosphorus (p_o), resulting from plant death and is lost through hydrolysis and settling (Equation (31)). Inorganic suspended solids (m_i) are lost through settling (Equation (33)). Loss and increase of phosphorus are represented by the equations,

$S_{p_0}=r_{Pa}\left(k_d{p}^{a}P\right)+\left(\frac{k_{bd}{a}_b}{H}\right)q_p-k_{np}{p}_0-\frac{v_{0}p}{H}{p}_0$ (31)

$\begin{array}{c}{S}_{pi}=\text{DOPHydr}+r_{pa}\text{PhytoResp}-r_{pa}\text{PhytoPhoto}-\text{IPSettl}\\ +\left[J_{\text{PO4}}-\text{BotAlgUptakeP}\left(\text{PUpWCfrac}\right)\right]\frac{A_{st,i}}{V_i}\end{array}$ (32)

$S_{mi}=\frac{v_i}{H}{m}_i$ (33)

where r_pa = ratio of phosphorus to chlorophyll a (mgP/mgA), k_hp = organic phosphorus hydrolysis rate (/day), k_dp = phytoplankton death rate (/day), q_P = actual cell quotas of phosphorous (mgP/gD), v_op = organic phosphorus settling velocity (m/day), v_ip = inorganic phosphorus settling velocity (m/day), v_i = inorganic suspended solids settling velocity (m/day), PUpWCfrac = fraction of P uptake from the water column by bottom plants, J_PO4 = the sediment flux of inorganic P (mgP/m²/day) and BotAlgUptakeP = uptake rate for phosphorous in bottom algae (mgP/m²/day) as defined in Equation (35g).

Phytoplankton (a_p) increases due to photosynthesis and is lost via respiration, death and settling (Equations (34a)-(34e)). Michaelis-Menten equations (Equation (34b)) are used to compute the nutrient attenuation factor φ_Np. Three models (Equations (34c)-(34e)): half-saturation, Smith’s function and Steele’s equations are combined with Beer-Lambert law to characterize the impact of light on phytoplankton photosynthesis (to estimate’ Lp) as,

$S_{ap}=k{g}_P{\varphi }_{N_P}{\varphi }_{LP}{a}_p-F_{0xp}{k}_{\gamma p}{a}_p-k_{dP}{a}_p-\frac{v_{a_P}}{H}{a}_P$ (34a)

${\varphi }_{Np}=\min \left(\frac{n_a+n_n}{k_{sNp}+n_a+n_n},\frac{p_i}{k_{sPp}+p_i},\frac{\left[{\text{H}}_{\text{2}}{\text{CO}}_3^{*}\right]+\left[{\text{HCO}}_3^{-}\right]}{k_{sCp}+\left[{\text{H}}_{\text{2}}{\text{CO}}_3^{*}\right]+\left[{\text{HCO}}_3^{-}\right]}\right)$ (34b)

${\varphi }_{Lp}=\frac{1}{k_{e}H}\ln \left(\frac{K_{Lp}+I\left(0\right)}{K_{Lp}+I\left(0\right){\text{e}}^{-k_{e}H}}\right)$ (34c)

${\varphi }_{Lp}=\frac{I\left(0\right){\text{e}}^{-k_{e}H}}{\sqrt{K_{Lb}^2+{\left(I\left(0\right){\text{e}}^{-k_{e}H}\right)}^2}}$ (smith 1932) (34d)

${\varphi }_{Lb}=\frac{I\left(0\right){\text{e}}^{-k_{e}H}}{K_{Lb}}{\text{e}}^{1+\frac{I\left(0\right){\text{e}}^{-k_{e}H}}{K_{Lb}}}$ (steel 1962) (34e)

where v_ap = phytoplankton settling velocity (m/day), k_sNp = nitrogen half-saturation constant for phytoplankton (μg N/L), k_sPp = phosphorus half-saturation constant for phytoplankton (μg N/L), k_sCp = inorganic carbon half-saturation constant for phytoplankton (mole/L), K_Lp = light constant for phytoplankton (ly/day), and K_e = extinction coefficient, $\left[{\text{H}}_{\text{2}}{\text{CO}}_3^{*}\right]$ = sum of dissolved CO₂ and HCO₃ (mg/m³/day).

Bottom algae (a_b) increase due to photosynthesis and are lost via respiration and death (Equation (35a)). Attenuation due to low oxygen F_oxb is modeled by Equations (26a)-(26c) with the oxygen inhibition parameter k_sob. A logistic model (Equation (35b)) is used to estimate space limitation attenuation factor φ_Sb. The effect of nutrient limitation on bottom algae photosynthesis is modeled using a formulation, according to which the photosynthesis rate is dependent on intracellular nutrients (Equations (35c)-(35g)). The three models Equations (36a)-(36c), as in phytoplankton, are combined with Beer-Lambert law and integrated to yield light attenuation factor φ_Lb as given by equations,

$\begin{array}{c}{S}_o=r_{oa}\left(\text{PhytoPhoto-PhytoResp}\right)+r_{od}\left(\text{BotAlgPhoto-BotAlgResp}\right)\frac{A_{st,i}}{V_i}\\ -r_{oc}\text{FastCOxid}-r_{oc}\text{SlowCOxid}-r_{on}\text{NH4Nitr}+\text{OxReaer}\\ -\text{CODoxid}-\text{SOD}\frac{A_{st,i}}{V_i}\end{array}$ (35a)

where bottom algae photosynthesis (BotALgPhoto) can be modeled by two options: 1) BotAlgPhoto 1/4 CgbfNbfLbfSbap (zero-order growth rate), where C_gb is the zero-order temperature/dependent maximum photosynthesis rate (gD/m²/day) and 2) BotAlgPhoto 1/4 CgbfNbfLbfSbab (first-order growth rate), where C_gb = the first-order temperature-dependent maximum photosynthesis rate (day⁻¹)

${\varphi }_{Sb}=1-\frac{a_b}{a_{b,\max }}$ (35b)

${\varphi }_{N_b}=\min \left(1-\frac{qN}{9N}\right)\left(1-\frac{q_{0}p}{qp}\right)\left[\frac{\left[{\text{H}}_{\text{2}}{\text{CO}}_3^{*}\right]+\left[{\text{HCO}}_3^{-}\right]}{k_5{c}_b+\left[{\text{H}}_{\text{2}}{\text{CO}}_3^{*}\right]+\left[{\text{HCO}}_3^{-}\right]}\right]$ (35c)

$q_N=\frac{I{N}_b}{a_b}$ (35d)

$q_P=\frac{I{P}_b}{a_b}$ (35e)

$\text{BotAlgUptakeN}={\rho }_{mN}\frac{n_a+n_n}{k_{sNb}+n_a+n_n}\frac{K_{qN}}{K_{qN}+\left(q_N-q_{0N}\right)}{a}_b$ (35f)

$\text{BotAlgUptakeP}={\rho }_{mP}\frac{p_i}{k_{sPb}+p_i}\frac{K_{qP}}{K_{qP}+\left(q_P-q_{0P}\right)}{a}_b$ (35g)

The change in intracellular nitrogen and phosphorus in bottom algal cells are calculated respectively from Equations (35h) and (35i) as,

$S_{bN}=\text{BotAlgUptakeN}-q_N\text{BotAlgDeath}-q_N{k}_{excb}\left(T\right)a_b$ (35h)

$S_{bP}=\text{BotAlgUptakeP}-q_P\text{BotAlgDeath}-q_P{k}_{excb}\left(T\right)a_b$ (35i)

Following formulas are used for the bottom algae light attenuation coefficient:

${\varphi }_{Lb}=\frac{I\left(0\right){\text{e}}^{-k_{e}H}}{K_{Lb}+I\left(0\right){\text{e}}^{-k_{e}H}}$ (36a)

${\varphi }_{Lp}=\frac{I\left(0\right){\text{e}}^{-k_{e}H}}{\sqrt{K_{Lb}^2+{\left(I\left(0\right){\text{e}}^{-k_{e}H}\right)}^2}}$ (36b)

${\varphi }_{Lb}=\frac{I\left(0\right){\text{e}}^{-k_{e}H}}{K_{Lb}}{\text{e}}^{1+\frac{I\left(0\right){\text{e}}^{-k_{e}H}}{K_{Lb}}}$ (36c)

where k_rb1 = temperature-dependent bottom algae basal respiration rate [/day], k_rb2 = bottom algae photo-respiration rate constant [dimensionless], F_oxb = attenuation due to low oxygen on bottom algae respiration, k_rb Bottom algae respiration rate (/day), Φ_Lb = bottom algae light attenuation factor (dimensionless), Φ_Nb = bottom algae nutrient attenuation factor (dimensionless), Φ_Sp = bottom algae space limitation factor (dimensionless), q_N = cell quotas of nitrogen (mgN/gD), q_P = cell quotas of phosphorus (mgP/gD), q_0P = minimum cell quotas of phosphorus (mgP/gD), K_Lb = light constant for bottom algae (Ly/day), IN_b = intracellular nitrogen concentration (mgN/m²) and IP_b = intracellular phosphorus concentration (mgP/m²). k_sNb = nitrogen half-saturation constant for bottom algae (μg N/L), k_sPb = phosphorus half-saturation constant for bottom algae (μg N/L), k_scb = inorganic carbon halfsaturation constant for the bottom algae (mole/L), K_qN = half-saturation constants for intracellular nitrogen (mgN/gD), K_qP = half-saturation constants for intracellular phosphorus (mgP/gD), ρ_mN = maximum uptake rates for nitrogen in bottom algae (mgN/gD/day), ρ_mP = maximum uptake rates phosphorus in bottom algae (mgP/gD/day) and a_b,max = carrying capacity (gA/m²).

Pathogens (x) are lost due to death and settling (Equation (37)). Pathogens death is due to natural die-off and light and their death is modeled as a first-order temperature dependent decay. The death rate due to light is based on the Beer–Lambert law. The representing equation for pathogens is

$Sx=k_{dx}\left(T\right)x+{\alpha }_{path}\frac{I\left(0\right)/24}{k_{e}H}\left(1-{\text{e}}^{-k_{e}H}\right)x$ (37)

where x = pathogens concentrations (cfu/100 mL), k_dx = temperature dependent pathogen die-off rate (/day), H = water depth (m), v_x = pathogen settling velocity (m/day), I(0) = light intensity just below the water surface (Ly/day), α_path = a light efficiency factor (dimensionless), ke is the light extinction coefficient (/day). The model simulates a generic constituent and has option of either assuming no interaction with other state variables, or it can be treated as a non-carbonaceous non-nitrogenous chemical oxygen demand [COD]. The generic constituent is subject to firstorder decay and settling as shown (Equation (38)).

$S_{\text{COD}}=-K_{\text{COD}}\left[\text{COD}\right]-\left[\frac{v_{\text{COD}}}{H}\right]\left(\text{COD}\right)$ (38)

where k_COD = COD oxidation rate (/day), v_COD = settling velocity of COD (m/day). Dissolved oxygen (O) increases due to photosynthesis and reaeration (Equations (39a)-(39c)). It is lost via CBOD oxidations, nitrification, plant respiration, and oxidation of COD. The gain or loss of oxygen from atmosphere depends upon whether the water is under saturated or over saturated. The saturation concentration depends upon local temperature and elevation above mean sea level (Equation (39c)). The reaeration rate ka is calculated using eight equations available in the model which include: O’Connor-Dobbins, Churchill, Owens-Gibbs, Tsivoglou-Neal, Thackston-Dawson, USGS (pool-riffle), USGS (Channel control) and internal calculation. The effects of control structures are also included in the model. The equations involving dissolved oxygen (O) are:

$\begin{array}{c}{S}_o=r_{oa}\left(\text{PhytoPhoto-PhytoResp}\right)+r_{od}\left(\text{BotAlgPhoto-BotAlgResp}\right)\frac{A_{st,i}}{V_i}\\ -r_{oc}\text{FastCOxid}-r_{oc}\text{SlowCOxid}-r_{on}\text{NH4Nitr}+\text{OxReaer}\\ -\text{CODoxid}-\text{SOD}\frac{A_{st,i}}{V_i}\end{array}$ (39a)

$\text{Reaeration}=k_a\left[\left({\text{e}}^{\ln 0s}\left(1-0.0001148\cdot elev\right)-o\right)\right]$ (39b)

$\begin{array}{c}\ln o_s\left(T,0\right)=-139.34411+\frac{1.575701\times {10}^5}{T_a}-\frac{6.642308\times {10}^7}{T_a^2}\\ +\frac{1.243800\times {10}^{10}}{T_a^3}-\frac{8.621949\times {10}^{11}}{T_a^4}\end{array}$ (39c)

where k_a = reaeration rate (1/day), O_s = saturation concentration of dissolved oxygen (mgO₂/L), T_a = absolute temperature (˚K), elev = elevation of the area (m), r_od = ratio of oxygen consumed to detritus (mgO₂/mgD) during bottom algae respiration, r_on = ratio of oxygen consumed to nitrogen during nitrification (mgO₂/mgN), r_oa = ratio of oxygen generated to phytoplankton growth (mgO₂/mgA), r_oc = ratio of oxygen consumed during CBOD oxidation (mgO₂/mgO₂), Phytophoto = phytoplankton photosynthesis (gO₂/m²/day), PhytoResp = Phytoplankton respiration (gO₂/m²/day), BotAlgPhoto = bottom algae photosynthesis (gO₂/m²/day), BotAlResp = bottom algae respiration (gO₂/m²/day), FastCOxid = fast CBOD oxidation (gO₂/m²/day), SlowCOxid = slow CBOD oxidation (gO₂/m²/day), NH4nitr = ammonium nitrification (gO₂/m²/day), Reaeration = (gO₂/m²/day) and CODoxid = oxidation of non-carbonaceous non-nitrogenous chemical oxygen demand (gO₂/m²/day). The temperature effect for all first-order reactions used in the model is represented by:

$k\left(T\right)=k\left(20\right){\theta }^{T-20}$ (40)

where k(T) = the reaction rate [/day] at temperature T [˚C] and θ = the temperature coefficient for the reaction. The pH is estimated by equilibrium, mass balance and electro-neutrality equations as follows: (Equations (41)-(43))

$K_1=\frac{\left[{\text{HCO}}_{\text{3}}^{-}\right]\left[{\text{H}}^{+}\right]}{\left[{\text{H}}_{\text{2}}{\text{CO}}_{\text{3}}^{\text{*}}\right]}$ (41a)

$K_2=\frac{\left[{\text{CO}}_{\text{3}}^{2-}\right]\left[{\text{H}}^{+}\right]}{\left[{\text{HCO}}_{\text{3}}^{-}\right]}$ (41b)

$K_w=\left[{\text{H}}^{+}\right]\left[{\text{OH}}^{-}\right]$ (41c)

$c_T=\left[{\text{H}}_{\text{2}}{\text{CO}}_{\text{3}}^{\text{*}}\right]+\left[{\text{HCO}}_{\text{3}}^{-}\right]+\left[{\text{CO}}_{\text{3}}^{2-}\right]$ (42)

$\text{pH}=-{\log }_{10}\left[{\text{H}}^{+}\right]$ (43)

where K₁, K₂, and K_w are acidity constants, Alk = alkalinity [eq L-1], H₂CO₃* = the sum of dissolved carbon dioxide and carbonic acid, ${\text{HCO}}_{\text{3}}^{-}$ = bicarbonate ion, ${\text{CO}}_{\text{3}}^{2-}$ = carbonate ion, H⁺ = hydronium ion, OH⁻ = hydroxyl ion, and c_T = total inorganic carbon concentration [mol/L]. The model accounts changes in alkalinity considering plant photosynthesis, respiration, nutrients hydrolysis, nitrification, denitrification, and bottom algae uptakes. It estimates total inorganic carbon, considering fast carbon oxidation, plant respiration, plant photosynthesis, and atmospheric aeration.

Temperature is modeled by performing a heat balance on each element i (Equations (44a)-(44b)). The heat balance takes into account heat transfers from adjacent reaches, loads, abstractions, the atmosphere, and the sediments and includes the influences of radiation, convection, and evaporation. The surface heat exchange J_h is modeled as a combination of five processes: solar shortwave radiation at the water surface I(0), atmospheric long wave radiation J_an, long-wave back radiation from the water J_br, conduction J_c, and evaporation J_e (fluxes are expressed as cal/cm²/day) as represent by Equations (44a)-(44b)

$J_h=I\left(o\right)+J_{un}-J_{br}-J_c-J_e$ (44a)

$\begin{array}{c}\frac{\text{d}{T}_i}{\text{d}t}=\frac{Q_{i-1}}{V_i}{T}_{i-1}-\frac{Q_i}{V_i}{T}_i-\frac{Q_{ab,i}}{V_i}{T}_i+\frac{{E^{\prime }}_{i-1}}{V_i}\left(T_{i-1}-T_i\right)+\frac{{E^{\prime }}_i}{V_i}\left(T_{i+1}-T_i\right)\\ +\frac{W_{h,i}}{{\rho }_w{C}_{pw}{V}_i}\left(\frac{{\text{m}}^3}{{10}^6{\text{cm}}^3}\right)+\frac{J_{h,i}}{{\rho }_w{C}_{pw}{H}_i}\left(\frac{\text{m}}{100\text{cm}}\right)+\frac{J_{s,i}}{{\rho }_w{C}_{pw}{H}_i}\left(\frac{\text{m}}{100\text{cm}}\right)\end{array}$ (44b)

where T_i = temperature in reach i (˚C), t = time (day), ${E^{\prime }}_i$ = the bulk dispersion coefficient between reaches i and i + 1 (m³/day), W_h,I = the net heat load from point and non-point sources into reach i (cal/day), ρ_w = the density of water (g/cm³), C_pw = the specific heat of water (cal/g˚C), J_h,i = the air-water heat flux (cal/cm²/day), J_s,i = the sediment-water heat flux (cal/cm²/day), J_s,i = sediment-water heat flux (cal/cm²/day) and ρ_w = density of water (g/cm³). A complete discussion of the model theory is described by Pelletier and Chapra [26].

For auto-calibration, the model uses genetic algorithm to maximize the goodness of fit of the model results compared with measured data by adjusting a large number of parameters [27] [28]. The fitness is determined as the reciprocal of the weighted average of the normalized root mean squared error of the difference between the model predictions and the observed data for water quality constituents. The genetic algorithm maximizes the fitness function f(x) as (Equation (45)),

$f\left(x\right)=\left[{\displaystyle {\sum }_{i=1}^m{w}_i}\right]\left[{\displaystyle {\sum }_{i=1}^m\left[\frac{\frac{{\displaystyle {\sum }_{j=1}^m{O}_{ij}}}{m}}{{\left[\frac{{{\displaystyle \sum }}^{\text{}}{\left(P_{ij}-O_{ij}\right)}^2}{2}\right]}^{0.5}}\right]}\right]$ (45)

where O_i,j = observed values, P_i,j = predicted values, m = number of pairs of predicted and observed values, w_i = weighting factors, and n = number of different state variables included in the reciprocal of the weighted normalized root mean squared error.

2.4.3. Input Data

Date, location, boundary conditions of flow, numerical integration control options, headwater boundary flow and concentrations, concentration for tributary point sources and diffuse sources, air temperature, hydraulic geometry (rating curve or Manning’s equation inputs for depth and velocity), shade, light attenuation parameters, reach segment lengths, options of solar radiation, cloud cover, dew point temperature, long wave radiation, wind speed, evaporation, parameters to control the genetic algorithm for optional automatic calibration of water quality kinetics rates and constants, parameters for water quality kinetics rates and constants.

2.4.4. Model Capability

pH, temperature, conductivity, fast reacting CBOD, slow reacting CBOD, inorganic suspended solids, dissolved oxygen, organic phosphorus, inorganic phosphorus, ammonia nitrogen, organic nitrogen, nitrate nitrogen, pathogen, total organic carbon, bottom algae ( periphyton) nitrogen, bottom algae (periphyton) phosphorus, bottom algae (periphyton) biomass, detritus, phytoplankton, total inorganic carbon, alkalinity are all simulated by this model. It simulates only the main stream of a river and not its branches. Currently, it does not consistent an uncertainty component. It is a steady-flow, 1D model and is unable to assess variable flow.

2.4.5. Strengths and Its Application

The QUAL2KW is able to derive carbonaceous Biochemical Oxygen demand, via the conversion of algal death [19] [20]. Therefore, where macrophyte (rooted aquatic plants) serve a unique purpose, this model is suitable. It is well detailed and available. It is made up of automatic calibration system. This model has been used on the south Umpqua river, Oregon, USA to assess the dissolved oxygen [29], and used to assess dissolved oxygen in the Bagmati river, Nepal [30], and it has been used to calculate total maximum daily load studies for dissolved oxygen, temperature, pH and nutrients in the Wenatchee River, Washington state [31].

2.5. WASP7

WASP7 [32] is a modernized version of the original WASP [33] [34] [35] [36], which is found free in the United States Environmental Protection Agency’s website http://www.epa.gov/athens/wwqtsc/. The development of this model began in the 1970’s [37]. The WASP can be seen as a dynamic compartment-modeling program for aquatic systems, the water column and underlying benthos both inclusive. The following processes change with time, and they include: dispersion, point and diffuse mass loading, advection-dispersion and boundary exchange are simulated in the basic program. 1, 2, and 3 dimensional systems and different types of pollutant can be assessed using WASP7. It is useful in assessing different problems of water quality in various water bodies such as, streams, rivers, coastal waters, ponds, estuaries, lakes, and reservoirs. It is associated with sediment transport models and hydrodynamic models (like EFDC, DYNHYD; environmental fluid Dynamics code), which provides salinity, flows, sediment fluxes, temperature and depth velocities. It contains

· High speed WASP7 sub-model processors.

· Easily accessible Windows-based interface.

· A pre-processor to aid in processing of data into a format which can be used in WASP7.

· Graphical postprocessor to enable viewing of the WASP7 results.

The old versions of WASP have two general kinetic modules [33] : EUTRO for conventional water quality and TOXU for toxicants. Advanced sub-models of EUTRO (named periphyton), HEAT and MECURY were added for particular modeling needs. This provided a medium for Simulation of periphyton (attached bottom algae) with nutrients uptake. Routines in the QUAL2K model contain algorithms which assess changes in concentration of periphyton and detritus [27]. The kinetic sub-models TOXI and EUTRO can aid in simulation of two major types of water quality problems: toxic pollution (with sediments, organic chemicals and metals) and conventional pollution (with Biochemical oxygen demand, dissolved oxygen, eutrophication and nutrients). The EUTRO can be utilized at different complex levels to simulate the variables and interactions. About three types of particulate material and about three chemicals have been simulated by TOXI. The chemicals may be associated with reaction yields or may not. Any of the chemicals being simulated, may present in five ways:

1) Doubly charged cations.

2) Singly charged cations.

3) Doubly charged anions.

4) Singly charged anions.

5) Neutral molecules.

The ionic or neutral molecules can also occur in five ways:

1) Dissolved phase.

2) Sorbed to dissolved organic carbon (DOC).

3) Sorbed to the three different solids types.

Therefore, 25 forms of each chemical can be modeled in TOXI. Ambrose et Al [32] and Wool et al. [33] contain a well detailed documentation.

2.5.1. Conceptualization

This model is based on compartmentalization principle. For each equation, a mass balance equation is stated out. Within each compartment, there us rapid and complete mixing. WASP7 solves the equations according to the conservation of mass principle. The three main classes of water quality processes (loading, transformation and transportation), are represented by these equations. Out of these three processes, about three components play a major role with concentration variability along the river reach, and these processes include: dispersion, advection-dispersion and kinetic transformation (biological, physical and chemical transformation). WASP7 assessed every component of water quality from the point of temporal and spatial input to the point of export and conserving mass in space and time, which is its final point. A finite-difference equation is used for every segment, when temporal and spatial variations in concentration of the constituent are being accounted for. For easy use, derivation of the finite difference form of the equation of mass balance is for a 1D reach. For every segment, the concentration is calculated. The final concentration which is derived from the previous segment is the initial value used for each segment at time zero. In a body of water, the dissolved components representing the three main grades of water quality processes-transport (term 1), loading (term 2) and transformation (term 3), the equation of mass balance is (Equation (46)):

$\frac{\partial \left(AC\right)}{\partial t}=\frac{\partial }{\partial x}\left(-u_{x}AC+E_{x}A\frac{\partial \left(C\right)}{\partial t}\right)+A\left(S_L+S_B\right)+A{S}_k$ (46)

where C is concentration of water quality (g/m³), t is time (day), U_x is longitudinal velocity (m/day), E_x is longitudinal diffusion coefficient (m²/day), S_L is diffusion loading rate (g/m³ day), S_B is boundary loading rate including upstream, downstream, benthic, atmospheric (g/m³ day), S_k = transformation term (total kinetic transformation rate; positive is source, negative is sink, g/m³ day for variable i in a segment) and A is cross-sectional area (m²). The physical and chemical processes, that affect the transport and interaction among the nutrients, phytoplankton, carbonaceous material, sediment, atmosphere and dissolved oxygen in the aquatic environment, is shown in Figure 4.

2.5.2. Processes

Various physico-chemical processes have effect on periphyton, carbonaceous material, phytoplankton, transport and interactions between nutrients and dissolved oxygen in the aquatic environment [33]. The major kinetic interactions for dissolved oxygen and nutrient cycling as modeled by WASP7, is represented by Figure 4. The sub-model of EUTRO consists of phytoplankton, nutrients, Biochemical oxygen demand and dissolved oxygen.

Dissolved oxygen (DO): anaerobic processes in the underlying sediments and aerobic respiratory processes in the water column results in a decrease in DO. Phytoplankton growth, re-aeration result in increase in DO and sediment oxygen demand, phytoplankton respiration and oxidation of CBOD result in loss of DO (Equation (47)):

Figure 4. State variable interactions in advanced eutrophication submodel in WASP7 (Phyto is phytoplankton as carbon, NO₃ is nitrate, NH₄ is ammonium, PO₄ is ortho-phosphorus, CBOD is carbonaceous biochemical oxygen demand, DO is dissolved oxygen ON is organic nitrogen, OP is organic phosphorous, DOM is dissolved organic matter, SS is inorganic suspended solids) [36].

$\begin{array}{c}\frac{\partial c_b}{\partial t}=k_a\left(c_{sat}-c_6\right)-k_d\left[\frac{C_6}{k_{BOD}+C_6}\right]c_1-\frac{SOD}{D}\\ +G_{P_1}\left(\frac{32}{12}+\frac{48}{14}\frac{14}{12}\left(1-P_{N{H}_3}\right)\right)c_4-\frac{32}{12}{k}_{1^R}{c}_4\end{array}$ (47)

where c₆ = DO (mg/L), c₅ = CBOD (mg/L), c_sat = saturated concentration of DO (mg/L), c₅ = biochemical oxygen demand (mgO₂/L), c₁ = ammonia-nitrogen (mgN/L), k_d = de-oxygenation rate (1/day), k_BOD = half saturation constant for oxygen limitation for CBOD (mg/O₂/L), k₁₂ = nitrification rate (1/day), k_NIT = half saturation constant for oxygen limitation for nitrification (mgN/L), k_a = re-aeration rate (1/day), D = depth of water (m), SOD = sediment oxygen demand (mg/m²/day), $P_{N{H}_3}$ = preference for ammonia term, G_p1 = phytoplankton growth rate (1/day), c₄ = the phytoplankton biomass in carbon units (mgC/L) and k_1R = phytoplankton respiration rate (1/day).

The covar method’s flow-induced re-aeration is calculated by EUTRO. This method used one of three formulas (Owens, Churchill, O’Connor-Dobbins) to calculate reaeration, as a function of depth and velocity. Churchill or O’Connor-Dobbins is used for segments that are more than 2 feet. For segments less than two feet, the Owen’s formula is used. Carbonaceous Biochemical Oxygen demand (c₅): detritus phytoplankton carbon derived from the death of algae, is a major source of CBOD apart from man-made sources. Therefore, CBOD is lost through settling, denitrification and oxidation and increases through death of phytoplankton (Equation (48)):

$\begin{array}{c}\frac{\partial c_5}{\partial t}=a_{Oc}{k}_{1d}{c}_4-k_d\left(\frac{C_6}{k_{B0D}+C_6}\right)C_5-\frac{v_{s_3}\left(1-F_{d_5}\right)}{D}{C}_5\\ -\frac{5}{4}\frac{32}{14}{k}_{2d}+\left(\frac{k_{N{O}_3}}{k_{N{O}_3}}+C_6\right)C_2\end{array}$ (48)

where a_oc = oxygen to carbon ratio = 32/12 (mgO₂/mgC), k_1d = phytoplankton death rate (1/day), k_d = CBOD deoxygenation rate (1/day), k_2d = denitrification rate (1/day), $k_{N{O}_3}$ = half saturation constant for oxygen limitation for denitrification (mgN/L), f_d5 = dissolved fraction of CBOD, V_s3 = settling velocity of organic matters (m/day) and c₂ = nitrate nitrogen (mg/L).

The internal CBOD documented by EUTRO is unable to make direct comparisons between model output and BOD₅ data, as measurements from fields may be affected by the decay of algal carbon and algal respiration. Thus, this requires corrective measures taken towards the internally computed CBOD model. This gives rise to a new variable called the bottle BOD₅, and is given by,

$BO{D}_5=c_5\left(1-{\text{e}}^{-k_{dbot}}\right)+\frac{64}{14}{c}_1\left(1-{\text{e}}^{-5k_{nbot}}\right)+a_{0c}{c}_4\left(1-{\text{e}}^{-5k_{1R}}\right)$ (49)

where c₅ = the internally computed CBOD (mg/L), C₁ = the internally computed NH₃, (mg/L), C₄ = the phytoplankton biomass in carbon units (mg/L), a_oc = the oxygen to carbon ratio = 32/12 (mgO₂/mgC), k_dbot = the laboratory “bottle” deoxygenation rate constant (1/day), k_nbot = the laboratory “bottle” nitrification rate constant (1/day) and k_1R = the algal respiration rate at 20˚C (1/day)

Phytoplankton phosphorus (c₄): dissolved inorganic phosphorus is taken up, stored and converted to biomass during the growth of the phytoplankton. As each mg of phytoplankton carbon is produced, the amount of inorganic phosphorus taken up is a_pc mg. Biomass is converted to both non-living inorganic and organic matter, as phytoplanktons undergo respiration and death. There is a release of a_pc mg of phosphorus, for each mg of phytoplankton carbon lost. Phytoplankton phosphorus equation is given by (Equation (50)):

$\frac{\partial \left(C_4{a}_{PC}\right)}{\partial t}=\left(G_{p\Gamma }D{p}_1-D_{P_1}-\frac{v_{54}}{D}\right)a_{Pc}{c}_4$ (50)

where a_pc = phosphorus to carbon ratio (mgP/mgC), G_p1 = phytoplankton growth rate (1/day), D_p1 = phytoplankton death rate (1/day), V_s4 = phytoplankton settling velocity (m/day).

Inorganic and organic phosphorus: there is an interaction between particulate inorganic phosphorus and dissolved inorganic phosphorus (c₃) through a sorption-desorption mechanism. There is a return of phosphorus to particulate organic form and Dissolved inorganic form, from pool of phytoplankton biomass via non-predatory mortality and endogenous respiration. At a temperature dependent mineralization rate, dissolved inorganic phosphorus is derived from organic phosphorus. For phytoplankton yo utilize it, non-living organic phosphorus must ho through bacterial decomposition or mineralization into inorganic phosphorus. EUTRO utilizes a saturating recycle mechanism; which is a compromise between the usual first-order kinetics and a second-order mechanism for recycle, in which the rate of recycle is directly proportional to the phytoplankton biomass present. Saturating recycle allows for second-order dependency, at reduced phytoplankton concentrations, when $P_c\ll K_{mPc}$ (where K_mPc is the half-saturation constant for a cycle) and allows for first-order recycle when there is an increase in phytoplankton, to an amount greater than the half-saturation constant. If the population of the phytoplankton is small, the mechanism reduces the rate of recycle, and doesn’t allow for simultaneous increase in rate with increase in phytoplankton. The hypothesis is that, recycle kinetics occur at the maximum first-order rate with higher levels of population. Zero, is the default value for K_mPc, and this is the point at which it causes mineralization to occur in first-order rate, at all classes of phytoplanktons. Phytoplankton growth results in loss of inorganic phosphorus and mineralization, phytoplankton death result in gain of inorganic phosphorus. The equations for inorganic and organic phosphorus (Equation (52)) are:

$\frac{\partial C_8}{\partial t}=\left(D_{P_1}{a}_{pc}{f}_{0}p\right)C_4-k_{83}\left[\frac{C_4}{k_{mPC}+C_4}\right]C_8$ (51)

$\frac{\partial C_3}{\partial t}=\left(D_{P_1}{a}_{pc}\left(1-f_{0}p\right)\right)C_4-k_{83}\left[\frac{C_4}{k_{mPC}+C_4}\right]C_8-G_{p_1}{a}_{p_c}{c}_4$ (52)

where c₃ = phosphate phosphorus (mg/L), c₈ = organic phosphorus (mg/L), f_op = fraction of dead and respired phytoplankton recycled to the organic phosphorus pool, f_d8 = fraction dissolved organic phosphorus, k₈₃ = dissolved organic phosphorus mineralization (1/day), k_mPc = half saturation constant for phytoplankton limitation of phosphorus recycle (mgC/L) and V_s3 = organic matter settling velocity (m/day).

Phytoplankton nitrogen: as phytoplanktons grow, dissolved inorganic nitrogen is taken up and incorporated into biomass. For every mg of phytoplankton carbon produced, a_nc mg of inorganic nitrogen is taken up. Both ammonia and nitrate are available for uptake and use in cell growth by phytoplankton; however, for physiological reasons, the preferred form is ammonia nitrogen. The phytoplankton nitrogen is gained during growth and lost during death and settling (Equation (53); V_s4 = settling velocity (m/day) and a_nc = nitrogen to carbon ratio = 0.25 mgN/mgC):

$\frac{\partial \left(C_4{a}_{nC}\right)}{\partial t}=\left(G_{p_1}-D_{PT}-\frac{v_{s_4}}{D}\right)c_4{a}_{nc}$ (53)

Organic and inorganic nitrogen: as phytoplankton respires and dies, living organic material is recycled to nonliving organic and inorganic matter. For every mg of phytoplankton carbon consumed or lost, a_nc mg of nitrogen is released. During phytoplankton respiration and death, a fraction of the cellular nitrogen f_on is organic, while (1-f_on) is in the inorganic form of ammonia nitrogen

Nonliving organic nitrogen must undergo mineralization or bacterial decomposition into ammonia nitrogen before utilization by phytoplankton. The first-order rate constant is modified by a saturated recycle mechanism, a compromise between conventional first-order kinetics and a second order recycle mechanism, wherein the recycle rate is directly proportional to the phytoplankton biomass present. This mechanism slows the mineralization rate if the phytoplankton population is small, but does not permit the rate to increase continuously as phytoplankton increase. Ammonia nitrogen, in the presence of nitrifying bacteria and oxygen, is converted to nitrate nitrogen. The nitrification in natural waters is a two-step process with nitrosomonas bacteria responsible for the conversion of ammonia to nitrite and nitrobacter responsible for the conversion of nitrite to nitrate. The nitrification in natural waters depends upon dissolved oxygen, pH, and flow conditions. This process appears to be affected by high or low values of pH that inhibit nitrosomonas growth, particularly for pH below 7 and greater than 9. The nitrifying bacterial populations are sensitive to flow. During periods of high flow, upstream bacteria may be advected downstream with some lag time after a flow transient, before they can build up to significant levels again. Organic nitrogen (c₇) is gained by phytoplankton death and lost by mineralization and settling (Equation (54)). Ammonia nitrogen (c₁) is gained by phytoplankton death and mineralization and lost by nitrification & phytoplankton growth (Equation (55)). The nitrate nitrogen (c₂) is gained via nitrification and lost through phytoplankton growth and de-nitrification (Equation (56)). The reaction equations involving organic nitrogen (c₇), ammonia nitrogen (c₁) and nitrate nitrogen (c₂) are given by:

$\frac{\partial C_7}{\partial t}=D_{P_1}\left(a_{nC}{f}_{on}\right)C_4-K_{71}\left[\frac{C_4}{k_{mpC}+C_4}\right]-\frac{v_{s_3}\left(1-f_{d_7}\right)}{D}{C}_7$ (54)

$\begin{array}{c}\frac{\partial C_1}{\partial t}=D_{P_1}{a}_{nC}\left(1-f_{on}\right)C_4-K_{71}\left[\frac{C_4}{k_{mpC}+C_4}\right]C_7\\ -k_1\partial \left[\frac{c_6}{k_{NIT}+c_6}\right]c_1-G_{p_1}{a}_{nc}{p}_{N{H}_3}{C}_4\end{array}$ (55)

$\frac{\partial C_2}{\partial t}=k_{12}\partial \left[\frac{c_6}{k_{NIT}+c_6}\right]c_1-G_{p_1}{a}_{nc}{p}_{N{H}_3}{C}_4-k_{2}D\left(\frac{k_{N{O}_3}}{k_{N{O}_3}+C_6}\right)C_2$ (56)

$P_{N{H}_3}=\left[\frac{C_2}{\left(k_{m_N}+C_1\right)\left(k_{m_N}+C_2\right)}\right]C_1+\left[\frac{C_2}{\left(C_2+C_1\right)\left(k_{m_N}+C_2\right)}\right]C_2$ (57)

where c₇ = organic nitrogen (mg/L), c₁ = ammonia nitrogen (mg/L), c₂ = nitrate nitrogen (mg/L), f_on = fraction of dead and respired phytoplankton recycled to the organic nitrogen pool, k₇₁ = organic nitrogen mineralization rate (1/day), f_d7 = fraction of dissolved organic nitrogen, k_mN = Michaelis value for ammonia preference (μgN/L), k_mPc = half saturation constant for phytoplankton limitation of phosphorous recycle (mgC/L), $k_{N{O}_3}$ = half saturation constant for oxygen limitation of de-nitrification (mgO₂/L), k₁₂ = nitrification rate (/day) and k_2D = de-nitrification rate (/day). The temperature effect for all first-order reactions used in the model is represented by:

$k\left(T\right)=k\left(20\right){\theta }^{T-20}$ (58)

where k(T) = the reaction rate (/day) at temperature T (˚C) and θ = the temperature coefficient for the reaction.

2.5.3. Input Data

To run WASP7 the user must first insert the data, which includes: output control; model segmentations; boundary concentrations; point and diffuse source waste loads; kinetic parameters; constants, time series flow and initial concentrations.

2.5.4. Model Capability

One of the key abilities of WASP7 involves the analysis of various water quality problems in different water bodies. It works by stimulating the transport and transformation reactions of the state variables identified. With the WASP7, these state variables range from 10 to 14, unlike previous versions only worked with 8 state variables [33] [38]. The state variables for WASP7 are; CBOD (fast, intermediate, slow), salinity, temperature, pesticides, heavy metals, organic chemicals, inorganic solids, conservative tracer, diagenesis, coliform bacteria, cohesive sediments, non-cohesive sediment, silica, and inorganic solids, P(organic, inorganic), phytoplankton and periphyton (bottom algae C,N.P), particulate detritus (N,P,C) [33] [36].

2.5.5. Limitations of Model

There are a couple of limitations of this model in relation to other models out there. First, the model does not take into consideration mixing zones or near field effects. Additionally, it does not employ the use of sinkable/floatable materials. It's calibration and verification process requires a massive amount of data, and it also involves the use of vast site-specific linkage efforts along with multidimensional hydrodynamic models.

2.5.6. Strengths and Its Application

One of the model’s strength is its versatility when it comes to modeling organic chemicals. It can model 25 forms of each chemical. The model can be used for 1D, 2D, and 3D analysis. Over time, it has been used to stimulate a broad array of nutrients, organic compounds, PCBs, and heavy metals, in various water bodies such as the Gulf of Greece [33] [37] [39] [40].

2.6. QUASAR

The primary function of the QUASAR is to analyze the environmental impact of pollutants on river water quality. The QUASAR was an offspring of the Bedford Ouse Study, which was aimed at stimulating the dynamic behavior of river water quality or flow [4] and [41]. The model was first used to along the river to collect telemetered data and present forecasts at specific abstraction sites [42].

In addition, the model was applied within a stochastic or Monte Carlo framework to garner information on the distribution of water quality in rivers exposed to a high level of effluent discharges. To determine the water quality and flow of specific downstream points of non-tidal rivers, the model makes use of point and diffuse elements, as well as upstream inputs from tributaries. To obtain values, it takes into consideration inputs from tributaries downstream, mass balance of flow and water qualities downstream and effluent discharges. It also takes into account chemical decay processes and biological behavior as well as abstractions.

Access to the QUASAR is free, and it is available for PC users on http://www.ceh.ac.uk/. Because of its performance levels, and quality of results, most river regulators use it to set effluent consent levels to meet the various quality standards. It can also be used in estimating the river quality and flow at various points across a defined timeline. With the results, the quality or flow can be analyzed, and also any proposed changes in the river can be carried out.

The versatility of the model is one of its significant strengths, and it can be used in both planning and dynamic mode. The planning mode employs the Monte Carlo simulation method cumulative frequency distribution of selected water quality variables from a given set of hydrological inputs and operating conditions. The data obtained from this mode, ensure standards are upheld in relation to river quality targets. On the other hand, in the dynamic mode, time-series data are used. These data are put in the model to generate estimates of flow and quality at specific boundaries over a given period of time. This mode generates two output data; user-edited values, and observed values.

2.6.1. Conceptualization

The working principles of QUASAR are based on the mass balance on each reach of the river that clearly explained in shown in Figure 5. Each river reach is modeled as a well-mixed tank connected in series. QUASAR shares some similarities and differences with other models like QUAL2EU. Some of the similarities include splitting of the reach into a various number of mixed elements with equal sizes. On the other hand, one significant difference is that in QUASAR, influences are only added to the top element. The duration of solutes in each reach

Figure 5. Interacting water quality state variables in QUASAR model [4] [45].

is calculated using the stream velocity [4].

$\frac{\text{d}c}{\text{d}t}=\frac{1}{T}\left[C_i-C\right]+\Delta C$ (59)

where ΔC is a term representing the net accumulation in the reach of that determinant due to internal transformations, i.e., the “sources” minus the “sinks”, C is the concentration of the solute in the reach, C_i is the concentration entering the reach and τ is the residence time

$\frac{\text{d}\left[{\text{NO}}_{\text{3}}\text{N}\right]}{\text{d}t}=\frac{\left[{\left({\text{NO}}_{\text{3}}\text{N}\right)}^1-\left({\text{NO}}_{\text{3}}\text{N}\right)\right]}{T}+k_1\left[{\text{NH}}_{\text{4}}\text{N}\right]-K_2\left[{\text{NO}}_{\text{3}}\text{N}\right]$ (60)

where [NO₃N]’ and [NO₃N] are the input and output concentration of nitrate nitrogen, [NH₄N] is the concentration of ammonium nitrogen, K₁ is the nitrification rate coefficient and K₂ is the de-nitrification rate coefficient and τ is the residence time (reach length/water velocity). The ammonium nitrate concentration is reduced by nitrification [45] as

$\frac{\text{d}\left[{\text{NH}}_{\text{4}}\text{N}\right]}{\text{d}t}={\left[{\text{NH}}_{\text{4}}\text{N}\right]}^{\prime }-\left[{\text{NH}}_{\text{4}}\text{N}\right]$ (61)

where [NH₄N]’ and [NH₄N] are input and output concentrations of ammonia nitrogen. The concentration of DO in the river is simulated (Equation (62)) as being affected by algal photosynthesis and respiration, a sediment oxygen demand, re-aeration, nitrification and BOD [45] as,

$\frac{\text{d}c}{\text{d}t}=\left(c^{\prime }-c+c_{\omega e_i\gamma }\right)+p-R-k_{3}c+k_4\left(c_5-C\right)-4.57k_1\left[{\text{NH}}_{\text{4}}\text{-N}\right]-K_{5}L$ (62)

where C’ and C are the input and output dissolved oxygen concentrations and CWEIR is the increase in DO concentration due to a weir, P is the rate of photosynthetic oxygen production, R is the rate of oxygen uptake due to algal respiration, K₃ is the rate coefficient for the sediment oxygen demand, K₄ is the re-aeration rate coefficient, C_s is the DO saturation concentration, L is the BOD concentration and K₅ is the rate parameter for BOD oxidation. BOD in the water column is simulated (Equation (63)) as being removed by biological oxidation and sedimentation and there is a contribution to the BOD due to dead algae [4] :

$\frac{\text{d}L}{\text{d}t}=\frac{\left(L^{\prime }-L\right)}{T}-k_5\cdot L-k_{6}L-k_{7}chla$ (63)

where L’ is the input of BOD concentration, K₆ is rate coefficient for sedimentation, K₇ is the rate coefficient for the increase in BOD due to algal death, and Chla is the concentration of chlorophyll-a representing the algal (biomass) concentration. Temperature (T) dependability of rate constants (1/day) k₁, k₂, k₃, k₄ and k₅ are given by

(64a)

(64b)

(64c)

(64d)

(64e)

2.6.2. Input Data

QUASAR operates on four sets of input data. They are; (a) A catchment structure consisting of a river map (b) Boundary conditions that define the water quality and flow of the tributaries and the water at the top of the river (c) Two independently measured water quality data sets at points down the river network to allow model calibration and validation (d) Reach parameters consisting of values of the rate coefficients specific to each reach.

In dynamic mode, the QUASAR inputs are time series of flow and quality. Whereas in the stochastic mode, the user supplies flow and quality data as means and standard deviations of a normal or lognormal distribution, or as lower and upper boundaries of a rectangular distribution.

2.6.3. Model Capability

Some of its modeling capabilities include; temperature, nitrate, dissolved oxygen (DO), pH, algae, E. coli, biochemical oxygen demand (BOD) and conservative pollutant. Its modeling abilities range from dead algae to biochemical oxidation demand.

2.6.4. Limitations of Model

The model is more suitable for dynamic simulations than stochastic simulations. As a stochastic model, QUASAR is very limited, and it does not include correlations and distribution types as seen in other models.

2.6.5. Strengths and Its Application

One of the strengths of the model is its versatility. It is suitable for both dynamic and planning purposes. For the former it can be used to model sizeable freshwater river systems as long as there are sufficient data available [3]. It has been put to use over the years in several areas such as the river Pelenna (in Wales, UK) to assess heavy metal pollution [46], the river Thames (UK) to assess the movement and distribution of nitrates and algae along the Fiver system [47], for modeling impacts of land use and climate change on nitrate-nitrogen in UK [44] and land ocean interaction study [48].

3. Summary and Conclusion

Several parameters differentiate the public domain water quality models used for rivers and streams. These differences can be seen in terms of modeling capability, data input requirements, processes represented, assumptions, as well as strengths and weaknesses. Used appropriately, these models can deliver exceptional results. The six significant public domain water quality models chosen we have looked at in this review are; SIMCAT, TOMCAT, QUAL2EU, QUAL2Kw, WASP7, and QUASAR. From the models reviewed, SIMCAT and TOMCAT are basic. They are used when certain factors like photosynthesis, respiration, and sediment oxygen demand are not essential. They are also used when there is a need to quickly assess impact of point sources. The QUAL2EU has lack of provision for conversion of algal death to CBOD. It falls short in areas where macrophytes (rooted aquatic plants) play a vital part. The QUAL2Kw has an automatic calibration system and new constituent interactions, such as denitrification, algal BOD, DO change caused by fixed plants, hyporheic exchange, and sediment pore-water quality. The current version of the model is suited to simulate branches of the river systems. Finally, both WASP7 and QUASAR have their shortcomings. They require an extensive array of data which requires lots of time and considerable cost. However, used rightly, they can deliver exceptional results. The review of available model indicates not any one model can provide adequate service needed for each evaluation, and as such selecting a model to use will depend on time, cost, and area of application.

Appendix

Table A1. Comparison of models.

CSRTS continiously stirred tank reactors in the series, CBOD carbonaceous biochemical oxygen demand, SOD seediment oxygen demand, OP organic phosphorous, ON organic nitrogen, DO dissolved oxygen, temp water temperature, ADE advection disperion equation, TIC total inorganic carbon, organic chemicals, NCPAR non-conservative parameter.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Agoshkov, V.I. (2002) Mathematical Models of Life Support Systems. In: Agoshkov, V.I., Ed., Knowledge for Sustainable Development: An Insight into the Encyclopaedia of Life Support Systems, Vol. 1, UNESCO/EOLSS, Paris, 335-281.
[2]	USEPA (2009) Water Quality Models and Tools. United States Environmental Protection Agency (USEPA), Washington DC. http://www.epa.gov/waterscience/models
[3]	Chapra, S.C. (1997) Surface Water-Quality Modeling. McGraw-Hill International Edition, New York, 844.
[4]	Cox, B.A. (2003) A Review of Currently Available In-Stream Water-Quality Models and Their Applicability for Simulating Dissolved Oxygen in Lowland Rivers. The Science of the Total Environment, 314, 335-377. https://doi.org/10.1016/S0048-9697(03)00063-9
[5]	Reckhow, K.H. (1994) Water-Quality Simulation Modeling and Uncertainty Analysis for Risk Assessment and Decision-Making. Ecological Modelling, 72, 1-20. https://doi.org/10.1016/0304-3800(94)90143-0
[6]	Warn, A.E. (1987) SIMCAT—A Catchment Simulation Model for Planning Investment for River Quality. IAWPRC, Pergamon, Oxford, 211-218.
[7]	Crabtree, B., Seward, A.J. and Thompson, L. (2006) A Case Study of Regional Catchment Water Quality Modelling to Identify Pollution Control Requirements. Water Science and Technology, 53, 47-54. https://doi.org/10.2166/wst.2006.296
[8]	Elmore, H.L. and Hayes, T.W. (1960) Solubility of Atmospheric Oxygen in Water, Twenty-Ninth Progress Report of the Committee on Sanitary Engineering Research. Journal of the Sanitary Engineering Division, 86, 41-53.
[9]	Crabtree, B., Hickman, M. and Martin, D. (1999) Integrated Water Quality and Environmental Cost-Benefit Modelling for the Management of the River Tame. Water Science and Technology, 39, 213-220. https://doi.org/10.2166/wst.1999.0208
[10]	Bowden, K. and Brown, S.R. (1984) Relating Effluent Control Parameters to River Quality Objectives Using a Generalized Catchment Simulation Model. Water Science and Technology, 16, 197-206. https://doi.org/10.2166/wst.1984.0132
[11]	Keller, V. (2005) Risk Assessment of “Down-the-Drain” Chemicals: Search for a Suitable Model. The Science of the Total Environment, 360, 305-318. https://doi.org/10.1016/j.scitotenv.2005.08.042
[12]	Kinniburgh, J.H., Tinsley, M.R. and Bennett, J. (1997) Orthophosphate Concentrations in the River Thames. Journal of the Chartered Institution of Water and Environmental Management, 11, 178-185. https://doi.org/10.1111/j.1747-6593.1997.tb00113.x
[13]	Brown, L.C. and Barnwell, T.O. (1987) The Enhanced Stream Water Quality Models QUAL2E and QUAL2E-UNCAS: Documentation and User Manual. Env. Res. Lab. USEPA, EPA/600/3-87/007, Athens, 204.
[14]	Rauch, W., Henze, M., Koncsos, L., Reichert, P., Shanahan, P., Somlyody, L., et al. (1998) River Water Quality Modelling: I. State of the Art. Water Science and Technology, 38, 237-244. https://doi.org/10.2166/wst.1998.0473
[15]	Streeter, W.H. and Phelps, E.B. (1925) A Study of the Pollution and Natural Purification of the Ohio River. Public Health Service, USA Public Health Bulletin 146, Washington DC.
[16]	Orlob, G.T.E. (1982) Mathematical Modelling of Water Quality: Streams, Lakes and Reservoirs. Wiley, Chichester, 542.
[17]	Walton, R. and Webb, M. (1994) Qual2e Simulations of Pulse Loads. Journal of Environmental Engineering, 120, 1017-1031. https://doi.org/10.1061/(ASCE)0733-9372(1994)120:5(1017)
[18]	USEPA (1995) QUAL2E Windows Interface User’s Guide. United States Environmental Protection Agency, EPA/823/B/95/003, Athens, 66.
[19]	Ambrose, R.B., Wool, T.A., Connolly, J.P. and Shanz, R.W. (1987) WASP5, a Hydrodynamic and Water Quality Model. U.S. Environmental Protection Agency, EPA/600/3-87/039, Athens.
[20]	Park, S.S. and Uchrin, C.G. (1997) A Stoichiometric Model for Water Quality Interactions in Macrophyte Dominated Water Bodies. Ecological Modelling, 96, 165-174. https://doi.org/10.1016/S0304-3800(96)00064-6
[21]	Park, S.S. and Lee, Y.S. (2002) A Water Quality Modeling Study of the Nakdong River, Korea. Ecological Modelling, 152, 65-75. https://doi.org/10.1016/S0304-3800(01)00489-6
[22]	Barnwell, T.O., Brown, L.C. and Whittemore, R.C. (2004) Importance of Field Data in Stream Water Quality Modeling Using QUAL2E-UNCAS. Journal of Environmental Engineering, 130, 643-647. https://doi.org/10.1061/(ASCE)0733-9372(2004)130:6(643)
[23]	Cubillo, F., Rodriguez, B. and Barnwell, T.O. (1992) A System for Control of River Water-Quality for the Community of Madrid Using Qual2e. Water Science and Technology, 26, 1867-1873. https://doi.org/10.2166/wst.1992.0631
[24]	Ning, S.K., Chang, N.B., Yang, L., Chen, H.W. and Hsu, H.Y. (2001) Assessing Pollution Prevention Program by QUAL2E Simulation Analysis for the Kao-Ping River Basin, Taiwan. Journal of Environmental Management, 61, 61-76. https://doi.org/10.1006/jema.2000.0397
[25]	Chaudhury, R.R., Sobrinho, J.A.H., Wright, R.M. and Sreenivas, M. (1998) Dissolved Oxygen Modeling of the Blackstone River (Northeastern United States). Water Research, 32, 2400-2412. https://doi.org/10.1016/S0043-1354(98)00004-9
[26]	Pelletier, G.J. and Chapra, S.C. (2005) QUAL2Kw Theory and Documentation (Version 5.1), a Modeling Framework for Simulating River and Stream Water Quality. Department of Ecology, Washington DC.
[27]	Chapra, S.C. and Pelletier, G.J. (2003) QUAL2K: A Modeling Framework for Simulating River and Stream Water Quality (Beta Version): Documentation and Users Manual. Civil and Environmental Engineering Department. Tufts University, Medford.
[28]	Pelletier, G.J., Chapra, S.C. and Tao, H. (2006) QUAL2Kw—A Framework for Modeling Water Quality in Streams and Rivers Using a Genetic Algorithm for Calibration. Environmental Modelling and Software, 21, 419-425. https://doi.org/10.1016/j.envsoft.2005.07.002
[29]	Turner, D.F., Pelletier, G.J. and Kasper, B. (2009) Dissolved Oxygen and pH Modeling of a Periphyton Dominated, Nutrient Enriched River. Journal of Environmental Engineering, 135, 645-652. https://doi.org/10.1061/(ASCE)0733-9372(2009)135:8(645)
[30]	Kannel, P.R., Lee, S., Lee, Y.S., Kanel, S.R. and Pelletier, G.J. (2007) Application of Automated QUAL2Kw for Water Quality Modeling and Management in the Bagmati River, Nepal. Ecological Modelling, 202, 503-517. https://doi.org/10.1016/j.ecolmodel.2006.12.033
[31]	Cristea, N. and Pelletier, G. (2005) Wenatchee River Temperature Total Maximum Daily Load Study. Washington State Department of Ecology, Olympia.
[32]	Ambrose, R.B. and Wool, T.A. (2009) WASP7 Stream Transport Model Theory and User’s Guide. U.S. Environmental Protection Agency, EPA/600/R-09/100, Athens.
[33]	Wool, T.A., Ambrose, R.B., Martin, J.L. and Comer, E.A. (2001) Water Quality Analysis Simulation Program (WASP) Version 6.0: User’s Manual. U.S. Environmental Protection Agency, Athens.
[34]	Ditoro, D.M., Fitzpatrick, J.J. and Thomann, R.V. (1983) Documentation for Water Analysis Simulation Program (WASP) and Model Verification Program (MVP) Westwood: Hydroscience. USEPA Contract No. 68-01-3872.
[35]	Ambrose, R.B., Wool, T.A. and Connolly, J.P. (1988) WASP, Version 4, a Hydrodynamic and Water Quality Model—Model Theory, User’s Manual, and Programmer’s Guide. U.S. Environmental Protection Agency, Athens.
[36]	Ambrose, R.B., Martin, J.L. and Wool, T.A. (2006) WASP7 Benthic Algae—Model Theory and User’s Guide. U.S. Environmental Protection Agency, EPA/600/R-06/106 (NTIS PB2007-100139), Washington DC.
[37]	Nikolaidis, N.P., Karageorgis, A.P., Kapsimalis, V., Marconis, G., Drakopoulou, P., Kontoyiannis, H., et al. (2006) Circulation and Nutrient Modeling of Thermaikos Gulf, Greece. Journal of Marine Systems, 60, 51-62. https://doi.org/10.1016/j.jmarsys.2005.11.007
[38]	Ambrose, R.B., Wool, T.A. and Martin, J.L. (1993) The Water Quality Simulation Program (WASP), Version 5, Part A, Model Documentation. Environmental Research Lab, Athens.
[39]	Rygwelski, K.R., Richardson, W.L. and Endicott, D.D. (1999) A Screening-Level Model Evaluation of Atrazine in the Lake Michigan Basin. Journal of Great Lakes Research, 25, 94-106. https://doi.org/10.1016/S0380-1330(99)70719-7
[40]	Stansbury, J. and Admiraal, D.M. (2004) Modeling to Evaluate Macrophyte Induced Impacts to Dissolved Oxygen in a Tailwater Reservoir. Journal of the American Water Resources Association, 40, 1483-1497. https://doi.org/10.1111/j.1752-1688.2004.tb01600.x
[41]	Whitehead, P., Beck, B. and Oconnell, E. (1981) A Systems-Model of Streamflow and Water-Quality in the Bedford Ouse River System. 2. Water-Quality Modeling. Water Research, 15, 1157-1171. https://doi.org/10.1016/0043-1354(81)90091-9
[42]	Whitehead, P.G., Caddy, D.E. and Templeman, R.F. (1984) An Online Monitoring, Data Management and Forecasting System for the Bedford Ouse River Basin. Water Science and Technology, 16, 295-314. https://doi.org/10.2166/wst.1984.0139
[43]	Whitehead, P. and Young, P. (1979) Water-Quality in River Systems—Monte-Carlo Analysis. Water Resources Research, 15, 451-459. https://doi.org/10.1029/WR015i002p00451
[44]	Ferrier, R.C., Whitehead, P.G., Sefton, C., Edwards, A.C. and Pugh, K. (1995) Modeling Impacts of Land-Use Change and Climate-Change on Nitrate-Nitrogen in the River Don, North-East Scotland. Water Research, 29, 1950-1956. https://doi.org/10.1016/0043-1354(95)00004-5
[45]	Whitehead, P.G., Williams, R.J. and Lewis, D.R. (1997) Quality Simulation along River Systems (QUASAR): Model Theory and Development. The Science of the Total Environment, 194, 447-456. https://doi.org/10.1016/S0048-9697(96)05382-X
[46]	Whitehead, P.G., Mccartney, M.P., Williams, R.J., Ishemo, C.A.L. and Thomas, R. (1995) A Method to Simulate the Impact of Acid-Mine Drainage on River Systems. Journal of the Chartered Institution of Water and Environmental Management, 9, 119-131. https://doi.org/10.1111/j.1747-6593.1995.tb01601.x
[47]	Whitehead, P.G. and Williams, R.J. (1982) A Dynamic Nitrate Balance Model for Fiver Basins, IAHS Exeter Conference Proceedings IAHS Publication.
[48]	Lewis, D.R., Williams, R.J. and Whitehead, P.G. (1997) Quality Simulation along Rivers (QUASAR): An Application to the Yorkshire Ouse. The Science of the Total Environment, 194, 399-418. https://doi.org/10.1016/S0048-9697(96)05379-X