Pollutant Dispersion Modeling in Natural Streams Using the Transmission Line Matrix Method

Numerical modeling has become an indispensable tool for solving various physical problems. In this context, we present a model of pollutant dispersion in natural streams for the far field case where dispersion is considered longitudinal and one-dimensional in the flow direction. The Transmission Line Matrix (TLM), which has earned a reputation as powerful and efficient numerical method, is used. The presented one-dimensional TLM model requires a minimum input data and provides a significant gain in computing time. To validate our model, the results are compared with observations and experimental data from the river Severn (UK). The results show a good agreement with experimental data. The model can be used to predict the spatiotemporal evolution of a pollutant in natural streams for effective and rapid decision-making in a case of emergency, such as accidental discharges in a stream with a dynamic similar to that of the river Severn (UK).


Introduction
The evolution of a pollutant in both natural and artificial streams may be subject to different phenomena such as dispersion, diffusion, advection, sedimentation, adsorption, desorption, etc. Predicting its development and its spread is important for the environmental protection.In the field of water quality modeling, several researchers [1][2][3][4][5] presented different approaches to understand and interpret the basic concept of water quality problems.In a case of emergency, such as accidental discharges in a stream, the prediction of the pollutant evolution is crucial in effective and rapid decision-making.Such a help in decision-making must be simple and precise, needs a reduced computational time and a minimum input data.Taking into account most dispersion phenomena found in streams complicates the model by increasing the number of input data which makes simulations more difficult, slow, and inadequate for emergency decision-making, where simplicity and speed are the two critical parameters [6].
Dispersion is a combination of diffusion, advection, and shear.It is created by the non-uniformity of velocity fields related to the different characteristics of the stream such as geometry, roughness, and kinematics.The dispersion area is usually composed of three distinct zones: The initial mixing zone, the full mixing zone and the far field zone [6].In this work, we focus on the third zone where dispersion is considered longitudinal and one-dimensional in the flow direction.The one-dimensional differential Advection-Dispersion Equation (ADE) describing the phenomenon is given by Equation (1) [6][7][8][9][10][11], where the following assumptions are considered: • Vertical and transversal dispersions are very small; • The pollutant is completely miscible in water; • Chemical reactions between the pollutant and its environment are absent; and • The overall mass of pollutant is maintained during transport.
where: c is the concentration of pollutant (g/l), U is the flow velocity (m/s), DL is the longitudinal dispersion coefficient (m 2 /s), x is the distance (m), and t is the time (s).Since velocity and dispersion coefficient are constant, Equation (1) becomes Equation (2) or Equation (3).

Methods
The setup of the modeling method is made of the following steps: (1) get the best description of the physical phenomenon (pollutant dispersion in stream); (2) define the characteristic parameters of the profile temporal evolution of the pollutant dispersion; (3) find the experimental parameters of the pollutant dispersion; (4) set the hydraulic conditions corresponding to the phenomenon; (5) develop the TLM (Transmission Line Matrix) method and the TLM model to simulate the phenomenon; (6) optimize the TLM parameters by comparing the obtained TLM results to the experimental data; (7) validate statistically and interpret physically the differences between the TLM results and the experimental data; and (8) deduce the final TLM dispersion model.

Physical Phenomenon Description
To study the pollutant dispersion phenomenon in a channel, the experimental data of Atkinson and Davis [12] have been chosen.The experiment was conducted over a length of approximately 14 km on the river Severn in Wales (UK).A volume of 5 L of tracer (Rhodamine WT 40 g/l) is poured for 107 s while moving on Llanidloes Bridge.At this level, the water depth is equal to 0.30 m while the average width is about 24 m.Samples were taken along the test distance and analyzed by molecular fluorescence to find the temporal evolution of the pollutant concentration at seven locations (stations) along that distance.Table 1 shows the characteristic parameters of the temporal evolution that are the start time of pollution (t0), the maximum concentration (Cmax), the peak time (tp), and the end time of pollution (tf).The most freely-meandering part of the test channel is downstream of the station E where the flood plain broadens.The channel width increases in the downstream direction, averaging 24 m.Upstream of the station E, the channel is narrower, with an average width of 21.4 m (standard deviation: 4.7 m), whereas downstream it is wider with a greater variability of 28.3 m (standard deviation: 7.5 m).The channel depth varies considerably because the long profile of the riverbed consists of a series of riffles and pools (topography), the average depth is 0.5 m (standard deviation: 0.18 m).Between 6.5 km from the injection point and station E, the channel is mostly much deeper.There is a relatively sudden increase in average section at about 6.5 km downstream of the injection point, upstream of this distance, the average section is 10.7 m 2 (standard deviation: 2.8 m 2 ) while below this, it is 13.8 m 2 (standard deviation: 4.1 m 2 ) [12].

