Impact of the Finite Element Mesh Structure on the Solution Accuracy of a Two-Dimensional Kinematic Wave Equation

: This paper presents the inﬂuence of the ﬁnite element mesh structure on the accuracy of the numerical solution of a two-dimensional linear kinematic wave equation. This equation was solved using a two-level scheme for time integration and a modiﬁed ﬁnite element method with triangular elements for space discretization. The accuracy analysis of the applied scheme was performed using a modiﬁed equation method for three different uniform triangular meshes with the same resolution, but with a different structure. The modiﬁed equation approach based on the Taylor series truncation allowed the numerical diffusivity and dispersivity tensors to be derived, which are directly associated with numerical errors. The derived tensors depend on parameters such as the space and time interval, ﬂow velocity, and weighting coefﬁcients. A detailed analysis carried out for the particular values of these parameters enabled an assessment of the numerical errors that may be generated in the solution for the assumed mesh structure. The theoretical analysis was conﬁrmed by using numerical simulations carried out for an arbitrary domain and auxiliary conditions. According to the obtained results, it appears that it is possible to improve the accuracy of the numerical solution by choosing the proper mesh structure and numerical parameters for the applied algorithm.


Introduction
The kinematic wave model, being a simplification of the dynamic wave model, is widely used in hydrology and environmental engineering. This model, in the form of a one-dimensional equation, usually describes a simulation of unsteady flow in an open channel [1,2], whereas the two-dimensional form of the kinematic wave equation was developed for overland flow routing [3][4][5]. Moreover, the kinematic wave equation is often considered as a component of more complex rainfall-runoff modeling systems [6][7][8][9]. The properties, applicability and solution of the kinematic wave equation have been described in many publications, and a comprehensive discussion on these problems is presented by Szymkiewicz [2], Miler and Cunge [10], and Singh [11].
It is worth noting that the kinematic wave equation does not include a term representing diffusion processes. Therefore, this equation can be considered as a pure advective transport equation with a numerical solution that presents a significant problem. In order to solve this equation, any algorithm of the finite difference method (FDM), finite element method (FEM), or finite volume method (FVM) can be used. In the two-dimensional case, it is convenient to use FEM or FVM because of the better fit to the boundaries of the area under consideration.
It is known that the numerical algorithms used to solve a differential equation are frequently inaccurate due to the numerical errors present in the solution. Depending on the approximation order of derivatives, numerical schemes generate errors changing the amplitude and phase celerity of waves representing the solution. It should be noted that the attenuation of the wave amplitude is caused by dissipative (diffusive) schemes that generate errors of numerical diffusion, while dispersive schemes introduce dispersion errors that change the phase velocity manifested by non-physical oscillations [12][13][14][15][16].
The solution of the advection equation is a particular challenge due to the fact that the physical dissipation process is not represented in this equation. Therefore, it is not possible to suppress the oscillations caused by the numerical dispersion resulting from the applied scheme [12,17]. On the other hand, the effect of dispersivity can be mitigated by using a scheme generating numerical diffusion [13,18]. However, it should be remembered that, in this case, the introduced diffusion will cause excessive wave attenuation.
The numerical errors depend on the applied algorithm, the size of the adopted elements (or volumes), and the numerical parameters such as the time step or the weighting parameters of the method. While analyzing such algorithms, the quality of the obtained solution is most often examined due to: (1) the shape of the elements (triangular, square, or other) and (2) the resolution of the mesh, which determines the size of both the spatial step and the time step according to the stability and accuracy of the solution. An important element of the numerical analysis is also the influence of the mesh structure on the accuracy of the obtained solution. This aspect was discussed by Gresho and Sani [6] in the context of the linear advection-diffusion equation. Numerical experiments show that even structured meshes with the same resolution but with a different arrangement are responsible for generating solutions with numerical errors of various sizes. Consequently, a numerical solution obtained with the same algorithm may be unstable or stable with variable accuracy. It should be noted that for a given mesh structure, the accuracy of the solution can also depend on the preferential flow direction. Taking into account this additional aspect means that the correct selection of the mesh with a given resolution and structure is not a trivial task.
The accuracy analysis of commonly used algorithms for the solution of a 1D advection equation is well established. However, in the case of 2D problems, due to the more complicated structure of truncation errors, accuracy and stability analyses are most often performed only by use of numerical experiments. In this approach, the numerical algorithm is inferred by analyzing the behavior of the solution error when the grid is refined. In other words, by reducing the spatial and temporal step, it is checked whether the obtained approximate solution coincides with the corresponding exact or reference solution. Unfortunately, such an approach does not provide complete information about the properties of the numerical scheme.
The estimation of numerical errors can be made with regard to a phase error analysis using the von Neumann method [5]. Such an approach was applied for a 2D advection equation [19] as well as for a 2D advection-diffusion problem [20], where equations were solved using the finite difference scheme on the staggered grid. However, an analysis carried out in this way directly reveals only the properties related to the stability and dispersivity of the applied scheme.
An alternative analysis of numerical errors can be performed with the modified equation method proposed by Warming and Hyett [21]. This approach, based on the Taylor series truncation, provides full information about the numerical dissipation and dispersion errors that can be generated in the solution [17]. In the context of a 2D advection equation, such an accuracy analysis was applied for the finite difference scheme [22] as well as for the dimensional decomposition using the 1D finite element method with a rectangular mesh [23]. An accuracy analysis performed using the modified equation approach was also presented by Szymkiwicz and Gąsiorowski [24] for a 2D diffusion wave equation solved using the finite element method for the floodplain inundation problem.
This paper presents the influence of the mesh structure on the accuracy of the numerical solution of a 2D kinematic wave equation, which can be considered as a pure 2D advection equation. The accuracy of the solution was assessed for three different uniform triangular meshes with the same resolution but with a different structure. Such an analysis using the modified equation approach has not yet been published. The governing equation was solved by means of the two-level scheme and the modified FEM with triangular elements. According to the performed accuracy analysis, the derived diffusivity and dispersivity tensors enable the properties of the algorithm to be explained in detail with regard to the generated numerical errors.
Many works provide only information about the approximation order used in an applied algorithm, which is not complete information to properly assess the obtained solution with regard to the solution's accuracy. Therefore, in this paper it was proved that a more precise numerical analysis is possible that does not require many time-consuming numerical tests.
In view of modeling practice for two-dimensional flow problems, the study shows that the accuracy of the solution can be significantly improved when a proper mesh structure is applied. Moreover, it is also possible to minimize the numerical errors by appropriately orienting the coordinate system with regard to the main flow direction.

