# A Mathematical Model for Variable Chlorine Decay Rates in Water Distribution Systems

Authors: Ababu T. Tiruneh; Tesfamariam Y. Debessai; Gabriel C. Bwembya; Stanley J. Nkambule
Journal: Modelling and Simulation in Engineering
 (2019)

---

## Abstract

A model for relating the time-dependent variable rate of reaction to the decay of chlorine residual in water is developed based on the initial chlorine dose, molar concentrations of reactants, and the rate constant itself. The mathematical model, while retaining its second order nature, simplifies the solution as the residual chlorine and aggregate parameters such as molar concentration of reactants can be estimated. The model is based on molar-averaged reaction rates involving arithmetic and harmonic means of reactants that eliminate the individual reaction rates that are difficult to determine. Part of the mathematical assumption used in the derivation of the equations using molar averaging is tested for its validity through theoretical as well as Monte Carlo simulation of the error term over wide ranges of assumed reaction rates and molar concentration of reactants. The second-order variation of the rate of reaction with respect to the initial chlorine concentration has been verified through experimental tests of bulk chlorine decay carried out at different chlorine doses.


---

## Body

## 1. Introduction

Disinfection of water using chlorine is an essential water treatment process step that renders the water supplied to consumers bacteriologically safe [1]. The amount of chlorine added for disinfection is controlled in such a way that there is an adequate residual chlorine present as the water flows through the distribution system until it reaches the consumers. The presence of a minimum residual chlorine guarantees that recontamination does not occur as well as deterioration in the aesthetic quality of water due to growth of organisms within the distribution system does not occur. Excess chlorination is chlorination of water over and above the minimum dose that is required to ensure potable water quality. Practices of excess chlorination give rise to the formation of disinfection by-products (DBPs) in water, which are compounds known to be associated with health risks related to cancer [2].The problem related to maintenance of residual chlorine is increased by water supply interruptions that lead to formation of stagnant water, which, during resumption of supply, may be drawn with no chlorine residual left as the water reaches the consumers [3]. Old water supply systems with pipes that have deteriorated linings encourage microbial growth in the distribution pipe, which in turn results in rapid loss of chlorine due to the wall decay reaction [4].Chlorine is added for disinfection of water in an optimal way depending on the water quality of the treated water prior to disinfection due to the presence of reactants that exert chlorine demand in the bulk water as the water travels through the distribution system. In addition, the dosage of chlorine is also influenced by chlorine consuming reactants present on the walls of the water supply distribution pipes. An optimum dosage is desirable that addresses these consumptions of chlorine in such a way that a low dosage is avoided that does not ensure adequate disinfection of the water supplied while at the same time minimising the instances of excess chlorine that encourage the formation of undesirable disinfection by-products such as trihalomethanes [5].The reduction of chlorine residual in water supply system take place in the bulk water as well as the walls of pipes and surfaces such as storage tanks present along the distribution system. The bulk decay is a volume-based reduction process whereas the wall decay is a surface area based process. In addition, the nature of reaction and the types of reactants involved in these two processes are different. Water quality models for the reduction of chlorine, therefore, require separate steps for determination of the wall and bulk decay coefficients [5]. The bulk decay is determined by the laboratory batch test on sample of water taken from the water treatment system ready for disinfection. The wall decay rate is often determined by a calibration process as a difference between the chlorine consumption observed in the distribution system and the chlorine consumption due to bulk decay alone determined by laboratory tests [6]. The wall decay depends on pipe conditions including the materials from which pipes are made [7]. In general, laboratory and pilot tests alone cannot adequately represent the chlorine decay process in the actual distribution system. Therefore, the decay model parameters in water quality modeling programs such as the EPANET need to be calibrated against actual observation of chlorine residuals within the distribution systems [8].Because of the greater domain of space and time through which chlorine residual need to be determined within the distribution system, monitoring of chlorine residuals based on laboratory determination alone is impractical and cost-ineffective [9]. Mathematical models have been developed that trace the decay of chlorine using conservation of mass equation along the travel path of the water in the distribution system. The conservation of mass equation in space and time takes into account the transport of chlorine with the bulk water (advection) and its reduction in the bulk water as well as along the pipe walls as the water travels through the distribution system pipes [10].Traditional chlorine decay models are simple first-order reactions [2, 3]. The reaction rate constant in the traditional first-order model is an apparent rate constant in which the molar concentration of the reactants exerting chlorine demand is implicitly represented. Because of this, the reaction rate constant changes when the nature of reactants such as the concentration of dissolved organic matter changes [11]. The rate of reaction is also variable based on the concentration of chlorine that is initially applied [12]. The rate of reaction is also influenced by the presence of different reactants with heterogeneous kinetics of reaction. Initially, chlorine is consumed faster by fast reacting species such as organic matter and by species having greater molarity. This is followed by slowly reacting species and species with lower molar concentration present in the chlorinated water.In an attempt to take into account the effect of different reactants as well as concentration of chlorine applied, a number of different models have been proposed that more or less deviate from the implicit first-order decay rate model. One such model is a parallel two-reaction decay rate involving fast and slow reactants [13]. A number of other variations of such models are also proposed that are second-order models or a hybrid of first- and second-order models [13–17]. However, few researchers played down the importance deviating from the traditional pseudo-first-order reaction in which they argued that the difference in model results in terms of the free residual chlorine present is not of great practical significance [18–20].The dependence of the rate of reaction for chlorine decay on the initial concentration of chlorine applied has been cited by several researchers [7, 12, 13, 18, 21]. In addition, researchers have addressed several relevant factors such as temperature and oxidisable organic as well as inorganic matter. The model by Hua et al. [13] integrates these several factors in determining the chlorine residual. Hallam et al. [7] found inverse relation between the initial chlorine dose and the rate of reaction. This relationship has also been confirmed by a number of other researchers [13, 18, 22, 23].Different researchers have approached the effect of the initial chlorine present and of the different reactants exerting chlorine demand differently. Some researchers such as Fisher et al. [24] assumed two-reaction model in which the two reactions involving fast and slow reactants act in parallel. The model is calibrated against ranges of minimum and maximum initial chlorine concentrations expected. Another model by Kastl et al. [25] tried to determine the model parameters for initial chlorine varying between 1 and 4 mg/L. According to these researchers, the model parameters remain constant in this calibrated range.The model developed by Phillip et al. [5] is a more explicit second-order time-dependent model. The rate of reaction varies second order in time and is also varying directly with the chlorine concentration present. The proposed solution is a trial and error procedure involving solving two differential equations, namely, the equation involving reduction of chlorine and the time-dependent rate of reaction equation. The boundary value problem is solved using Euler’s method. Four parameters need calibration for this model as demanded by the two differential equations. This is done through the bulk decay data in which the optimal values of the parameters are determined by minimising the error between the model result and actual chlorine residual experimentally observed through evolutionary algorithms.


## 2. Materials and Methods

