Shallow Water 1D Model for Pollution River Study

In this paper a finite element 1D model for shallow water flows with distribution of chemical substances is presented. The deterministic model, based on unsteady flow and convection-diffusion-decay of the pollutants, allows for evaluating in any point of the space-time domain the concentration values of the chemical compounds. The numerical approach followed is computationally cost-effectiveness respect both the stability and the accuracy, and by means of it is possible to foresee the evolution of the concentrations.


Introduction
When considering river flows with an essential lack of transversal velocity component, the obvious choice to calculate the discharge and the pollutant concentrations is to use a 1D shallow water (SW) model. One dimensional schemes are still a practical tool for modeling river flows, in particular for simulations over extended time [23]. For flow essentially 2D, of course, a fully 2D shallow water model should be considered [2]. In literature, exists a large amount of numerical models used for the solution of the systems of partial differential equations, they are traditionally based on finite element, finite difference or finite volume approximations in space and on implicit or explicit schemes in time [1]. We adopt for the solution of the SW system a fractional step scheme in time and a finite element (FE) approximation in space [2,15]. In particular, the choice of the fractional step afford, at every time step, to decouple the physical contributions so that waves traveling at speed of √gh can be calculated implicitly, with a low computational cost. In order to save computational time, moreover, the non-linear part of the momentum equation (and also the total derivative of the chemical pollutant transport equations) is solved using the characteristics method, so that the most onerous computational kernels are reduced to elliptic like ones. In particular, for the solution of the algebraic systems we use a direct Thomas method for the tri-diagonal matrix of elevation and an iterative method, preconditioned with an additive scalable Schwarz preconditioner, for the matrices of the concentrations of pollutants. By the numerical model it is possible, assigned the suitable initial and boundary conditions, to calculate the unit width discharge and the concentration of pollutants in each point of the space-time domain at a very low computational cost, running on a common personal computer the informatics code in which the model is implemented. Several approaches have been developed for predicting the evolution of chemical concentrations of pollutants along rivers, in particular relatively of plant protection products due to their potential adverse effects on human health and environment [3,7,9,17,22]. In this paper, we present a study relevant to the distribution of three chemical compounds: Pirimicarb (P), Dithianon (D), Methoxyfenozide (M), along 3500 m of the Novella River located in the north-east part of Italy, during a transient of 28 hours occurred in the autumn of 2009. The distributions are characterized by transport and decay phenomena and for the study has been adopted the new numeric model based on the shallow water equations for the flow and three convection-diffusion equations and decay lows (without air and sediments interactions) for the pollutant concentrations.
The remainder part of the paper is organized as follows: in the next section the numerical model is presented in detail; then follows two sections reporting two numerical tests with known analytical solution; in the fourth section an application to a real problem is presented (the Novella river problem) reporting the results of a transient of 28 hours occurred in the autumn of 2009; finally some conclusions are drawn in the last section.

Governing Equations
To develop the numerical model, we start from the 1D classical shallow water equations represented by the Saint-Venant system, it is formed by the momentum equation (with dependent variable the unit width discharge q) and by the water height h equation (in our formulation replaced by the elevation ξ of the free surface respect to a reference level). The 1D approach presented here is suitable for situations in which the river presents an almost constant width and meanders are absent, it is based on the assumptions of constant flow in the vertical direction and hydro-static pressure. The Saint-Venant system of partial differential equations, in conservative form reads: and the transport and decay equations for concentrations read: = " − # , = , , !
where: x is the spatial coordinate, t is the time, h(x,t) is the total depth water measured from the bottom of the river, is the elevation over the reference plane, f(x,t) is the source term, g is the gravitational acceleration, µ is the diffusion coefficient, α is the bottom slope of the river. Marking by i in turn P, D, M, a i (x,t) is the reaction function, c i (x,t) is the concentration of the pollutant, s i (x,t) is the source term for transport equation, Γ i is the diffusion coefficient, S i (x,t) is the source term for decay equation and r i is a decay ratio. We recall that (1) represents the conservation of the momentum, while (2) ensure the conservation of mass, (3) is the convection dominated transport of the pollutants along the river and, finally, (4) keeps into account the decay effects. Suitable initial and boundary conditions have to be added in order to make the overall system solvable.