Two-Dimensional Kinematic Wave Equation
The kinematic wave equation for overland flow in the 2D case can be described by the following equation [25]:

∂H ∂t
where x, y-space coordinates; t-time; H-water depth; n m -manning roughness coefficient; s-bottom slope along the flow direction; s x and s y -bottom slope components along the xand y-directions, respectively; φ-source term (φ = P· cosβ − f ) involving rainfall (P) and infiltration (f ); β-inclination angle of the slope (in degrees). The kinematic wave equation (Equation (1)), describing the advective transport, is a partial differential equation of the hyperbolic type, written in so-called conservative form [26,27]. Moreover, this equation is nonlinear with regard to the unknown function H. However, it should be noted that the accuracy analysis of the scheme used to solve this equation can only be performed for its linear form without the source term φ. In order to linearize Equation (1), it can be rewritten as follows: where u and v represent the kinematic wave celerity in the xand y-directions, respectively. These parameters are expressed as follows: For the constant values u and v (assuming that H is the average water depth) the governing equation becomes the 2D linear kinematic wave equation, which can also be considered as a 2D linear advection equation.

Solution of the Kinematic Wave Equation Using the Finite Element Method
Assume that the continuous domain S is covered by a mesh with M triangular elements ( Figure 1a). The finite elements with area A are joined in a finite number of points N, which are also computational nodes where an approximate solution will be obtained. For a solution to Equation (2), an algorithm of the modified FEM with a linear base function in the Galerkin formulation was applied [24]. The modification of the FEM provides a more general approximation formula, in which, by an appropriate selection of the parameter values, the instabilities occurring in the case of the standard FEM can be eliminated [2]. Moreover, owing to the involved weighting parameter, such an algorithm is a more flexible and makes possible an analysis of the influence of the numerical errors on the solution.
According to the modified FEM algorithm, the solution of the governing equation satisfies the following equation:  Subscript "c" means that the function is constant in the element and equal to the average weighting value. When the linear basis functions for approximation are assumed, then in element e, only three non-zero integrals occur: with the average weighting values H c (t) defined as follows: • For Equation (7): • For Equation (8): • For Equation (9): where α is a weighting parameter from the range [0, 1].
The above approximations H a and H c applied for one term of the sum in Equation (5), after its integration, yield the following formulas: where the coefficients and area of the element are calculated according to the following formulas [28]: The calculation of all integrals in Equation (5) and the assembling of the obtained relations for all elements (e = 1, 2, . . . , M) leads to a system of ordinary differential equations (ODEs): where A-constant, banded, and symmetrical matrix; B-banded and variable matrix; H-vector of the nodal values of H(t). Such a system can be solved using a two-level difference scheme [2]: where θ-weighting parameter from the range [0, 1]; n-time level index; ∆t-time step.
Formula (24) applied for the obtained ODE (Equation (23)) finally yields the following system of algebraic equations: which must be closed by the imposed boundary conditions. The presented scheme includes two weighting parameters determining the numerical properties of the obtained solution. The parameter α is associated with the averaging of time derivatives in space. For α = 1/2, the scheme corresponds to the standard FEM with linear basis functions, whereas for α = 1/3 an approximation with an arithmetical average is obtained. The weighting parameter θ defines the averaging of spatial derivatives in time leading to commonly used algorithms for the solution of ODE. For θ = 1/2, the scheme (Equation (24)) becomes the implicit trapezoidal method, whereas for θ = 1, the implicit Euler formula is obtained. Moreover, the scheme is absolutely stable for θ ≥ 1/2 and α ≥ 1/3 [24]. Additional information associated with the influence of both parameters on the numerical properties of the applied scheme can be obtained from the accuracy analysis carried out using the modified equation approach.