This research paper provides outline of the mathematical model for variable rate of decay of chlorine that is also accompanied by laboratory experiment for verification of the second order model with respect to the initial chlorine dose used. The mathematical model considers two cases. In the first case, the molar-averaged reaction rate constant variation with time is modeled. In the second case, the overall reaction rate constant in which the molar concentration of reactants are in-built in the rate constant is modeled.


### 2.1. Experimental Determination of Chlorine Residual

For the experiment involving determination of chlorine residual at different times and under different initial chlorine dosages used, samples of water that underwent conventional treatment up to and including rapid sand filtration and just prior to disinfection were collected form the Matsapha water treatment plant and were subsequently used in the experimental trials. The samples collected were brought to the Chemistry Laboratory of the University of Eswatini. Different dosages of chlorine were added to multiple samples, and the chlorine residuals were determined at predetermined intervals of time. The method used for determination of residual chlorine was the iodometric titration as stated in Standard Methods for the Examination of Water and Wastewater [26]. The method is based on the principle that when potassium iodide is added to a sample of water containing residual chlorine at pH of 8 or less, the residual chlorine liberates iodine from the potassium iodide and the liberated iodine is titrated with sodium thiosulphate (Na2S2O3) titrant.In order to make the titration stoichiometric, the pH of the sample is reduced to between 3 and 4 by adding acetic acid. The method has expected a detection limit of 0.04 mg/L if 0.001 N sodium thiosulphate titrant is used. All chemicals used had reagent-grade quality, and fresh solutions were prepared every time the experiment was repeated at different time periods.Depending upon the concentration anticipated, suitable sample volume was taken so that the titrant sodium thiosulphate (Na2S2O3) volume would not exceed 20 mL and the starch indicator volume also is above 0.2 mL. To the sample volume, 5 mL acetic acid was added followed by 1 g of potassium iodide powder measured on weighing balance. After this, the sample was titrated with sodium thiosulphate (Na2S2O3) of appropriate normality. For samples with low anticipated chlorine concentrations, the titrant was prepared to be of low normality. The titration was continued until the yellow colour was almost disappearing. Then, 1 mL of the starch solution was added until the blue colour disappeared. In order to compensate for the method error, blank titration was also performed using distilled-deionised water and that passed through the same procedure as the one used for the actual sample. A minimum of three repetitions were performed for each determination. After titration, mg/L of chlorine residual present in the sample was determined using the following formula:(1)mgLof chlorine=A±B∗N∗35450mL of sample,where A is the volume of the titrant used for the sample, B is the volume of the titrant used for blank, and N is the normality of the Na2S2O3 titrant.


### 2.2. Basis of Mathematical Model