Numerical Scheme
Some features of the numerical approach developed have to be remarked, in fact the non-linear convective term of momentum equation is approximated by the characteristics method, so that the spurious oscillations due to centered approximations are avoided under mild restrictions on the time step and, moreover, a fractional step procedure is adopted for advancing in time. This last choice is made in order to decouple, at every time step, the physical contributions [5] obtaining the important advantage that the critical time step is in term of the flow velocity instead of the wave celerity. The spatial discretization of the equations is based on a Galerkin FE method, with two different spaces for the unknowns: polynomials of degree two for the unit width discharge and concentrations, polynomials of degree one for the elevation. It is notable to underline that the following advantages are guaranteed from our FE choice: the possibility to use h-adaptive techniques (used in this paper for the solution of the second numerical test), the absence of numerical viscosity inherent to the low order interpolation at the foot of characteristics and the elimination of the spurious oscillations that can arise when an equal order approximation is used [26]. Let Ω a domain in R and [0, T] the time interval of study. In detail, the fractional step we use for advancing in time is: and (x f is the foot of characteristic relevant to the node with abscissa x) 1 $% 0 2 = 1 $% 3 2 − ∆, ℎ $ + ∆, 5 0 1 $% 3 where r i = r i ∆t and a i ,, s i ,, S i are null. Remark 2.2.1. Equations (15) and (16) make completed the phenomenon of convection, diffusion and decay.
Remark 2.2.2. The use of characteristics for the approximation of the convective terms (both for the momentum equation and for the transport equations of chemical substances) requires an interpolation of high order; this is easily obtained because of our P2 choice for q and c i variables (we recall that we use P1 elements only for ξ variable).

Algebraic Formulation
Since the spatial approximation is based on Galerkin FE method, the algebraic formulation reads: for the provisional discharge q n+2/3 , indicating by M L q the lumped mass matrix, by A q the stiffness matrix, by For the elevation ξ n+1 , indicating by M L and by A the lumped mass matrix and stiffness matrix respectively, by B T the matrix of the integral of the product between the basis function and its derivative, by G 2 the diagonal matrix with g 2i,i =g(h n ) i , by G 3 the diagonal matrix with g 3i,i = (q n+2/3 ) i /(h n ) i : In equation (18), the q values considered are those associated with the boundary nodes of the elements.
For the update value of the q n+1 discharge, indicating by M L q the lumped mass matrix of discharge, by C T the matrix of the integral of the product between the basis function and its derivative, by G 2 and G 3 the above defined diagonal matrices: In (19), the ξ values considered are those calculated in the boundary nodes and in the middle node of each element.
For the chemical pollutant c i , indicating by M L ci the lumped mass matrix and by A ci the stiffness matrix:  (18) is tridiagonal, the Thomas algorithm has been used for its solution, while the Schwarz preconditioned iterative Bi-CGSTAB solver has been used for the systems of (20). In the next section a short presentation of the Schwarz preconditioner is presented.

Schwarz Scalable Additive Preconditioner
We give a brief description of the Schwarz preconditioner for elliptic problems introduced in [21] and adopted in this paper. In order to fix the ideas, we refer to a Poisson problem with homogeneous boundary conditions: The FE approximation to the Poisson problem consists in If we consider, for example, the domain as decomposed in two overlapping sub-domains Ω 1 and Ω 2 , let be ∂Ω, ∂Ω 1 , ∂Ω 2 the boundaries of Ω, Ω 1 , Ω 2 respectively, and Γ 1, Γ 2 the internal boundaries of Ω 1 , Ω 2 respectively. The Schwarz approach consists of defining two sequences u 1 k+1 and u 2 k+1 with k ≥ 0, and solving iteratively the two boundary-value problems: with u 1 k+1 = 0 on ∂Ω U ∂Ω 1 and u 1 k+1 = u 2 k on Γ 1 (22) and For the algebraic formulation of the Schwarz method one can think that the domain Ω is subdivided by M sub-domains Ω i with i = 1,…, M. We define the prolongation operator R i T , that extends by zeros the degree of freedom outside Ω i, and the restriction operator R i like its transposed. By denoting A h and A i the stiffness matrices with respect to the global and local bi-linear forms a h (·, ·) and a i (·, ·) and defining the operator Q i like the matrix Q i = R i T A i -1 R i , then after some algebraic manipulations the one level additive Schwarz iterative method reads for M ≥ 2 sub-domains: with k ≥ 0 index of iteration. Regarding (24) as a fixed-point relation, we can define the Schwarz method simply as a Richardson iterative method with preconditioner P as Unfortunately, the convergence of the one level Schwarz method deteriorates when the number M of sub-domains becomes larger because of the lack of information between sub-domains Ω i . In order to overcome this drawback, we introduce a global coarse mesh as a particular sub-domain Ω 1 , so that communications among all the sub-domains are guaranteed. Therefore, the final form of the preconditioner is P asc Since the algebraic system of (21) is A h u h = b h , the preconditioned two level additive Schwarz system becomes: The number M of the sub-domains in which Ω is divided, decides the efficiency of the Schwarz preconditioner. In fact, the number of nodes belonging to each sub-domain depends on M and also on the sub-domains length; both the choices should be as shrewd as possible in order to reach a good balance among the dimension of the local stiffness matrices A i and the optimal computational time. Another choice influences the efficiency of the P as and P asc preconditioners, namely the level of fill-in allowed in the process of LU factorization to calculate the matrices A i -1 , the so-called Incomplete LU factorization (ILU) process. The algorithm of the implementation of the Schwarz scalable additive preconditioner can be seen in [15].

