B-Spline Method of Lines for Simulation of Contaminant Transport in Groundwater

: In this study, we propose a new numerical method, which can be e ﬀ ectively applied to the advection-dispersion equation, based on B-spline functions and method of lines approach. In the proposed approach, spatial derivatives are calculated using quintic B-spline functions. Thanks to the method of lines approach, the partial di ﬀ erential equation governing the contaminant transport in groundwater is converted into time-dependent ordinary di ﬀ erential equations. After this transformation, the time-integration of this system is realized by using an adaptive Runge–Kutta formula. In order to test the accuracy of the proposed method, four numerical examples were solved and the obtained results compared with various analytical and numerical solutions given in the literature. It is proven that the proposed method is faster and more reliable than other methods referenced herein and is a good alternative for simulation of contaminant transport problems as a result of these comparisons.


Introduction
Groundwater is one of the most important resources in the world. This resource is used as clean water in many different fields such as domestic, industrial and agricultural. Especially in rural areas, most of the people's need for clean water is met by groundwater. The need for clean water is increasing with the increase in the world population and the development of industrial and agricultural technologies. Therefore, both the quantity and the quality of groundwater are very important, because contaminated groundwater loses its intended use. In such cases, remediation works should be carried out for contaminated groundwater. In these studies, modeling of contaminant transport in the relevant region is very important as it will provide an analysis of the current situation and provide information about the future status of the current pollution [1].
Groundwater contaminant transport modeling is carried out by solving the partial differential equation which represents the physical processes occur in nature, under certain boundary and initial conditions using analytical or numerical methods. Analytical models produce continuous results in the solution space, but these models are available only for simplified cases with simple geometry and constant parameters. However, the geometry of contamination problems encountered in the real world is complex, and its parameters vary depending on location and time. For this reason, it is mandatory to use numerical models in real-world applications. Although numerical models produce discrete solutions at the nodes, they are preferred because they can provide solutions when the geometry is complex and the parameters are variable [2].
In last two decades many different numerical models have been developed for the solution of contaminant transport problem. These models are based on different approaches such as finite difference method [3][4][5][6], finite element method [7], finite volume method [8], boundary element method [9],

Governing Equations
The partial differential equation describes the transport of a single chemical constituent in saturated porous media in two-dimensions, considering the advection-dispersion process is given by [2].
where V x and V y are seepage velocities in x and y directions (LT −1 ); D xx , D xy , D yx and D yy are the elements of the dispersion coefficient tensor (L 2 T −1 ); C is dissolved concentration (MT −3 ); x and y are cartesian coordinates (L); t is time (T). The elements of the dispersion coefficient tensor D xx , D xy , D yx and D yy can be calculated by using seepage velocities and longitudinal (α L ) and transverse (α T ) dispersivities from relations given in Equation (2) [28].
where D * is the effective molecular diffusion coefficient. The molecular diffusion coefficient is extremely small compared to the mechanical dispersion coefficient. For this reason, it can be neglected in calculations.

B-Spline Method of Lines
In the proposed method, MOL is used with the QBS functions for numerical discretization of the Equation (1). Firstly, space derivatives are approximated with the help of the QBS functions. Then the approximated space derivatives are written in Equation (1). Thus, the problem transforms to time-dependent system of ordinary differential equations. This system is solved with the help of DOPRI5. The formulation of this solution is presented in the following sections.

