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

Apart from the transport of contaminations, similar equation describes other flow phenomena as flood wave propagation (Gąsiorowski 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’.

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-to-date 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 high-order 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 & Gąsiorowski 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.

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).

Detailed information on the FEM is provided, e.g. by Zienkiewicz (1972), Fletcher (1991) and Gresho & Sani (1998).

The initial-value problem for the system (12) is solved using the implicit TR ⁠, where n is the index time level).

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.

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.

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.

In its standard form, the Strang splitting process consists of the three following steps:

• Ist step – time integration in x-direction over half of the time step Δt/2:

  

  (28)
• IInd step – time integration in y-direction over the time step Δt:

  

  (29)
• IIIrd step – time integration in x-direction over half of the time step Δt/2:

  

  (30)

The solution obtained in the IIIrd step is the sought solution of the governing equation at the next time level, i.e. c(t + Δt) = c⁽³⁾(t + Δt).

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).

Equations (28)–(30) are solved successively along all rows parallel to the x-axis (1 ≤ i ≤ N), all columns parallel to the y-axis (1 ≤ j ≤ M) and once more along all rows (1 ≤ i ≤ N). To this order, the method of solution 1D transport equation presented previously can be used. If in the considered flow case variable flow velocities are expected, then for the solution of Equations (28)–(30) the equations of type of Equation (12) with variable weighting parameter ω should be used. All final systems of algebraic equations have the tri-diagonal matrices. Their dimensions are M × M for the consecutive rows and N × N for the consecutive columns.

For an illustration of how the presented approaches work, a couple of examples dealing with various initial and boundary conditions are presented below. In all performed examples, the numerical solutions are compared with the analytical solutions available for the considered cases of contaminant transport. Moreover, the Theis equation describing groundwater flowing towards a well is also considered.

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²/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.

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²/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.

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.

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⁻¹).

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²/s. The imposed pulse (Figure 8) is defined by the following value of its parameters: t₁ = 5 s and t₂ = 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⁻¹. 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.⧫

The so-called well function W(u), where u = r²·S/(4T·t), can be approximated by the following formulas (McElwee 1980):

• W(u) = −ln(u) − 0.57721566 + 0.99999193u − 0.24991055u² + 0.05519968u³ + 0.00976004u⁴ + 0.00107857u⁵ for 0 ≤ u ≤ 1.

• W(u) = (exp(−u)/u)·(u² + 2.334733u + 0.250621)/(u² + 3.330657u + 1.681534) for u > 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₀ where r₀ 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.

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²/s is constant, whereas the diffusive number is C_dx = D_x·Δt/Δx² = 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.⧫

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 known – the 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.

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 φ = 0). The results obtained for Δt = 1.0 s are displayed in Figure 11.

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 φ = 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.

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²/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 13 – curve a).

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.

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