Accuracy Analysis for FEM Using Various Meshes of Triangular Elements
The influence of the mesh structure on the accuracy of the numerical solution was assessed for the three different uniform triangular meshes shown in Figure 2. The first of the meshes was based on a pattern in which the inner point P connects to the adjacent points resulting in a six-node structure (6N) (Figure 2a). The second mesh also used the six-node structure, but it had an inverse orientation (6NI) (Figure 2b). The last mesh was based on eight external nodes (8N) (Figure 2c). Using the meshes defined in this way, the 2D linear kinematic wave equation (Equation (2)) was solved using the modified FEM with a two-level difference scheme. Let us consider a single rearranged equation involved in the system (Equation (23)) corresponding to an internal node of the 6N mesh structure presented in Figure 2a. The spatial approximation of this equation was performed at the node P, and then the obtained ODE was integrated using the Equation (24). Finally, as a result of the spatial and temporal approximation, the following algebraic equation was obtained: where: with ∆x and ∆y-spatial step in x and y, respectively ( Figure 2). The accuracy analysis of the applied scheme can be performed using the modified equation approach proposed by Warming and Hayett [21]. Although this technique is rarely used, it provides very detailed information on the numerical solution with respect to its accuracy and stability. The modified equation procedure is a process inverse to the approximation of the governing equation. In the case under consideration, all nodal values of the function H involved in Equation (26) were expanded in the Taylor series around the node P. It was assumed that in the Taylor series expansion the derivatives up to the 3rd order are taken into account. Then, time derivatives higher than the first order were eliminated by the appropriate substitution and algebraic conversions. The applied approach led to the following equation valid in each node: where B = 1 − 6θ + 6θ 2 , and Cx and C y are the Courant numbers in the xand y-directions, respectively, defined as follows: Components of the right-hand side of Equation (28) can be considered as elements of tensors, and then this equation can be rewritten in a more compact form: where i, j = 1,2; X 1 = x; X 2 = y; w 1 = u; w 2 = v; D is the numerical diffusivity tensor: whereas E is the numerical dispersivity tensor: Equation (31) represents the governing equation which was modified by the applied numerical method of the solution. The left side of this equation corresponds to the 2D linear kinematic wave equation (Equation (2)), while the right side shows all the Taylor series expansion terms that were omitted during the approximation. When the mesh dimensions tend to zero, then all of those terms disappear, and the approximating equation is consistent with the solved differential equation. It is also worth noting that the modified equation (Equation (31)) with E = 0 becomes the advection-diffusion equation.
In Equation (31), the diffusive and dispersive terms represent the numerical errors which can be observed in the solution. The term associated with the 2nd order derivative (tensor D) is related to a dissipation error, whereas the term associated with the 3rd order derivative (tensor E) is considered as a dispersion error. These truncation errors are generated in the numerical solution as the non-physical damping of the wave or as non-physical oscillations, respectively [5,12]. The type of error that will dominate the solution is determined by the first term of the truncated part of the Taylor series expansion (the term on the right side of Equation (31)). If the first term contains an even derivative, then the scheme will be dissipative, i.e., the solution will be overly smoothed without non-physical oscillations. On the other hand, if the first term contains an odd derivative, then oscillations will appear in the solution, proving the dispersivity of the scheme. In the case of the presented FEM scheme, the occurrence of possible oscillations or the smoothing of the solution will depend on the values of the weighting parameters θ and α, time step, spatial step as well as flow velocities u and v (see Equations (32) and (33)). Thus, depending on the values of these parameters, the applied method can be an approximation of the 1st, 2nd, and 3rd order, with regard to x and y as well as t. For values θ > 1/2, the scheme is an approximation of the 1st order, and it introduces the numerical diffusion into the solution. According to the tensor D (Equation (32)), the amount of numerical diffusion, apart from the value of parameter θ, depends also on the time step ∆t and the components u and v of velocity vector. When parameter θ = 1/2, then all elements of the diffusivity tensor D are deleted, and the scheme is an approximation of the 2nd order. Consequently, it does not generate numerical diffusion. However, the elimination of diffusion causes the truncation error of the Taylor series to be dominated by terms with 3rd order derivatives that are related to the numerical dispersion of wave speeds. In such a case, an oscillating solution according to the tensor E (Equation (33)) can be generated.
To obtain additional information about the properties of the applied scheme with the associated mesh structure, the rotation of the x-y coordinate system by the angle β, with respect to the origin of the coordinate system, is performed [22]. In the obtained new t-n coordinate system, the t-axis is tangent, and the n-axis is perpendicular to the flow direction indicated by the velocity vector w (Figure 3a). Consequently, the tensor elements lying on the main diagonal D 11 (or E 11 ) and D 22 (or E 22 ) correspond to the tangent and normal direction, respectively (Figure 3b). It was assumed that the sign convention for rotation is positive counterclockwise. The rotation of the coordinate system causes the tensors D and E to be transformed into the new system as follows [29]: where D and E-tensors in the x-y co-system; D r and E r -tensors in the new t-n co-system; R-rotation tensor in two dimensions: According to Equation (34), the numerical diffusivity tensor in the t-n coordinate system has the following elements: The numerical diffusion, apart from the parameter θ, depends also on the flow direction determined by the value of the angle β. While analyzing Equations (37)-(39), it can be noticed that the maximum numerical diffusion takes place when u = v, i.e., for the angle β equal to 45 o . In this case, the flow velocity is oriented along the diagonal in the tangent direction t. Consequently, using β = 45 o the numerical diffusivity tensor for the 6N mesh structure has the following form: where w is a modulus of the velocity vector (w = √ u 2 + v 2 ). According to the obtained relation only the tensor element in the tangent direction D 11 is non-zero. Thus, in this direction we should expect the maximum numerical diffusion. At the same time, the diffusion in the normal direction n (D 22 ) as well as the tensor components lying beyond the main diagonal (D 12 and D 21 ) disappear. This dependence is schematically shown in Figure 3b. Similarly, the dispersivity tensor can be transformed to the new system according to Equation (35) for β = 45 o and, finally, we obtain: where ∆s = ∆x = ∆y; C = C x = C y . As a result of the rotation (Equation (41)), the non-zero values of the tensor components are to be found only on the main diagonal, i.e., in the tangential (t) and perpendicular (n) directions to the flow direction. The same procedure was used to derive the numerical diffusivity and dispersivity tensors for the remaining mesh structures of triangular elements presented in Figure 2b,c. Numerical diffusivity and dispersivity tensors in the new t-n coordinate system after rotation with β = 45 o for various mesh structures were as follows: • For a six-node structure (6N): • For a six-node structure-inverse orientation (6NI): • For a eight-node structure (8N): As can be seen, for all mesh structures, the diffusion tensors D r take the same form (Equations (42), (44), and (46)). Since the component D 11 of the tensor D r for θ > 1/2 is a non-zero value, whereas the other components are equal to zero, then the applied scheme generates numerical diffusion only in the tangential direction. On the other hand, depending on the mesh structure, the dispersivity tensors E r have different components lying on the main diagonal. Thus, although meshes with the same triangular elements were used, the various orientations of triangular elements lead to schemes with different numerical properties.
For the value θ = 1/2, in the diffusivity tensor D all components are zero. Consequently, the scheme does not generate any numerical diffusion. In this case, the elimination of diffusion causes the solution to be dominated by numerical dispersion. It is worth recalling that then non-physical oscillations will be generated in the solution. According to the form of the dispersion tensors E r (Equations (43), (45), and (47)), these oscillations depend on numerical parameters such as weighting parameters, spatial step, time step, and flow velocity. Figures 4-6 present the values of the dispersivity tensor components E 11 and E 22 depending on the Courant number, C, which were calculated for θ = 1/2 using selected specific values of the weighting parameters α = 1/3 and α = 1/2. Such values of the parameter α adopted in the applied scheme correspond to the modified (α = 1/3) and the standard FEM (α = 1/2), respectively.    For the value of the weighing parameter α = 1/2 (Figure 4b), the dispersion error will converge to zero with decreasing values of the Courant number. Thus, by assuming low values of this number, it is also possible to significantly minimize the dispersion error. However, in this case, with a fixed value of the space interval ∆x, very small values of C also correspond to very small values of the time step ∆t. Consequently, the computation time is significantly extended, which makes such a mesh with α = 1/2 less useful in practical applications.
In Figure 5, the relationship E = f (C) for the mesh structure with inverse orientation (6NI) is shown. When analyzing these relations, it can be seen that for both α = 1/3 and α = 1/2, the dispersion error can be systematically minimized as the Courant number is reduced. In the case of α = 1/2 presented in Figure 5b, the same tendency was obtained as for the previous 6N mesh structure (Figure 4b), where for the smaller numbers of C, the dispersion error decreases taking the value of zero if C = 0, that is ∆t = 0. This property can be explained by the tensors of numerical dispersivity derived using a corresponding mesh structure. For the accepted values of θ = 1/2 and α = 1/2, the elements of these tensors (Equations (43) and (45)) for both the 6N and 6NI mesh structures take the same form as follows: Thus, the components E 11 and E 22 depend on the Courant number only and their absolute values decrease when C decreases as indicated in Figures 4b and 5b. On the other hand, for α = 1/3 (Figure 5a), the components E 11 and E 22 are different from zero for the entire range of variability of C. Only for the one value of C ≈ 0.41 is the component E 11 equal to zero. Thus, in both perpendicular and tangential directions, a numerical dispersion error will always be generated. Even for very small values of the time step ∆t, this error will reach certain constant values causing oscillations.
The relation E = f (C) obtained for the 8N mesh structure is shown in Figure 6. As can be seen in this case, regardless of the value of the Courant number, the components E 11 or E 22 of the dispersivity tensor take non-zero values for both variants of the scheme with α = 1/3 and α = 1/2. This error will decrease if at the same time the value of the Courant number decreases. However, even for very small values of C, this error cannot be eliminated. A detailed explanation of this feature is provided by the analysis of the dispersivity component E 22 (Equation (47)) for C = 0: The above relation shows that the dispersion error takes a constant value regardless of the weighting parameter α and the Courant number C. Consequently, when the numerical diffusion is eliminated then non-physical oscillations will always be generated in the solution, even for C tending to zero.
Comparing these results, it can be seen that for the same values of numerical parameters but with a different mesh arrangement, one should expect a different accuracy of the solution. It is worth noting that although the dispersivity of the applied scheme cannot be eliminated, its reduction is possible by the appropriate selection of the numerical parameters α as well as C. The presented analysis shows that the best numerical properties are obtained by using the modified FEM with α = 1/3 applied on the 6N mesh structure. On the other hand, the worst effects in the form of numerical dispersion errors always causing oscillations should be expected in the case of the 8N mesh structure, regardless of the values α and C. Numerical simulations presented in the next chapter confirm the results of the theoretical analysis.