Spatial Discretization
Let the two-dimensional solution domain partitioned into a mesh of uniform lengths ∆x = x i+1 − x i and ∆y = y j+1 − y j , by the nodes x i , y j where i = 0, 1, . . . , N and j = 0, 1, . . . , M. Our goal is to construct QBS functions for every vector in both directions. Then, space derivatives at the nodes can be calculated easily without needing complicated formulations. It is clear that solution domain consists of M + 1 row vectors and N + 1 column vectors. Therefore, (M + 1) × (N + 1) times numerical approximations for contaminant concentrations are established with QBS functions for both xand y-directions.
The concentration values C(x, y, t) can be approximated by C(x, y, t) with QBS functions. The mathematical formulation of QBS functions is given below considering the row and column vectors in xand y-directions.
where δ j,i (t) and γ j,i (t) are time-dependent parameter matrices that are needed to be calculated by using initial and boundary conditions. B i (x) and Q j (y) are the QBS basis functions for xand y-directions given in Equations (5) and (6), respectively [29].
Water 2020, 12, 1607 The values of QBS basis functions and its first and second order derivatives at nodal points are given in Tables 1 and 2. Obviously, third and fourth order derivatives can be calculated as well but as it is clear from Equation (1) that only first and second order derivatives are needed. Table 1. The values of B i (x) and its derivatives at nodal points.
0 20/∆x 2 40/∆x 2 −120/∆x 2 40/∆x 2 20/∆x 2 0 Approximate concentration values and its first and second order derivatives at the nodes in both directions can be calculated in terms of δ j,i (t) and γ j,i (t) by using approximate functions in Equations (3) and (4) with the values of QBS functions given in Tables 1 and 2. The steps required to find the δ j,i (t) and γ j,i (t) values are presented in detail in Appendix A.

Temporal Discretization
MOL is a powerful approach in which partial differential equation transforms to ordinary differential equation system. The main advantage of this approach is that in the integration of ordinary differential equation system well-established, robust, thoroughly tested integrators can be used [30]. When the spatial derivatives of the contaminant transport equation given in Equation (1) are discretized by QBS functions, the time-dependent ordinary differential equation system given below is obtained.
where L indicates a spatial differential operator. A fifth-order adaptive Runge-Kutta formula, known as DOPRI5 [27], is used for the time integration of the Equation (7) which consists of (M + 1) × (N + 1) ordinary differential equations. The DOPRI5 is an adaptive integration method used in the solution of ordinary differential equations and is a member of the Runge-Kutta family [27]. It has seven stages, but it uses only six function evaluations per step because it has the FSAL (First Same as Last) property: the last stage is evaluated at the same point as the first stage of the next step. The coefficients in the DOPRI5 are chosen to minimize the error of the fifth order solution. This is the main difference with the Fehlberg method [31], which was constructed so that the fourth-order solution has a small error. For this reason, the DOPRI5 method is more suitable when the higher-order solution is used to continue the integration, a practice known as local extrapolation [32]. Such adaptive time-integration methods always keep the truncation errors below a certain value for each step and prevent the errors from growing uncontrollably. In other words, the numerical errors are prevented from overgrowth, and the stability of the numerical method is maintained. The numerical concentration values for successive time steps are calculated with the DOPRI5 scheme as follows: where p and ν represent time and stage indexes, k ν is the approximated slope matrix calculated at each stage, ∆t p is the dynamically determined time step at time t p , the coefficients ω ν , φ ν,ξ , ψ ν are elements of Butcher array which is given in Table 3 [33].  The local truncation error for the DOPRI5 scheme is obtained approximately by the formula given in Equation (11). The infinity norm of this error matrix is compared to a user-defined maximum allowable error tolerance. If the e p+1 ∞ ≤ e tol condition is satisfied, the time step is calculated according to the formulation given in Equation (12). Otherwise, ∆t p+1 is iteratively updated using Equation (12) until e p+1 ∞ ≤ e tol condition is satisfied. In this study, the ∆t p+1 /∆t p ratio was limited to [0.1, 10] and e tol = 1 × 10 −6 .

Numerical Applications
In this section, the success of the proposed method (QBS-DOPRI5) will be tested on four different problems, which are widely used in the literature. All problems were solved with QBS-DOPRI5 and their results were compared with analytical or numerical solutions given in the literature. In the following section, detailed information about this comparison will be presented. All simulations are performed on a computer with an i5-4460 processor with 3.20 GHz speed and 8 GB RAM.