The time varying reaction rate constant equation as is suggested by Phillip et al. [5] is used for the mathematical formulation. According to this equation, the concentration-averaged rate of reaction variation with time is given by the following expression:(2)dktdt=Ctkt2−∑i=1Nki2XitXT.The original second-order variation of chlorine decay rate with molar concentration of chlorine and reactants and from which equation (2) was derived is given by(3)dCtdt=−ktCtXt=dXTdt.In Equations (2) and (3), kt is the molar weighed reaction rate, Ct is the chlorine residual measured at time t, Xit is the molar concentration of reactant i measured at time t and reacting with chlorine, XT is the total molarity of the reactants exerting chlorine demand, and N is the total number of chlorine consuming reactants present in the water.Meanwhile, the molar-averaged reaction ratekt is defined as(4)kt=∑i=1NkiXit∑i=1NXit=∑i=1NkiXitXT,where ki is the individual rate of reaction of reactant i having molar concentration xi and is exerting chlorine demand. Other symbols are as defined for Equation (2). Equation (2) indicates the second-order nature of the variation of the molar-averaged reaction rate constant kt with time. It is also clear from Equation (2) that other factors staying constant, the reaction rate constant kt varies as a first-order reaction with respect to the chlorine residual present. The aggregate effect of the presence of chlorine residual and this second-order variation of the rate constant are because as the chlorine concentration increases the rate constant tends to decrease. Similar observations were noted by several researchers through empirical experimental determination of the relationship between the reaction rate constant and the initial chlorine used in the experiments [13, 18, 22, 23]. The rate constant appears to decrease inversely with increase in the chlorine dose used in the experiments. Phillip et al. [5] also argue that the rate constant variation with time always stays negative in Equation (2) as the expression in bracket in this equation can be proven either negative or zero only. In other words, it cannot take positive values.Equation (2) as developed by Phillip et al. [5] is extended in this paper to provide a formula that relates the time variation of rate of reaction kt with chlorine residual present and the molar concentration of reactants. The difficulty of dealing with the indeterminate term in the bracket on the right-hand side of Equation (2) and involving the individual reactants can be resolved through mathematical averaging of the product of the individual reaction rates and the molar concentration of reactants. This averaging is later verified by theoretical formulation of the error term and Monte Carlo simulation of the error of averaging under different values of these product terms, i.e., molar concentration of reactants and individual reaction rates.Starting with the expression in bracket in Equation (2) and reformulating it such that(5)∑i=1Nki2Xit=∑i=1Nki2Xit2Xit=∑i=1NkiXit2Xit.Rearranging the right-hand side of the above expression, we get(6)∑i=1NkiXit2Xit=∑kiXit.kiXitXit.Now defining the average ofkiXit such that(7)kiXitav=∑kiXitN.Since it was already defined in Equation (4) that(8)kt=∑kiXitXT.Equations (4) and (7) are combined to give(9)kiXitav=∑kiXitN=ktXTN.Now returning to Equation (5) once again, we get(10)∑i=1NkiXit2Xit=∑kiXit·kiXitXit=kiXitav·∑kiXitXit.In the above expression, (kiXit/Xit) is used as a weighing factor. This is also analogous to the weighing of a continuous function, i.e.,(11)∫x.fxdx=x¯∫fxdx.In addition, by using similar averaging method, we get(12)∑kiXitXit=kiXitav·∑1Xit.In this case, 1/Xit are used as weighing factors.Overall,(13)∑i=1NkiXit2Xit=kiXitav2.∑1Xit=kt2XT2N2.∑1Xit.The use of weighted average in Equation (13) above for discrete values can be verified theoretically as well as practically using Monte Carlo simulation using random variation of individual reactants’ rate constants and reactants’ molar concentrations as described below. First, the theoretical basis is explained in the following:(14)∑i=1Nkx2≈kxav·∑i=1NX=1N∑i=1Nkx∑i=1NX.Working a little with the right-hand side expression, we get(15)1N∑i=1Nkx∑i=1NX=1N∑i=1Nkx2+∑i=1N∑j=1,j≠iNkixixj.Now defining the difference term in molar expression of reactants Δxij, we get(16)Δxij=xi−xj,so thatxj=xi−Δxij.Substituting the above expression for thexj term, we get(17)1N∑i=1Nkx∑i=1NX=1N∑i=1Nkx2+∑i=1N∑j=1,j≠iNkixixi−Δxij.The right-hand side of the above expression changes to(18)1N∑i=1Nkx∑i=1NX=1N∑i=1Nkx2+N−1∑i=1Nkx2−∑i=1N∑j=1,j≠iNkixiΔxij,1N∑i=1Nkx∑i=1NX=∑i=1Nkx2−1N∑i=1N∑j=1,j≠iNkixiΔxij.Because of symmetry, it is easy to show that(19)Δxij=−Δxji.Using this symmetry, the expression is reduced further to(20)1N∑i=1Nkx∑i=1NX=∑i=1Nkx2−1N∑i=1N∑j=i+1Nkixi−kjxjΔxij.Defining the difference terms Δkx such that(21)Δkxij=kixi−kjxj,and using this term in the above equation, we get(22)1N∑i=1Nkx∑i=1NX=∑i=1Nkx2−1N∑i=1N∑j=i+1NΔkxijΔxij.Therefore,(23)∑i=1Nkx2=1N∑i=1Nkx∑i=1NX+E,where the error term E is defined as(24)E=1N∑i=1N∑j=i+1NΔkxijΔxij.The above error term is summation over second-order difference (i.e., Δkx times Δx) and tends to be small compared with other terms in the equation containing E:The second practical verification is using Monte Carlo simulation of the error term and finding the average as well as bracketing the range of variation of the error termE. For this purpose, it is assumed that there were five chlorine demanding reactants with individual reaction rates ki varying by a minimum factor of 10 and maximum factor of 1000 and with molar concentrations of reactants xi also varying by a minimum factor of 10 and maximum factor of 1000. The Monte Carlo simulation is done for each range of variation 10000 times, with the following results shown in Table 1.Table 1
Results of Monte Carlo simulation: error of estimationE for different ranges of molar concentration of reactants and reaction rates.
K
10, X10
K
20, X20
K
50, X50
K
100, X100
K
200, X200
K
500, X500
K
1000, X1000
Average error,E
12.68
12.28
12.53
12.53
12.63
12.44
12.61
Standard deviation
8.68
8.67
8.68
8.76
8.79
8.73
8.84
Geometric mean of errors,E
8.98
NA∗
NA∗
NA∗
8.85
8.79
NA
Harmonic mean of errors,E
3.25
1.54
3.12
0.65
0.03
2.83
2.64
Maximum error,E
54.99
54.73
56.69
56.73
59.75
50.56
59.70
Minimum error,E
0.02
0.00∗
0.02
0.00∗
0.00∗
0.01
0.01
∗Geometric mean calculation invalidated by presence of zero in the data as is noticed in the last row of this table.As is shown in Table1, the mean error E is about 12% with 95% confidence range lying between 12.11 and 12.88 which is quite narrow range. The error term is therefore comparatively small, and the use of the weighing factor in the derivation of the rate equation is justified. The same can be said of the geometric and harmonic means of the error term E which are in the order of 8% and 3%, respectively, as shown in Table 1. These data support the heuristic use of the averaging used in deriving Equation (13) in terms of the ranges of variations of the kinetic rate and molar concentration of reactants.Now returning back to Equation (13) and using the arithmetic average (XA) and harmonic average (XH) of the molar concentration of the reactants exerting chlorine demand, we get(25)X¯A=XTN,X¯H=N∑i=1N1/Xit.The expression in Equation (13) can now be written as(26)∑i=1NkiXit2Xit=kt2XT2N2·∑1Xit=kt2XTXTN·1N·∑1Xit=kt2XTX¯ANNX¯H.Simplifying further, we get(27)∑i=1Nki2Xit=∑i=1NkiXit2Xit=kt2XTX¯AX¯H.Now, the original expression given in Equation (2) is(28)dktdt=Ctkt2−∑i=1Nki2XitXT,which can be written as(29)dktdt=Ctkt2−∑i=1Nki2XitXT=Ctkt2−kt2XTX¯A/X¯HXT,which is the same as(30)dktdt==Ctkt21−X¯AX¯H.Changing the order of the expression in bracket, we obtain(31)dktdt=−Ctkt2X¯AX¯H−1.The above expression in Equation (31) shows that the reaction rate is second order with respect to kt and first order with Ct. The reaction rate decreases with time being always negative since the expression in bracket is always positive (the arithmetic mean is always greater than the harmonic mean).Equation (31) as such replaces the individual reaction rates by the concentration-weighed aggregate reaction rate kt as it is the one that can be easily measured. The concentration-weighed aggregate reaction rate kt is influenced by the ratio of arithmetic to harmonic mean of reactants. If the molar concentration of reactants is close in magnitude to each other, then(32)X¯AX¯H≈1,X¯AX¯H−1⟶0.The ratekt becomes constant (dkt/dt = 0 in Equation (31)), and the chlorine decay model approaches first order rate reaction. By contrast, if there is significant variation in the molar concentration among the reactants, the ratio(33)X¯AX¯H−1≫0,holds true. The chlorine decay follows second-order reaction and is highly influenced by the reactant with greater molar concentration. Because of the validity of the expression(34)X¯AX¯H−1≥0,as the arithmetic mean is always greater than the harmonic mean, the time rate of decay of the concentration-weighed reaction rate is always negative and is decreasing with time. This is also observed in the literature through empirical inverse relationship between the chlorine decay rate and initial chlorine dose obtained experimentally. It is also clear from Equation (31) that the rate kt decreases faster with increase in initial chlorine dose. The implication is that the chlorine decay curve becomes flatter (second order) at higher initial chlorine dosages. In other words, within the distribution system, it takes longer time for the chlorine to decay compared to first-order rate decay and the residuals tend to be greater.


### 2.3. Procedure for Determination of the Chlorine Decay Rate Based on Equations (31) and (36)

In order to determine the chlorine decay rate, it is necessary to solve the differential equations involving equation (31) together with equation (36), namely,(35)dktdt=−Ctkt2X¯AX¯H−1,(36)dCtdt=−ktCtXt=dXTdt.To determine the decay rate from the above two equations, it is necessary to determine theXt and kt values by alternating between these two equations. In other words, the Xt and Ct values are determined from equation (36) and the kt value from equation (31). To begin with, however, the initial rates k0 and X0 as well as the ratio of arithmetic to harmonic mean of the reactants need to be determined. It is proposed that the initial k0 value as well as the expression in equation (31), i.e.,(37)X¯AX¯H−1,can be determined by calibration. This means assuming different values of k0 and X¯A/X¯H and choosing the values that give minimum sum of squares of error between the observed and calculated chlorine residual concentrations measured at different times. It is possible that the combination of k0 and X¯A/X¯H can be determined using the evolutionary algorithm or Monte Carlo method.It is also possible to relate theX¯A/X¯H value to reactant chemical species parameters involving, for example, TOC and UV absorption values. One possibility is using the TOC/UV ratio to relate with X¯A/X¯H and use this relationship to model the variable decay rate kt as given by equation (31).It is also proposed that the initial reactant concentrationX0 be the molar equivalent of the total organic carbon (TOC) concentration. The actual value of X0 can also be ascertained through dilution experiment by running a chlorine decay experiment at successive dilution and extrapolating the first order rate constant corresponding to the undiluted sample. The X0 value in molar units is greater than or equal to the maximum chlorine demand exerted by the sample.The procedure for solving Equations (31) and (36) simultaneously follows using the Euler forward method. The procedure is illustrated for the first two steps as follows. Let the initial assumed values of k at t = 0 X¯A/X¯H be k0 and X¯A/X¯H=a, respectively. In addition, let the initial reactant concentration be X0 as estimated from the TOC value of the sample, the actual molar value of X0 determined through chlorine decay experiments carried out on a series of dilution of the original sample, and the X0 value extrapolated to the undiluted sample. Using a time step Δt and the forward Euler method of solving first-order differential equations, the values at the next time step t1 = t0 + Δt are determined as follows.K
1 is determined using equation (31), i.e.,(38)k1−k0Δt=−C0k02a−1.NextC1 and X1 are determined from equation (36), i.e.,(39)C1−C0Δt=−kavC0X0=X1−X0Δt.The average decay ratekav in the above equation is determined from the following formula:(40)kav=k0+k12.The procedure is continued to the next time stept2 = t1 + Δt using the newly computed values of X1, C1 and k1 to compute the values X2, C2 and k2 corresponding to the time step t2. This procedure continues further until the reactants and/or the chlorine available is exhausted. As discussed above, the values of k0 and a=X¯A/X¯H are determined by calibration as the values that give the minimum sum squares of error between the experimentally determined residual chlorine concentrations and the calculated chlorine residuals according to the steps described above:


