#### Abstract

Migration of pollutant particles into subsurface water reservoirs through point sources is largely involved mixing processes within the system of water flow. Possible potential sources of pollution to these point sources include municipal wastes, septic loads, landfills, uncontrolled hazardous wastes, and sewage storage tanks. The mixing processes of pollutant significantly alter their predictive rate of flow in the water reservoirs, and therefore the time inherent in mixing processes need to be accounted for. In this study, pollution of subsurface water reservoirs mainly rivers and streams through contaminated water point sources (CWPS) was studied through a conceptual perspective of mixing problem processes in water tanks. The objective was to formulate a discrete time delay mathematical model which describes the dynamics of water reservoir pollution that involve single species contaminants such as nitrates, phosphorous, and detergents injecting from a point source. The concentration of pollutants was expressed as a function of the inflow and outflow rates using the principle for the conservation of mass. The major assumption made in modeling of mixing problems using tanks is that mixing is instantaneous. Practical realities dictate that mixing cannot occur instantaneously throughout the tank. So as to accommodate these realities, the study refined the systems of ordinary differential equations (ODEs) generated from principles of mixing problems in cascading tanks, into a system of delayed differential equations (DDEs) so that the concentration of pollutant leaving the reservoir at time would be equal to the average concentration at some earlier instant, for the delay The formulated model is a mathematical discrete time delay model which can be used to describe the dynamics of subsurface water reservoir pollution through a point source. The model was simulated on municipal River Nyakomisaro in Kisii County, Kenya. Physical and kinematic parameters of the river (cross-sectional lengths, depths, flow velocities) at three river sectional reservoirs were measured and the obtained parameter values were then used to evaluate coefficients of the formulated model equation. The system of DDEs from this simulation was solved numerically on MATLAB using dde23 software. From the graphical views generated for concentration of pollutant versus time it was established that the developed DDEs cover longer time series solutions (characteristic curves) than that from the corresponding ODEs in the same reservoir indicating that time necessary for particle flow through water reservoirs is underestimated if ODEs are used to describe particle flow. Also, the graphical views indicated similar tendencies (characteristics) in particle flow with time elapse even though initial values of concentration were different for every potentially recognized single species pollutant considered in each river reservoir. Hence, longer values of time will imply more pollution in the water reservoir and vice versa. By introducing time delays due to constituent mixing processes in water quality simulation models that make use of advection-diffusion equation such as Qual2kw, the findings of this study can help for better understanding of the contaminant’s accumulation levels and their rate of transport in water resource. These will assist, for example, water-quality protection agencies such as Environmental Protection Agency (EPA), World Health Organization (WHO), and National Environmental Management Authority (NEMA) for the need to generate efficient and effective remedial strategies to control or mitigate hazardous or risks arising from water pollution.

#### 1. Introduction

In addition to groundwater, subsurface water systems that include rivers, ponds, dams, and lakes are important water resource for municipalities, domestic, agriculture, and industry. Therefore, the capability to predict the levels of concentration and rate of movement of contaminants in the water resource is of vital importance for the reliable assessment of hazardous or risks arising from water contamination problems, and for the decision of efficient and effective techniques to mitigate them [1].

Reduction of subsurface water pollution largely involves the scientific understanding of these processes which control the fate and movement of pollutants in the subsurface environment. First, the advection transport involves movement of dissolved chemicals with the flowing water. Second, hydrodynamic dispersion, involve molecular and ionic diffusion together with small-scale variations in the flow velocity through the porous media causing the paths of dissolved molecules and ions to diverge or spread from the average direction of water flow. Third, reactions in which some amount of a particular dissolved chemical species (such as herbicides) may be removed from the water as a result of biological, chemical, and physical reactions in water or between the water and the solid aquifer materials or other liquid phases. Lastly, fluid point sources, where water of one constituent level is introduced into and mixed with water of a different constituent level. In addition, it is important to understand the most challenging problems associated with pollution of water systems which include prevention of the introduction of contaminants in water systems, prediction of their movement if they are introduced, and methods of their removal, to some extent, in order to protect the biosphere effectively [2].