Example 1
Consider Equation (1) with constant parameters such as V x = V y = 0.8 m/s and D xx = D yy = 0.01 m 2 /s, D xy = D yx = 0 m 2 /s and the solution domain of 0 < x, y < 2. Analytical solution of this problem was given by Kalita et al. [34].
Water 2020, 12, 1607 6 of 19 In this problem, a Gaussian pulse is considered as an initial condition with the peak value of unity. Location of the peak value is indicated by (x,ŷ), which is taken as (0.5, 0.5) in the experiments. The initial and boundary conditions can be taken from the analytical solution. To compare our results with the other proposed numerical methods in the literature, the same simulation parameters such as flow velocities, dispersion coefficients and grid sizes are used. Example 1 has been studied by many researchers to test their proposed methods. Comparisons were made in terms of average and maximum absolute errors. Average absolute error values at the time t = 1 s over the domain [1,2] × [1,2] with V x = V y = 0.8 m/s and D xx = D yy = 0.005 m 2 /s are presented in Table 4. QBS-DOPRI5 is produced better average absolute error values than the other methods given in Table 4.  The parameters used in Table 5 are V x = V y = 0.8 m/s, D xx = D yy = 0.01 m 2 /s, D xx = D yy = 0.001 m 2 /s and ∆x = ∆y = 0.025 m. In Table 5, both average absolute error and maximum error values are presented and compared with the results of the other available numerical methods in the literature. It is clear that QBS-DOPRI5 produces error values similar to or lower than other methods. Moreover, the results obtained with QBS-DOPRI5 at time t = 1.25 s for Figure 1 in which Gaussian shape of the concentration can be seen clearly. The contour lines of these results are plotted in Figure 2.
It should be noted that the smallest time step chosen in the methods presented in Table 5 (except QBS-DOPRI5) is ∆t = 0.00625 s. In the QBS-Euler method, the time step is reduced by 50 times and taken into consideration as ∆t = 0.000125 s. The average absolute error and maximum absolute error values obtained with the QBS-Euler method for D xx = D yy = 0.01 m 2 /s are 2.435 × 10 −5 and 5.572 × 10 −4 , respectively. It is clear that the error values given in Table 5 for QBS-DOPRI5 are considerably smaller than the error values obtained with QBS-Euler. CPU times required for QBS-Euler and QBS-DOPRI5 are obtained as 335.57 s and 17.48 s, respectively. As a result of this analysis, it is necessary to use high-order schemes together in space and time, both in terms of the quality of the solutions and the reduction of the CPU time.

Example 2
In this example, migration of chloride ion in landfill leachate through a narrow, relatively thin, valley-fill aquifer is considered in which the simulation is governed by Equation (1). This example is

Example 2
In this example, migration of chloride ion in landfill leachate through a narrow, relatively thin, valley-fill aquifer is considered in which the simulation is governed by Equation (1). This example is

Example 2
In this example, migration of chloride ion in landfill leachate through a narrow, relatively thin, valley-fill aquifer is considered in which the simulation is governed by Equation (1). This example is taken from the work of Wexler [38], the sample problem 6. Boundary conditions of the problem are as follows: 0,x = x 0 and y ŷ 1 or y ŷ 2 (14) Water 2020, 12, 1607 and where, y M = 1200 m is the aquifer width,ŷ 1 = 300 m,ŷ 2 = 800 m, x N = 1500 m is the aquifer length,Ĉ = 1 kg/m 2 , flow velocities are V x = 3 × 10 −6 m/s, V y = 0 and dispersion coefficients are D xx = 2 × 10 −4 m 2 /s, D yy = 6 × 10 −5 m 2 /s. According to the assumptions of Wexler [38], D xy and D yx dispersion coefficients are taken as zero. Under these conditions, analytical solution of the problem given by Wexler [38] C(x, y, t) =Ĉ ∞ n=0 L n P n cos(ηy) exp Simulations are carried out for different times, i.e., t = 1500 days, t = 3000 days. Contour lines of the approximated concentration results for different simulation times are plotted in Figures 3 and 4 together with the analytical solution, showing very good agreement. Moreover, central concentration profiles are compared and plotted in Figure 5. An excellent agreement was observed.
Water 2020, 12, x FOR PEER REVIEW 8 of 21 taken from the work of Wexler [38], the sample problem 6. Boundary conditions of the problem are as follows: ( , , ) = , = and ≤ ≤ 0, = and < or > and where, = 1200 m is the aquifer width, = 300 m, = 800 m, = 1500 m is the aquifer length, = 1 kg/m , flow velocities are = 3 × 10 m/s, = 0 and dispersion coefficients are = 2 × 10 m /s, = 6 × 10 m /s. According to the assumptions of Wexler [38], and dispersion coefficients are taken as zero. Under these conditions, analytical solution of the problem given by Wexler [38] in which