### 2.4. Expression for the Time Variation of the Overall Reaction RateK

The chlorine residual decay rate can be expressed as a first-order reaction with the overall reaction rateK, whereby(41)K=ktXT.The decay rate is given by the following expression (first order):(42)dCdt=dXTdt=−ktXTC=−KC.Using the time derivative, we get(43)dktdt=dK/XTdt=dK/dTXT−dXT/dtKXT2,dktdt=dK/XTdt=dK/dTXT−−KCKXT2.Rearranging further, we get(44)dKdt=dktdtXT−K2CXT.Using Equation (31) for dkt/dt, we get(45)dKdt=−Ckt2X¯AX¯H−1XT−K2CXT,=−CK2XT2X¯AX¯H−1XT−K2CXT,=−CK2XTX¯AX¯H+K2CXT−K2CXT.Eventually,(46)dKdt=−CK2XTX¯AX¯H.The overall rate contains similar factors as in Equation (31) in terms of the initial chlorine dose and the ratio of the arithmetic to harmonic mean of the molar concentration of the reactants. The additional factor is the total molar concentration of reactants XT. Higher molar concentration reduces the overall rate resulting in steeper curve approaching first-order reaction. Lower overall molar concentration XT increases the decay rate approaching second-order eventually resulting in flatter overall decay curve typical of second-order decay rate. Overall, the chlorine decay rate is steeper at low chlorine dose and higher molar concentration of reactants and is flatter on the contrary at higher initial chlorine dose and lower molar concentration of reactants. Since the overall rate K is molar averaged, these values can be determined initially similar to the initial chlorine dose.


### 2.5. Modeling of Variation of Overall Reaction Rate with respect to Initial Chlorine

Referring to the first-order rate of decay of chlorine formula provided by Equation (42), we get(47)dCdt=KC,where C is the concentration at t, t is the time usually converted to unit of days, and K is the overall reaction rate constant.From Equation (46), we obtain(48)dKdt=−CK2XTX¯AX¯H.Since the molar-averaged overall rate of reaction is determined by the initial molar concentrations of reactants which means considering only the partial derivative of the overall rate with the chlorine concentrationC,(49)∂K∂t=−αCK2,where(50)α=X¯AX¯H1XT=constant.Equation (49) can also be written as(51)∂1/K∂t=αC.Taking the first derivative with respect to the chlorine concentration, we get(52)∂∂C∂1/K∂t=α.Reversing the order of differentiation, we get(53)∂∂t∂1/K∂C=α.At any given timet, the above expression after integration becomes(54)∂1/K∂C=αt=constant=K0,or(55)1K2∂K∂C=−X¯AX¯HtXT=−K0,since(56)α=X¯AX¯H1XT=constant.Finally,(57)∂K∂C=−K0K2.Integrating between the initial rate atC0 = 0 and at any given initial concentration C0 gives(58)∫βKdKK2=∫0C0−K0dC0,where β is the initial reaction rate constant when the initial concentration of chlorine approaches zero.After integration, the expression becomes(59)−1K−1β=−K0C0.Finally,(60)1K=1β+K0C0.The regression-based modeling based on Equation (59) is carried out by linear regression of (1/K) against the initial concentration for a number of chlorine decay tests carried out at different initial concentrations of chlorine. The regression parameters β and K0 are determined from this step.The expression for the concentration-based reaction rate constant after the regression parameters has been determined which then becomes(61)K=β1+βK0C0.The overall decay rate modeling is then obtained by combining the traditional first -order decay rate with the concentration-based reaction rate constant:(62)dCdt=−KC,Ct=C0e−Kt.Substituting the expression for the concentration basedK value in the above equation, we get(63)Ct=C0e−β/1+βK0C0t.Equation (63) can be used to develop the bulk decay of chlorine in distribution systems. Water quality modeling programs such as the one used in EPANET can be adapted with the change in overall rate with respect to the initial chlorine dose taken into account with the use of equation instead of the constant values used in the program.


## 2.1. Experimental Determination of Chlorine Residual

For the experiment involving determination of chlorine residual at different times and under different initial chlorine dosages used, samples of water that underwent conventional treatment up to and including rapid sand filtration and just prior to disinfection were collected form the Matsapha water treatment plant and were subsequently used in the experimental trials. The samples collected were brought to the Chemistry Laboratory of the University of Eswatini. Different dosages of chlorine were added to multiple samples, and the chlorine residuals were determined at predetermined intervals of time. The method used for determination of residual chlorine was the iodometric titration as stated in Standard Methods for the Examination of Water and Wastewater [26]. The method is based on the principle that when potassium iodide is added to a sample of water containing residual chlorine at pH of 8 or less, the residual chlorine liberates iodine from the potassium iodide and the liberated iodine is titrated with sodium thiosulphate (Na2S2O3) titrant.In order to make the titration stoichiometric, the pH of the sample is reduced to between 3 and 4 by adding acetic acid. The method has expected a detection limit of 0.04 mg/L if 0.001 N sodium thiosulphate titrant is used. All chemicals used had reagent-grade quality, and fresh solutions were prepared every time the experiment was repeated at different time periods.Depending upon the concentration anticipated, suitable sample volume was taken so that the titrant sodium thiosulphate (Na2S2O3) volume would not exceed 20 mL and the starch indicator volume also is above 0.2 mL. To the sample volume, 5 mL acetic acid was added followed by 1 g of potassium iodide powder measured on weighing balance. After this, the sample was titrated with sodium thiosulphate (Na2S2O3) of appropriate normality. For samples with low anticipated chlorine concentrations, the titrant was prepared to be of low normality. The titration was continued until the yellow colour was almost disappearing. Then, 1 mL of the starch solution was added until the blue colour disappeared. In order to compensate for the method error, blank titration was also performed using distilled-deionised water and that passed through the same procedure as the one used for the actual sample. A minimum of three repetitions were performed for each determination. After titration, mg/L of chlorine residual present in the sample was determined using the following formula:(1)mgLof chlorine=A±B∗N∗35450mL of sample,where A is the volume of the titrant used for the sample, B is the volume of the titrant used for blank, and N is the normality of the Na2S2O3 titrant.


