Adaptive method for the solution of 1D and 2D advection–diffusion equations used in environmental engineering

The paper concerns the numerical solution of one-dimensional (1D) and two-dimensional (2D) advection–diffusion equations. For the numerical solution of the 1D advection–diffusion equation, a method, originally proposed for the solution of the 1D pure advection equation, has been developed. A modified equation analysis carried out for the proposed method allowed increasing of the resulting solution accuracy and, consequently, to reduce the numerical dissipation and dispersion. This is achieved by proper choice of the involved weighting parameter being a function of the Courant number and the diffusive number. The method is adaptive because for uniform grid point and for uniform flow velocity, the weighting parameter takes a constant value, whereas for non-uniform grid and for varying flow velocity, its value varies in the region of solution. For the solution of the 2D transport equation, the dimensional decomposition in the form of Strang splitting technique is used. Consequently, this equation is reduced to a series of the 1D equations with regard to xand y-directions which next are solved using the aforementioned method. The results of computational experiments compared with the exact solutions confirmed that the proposed approaches ensure high solution accuracy of the transport equations.


INTRODUCTION
Consider the two-dimensional (2D) advection-diffusion transport equation with a source term: where t is the time, x and y are the space co-ordinates, c is the scalar function, U and V are the flow velocities in x-and ydirections, respectively, D x and D y are the coefficients of diffusion and w is the source term. Equation (1) is widely applied in many scientific and technique domains. It has a particularly important meaning in environmental engineering. The 2D equation describes most of the cases of transport of pollutants dissolved in water flowing in the open channels and in the shallow reservoirs. Note that the 3D advection-diffusion transport equation is applied just below the pollution discharge point, whereas when the contaminants are well mixed in the water over a channel reach then Equation (1) is usually simplified to the form of one-dimensional (1D) advection-diffusion transport equation: @c @t þ U @c @x À D x @ 2 c @x 2 ¼ f: Apart from the transport of contaminations, similar equation describes other flow phenomena as flood wave propagation (Gasiorowski 2014) and unsteady ground water flow towards a well (the Theis equation). In the case of pollutant transport, c(x, t) signifies a concentration, whereas in the case of unsteady flow, it means the discharge rate or the hydraulic head, respectively.
Usually, in environmental engineering, Equations (1) and (2) are called 'the advection-dispersion equation' (see for instance van Genuchten et al. 2013). However, in hydromechanics, the term 'dispersion' has a triple meaning and it is used in three different situations. The first meaning is related to the process of averaging the flow and transport parameters over the water depth or a channel cross-section, when the 2D or 1D equation is derived from the general form of the 3D transport equation. The second meaning is related to the transport of a scalar quantity in the ground waters, whereas the third meaning is related to the wave propagation process. As in the paper this term is used in the last sense, then to avoid any misunderstanding, the 2nd order terms in both transport equations will be considered as 'the diffusive terms'.
Both Equations (1) and (2) are the subject of a vast number of papers in which one can find both exact analytical and numerical solutions obtained for the appropriately formulated initial and boundary conditions. However, their numerical solution still attracts the attention of researchers because the provided solution very often is seriously affected by the numerical errors. The difficulties related to the numerical solution are caused by the properties of transport equations. Equations (1) and (2) represent the superposition of two different physical mechanisms of transport: advection having a hyperbolic character and diffusion having a parabolic character. The numerical problems arise when the transport is dominated by advection (Patankar 1980;Fletcher 1991). In such a case, the numerical solution of both Equations (1) and (2) becomes a challenging problem. For instance, if for the solution of Equation (2) the used method is of 2nd order, then to avoid spatial oscillations the so-called grid Péclet number must be not greater than 2 (Fletcher 1991): where Δx is the spatial mesh dimension. To obtain smooth solution of Equation (2) with P e . 2, more accurate methods should be used.
For the numerical solution of the 1D advection-diffusion transport Equation (2), a variety of numerical methods and techniques have been proposed. They use the different schemes of the finite difference method (FDM), the finite element method (FEM) or the finite volume method (FVM). A comprehensive discussion on the mentioned methods is presented, for instance, by Cunge Holly & Verwey (1980), Fletcher (1991, Gresho & Sani (1998), Quarteroni Sacco & Saleri (2000, LeVeque (2002) and Morton & Mayers (2005), as well as by many other authors of books and papers. For instance, all standard methods as the Lax-Wendroff scheme, the Crank-Nicolson (C-N) FDM and linear FEM scheme and the three-level fully implicit scheme represent the 2nd order of accuracy (Fletcher 1991). Consequently, they do not generate the numerical diffusion. For this reason, the solution without spatial oscillations is possible only on condition that the relation (3) is satisfied. Therefore, the implementation of aforementioned methods for advection-dominated transport requires decreasing of the space interval Δx, usually to an unacceptable small value. For this reason, the methods of higher accuracy are more interesting. As an example, the QUICKEST scheme (Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms) can be considered. It approximates the advection-diffusion equation with the accuracy of 3rd order involving four nodes from the preceding time level (Leonard 1979). The scheme does not generate numerical diffusion, and it does not produce the spatial oscillations for P e 8/3. Another idea frequently applied for the numerical solution of the advection-diffusion equation bases on the splitting process. The concept of splitting technique (or decomposition method) is to split in each time step the governing equation into a sequence of simpler equations. Such approach saves the computational time and memory. The splitting process can be carried out in different ways. Theoretical background of this technique was given by Yanenko (1971). Comprehensive and up-todate information on the splitting technique for hyperbolic problems is given by LeVeque (2002). According to this technique, the 1D advection-diffusion transport equation is split into advection equation and diffusion equation in each time step. Next, both equations are solved separately. For instance, Cunge et al. (1980) suggest the approach proposed by Holly & Preissmann (1977) for the solution of the advective part and an explicit difference scheme for the diffusive part with the initial condition being the solution obtained in the previous stage. Other methods can also provide satisfying results (Szymkiewicz 2010). However, it should be added that the splitting process introduces additional error. In the case of 2D advection-diffusion transport equation, usually dimensional splitting is used. Consequently, one obtains two series of 1D convection-diffusion problems with regard to both directions x and y, which are easier to solve numerically.
There are also known other methods of higher accuracy. Such methods were developed via the approximation of the highorder derivatives occurring in the terms representing the error of approximation. However, if higher order terms are included, then additional nodes required for the approximation of the derivatives must be involved. As examples of such an approach, one can mention the methods proposed, for instance, by Kim (2003) as well as by Man & Tsai (2008), Appadu et al. (2016) and by many others. For increasing the solution accuracy typically, a two-level stencil is developed in time, in space or simultaneously, in time and space. Instead of simple grid points allowing the standard approximation of space and time derivatives, more complex grid points and, consequently, more complex algorithm, decreasing the effectiveness of computations, must be used. Tian & Ge (2007) proposed a method which can be considered as a combination of both above-mentioned approaches. For the solution of the unsteady 2D convection-diffusion equation, the alternating direction implicit (ADI) method is used. The ADI method is based on splitting in time carried out in each time step. This process leads to 1D convection-diffusion problems which are solved by Tian & Ge (2007) using the 2nd order formula for the time discretization and the 4th order difference formula for the spatial discretization. The method uses a regular five-point stencil typical for 2D problems, and it requires solving the systems of algebraic equations with tri-diagonal matrices.
It seems that it is possible to develop a numerical method for the solution of the 1D advection-diffusion equation which is free from previously mentioned disadvantages. The method proposed here is an unsplit one. It uses the standard two-level stencil involving only three nodes at each level, whereas increasing of solution accuracy is reached via vanishing the terms containing the higher-order derivatives occurring in the error of approximation. To this order, the modified FEM with the linear shape functions for the spatial approximation and the implicit trapezoidal rule (TR) for the time integration are used. This approach, originally applied for the 1D pure advection equation and for the systems of 1D hyperbolic equations describing an unsteady open channel flow (Szymkiewicz 1995), in this paper is developed for the 1D advective-diffusive transport equation. The developed algorithm possesses an important feature. Although both applied methods, i.e. FEM and TR, represent an approximation accuracy of the 2nd order, the final method's accuracy can be increased up to the 4th order with regard to both the independent variables x and t. This could be achieved owing to the modified equation analysis proposed by Warming & Hyett (1974). This analysis provides explicit formulas for the coefficients of numerical dissipation and dispersion in which the weighting parameter is involved. Vanishing the derived coefficients via a proper choice of their value increases the solution accuracy. It is important that the value of weighting parameter is able to adopt. Because it depends on the Courant and diffusive numbers, it is a function of flow velocity and, consequently, it varies for non-uniform grid points or non-uniform flow.
As far as the 2D advection-diffusion transport Equation (1) is considered, it can be solved using the unsplit techniques as FDM, FEM or FVM. In such approach, Equation (1) is directly approximated in a way typical for the used method of solution (see for instance Fletcher 1991; Gresho & Sani 1998;Quarteroni Sacco & Saleri 2000). However, such an approach usually leads to the time-consuming algorithm with large systems of algebraic equations. To avoid this disadvantage, the aforementioned splitting technique can be applied also for the solution of Equation (1). This equation can be transformed into simpler equations via splitting with regard to the dimensions or to the involved physical processes (Fletcher 1991;LeVeque 2002). The computational advantage of such an approach is obvious. The solution of a series of 1D problems is less costly, particularly when the parallel computing is applied (Artichowicz & Gasiorowski 2019).
One of the most accurate and effective variants of decomposition seems to be the method proposed by Strang (1968) belonging to the family of operator splitting techniques. The Strang method is 2nd order accurate with respect to the unsplit or exact solution. Consequently, this method does not introduce or minimizes the so-called splitting error if the governing equation is non-linear (LeVeque 2002). For this reason, the Strang method is used in the paper to perform the dimensional splitting of Equation (1).
The applied numerical methods should ensure the physically justified shape of the propagating wavefront. This means that the solution should not be artificially smoothed and should be oscillation-free. Both effects are caused by the errors of numerical diffusion and numerical dispersion generated in the numerical solution. These issues concern 1D problems as well as 2D and 3D ones. Since water resources management requires reliable computational tools not only for the analysis of the available water quantity but also for the analysis of its quality, then effective and accurate numerical methods, which are able to ensure that the results are not affected by the aforementioned spurious effects, are still sought. This article attempts to develop the numerical methods for solving the 1D and 2D advection-diffusion equations that are effective and provide the solutions with high accuracy.

Numerical solution of the 1D advection-diffusion equation
Consider the 1D advection-diffusion transport equation (2) with a variable flow velocity (U ¼ U(x)) and a constant coefficient of diffusion (D x, ¼ const) which will be solved for 0 x L x and t ! 0. Assuming that a non-uniform grid constituted by M nodes covers a space domain of length L x being the length of considered open-channel reach, Equation (2) will be solved numerically using the modified FEM. For its solution, the modified Galerkin FEM is applied. The considered channel reach of length L x is divided into M À 1 elements of length Δx j ( j ¼ 1, 2, …, M À 1) (Figure 1).
According to the Galerkin procedure, the numerical solution of Equation (2) must satisfy the following condition (Zienkiewicz 1972): T is the vector of linear shape functions, j is the index of node, M is the total number of grid points and c a is the aproximation of function c(x,t) defined as follows: Certain terms in formula (4) can be approximated using another approach: where the subscript c denotes the approximation of function according to the modified formula. Let us consider the first term of the above sum, corresponding to the element j which is bounded by the nodes j and j þ 1 (Figure 1). Since in the element j only two components of the vector N(x), i.e. N j (x) and N jþ1 (x), are non-zero, then in Equation (6), only two following non-zero products will occur: where c a is the approximation of function c corresponding to the standard approach (5): where c j represents the nodal value of function c(x, t). Integration of the first terms of Equations (7a and 7b) yields: Detailed information on the FEM is provided, e.g. by Zienkiewicz (1972), Fletcher (1991 and Gresho & Sani (1998). More general formulas can be used instead of standard approximations (Equations (10a) and 10(b)): where ω j is a weighting parameter ranging from 0 to 1 (Szymkiewicz 2010). It is important that this parameter is related to the element j but not to the node j. Note also that for ω j ¼ 2/3, Equations (11a) and 11(b) corresponding to Equations (10a) and 10(b) are given by the standard FEM. The calculation of the integrals (7a) and (7b) in any single element is carried out step by step assuming that the weighting parameter ω varies in the consecutive elements. The equations obtained for all elements ( j ¼ 1, 2, …, M À 1) are assembled according to Equation (4) giving the system of ordinary differential equations: The initial-value problem for the system (12) is solved using the implicit TR ( For further analysis, assume that a uniform grid with equally spaced nodes (Δx j ¼ Δx ¼ const) covers considered open-channel reach of length L x in which flow velocity is constant (U ¼ const) and the source term ϕ is neglected. In such a case, the solution of Equation (12) using the implicit TR yields the following system of algebraic equations: where c n j is the nodal value of function c at the time level n, ω is a weighting parameter considered as constant. In each time step, the system (13), completed by the imposed boundary conditions, is solved using the Thomas algorithm (Fletcher 1991;Quarteroni Sacco & Saleri 2000) dedicated to the systems with tri-diagonal matrices.
The involved weighting parameter ω determines the type of approximating formula applied. For ω ¼ 1, the formula corresponding to the C-N FDM is obtained, whereas for ω ¼ 2/3, Equation (13) coincides with the C-N FEM with linear basis functions. Consider a single rearranged equation involved in the system (13), corresponding to any internal node j ¼ 2, 3, …, M À 1: where are the Courant number and the diffusive number, respectively. A stability analysis carried out using the Neumann method shows that the numerical method represented by Equation (14) is absolutely stable on condition that the weighting parameter ω ranges from the interval 〈0.50, 1.0〉 (Szymkiewicz 1995). As far as the solution accuracy is considered, it can be estimated via the modified equation analysis (Warming & Hyett 1974;Morton & Mayers 2005). To this order, let us replace in Equation (14) all nodal values of the function c by their Taylor series expansions around the point in which the solved equation was approximated. After appropriate arrangement at each node ( j, n), one obtains the following modified advection-diffusion equation: where E n2 , E n3 , E n4 and E n5 are the coefficients of numerical dissipation if they are related to the even order derivatives, or the coefficients of numerical dispersion if they are related to the odd order derivatives. They are given by the following formulas: The derived Equation (17) is the advection-diffusion equation (2) modified by the used numerical method. Its left-hand side is consistent with the left-hand side of the considered equation (2), whereas its right-hand side contains all terms of the Taylor series expansion neglected during the approximation of the derivatives. Since at the right-hand side of Equation (17) the term of the 2nd order does not exist (E n2 ¼ 0), then the method does not produce numerical diffusion.
The accuracy of the proposed method can be increased by eliminating also the second term occurring at the right-hand side of the modified equation (17) related to the derivative of the 3rd order. In this case, for the following value of ω: the term with the derivative of the 3rd order in Equation (14) is cancelled because E n3 ¼ 0. However, for this value of ω, the term with the 4th order derivative exists because E n4 ≠ 0. Therefore, with Equation (22), the order of approximation is increased up to the 3rd one with regard to both variables t and x. For other values of ω, the method's accuracy is of the 2nd order because E n3 ≠ 0. When C dx ¼ 0, the value of weighting parameter (22) corresponds to its value derived for the pure advection transport equation (Szymkiewicz 1995). Note that for the pure advection equation (C dx ¼ 0), the parameter ω given by Equation (22) cancels also the coefficient E n4 ensuring the accuracy of the 4th order. Moreover, for the Courant number equal to 1, all terms at the right-hand side are vanished, which means that the method provides the exact solution of the pure advection equation.
In view of numerical stability, the assumed value of the weighting parameter ω has to belong to the interval 〈1/2, 1.0〉. According to Equation (22), this parameter depends on the Courant number C ax and the diffusive number C dx . Therefore, this condition can be satisfied for various combinations of C ax and C dx , including also a value of C ax greater than 1. In Figure 2, a map of the function ω ¼ f(C ax , C dx ) is shown. Any value of the weighting parameter ω belonging to the shadowed area provides a numerically stable solution.
If in Equation (12) the source term ϕ is taken into account, then instead of Equation (13) the following system of algebraic equations should be used: Journal of Hydroinformatics Vol 23 No 6, 1296 Decomposition of the 2D advection-diffusion transport equation using the Strang method Consider the 2D advection-diffusion transport equation with a source term (1). Assume that this equation will be solved using the Strang splitting method (Strang 1968). This method is used to perform the dimensional splitting.
Equation (1) can be rewritten in the following form: where In its standard form, the Strang splitting process consists of the three following steps: • Ist steptime integration in x-direction over half of the time step Δt/2: • IInd steptime integration in y-direction over the time step Δt: ω belonging to the shadowed area ensures numerical stability.
Journal of Hydroinformatics Vol 23 No 6, 1297 • IIIrd steptime integration in x-direction over half of the time step Δt/2: The solution obtained in the IIIrd step is the sought solution of the governing equation at the next time level, i.e.
The idea of the Strang splitting technique can be easily interpreted using the properties of integral. To this order, let us integrate Equation (1) over the time step Δt: The resulting relation: can be rewritten as follows: Substituting one obtains: Next, substitution leads to the following formula: The formulas (34), (36) and (37) correspond to the 1D differential ordinary equations (28), (29) and (30), respectively. Assume that Equation (1) is solved numerically in the rectangular region of the dimensions L x Â L y with the appropriate auxiliary conditions. The considered region is covered with a regular mesh having N rows and M columns (Figure 3).

RESULTS AND DISCUSSION
where erfc(·) is the complementary function to the error function erf(·) (McQuarrie 2003). Equation (38) illustrates the propagation of a discontinuity in space and time with a constant advection speed which simultaneously is subjected to physical diffusion ( Figure 4). Assume also the following arbitrarily chosen set of data: constant flow velocity U ¼ 0.50 m/s, space interval Δx ¼ 0.50 m, time step Δt ¼ 0.75 s and coefficient of diffusion D x ¼ 0.0075 m 2 /s. The integration is carried out consecutively up to t ¼ 60 and 120 s. Since the discontinuity travelling with a constant speed is simultaneously subjected to diffusion, then the initially steep front is systematically smoothed as time goes by. For the assumed data, the following values of the numerical parameters are equal to C ax ¼ 0.75, C dx ¼ 0.00225, P e ¼ 33.33 and ω ¼ 0.595 (given by Equation (22)). The analytical and numerical solutions are compared in Figure 4. As one can see, an increase in the solution accuracy ensures a satisfying agreement between the exact and numerical solutions. Parameter ω chosen according to Equation (22) causes that E n3 ¼ 0. This means that in this case, the numerical method represents the 3rd order of accuracy. Additionally, this value of ω decreases the next term E n4 at the right-hand side of Equation (17) associated with the derivative of the 4th order. Consequently, the increased solution accuracy allows a very good agreement to be obtained for a relatively high value of the Péclet number. For comparison, the same problem was solved using the finite difference C-N method. The C-N method, commonly used for the solution of the pure diffusion equation, is useless for the pure advection equation. Note also that the C-N method corresponds to the proposed method with ω ¼ 1. The results of solution obtained with exactly the same set of data as previously are shown in Figure 5.
The comparison of the graphs presented in Figures 4 and 5(a) shows that the C-N method is remarkably less accurate than the proposed one. The spurious oscillations affecting the travelling steep front of concentration are still present. It is well known that for the convergent numerical method, the accuracy of the provided solution can be increased via the reduction of the used mesh dimensions. In the next test, let us decrease the time step assuming a 10 times smaller value than previously (Δx ¼ 0.50 m and Δt ¼ 0.075 s). The results of computations presented in Figure 5(b) are very similar to the ones displayed in Figure 5(a) so that any improvement of solution accuracy is not observed. However, if simultaneously apart from the time step also the space interval is decreased, the solution accuracy increases. This effect is observed in Figure 5(c) where the results obtained for Δx ¼ 0.050 m and Δt ¼ 0.075 s are displayed. As one can see, in this case the perfect agreement between the exact and numerical solutions is obtained. Unfortunately, this approach increases remarkably both the dimension of the final system of algebraic equations and the total number of the time steps. Consequently, it can be applied while the theoretical consideration but it rather cannot be used in engineering practice.⧫ One can find out that there is a remarkable difference between the proposed approach and other well-known methods, which usually provides satisfying results only for a smaller value of the Péclet number and rather for smooth functions. For instance, the method of solution of the advection-diffusion equation proposed by Nazir et al. (2016), in which the cubic trigonometric B-spline functions are used for interpolation, provides satisfying results for the Péclet number not greater than 5 (P e 5), whereas using another kind of the B-spline functions and the collocation method Mittal & Jain (2012) reached the satisfying results for P e 10. It is worth adding that the very famous QUICKEST scheme requires P e 8/3. Note that in Example 1, we obtained very good agreement with the exact solution for P e ¼ 33 although in the imposed upstream boundary condition a discontinuity was present. It seems that this makes a great difference between the methods.

Example 2
The advection-diffusion equation (2) without the source term is solved for the initial condition in the form of a Gaussian distribution ( Figure 6, t ¼ 0):  where m p is the mass of the pollutant introduced into the water, whereas σ is the standard deviation. At the ends of the considered channel's reach, the following boundary conditions: c(0, t) ¼ 0 and c(25, t) ¼ 0 are imposed, respectively. For these initial and boundary conditions, the advection-diffusion equation has the following exact solution (Massel 2010): Assume that the mass of the introduced pollutant is equal to m p ¼ 1 kg/s, the flow velocity is equal to U ¼ 1 m/s and the coefficient of diffusion is equal to D x ¼ 0.020 m 2 /s, whereas the standard deviation is σ ¼ 0.25 m. The analytical solution (40) corresponding to the chosen values of parameters calculated for the selected times (t ¼ 5, 10, 15 s) is represented in Figure 6 by the solid lines.
As the initial distribution of the pollutant concentration is simultaneously subjected to advection and diffusion, then while it travels along the channel axis with a constant speed U ¼ 1 m/s, it is systematically smoothed as time goes by. In the same Figure 6, the analytical solution is compared with the numerical one provided by the proposed method. In this case, the following values of the numerical parameters are assumed: the space interval Δx ¼ 0.10 m and the time step Δt ¼ 0.05 s. For these mesh dimension values, the Courant number is equal to C ax ¼ 0.50, the diffusive number is equal to C dx ¼ 0.10 and the Péclet number is P e ¼ 5.0, whereas the weighting parameter given by Equation (22) is ω ¼ 0.725. The results of computations represented in Figure 6 by the dotted lines confirm the satisfying accuracy of the proposed method.
The obvious advantage of the presented method is illustrated by the extreme case of the transport equation (2), i.e. when diffusion does not exist. For D x ¼ 0, Equation (40) becomes the analytical solution of the pure advection equation: which describes a translation of the initial distribution of the pollutant c(x,0) along the channel axis without its shape deformation. This fact is confirmed in Figure 7, where the analytical solution (solid lines) is compared with the numerical solution (dotted lines). The presented results were obtained for the same set of data as assumed previously except for the value of the diffusion coefficient, which was equal to D x ¼ 0. The centre of gravity of the initial distribution being initially in the position x ¼ 0, after t ¼ 5, 10 and 15 s occurs at the positions: x ¼ 5, 10 and 15 m, respectively (Figure 7). Since E n2 ¼ 0, then numerical diffusion is not generated in the solution. Parameter ω chosen according to Equation (22) with D x ¼ 0 means that C ax 1.0 results in E n3 ¼ 0 and E n4 ¼ 0. As far as the coefficient of numerical dispersion E n5 is concerned, its value is different from zero for 0 , C ax , 1.0. In spite of this, the numerical solution of the pure advection equation obtained for C ax ¼ 0.5 coincides with the analytical one ( Figure 7). Bear in mind that the Péclet number corresponding to the pure advection equation is P e ¼ ∞.⧫ Note that for the advection equation, the value of parameter ω chosen according to Equation (22) and C ax ¼ 1.0 causes all terms at the right side of the modified equation (14) to disappear. In such a case, all coefficients E ni related to the consecutive derivatives become equal to zero so that the proposed method ensures the exact solution. However, the results presented in Figure 7, obtained for C ax ¼ 0.5, are also in excellent agreement with the exact ones. This fact is obvious because the approximation of the derivatives involved in the advection transport equation gives better results for a smooth function than the approximation of the derivatives carried out for a function containing discontinuity. The approximation of the derivatives is based on the Taylor series expansion, valid for smooth functions.

Example 3
In the semi-infinite open channel, initially clear water (c(x, 0) ¼ 0) flows with a constant velocity U. Consider the transport of a pollutant introduced at the upstream end (x ¼ 0) in the form of the following pulse (Figure 8): Assume also that the coefficient of diffusion D x is constant, whereas the source term takes the form ϕ ¼ Àk·c (the first-order reaction) with the constant coefficient of decay k (T À1 ).
For the assumed data, the analytical solution of Equation (2) has the following form (Mazaheri et al. 2013): The following data are arbitrarily assumed: flow velocity U ¼ 1 m/s, space interval Δx ¼ 0.50 m, time step Δt ¼ 0.20 s and coefficient of diffusion D x ¼ 0.020 m 2 /s. The imposed pulse (Figure 8) is defined by the following value of its parameters: t 1 ¼ 5 s and t 2 ¼ 20 s. Travelling along the channel's axis, it is simultaneously subjected to both diffusion and a reaction characterized by the decay parameter equal to k ¼ 0.0025 s À1 . Consequently, the pulse height is continuously reduced and the discontinuities present at both its sides are systematically smoothed. The solutions corresponding to time t ¼ 45 s, obtained analytically and numerically with ω given by Equation (22), are compared in Figure 9.
For the assumed data, the value of the Péclet number is equal to P e ¼ 25, which means that in this case, the transport is dominated by advection. In spite of this, both solutions agree perfectly. The value of the weighting parameter ω ¼ 0.656, chosen according to Equation (22), ensures the 3rd order of accuracy. Moreover, this value of ω reduces the value of the coefficient E n4 associated with the term containing the derivative of the 4th order present at the right-hand side of Equation (17). Consequently, a very good accuracy of the numerical solution was obtained. This result was achieved despite assuming a relatively large value of the Péclet number.⧫ This example deals with the Theis equation describing the radial groundwater flow in a confined aquifer towards a well. Assuming one directional flow uniformly distributed over depth in infinite horizontal layer with constant thickness, Theis (1935) proposed the following equation: where t is the time, r is the space co-ordinates taken from the axis of well, h is the hydraulic head, U is the flow velocity in r directions, S is the storage coefficient, T ¼ k·b is the transmissivity of aquifer, k is the hydraulic conductivity and b is the thickness of layer. With the following initial and boundary conditions: h(r, 0) ¼ h 0 for 0 r R, h(R, t) ¼ h 0 for t . 0 and for the imposed well discharge rate Q, the analytical solution of Equation (43) expressed in terms of drawdown s ¼ h 0 À h(r, t) is as follows: The so-called well function W(u), where u ¼ r 2 ·S/(4T·t), can be approximated by the following formulas (McElwee 1980): • W(u) ¼ Àln(u) À 0.57721566 þ 0.99999193u À 0.24991055u 2 þ 0.05519968u 3 þ 0.00976004u 4 þ 0.00107857u 5 for 0 u 1.
The Theis equation (43) can be solved also numerically using the proposed method. To this order, it is rearranged to the standard form of advection-diffusion equation: in which U ¼ ÀT/(r·S) and D r ¼ T/S, where D r is the coefficient of diffusion. Computations were carried out for the following arbitral set of data: thickness of ground layer is b ¼ 15 m, its total length layer is R ¼ 100 m, hydraulic conductivity is k ¼ 0.0001 m/s, initial head is h 0 ¼ 30 m, storage coefficient is S ¼ 0.002 and well discharge rate is Q ¼ 30 m 3 /h ¼ 0.0083 m 3 /s. The space interval was assumed to be Δr ¼ 1.0 m, while the time step was equal to Δt ¼ 10.0 s. The obtained results are compared in Figure 10. Note that in the considered case, the flow velocity U varies from infinity at the well axis (r ¼ 0) to zero for r→∞. Consequently, the weighting parameter ω, which according to Equation (22) depends on the Courant number (C ax ¼ U·Δt/Δr), also varies in time and in space. As it is related to the consecutive finite elements (see Figure 1), then for element j the Courant number is evaluated using the average flow velocity taken from the nodes j and j þ 1: This means that the value of ω varies from node to node. Therefore, the equation of type (12) is applied. Summarizing, one can find out that the proposed method is adaptive one because while computing the value of weighting parameter ω (r, t) varies for variable flow velocity and consequently, the same higher order of method accuracy is still kept.
The hydraulic head being the analytical solution of Theis equation at time t ¼ 1,200 s is represented by the continuous line, whereas its numerical solution provided by the proposed method is represented by the dotted line. As one can see, both curves are very similar to one another. Because for r ¼ 0 a singularity appears in Equation (45), the computations must start from r ¼ r 0 where r 0 is a radius of well assumed to be equal to 0.30 m.
It is interesting that if ω is taken as for the pure advection equation (Equation (22) with C dx ¼ 0), the results (cross line in Figure 10) are remarkably less accurate than obtained previously. This fact confirms that the weighting parameter ω has important meaning. Its values resulting from the modified advection-diffusion equation analysis ensure more accurate solution than ones derived for the pure advection equation.

Impact of mesh dimensions on the solution error
The order of the numerical method used to solve the partial differential equation can be inferred by analyzing the behaviour of the solution error when the grid is refined. To this end, reconsider Equation (2) without source term (ϕ ¼ 0). At first, assume that it is solved for the same auxiliary conditions as in Example 1 (i.e. c(x,0) ¼ 0 for x ! 0 and c(0, t) ¼ 1.0 for t ! 0) so that Equation (38) is its exact solution. The aim of numerical tests is to check the effect caused by successively refined mesh. It can be evaluated by comparison of the numerical solution with the exact one. To this order, the following definition of solution accuracy is used: where L is the length of considered channel reach, T is the assumed limited time, c(x,T ) is the exact solution of Equation (2) given by Equation (38) and c(x, T) is the numerical solution, both corresponding to t ¼ T. Starting with the assumed initial values of Δx and Δt, one obtains the corresponding value of the Courant number C ax and the diffusive number C dx involved in Equation (2). The grid refinement should be carried out in such a way that constant values of both dimensionless numbers are kept. However, according to Equations (15) and (16), such a value of time step Δt, which simultaneously ensures the constant value of both C ax and C dx , does not exist. Assuming the value of space interval Δx, the selected value of time step Δt can ensure constant value either the Courant number or the diffusive number. Consequently, it is impossible to show the effect of mesh refinement using the standard approach for the numerical examination of the order of accuracy. It seems that the evaluation of the order of proposed method can be carried out indirectly. To this order, one can analyse separately both advective and diffusive parts of Equation (2).

Example 5
Putting U ¼ 0, Equation (2) is reduced to the diffusion equation. Note that its exact solution for imposed initial and boundary conditions is given also by Equation (38) with U ¼ 0. Diffusion equation is numerically solved over an interval 〈0.0, 100.0 m〉, whereas the time integration is carried out up to t ¼ T ¼ 120 s. Assume also that the coefficient of diffusion D x ¼ 0.10 m 2 /s is constant, whereas the diffusive number is C dx ¼ D x ·Δt/Δx 2 ¼ 0.25. Since U ¼ 0, the Courant number C ax is also equal to zero and formula (22) is reduced to the following one: ω ¼ 2/3 þ C dx . The results of calculations carried out for systematically reduced values of Δx and for the corresponding values of Δt are presented in Table 1.
As one can see, all variants of the method, including the proposed one, represent accuracy of the 2nd order because the δ error changes with the square of Δx. Moreover, all the proposed variants produce the solution error having similar value. This fact is understandable because all three considered variants of the method are the methods of C-N.⧫

Example 6
The second particular case of the advection-diffusion equation (2) is the pure advection equation. When diffusion does not exist (D x ¼ 0), the diffusive number is also equal to zero (C dx ¼ 0), whereas the Péclet number is equal to infinity (P e ¼ ∞). The advection equation (2) without the source term is solved for the initial condition in the form of a Gaussian distribution and for the same boundary conditions as used in Example 2. Its exact solution is well knownthe initial distribution of the pollutant c(x,0) travels along the channel axis without shape deformation. Note that for C ax ¼ 1.0, the method provides exact solution.
The flow with velocity equal to U ¼ 1 m/s takes place in a channel reach of length L ¼ 25 m, whereas the mass of the introduced pollutant is m p ¼ 1 kg/s. The standard deviation equal to σ ¼ 0.50 m causes that the initial distribution of the pollutant is relatively steep because it is spread over approximately 4.0 m. The equation is integrated up to t ¼ T ¼ 15 s. The computations were carried out for two values of space interval Δx, whereas the value of time step Δt corresponds to assumed values of the Courant number equal to C ax ¼ 0.25 and 0.50. The results of computations are displayed in Table 2. As one can see, the relation between the errors corresponding to both assumed values of the Courant number shows that doubling of Δx causes an approximately 16-fold increase in the solution error δ. This suggests that the method represents 4th order accuracy.
As initial distribution of the pollutant in the form of a Gaussian distribution occurs over short distance of channel (about 4.0 m only), then larger Δx values were not used.
The modified equation (17) confirms this conclusion. For both used values of the Courant number, the term with derivative of the 4th order disappears because E n4 ¼ 0 (Equation (20)). Moreover, both used values of ω decrease the value of next term at the right-hand side of Equation (17).⧫  Unfortunately, a similar analysis as performed above for both diffusive and advection equations separately cannot be carried out for the advection-diffusion equation. This is because in the approximating algebraic equations, both the Courant number and the diffusive number are involved. Consequently, the time step Δt chosen for assumed space interval Δx ensures that only one mentioned number, C ax or C dx , can preserve its constant value during mesh refinement. Moreover, increasing of the space interval Δx increases the Péclet number which is not related to the time step Δt. This means that spurious effects will be generated in the solution. Nevertheless, our numerical experiments showed that the proposed method ensures very good agreement with exact solution even for a relatively high value of the Péclet number.

Examples of the numerical solution of the 2D advection-diffusion equation
Example 7 The transport of matter dissolved in water takes place in a large shallow reservoir. The horizontal dimensions of the considered part of the reservoir are assumed to be equal to L x ¼ 100 m and L y ¼ 100 m (Figure 11). The region of the solution is covered with a uniform grid point having constant dimensions Δx ¼ Δy ¼ 1.0 m. Consequently, M ¼ 101 and N ¼ 101. Both components U and V of the nodal flow velocity vector are assumed to be constant and equal to 0.5 m/s. The imposed initial condition has the form of the 2D Gaussian distribution with the standard deviation σ ¼ 4 m. Its extreme concentration equal to c max ¼ 1.00 occurs at the point (x g ¼ 20 m, y g ¼ 20 m) (see Figure 11, t ¼ 0). The assumed data ensures that the initial distribution of the concentration (t ¼ 0) travels along the diagonal of considered region.
The first test deals with the solution of the pure advection equation (in Equation (1), D x ¼ D y ¼ 0 and w ¼ 0). The results obtained for Δt ¼ 1.0 s are displayed in Figure 11.
Note that in this case, the exact solution of the 2D advection equation is known. The concentration distribution imposed initially does not change its shape as time goes on. Therefore, the extreme value of the 2D Gaussian distribution imposed as the initial condition occurred after t ¼ 60 s at the point (x ¼ 50 m, y ¼ 50 m), whereas after t ¼ 120 s, it reached the point (x ¼ 80 m, y ¼ 80 m). It is very important that the extreme values of the Gaussian distribution after t ¼ 60 s and t ¼ 120 s are equal to c max ¼ 1.00, which means that any numerical dissipation does not exist. Moreover, the computed distributions of the concentration are not affected by any spurious oscillations. Note that the Courant numbers were equal to  whereas the weighting parameter taken according to Equation (22) was constant and equal to ω ¼ 0.625. In the next test, the non-uniform field of flow velocity was assumed. Water rotating around the centre of the reservoir (x ¼ 50 m and y ¼ 50 m) with an angular velocity equal to f ¼ 0.01 1/s requires T ¼ 2π/f ¼ 628 s to carry out one cycle. The imposed initial condition has the same form as the 2D Gaussian distribution but its extreme concentration equal to c max-¼ 1.00 occurs at the point (x g ¼ 20 m, y g ¼ 50 m) (see Figure 12, t ¼ 0). In this case, the initial distribution of the concentration (t ¼ 0) travels around the axis of rotation, and after t ¼ T ¼ 628 s, it returns to its starting position. Both components U and V of the nodal linear flow velocity vector, corresponding to the assumed angular velocity f, vary from À0.50 to 0.50 m/s. At first, the pure advection equation (Equation (1) with D x ¼ D y ¼ 0 and w ¼ 0) was solved. The results obtained for Δt ¼ 0.5 s are displayed in Figure 12. Note that the concentration distribution imposed initially does not change its shape during rotation. Therefore, after t ¼ T ¼ 628 s, the assumed 2D Gaussian distribution is expected exactly at its initial position corresponding to t ¼ 0. As one can see, the Strang decomposition process assembled with the proposed method provides a very good final effect. Although the horizontal rotation of the Gaussian distribution is a very challenging issue for numerical solution methods, the results of computation displayed in Figure 12 represent very high accuracy. After one cycle (t ¼ T ), the numerically computed concentration is practically the same as its initial distribution assumed at t ¼ 0. This statement is proved by the following facts.
After t ¼ T ¼ 628 s, the extreme value of the concentration was reduced from 1.000 to 0.999, whereas the total mass of the dissolved matter dispersed in the water, estimated using the formula: which initially was equal to Total in ¼ 100.531, after t ¼ T is also equal to Total T ¼ 100.531. This means that any serious mass balance error was not generated. Moreover, the isolines displayed in Figure 12 confirm that the used method neither smoothes the provided solution nor produces spurious oscillations. It is important that this perfect result has been obtained for the spatially variable vector of flow velocity. The advective Courant numbers with regard to x-and y-directions, defined by Equations (48) and (49), take a value ranging from 0 to 0.25 over the considered area. Consequently, both variable Courant numbers cause the weighting parameter ω to vary as well in space but not in time. Therefore, depending on the solved equation, one has ω ¼ ω(x) or ω ¼ ω(y) which are calculated according to Equation (46). Owing to this adaptive property, although the flow velocities vary, the method still maintains the 4th order of accuracy. It should be added that this is valid as long as the greatest values of both Courant numbers do not exceed unity (C ax 1 and C ay 1). The computations confirmed that, indeed, similarly satisfying results are provided for a greater value of Δt than used previously, whereas spurious effects occur in the solution for C ax . 1 or C ay . 1 when Equation (46), determining the value of the weighting parameter ω, is no longer valid.
In the next tests, the 2D advection-diffusion equation and the advection-diffusion equation with a source term are considered. In both cases, the coefficient of diffusion D x ¼ D y ¼ 0.010 m 2 /s is assumed, whereas the decay parameter is assumed to be equal to k ¼ 0.0005 1/s. The computations were carried out for the same data as used previously. The results are presented in Figure 13 in the form of concentration distributions along a row corresponding to y ¼ 50 m (see Figure 12) and t ¼ T ¼ 628 s. For comparison, an additional concentration distribution, being the solution of the previously considered 2D pure advection equation, is also displayed (Figure 13curve a).
As one can see, also for the advection-diffusion equation, the computed results are not affected by any spurious effects. The balance of total matter dissolved in water is satisfied because after t ¼ 628 s the integral (49) is equal to Total T ¼ 100.516, while initially it was equal to Total in ¼ 100.531. In this case, the diffusive numbers, take values equal to C dx ¼ C dy ¼ 0.005. The variable weighting parameter ω, given by Equation (22), ensures an oscillation-free solution also in this case ( Figure 13curve b). As far as the solution of the advection-diffusion transport equation with a source term is considered, it is illustrated in Figure 13 by curve c. From the mass balance of matter dissolved in water, it results that after t ¼ T ¼ 628 s, the integral (50) is equal to Total T ¼ 62.760, while initially it was equal to Total in ¼ 100.531. The mass of dissolved matter, which disappeared during rotation, depends on the assumed value of the decay coefficient k. Consequently, compared with the solution of the advection-diffusion equation, the extreme concentration is decreased, whereas the borders of spatial distribution of the concentration are the same.⧫ In the paper, it was shown that the proposed method, which bases on the modified FEM and the implicit TR for time integration, is a very effective tool for the numerical solution of the 1D advection-diffusion transport equation. It works on the standard two-level stencil. The modified equation analysis allowed deriving explicit formulas for the numerical dissipation and dispersion. As they are the functions of the involved weighting parameter ω, its proper choice vanishes some terms of approximation error increasing the method's accuracy. This method possesses a very important feature. It is adaptive because for the varying flow velocities, its weighting parameter also varies in space and in time ω ¼ ω (x,t) (or ω ¼ ω (y,t)). It was shown that the method ensures very good agreement between the numerical and analytical solutions of the 1D advection-diffusion equation with and without a source term. It seems that in comparison with other known methods, the method of solution presented in this paper is relatively simple and they significantly reduce the numerical errors affecting the provided solution. This fact has been confirmed via a comparison of the computational results with the available exact solutions.
Regarding the 2D advection-diffusion transport equation, the Strang technique has been used for its splitting into a series of the 1D advection-diffusion transport equations. Next for the solution of all obtained equations, the method proposed for the 1D equation was implemented as a solver. The combination of both numerical techniques provided also a robust tool for the numerical solution of the 2D advection-diffusion transport equation even in the case when the flow velocities vary remarkably in the region of solution. The effectiveness of the proposed approach confirmed for the 2D advection-diffusion transport equation suggests that the presented approach can also be easily developed for the 3D case.
The results provided by the proposed method consistent with the exact solutions suggests that probably it will be effective also in the cases dealing with the real transport phenomena in environmental engineering. To this order, the method should be associated with the model of transient flow or steady non-uniform flow in order to calculate the velocity field. The positive properties of the method shown in the paper will be preserved if while computation the weighting parameter ω will be in the range from 0.5 to 1.0 for 1D flow or its value will belong to the shaded area presented in Figure 2 for 2D flow. However, to recognize the behaviour of the method in realistic cases, additional studies are needed.

DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.