Example 3
In this example, an unsymmetrical continuous solute line source at the left boundary (x = x 0 ) is considered. The concentrations of the solute source at the boundary are calculated by the following formulation [9,11].      C(x 0 , y, t) =Ĉ

Example 3
where,Ĉ is the peak value of concentration at the line source and assumed to be unity. This way obtained concentration values can be considered as relative concentrations.ŷ = 125 m is the location of the peak concentration value at solute line source. The other boundary conditions are assumed to be no flux and mathematical representation of them are as follows: where y M = 300 m and x N = 600 m are the width and length of the aquifer field, respectively. Two-dimensional velocity field is considered. The simulations are conducted by taking velocities as V x = 1.1784 m/d and V y = 0.3157 m/d. Dispersion coefficients are calculated from Equation (2) by taking longitudinal dispersivity, α L = 6.248 m, and transverse dispersivity, α T = 0.393 m. As can be seen from the parameters, the aquifer is homogeneous and anisotropic. For this reason, the ∂ ∂x D xy

∂C ∂y
and ∂ ∂y D yx ∂C ∂x terms in Equation (1) must be included in the solution. However, in the works of Eldho and Rao [9] and Boddula and Eldho [11], it is seen that these terms are not included in the solution. Two different scenarios are examined for this example. In the first scenario (Scenario A), these terms are neglected in order to compare the solutions under the same conditions. In the second scenario (Scenario B), these terms are taken into account in order to see the real behavior of the problem.
Since there is no available analytical solution for this problem, FDM-DOPRI5 was developed and used to obtain a reference solution by taking space grid sizes considerably small such as h = ∆x = ∆y = 0.5 m for scenario A and h = ∆x = ∆y = 0.25 m for scenario B. These grid sizes mean 1201 × 601 nodes for scenario A and 2401 × 1201 nodes for scenario B are presented in the solution domain. The solutions obtained with these many nodes can be considered as almost exact solutions.
The obtained results with the QBS-DOPRI5 considering Scenario A are compared with the other solutions taken from the literature in Table 6 for different grid sizes. QBS-DOPRI5 produced very good results even when the grid size h = 25 m, which means 25 × 13 nodes in the solution domain. For h = 5 m, meaning 121 × 61 nodes, the obtained concentration values with QBS-DOPRI5 are almost identical with reference solution. Contour lines of the approximated solution obtained by QBS-DOPRI5 and reference solution are plotted in Figure 6. As clearly seen from Figure 6, the contour lines of the solutions are identical.  Table 7 was created to show that the proposed method achieved the same accuracy in a shorter time than standard numerical methods. As can be seen from Table 7, solutions for QBS-DOPRI5 and FDM-DOPRI5 for different h values were obtained, and their concentration values at different locations were given. CPU times spent for these solutions are also presented in the same table. When Table 7 is examined carefully, the solutions obtained for QBS-DOPRI5 with h = 25/4 m and FDM-DOPRI5 with h = 25/32 m appear to be the same as the reference solution. Comparing the CPU times (19.57 s and 1302.46 s) obtained in this case, it is obvious that the QBS-DOPRI5 method is 66.55 times faster than the FDM-DOPRI5 method.  Table 7 was created to show that the proposed method achieved the same accuracy in a shorter time than standard numerical methods. As can be seen from Table 7, solutions for QBS-DOPRI5 and FDM-DOPRI5 for different h values were obtained, and their concentration values at different locations were given. CPU times spent for these solutions are also presented in the same table. When Table 7 is examined carefully, the solutions obtained for QBS-DOPRI5 with ℎ = 25/4 m and FDM-DOPRI5 with ℎ = 25/32 m appear to be the same as the reference solution. Comparing the CPU times (19.57 s and 1302.46 s) obtained in this case, it is obvious that the QBS-DOPRI5 method is 66.55 times faster than the FDM-DOPRI5 method.