Moment Method and Experimental Parameters Determination
The method of moments (Equation ( 4)) [13][14][15] is used to determine, experimentally, the longitudinal dispersion coefficient DL: where, σ (s ) is the temporal variance given by Equation (5) [16]; ̅ (s) is the time to the centroid of the concentration distribution, determined by Equation ( 6) [15] and U(m/s) is the mean flow speed upstream of each station i, calculated by Equation ( 7) [13,15,17].In analyzing tracer dispersion data, some models may require upstream average values for velocity which is, indeed, more real [12].
where, xi is the distance between station i and dumping point (m).
In the case study, at each station, the concentration is a discrete series of N values sampled at a certain time interval, so the equivalent expressions to Equation ( 5) and Equation ( 6) are respectively: Thus, DL and U are essentially the two experimental parameters of the pollutant dispersion.

General Assumptions and Hydraulic Conditions
This study assumed that the tracer is conservative and miscible in water without chemical reactions with the medium and that the test channel is divided into seven reaches, each one of them is regular, uniform, and is characterized by a stable flow with a constant velocity and a constant dispersion coefficient.It focuses on the channel far field where dispersion is considered longitudinal and one-dimensional in the flow direction.The effect of dead zones, skewness, roughness, expansions, and contractions are neglected.
To determine the far field area where longitudinal dispersion is predominant, Fischer proposed Equation (10) to calculate the length of the mixture, which was found to be 880 m [6,18].Thus, station A is not considered since its distance from dumping point is below this value.

=
where, Lm is the length of the mixture (m), K2 is the constant depends on the manner of tracer pouring, U is the flow speed (m/s), B is the average width between the dumping and the sampling points (m), and EY is the lateral dispersion coefficient (m 2 /s).

Modeling Method
The study of propagation, diffusion, and dispersion phenomena by analytical methods is possible only for few cases.Therefore, several researchers have proposed new approaches to better represent these phenomena.For several decades, a great number of numerical methods were initiated but found real use only in recent years with the rapid evolution of computing.Thus, numerical modeling has become an indispensable tool for a wide range of applications.Among these methods, Transmission Line Matrix (TLM) method is being increasingly used since seventies [11], and its application field continues to expand [19][20][21].