The mathematical representation of the transport models to predict the contaminant migration in subsurface water systems is one major area to supplement to the scientific knowledge on pollution processes. Formulation of the transport models involves determination of different parameters in the field and the laboratory using different methods and consequent applications of the established conceptual models for development of management measures to control and/or prevent the introduction of contaminants in water systems [3].

Previous studies on mixing problems that exist in literature are experimentally covered in Ukwu [4], using variation of concentration of brine with time through cascading tanks. The common assumption in these cited literatures while formulating their mixing-problem models is that mixing in the tanks is instantaneous. Practical realities dictate that mixing processes cannot occur throughout the tank at once. Thus, the concentration of salt flowing out of the tank at time, will be equal to the average concentration at some earlier time , where. Lead this way, the huddle with the assumption of instantaneous mixing is settled. Mathematically, however, these criteria involve refinement of ODE models to DDE models. The novelty of this study is the extension of this experimental approach to mathematically model the rate and concentration of particle flow (from point source) of single species pollutants in subsurface water-flow systems.

The sketch in Figure 1 below introduces the general problem involving pollution of subsurface water systems that was covered under this study. In this flow diagram, water and pollutant meet at a point source and the mixture flows downstream in discretized cascades as shown in Figure 2. The pollutant undergoes processes of sedimentation, filtration, and leakages, and finally, the pollutant either accumulates on the water surface or in the soil. The concern of this research is to model the process leading to surface accumulation of pollutant within the water flow beginning from the point source.

This paper is organized as follows. First, the background to the study is described. Then, the governing mathematical principles and equations are given in section 2, model formulation and applications are presented in section 3, summary and conclusion in section 4, and recommendations in section 5.

#### 2. Background

Migration of dissolved (or undissolved) organic and inorganic substances into subsurface water reservoirs through point sources, largely involve mixing processes within the system of flow. The mixing processes significantly alter the predictive rate and concentration of pollutants in the water-flow systems, and therefore, the time inherent in these processes need to be accounted for in modeling rate of pollutant flow. This study considered changes in contaminant concentration which occur within dynamic subsurface water systems of cascading water reservoirs due to fluid point sources, where water of one composition of pollutant is introduced into and mixed with water of a different composition. A numerical model for the process was formulated using mixing-problem analyzes in water tanks [6]. The objective was to apply the formulated model to predict pollution concentration attenuation in water systems and this is helpful to approximate time inherent in the process of pollutants accumulation in water reservoirs.

The following governing equations and principles were found resourceful to this study.

##### 2.1. Continuity Equation

In water dynamical systems, development of mathematical equations to describe the groundwater flow and transport processes is mainly achieved through the fundamental principle of mass conservation of fluid or solute. The general equation for conservation of mass may be represented in global variables, using a representative volume of porous medium as where is the mass of fluid or solute. This is the equation for conservation of mass (or continuity equation) which may be combined with a mathematical expression of the relevant dynamic process to obtain a differential equation describing flow or transport [7].

##### 2.2. River Water Quality Simulation Model

Qual2kw model is one commonly used water quality simulation models. The Qual2kw model is a one-dimensional model for simulation of river and stream water quality variations. This model uses the advection-diffusion equation to simulate the transport and fate of several water constituents. The model deploys the finite-difference numerical scheme for the solution of the adjective-dispersion mass transport and reaction ordinary differential equations.

For a water column, the variation of constituent concentration based on a general mass balance equation (Equation (1)) is given as (Maridi [8]) where is the constituent concentration; is flux at reach ; is abstraction flux at reach ; is bulk coefficient between adjacent reaches; are consecutive cascading reaches; are external loading of constituents; and are sources and sinks of constituents due to reactions and mass transfer mechanisms. The detail description of interacting water quality state variables processes that are effective on equation (2) is described in Pelletier and Chapra [9].