Method
( )   When the same problem is solved considering scenario B, the results are presented in Table 8. As can be seen from Table 8, QBS-DOPRI5 with h = 25/4 m and FDM-DOPRI5 with h = 25/64 m reach the same accuracy as the reference solution. Comparing the CPU times spent under these conditions, QBS-DOPRI5 was found to be 1061.82 times faster than FDM-DOPRI5. Moreover, this example is solved with COMSOL Multiphysics, the well-known FEM commercial code, under scenario B. This solution was obtained by activating the adaptive mesh refinement property in the COMSOL Multiphysics program and using 30763 triangular elements. It should be noted that QBS-DOPRI5 and COMSOL Multiphysics results are exactly the same. These results are illustrated in Figure 7. Also, COMSOL Multiphysics output results used to plot Figure 7 are presented in Table S1 file in Supplementary Materials. COMSOL Multiphysics program and using 30763 triangular elements. It should be noted that QBS-DOPRI5 and COMSOL Multiphysics results are exactly the same. These results are illustrated in Figure 7. Also, COMSOL Multiphysics output results used to plot Figure 7 are presented in Table S1 file in Supplementary Materials.

Example 4
In reality, the aquifers are complex and can vary significantly in terms of physical parameters, such as permeability and diffusivity. For this reason, an example that is given in the work of Perko [17] is studied. In this example, anisotropic dispersion with spatially varying velocity field is considered for testing the proposed method's robustness. The problem is invariant in the zdirection. Hence, it can be considered as a 2D problem. The velocity components in xand ydirections are as follows: A rectangular flow domain is considered with the dimensions −1 ≤ x ≤ 1 and 0 ≤ y ≤ 1. Dirichlet boundary condition is considered along the left bottom boundary (−1 ≤ x ≤ 0, y = 0). The concentration profile at the Dirichlet boundary varies from 0 to 1 as a sharp pulse-like shape. Open boundary condition is defined along the right bottom boundary (0 ≤ x ≤ 1, y = 0). The Dirichlet boundary concentration profile is defined as where the parameter Υ determines the sharpness of the concentration profile. The value of Υ is taken as 100, which results in a sharp profile. Neumann boundary condition with zero flux is defined at the left, right and upper boundaries.
Since this example also does not have an analytical solution, a reference solution is needed to compare QBS-DOPRI5 results. For this purpose, a reference solution is obtained by activating the adaptive mesh refinement property within the COMSOL Multiphysics software. In this example, two different scenarios are created by changing the ratio of the longitudinal dispersivity to transverse dispersivity coefficients. In the first scenario (scenario A), a typical ratio of 10 is considered by taking α L = 0.1 m and α T = 0.01 m. In the second scenario (scenario B), this ratio is taken as 1000 in order to test the real potential of the proposed method. In this scenario, dispersivities are selected as α L = 1 m and α T = 0.001 m. COMSOL Multiphysics models contain 25,806 and 106,268 triangular elements in the solution domain for Scenarios A and B, respectively.
The contour plot of the solutions obtained with QBS-DOPRI5 and COMSOL Multiphysics under Scenarios A and B are presented in Figures 8 and 9, respectively. When Figures 8 and 9 are examined together, it is clear that the dissipation of the concentration profile is greater in Scenario A. In Scenario B, the sharp concentration profile at the inflow tends to be preserved.    Moreover, the concentration profiles at the inlet and outlet boundaries, which are taken from the obtained solutions, are illustrated in Figures 10 and 11, respectively. The effects of scenarios A and B conditions on contaminant transport are clearly seen in Figures 10 and 11. It is clear that the QBS-DOPRI5 and COMSOL Multiphysics solutions are in good agreement for both scenarios. Also, Moreover, the concentration profiles at the inlet and outlet boundaries, which are taken from the obtained solutions, are illustrated in Figures 10 and 11, respectively. The effects of scenarios A and B conditions on contaminant transport are clearly seen in Figures 10 and 11. It is clear that the QBS-DOPRI5 and COMSOL Multiphysics solutions are in good agreement for both scenarios. Also, COMSOL Multiphysics output results used to plot Figures 8 and 9 are presented in Tables S2 and S3 file in Supplementary Materials, respectively. Moreover, the numerical values of the concentration profiles at the open boundary are presented separately in Table 9 for both scenarios.  Table 9 for both scenarios.