Numerical Tests
The FEM algorithm using the considered mesh structures was applied to solve the 2D linear kinematic wave equation. To this order, the advection transport with a uniform velocity field u = v = const = 0.1 m/s was assumed. A domain of the dimensions 50 × 50 m was divided into triangular elements of the dimensions ∆x = ∆y = 1 m using the 6N, 6NI, and 8N mesh structures (Figure 2). Consequently, the domain was covered with 2601 nodes (N x × N y = 51 × 51 = 2601). At the boundaries of the domain, the Dirichlet boundary conditions were imposed with H(x = 0,y) = H(x = 50 m,y) = H(x,y = 0) = H(x,y = 50 m) = 0 for 0 ≤ t ≤ T. The initial condition H(t = 0,x,y) is described by the Gaussian distribution as follows: with the parameters σ = 2.5, µ x = µ y = 10 m, and H s = 1 m (Figure 7). The simulation was carried out up to the time t = T = 250 s. The imposed initial condition for the assumed uniform velocity field with u = v = 0.1 m/s means the initial distribution of the quantity H is advected with constant velocity along the diagonal direction. Moreover, if the solution is not affected by any dispersion and diffusion errors, then the distribution should remain unchanged at the final simulation time T.
All numerical tests were performed with software developed by the author using the Fortran programming language, while the results of the simulations were visualized using the scientific graphics software Surfer (version 9.0).