Based on equation (2), this study proposes a criterion for more accurate mathematical evaluation of constituent concentration, by considering the time delays inherent in the mixing processes of the constituents (pollutant loads) with flowing water. These mixing processes greatly contribute to rate of variations of constituent concentration, in the flowing water. On principle, the ordinary differential equation (2) in the given variables and constants is transformed to a corresponding delay differential equation with the concentration where is the time delayed in the mixing process at a point as the flow of pollutant load proceeds downstream. Consequently, is the concentration variation rate in time including the temporal variation of concentration based on mixing, diffusion, advection, reaction, and rate of changes due to pollution load discharged into the river.

This study propagates and applies a wide range of known mathematical theories and properties of delay differential equations (DDEs) to demonstrate the essentiality of the DDEs in describing dynamical systems over their counterpart, the ordinary differential equations (ODEs). In application, the outcome of this analysis is ventured towards the use of DDEs in place of ODEs in water quality simulation models that make use of numerical schemes such as Qual2kw and Ce Qualw2. This mathematical contribution is geared towards advancement of river water assessment and management.

##### 2.3. Equations for Mixing-Problem Processes in Cascade Reservoirs

In modeling population dynamics, Antonin [6] derives the following system of ODE from cascade arrangement of water tanks: where is concentration of pollutant in the cascade; is number of reservoirs in the cascade; and . Equation (3) indicates that fresh water enters the initial water reservoir so that the particle concentration due to a single particle source is negative. Equation (4) is a global ordinary differential equation which describes the change in concentration between two consecutive reservoirsand in the cascade through the 𝑛 reservoirs in the flow. It should be noted that in its application, such as on the dynamics of river flow involving pollution, the system of ordinary differential equation (equation (4)) provides net concentration of particles only between two progressive adjacent cascade reservoirs, the and reservoirs.

##### 2.4. Equation for Single Species Population Model with Delay

In modeling population dynamics, Forde [10] provides DDE model for a single species population as

If , then equation (5)will have nontrivial steady states. The study considered the case of equation (5) when the average death rate is constant, and the average birth rate has an exponential decaying factor. Examples where this type of model has been applied include “Nicholson’s blowfly data” on laboratory populations of blowflies, in which the essential features of the model are two functions relating fecundity and egg-to-adult survival to the amount of food available to each individual. The model describes thresholds representing maintenance requirements in the two functions, and , together with a time delay due to the development that cause cycling in the blowfly populations. Similar thresholds and time delays occur in more complex biological systems with cycling in closely interacting populations of predators and prey [11].

#### 3. Results and Discussions

##### 3.1. Formulated Discrete Time-Delay Mathematical Model

Steady, uniform and three-directional flow in a river, stream or generally subsurface water flow system with cascade geological setting was considered. First, equations (3) and (4) derived from mixing-problem processes in controlled set-up of water tanks [6] were used to provide a conceptual framework and methodology for natural subsurface water reservoirs of cascading geological setting. Second, since pollution concentration in the reservoirs is never instantaneous, the equations were modified to include lag or delay time . This modification considered concentration of pollutant species as a reducing or decaying factor so that the final model equation for particle transport would conform to equation (5) from population dynamics. To achieve this, the changes in concentration of particles in a particular reservoir were analogized to the dynamics of population change of a single species in an ecological niche, and therefore the following derivation was resourceful:

Considering a simple population growth model with unlimited resources, the changes in population size at any given time can be represented in global variables as

The constants and are the birth and death rates, respectively. This model assumes that birth and death rates are directly commensurate! Natural realities in biological populations indicate that these rates are never commensurate. Can this model, therefore, be improved to simulate such realities in a population? Let be the number of adults. The babies of the population species should take time to grow up and with probability to make it to full adults. Then,

This is a delay differential equation (DDE), named Hayes Equation [10]. In particular, by considering the independent variables and equation (7) is similar to equation (4) in applications.

##### 3.2. Model Formulation

The model DDE equation was provided in the form of equation (5) from population dynamics through equation (4) and then parameter setting mathematical approach was used to investigate the stable states of the resulting model equation. This approach involved combinations of a range of chosen values of the time lags and coefficients in the DDE equation so as to obtain a steady stable state.