## 2.2. Basis of Mathematical Model

The time varying reaction rate constant equation as is suggested by Phillip et al. [5] is used for the mathematical formulation. According to this equation, the concentration-averaged rate of reaction variation with time is given by the following expression:(2)dktdt=Ctkt2−∑i=1Nki2XitXT.The original second-order variation of chlorine decay rate with molar concentration of chlorine and reactants and from which equation (2) was derived is given by(3)dCtdt=−ktCtXt=dXTdt.In Equations (2) and (3), kt is the molar weighed reaction rate, Ct is the chlorine residual measured at time t, Xit is the molar concentration of reactant i measured at time t and reacting with chlorine, XT is the total molarity of the reactants exerting chlorine demand, and N is the total number of chlorine consuming reactants present in the water.Meanwhile, the molar-averaged reaction ratekt is defined as(4)kt=∑i=1NkiXit∑i=1NXit=∑i=1NkiXitXT,where ki is the individual rate of reaction of reactant i having molar concentration xi and is exerting chlorine demand. Other symbols are as defined for Equation (2). Equation (2) indicates the second-order nature of the variation of the molar-averaged reaction rate constant kt with time. It is also clear from Equation (2) that other factors staying constant, the reaction rate constant kt varies as a first-order reaction with respect to the chlorine residual present. The aggregate effect of the presence of chlorine residual and this second-order variation of the rate constant are because as the chlorine concentration increases the rate constant tends to decrease. Similar observations were noted by several researchers through empirical experimental determination of the relationship between the reaction rate constant and the initial chlorine used in the experiments [13, 18, 22, 23]. The rate constant appears to decrease inversely with increase in the chlorine dose used in the experiments. Phillip et al. [5] also argue that the rate constant variation with time always stays negative in Equation (2) as the expression in bracket in this equation can be proven either negative or zero only. In other words, it cannot take positive values.Equation (2) as developed by Phillip et al. [5] is extended in this paper to provide a formula that relates the time variation of rate of reaction kt with chlorine residual present and the molar concentration of reactants. The difficulty of dealing with the indeterminate term in the bracket on the right-hand side of Equation (2) and involving the individual reactants can be resolved through mathematical averaging of the product of the individual reaction rates and the molar concentration of reactants. This averaging is later verified by theoretical formulation of the error term and Monte Carlo simulation of the error of averaging under different values of these product terms, i.e., molar concentration of reactants and individual reaction rates.Starting with the expression in bracket in Equation (2) and reformulating it such that(5)∑i=1Nki2Xit=∑i=1Nki2Xit2Xit=∑i=1NkiXit2Xit.Rearranging the right-hand side of the above expression, we get(6)∑i=1NkiXit2Xit=∑kiXit.kiXitXit.Now defining the average ofkiXit such that(7)kiXitav=∑kiXitN.Since it was already defined in Equation (4) that(8)kt=∑kiXitXT.Equations (4) and (7) are combined to give(9)kiXitav=∑kiXitN=ktXTN.Now returning to Equation (5) once again, we get(10)∑i=1NkiXit2Xit=∑kiXit·kiXitXit=kiXitav·∑kiXitXit.In the above expression, (kiXit/Xit) is used as a weighing factor. This is also analogous to the weighing of a continuous function, i.e.,(11)∫x.fxdx=x¯∫fxdx.In addition, by using similar averaging method, we get(12)∑kiXitXit=kiXitav·∑1Xit.In this case, 1/Xit are used as weighing factors.Overall,(13)∑i=1NkiXit2Xit=kiXitav2.∑1Xit=kt2XT2N2.∑1Xit.The use of weighted average in Equation (13) above for discrete values can be verified theoretically as well as practically using Monte Carlo simulation using random variation of individual reactants’ rate constants and reactants’ molar concentrations as described below. First, the theoretical basis is explained in the following:(14)∑i=1Nkx2≈kxav·∑i=1NX=1N∑i=1Nkx∑i=1NX.Working a little with the right-hand side expression, we get(15)1N∑i=1Nkx∑i=1NX=1N∑i=1Nkx2+∑i=1N∑j=1,j≠iNkixixj.Now defining the difference term in molar expression of reactants Δxij, we get(16)Δxij=xi−xj,so thatxj=xi−Δxij.Substituting the above expression for thexj term, we get(17)1N∑i=1Nkx∑i=1NX=1N∑i=1Nkx2+∑i=1N∑j=1,j≠iNkixixi−Δxij.The right-hand side of the above expression changes to(18)1N∑i=1Nkx∑i=1NX=1N∑i=1Nkx2+N−1∑i=1Nkx2−∑i=1N∑j=1,j≠iNkixiΔxij,1N∑i=1Nkx∑i=1NX=∑i=1Nkx2−1N∑i=1N∑j=1,j≠iNkixiΔxij.Because of symmetry, it is easy to show that(19)Δxij=−Δxji.Using this symmetry, the expression is reduced further to(20)1N∑i=1Nkx∑i=1NX=∑i=1Nkx2−1N∑i=1N∑j=i+1Nkixi−kjxjΔxij.Defining the difference terms Δkx such that(21)Δkxij=kixi−kjxj,and using this term in the above equation, we get(22)1N∑i=1Nkx∑i=1NX=∑i=1Nkx2−1N∑i=1N∑j=i+1NΔkxijΔxij.Therefore,(23)∑i=1Nkx2=1N∑i=1Nkx∑i=1NX+E,where the error term E is defined as(24)E=1N∑i=1N∑j=i+1NΔkxijΔxij.The above error term is summation over second-order difference (i.e., Δkx times Δx) and tends to be small compared with other terms in the equation containing E:The second practical verification is using Monte Carlo simulation of the error term and finding the average as well as bracketing the range of variation of the error termE. For this purpose, it is assumed that there were five chlorine demanding reactants with individual reaction rates ki varying by a minimum factor of 10 and maximum factor of 1000 and with molar concentrations of reactants xi also varying by a minimum factor of 10 and maximum factor of 1000. The Monte Carlo simulation is done for each range of variation 10000 times, with the following results shown in Table 1.Table 1
Results of Monte Carlo simulation: error of estimationE for different ranges of molar concentration of reactants and reaction rates.
K
10, X10
K
20, X20
K
50, X50
K
100, X100
K
200, X200
K
500, X500
K
1000, X1000
Average error,E
12.68
12.28
12.53
12.53
12.63
12.44
12.61
Standard deviation
8.68
8.67
8.68
8.76
8.79
8.73
8.84
Geometric mean of errors,E
8.98
NA∗
NA∗
NA∗
8.85
8.79
NA
Harmonic mean of errors,E
3.25
1.54
3.12
0.65
0.03
2.83
2.64
Maximum error,E
54.99
54.73
56.69
56.73
59.75
50.56
59.70
Minimum error,E
0.02
0.00∗
0.02
0.00∗
0.00∗
0.01
0.01
∗Geometric mean calculation invalidated by presence of zero in the data as is noticed in the last row of this table.As is shown in Table1, the mean error E is about 12% with 95% confidence range lying between 12.11 and 12.88 which is quite narrow range. The error term is therefore comparatively small, and the use of the weighing factor in the derivation of the rate equation is justified. The same can be said of the geometric and harmonic means of the error term E which are in the order of 8% and 3%, respectively, as shown in Table 1. These data support the heuristic use of the averaging used in deriving Equation (13) in terms of the ranges of variations of the kinetic rate and molar concentration of reactants.Now returning back to Equation (13) and using the arithmetic average (XA) and harmonic average (XH) of the molar concentration of the reactants exerting chlorine demand, we get(25)X¯A=XTN,X¯H=N∑i=1N1/Xit.The expression in Equation (13) can now be written as(26)∑i=1NkiXit2Xit=kt2XT2N2·∑1Xit=kt2XTXTN·1N·∑1Xit=kt2XTX¯ANNX¯H.Simplifying further, we get(27)∑i=1Nki2Xit=∑i=1NkiXit2Xit=kt2XTX¯AX¯H.Now, the original expression given in Equation (2) is(28)dktdt=Ctkt2−∑i=1Nki2XitXT,which can be written as(29)dktdt=Ctkt2−∑i=1Nki2XitXT=Ctkt2−kt2XTX¯A/X¯HXT,which is the same as(30)dktdt==Ctkt21−X¯AX¯H.Changing the order of the expression in bracket, we obtain(31)dktdt=−Ctkt2X¯AX¯H−1.The above expression in Equation (31) shows that the reaction rate is second order with respect to kt and first order with Ct. The reaction rate decreases with time being always negative since the expression in bracket is always positive (the arithmetic mean is always greater than the harmonic mean).Equation (31) as such replaces the individual reaction rates by the concentration-weighed aggregate reaction rate kt as it is the one that can be easily measured. The concentration-weighed aggregate reaction rate kt is influenced by the ratio of arithmetic to harmonic mean of reactants. If the molar concentration of reactants is close in magnitude to each other, then(32)X¯AX¯H≈1,X¯AX¯H−1⟶0.The ratekt becomes constant (dkt/dt = 0 in Equation (31)), and the chlorine decay model approaches first order rate reaction. By contrast, if there is significant variation in the molar concentration among the reactants, the ratio(33)X¯AX¯H−1≫0,holds true. The chlorine decay follows second-order reaction and is highly influenced by the reactant with greater molar concentration. Because of the validity of the expression(34)X¯AX¯H−1≥0,as the arithmetic mean is always greater than the harmonic mean, the time rate of decay of the concentration-weighed reaction rate is always negative and is decreasing with time. This is also observed in the literature through empirical inverse relationship between the chlorine decay rate and initial chlorine dose obtained experimentally. It is also clear from Equation (31) that the rate kt decreases faster with increase in initial chlorine dose. The implication is that the chlorine decay curve becomes flatter (second order) at higher initial chlorine dosages. In other words, within the distribution system, it takes longer time for the chlorine to decay compared to first-order rate decay and the residuals tend to be greater.