Numerical Procedure
In summary, the fractional step works following these steps: 1) The unit width discharge q n , the depth h n , the elevation ξ n and the concentrations c i n are given at time t = t n 2) Compute the velocity u n by (5) 3) Compute the velocity u n+1/3 by (6) 4) Compute the unit width discharge q n+1/3 by (7) 5) Compute the unit width discharge q n+2/3 by (17)  6) Compute the adjourned value of the elevation ξ n+1 by (18)  7) Compute the adjourned value of the unit width discharge q n+1 by (19)  8

Numerical Test with Analytical Solution
The correctness and efficiency of the developed model have been checked solving the following problem with analytical solution. In fact, the expression of the elevation and the discharge, similar to that of [25], are the following: giving rise to the source function: where: A = 2\ 3800 , e = cos(d)), f cos(d,) = sin(d)) , q = sin(d,) Two analytical solutions were also taken for the concentrations: therefore the source functions for the concentrations were respectively: 3 (), ,) = −u sin(u)) sin(u,) + u cos(u)) cos(u,) + Γ 3 u 0 sin(u)) cos(u,) + 3 (1 + sin(u)) cos(u,)) (31) 0 (), ,) = −u cos(u)) cos(u,) + u sin(u)) sin(u,) − Γ 0 u 0 cos(u)) sin(u,) + 0 (1 − cos(u)) sin(u,)) with u=q/h and F=4π/3800 For the analytical test we fixed α=0; moreover we choose µ=0, g=9.81, Γ 1 = Γ 2 =0.001, a 1 =0.0000174, a 2 =0.0115 and did not take into account decay effects. The domain considered was Ω=[0, 3800], partitioned in 950 elements of equal length; the number of nodes was 1901 for the discharge q and the pollutants c 1 and c 2 ; 951 for the elevation ξ. The transient studied was 10800 s and ∆t=100 s. The initial and boundary conditions were obtained by the analytical solutions and we imposed Dirichlet at inflow and Neumann at outflow for ξ, c 1 and c 2 ; while we imposed only Dirichlet at inflow for q. In figure 1 we report the distribution of analytical and numerical solutions of concentration c 2 at time 10800 s and in figure 2 the time evolution of analytical and numerical solutions of concentration c 2 at node 238 (x=237). For the other results and comments, we refers to [14]. concentration c2 at t = 10800s.