Following the derivation of Hayes equation in equation (7), equation (5) was modified to the following model equation: where is the position of reservoir in the cascades, is the concentration of pollutant species in subsurface water reservoir, is the time lag of pollutant flow in the cascade reservoir and is considered to be of unique value for a particular reservoir so that for reservoirs . This specifies equation (8) as DDE with multiple discrete delays. The parameters are time-dependent coefficients so that equation (8) is nonautonomous.

The coefficient is termed as the incursion term and represents pollutant particles that enter the reservoir with a rate of flow , which is the rate of flow in the previous adjacent reservoir. is considered a continuous, positive decreasing function, so that the average concentration rate of pollutant decreases with increasing number cascade water reservoirs. Similar consideration on this coefficient has been done for density-limited growth and the logistic models in population dynamics [10].

The coefficient is called the dilution term which accounts for the spreading rate (flux divergence rate) of pollutant particles and therefore reducing the average concentration per unit volume of particles flowing in the ith reservoir. This term is comparable to the death rate of a single species in population dynamics, which is increased by factors like intraspecific competitions for survival. For example, in Nicholson’s blowfly model [10], is considered constant and represents constant per capita death rate. Similarly, in applications of equation (8), 𝜕i in the ith reservoir is considered constant in that particular reservoir.

The invading particles in the reservoir do not mix up instantaneously with water, but rather takes time in the mixing process. Some of these particles will be harbored in the reservoir due to chemical, biological, and physical processes and therefore may or may not eventually flow into the reservoir. With reference to equation (7), Hayes’ equation, this possibility is represented in the model equation (8) by the mathematical probability This term is attached to the incursion term to mathematically account for the progressive exponentially reducing (decaying) incursion rate of pollutant particles in the cascades. The constant is mathematically used as a scaling factor and is assigned a negative value so that the exponential term can reciprocally provide the decaying factor in the incursion term.

The following mathematical definition was provided for particle concentration so as to specify equation (8) for a solution which is bounded for each water reservoir in the cascade.

Definition 1. If , where represent a trajectory of the solution in the past, then represent a functional operator from to . Further, for , define as functional DDE with discrete delays [12].
The consequence of this definition is to specify the model equation (4.3) as a functional delay differential equation with the following properties:
The delays 𝜏𝑖 are discrete with 𝜏1 < 𝜏2 < ⋯ < 𝜏𝑛 for every
Concentration. Initial conditions:
The characteristic equation is transcendental.

###### 3.2.1. Linearized Stability of Model Equation (8)

This subsection was aimed at setting conditions for the determination of values of the coefficients and as a pair in Equation (8) that will always give rise to a stable state on application of the DDE. This was achieved through a criterion derived to evaluate a critical delay, from a suitable combination of values of and . The delays, 𝜏𝑖 will then be determined based on the value of 𝜏𝑐 so as to sustain stability of resulting DDE.

The study derived a lot of insight from stability analysis of equilibria. An equilibrium point in the state space was considered for which is a solution of equation (8) for all , so that the equilibrium points satisfy

In principle, when a solution or solutions of the characteristic equation has positive real part, then the equilibrium point is unstable. If they all have negative real parts, the equilibrium point is stable. When the leading characteristic values are zero, then the stability is irreducible to linear order.

It should be realized that the basic properties of this model are the type of functions and which might lead to the existence of periodic solutions of equation (8). The study specified the use of the case and constant in each of the cascade reservoirs so as to provide periodic solutions of the model equation.

###### 3.2.2. Preliminary Analysis of Model Equation

An analytical approach was provided to establish basic properties of solutions to the equation (8) under the following projections: (i)Given positive initial data of concentration , solutions remain positive in all time . This positivity of solutions is important in modeling physical and biological phenomena because negative solutions (values) will be meaningless(ii)Solutions are bounded and eventually uniform bounded regardless of initial data. Pollution in a flow is characterized with particles moving along a concentration gradient. The rate of the mixing process of pollutant particles is affected by initial concentration. Such rates are expressed in terms of time delays in any dynamics of flow. Eventually a uniform particle flow is attained for every initial concentration when particles concentration is minimum constant or otherwise zero in the long run, or maximum levels(iii)Determine steady state solutions and their stability:(a)Identify delays that cause instability to nontrivial steady states of equation (8). If the nontrivial steady states exist for a solution, then there is a critical delay corresponding to every such nontrivial steady state.(b)Propagate the existence of periodic solutions. DDEs that model physical dynamics exhibit oscillatory phenomena whose analytical results provide periodic solutions [10].