## 2.3. Procedure for Determination of the Chlorine Decay Rate Based on Equations (31) and (36)

In order to determine the chlorine decay rate, it is necessary to solve the differential equations involving equation (31) together with equation (36), namely,(35)dktdt=−Ctkt2X¯AX¯H−1,(36)dCtdt=−ktCtXt=dXTdt.To determine the decay rate from the above two equations, it is necessary to determine theXt and kt values by alternating between these two equations. In other words, the Xt and Ct values are determined from equation (36) and the kt value from equation (31). To begin with, however, the initial rates k0 and X0 as well as the ratio of arithmetic to harmonic mean of the reactants need to be determined. It is proposed that the initial k0 value as well as the expression in equation (31), i.e.,(37)X¯AX¯H−1,can be determined by calibration. This means assuming different values of k0 and X¯A/X¯H and choosing the values that give minimum sum of squares of error between the observed and calculated chlorine residual concentrations measured at different times. It is possible that the combination of k0 and X¯A/X¯H can be determined using the evolutionary algorithm or Monte Carlo method.It is also possible to relate theX¯A/X¯H value to reactant chemical species parameters involving, for example, TOC and UV absorption values. One possibility is using the TOC/UV ratio to relate with X¯A/X¯H and use this relationship to model the variable decay rate kt as given by equation (31).It is also proposed that the initial reactant concentrationX0 be the molar equivalent of the total organic carbon (TOC) concentration. The actual value of X0 can also be ascertained through dilution experiment by running a chlorine decay experiment at successive dilution and extrapolating the first order rate constant corresponding to the undiluted sample. The X0 value in molar units is greater than or equal to the maximum chlorine demand exerted by the sample.The procedure for solving Equations (31) and (36) simultaneously follows using the Euler forward method. The procedure is illustrated for the first two steps as follows. Let the initial assumed values of k at t = 0 X¯A/X¯H be k0 and X¯A/X¯H=a, respectively. In addition, let the initial reactant concentration be X0 as estimated from the TOC value of the sample, the actual molar value of X0 determined through chlorine decay experiments carried out on a series of dilution of the original sample, and the X0 value extrapolated to the undiluted sample. Using a time step Δt and the forward Euler method of solving first-order differential equations, the values at the next time step t1 = t0 + Δt are determined as follows.K
1 is determined using equation (31), i.e.,(38)k1−k0Δt=−C0k02a−1.NextC1 and X1 are determined from equation (36), i.e.,(39)C1−C0Δt=−kavC0X0=X1−X0Δt.The average decay ratekav in the above equation is determined from the following formula:(40)kav=k0+k12.The procedure is continued to the next time stept2 = t1 + Δt using the newly computed values of X1, C1 and k1 to compute the values X2, C2 and k2 corresponding to the time step t2. This procedure continues further until the reactants and/or the chlorine available is exhausted. As discussed above, the values of k0 and a=X¯A/X¯H are determined by calibration as the values that give the minimum sum squares of error between the experimentally determined residual chlorine concentrations and the calculated chlorine residuals according to the steps described above:


## 2.4. Expression for the Time Variation of the Overall Reaction RateK