Sub-Critic Flow in a Channel
In order to solve problems defined in channels with variating cross section the equations (1) and (2) have to be updated: with Q total discharge, h depth of water, A f and B f cross section and width channel respectively. Following the fractional step seen in the section 2, the modified algebraic system reads: given at t n : Q n , A f n , B f n , h n = h 0 +ξ n , u n =Q n /A f n , u n+1/3 =u n (x f ) with x f =x-∆tu n (x I ) and x I =x-(∆t/2)u n (x) (35) Indicating by M L the lumped mass matrix, by A the stiffness matrix, by H 1 the diagonal matrix with h 1i,i =1- For the elevation ξ n+1 , indicating by M L and by A the lumped mass matrix and stiffness matrix respectively, by C T the matrix of the integral of the product between the basis function and its derivative, by G 2 the diagonal matrix with For the update value of the Q n+1 discharge, indicating by M L the lumped mass matrix of discharge, by C T the matrix of the integral of the product between the basis function and its derivative, by H 5 the diagonal matrix with h 5i,i = g(A f n ) i , and H 6 the diagonal matrix with h 6i,i = (Q n+2/3 ) i /(h n ) i : By means of the above method we solve the steady problem reported in [8,12]. The domain is Ω=[0, 3], the height is h = 1-z b -ξ 0 , with ξ 0 as figure 3, and bottom z b as: The channel width is given by: . otherwise  (38) and in order to save computational time, we use two different time intervals for Q and ξ, that is ∆t1=0.0002 for Q and ∆t2=0.005 for ξ. The system (38) was solved by means of Bi-CGSTAB iterative solver, with about 60 iterations for each time step. For this test we don't use preconditioning and the iterations are stopped when dif < ε, with dif=difference of two consecutive iterations and ε = 10 -12 .
In figures 3 and 4 are reported the graphics of analytical and numerical solutions for ξ and Q respectively and in figures 5 and 6 are indicated the time history of the errors relevant to ξ and Q respectively. The maximum differences err between the analytical and numerical solution, at the end of the transient, are: for the elevation ξ, err = 3.04*10 -5 and for the discharge Q, err = 4.73*10 -6 .    Some comments regarding the results are appropriate. The stability and the accuracy of the presented method are confirmed, even if an equal order scheme for spatial discretization has been used. The Q graphic of numerical solution shows the oscillations typical of the methods lacking of an exact balance between the flux gradients and the source terms. For what concerns the elevation ξ, the graphic of numerical solution is coherent with that of analytical solution, even if the error is worse than the previous one. We ascribe this minor accuracy to the value of ∆t2 used in (38). From both the error time histories, it is evident that the stationery flow has been reached after 100s (in fact both the errors have reached their asymptotic value).
Remark 4.1. The flow studied in this test is relevant to a sub-critic case, but it is well known that characteristics could meet difficulties in presence of hyper-critic flows (e.g. hydraulic jumps) [10]. In order to overcome this drawback, recently some FE methods based on discontinuous Galerkin total variation diminishing Runge-Kutta schemes [16] and Streamline Upwind Petrov-Galerkin schemes [4] have been published. However, in modeling SW equations the most difficulties come from the approximation of the convective flux and source terms so that many conservative numerical schemes have been developed, particularly in the framework of finite volumes and appropriate Riemann solvers [8,12].