To achieve the above set projections the following analytical approach provided in Forde, [10] was applied. The study relied on the condition that so as to be led by the case of Theorem (2) below, where the nontrivial steady state exists. Besides, for particles to flow in a water reservoir, their rate of “production” must be higher than their “death” rates.

Theorem 2. Let and be positive functions. Suppose that there exists an such that
and Then, is a positive steady state, and the trivial steady state is unstable. If then is linearly stable for all Otherwise, there exists a such that is stable for and unstable for [10].
Using Theorem (2) above, the particulars of the case were analyzed, which is a necessary condition for particle flow to occur and for completeness of the mixing process within a water reservoir. The analysis leads to a condition for the boundedness of the delays in SWR, and the proof of the theorem follows from the application.
From Theorem (2) and considering the case the nontrivial steady state occurs whenever, or Then, according to Theorem (2), is stable for all if and only if This is equivalent to the condition that , which is a contradiction to the necessary requirement that for particles to flow in solution. This condition was used to validate model equation (8) by determining a pair of measured values from the dynamics of water flow involving change of concentration of pollutant particles, so that for unstable solution that is essential for the periodic property of the DDEs. So as to determine the range of values of delay that would give rise to a stable or unstable state solution of equation (8) in its application, the following criterion based on Theorem (2) was applied [10].
The condition 𝛽 > 𝜕𝑒2 was considered for a steady state of equation (8). Laplace transforms show that all analytic solutions of autonomous constant delay DDEs are exponentials [12]. To obtain a characteristic equation, therefore, the simple exponential may be chosen and posted in equation (8) to get, Dividing by and rearranging gives a transcendental exponential, Taking then The roots have two cases to address on the application of equation (8). When there is no delay, such as in the case of flow without mixing processes, resulting to simple polynomial, , which is insignificant as far as pollution is concerned. But, if , then , which is a complex root. Taking the real and imaginary parts, Squaring and summing up yields or, Rewriting the real and imaginary parts of the roots of the characteristic equations, Thus, there exists a critical delay, at which an eigen-value crosses into the right half-plane, and the delay is given by

###### 3.2.3. Time Necessary for Transport of Pollutants through Cascaded SWR

The formulated model equation (8) was solved in Matlab dde23 software to obtain time necessary for transport of single species pollutant in cascading sub-surface water reservoirs (SWR). It should be noted that the basic properties (such as linearity, stability, and oscillatory) of equation (8) are determined by the nature of functions 𝛽 and 𝜕.

The time series solutions for the system of equation (8) were generated graphically through Matlab by assigning paired values of the coefficients and in a particular reservoir that were predetermined for the boundedness of the corresponding delay using equation (21). These time series curves displayed similar tendencies downstream from the point source to the last reservoir in the cascade, and therefore they constitute characteristic curves. Curves resulting from the dynamics in reservoirs nearer the point source have longer times to cover than those far from the point source. The intensity (magnitude) of the point source pollutant also affects the duration of the curves. Higher levels of the pollutant concentration at point source cause longer amplitudes, and vice versa. One remarkable result of the DDE model was obtained on comparing its results to the results of ODEs it modifies under the same dynamic system of flow. The time series solutions obtained for the DDE model displayed longer curves than the ODE counterparts in the same reservoir for the same single species pollutant.

Therefore, the time necessary for transport of pollutant in cascading SWR can be evaluated numerically using the DDE model. Further, since the process of water reservoir pollution is not instantaneous, but rather takes time to occur, the model provides the time inherent for this process as .

###### 3.2.4. Assumptions for the Model