Solution Dominated by a Diffusion Error
In the first test, calculations were performed for flow dominated by a diffusion error. This situation occurs when in the applied scheme θ > 1/2 is assumed. Then, the diffusivity tensor component D 11 is different from zero and in the solution a diffusion error is generated according to the following relation: Numerical calculations were performed for various values of the parameter θ using a fixed time step ∆t equal to 10 s, which corresponds to the Courant number C = 1. In Figure 8 the results of a simulation using the 6N mesh structure with α = 1/3 are presented.
The results of the computations are presented in the form of the water surface elevation (left part of the figure) as well as isolines (right part of the figure). As can be seen, for θ = 0.55 ( Figure 8a) the initial symmetrical distribution of H was deformed during its propagation with a constant speed w = (0.1 2 + 0.1 2 ) 0.5 m/s. After the time t = 250 s, the introduced numerical diffusion causes a reduction in the wave amplitude and extension of the isolines. In other words, the initially symmetrical distribution was systematically smoothed over time. Assuming a constant time step and velocity flow, this deformation becomes greater as the weighting parameter θ increased (Figure 8b), and a maximum diffusivity occurred when the parameter θ was equal to 1 (Figure 8c). It is worth noting that the isolines were lengthened and smoothed only in the flow direction, which results from the fact that the diffusivity tensor D has only one non-zero component D 11 , which corresponds to the tangential direction (see Figure 3b and Equation (40)). According to Equation (52), a similar solution with a diffusion error can be obtained when the Courant number C or the time step ∆t is increased for the assumed constant value of the parameter θ > 1/2. Moreover, the same solutions were obtained for the other two mesh structures 6NI and 8N, because, as previously shown, the numerical diffusion tensors D have an identical form for each mesh (see Equations (42), (44), and(46)).

