# Formulated Mathematical Model for Delayed Particle Flow in Cascaded Subsurface Water Reservoirs with Validation on River Flow

Authors: Richard Ombaki; Joash Kerongo
Journal: Journal of Applied Mathematics
 (2022)

---

## 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 concentrationxt 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 t would be equal to the average concentration at some earlier instant, t−τ for the delay τ>0. 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 xt versus time t, 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 xt were different for every potentially recognized single species pollutant considered in each river reservoir. Hence, longer values of timet 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.


---

## Body

## 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, twill be equal to the average concentration at some earlier time t−τ, whereτ>0. 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 Figure1 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.Figure 1
Scheme of pollutant propagation in water reservoir from point source (adopted from [5]).Figure 2
Setup for discretization of River flow channel into cascading water reservoirsi=0,1,2,⋯,n is the position of cascade reservoir). Flowing water meets with the pollutant point source of predicted pollutant in the first reservoir, and the mixing process proceeds downstream through the discretized cascades.This paper is organized as follows. First, the background to the study is described. Then, the governing mathematical principles and equations are given in section2, 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(1)Qin−Qout−Qproduction=Qaccumulation.where Q 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])
(2)dxidt=Qi−1Vixi−1−QiVixi−Qab,iVixi+Ei−1Aixi−1−xi+EiAixi−1−xi+WiVi+Si, where xi is the constituent concentration; Qi is flux at reach i; Qab,i is abstraction flux at reach i; Ei is bulk coefficient between adjacent reaches; i−1,i,i+1 are consecutive cascading reaches; Wi are external loading of constituents; and Siare 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, xit 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, xit 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 xit−τ,where τ is the time delayed in the mixing process at a point as the flow of pollutant load proceeds downstream. Consequently, xit−τ/dt 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:
(3)dx1tdt=−k1x1t,(4)dxidt=−ki−1xi−1t−kixit,2≤i≤n,where x is concentration of pollutant in the cascade; n is number of reservoirs in the cascade; and k=f/Vi. Equation (3) indicates that fresh water enters the initial water reservoir so that the particle concentration x1t 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 reservoirs,i−1thand ith 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 xtof particles only between two progressive adjacent cascade reservoirs, the i−1thand ith 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
(5)dxtdt=xt−he−bxt−τ−dxt.Ifb>d, then equation (5)will have nontrivial steady states. The study considered the case of equation (5) when the average death rate d is constant, and the average birth rateb 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, band d, 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].


## 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(1)Qin−Qout−Qproduction=Qaccumulation.where Q 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])
(2)dxidt=Qi−1Vixi−1−QiVixi−Qab,iVixi+Ei−1Aixi−1−xi+EiAixi−1−xi+WiVi+Si, where xi is the constituent concentration; Qi is flux at reach i; Qab,i is abstraction flux at reach i; Ei is bulk coefficient between adjacent reaches; i−1,i,i+1 are consecutive cascading reaches; Wi are external loading of constituents; and Siare 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, xit 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, xit 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 xit−τ,where τ is the time delayed in the mixing process at a point as the flow of pollutant load proceeds downstream. Consequently, xit−τ/dt 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:
(3)dx1tdt=−k1x1t,(4)dxidt=−ki−1xi−1t−kixit,2≤i≤n,where x is concentration of pollutant in the cascade; n is number of reservoirs in the cascade; and k=f/Vi. Equation (3) indicates that fresh water enters the initial water reservoir so that the particle concentration x1t 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 reservoirs,i−1thand ith 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 xtof particles only between two progressive adjacent cascade reservoirs, the i−1thand ith 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
(5)dxtdt=xt−he−bxt−τ−dxt.Ifb>d, then equation (5)will have nontrivial steady states. The study considered the case of equation (5) when the average death rate d is constant, and the average birth rateb 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, band d, 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 xt 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 sizeNt at any given time can be represented in global variables as
(6)dNdt=B−Dt,B,D>0.The constantsB andD 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 Nt 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,
(7)dNdt=BρNt−τ−DNt.This is a delay differential equation(DDE), named Hayes Equation [10]. In particular, by considering the independent variables B,D,ρ and τ, equation (7) is similar to equation (4) in applications.


### 3.2. Model Formulation

