Deterministic Advection-Diffusion Model Based on Markov Processes

Deterministic Advection-Diffusion Model Based on Markov Processes João S. Ferreira1 and Manuel Costa2 Abstract: A new deterministic numerical formulation named DisPar-k based on particle displacement probability distribution for Markov processes was developed to solve advection-diffusion problems in a one-dimensional discrete spatial grid. DisPar-k is an extension of DisPar, and the major difference is the possibility of establishing a number of consecutive particle destination nodes. This was achieved by solving an algebraic linear system where the particle displacement distribution moments are known parameters taken from the Gaussian distribution. The average was evaluated by an analogy between the Fokker-Planck and the transport equations, being the variance Fickian. The particle displacement distribution is used to predict deterministic mass transfers between domain nodes. Mass conservation was guaranteed by the distribution concept. It was shown that, for linear conditions, the accuracy order is proportional to the number of particle destination nodes. DisPar-k showed to be very sensible to physical discontinuities in the transport parameters ~water depth, dispersion, and velocity!, showing that this type of problem can only be disguised by introducing numerical dispersion ~i.e., changing the Fickian variance!. DOI: 10.1061/~ASCE!0733-9429~2002!128:4~399! CE Database keywords: Markov process; Diffusion; Probability distribution. Introduction The numerical simulation of dissolved substances transport in natural systems has become an increasingly important tool used, for example, in water quality management and environmental impact assess of engineered facilities. The solution of these models has been based on the deterministic solution of the advectiondiffusion equation, by means of Eulerian models-EMs ~Hoffman 1992; Chapra 1997! and Eulerian-Lagrangian models ~ELMs! ~Neuman 1984; Baptista 1987; Oliveira et al. 2000!. However, these methods do not make explicit use of stochastic concepts, which can be seen as a disadvantage in the comprehension of physical processes involving randomness. In fact, the complexity involving the transport of a particle in a turbulent water body is so great that only its statistical consequences can be measured. This stochastic nature of particle motion in some systems encouraged the development of several numerical formulations in statistical physics. Some of these systems—e.g., Brownian motion, birthdeath processes, noise in electronic systems ~Gardiner 1985!— follow the principles of Markov processes. The master equation, which describes the motion of a particle by transition probabilities, is often used to represent these processes ~Gardiner 1985; 1 Environmental Systems Analysis Group ~GASA!, Univ. Nova de Lisboa, Quinta da Torre, 2825, Monte da Caparica, Portugal. E-mail: [email protected] 2 Environmental Systems Analysis Group ~GASA!, Univ. Nova de Lisboa, Quinta da Torre, 2825, Monte da Caparica, Portugal. E-mail: [email protected] Note. Discussion open until September 1, 2002. Separate discussions must be submitted for individual papers. To extend the closing date by one month, a written request must be filed with the ASCE Managing Editor. The manuscript for this paper was submitted for review and possible publication on May 9, 2000; approved on October 31, 2001. This paper is part of the Journal of Hydraulic Engineering, Vol. 128, No. 4, April 1, 2002. ©ASCE, ISSN 0733-9429/2002/4-399– 411/$8.001$.50 per page. Risken 1989; Van Kampen 1992!. The Fokker-Planck equation is a simplified version of the master equation and authors like Heemink ~1990! and Dimou and Adams ~1993! used it to develop stochastic particle models, assuming that the concentration in the transport equation can also represent a particle transition probability. These formulations are a good alternative to EMs and ELMs, since they do not require a spatial grid for transport simulations. However, when applied to large scale modeling the required computational effort makes their usage unsuitable for ordinary implementation. Costa and Ferreira ~2000! proposed a new numerical formulation called DisPar, where a discrete probability distribution for particle displacement was developed. The model concept was developed to a spatial discretization on cells instead of nodes, but the mathematical treatment is similar for both representations. It was considered that a particle had three destination nodes corresponding to a distribution with three probabilities. To evaluate these probabilities the authors used the particle displacement average and variance, which were developed as a function of modeling coefficients. These probabilities were used as deterministic mass transfer prediction between neighboring nodes. However, the Courant number represents a stability restriction since only origin and two neighboring nodes can be considered as particle destination. Also the use of only three destination nodes imposes numerical restrictions in the dispersion term. In this paper an extension of DisPar, called DisPar-k is developed, also based on the discrete particle displacement distribution over a time step. The major difference is the possibility of using a user-specified number of consecutive particle destination nodes, corresponding to a distribution with an identical number of probabilities. These transition probabilities are obtained by solving an algebraic linear system taking the particle displacement distribution moments as known parameters. For a specified number of destination nodes, it is necessary to use the same number of moments, choosing them by ascending order and starting at zero. The particle displacement distribution moments can be evaluated asJOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 / 399 suming a Gaussian behavior for the transition probabilities. The average is obtained from the random walk particle models ~Heemink 1990; Dimou and Adams, 1993! and the variance is considered Fickian. So, all the moments used in the formulation are computed as a function of these two statistical parameters. Another important feature of DisPar-k is the possibility of considering any group of consecutive domain nodes. Therefore, the particle displacement average is used to choose the computation nodes, avoiding stability issues related to Courant restriction. A similar process is made in the advection treatment of ELMs, with the important difference of resorting to spatial points noncoincident with the grid nodes, which is avoided in DisPar-k. As in Costa and Ferreira ~2000!, the development of DisPar-k aims to show that transport models based on Markov processes may represent an alternative to EMs and ELMs, since some underlying concepts may become more visible by means of particle individual analyses. The present paper has six sections besides this introduction. First, the DisPar-k concept is presented in detail, which is followed by truncation error and convergence analysis. Then, three numerical tests are performed to compare DisPar-k with analytical solutions and three other tests are carried out to assess the influence of nonlinearities in the methodology performance. A real data application is made using a Dutch-Rhine branch hydrodynamic model. The paper ends with some concluding remarks. Model Development This section presents the concept of particle movement in a discrete space and over a time interval. The concept is developed, leading to the mass transfer predictions between nodes, and therefore it is possible to obtain the concentration of particles in a generic node after a time step. Background A Markov process is defined as a stochastic process that has the following property: P ~ x n ,t n u x 1 ,t 1 ,...,x n21 ,t n21 ! 5 P ~ x n ,t n u x n21 ,t n21 ! (1) where P(x n ,t n u x 1 ,t 1 ;...;x n21 ,t n21 )5the transition probability of a particle to be in position x n at time t n if it was in all the positions x at some specific time ~i.e., if it was in x 1 at time t 1 and ... and in x n21 at time t n21 !; P(x n ,t n u x n21 ,t n21 ) represents the transition probability conditioned only by the particle spatial position at the previous time. Thus, in a Markov process the transition probability is solely dependent on the previous spatial position. The motion of a particle obeying this condition can be expressed by the master equation, a form often used in statistical physics. This equation represents a differential form of the Chapman-Kolmogorov equation, which expresses the fact that a particle initially positioned in x 1 at time t 1 will get to position x 3 at time t 3 via any middle position x 2 at time t 2 ~Van Kampen 1992!. Any transition probability for a Markov process obeys this equation. A possible way of writing the master equation is through the Kramers-Moyal expansion ~Risken 1989!: S ] P ~ x,t ! ~ 21 ! a ] a E d ~ x a ! 5 P ~ x,t ! ]t a! ]x a dt a51 ( D 400 / JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 (2) where P(x,t)5the probability of a particle to be in x at time t and E d (x)5the particle displacement expectation. This expression was meant to transition probabilities and, for that case, P 5conditional probability ~Van Kampen 1992!. Nevertheless, expression ~2! represents a valid relationship for random Markovian variables. The Fokker-Planck is a special case of Eq. ~2!, which assumes that all terms bigger than 2 are negligible. This equation has been used as the basis for particle-tracking formulations in transport models ~Heemink 1990; Dimou and Adams 1993!. These formulations track each particle individually, and use an analogy between the transport equation and the Fokker-Planck equation to obtain parameters E d (x) and E d (x 2 ) S E d ~ x ! 5 u1 D ]D D ]A 1 dt ]x A ]x E d ~ x 2 ! 52D dt (3) (4) where u5velocity; D5Fickian coefficient; and A5section area. In the mentioned random walk models, particle displacement is caused by an advective deterministic component E d (x) and by an independent, random component statistically close to the random and/or chaotic nature of time-averaged mixing. This random component characterizes the random forces and the mentioned authors used a temporal discretization of E d (x 2 ) to produce results. As it will be shown, in DisPar-k the particle displacement is implemented as a Markov process, transposing the particle motion in a continuous space to a discrete space. The discretization scheme allows us to build a deterministic model in opposition to the conventional particle tracking ones. Concept In the present paper, space is divided in a one-dimensional ~1D! uniform grid and a particle located in node i can move to a node x, over a time step Dt ~Fig. 1!. This particle displacement depends on u, D, and A, which are known parameters in all spatial nodes. To establish that relationship, two important particle displacement statistical parameters, average and variance, are defined with a spatial and temporal discrete nature. The average is obtained using expression ~3! S a i5 u i1 D ]D D i ]A Dt 1 ]x A i ]x Dx (5) and variance is assumed to be Fickian s i2 5 2D i Dt Dx 2 (6) where a i 5average particle displacement; s 2i 5variance particle displacement; in the discrete space and over a time step; and u i , D i , and A i correspond to the coefficient values at the particle origin node i. The integer part of parameter a i , defined by b i , is embedded in the formulation to compute the particle displacement distribution. For that purpose, b i is used as the central node of an userspecified particle possible destination, which corresponds to k neighboring nodes for each side as it can be seen in Fig. 1. Thus, there are 2k11 possible events in this scheme, each one having an associated probability. The particle displacement distribution is Fig. 1. Possible events for particle in time step Dt defined by the probability that a particle will move from node i to a destination node x, with xP $ i1b i 2k,...,i1b i ,...,i1b i 1k % . So, the array of 2k11 probabilities can be written in matrix notation as F P ~ i1b i 2k,n11 u i,n ! P ~ i1b i 2k11,n11 u i,n ! ] W i5 P ~ i1b i 1k21,n11 u i,n ! P ~ i1b i 1k,n11 u i,n ! G (7) ~ 2k11 ! where P(x,n11 u i,n)5probability that a particle will move from node i to node x ~i.e., probability that a particle located in i at t 5n, will move to x at t5n11!. The use of parameter b i to center the particle destination nodes is similar to the ELMs particle tracking. However, in DisPar-k only the grid nodes are used ~i.e., the grid does not move with the flow as in typical Lagrangian advection treatment!, which avoids mass errors that can occur due to interpolations and/or integration between domain nodes in ELMs ~Oliveira et al. 2000!. The particle displacement in a discrete space and over a time interval has a discrete probability distribution, which is mathematically defined by the moments centered at a generic spatial point. Considering point i ~particle origin node! as the spatial reference origin, the moments centered at i are expressed as E i (x v ), with v 51,2,3, . . . , where v represents the order of moments. Another important moment class is centered at the average distribution, which equals a i for the particle displacement. This class is expressed by E i @ (x2a i ) v # and the second-order ( v 52) corresponds to the particle displacement variance (s 2i ). Note that every zero-order moment E i (x 0 ) equals one. These discrete statistical parameters can all be evaluated by means of theoretical probability distributions such as the Gaussian distribution, which is further used in the present paper. M5 3 Thus, after the computation of the first 2k11-order moments ~including the zero order!, the DisPar-k methodology allows us to evaluate the 2k11 probabilities associated with the particle possible destination nodes, as it is described afterwards ~i.e., to obtain 2k11 probabilities in a discrete distribution it is necessary to compute 2k11 moments!. To respect the Courant restriction for any k, the first-order moment must be lower than one ~considering this displacement always positive!. That purpose can be achieved by centering the particle destination nodes in i1b i . Therefore, the relationship between the discrete distribution moments centered at b i (E i @ (x2b i ) v # ) and the W i probabilities is, by definition, given as b i 1k E i @~ x2b i ! # 5 ( v x5b i 2k F G E ~ x2b i ! 0 E ~ x2b i ! 1 ] E i5 E ~ x2b i ! 2k21 E ~ x2b i ! 2k E i 5M W i ~ k21 ! 0 k0 2k 1 ~ 2k11 ! 1 ¯ 0 ¯ ~ k21 ! 1 k1 ] ]    ] ] ]    ] 2k 2k21 ~ 2k11 ! 2k21 ¯ 0 ¯ ~ k21 ! 2k21 ~ 2k11 ! ~ 2k11 ! (10) where M 5the square matrix, with (2k11)3(2k11) elements, given by ¯ 0 ¯ 2k (9) Thus remembering expression ~8!, it is possible to establish the following relationship between the matrices E i and W i : ~ 2k11 ! 0 2k (8) To evaluate the probabilities as E i @ (x2b i ) v # function, one can use matrix E i containing the first 2k11-order moments centered at b i , including the zero-order moment 2k 0 2k @~ x2b i ! v P ~ i1x,n11 u i,n !# ¯ 0 ¯ ~ k21 ! 2k k 2k21 k 2k 4 (11) ~ 2k11 !~ 2k11 ! JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 / 401 So W i , and therefore each probability, can be written as a function of the moment matrix E i W i 5M 21 E i The DisPar-k can be easily adapted to build distributions with an even number of particle destination nodes, since one can use the desired number of moments. This features was not included in this paper to facilitate the comprehension of the proposed method. (12) The distribution moments can be evaluated by means of a theoretical distribution. In the present situation, the particle displacement is approximate to the spatial continuous Gaussian ~or normal! distribution with average a i Dx and variance s 2i Dx 2 , over a time lapse Dt P ~ x,n11 u i,n ! 5 1 A2ps i2 Dx F exp 2 2 ~ x2a i Dx ! 2 2s i2 Dx 2 G Truncation Error Analysis In this section, the truncation error analysis was made for the linear problem, to simplify its numerical treatment. This analysis was carried out considering the probability of a particle to be at time n11 in a specific position as the model’s formulation base. The probability is expressed as the sum of the product of the possible origin probabilities by the transition ones. Both sides of this equation are decomposed into a Taylor series and its result is then analyzed. This section ends with a subsection that aims to strengthen the results by studying an explicit finite difference method embedded in the proposed formulation. (13) Thus, the particle displacement distribution in a discrete space is characterized by a i and s 2i , from which it is possible to compute all higher-order moments presented in matrix E i ~expression ~9! r21 E i @~ x2b i ! # 5 v ( j50 E i ~ x2b i ! 5d i (14) E i ⌊ ~ x2b i ! 2 ⌋ 5s i2 1d i2 (15) Linear Formulation In the linear case if a particle is found in a spatial position at time n11, the described formulation implies that there are only 2k 11 particle possible origins. Thus, the model’s state equation can be seen as v! @ E ~ x2b i ! 2 2E i ~ x2b l !# j 2 j j! ~ v 22 j ! ! i 3E i ~ x2b i ! v 22 j (16) k where r5( v 12)/2, if z is even or r5( v 11)/2, if z is odd. Expression ~16! is proved in the Appendix, Theorem 2. Note that these statistical parameters only depend on the modeling coefficients Dt, Dx, u, D, and A. Expression ~16! can now be used to obtain the 2k11 probabilities, from which the mass transfer between nodes over a time step is directly evaluated. Thus, the mass transfer from i to x is simply given by the product of the node i particle mass at time n by P(x,n11 u i,n), which are variables that only depend on the model conditions at time n. S5 F P ~ i1b,n11 ! 5 Note that the index i is omitted in b, since the particle displacement average is constant for all nodes in the linear situation. In this section, W i matrix will be constituted by the moments centered at the origin i and not at i1b i , as was previously made. The matrix with the spatial coefficients earlier known as M is now replaced by a new one called S, which can be yielded as ¯ b0 ¯ ~ b1k ! 0 ~ b2k ! 1 ¯ b1 ¯ ~ b1k ! 1 ]  ] ~ b2k ! ~ b2k ! 2k ¯ ¯  ] 2k21 ¯ ~ b1k ! 2k21 b 2k ¯ ~ b1k ! 2k b Hence, the matrix composed by the moments centered at i (E 8i ) can be represented by E i8 5SW i (19) which means that W i matrix is given by W i 5S 21 E i8 (20) With linear conditions the probability for particle displacement will only depend on the distance between the origin node and the destination one, and therefore it is possible to establish the generic equality P ~ i1b,n11 u i2q,n ! 5 P ~ i1b1q,n11 u i,n ! , qP $ 2k,..., 0,...,k % which means that W i can be rewritten as 402 / JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 (21) P ~ i1b,n11 u i2q,n ! P ~ i2q,n ! (17) ~ b2k ! 0 2k21 ( q52k G (18) ~ 2k11 !~ 2k11 ! F F P ~ i1b2k,n11 u i,n ! ] P ~ i1b,n11 u i,n ! W i5 ] P ~ i1b1k,n11 u i,n ! P ~ i1b,n11 u i1k,n ! ] P ~ i1b,n11 u i,n ! 5 ] P ~ i1b,n11 u i2k,n ! G G ~ 2k11 ! (22) ~ 2k11 ! Let C i be the matrix of probabilities c i 5 @ P ~ i1k,n ! ¯ P ~ i,n ! ¯ P ~ i2k,n !# ~ 2k11 ! (23) Thus, it is possible to express the Eq. ~17! as a function of C i and W i P ~ i1b,n11 ! 5c i W i where h x 5the first 2k11 spatial derivative orders, including the zero order; the coefficient matrix L (24) In the next two subsections, both sides of this equation will be developed into Taylor series relative to point (i1b,n), and they will be both truncated after the 2knd derivative term to show that they are equal. L5 Right-hand Side Development In this subsection all the terms present on the right side of Eq. ~15! will be developed into Taylor series relative to point (i 1b,n) and truncated after the 2knd spatial derivative. To perform this decomposition, one can consider the following matrices: h x5 F G ] 2k P ~ i1b,n ! ]x 2k 0 ¯ 0 0 0 1 1! ¯ 0 0 ] ]  ] ] ¯ 1 ~ 2k21 ! ! 0 0 1 ~ 2k ! ! 0 0 0 (25) F ~ 2b1k ! 0 ¯ ~ 2b ! 0 ¯ ~ 2b2k ! 0 ~ 2b1k ! 1 ¯ ~ 2b ! 1 ¯ ~ 2b2k ! 1 ]    ] ¯ ~ 2b2k ! 2k21 ¯ ~ 2b2k ! 2k ~ 2b1k ! ¯ 2k21 ¯ ~ 2b1k ! 2k ~ 2b ! c i 5h x LZ 2k21 ~ 2b ! 2k Thus, the c i matrix can now be written as (28) R v, j5 P ~ i1b,n11 ! 5h x LZW i 5h x LZS 21 E i8 (29) Hence, and by Theorem 4 expressed in the Appendix, it is possible to write Eq. ~29! as follows: ( v 50 ]vP ~ 21 ! v E d ~ x v ! v ~ i1b,n ! ]x v! H (27) ~ 2k11 ! 3 ~ 2k11 ! S D v P @ j,2j # 0 v ¹ @ j,2j # j D v 2 j ~ 2u ! j2 ~ v 2 j ! v2 j (33) The conversion from temporal to spatial derivatives is proved in Theorem 3’s demonstration from the Appendix and its general expression, written in a matrix format, can be expressed as ] jP 5h x R j ]t j (30) (34) Let h t 5the matrix of P temporal derivatives Left-hand Side Development F ]0P ]1P ] 2k21 P ~ i1b,n ! 0 ~ i1b,n ! 1 ~ i1b,n ! ¯ ]t ]t ]t 2k21 To prove that the left-hand side of the equation is equal to the right-hand side if truncated after the 2k derivative term, it is necessary to decompose it into a Taylor series relative to point (i1b,n). To achieve that, let us assume that the P variable can be represented by the linear Fokker-Planck equation h t5 ] P ~ x,t ! ] P ~ x,t ! ] 2 P ~ x,t ! 52u 1D ]t ]x ]x 2 and T the matrix (31) Let R j be the matrix F SD G where the first line is referenced by 0 and the nonzero terms begin at line j and end at line 2 j. R j ’s general term belonging to line v can be expressed as Replacing c i in the model’s Eq. ~24!, it is possible to write it as 2k ~ 2k11 ! 3 ~ 2k11 ! ~ 2k11 ! Z5 P ~ i1b,n11 ! 5 ¯ 0 4 (26) and Z matrix is expressed as ]1P ]0P ] 2k21 P i1b,n i1b,n ¯ ! ! ~ ~ ~ i1b,n ! ]x 0 ]x 1 ]x 2k21 3 3 1 0! SD j j R j 5 0¯ 0 D 0 ~ 2u ! j20 ¯ j D j ~ 2u ! 0 ¯0 G T (32) ~ 2k11 ! 3 ] 2k P ~ i1b,n ! ]t 2k G (35) ~ 2k11 ! F G Dt 0 Dt 1 ] T5 Dt 2k21 Dt 2k (36) ~ 2k11 ! JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 / 403 The P(i1b,n11) development into Taylor Series truncated after the 2k term and relative to point (i1b,n) leads to P ~ i1b,n11 ! 5 F ]0P ]1P i1b,n ! ~ ~ i1b,n ! ]t 0 ]t 1 3¯ ] 2k P ~ i1b,n ! ]t 2k G LT (37) ~ 2k11 ! Replacing the derivatives in expression ~37! using expression ~34! P ~ i1b,n11 ! 5n x @ R 0 R 1 ¯ R 2k # ~ 2k11 !~ 2k11 ! LT F Dt j S 2 j D D v 2 j ~ 2u ! j2 ~ v 2 j ! ( v j5 v 2 ~ r21 ! j! 1 j G (39) Multiplying by v ! and rearranging the sum from expression ~39! in inverse order, a new one can be yielded expressed as a function of the particle displacement average and variance, being also possible to verify that it is equal to (21) v @ E d (x v ) # . Therefore, it was proved that the formulation for the linear case respects the 2k terms from the Taylor series developments. These developments showed the evenness between both sides of Eq. ~17! if decomposed into Taylor series developments up to order 2k. This means that for 2k11 nodes used in the model’s formulation, the results obtained are very close to the best ones for an explicit numerical formulation. Spatial error cannot be considered the best spatial error for a specific number of destination nodes, since all higher terms from Taylor expansion introduce an error. Changing the first term of the sum from the right-hand side of the equation to the left-hand side and dividing both by Dt, one gets expression ~2! when Dt→0. However, in that expression all terms bigger than two vanish, and therefore only two terms from the Taylor expansion are necessary in this convergence situation. Particle Formulation as Truncation Error Evaluation The result shown for the left-hand side of Eq. ~15! expresses an important issue since it is a function of the particle displacement moments. Each derivative term has a corresponding moment associated with it. From the particle’s perspective, numerical errors represent the increase or decrease of the displacement moments. If P represents concentration, enlargements in these moments represent errors such as numerical dispersion or phase. For example, if the model is running with the smallest possible value that can be applied to k (k51) ~representing the linear DisPar 404 / JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 F G P ~ i,n11 ! 5 @ P ~ i11,n ! P ~ i,n ! P ~ i21,n !# DDt 1 uDt 2 Dx 2 2 Dx (38) Now, it is necessary to evaluate the number of nonzero terms present in a R matrix line ~i.e., the matrix with all submatrices R j !. To accomplish that, one must look at Rj’s expression and verify that the first nonzero term begins at j. This means that the last nonzero value in line v will be in column v , which is the first from this column. Assuming that r represents the amount of terms from line v not equal to zero the first entry can be given by v 2(r21). Therefore, so that a line v entry from matrix R may be different from zero, it must obey the condition: v 2(r21)< v <2 @v 2(r 21) # , which means that: r>1 and r<( v 12)/2. The first condition is universal and the second one imposes that the number of nonzero terms in line v is given by r5( v 12)/2 if v is even and r5( v 11)/2 if v is odd. Thus line v obtained from the product RLT can now be represented by v model!, there will be no numerical error up to the second derivative term ~i.e., the 2k derivative term!. Therefore, this three node concentration model has no numerical dispersion ~Costa and Ferreira 2000!. On the other hand, if one considers the forward time centered space model ~FTCS! ~see, for example, Hoffmam 1992 for details! it is possible to see that the model has numerical dispersion. Thus, writing the FTCS model as 12 3 DDt Dx 2 (40) DDt 1 uDt 1 Dx 2 2 Dx The numerical error can be analyzed calculating the parameters E(x) and E(x 2 ) from the probability matrix. E(x) is equal to the one picked up from the Gaussian curve, but E(x 2 ) has a decrease of E num52 S D uDt Dx 2 (41) where E num represents the numerical error associated with the second derivative term. This scheme can be used to calculate the errors coupled to other numerical models, but in the present paper it is not a goal to develop this issue deeply. It was used, as mentioned before, just to emphasize the importance of the proposed formulation. Convergence Analysis The convergence analysis was carried out considering the parameter b equal to zero for all the nodes in the domain. It was considered that if a particle is in position i at time n11, there are k possible origin-neighboring nodes for each side, besides the origin one. Thus, the model’s equation, for this specific case, can be written as k P ~ i,n11 ! 5 ( P ~ i,n11 u i2q,n ! P ~ i2q,n ! (42) q52k The 1D Fokker-Plank equation for nonlinear situations can be written as ]P ] 1 ]2 1 ~ v P !5 ~ 2D P ! ]t ]x 2 ]x 2 (43) The Fokker-Planck equation only has two parameters affecting the state variable P variation, which are the particle displacement average (vdt) and second-order moment, and therefore the convergence of the proposed formulation is shown only for k51 ~three nodes model!. It was assumed that if this situation is convergent, all the other formulations (k.1) are also convergent since they accomplish both moments and a few more, depending on k’s value. Increasing k will speed up the model convergence on the spatial component. Decomposing the probabilities ~P! from the right-hand side of Eq. ~42! into Taylor series relative to point ~i,n! and truncating them after the second derivative term, one can rewrite it as Table 1. Parameters and Conditions Adopted in Tests Parameters Dt Dx Total points x8 u(x) D(x) Time step number Initial condition C(x,0) C(0,t) c @ (s21) Dx,t # Maximum courant (uDt/Dx) Maximum dispersion coefficient (2DDt/Dx 2 ) Maximum Péclet number @ 2E i (x)/V i (x) # P ~ i,n11 ! 2 P ~ i,n ! 5 Test 1 linear Test 2 linear Test 3 nonlinear 24 200 64 0 10 0 25 Gauss hill, x 0 52000, d 0 5264 0 0 1.2 9.6 200 64 0 50 2500 10 Gauss hill, x 0 52000, d 0 5264 0 0 2.4 0.05 0.1 66 1x,u 0 50.1 0.003x 2 ,D 0 50.003 40 0 0 1.2 1.7 ` 4 33.5 D D D S D i11 22D i 1D i21 v i21 2v i11 Dt 2 1 2 2 v 1v Dt1 Dt ! ~ 2 21 i11 i21 Dx 2 Dx 2Dx 1 1 1! 1 1 2! FS FS S S D S D S 2 2 D i11 2D i21 Dt 1 v i11 2v i21 Dt 2 1 v i11 1v l21 ]P Dt 1 2 Dx Dx Dx 2 Dx Dx 2 Dx ]x D G 2 2 1 v i11 1v i21 D i11 1D i21 v i11 2v i21 ]2P Dt1 Dt 2 2 Dt Dx 2 1¯ 2 2 Dx 2 Dx Dx ]x 2 D S D dt 2 ]v ] 2D ]P ]D 2 2 dt P1 2 dt2vdt dt12v 2 ]x dx ]x ]x ]x 1D D G D G FS When both parameter Dx and Dt converge to zero, Eq. ~44! can be written as ] P5 S 1 0.062 3.8 ]2P dt ]x 2 (45) Since 2v 2 dt 2 /dx is one order higher (dt 2 !dt), this term vanishes from the equation. Dividing both sides by dt and rearranging Eq. ~45! the Fokker-Planck Eq. ~43! is obtained proving the model’s convergence for this particular case. (44) ditions are imposed: C(x,0)50 for x.x 8 , C(x 8 ,t)5C 0 for x <x 8 and C(`,t)50. The velocity field and the diffusion coefficient vary, respectively, linearly and quadratically with distance, i.e., u(x)5u 0 x and D(x)5D 0 x 2 and the section area is constant. So, the analytical solution becomes C ~ x,t ! 5 F S ln~ x/x 8 ! t ~ u 0 1D 0 ! C0 x8 erfc 2 x 2 AD 0 t 1exp S D S D u 0 ln~ x/x 8 ! ln~ x/x 8 ! 1t ~ u 0 1D 0 ! erfc D0 2 AD 0 t DG (47) Comparison with Analytical Solution The accuracy of the developed method was tested by two wellknown problems. The first problem, a linear situation, is a transport with the initial condition of a Gaussian profile, which has an average of x 0 and a standard deviation of d 0 . The boundary conditions imposed are C(0,t)5C(},0)50 and the analytical solution for this problem is given by C ~ x,t ! 5 d0 Ad 0 12Dt F exp 2 ~ x2x 0 2ut ! 2 2d 20 14Dt G (46) The second problem is a conservative transport of continuous injection where u and D have spatial variability. A methodology provided by Zoppou and Knight ~1997! is used to obtain this analytical solution where the following initial and boundary con- Two tests were done for the linear case with the Gaussian profile and a third one was carried out for the continuous injection with nonlinear conditions. Each of these tests has two models, showing their growing power in prediction capability. The values used in each test are summarized in Table 1. The first test was done to show the importance the number of nodes used (2k11) has on the accuracy of the results on a pure advective situation. Observing Fig. 2, it is possible to verify that the model with k56 produces more accurate results than the one with k51. This result corresponds to what was theoretically predicted, since the increase of nodes reduces the spatial error, which is the most important one introduced by the drift term. Increasing JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 / 405 Fig. 4. Results from DisPar-k in nonlinear situation ~Test 3! Fig. 2. Results from DisPar-k in pure advection situation ~Test 1! Nonlinear Water Depth Tests the Courant number is not a restriction since the spatial error will depend exclusively on the fractional part of the particle displacement average. On the other hand, the diffusion term is really dependent on the time step, meaning that temporal discretization can represent the most important issue in terms of accuracy. However, by increasing the number of nodes, this problem is expected to disappear as it can be seen in the second test ~Fig. 3!. A test closer to reality will be done now to better evaluate the formulation ~Fig. 4!. For the boundary treatment it was considered that b1k21 nodes to each side of the upstream and downstream boundaries influence the domain. This means that there are 2 (b1k21) hypothetical nodes with possible influence on the computational domain according to the boundary parameters. The values used in these possible mass origins are equal to the corresponding boundary and they were treated exactly in the same way as the domain nodes. The highest Péclet number can be found in the upstream node decreasing progressively to downstream. The results near this advection-dominated region are accurate in both models, reflecting the DisPar-k power to treat the advective term. However, downstream, it is possible to verify the instability produced by the three-node model. As it happens on the second test, the temporal error introduced by the diffusion term is extremely visible in this part of the computational domain. Once again the increase of the number of nodes used to compute the model at each time solved the problem and the results produced are remarkably accurate. Fig. 3. Results from DisPar-k in diffusive-dominated situation ~Test 2! 406 / JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 To test the DisPar-k model with spatial variable water depth, three tests are presented, where u50 and D5constant for all spatial points. In these conditions, a uniform initial concentration field should maintain the same values over any simulation time. The three tests represent three different water height profiles: the first is a function that represents a physical discontinuity ~if x<53 then y52; if x.53 then y56, Fig. 5!, where the derivative significantly changes. The second situation corresponds to a continuum function @y5 3 A (x250)15, Fig. 6# with an impossible derivative at a specific point (x550). The last situation is a fourth-order polynomial ~Fig. 7!, derivable at all points. In all three situations Dx51, D50.01, k51, the total simulation time is equal to 100 and the boundaries do not influence the results at the regions presented in the figures. The water depth spatial derivatives are approximated with a centered difference, since higher orders in the derivative calculation do not improve the results. As it was already pointed out, the conditions described above imply the theoretical conservation of the initial concentration field. If these concentration values change, it is possible to understand the influence of the water depth spatial variability in the DisPar-k numerical errors. The two first situations ~Figs. 5 and 6! show that concentration changes are not only dependent on temporal error, but also on the Fig. 5. Results for water depth function representing physical discontinuity provoke the spatial error presented in the other situations. Therefore, the decrease in the time step is sufficient to obtain accurate results. Real Data Application Fig. 6. Results for continuum water height function with nonderivable point water depth profile, since the time step decrease from 1 to 0,1 does not produce any improvement in the results. This problem is handled in Costa and Ferreira ~2000! due to a specific balancing for the dispersion flow ~i.e., the dispersion flow from point i to i11 is equal to the dispersion flow from point i11 to i!, which is achieved by a numerical diffusion introduction in the secondorder term. In the present paper, the average and variance imposed, respectively, in expressions ~5! and ~6! lead to the sort of errors presented in Figs. 5 and 6 when discontinuities and impossible derivatives are presented in the water depth data. Nevertheless, these hydrodynamic features are generally associated with mass imbalance errors ~Oliveira et al. 2000!, and therefore accurate solutions in transport simulations are very difficult to obtain. In random walk particle tracking models ~Heemink 1990; Dimou and Adams 1993! these problems are expected since the advective term includes the water depth spatial derivative. From a practical point of view, one can conclude the need of another spatial dimension, in this case the vertical dimension, in order to model the transport process correctly. In Fig. 7 it is possible to verify that the water depth profile with the polynomial function, which is always derivable, does not In this section two tests will be carried out using a hydrodynamic model with real data in a steady-state situation. By doing so, it will be possible to test the model in a practical case and therefore compare what was predicted theoretically in the previous sections with what will happen in practical cases. The first test aims to evaluate possible mass imbalances in the transport model. Thus, a uniform concentration field will be applied to the entire domain and it will be evaluated after some time. The goal of the second test is to appraise and reinforce the idea that numerical dispersion must be added in the model formulation, so as to correct mass imbalances caused by discontinuities. For that purpose, a comparison with the first version of DisPar ~Costa and Ferreira 2000! will be made by an instantaneous spill of mass. Both tests will run with a small Dt since the goal is to show problems due to the discontinuities in the particle displacement average derivatives. As it was done in the previous section, each derivative from the average term @Eq. ~5!# was calculated by a centered differences scheme. In the first test the parameters used were Dt 50.01 s and time step number5100 and in the second one Dt 51 s and time step number51,000. The number of destination cells was 11 (k55) for both tests. The case study was applied to a Dutch-Rhine branch called the River Waal. The Waal part in the study is located between 900 and 910 km relative to the Rhine datum. The hydrodynamic results were obtained from SOBEK, a computational 1D river model developed by the Delft Hydraulics and the Institute of Inland Water Management and Waste Water Treatment ~RIZA! of the Dutch government. The results were obtained with a dominant discharge of 1600 m3 s21 with a constant Dx of 99.58 m and can be seen in Fig. 8. The data used for the model calibration was recollected in the years 1995 and 1996. The hydrodynamic simulation was performed with a constant section width of 271.00 m except for section 15, where the value was 298.00 m. The dispersion term was calculated using the well-known Fischer’s formula ~Fischer et al. 1979! D50.011 Fig. 7. Results for continuum water height function with all points derivable U 2B 2 HU * (48) where D5dispersion coefficient; U5velocity; B5width; H 5mean depth; and U ! 5is the shear velocity @U ! 5(gHS) 0.5; g 5acceleration due to gravity; S5channel slope#. As it was explained in the previous section, DisPar-k is very sensible to noncontinuous derivatives in the average term @Eq. ~5!#, which happens with the dispersion derivative. To assess its importance the dispersion variability is assuaged by redefining each point value as the average of its own and the two neighbors. In Fig. 9 it is possible to observe that the dispersion peak decreases from 2,401 to 1,722 m2 s21, which is a fall of almost 30%. The smoothed dispersion variability is clearly much softer. The first test with constant concentration undoubtedly shows the mass transfer imbalances in the region where the parameters have more spatial variability ~Fig. 10!. The second test ~with smoothed dispersion! also shows this type of unsteadiness, but in a much thinner scale ~Fig. 10!, which shows the importance of noncontinuities in this type of model. JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 / 407 Fig. 8. River Waal profile ~water level, bed level, and velocity! As it was strengthened in the previous section, this unsteadiness can only be disguised by introducing numerical error in the particle displacement variance @i.e., changing the Fickian variance imposed on Eq. ~6!#. The spill of mass test ~Fig. 11! clearly shows this issue, since DisPar ~Costa and Ferreira 2000! has a higher peak than the two tests with DisPar-k. These differences in the tests occur in the small imbalances near the peak of the distribution where the dispersion coefficient was not smoothed. Conclusion This paper described DisPar-k, a new deterministic numerical formulation based on Markov processes, consisting in the development of a particle displacement probability distribution in a discrete space. Therefore, DisPar-k is an explicit scheme with a userspecified number of particle destination nodes allowing, at least for linear situations, to obtain the desired spatial numerical error. The overlap of temporal Courant restrictions and the control of spatial accuracy lead to excellent results in linear advectiondominated situations. In the numerical tests the spatial accuracy is achieved with a few particle destination nodes, since the spatial error can be corrected up to a very significant order. Thus, the typical problems in EMs related with numerical dispersion or instability are overtaken by the DisPar-k formulation. The diffusion component is strongly dominated by the temporal error and Fig. 9. Dispersion coefficient profile for two situations: directly obtained from expression ~48!; averaged dispersion 408 / JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 as it was strengthened in the numerical tests, this issue can only be solved by increasing the number of particle destination nodes. However, the discontinuities in the physical parameters ~velocity, Fickian number, and section area! lead to numerical errors that can only be accurately handled by studying two and/or three spatial dimensions. Mass conservation is guaranteed, since only grid nodes are incorporated in the computations. Thus mass errors that occur, for example, in the ELMs due to interpolations and/or integrations are avoided. This particle concept can be applied to unstructured meshes and also extended to two and three dimensions, without loss of perception of the physical phenomena analyzed. As it was exemplified, the DisPar-k formulation can be applied to evaluate truncation errors from other numerical methods with a spatial discrete nature. Finally, it is possible to conclude that the explicit use of stochastic concepts can help to understand and solve numerical problems in transport modeling. Acknowledgments The writers wish to acknowledge the contribution, comments, and suggestions of António Câmara, André Fortunato, António Carmona Rodrigues, and M. Rosário Cristóvão. A special acknowledgement goes to Mário Franca who provided the SOBEK hydro- Fig. 10. Results obtained with initial concentration of 1 in entire domain ~Dt50.01; time steps5100! E~ xv!5 E xv 1 A2V ~ x ! p S exp 2 D @ x2E ~ x !# 2 dx 2V ~ x ! (53) Integrating this expression by parts it is possible to express the moment of order v 12 as function of the two earlier ones and thus E ~ x v 12 ! 5E ~ x ! E ~ x v 11 ! 1 ~ v 11 ! V ~ x ! E ~ x v ! Fig. 11. Results obtained for spill of mass in cell 11 ~Dt51; time steps51,000! (54) To demonstrate the equality ~54! let us assume, for example, that v is zero in expression ~54! and replace the independent variable x by a new one centered on average. In this case the first product of the right-hand side is always zero, since the independent variable is of odd order and is centered on average. Now, using expression ~52! to get both expectations on both sides of Eq. ~54! it is possible to verify that both are in fact equal, which means that Eq. ~52! is true for that case and therefore for all the others, proving the theorem by induction. Theorem 2 dynamic model with data supplied by the Rijkswaterstaat/RIZA. This work is supported by Fundação para a Ciência e a Tecnologia from the Portuguese Ministry of Science and Technology under research Contract No. MGS/33998/99-00. If x is a random variable with Gaussian distribution the Gaussian moment of order n can be yielded as a function of average and variance replacing expression ~52! in expression ~51! r21 E~ xv Appendix: Proof of Theorems This Appendix shows some developments, which made it possible to achieve some of the stated results made in the present paper. These developments include the formulation and demonstration of four theorems. Gaussian Distribution For any distribution it is possible to express its moments of order n as follows: E ~ x v ! 5E @ $ x2E ~ x ! 1E ~ x ! % v # (S D j50 Theorem 3 The linear Fokker-Plank equation can be expressed as ]P ]P ]2P 52u 1D 2 ]t ]x ]x (50) All odd terms from E @ (x2E(x)) i # 50 which means that expression ~50! can be rewritten as r21 E~ xv!5 ( j50 S D (S D ]vP ] v1 j P v j v2 j 5 2u D ! ~ v ]t ]x v 1 j j50 j (51) where r5( v 12)/2 if v is even or r5( v 11)/2 if n is odd. To get the Gaussian moment of order v expressed only as a function of average and variance, two theorems will be formulated and shown. (56) which means that the temporal derivative of order v converted to spatial derivatives can be expressed as v v 2j v 22 j 2 j E @ $ x2E ~ x ! % # $ E ~ x ! % (55) Demonstration: To demonstrate this theorem it is possible to use again the expectation relationship ~54!. Thus, to perform the demonstration by induction let us assume that v 51 and use expression ~55! to get the right-hand side as a function of average and variance. The result obtained is formally equal to the result produced by expression ~55! for the third-order moment. An even order can also be applied to the left-hand side of Eq. ~54! and verify that both sides are equal. Thus, it was proved by induction that E(x n ) can be expressed as a function of average and variance like in expression ~55!. v v j v2 j j E @ $ x2E ~ x ! % # $ E ~ x ! % v! @ V ~ x !# j @ E ~ x !# v 22 j ( j j50 2 j! ~ v 22 j ! ! (49) Decomposing expression ~49! according to the binomial theorem yields a new expression as E~ xv!5 !5 (57) Demonstration: To demonstrate this theorem the derivative of order v 11 will be obtained from the one of order v and therefore S D F( S D ] ]vP ] 5 v ]t ]t ]t v j50 ] v1 j P v j v2 j 2u D ! ~ j ]x v 1 j G (58) Calculating this derivative expression ~58! can be yielded as Theorem 1 S D E @ $ x2E ~ x ! % 2 j # 5 ~ 2 j !! $V~ x !% j 2 j j! (52) Demonstration: The Gaussian moment of order v is written by definition as S D ] v 11 P ]2 ]vP ] ]vP 1D 2 2 11 52u v v ]t ]x ]t ]x ]t v If x is a random variable with Gaussian distribution it is possible to establish the following relationship: (59) For example, if v is equal to 1, expression ~59! can be written as ]2P ]2P ]3P ]4P 2 2 2 5u 2 22uD 3 1D ]t ]x ]x ]x 4 (60) JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 / 409 what is true, proving the theorem. Theorem 4 Let l be the diagonal matrix l5 3 1 0 ¯ ¯ 0 0 0 21 ¯ ¯ 0 0 A A   A A A A   A A 0 0 ¯ ¯ 21 0 0 0 ¯ ¯ 0 1 4 (61) ~ 2k11 ! 3 ~ 2k11 ! If Z and S are the matrices presented in the truncation error analysis section then l5ZS 21 (62) Demonstration: Let B be the matrix B5lS (63) (64) where l ii 5diagonal entry from matrix l. Entry z i j is equal to b i j , and so Z5B. Writing this equality as lS5Z ai 5 di 5 Dt 5 Dx 5 « 5 hx 5 (65) l 5 r 5 and multiplying both sides of the equation by S 21 , it is possible to verify that l5ZS 21 proving the theorem. Notation The following symbols are used in this paper: A 5 section area; B 5 matrix product; C(x,t) 5 particle concentration field; Ct 5 computational time; D 5 Fickian coefficient; D 0 5 constant diffusion coefficient; d 0 5 standard deviation of initial Gaussian profile; E d (x) 5 particle displacement expectation; E i 5 moments centered at i1b i node matrix; E i (x v ) 5 v th-order moment centered at origin node for particle displacement distribution; E i8 5 moments centered at i node matrix; E num 5 numerical error associated with second derivative term; i 5 particle origin node; k 5 constant characterizing the number of particle possible destination nodes; L 5 coefficient matrix; M 5 coefficient matrix; P(x,n11 u i,n) 5 probability that particle will move from node i to node x over time step; P(x,t) 5 probability of particle to be in x at time t; 410 / JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 5 5 5 5 5 5 5 5 5 5 5 5 5 5 bi 5 Multiplying the column j from matrix S by the line i from matrix l, it is possible to write the entry b i j from matrix B as b i j 5l i j s i j 5 ~ 21 ! i21 s i j 5z i j Rj S t tn u u0 v Wi x x8 xn x0 Y Z (66) s 2i 5 Ci 5 v 5 coefficient matrix; coefficient matrix; time; generic temporal point; velocity; constant fluid velocity; moment order; transition probability matrix; spatial independent variable; allocation of boundary condition; generic spatial point; center of mass of initial Gaussian profile; time step 2k order matrix; (2k11)(2k11) element matrix centered at b i ; average particle displacement over time step; integer part of average particle displacement; fractional part of average particle displacement; time step; spatial resolution; absolute sum of differences between numerical models and analytical solutions, matrix with 2k spatial derivatives for P(x,n); coefficient matrix; coefficient associated to displacement moments; variance particle displacement; probability matrix; and average particle displacement. Subscripts i 5 discrete space index; and n 5 discrete time index. References Baptista, A. M. ~1987!. ‘‘Solution of advection-dominated transport by eulerian-lagrangian methods using the backwards method of characteristics.’’ PhD dissertation, Massachusetts Institute of Technology, Cambridge, Mass. Chapra, S. ~1997!. Surface water quality modeling, McGraw-Hill, N.Y. Costa, M., and Ferreira, J. S. ~2000!. ‘‘Discrete particle distribution model for advection-diffusion transport.’’ J. Hydraul. Eng., 126~7!, 525–532. Dimou, K., and Adams, E. ~1993!. ‘‘A Random-Walk particle trackin model for well-mixed estuaries and coastal waters.’’ Estuarine Coastal Mar. Sci., 37, 99–110. Fischer, H. B., List, E. J., Koh, R. C. Y., Imberger, J., and Brooks, N. H. ~1979!. Mixing in inland and coastal waters, Academic, N.Y. Gardiner, C. W. ~1985!. Handbook of stochastic methods for physics, chemistry and the natural sciences, 2nd Ed., Springer- Verlag, Berlin. Heemink, A. ~1990!. ‘‘Stochastic modelling of dispersion in shallow water.’’ Stochastic Hydrol. Hydraul., 4, 161–174. Hoffman, J. ~1992!. Numerical methods for engineers and scientists. McGraw-Hill, N.Y. Neuman, S. P. ~1984!. ‘‘Adaptative eulerian-lagrangian finite element method for advection-dispersion.’’ Int. J. Numer. Methods Fluids, 20, 321–337. Oliveira, A., Fortunato, A. B., and Baptista, A. M. ~2000!. ‘‘Mass balance in eularian-lagrangian transport simulations in estuaries.’’ J. Hydraul. Eng., 126~8!, 605– 614. Risken, H. ~1989!. The Fokker-Plank equation—Methods of solution and application, 2nd Ed., Springer-Verlag, Berlin. Van Kampen, N. G. ~1992!. Stochastic processes in physics and chemistry, 2nd Ed., North-Holland, Amsterdam, The Netherlands. Zoppou, C., and Knight, J. H. ~1997!. ‘‘Analytical solution for advection and advection-diffusion equations with spatially variable coefficients.’’ J. Hydraul. Eng., 123~2!, 144 –148. JOURNAL OF HYDRAULIC ENGINEERING / APRIL 2002 / 411