Conclusions
In this paper, a new method of line approach based on B-spline functions is proposed for accurate simulation of two-dimensional contaminant transport problems. In the proposed method, quintic B-spline functions were used to approximate spatial derivatives. After the approximation of spatial derivatives, the partial differential equation representing the contaminant transport problem was converted to a time-dependent ordinary differential equation at each grid point. The system of ordinary differential equation system obtained for all grid points was solved with an adaptive time-integration scheme called DOPRI5. The biggest advantage of the proposed method is to convert the governing partial differential equation to the ordinary differential equation system and to obtain accurate solutions in an easier way and in a shorter time than conventional numerical methods such as finite difference, finite element, etc.
In this study, four numerical applications were carried out. The first two examples have analytical solutions, and the success of the proposed method on these examples is compared with the results of the studies available in the literature. As a result of this comparison, it was seen that the results obtained with QBS-DOPRI5 were compatible with analytical solutions and produced better results than other studies except for the compact finite difference method. In the third example, a hypothetical aquifer model with homogeneous and anisotropic dispersion coefficients without an analytical solution is solved. In order to question the quality of the results of the QBS-DOPRI5 model, the FDM-DOPRI5 model, which has an extremely high number of grids and can be considered almost an analytical solution, was used as a reference solution. In this example, the results of the proposed method are compared with the results obtained with the BEM and MLPG-MLS methods. As a result of this comparison, it has been determined that the proposed method produces the same results as the reference solution, whereas the solutions obtained by the BEM and MLPG-MLS methods are far from the reference solution. It is observed that the proposed method is approximately 1000 times faster than the FDM-DOPRI5 method, provided that the solution with the same accuracy is obtained. In addition, in order to test the reliability of the reference solution, this problem was solved with the COMSOL Multiphysics program. It should be noted that the adaptive mesh refinement property has been activated in the COMSOL Multiphysics program, and a large number of triangular elements have been used.
As is well known by groundwater hydrologists, aquifer parameters such as permeability and dispersivity can vary significantly in field conditions. Therefore, anisotropic dispersion with spatially varying velocity field was taken into account in the last example to question the usability of the proposed method in real applications. Since this example also does not have an analytical solution, a reference solution for determining the quality of QBS-DOPRI5 results was obtained with the help of the COMSOL Multiphysics program. For this example, two different scenarios have been created, and the results obtained with the proposed method under both scenarios are considered to be quite accurate and can be used safely in real applications. The proposed method of lines approach uses uniform mesh in space and adaptive mesh in time. Especially in the case of contaminant transport problems with sharp front, non-uniform mesh or adaptive mesh should be preferred in the space. By making a small modification in the code, these mesh structures can be easily created. The proposed method can also be used easily in the solution of three-dimensional contaminant transport problems and does not cause any changes in the solution algorithm. In addition, it can be safely preferred for the solution of engineering problems involving many partial differential equations due to the simplicity of the solution algorithm. The authors are planning to use the proposed method in several groundwater management problems.  Acknowledgments: In this study, a part of the first author's PhD thesis study is presented.

Conflicts of Interest:
The authors declare no conflict of interest.
where f, F, g, G are boundary matrices consisting of first and second derivatives as scalar values, either calculated from FD schemes or known from boundary conditions. From Equations (A2-A3) and (A7-A10) following relations can be written.
Consequently, the first and second derivatives in xand y-directions at each nodal point can be calculated by using Equations (A2-A3) and (A5-A6), respectively since the parameter matrices are known.