Mathematical Modeling of Solutes Migration under the Conditions of Groundwater Filtration by the Model with the k-Caputo Fractional Derivative

Within the framework of a new mathematical model of convective diffusion with the k-Caputo derivative, we simulate the dynamics of anomalous soluble substances migration under the conditions of two-dimensional steady-state plane-vertical filtration with a free surface. As a corresponding filtration scheme, we consider the scheme for the spread of pollution from rivers, canals, or storages of industrial wastes. On the base of a locally one-dimensional finite-difference scheme, we develop a numerical method for obtaining solutions of boundary value problem for fractional differential equation with k-Caputo derivative with respect to the time variable that describes the convective diffusion of salt solution. The results of numerical experiments on modeling the dynamics of the considered process are presented. The results that show an existence of a time lag in the process of diffusion field formation are presented.


Introduction
In this paper, we study the problem of modeling convective diffusion of solutes under the conditions of plane-vertical steady-state filtration with a free surface.Such problems arise, for example, when dealing with desalinization of soils during land reclamation or desalinization of groundwater polluted by industrial or domestic sewage [1,2].The theory and practice of mathematical modeling of such problems within the framework of classical differential models are now well-developed and verified [1][2][3][4].It should be noted that mathematically analogous statements of boundary value problems are also encountered within the theory of heat and mass transfer of solutes and gas mixtures.
At present, the challenging issue is to increase the adequacy of classical quantitative models of heat and mass transfer in systems with a complex space-time structure characterized by memory effects, spatial non-locality, and self-organization.There is a tendency to revise the main principles of the classical theory of heat and mass transfer.In particular, significant progress while modeling convective diffusion in anomalous conditions and several other processes in the continuum mechanics was achieved using the fractional-order integro-differentiation formalism [5][6][7][8][9][10][11][12][13].These transport processes are strongly non-local in time and (or) space, hence mathematical models have the form of differential equations of fractional order.
While modeling anomalous diffusion, parameters of fractional derivatives represent the properties of interactions between particles and medium.In practice, they are often determined from experimental data.As an increased number of parameters in generalizations of derivatives can increase the accuracy of such fitting, studies of diffusion equation with multiparameter derivatives, in particular the k-Caputo derivative, are topical.
Considering a multi-dimensional problem, we use numerical techniques that allow obtaining solutions faster without significant reduction of accuracy.The methods most used in such case are finite-difference splitting schemes that reduce high-dimensional problems to a set of one-dimensional ones.Such schemes for the classical models of diffusion were thoroughly studied in [14].The unique solvability, stability and convergence of a splitting scheme for the fractional diffusion equation with the Caputo derivative with respect to the time variable was proved in [15].Compared with the conventional methods (e.g., [16]) that discretize fractional differential equation to a single system of linear or non-linear algebraic equations, splitting schemes result in a series of independent equation systems of lower size that allows developing simple and efficient computational algorithms [17].
In Section 2 we propose a new two-dimensional mathematical model of anomalous convective diffusion under the conditions of steady-state filtration with a free surface and pose the corresponding boundary value problem.As a filtration scheme, we use a scheme for the spread of pollution from rivers, canals, or storages of industrial wastes [1].In Section 3 a numerical method based on the locally one-dimensional finite difference scheme of A.A. Samarsky [14] is described.In Section 4 we give the results of convective diffusion process simulation and compare them with the results obtained using the classical mathematical model and the model with the standard Caputo derivative.