Solution Dominated by a Dispersion Error
In the next tests, calculations were performed for θ = 1/2, which eliminates numerical diffusion. Consequently, the solution will be dominated by numerical dispersion. As results from the previously presented analysis of the dispersivity tensor (Figures 4-6), the accuracy of the solution depends on the weighting parameter α as well as on the Courant number C.

Solution Using the 6N Mesh Structure
The solution obtained using the modified FEM with α = 1/3 for C = 1.5 had a strong oscillating character (Figure 9a). It is worth noting that the oscillations appeared behind the propagating wave and were mainly oriented along the velocity vector w, i.e., along the tangent direction in the transformed coordinate system. This effect is due to the fact that the tensor component in the tangential direction E 11 is twice as large as the component in the normal direction E 22 .
A significant improvement in the accuracy was achieved by decreasing the value of the Courant number. The best solution was obtained for C = 0.65 (Figure 9b). In this case, the propagating wave was symmetrical and smooth without any oscillations. However, a further successive reduction of the Courant number led to an increasingly less accurate solution. For example, with C equal to 0.2 ( Figure 9c) the oscillations were generated again but with a smaller magnitude than in the solution obtained with C = 1.5. Moreover, the oscillations appeared before the wave front, because in this case for, C < 0.71 (see Figure 4a), the value of the dispersion tensor component E 11 changed the sign to the opposite, i.e., E 11 takes a negative value. It is also worth noting that smaller Courant numbers correspond to a smaller time step ∆t. This untypical numerical property has its justification in the previously performed accuracy analysis.
Similar calculations were performed with α = 1/2 ( Figure 10). It can be noticed that in this case, the applied scheme also provided oscillatory results with the Courant number equal to 1.5 (Figure 10a). Only a decrease in the value of the Courant number leads to a better result with smaller oscillations (Figure 10b). Finally, an acceptable smooth solution without any oscillations is obtained for a very small value of C equal to 0.2 (Figure 10c).