The following assumptions were found necessary for the formulation and application of the model: (i)Flow in subsurface water reservoirs was assumed steady, uniform, and unidirectional(ii)Water free from pollutant-species flows into the first reservoir in the cascade and pollutant in solution flows out of every reservoir at the same rate, ensuring that the amount of fluid flow remains constant in the cascade. This way, the time 𝑡 was determined that accounted for movement of pollutant through the cascades since its injection into water reservoirs through a point source until maximum concentration is reached in the whole system of water flow(iii)Only a single point source of pollutant species into water system was considered(iv)It was assumed that mixing of pollutant particles in water (in every reservoir) was perfect and complete and that the flow rates (in and out) of reservoir were well observed (measurable)

##### 3.3. Model Implementation

The described mathematical model was implemented on the dynamics of flow involving pollution on river Nyakomisaro in Kisii town (Kenya). The river is of shallow water with average depth metres (at minimum during dry season) and 1–2 metres (at maximum during high rainfall) and average width below meters. The river transverses highly populated Kisii town and therefore it is characterized with several point sources of pollution from the municipal solid and dissolved wastes flowing into the river through its banks. Concentration, for pollutant species were selected at particular point sources along the river banks.

From the studies conducted on the river involving pollution [13], the initial pollutant concentrations at the identified point sources for SECTIONS A, B, and C were taken as ppm, ppm, and ppm, respectively.

###### 3.3.1. Determination of Flow Volume

Each reservoir volume was calculated from the length and average cross-sectional area for every reservoir as where Measurements of depth, width, and length were taken downstream at every discretized endpoint of a reservoir using meter-rule and string.

###### 3.3.2. Discretization of Time-Lags

In the physical applications of the DDE equation (8), the lags are parameters distinctively calculated as a ratio of the cascade length to total length of the cascade zone where These proportions were necessary to maintain that so that the time lags are distinctively discrete and in series. The discrete delays were calculated and bound as shown in Table 1 below.

The coefficient (𝑡) for every cascade reservoir was determined using the method of “serial dilution”. The dilution factor (DF) for each step was computed as the ratio of the total volume, (obtained by summing up for where is the number of zones in the cascade) to , volume of water (assumed to be free of single species pollutant) entering the cascade zone. Thus, where is the volume of precascade reservoir which is determined just before entering the first cascade reservoirs. Using this approach, and with the assumption that “freshwater” enters the contaminated cascade through the first reservoir (at the point source), the dilution term was computed as the DF of the system. Thus, the DF up to sixth reservoir was reached as:

###### 3.3.3. Coefficients for Model Equation

Using the results in Table 2, the incursion term for every cascade was calculated from the relation

The computed results for the six cascade zones were recorded in Table 3. This computation was necessary to show that the concentration of the pollutant reduces downstream leading up to sixth reservoir where it is 4096 times less than the original “undiluted” solution inside the first reservoir.

###### 3.3.4. Particle Flow for

The model equation (8) was transformed into its ordinary differential equation counterpart by setting the delays 𝜏𝑖 =0, which resulted to where represents the rate per unit volume (flow flux) at which pollutant particles are entering the reservoir at the boundary point. Because no history of particle flow is required to solve ODEs, a solution is for equation (8) guaranteed at the entry boundary point of each reservoir. The consequence of the choice of this boundary point for the ODEs is to determine using the relation, where 𝑓𝑖 and 𝑉𝑖 represents rate of flow (velocity) of pollutant particles and the water volume in the 𝑖𝑡ℎ reservoir, respectively. It turns out that for a given reservoir is a constant defined at the commencement of the flow process of pollutant particles in solution for that particular reservoir. This specification of executes the mixing process of pollutant particles in water flow as “instantaneous”, so that particle concentration in all other points in a reservoir is uniform without any time elapse. The coefficient remains the serial dilution factor in each reservoir and was computed using relation (24).

Graphical displays of time series solutions for each reservoir in the cascade were demonstrated in comparison for both model equation (8) and ODE equation (27) in the following section.

###### 3.3.5. Comparison of Results and Discussion

When the model equation was ran using MATLAB for the data sets for where obtained from River Nyakomisaro, it produced graphical views that varied from the paired parameter values for β, 𝜕, and 𝜏 in each reservoir, downstream. From the graphs 3.1(a), 3.1(b), and 3.1(c), it was observed that as the values of the parameters and become smaller, the curve is less clear and has shorter times to cover than when the parameter values are larger.