Mathematical Model and the Corresponding Boundary Value Problem
When modeling the peculiarities of anomalous transport dynamics, fractional-differential mathematical models, that is, the models with fractional derivatives with respect to the time and/or spatial variables are widely used [18][19][20][21][22].In the case of mathematical description of non-local convective diffusion dynamics, a model based on the following fractional-differential counterpart of the classical convective diffusion equation [17,[23][24][25][26] is very common: where C = C (ϕ, ψ, t) is the concentration, υ is the velocity of convection, d is the diffusion coefficient, σ is the porosity of a medium, D t is the operator of the Caputo fractional derivative [5][6][7][8] of order α(0 < α < 1), ∆ is the Laplace operator.
However, when studying non-local processes of convective diffusion in saturated porous media with a complex space-time structure, mathematical model based on Equation (1) with a parameter α may not be adequate.As it is noted in [27], in the case when a simulated medium consists of layers with substantially different properties (in particular, with different permeability parameters), the direct use of this one-parameter fractional model without sufficiently complete substantiation of the fulfillment of its applicability conditions may lead to unacceptable errors.For this case, it is proposed [27] to make a transition from the one-parameter model to a corresponding model with a larger number of parameters that allows in some cases increasing modeling accuracy by their proper selection.
Performing mathematical modeling of the dynamics of non-equilibrium in time convective diffusion process, we start from the following generalization of Equation ( 1): where k D α t is the operator of the fractional derivative with two parameters, the so-called k-Caputo derivative of order α (0 < α < 1) that is defined as [28] It should be noted that the widely used Equation ( 1) is a particular case of (2) and follows from the latter when k → 1.
Under the conditions of filtration in a potential velocity field, taking into account the continuity equation div υ = 0, from (2) we finally obtain where υ = (υ x , υ y ) = ∇ϕ, ϕ is the filtration velocity potential.
In particular, when α, k → 1 from (3) we obtain the classical equation of convective diffusion [1,2] Suppose that for the considered filtration scheme, the domain of complex flow potential ω = ϕ + iψ (ψ is the flow function) is known, as well as the solution of the corresponding filtration problem, i.e., the characteristic flow function z = f (ω) (for many practical filtration schemes, such a function is given, for example, in [1]).Then, transforming Equation ( 3) to new independent variables-points of geometrically simpler domain of complex flow potential (ϕ, ψ), we rewrite it in the form where Let us consider a filtration scheme corresponding to the problem of convective diffusion of pollution from rivers, canals, or surface storages of industrial effluents (Figure 1).For this scheme, the domain of complex flow potential G ω has the form of a half-strip (Figure 2) and the solution of the corresponding filtration problem is given as [1] x = He The problem of modeling pollutants migration in the right side G z of the symmetric filtration region (Figure 1) can be mathematically formulated as the problem of finding the solution of Equation ( 4) ((x, y, t) ∈ G z × (0, +∞)) under the following boundary conditions: = 0, C| t=0 = 0 where C 1 is the given concentration of soluble substances on filtration inflow, n is the outer normal to the corresponding line.
As the filtration domain G z is a domain with a partially unknown boundary, we make a transition to the variables (ϕ, ψ) that represent the points of the domain G ω = {(ϕ, ψ) : 0 < ϕ < +∞, 0 < ψ < Q} (Figure 2).Finally, the boundary value problem in the domain of complex flow potential is as follows: We solve the considered problem under the following restrictions on the parameters of the derivative: 0 < α < 1, k > 1 − α.It is easy to show that in this case we have is the operator of the standard Caputo derivative of order β with respect to the variable t [5][6][7][8], β = k+α−1 k (0 < β < 1).Then, Equation ( 5) can be rewritten as where υ2 (ϕ, ψ) = k 1−β υ 2 (ϕ, ψ) and the simulation of the dynamics of the considered anomalous convective diffusion process reduces to the solution of the boundary value problem for Equation (7) with the conditions (6) in the domain of complex flow potential.Let us introduce the following dimensionless variables and parameters: where L is the parameter of scale, υ 0 is the characteristic parameter of velocity.
Performing transition to the dimensionless variables ( 8) in ( 7), (6) (further dropping the prime sign over the dimensionless variables), we obtain in the domain G ω × (0, +∞) the boundary value problem where The transition from the domain of complex flow potential (Figure 2) to the physical domain (Figure 1) is carried out using the following formulas:

Approximate Solution of the Boundary Value Problem
In this section we present a finite-difference method for obtaining an approximate solution of the boundary value problem ( 9), (10).
Let us introduce the grid domain where h 1 = 2ϕ 0 2m+1 , h 2 = 1 n are the grid steps with respect to the geometric variables ϕ, ψ; τ is the grid step with respect to the time variable; ϕ 0 = const.
Restricting the domain of complex flow potential from the right by the line ϕ = ϕ 0 (ϕ 0 1) and posing an additional boundary condition on it (for example, the Neumann condition), we associate with the considered boundary value problem the following locally one-dimensional finite-difference scheme [14,29]: where Ĉ are the finite-difference counterpart of the fractional differentiation operator defined as is the Gamma function [30].Let us note that in the class of sufficiently smooth functions we have the following error estimate: Unrolling the finite-difference operators in (12) and collecting similar terms, we obtain on the half-integer step t j+ 1 2 the following series of systems of linear algebraic equations where The solution of the system ( 14) can be found by the Thomas algorithm as where We obtain the following starting values of the coefficients α, β using the finite-difference counterpart C j 0l = 1 of the boundary condition in the point ϕ = 0: m+1,l (l = 1, n; j = 0, N) of the second-order boundary condition in the point ϕ = ϕ 0 using (15) we have After obtaining the solution on a half-integer time step, on an integer time step from ( 13) we obtain where .
The solution of the system (16) can be found as Taking into account the finite-difference form of the boundary condition at ϕ = 0 we obtain the following starting values of the coefficients: These relationships make it possible to compute the solution on an integer time step.On each step of the solution procedure, we have series of three-diagonal linear systems that are independent and, in the case of large grids, can be efficiently solved in parallel.As we consider a time-fractional model, the biggest impact on computational complexity of the scheme has the computation of the right-hand side vectors Φ j il and Ω j il that can be further optimized and parallelized.The stability of the Thomas algorithm for ( 14) and ( 16) follows from the fact of the diagonal predominance in the matrices of the coefficients of the linear systems and can be controlled by a proper choice of grid steps.
The subsequent transition to the physical domain is carried out according to (11).

Results of Numerical Experiments
Numerical modeling of the dynamics of solute migration using the considered non-classical model of convective diffusion is performed with respect to the dimensionless variables (8) for the input data from [26].With the dimensionless time step equal to 0.1 we performed the simulation for the dimensionless time range T ≤ 100.Experimentally checking the convergence of the scheme, we solve the problem increasing the grid size starting from 50 × 50.Maximal local relative error started changing not more than by 2. During the simulation, the time spent on computations on one step, as it is expected, increases from t < 10 ms on the first step to t = 1820 ms on the 1000-th time step.In the additional experiments, for the grid size equal to 30 × 30 computation time on one step increase up to t = 90 and for the grid size of 65 × 65-up to t = 680.The increase was linear for time steps less than 300, but become slightly non-linear in further computations.This effect can be explained by non-optimal cache usage [17].Total computation time for the dimensionless time range T ≤ 30 quadratically depended on the grid size.For the grid size equal to 100 × 100, time spent on computing the right-hand side vectors became more than 90% of the total computation time starting from the 125-th time step, for the grid size of 65 × 65-starting from the 231-th time step.For the smallest tested grid of 30 × 30 nodes, time spent on computing right-hand side vector became more than 80% of the total computation time starting from the 451-th time step.Some of the obtained results are depicted in Figures 3-5.The analysis of the results of the numerical experiments allows us to draw the following conclusions about the dynamics of solute concentration field modeled using the considered convective diffusion model:

•
There is a delay in the formation of concentration field while modeling the diffusion process using the model with the k-Caputo derivative comparing both to the case of modeling this process using the classical model [1,2] and the model based on the standard definition of the Caputo derivative (Figure 3).

•
For fixed values of order α (0 < α < 1), when the value of the parameter k varies, this delay is the greater, the smaller the value of the parameter k is (Figures 3 and 4).

•
For fixed values of the parameter k, when the value of α varies, there is also a delay in the advance of concentration front for the diffusion model with the k-Caputo derivative compared with the case of the classical model.This delay is the greater, the smaller the value of the parameter α is (curves 2-4 in Figure 5).The performed computations also allow us to conclude that the dynamics of concentration front's delay is determined to a greater extent by a change in the magnitude of order α than by a change of the value of the parameter k.

Conclusions
Consideration of a temporal non-locality of the convective diffusion in porous media with a complex internal structure on the base of the proposed model with the k-Caputo derivative opens the possibility of mathematical modeling of the phenomenon of a significant time lag in the process of diffusion field formation compared with the case when this process is described by the classical mathematical model or the model with the standard Caputo fractional derivative.Ignoring the phenomenon of a temporal non-locality of anomalous convective diffusion while developing engineering solutions (for example, in the field of designing systems for the environmentally safe operation of surface sewage storages in complex geological conditions) can lead to errors in forecasts of the degree of specified objects safety.

Figure 1 .
Figure 1.Scheme of filtration in a physical domain.

Figure 2 .
Figure 2. Scheme of filtration in the domain of complex flow potential.

Figure 3
Figure 3 shows the changes in concentration field along the streamline ψ = 0.5 at a fixed time t = 5 depending on the value of the variable ϕ for the classical convective diffusion model (curve 1), the model with the standard Caputo derivative D (α) t C (curve 2, α = 0.7), and the model with the k-Caputo derivative k D α t C (curves 3 and 4) for various values of the parameter k and a fixed value of order α = 0.7 (3 -k = 0.8, 4 -k = 0.7).

Figure 4
Figure 4 shows the growth pattern of the solute's concentration at a fixed point of the filtration domain (ϕ = 2.44, ψ = 0.5) depending on the value of the dimensionless time t: curve 1 corresponds to the classical model of convective diffusion, curve 2 corresponds to the model with the standard Caputo derivative D (α) t C for α = 0.9, curves 3-5 correspond to the model with the k-Caputo derivative (3 -k = 0.9, 4 -k = 0.8, 5 -k = 0.7; 3, 4, 5α = 0.9).