The chlorine residual decay rate can be expressed as a first-order reaction with the overall reaction rateK, whereby(41)K=ktXT.The decay rate is given by the following expression (first order):(42)dCdt=dXTdt=−ktXTC=−KC.Using the time derivative, we get(43)dktdt=dK/XTdt=dK/dTXT−dXT/dtKXT2,dktdt=dK/XTdt=dK/dTXT−−KCKXT2.Rearranging further, we get(44)dKdt=dktdtXT−K2CXT.Using Equation (31) for dkt/dt, we get(45)dKdt=−Ckt2X¯AX¯H−1XT−K2CXT,=−CK2XT2X¯AX¯H−1XT−K2CXT,=−CK2XTX¯AX¯H+K2CXT−K2CXT.Eventually,(46)dKdt=−CK2XTX¯AX¯H.The overall rate contains similar factors as in Equation (31) in terms of the initial chlorine dose and the ratio of the arithmetic to harmonic mean of the molar concentration of the reactants. The additional factor is the total molar concentration of reactants XT. Higher molar concentration reduces the overall rate resulting in steeper curve approaching first-order reaction. Lower overall molar concentration XT increases the decay rate approaching second-order eventually resulting in flatter overall decay curve typical of second-order decay rate. Overall, the chlorine decay rate is steeper at low chlorine dose and higher molar concentration of reactants and is flatter on the contrary at higher initial chlorine dose and lower molar concentration of reactants. Since the overall rate K is molar averaged, these values can be determined initially similar to the initial chlorine dose.


## 2.5. Modeling of Variation of Overall Reaction Rate with respect to Initial Chlorine

Referring to the first-order rate of decay of chlorine formula provided by Equation (42), we get(47)dCdt=KC,where C is the concentration at t, t is the time usually converted to unit of days, and K is the overall reaction rate constant.From Equation (46), we obtain(48)dKdt=−CK2XTX¯AX¯H.Since the molar-averaged overall rate of reaction is determined by the initial molar concentrations of reactants which means considering only the partial derivative of the overall rate with the chlorine concentrationC,(49)∂K∂t=−αCK2,where(50)α=X¯AX¯H1XT=constant.Equation (49) can also be written as(51)∂1/K∂t=αC.Taking the first derivative with respect to the chlorine concentration, we get(52)∂∂C∂1/K∂t=α.Reversing the order of differentiation, we get(53)∂∂t∂1/K∂C=α.At any given timet, the above expression after integration becomes(54)∂1/K∂C=αt=constant=K0,or(55)1K2∂K∂C=−X¯AX¯HtXT=−K0,since(56)α=X¯AX¯H1XT=constant.Finally,(57)∂K∂C=−K0K2.Integrating between the initial rate atC0 = 0 and at any given initial concentration C0 gives(58)∫βKdKK2=∫0C0−K0dC0,where β is the initial reaction rate constant when the initial concentration of chlorine approaches zero.After integration, the expression becomes(59)−1K−1β=−K0C0.Finally,(60)1K=1β+K0C0.The regression-based modeling based on Equation (59) is carried out by linear regression of (1/K) against the initial concentration for a number of chlorine decay tests carried out at different initial concentrations of chlorine. The regression parameters β and K0 are determined from this step.The expression for the concentration-based reaction rate constant after the regression parameters has been determined which then becomes(61)K=β1+βK0C0.The overall decay rate modeling is then obtained by combining the traditional first -order decay rate with the concentration-based reaction rate constant:(62)dCdt=−KC,Ct=C0e−Kt.Substituting the expression for the concentration basedK value in the above equation, we get(63)Ct=C0e−β/1+βK0C0t.Equation (63) can be used to develop the bulk decay of chlorine in distribution systems. Water quality modeling programs such as the one used in EPANET can be adapted with the change in overall rate with respect to the initial chlorine dose taken into account with the use of equation instead of the constant values used in the program.


## 3. Results and Discussion

### 3.1. Modeling of Chlorine Decay Rate under Different Chlorine Dose Conditions

The bottle test experiments for the bulk decay of chlorine residual were run under different initial chlorine doses. The results of measurement of chlorine under different initial chlorine doses are plotted in Figures1–4. Within each batch experiment, the decay rate generally follows a first-order decay rate except for the very high initial chlorine (18.38 mg/L of initial chlorine dose) that shows deviation from linearity at the early stage of decay. The first-order model parameters were determined through linear regression of the exponential decay with coefficient of determination varying between 0.89 and 0.95. Figure 5 shows a summary of the experimental runs of the chlorine residual under different initial chlorine dosages used.Figure 1
Linearised plot of chlorine residual with the initial chlorine dose of 0.77 mg/L.Figure 2
Linearised plot of chlorine residual with the initial chlorine dose of 3.08 mg/L.Figure 3
Linearised plot of chlorine residual with the initial chlorine dose of 6.59 mg/L.Figure 4
Linearised plot of chlorine residual with the initial chlorine dose of 18.38 mg/L.Figure 5
Variation of residual chlorine under different initial chlorine doses used.


### 3.2. Variation of the Reaction RateK under Different Initial Chlorine Doses

The plot of the reaction rate constant obtained from each of the bulk decay experiment is plotted against the corresponding initial chlorine dose used and is shown in Figure6. It can be seen from the graph that there is an inverse relationship between the bulk decay reaction rate constant and the initial chlorine dose used in the experiment. This relationship is also similar to the second-order decay rate mathematical model that is developed and demonstrated in this paper and shown eventually in Equation (61), depicting the inverse relationship between the chlorine decay reaction rate constant K and the initial chlorine dose used.Figure 6
Plot of the bulk decay rate constantK against the initial chlorine dose used.The curve in Figure6 also shows that the reaction rate K reduces rapidly when the initial chlorine dosage is low. This a region where the chlorine dose is the rate limiting step or the chlorine decay reaction is chlorine limited. The steep slope of the curve in this region is an indication of the fact that the decay reaction is chlorine limited. This region is the first-order rate region because the reactant concentration does not change significantly as the decay reaction progresses being limited by the chlorine available. On the other hand, as the initial chlorine dose increases progressively, the decay rate curve tends to get flatter, demonstrating weak relationship between initial chlorine dose and reaction rate constant. In this case, the reaction is reactant limited. This region is second-order rate region since the reactant concentration changes significantly as the decay reaction progresses.The second-order decay rate model parameters, namely,β and K0 shown in Equation (61), were determined by regression after the expression in Equation (61) was linearised through the inverse plot, namely, 1/K versus the initial chlorine dose C. The resulting data are shown in Figure 7 plotted as a straight line with the coefficient of determination (R2) value of 0.99.Figure 7
Linearised plot of the reaction rate constant (1/K) against the initial chlorine dose used (C0).


### 3.3. Suggestion for Future Research

In order to be able to use the result of the mathematical model as stated in equation (31) and test its validity under different reaction conditions, it is helpful to run a series of experiments on different samples of water and determine the experimental constants. In doing so, it is further necessary to relate these experimental constants to reactant characteristics such as the total organic carbon concentration and UV absorption so that these parameters can help in setting the values of the experimental constants during mathematical modeling.The experimental results shown in this research were carried out on an identical sample of water collected from the Matsapha water treatment. As such, the influence of the reactivity of the reactant species on the model parameters of equation (63), namely, β and K0, cannot be fully ascertained unless the experiment is carried out over a wide range of samples. This is because the value of K0 depends on the reactivity of the most reactive reactant in the water while β depends on how quickly that and other reactants are used up by the chlorine. In order to be able to use equation (63) in water quality models such as the EPANET for wide range of water quality conditions, it is necessary to relate again the values of K0 and β with reactant chemical species parameter such as TOC and UV absorption so that appropriate values of these parameters can be selected for modeling purposes. To be able to do this, it is suggested that the initial chlorine based experiment is carried out over different types of water samples while at the same type the reactant species of the samples be characterised by measurement of appropriate reactant parameters such as TOC and UV absorption.