This also agreed with the pattern of the natural phenomenon in mixing-problems that the rate of particle flow decreases with decreased concentration, a condition described as “concentration attenuation”. Therefore, particle migration is faster and more feasible at early stages (front reservoirs) of pollutant incursion into water system than when particles must have stayed (rare reservoirs) in the system.

Further, it was observed from the graphs in Figures 3(a)3(c) that as the values of the parameters and were varied downstream in the river flow, the resulting curves change pattern. The dips (troughs) of the curves change gradually and eventually to linear with zero slopes.

From the graphs given in Figures 4(a)4(d), it was observed that for the same reservoir, under similar physical properties, an ODE curve has shorter times to cover than its DDE curve counterpart. It implies that if the time delays in the mixing processes are not accounted for in particle (pollutant) transport in water reservoirs, and since mixing processes are never instantaneous, then the time required for transport of a species of pollutant in the water reservoir is underestimated. The time series solutions that include delays generated by model equation (8) offers a better approximate value of the time necessary for particle flow in cascade subsurface water reservoirs than mathematical models without delays [6].

#### 4. Summary and Conclusions

(i)The numerical model formulated in this study describes the kind of time delays observed in mixing-problems that was applied to govern the dynamics of water pollution in subsurface water reservoirs. The results from the validation of the model, displayed in Figures (3) and (4), imply that if the average time elapsing for mixing-process of pollutant particles within a water reservoir are accounted for (in models with time delays), and if the duration of the particle movement within the reservoir up to their possible maximum accumulation in water systems are approximately negative exponential in distribution, then such a model is more accurate (than models without delays) to estimate time evolution of particle transport in water reservoirs(ii)The model was applied on three stretch sections on River Nyakomisaro in Kisii town, labeled “SECTION A, SECTION B and SECTION C”. These three sections were identified by recognizing possible potential pollutants entering the river at particular points (point sources). The three stretches were of equal lengths and each was partitioned into six cascade zones. Physical parameters of the river at these sections were measured and the coefficients of model equation (8) were evaluated. It turned out that the graphical views generated displayed similar tendencies (characteristics) even though initial values of concentration were different for every potentially recognized single species pollutant chosen for each section. The graphical comparisons for three sections with different starting points (but with constant delays) were given in Figures (3).(iii)The comparisons of results displayed in Figure 4 between the generated ODEs and the corresponding DDEs in modeling rate of particle flow in reservoirs of water flow systems, indicate the essentiality of invoking the time delays in water quality simulation models such as Qual2kw and CeQualw2 that make use of the finite-difference numerical methods for the solution of the adjective-dispersion mass transport and reaction equations. Time delays are inherent in water quality variations, and therefore, for accurate model simulation of the transport and fate of several constituents in water flow systems, time delays should be considered in the advection-diffusion equation.

#### 5. Recommendations

More validation of the formulated model (equation (8)) to different natural river flow systems is recommended.

#### Nomenclature

 : The gravitational constant (9.8 m/s2) : Rate of pollutant particle flow in water : Discrete Time delay (Dimensionless) : Dilution rate of concentration of pollutant particles in water flow (Dimensionless ratio) ODE: Ordinary Differential Equations PDE: Partial Differential Equations DDEs: Delayed Differential Equations FDE: Functional Differential Equation RCDS: Remote Control Dynamic System IVP: Initial-Value-Problem BVP: Boundary-Value-Problem MATLAB: Matrix Laboratory EPA: Environmental Protection Agency SWR: Subsurface Water Reservoirs ppm: Parts Per Million DF: Dilution Factor.

#### Data Availability

No data were used to support this study.

#### Disclosure

A version of this article was presented as manuscript on November 8, 2018, in http://www.preprints.org/ according to the link below: https://www.preprints.org/manuscript/201811.0201/v1.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors acknowledge Dr. Okwoyo M. James ([email protected], School of Mathematics, University of Nairobi) for his contributions as scientific advisor and critical view in the study proposal.