Solution Using the 6NI Mesh Structure
The next calculations were carried out for the 6NI mesh structure, in which the orientation of the triangular elements is opposite to the previous ones. Figure 11 shows the results obtained using the scheme with the weighting parameter α = 1/3. As previously, the solution for C = 1.5 was characterized by large oscillations. The improvement in the accuracy was achieved by decreasing the values of the Courant number. However, in this case, numerical dispersion cannot be completely eliminated and, consequently, oscillations were observed even for small values of the number C.
The obtained solution confirms the previous analysis (see Figure 5a), which showed that for the applied scheme with α = 1/3, both the components E 11 and E 22 of the dispersivity tensor take non-zero values, leading to oscillations in both the tangential and perpendicular directions with regard to the velocity vector w.
In the case of α = 1/2, the same solutions were obtained as for the previous mesh structure (6N) presented in Figure 10, where reducing the Courant number led to an increasingly accurate solution without any oscillations.

Solution Using the 8N Mesh Structure
In the last test, simulations were performed for the 8N mesh structure. Figures 12 and 13 present the results of calculations with α = 1/3 and α = 1/2, respectively. As can be seen, the obtained solutions were dispersive, which was manifested by significant oscillations generated in the solution.
As it was for the previous simulation for the 6NI mesh with α = 1/3, these oscillations can neither be eliminated by adopting very small values of the time step ∆t nor by selecting the appropriate value of the weighting parameter α. It is especially visible in the solution for α = 1/2. This property was explained by analyzing the component E 22 of the dispersivity tensor (Equation (50)). In this case, for the accepted value of the parameters, the dispersion will always be generated in the perpendicular direction, even for C equal to zero. This is a different situation compared to the previous simulations with the use of the 6N and 6NI mesh structures, where the dispersion could be significantly eliminated by the appropriate choice of a value of the weighting parameter α, the Courant number or time step. Therefore, the 8N mesh structure, among all analyzed meshes, had the least favorable numerical properties.