## 3.1. Modeling of Chlorine Decay Rate under Different Chlorine Dose Conditions

The bottle test experiments for the bulk decay of chlorine residual were run under different initial chlorine doses. The results of measurement of chlorine under different initial chlorine doses are plotted in Figures1–4. Within each batch experiment, the decay rate generally follows a first-order decay rate except for the very high initial chlorine (18.38 mg/L of initial chlorine dose) that shows deviation from linearity at the early stage of decay. The first-order model parameters were determined through linear regression of the exponential decay with coefficient of determination varying between 0.89 and 0.95. Figure 5 shows a summary of the experimental runs of the chlorine residual under different initial chlorine dosages used.Figure 1
Linearised plot of chlorine residual with the initial chlorine dose of 0.77 mg/L.Figure 2
Linearised plot of chlorine residual with the initial chlorine dose of 3.08 mg/L.Figure 3
Linearised plot of chlorine residual with the initial chlorine dose of 6.59 mg/L.Figure 4
Linearised plot of chlorine residual with the initial chlorine dose of 18.38 mg/L.Figure 5
Variation of residual chlorine under different initial chlorine doses used.


## 3.2. Variation of the Reaction RateK under Different Initial Chlorine Doses

The plot of the reaction rate constant obtained from each of the bulk decay experiment is plotted against the corresponding initial chlorine dose used and is shown in Figure6. It can be seen from the graph that there is an inverse relationship between the bulk decay reaction rate constant and the initial chlorine dose used in the experiment. This relationship is also similar to the second-order decay rate mathematical model that is developed and demonstrated in this paper and shown eventually in Equation (61), depicting the inverse relationship between the chlorine decay reaction rate constant K and the initial chlorine dose used.Figure 6
Plot of the bulk decay rate constantK against the initial chlorine dose used.The curve in Figure6 also shows that the reaction rate K reduces rapidly when the initial chlorine dosage is low. This a region where the chlorine dose is the rate limiting step or the chlorine decay reaction is chlorine limited. The steep slope of the curve in this region is an indication of the fact that the decay reaction is chlorine limited. This region is the first-order rate region because the reactant concentration does not change significantly as the decay reaction progresses being limited by the chlorine available. On the other hand, as the initial chlorine dose increases progressively, the decay rate curve tends to get flatter, demonstrating weak relationship between initial chlorine dose and reaction rate constant. In this case, the reaction is reactant limited. This region is second-order rate region since the reactant concentration changes significantly as the decay reaction progresses.The second-order decay rate model parameters, namely,β and K0 shown in Equation (61), were determined by regression after the expression in Equation (61) was linearised through the inverse plot, namely, 1/K versus the initial chlorine dose C. The resulting data are shown in Figure 7 plotted as a straight line with the coefficient of determination (R2) value of 0.99.Figure 7
Linearised plot of the reaction rate constant (1/K) against the initial chlorine dose used (C0).


## 3.3. Suggestion for Future Research

In order to be able to use the result of the mathematical model as stated in equation (31) and test its validity under different reaction conditions, it is helpful to run a series of experiments on different samples of water and determine the experimental constants. In doing so, it is further necessary to relate these experimental constants to reactant characteristics such as the total organic carbon concentration and UV absorption so that these parameters can help in setting the values of the experimental constants during mathematical modeling.The experimental results shown in this research were carried out on an identical sample of water collected from the Matsapha water treatment. As such, the influence of the reactivity of the reactant species on the model parameters of equation (63), namely, β and K0, cannot be fully ascertained unless the experiment is carried out over a wide range of samples. This is because the value of K0 depends on the reactivity of the most reactive reactant in the water while β depends on how quickly that and other reactants are used up by the chlorine. In order to be able to use equation (63) in water quality models such as the EPANET for wide range of water quality conditions, it is necessary to relate again the values of K0 and β with reactant chemical species parameter such as TOC and UV absorption so that appropriate values of these parameters can be selected for modeling purposes. To be able to do this, it is suggested that the initial chlorine based experiment is carried out over different types of water samples while at the same type the reactant species of the samples be characterised by measurement of appropriate reactant parameters such as TOC and UV absorption.


## 4. Conclusion

Prediction of chlorine decay rate is a valuable exercise for ascertaining chlorine residual that meets minimum requirements, ensuring safety of potable water supply coming from distribution systems. Traditional water quality modeling programs that use first-order decay rate are suggested to be modified to take into account the variable nature of the reaction rate constant for chlorine decay under different initial chlorine doses and with different molar concentrations of reactants each of which has variable rates of reaction with chlorine.Suggested models for handling variable chlorine dose and reactants participating in chlorine decay largely employ semiempirical approach, which limits their universal applicability. In addition, reliance on experimentally determined model parameters tends to make the models specific to the conditions under which the experiments were run. The mathematical model developed in this paper shows a greater move towards the analytical solution in line with the second-order decay rate of chlorine suggested by previous researchers. The mathematical models as shown in Equations (31), (46), and (61) have a more simplified analytical solution compared, for example, with the trial and error solution to a boundary value problem suggested by Phillip et al. [5]. Molar averaging of reaction rates partly verified by Monte Carlo simulation for average validity simplify the solution by eliminating the need for a trial and error assumption of reaction rates required, for example, in the solution proposed by Phillip et al. [5]. The solution becomes a simple solution to initial value problem of second-order decay rate variation of chlorine residual both for exploring the variation of decay rate with initial chlorine as well as molar concentration of reactants.The influence of the initial chlorine dose on the reaction rate constant for bulk decay rate of chlorine is inversely related as shown in Equation (61). This means that the reaction rate constant decreases when the initial chlorine increases. On the other hand, the influence of the different reactants present is explained by the ratio of the arithmetic mean to harmonic mean of reactants as shown in Equations (31) and (46). If the molar concentration of reactants is similar, the chlorine decay rate approaches first order and the influence of the molar concentration of reactants is minimal. By contrast, if there are reactants with low molar concentration, the reaction tends towards being second order with flatter decay rate curve similar to the effect of high initial chlorine dose as explained in Equations (31) and (46).The inverse relationship between initial chlorine dose used and the corresponding reaction rate constant for bulk chlorine decay, an experimentally observed reality quoted by several researchers, is shown in this paper as a mathematical solution of the second-order decay rate modeling when the only variable under consideration is the initial chlorine used (reactants and their molar concentration remain unchanged). The experimental trials for bulk decay of chlorine under the different chlorine doses used show that the inverse relationship and the regression model acceptably fit the predicted second-order rate variation.



---

*Source: 5863905-2019-05-23.xml*