The Novella River Pollution
The presented model was applied to a segment of the Novella River, flowing in the Non Valley, Northeast Italy. Given the assumptions inherent to 1D models, it was assumed that the selected segments was characterized by a constant width with stream path. The Novella River was chosen since data concerning the river depth and discharge for some periods and stations were available. The river is surrounded by apple orchards, to which some pesticides are applied in known periods of the year. For the simulations presented here, an almost linear 3500 m river segment was selected. Simulations were run for 28 hours (from November 2 at 00:03 to November 3 at 05:03, 2009) for which precise informations about the variation of river depth during a precipitation event were available from a monitoring station located in the Comune di Dambel. Measurements of depth, reported in figure 7, had a temporal resolution of 15 minutes [6]. River discharges, figure 8, were calculated using the statistically-derived formula: q=35.087h 2 -10.141h+1.1 (42) relating the discharge to river depth, measured at the monitoring station [6].     Three pesticides applied to apple orchards in the Non Valley and characterized by different physical-chemical properties (Table 1), Pirimicarb, Dithianon, Methoxfenozide, were selected for running the model. Chemical half-lives in water (33.3, 0.05 and 68 days, respectively) were used to update concentrations according to decay, as expressed by (16). Hourly chemical loading at the inflow deriving from water runoff in the Novella watershed were calculated for the three chemicals using SoilPlus [11], a dynamic multimedia model based on the fugacity concept developed to investigate the fate organic chemicals in the air/litter/soil system at the local scale. A 10 cm loamy soil compartment divided in 20 sub-layers (thickness: 0.5 cm) and starting from field-capacity conditions was simulated. Precipitation data for the simulation period were acquired from Meteotrentino web-side [18]. Despite realistic applications, rates were used for all chemicals and it was assumed that they were applied directly to the soil layer. It must be remarked that no attempt was made to depict the real agricultural scenario of the Non Valley. For example, only one soil type was used and slope was not considered in runoff calculations. The only intention was to create a realistic environmental scenario for model illustration. Chemical loading at inflow were obtained by multiplying the SoilPlus output values, i.e., hourly amounts of chemical per unit area, by the application area.
Remark 5.1. Reference temperature for vapor pressure is 25°C, while for the solubility is 20°C [24]; for Methoxyfenozide, given the lack of data, half-lives in sediment and water were assumed equal to the one in soil [13]. From these data, it was possible to calculate the hourly average river depth h and the concentrations c i , using the formula: where n i is the number of mole and w i is the atomic mass of each pollutant, while q t is the total discharge. The time evolution diagrams of concentrations at inflow (expressed in ng/L) are reported in figures 9, 10 and 11.
Finally, considering that the concentrations are reduced of 0.5 in 2877120 s, in 4320 s and in 5875200 s for Pirimicarb, Dithianon and Methoxyfenozide respectively, it was possible to find the compounds decaying ratios for ∆t=100s:  14), each one with about 135 nodes, moreover we add the domain Ω 0 with 80 nodes. The overlapping consisted of two elements (4 nodes). Keeping into account that the matrix of (20) is constant in time and that the matrix obtained from its multiplication with the preconditioner is a well-conditioned matrix, this last one was used in all the time steps. On the contrary, the right hand term of (20) changes at every time step so that it had to be multiplied at every step for the preconditioner. Consequently, the three preconditioned systems (20) were solved at each step by means of the Schwarz preconditioned iterative Bi-CGSTAB solver in a very efficient way, in fact the number of iterations for each of the three systems was less than 13.      By the graphics can be verified that the computed velocity field is sub-critical everywhere, since the Froude number Fr=0.23<1, and that the flow is coherent with the boundary conditions imposed. The approximation of the convective term by characteristics allowed a maximum number CFL=11. 25 Moreover, from figure 20, 21 and 22 it can be observed that the compounds are realistically spread out along the river and they are moving with respect of the river current. It can also see, for all the pollutants, the perfect coherence between the concentration value measured at a fixed point x of the spatial distribution relevant to instant t, and the concentration value measured at a fixed instant t in the time evolution relevant to the point x. The profiles of concentrations are closely related to the physical-chemical properties of the modeled chemicals, mainly log KOW and half-lives, which influences chemical runoff. For example, Pirimicarb shows a rapid decrease in concentration which is due to its extremely high solubility in water and low log KOW; during the first hours of the precipitation event, it is rapidly removed from soil by runoff water. In contrast, both Dithianon and Methoxyfenozide show a slowly decreasing concentration profile; this behavior can be ascribed to higher log KOW and low solubility. Moreover, Methoxyfenozide is characterized by a moderate persistence in the environment (Table 1).
It should be noted that the concentrations depicted in figures do not represent the actual pollution levels in the Novella River, since the reported values derive from a series of worst-case assumptions and have the only aim of illustrating the modeling approach presented in this work.

Conclusions
In this paper a deterministic numerical model for the prevision of the distribution of pollutants in rivers has been presented. The model is based, for the flow, on the 1D Saint Venant equations, i.e. on the hypothesis of constant vertical velocity and hydrostatic pressure, and for the compounds on convection-diffusion and decay equations. The most important numeric features of the model presented are: time advancing by a fractional step scheme that allows for decoupling the physical contributions, and consequently to solve the wave equation with a low computational cost; approximation of the convective terms by means of characteristics so that the computational kernels are reduced to elliptic like ones. The algebraic systems are solved by means of very efficient methods: that relevant to the elevation by a Thomas solver and those relevant to the concentrations by a Schwarz pre-conditioned iterative Bi-CGSTAB solver. Of course, studies relevant to rivers with significant transversal velocity shouldn't be afforded by this model, however for studies on channel flows or on flows of linear segments of rivers it appears to be a promising tool, especially in presence of very long transients. The new model has been tested on three cases which results confirm its accuracy and computational efficiency. The first and second tests have known analytic solutions and the third one is an application to the realistic case of chemical pollution of the Novella River. The expected computational efficiency has been confirmed, in fact the overall Novella transient took 11 minutes of elapsed time on a Toshiba Intel Core I5 Satellite Pro personal computer.