The modelDDE 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:
(8)dxtdt=βitxt−τe−Axt−τi−∂itxt,i=1,2,⋯,n,where i is the position of reservoir in the cascades, xt is the concentration of pollutant species in subsurface water reservoir, τi 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 τ1<τ2<,,,<τn, for reservoirs 1,2,3,⋯,n. This specifies equation (8) as DDE with multiple discrete delays. The parameters β,∂ are time-dependent coefficients so that equation (8) is nonautonomous.The coefficientβt is termed as the incursion term and represents pollutant particles that enter the ithreservoir with a rate of flow fi−1, which is the rate of flow in the previous adjacent reservoir. βt is considered a continuous, positive decreasing function, so that the average concentration rate ẋtof 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∂t 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 d 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], d 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 theith reservoir do not mix up instantaneously with water, but rather takes time τiin 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 i+1threservoir. With reference to equation (7), Hayes’ equation, this possibility is represented in the model equation (8) by the mathematical probability ρt=e−Axt−τi. This term is attached to the incursion term βtto mathematically account for the progressive exponentially reducing (decaying) incursion rate of pollutant particles in the cascades. The constant A 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 concentrationxt, so as to specify equation (8) for a solution which is bounded for each water reservoir in the cascade.Definition 1.
Ifdxt/dt=Ft,xt,xt∈ℝn, where xtt+θ,θ∈−∞,0 represent a trajectory of the solution in the past, then F represent a functional operator from ℝ×C1 to ℝ2. Further, dxt/dt=Fxt,xt+θ1,⋯,xt+θn, for θ1<⋯<θn, define F 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,xt∈ℝn.
(9)Xt=xτ:τ≤t.
Initial conditions:∅:−τ,0⟶ℝn.
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, τc 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 whichxt=x¯ is a solution of equation (8) for all t, so that the equilibrium points satisfy
(10)x¯,x¯,⋯,x¯=0.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 βt=βie−Axt−τiand ∂tconstant 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 concentrationxt, solutions remain positive in all time t. 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βtand ∂tbe positive functions. Suppose that there exists an x¯∈x such that
signβx−∂x=−signx−x¯,and β′x¯<∂′x¯. Then, x¯ is a positive steady state, and the trivial steady state is unstable. If β′x¯x¯>−2∂x¯−∂′x¯x¯, then x¯ is linearly stable for all τ. Otherwise, there exists a τi>0 such that x¯ is stable for τi<τc and unstable for τi>τc [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 τiin SWR, and the proof of the theorem follows from the application.
From Theorem (2) and considering the case β>∂, the nontrivial steady state occurs whenever,
(11)β=∂eAx̅,
or(12)x¯=1Alnβ∂.
Then, according to Theorem (2), x¯is stable for allβ<∂ if and only if
(13)ddxβe−Axx=x¯>−∂β.
This is equivalent to the condition thatβ<∂e2, 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 ekt may be chosen and posted in equation (8) to get,
(14)kekt=βekt−τe−Aekt−τ−∂ekt.
Dividing byekt and rearranging gives a transcendental exponential,
(15)k+∂lnβ/∂−1e−kτ+∂=0.
Takingσ=lnβ/∂−1, then
(16)k+∂σe−kτ+∂=0.
The rootsk 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, τ=0 resulting to simple polynomial, k=−∂σ−∂, which is insignificant as far as pollution is concerned. But, if τ>0, then k=k1+k2i, which is a complex root. Taking the real and imaginary parts,
(17)k1=−∂σcosk2τ,k2=−∂σsink2τ.
Squaring and summing up yields(18)k12+k22=∂2σ2,
or,(19)k1=k2σ2−11/2.
Rewriting the real and imaginary parts of the roots of the characteristic equations,(20)cosk2τ=−1/σ<0,sink2τ=σ2−11/2/σ>0.
Thus, there exists a critical delay,τc at which an eigen-value crosses into the right half-plane, k2τ2∈π/2,π and the delay is given by
(21)τc=1/∂σ2−1cos−1−1/σ.


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

The formulated model equation (8) was solved in Matlab dde23 software to obtain time t 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 βi and ∂i in a particular reservoir that were predetermined for the boundedness of the corresponding delay τi 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 timet 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 t−τ.


#### 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 depth0.5 metres (at minimum during dry season) and 1–2 metres (at maximum during high rainfall) and average width below 3.0 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 2.228 ppm, 4.50ppm, and 6.80ppm, respectively.


#### 3.3.1. Determination of Flow Volume

Each reservoir volumeVi was calculated from the length L and average cross-sectional area A for every ith reservoir as
(22)Vi=AiLi,where A=depth×width. 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
(23)τc=LnLT,where Ln=∑i=1nLi,i=1,2,⋯,n.These proportions were necessary to maintain thatτ1<τ2<…<τn, so that the time lags are distinctively discrete and in series. The discrete delays were calculated and bound as shown in Table 1 below.Table 1
Boundedness of delays,τi.
ithreservoirLength of reservoir,LmFor sections “A, B, and C”τi>τc0200.0001500.091>0.0642600.200>0.1663800.345>0.28641000.527>0.49751200.745>0.623The 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, VT (obtained by summing upVi for i=1,2,3,⋯,n,where n is the number of zones in the cascade) to V0, volume of water (assumed to be free of single species pollutant) entering the cascade zone. Thus,
(24)DF=VTV0=1630800400000≈4,where V0is 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 ∂i was computed as the DF of the system. Thus, the DF up to sixth reservoir was reached as:
(25)14,116,164,1256,11024,14096.


#### 3.3.3. Coefficients for Model Equation

Using the results in Table2, the incursion term βitfor every cascade was calculated from the relation
(26)βit=fi−1Vn,Vn=∑i=1nVi,i=1,2,⋯,n.Table 2
Kinematic parameters from River Nyakomisaro.
ithreservoirFlow ratefi−1m/s at sections:Reservoir volumeVim3 at sections:“A and B”“C”“A and B”“C”01.553.20400,000220,00011.502.96536,000145,00021.102.66250,000156,80030.922.35400,000172,00040.801.90381,700182,40050.601.76402,000188,00060.501.62322,000192,000The computed results for the six cascade zones were recorded in Table3. 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.Table 3
Values ofβ, ∂, and initial concentration, x0t.
ithreservoirβi×10−6m2s−1∂itABCx0t2.2284.5006.8000i2.8002.8000.2041/4i+11.3991.3990.1691/16i+20.7760.7760.1361/64i+30.5090.5090.1041/256i+40.3040.3040.0931/1024i+50.2180.2180.0841/4096


#### 3.3.4. Particle Flow forτ=0

The model equation (8) was transformed into its ordinary differential equation counterpart by setting the delays 𝜏𝑖 =0, which resulted to
(27)dxdt=β0xte−Axt−∂ix,where β0 represents the rate per unit volume (flow flux) at which pollutant particles are entering the ith 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β0 using the relation,
(28)β0=fiVi,where 𝑓𝑖 and 𝑉𝑖 represents rate of flow (velocity) of pollutant particles and the water volume in the 𝑖𝑡ℎ reservoir, respectively. It turns out that β0 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 β0 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 ∂i 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 equation3.3 was ran using MATLAB for the data sets for βi,∂i,τi where i=1,2,3,⋯,6 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 Figures3(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.Figure 3
Time series trajectories for particle flow in reservoirs. Figures3(a)–3(c) show plots of Concentration xtversus time t for SECTIONS A, B, and C, respectively, of river flow.
(a)(b)(c)From the graphs given in Figures4(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].Figure 4
Plots of concentrationxt versus time t for comparison DDEs and ODEs in modeling rate of particle flow in reservoirs of water flow systems.
(a)(b)(c)(d)(e)


## 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 xt 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 sizeNt at any given time can be represented in global variables as
(6)dNdt=B−Dt,B,D>0.The constantsB andD 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 Nt 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,
(7)dNdt=BρNt−τ−DNt.This is a delay differential equation(DDE), named Hayes Equation [10]. In particular, by considering the independent variables B,D,ρ and τ, equation (7) is similar to equation (4) in applications.


## 3.2. Model Formulation

The modelDDE 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:
(8)dxtdt=βitxt−τe−Axt−τi−∂itxt,i=1,2,⋯,n,where i is the position of reservoir in the cascades, xt is the concentration of pollutant species in subsurface water reservoir, τi 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 τ1<τ2<,,,<τn, for reservoirs 1,2,3,⋯,n. This specifies equation (8) as DDE with multiple discrete delays. The parameters β,∂ are time-dependent coefficients so that equation (8) is nonautonomous.The coefficientβt is termed as the incursion term and represents pollutant particles that enter the ithreservoir with a rate of flow fi−1, which is the rate of flow in the previous adjacent reservoir. βt is considered a continuous, positive decreasing function, so that the average concentration rate ẋtof 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∂t 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 d 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], d 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 theith reservoir do not mix up instantaneously with water, but rather takes time τiin 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 i+1threservoir. With reference to equation (7), Hayes’ equation, this possibility is represented in the model equation (8) by the mathematical probability ρt=e−Axt−τi. This term is attached to the incursion term βtto mathematically account for the progressive exponentially reducing (decaying) incursion rate of pollutant particles in the cascades. The constant A 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 concentrationxt, so as to specify equation (8) for a solution which is bounded for each water reservoir in the cascade.Definition 1.
Ifdxt/dt=Ft,xt,xt∈ℝn, where xtt+θ,θ∈−∞,0 represent a trajectory of the solution in the past, then F represent a functional operator from ℝ×C1 to ℝ2. Further, dxt/dt=Fxt,xt+θ1,⋯,xt+θn, for θ1<⋯<θn, define F 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,xt∈ℝn.
(9)Xt=xτ:τ≤t.
Initial conditions:∅:−τ,0⟶ℝn.
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, τc 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 whichxt=x¯ is a solution of equation (8) for all t, so that the equilibrium points satisfy
(10)x¯,x¯,⋯,x¯=0.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 βt=βie−Axt−τiand ∂tconstant 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 concentrationxt, solutions remain positive in all time t. 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βtand ∂tbe positive functions. Suppose that there exists an x¯∈x such that
signβx−∂x=−signx−x¯,and β′x¯<∂′x¯. Then, x¯ is a positive steady state, and the trivial steady state is unstable. If β′x¯x¯>−2∂x¯−∂′x¯x¯, then x¯ is linearly stable for all τ. Otherwise, there exists a τi>0 such that x¯ is stable for τi<τc and unstable for τi>τc [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 τiin SWR, and the proof of the theorem follows from the application.
From Theorem (2) and considering the case β>∂, the nontrivial steady state occurs whenever,
(11)β=∂eAx̅,
or(12)x¯=1Alnβ∂.
Then, according to Theorem (2), x¯is stable for allβ<∂ if and only if
(13)ddxβe−Axx=x¯>−∂β.
This is equivalent to the condition thatβ<∂e2, 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 ekt may be chosen and posted in equation (8) to get,
(14)kekt=βekt−τe−Aekt−τ−∂ekt.
Dividing byekt and rearranging gives a transcendental exponential,
(15)k+∂lnβ/∂−1e−kτ+∂=0.
Takingσ=lnβ/∂−1, then
(16)k+∂σe−kτ+∂=0.
The rootsk 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, τ=0 resulting to simple polynomial, k=−∂σ−∂, which is insignificant as far as pollution is concerned. But, if τ>0, then k=k1+k2i, which is a complex root. Taking the real and imaginary parts,
(17)k1=−∂σcosk2τ,k2=−∂σsink2τ.
Squaring and summing up yields(18)k12+k22=∂2σ2,
or,(19)k1=k2σ2−11/2.
Rewriting the real and imaginary parts of the roots of the characteristic equations,(20)cosk2τ=−1/σ<0,sink2τ=σ2−11/2/σ>0.
Thus, there exists a critical delay,τc at which an eigen-value crosses into the right half-plane, k2τ2∈π/2,π and the delay is given by
(21)τc=1/∂σ2−1cos−1−1/σ.


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

The formulated model equation (8) was solved in Matlab dde23 software to obtain time t 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 βi and ∂i in a particular reservoir that were predetermined for the boundedness of the corresponding delay τi 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 timet 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 t−τ.


### 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.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, τc 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 whichxt=x¯ is a solution of equation (8) for all t, so that the equilibrium points satisfy
(10)x¯,x¯,⋯,x¯=0.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 βt=βie−Axt−τiand ∂tconstant 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 concentrationxt, solutions remain positive in all time t. 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βtand ∂tbe positive functions. Suppose that there exists an x¯∈x such that
signβx−∂x=−signx−x¯,and β′x¯<∂′x¯. Then, x¯ is a positive steady state, and the trivial steady state is unstable. If β′x¯x¯>−2∂x¯−∂′x¯x¯, then x¯ is linearly stable for all τ. Otherwise, there exists a τi>0 such that x¯ is stable for τi<τc and unstable for τi>τc [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 τiin SWR, and the proof of the theorem follows from the application.
From Theorem (2) and considering the case β>∂, the nontrivial steady state occurs whenever,
(11)β=∂eAx̅,
or(12)x¯=1Alnβ∂.
Then, according to Theorem (2), x¯is stable for allβ<∂ if and only if
(13)ddxβe−Axx=x¯>−∂β.
This is equivalent to the condition thatβ<∂e2, 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 ekt may be chosen and posted in equation (8) to get,
(14)kekt=βekt−τe−Aekt−τ−∂ekt.
Dividing byekt and rearranging gives a transcendental exponential,
(15)k+∂lnβ/∂−1e−kτ+∂=0.
Takingσ=lnβ/∂−1, then
(16)k+∂σe−kτ+∂=0.
The rootsk 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, τ=0 resulting to simple polynomial, k=−∂σ−∂, which is insignificant as far as pollution is concerned. But, if τ>0, then k=k1+k2i, which is a complex root. Taking the real and imaginary parts,
(17)k1=−∂σcosk2τ,k2=−∂σsink2τ.
Squaring and summing up yields(18)k12+k22=∂2σ2,
or,(19)k1=k2σ2−11/2.
Rewriting the real and imaginary parts of the roots of the characteristic equations,(20)cosk2τ=−1/σ<0,sink2τ=σ2−11/2/σ>0.
Thus, there exists a critical delay,τc at which an eigen-value crosses into the right half-plane, k2τ2∈π/2,π and the delay is given by
(21)τc=1/∂σ2−1cos−1−1/σ.


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

The formulated model equation (8) was solved in Matlab dde23 software to obtain time t 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 βi and ∂i in a particular reservoir that were predetermined for the boundedness of the corresponding delay τi 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 timet 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 t−τ.


## 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 depth0.5 metres (at minimum during dry season) and 1–2 metres (at maximum during high rainfall) and average width below 3.0 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 2.228 ppm, 4.50ppm, and 6.80ppm, respectively.


### 3.3.1. Determination of Flow Volume

Each reservoir volumeVi was calculated from the length L and average cross-sectional area A for every ith reservoir as
(22)Vi=AiLi,where A=depth×width. 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
(23)τc=LnLT,where Ln=∑i=1nLi,i=1,2,⋯,n.These proportions were necessary to maintain thatτ1<τ2<…<τn, so that the time lags are distinctively discrete and in series. The discrete delays were calculated and bound as shown in Table 1 below.Table 1
Boundedness of delays,τi.
ithreservoirLength of reservoir,LmFor sections “A, B, and C”τi>τc0200.0001500.091>0.0642600.200>0.1663800.345>0.28641000.527>0.49751200.745>0.623The 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, VT (obtained by summing upVi for i=1,2,3,⋯,n,where n is the number of zones in the cascade) to V0, volume of water (assumed to be free of single species pollutant) entering the cascade zone. Thus,
(24)DF=VTV0=1630800400000≈4,where V0is 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 ∂i was computed as the DF of the system. Thus, the DF up to sixth reservoir was reached as:
(25)14,116,164,1256,11024,14096.


### 3.3.3. Coefficients for Model Equation

Using the results in Table2, the incursion term βitfor every cascade was calculated from the relation
(26)βit=fi−1Vn,Vn=∑i=1nVi,i=1,2,⋯,n.Table 2
Kinematic parameters from River Nyakomisaro.
ithreservoirFlow ratefi−1m/s at sections:Reservoir volumeVim3 at sections:“A and B”“C”“A and B”“C”01.553.20400,000220,00011.502.96536,000145,00021.102.66250,000156,80030.922.35400,000172,00040.801.90381,700182,40050.601.76402,000188,00060.501.62322,000192,000The computed results for the six cascade zones were recorded in Table3. 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.Table 3
Values ofβ, ∂, and initial concentration, x0t.
ithreservoirβi×10−6m2s−1∂itABCx0t2.2284.5006.8000i2.8002.8000.2041/4i+11.3991.3990.1691/16i+20.7760.7760.1361/64i+30.5090.5090.1041/256i+40.3040.3040.0931/1024i+50.2180.2180.0841/4096


### 3.3.4. Particle Flow forτ=0

The model equation (8) was transformed into its ordinary differential equation counterpart by setting the delays 𝜏𝑖 =0, which resulted to
(27)dxdt=β0xte−Axt−∂ix,where β0 represents the rate per unit volume (flow flux) at which pollutant particles are entering the ith 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β0 using the relation,
(28)β0=fiVi,where 𝑓𝑖 and 𝑉𝑖 represents rate of flow (velocity) of pollutant particles and the water volume in the 𝑖𝑡ℎ reservoir, respectively. It turns out that β0 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 β0 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 ∂i 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 equation3.3 was ran using MATLAB for the data sets for βi,∂i,τi where i=1,2,3,⋯,6 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 Figures3(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.Figure 3
Time series trajectories for particle flow in reservoirs. Figures3(a)–3(c) show plots of Concentration xtversus time t for SECTIONS A, B, and C, respectively, of river flow.
(a)(b)(c)From the graphs given in Figures4(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].Figure 4
Plots of concentrationxt versus time t for comparison DDEs and ODEs in modeling rate of particle flow in reservoirs of water flow systems.
(a)(b)(c)(d)(e)


## 3.3.1. Determination of Flow Volume

Each reservoir volumeVi was calculated from the length L and average cross-sectional area A for every ith reservoir as
(22)Vi=AiLi,where A=depth×width. 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
(23)τc=LnLT,where Ln=∑i=1nLi,i=1,2,⋯,n.These proportions were necessary to maintain thatτ1<τ2<…<τn, so that the time lags are distinctively discrete and in series. The discrete delays were calculated and bound as shown in Table 1 below.Table 1
Boundedness of delays,τi.
ithreservoirLength of reservoir,LmFor sections “A, B, and C”τi>τc0200.0001500.091>0.0642600.200>0.1663800.345>0.28641000.527>0.49751200.745>0.623The 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, VT (obtained by summing upVi for i=1,2,3,⋯,n,where n is the number of zones in the cascade) to V0, volume of water (assumed to be free of single species pollutant) entering the cascade zone. Thus,
(24)DF=VTV0=1630800400000≈4,where V0is 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 ∂i was computed as the DF of the system. Thus, the DF up to sixth reservoir was reached as:
(25)14,116,164,1256,11024,14096.


## 3.3.3. Coefficients for Model Equation

Using the results in Table2, the incursion term βitfor every cascade was calculated from the relation
(26)βit=fi−1Vn,Vn=∑i=1nVi,i=1,2,⋯,n.Table 2
Kinematic parameters from River Nyakomisaro.
ithreservoirFlow ratefi−1m/s at sections:Reservoir volumeVim3 at sections:“A and B”“C”“A and B”“C”01.553.20400,000220,00011.502.96536,000145,00021.102.66250,000156,80030.922.35400,000172,00040.801.90381,700182,40050.601.76402,000188,00060.501.62322,000192,000The computed results for the six cascade zones were recorded in Table3. 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.Table 3
Values ofβ, ∂, and initial concentration, x0t.
ithreservoirβi×10−6m2s−1∂itABCx0t2.2284.5006.8000i2.8002.8000.2041/4i+11.3991.3990.1691/16i+20.7760.7760.1361/64i+30.5090.5090.1041/256i+40.3040.3040.0931/1024i+50.2180.2180.0841/4096


## 3.3.4. Particle Flow forτ=0

The model equation (8) was transformed into its ordinary differential equation counterpart by setting the delays 𝜏𝑖 =0, which resulted to
(27)dxdt=β0xte−Axt−∂ix,where β0 represents the rate per unit volume (flow flux) at which pollutant particles are entering the ith 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β0 using the relation,
(28)β0=fiVi,where 𝑓𝑖 and 𝑉𝑖 represents rate of flow (velocity) of pollutant particles and the water volume in the 𝑖𝑡ℎ reservoir, respectively. It turns out that β0 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 β0 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 ∂i 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 equation3.3 was ran using MATLAB for the data sets for βi,∂i,τi where i=1,2,3,⋯,6 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 Figures3(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.Figure 3
Time series trajectories for particle flow in reservoirs. Figures3(a)–3(c) show plots of Concentration xtversus time t for SECTIONS A, B, and C, respectively, of river flow.
(a)(b)(c)From the graphs given in Figures4(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].Figure 4
Plots of concentrationxt versus time t for comparison DDEs and ODEs in modeling rate of particle flow in reservoirs of water flow systems.
(a)(b)(c)(d)(e)


## 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 xt 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 Figure4 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.



---

*Source: 3438200-2022-11-10.xml*