Transmission Line Matrix Method
TLM is a spatiotemporal numerical technique, explicit and stable, based on electrical networks [19][20][21][22].It is introduced by Professors Peter Johns and Raymond Beurle [23] at the electrical engineering department of Nottingham University (UK).This technique uses Huygens' principle and is based on Maxwell's equations [11,22,24,25].It operates on a mesh structure where each element is represented by a transmission line that acts as an analogy between the physical quantity and an electrical equivalent (voltage or current) [20,[24][25][26].
In TLM, each transmission line is excited with a pulse which will be monitored throughout the network.The injected pulse propagates along the line until it reaches a discontinuity (node) where it disperses [27].Each reflected pulse moves back along the line and becomes incident on the adjacent node after a time step Δt(s), and so on [19,21] (Figure 1).At each iteration k (Equation ( 11)), the TLM routine calculates the incoming pulses to determine the evolution of the physical quantity at a given point, then the reflected pulses in preparation for the next iteration [26,28] TLM can be one-dimensional, two-dimensional, or three-dimensional [11,28].
TLM has been compared with various other numerical methods such as the Finite Difference method (FDTD) [22,26,29].These comparisons define the power of this method to the studied case.

TLM Model for Pollutant Dispersion Phenomenon
Early TLM models for dispersion of a physical quantity, proposed the study of both phenomena (advection and diffusion) separately [30,31].They were based on the concept that transfers the physical content of each node to the node adjacent downstream along the flow direction.According to the results obtained by de Cogan and Henini [11,31], this approach is poor.An alternative method is proposed [11,32] where the TLM model of dispersion is identical to the TLM model for diffusion with an added current generator Im(A) at each node.Hence, a TLM model of a pollutant dispersion phenomenon is a perfectly insulated transmission line (conductance G(Ω) null) of length Δx(m), that is electrically represented by two resistors R(Ω), two inductors L(H), and a capacitor C(F) (Figure 2).The TLM node is given on Figure 3, where I(x), V(x), and Z are, respectively, the intensity of current (A), the voltage (Volt), and the characteristic impedance (Ω) of the line.The latter is defined by the relation from Equation (12).Rd, Cd, and Ld are, respectively, distributed resistance, capacitance, and inductance over the entire length Δx of the transmission line and are given by the following relationships: The variation of the voltage on a piece of transmission line is equal to the sum of the voltage variation due to the inductance and the voltage drop across the resistor as shown in Equation ( 16): The current intensity variation flowing through this piece of transmission line is equal to the current flowing through the capacitor to which we add a current generator Im at the node (n) as indicated by Equation ( 17): where, : Transconductance (Ω).If ∆x→0, Equation (17) becomes: If we neglect the influence of the inductance by taking a small enough time step, the combination of Equation ( 16) and Equation (19) gives: This equation is similar to Equation 3 and, therefore, we can model a phenomenon of pollutant dispersion by an electrical equivalent.
The addition of a current generator produces an additional voltage at node (n).Hence, the total pulse at iteration k, (n) becomes: ( ) and ( ) are, respectively, the incident pulses on left and right of node (n) at iteration k. ( + 1) and ( − 1) are the pulses at nodes adjacent (n+1) and (n-1) at iteration k-1.
Pulses reflected on left and right of node (n) at iteration k: ( ) and ( ) are given by Equation (23) and Equation ( 24): The reflected pulses become incident on adjacent nodes at next iteration.( ) and ( ) are incident pulses on left and right of node (n) at iteration k+1 and are given by Equation (25) and Equation ( 26): ( − 1) and ( + 1) are respectively reflected pulses on right of node (n-1) and left of node (n+1).Equations ( 21) to (26) form the TLM algorithm used to solve a wide variety of dispersion problems, such as chromatographic processes, heat exchanges, heat flow and diffusion of electric species in semiconductors [32,33].Analogy between Equation (3) and Equation (20) gives: ( , ) ⁄ = − ( , ) According to the Courant-Friedrichs-Lewy stability criterion [11,34], TLM Dispersion models are stable if: 1 ≥ 2 ≥ , where s is the convection number (Equation ( 31)) and r is the diffusion number (Equation ( 32)).

TLM Initial and Boundary Conditions
In the TLM method, the conditions cited in paragraph (2.1.3)are equivalent to: The first node (n = 1) is the tracer dumping point (x = 0) with the boundary condition: At the last node (n), the boundary condition is: The incident pulse returns with the same magnitude but opposite sign.In TLM, this boundary corresponds to a short circuit where the reflection coefficient Γ is equal to −1 [11,20,21].

TLM parameters Optimization
The TLM dispersion model is applied on experimental data of Atkinson and Davis [12] and uses the input data of Table 2.The injected total tension is equivalent to initial concentration at dumping point (Equation ( 29)).The tracer is poured into an elementary volume of the river ∆v = ∆x∆y∆z, where ∆y and ∆z are, respectively, the width (m) and depth (m) of the river at the dumping point.This dilution leads to an initial concentration of the tracer of 1110 μg/L at dumping point.The TLM parameters Rd, gm, and initial total tension V0 have been optimized to keep the modeling simple.In TLM technique, the time step ∆t must be less than the time constant (RC) leading to a choice of Rd > 0.8.The conductance is the resistance inverse, so gm < 1.25.The initial total tension at each node must be inferior to the injected total tension (first node), so V0 (n ≠ 1) < 1110.
The TLM model, with the input data of Table 2 and the above TLM parameters, is run.The obtained TLM results V = f(k) are compared to the experimental data of Atkinson and Davis C = f(t), i.e., the TLM characteristic parameters (k0, Vmax, kp, kf) are compared to the electrical equivalents of experimental characteristic parameters (t0, Cmax, tp, tf) (Table 1).The electrical equivalents (k'0, V'max, k'p, k'f) are deduced from Equation (11) and Equation (29).The optimized TLM parameters (Rd, gm, and V0) correspond to the minimal values (∆k0, ∆Vmax, ∆kp, ∆kf).
This optimization is acceptable if we can validate, statistically, and interpret, physically, the differences between the experimental parameters and the physical equivalents of TLM parameters which are calculated using Equation (27), Equation (28), and Equation (29).

Statistical Tests for TLM Model Validation
The model performance is analyzed by comparing, for each station, observed data (experimental data of Atkinson and Davis (Table 3)) with TLM model predictions (physical equivalents of optimal TLM results (Table 4)) for the following parameters: flow velocity (Utlm and U), longitudinal dispersion coefficient (Dtlm and DL), maximum concentration (Cmaxtlm and Cmax), peak time (tptlm and tp), start time of pollution (t0tlm and t0), and end time of pollution (tftlm and tf).
The statistical indices are the basic tools for analyzing the model quality [6,14,[35][36][37].The following tests have been chosen: The coefficient of divergence Rdiv with Rdiv >1 (Rdiv < 1) indicates an over-prediction (under-prediction) of a value; the mean percentage error E(%) representing a convenient model when it is small; the Mean Relative Square Error MRSE indicating a successful model if it is close to 0; the Factor Of EXcedence FOEX(%), signifying an overestimation (underestimation and/or correct estimation) of all values when it is equal to 100% (0%); the factor of two FA2(%) representing a reliable model if it is close to 100% and finally the scatter diagram indicating a situation of overestimation (underestimation) when the point is above (below) "y = x" line.

Determination of Experimental Parameters
The experimental concentration data of Atkinson and Davis [12] have been considered to calculate at each station: ̅ by Equation ( 9), σ by Equation ( 8), and U by Equation (7).DL is determined using Equation ( 4).The results are presented in Table 3 The flow velocity varies from one station to another, and so does the longitudinal dispersion coefficient with distance x.These two parameters depend on the roughness and the variable geometry of the river.The value of longitudinal dispersion coefficient of station E is very high because of an obvious change in geometry of the river.The value of longitudinal dispersion coefficient of station G is negative.According to Atkinson and Davis [12], the time was insufficient to carry out all the necessary measurements at station G; thus, the dispersion monitoring was incomplete.The following day, measurements were remade only at two stations: station F, to compare the results of the two experiments and describe the dispersive behavior of the test channel, and station G, to complete the concentration recording according to time.The comparison of results showed that the dispersive behavior of the test channel is reproducible.Finally, all the results were compiled by taking measurements of the first experiment at the stations A to F and adjusted measurements of the second experiment at the station G. Since the longitudinal dispersion coefficient value depends essentially on the spatio-temporal monitoring of concentrations, the lack of reliable experimental data at station G led to obtain a negative coefficient of longitudinal dispersion.Thus, station G will be ignored from our comparative study.

Determination of TLM Parameters
Figure 4 shows the optimized TLM results and Table 4 their physical equivalents.Globally, there is a good agreement between the TLM results and experimental data.For stations E and F, the TLM curve is slightly shifted from the experimental curve; this means that the TLM model underestimates the beginning and end of pollution and gives a polluting trail shorter than that observed with an overestimation of flow velocity.This is certainly due to the sudden and obvious change of the channel geometry at 6.5 km from the injection point (upstream of the station E); downstream from this distance, the channel section is bigger, and downstream from the station E, the channel is wider and more sinusoidal, but the TLM model does not take into account the channel geometry which strongly influences the dispersion evolution [10,17,38,39].Furthermore, between 6.5 km from the injection point and the station E, the channel is much deeper because its bed is composed of a series of riffles and pools storage areas or dead zones [9], which retain pollutant then release it after some time, with such effect stretching, longitudinally, the pollution slick leading to the slight difference between the experimental and modeling curves; but the TLM model does not take into consideration the channel topography which strongly influences the dispersive behavior of the channel [7,9,16,40].It's important to remark that Utlm, Dtlm, and C0tlm are, respectively, the flow velocity, the longitudinal dispersion coefficient, and the initial concentration predicted by the TLM model.
Our model does not take into account the variable geometry and topography of the channel in order to reduce the model input data, since a model is just required to be able to provide useful predictions in case of emergencies.Actually, the study of complex geomorphologic flows requires specific computational methods [10,41] that take into account: Geometry, curvature, meandering, flooding conditions, bed topography, riparian vegetation, and secondary flows that influence the meandering evolution [42].

TLM Model Performance
The statistical tests results are given in Table 5 and shown in Figure 5. Figure 5a presents the comparison between the predicted and experimental values of flow velocity.According to Figure 5a and Table 5, the scatter diagram and FOEX indicate that all the experimental values are over-predicted with a mean percentage error equal to 11.6%.The model seems to be able to simulate the observed values of flow velocity well at stations C and E, with a weak over-estimation at stations B and F, and with a significant over-prediction at station D. This may indicate that at station D, the predicted tracer cloud disperses more quickly than the observed one.The experimental flow velocity at station D is the smallest (Table 3) so the passage of the tracer could is slow.It can be slowed down by the different river morphologies: sinuosity, topography, dead zones, expansions and contractions; these aspects are known to affect the distribution profiles of velocity.Thus, the difference between TLM and experimental values can be attributed to the fact that the 1D-TLM model does not take into account these aspects.The values of Fa2, MRSE, and E indicate that the model representation of the experimental flow velocity is acceptable.On Figure 5b, the predicted longitudinal dispersion coefficients are plotted against the observed ones.According to Figure 5b and Table 5, the scatter diagram and FOEX indicate that there are situations of over and under-estimation with a mean percentage error of 40.78%.The model gives a good estimate of the longitudinal dispersion coefficient for station B, but slightly over-predicts the coefficients of stations C and F, over-predicts the coefficient for station D, and under-predicts that of station E. The differences between the predicted and observed values are remarkable for the two stations D and E. The comparison of these two stations to the others shows that station D has the smallest coefficient, while station E is characterized by the highest coefficient (Table 3).Looking at station D, its smallest flow velocity has led to a smaller coefficient since DL decreases with U. The overestimation of the flow velocity explains the overestimation of the longitudinal dispersion coefficient at this station.For station E, its coefficient is the highest because of an obvious and sudden change in river geometry upstream of the station E at 6.5 km from the injection point.In this region, the channel is much deeper and starts to widen; station E has the biggest section.The longitudinal dispersion coefficient depends on the width and depth of the channel; thus, it is strongly related to the channel geometry.Therefore, the difference between the predicted and observed longitudinal dispersion coefficients at station E is due to our model not taking into account the channel geometry.The values of Fa2, MRSE, and E indicate that the model representation of the experimental longitudinal dispersion coefficients is fairly acceptable.Moreover, the numerical estimation of this coefficient is very difficult; in the literature [7][8][9]14,17,38,[43][44][45][46] for natural rivers, there is a great disparity between the values of longitudinal dispersion coefficients calculated by different methods and empirical formulas.This shows that such a coefficient is difficult to estimate, its experimental determination is the fairest approach because it is based on spatio-temporal monitoring of pollutant concentrations in streams with the real effect of skewness, topography, and water storage areas.
Figure 5c shows the predicted maximum concentrations compared to the experimental ones and Figure 5d presents the comparison between predicted and experimental values of the maximum concentration time (peak time).The scatter diagram of Figure 5c and the statistical indices of Table 5 show that the TLM model correctly predicted all maximum concentrations with a null mean percentage error.The scatter diagram on Figure 5d shows that the maximum concentration time predicted by TLM model is in good agreement with that obtained experimentally.E is equal to 0.8%, MRSE is zero, FOEX is null, indicating that there is no overestimation, and Fa2 factor is equal to 100%.All the statistical tests indicate that the model provides a very well representation for the experimental maximum concentrations and their peak times.The principal role of a pollutant model is to correctly give the pollution degree (maximum concentration and its peak time) at a given distance; the hereby-presented original model satisfies this requirement.
Figure 5e,f present the comparison between predicted and experimental time values of the pollution start and pollution end, respectively.The scatter diagrams show that the model slightly underestimates the pollution start time and considerably underestimates the pollution end time.This may indicate that the predicted cloud passage across each sampling station is earlier than the observed cloud passage and the model predicts a trail of pollutant shorter than that observed experimentally.This can be explained by the fact that the long profile of the riverbed consists of a series of storage areas (dead zones and pools) which retained the tracer and released it after a certain duration; this retarded the cloud passage at each station and longitudinally stretched the pollutant slick.Thus, the difference between the TLM and the experimental values is due to the fact that the TLM model does not take into account the phenomena of roughness and dead zones.For the pollution start time, E is equal to 13.8%, MRSE is 0.04, FOEX is null, confirming that all the experimental values are under-predicted, and Fa2 factor is equal to 100%.For the pollution end time, E is equal to 36%, MRSE is 0.32, FOEX is null, and Fa2 factor is equal to 60%.These statistical tests show that the model provides a good estimate of pollution start time but less of an estimate of the pollution end time.It is preferable for a pollution model to predict early tracer cloud passage so one can react as quickly as possible to preserve the environment.The hereby-presented original model fulfills this demand.

TLM Model Stability
In the considered case study, the condition of the Courant-Friedrichs-Lewy stability criterion [29] is respected for every station (Table 6).

Final TLM Model
Differences between TLM parameters and experimental parameters are statistically validated and physically interpreted as due to neglecting some minor phenomena characterizing pollutant dispersion in a river.Therefore, we defined α, β, and γ as factors representing ones of these neglected phenomena with α = DL/Dtlm, β = U/Utlm, and γ = C0/C0tlm.These factors (Table 7) are different because of the channel geometry and topology variation from a station to another.
This new model can present the pollution evolution at each station by knowing only the channel length, poured pollutant quantity, longitudinal dispersion coefficient, and flow velocity.

Conclusions
The proposed one-dimensional TLM model provides a good estimate of pollutant longitudinal dispersion in the far field zone of a river, such as the river Severn (UK).Usually, this model overestimates flow velocities and longitudinal dispersion coefficients but describes quite well the evolution of maximum concentrations over time.It underestimates time of start and end of pollution.This overestimation or underestimation may be related to the natural aspects of the river such as roughness, topography, geometry, and dead zones that longitudinally stretch the pollutant slick.Statistical tests indicate that the model is quite efficient and reliable.Generally, a model representing a phenomenon of pollution should preferably underestimate pollution start time and peak time, and not the other way, so to stress the situation and, hence, respond as quickly as possible.According to the Courant-Friedrichs-Lewy stability criterion, our TLM Dispersion model is stable.Therefore, the hereby-presented original model with few input data can be used to predict the spatiotemporal evolution of a pollutant slick in natural streams with a dynamic similar to that of the river Severn (that is, the first considered case study).

Figure 1 .
Figure 1.Dispersion of an injected pulse.(a) pulse reflection; (b) the first iteration result; and (c) the second iteration result.

Figure 4 .
Figure 4. TLM modeling and experimental data of Atkinson and Davis.The total tension (tracer concentration) is traced according to the iteration number k at stations B, C, D, E, and F.

Table 1 .
Characteristic parameters of pollutant concentration temporal evolution.
Note: x: Distance from dumping point.

Table 3 .
: Experimental values of flow velocities and longitudinal dispersion coefficients.

Table 4 .
Optimized TLM parameters and their physical equivalents.

Table 5 .
Statistical tests for flow velocity, longitudinal dispersion coefficient and maximum concentration.

Table 6 .
Convection and diffusion numbers.