Conclusions
In the paper, it was presented that a mesh structure with triangular elements has a significant influence on the accuracy of the solution obtained. The accuracy analysis performed using the modified equation approach allowed the diffusivity and dispersivity tensors to be derived for the applied scheme based on three different mesh structures. It seems that both numerical tensors facilitate an analysis of the numerical properties of the applied scheme. Since the derived tensors are a function of the numerical parameters (space and time interval, Courant number), then an analysis of the particular values of these parameters for the assumed mesh structure enables a detailed assessment of the numerical errors that may be generated in the solution. The conducted analyses and numerical simulations showed that if a time integration with θ > 0.5 was applied, then a significant smoothing and attenuation of the propagating wave in the direction of the flow were observed. For the fixed both mesh dimensions and velocity, the damping increased with the increasing value of the weighting parameter θ as well as the time step. Maximal diffusion was generated for θ = 1, which corresponded to the implicit Euler method. When the numerical diffusion was eliminated by using the trapezoidal formula with θ = 0.5, then the decisive factor for the property of the numerical solution will be the dispersivity tensor.
In such a case, an oscillating solution according to the form of the dispersivity tensor can be generated. The performed computations demonstrate that in this case, the modified FEM with α = 1/3 ensured a higher accuracy of the solution than the standard FEM with α = 1/2. Moreover, the accuracy analysis proved that the most accurate solution for a given flow direction was provided using the 6N mesh structure with the Courant number equal to approximately 0.65. On the other hand, the worst results were obtained using the 8N mesh structure. In this case, regardless of the value of the weighting parameter α and the values of the Courant number, non-physical oscillations were always generated in the solution.
From a practical point of view the presented study allows the following conclusions to be formulated:

•
The accuracy of the solution is influenced not only by the resolution of the mesh but also by its structure. Thus, by choosing an appropriate mesh structure, it is possible to minimize the numerical errors generated in the solution; • It is important to notice that maximum numerical diffusion (as well as numerical dispersion) occurs when the flow is oriented at a 45 o angle with respect to the axis of the coordinate system. Therefore, it is possible to minimize the numerical errors by an orienting the system so that the flow direction is parallel to the axis of this system; • It is commonly believed that by decreasing the time step value (or the Courant number value) a more accurate solution should be obtained. Unfortunately, such a general statement is not always true. As it was shown, while using the 6N mesh structure, the highest accuracy of the solution was obtained for the Courant number equal to 0.65 and successive decreasing the Courant number (also the time step) led to an increasingly less accurate solution. It is interesting that an explanation for this paradox was provided by the application of the modified equation approach, which is a robust tool in numerical analysis.
The accuracy analysis using the modified equation approach can be applied only to the linear equations, while the governing kinematic wave equation is a nonlinear. Despite this fact, the conclusions resulting from this analysis are of a more general nature and can also be applied to the numerical schemes used to solve the nonlinear problems. However, in some special cases, the discrepancies may arise between the results of the theoretical analysis and the simulations. Therefore, future studies should consider more real cases of the overland flow described by the 2D kinematic wave equation in its nonlinear form.