Simpliﬁed Mathematical Model for the Description of Anomalous Migration of Soluble Substances in Vertical Filtration Flow

: Since the use of the fractional-differential mathematical model of anomalous geomigration process based on the MIM (mobile–immoble media) approach in engineering practice signiﬁcantly complicates simulations, a corresponding simpliﬁed mathematical model is constructed. For this model, we state a two-dimensional initial-boundary value problem of convective diffusion of soluble substances under the conditions of vertical steady-state ﬁltration of groundwater with free surface from a reservoir to a coastal drain. To simplify the domain of simulation, we use the technique of transition into the domain of complex ﬂow potential through a conformal mapping. In the case of averaging ﬁltration velocity over the domain of complex ﬂow potential, an analytical solution of the considered problem is obtained. In the general case of a variable ﬁltration velocity, an algorithm has been developed to obtain numerical solutions. The results of process simulation using the presented algorithm shows that the constructed mathematical model can be efﬁciently used to simplify and accelerate modeling process.


Introduction
We study the problem of mathematical modeling of anomalous dynamics of soluble substances' convective diffusion under the conditions of vertical steady-state filtration of groundwater. Such problems arise while solving many problems related to the protection of water resources from pollution by industrial and domestic wastewater, as well as due to the need for desalination, leaching of soils during land reclamation, etc. It should be noted that an extensive amount of literature is devoted to these problems considering them within the framework of classical models of mass transfer (e.g., [1][2][3][4][5]). Empirical studies of these processes along with their investigation using differential models are presented in [6][7][8]. Soft computing techniques are also but rarely used for their analysis and prediction [9]. Regarding the issues of modeling anomalous migration of soluble substances in underground filtration flows, several problems have been posed and solved using the fractional-differential approach ( [10][11][12][13][14]) that allows, in particular, taking into account the effects of memory and spatial correlations [15][16][17].
We consider the mathematical model of geomigration process based on the MIM (mobile-immobile media) approach [18][19][20]. If we take into account immobilization (entrapment of particles by soil's skeleton or their penetration into the volume of a bound fluid), then the equation of convective diffusion takes the form [18][19][20] ∂C ∂t + β ∂C I ∂t where C, C I are the volume concentrations of particles in mobile and immobile phases, β = θ I /θ, θ, θ I are the porosities in mobile and immobile zones, L(C) = d∆C − υ∇C, υ is the filtration velocity, d is the convective diffusion coefficient, ∇ is the Hamilton operator, and ∆ is the Laplace operator with respect to the geometric variables. According to Shumer et al. [20], the dynamics of particles' outflow in the immobile phase can be described by the equation where C D γ t is the operator of Caputo-Gerasimov fractional derivative of the order γ(0 < γ < 1) with respect to the variable t [21,22]. Then, considering Equation (2), Equation (1) of the MIM model takes the form [18][19][20] ∂C ∂t with the total concentration calculated as C tot = θC + θ I C I . Considering that medium's memory-related properties manifest themselves in the process of migration in mobile phase, Equation (3) takes the form where C D α t is the operator of Caputo-Gerasimov fractional derivative of the order α(0 < α < 1) with respect to the variable t.
Let us note, that Equation (4) is formally similar to the equation of water exchange in mobile-immobile zones of swelling soils obtained in [23].
The use of the mathematical model of anomalous migration process that is based on the integro-differential Equation (4) in engineering practice with the aim of developing reliable methods for predicting process dynamics is associated with mathematical difficulties and complicates the simulation process. In some cases, one of the effective approaches to accelerate and simplify estimative engineering calculations when studying migration dynamics under the conditions of geofiltration consists in the simplification of the original mathematical model by an appropriate approximation of fractional derivatives and transition to a new, simplified mathematical model based on the classical convective diffusion equation. Thus, approximating the fractional derivatives in Equation (4) using the relation [24] C D α t u(t) ≈ αu (t) we obtain from Equation (4) the equation of a simplified mathematical model of the considered process of convective diffusion with particles' immobilization in the form where Further, within the framework of the simplified mathematical model based on Equation (6), we state a two-dimensional initial-boundary value problem of convective diffusion of soluble substances under the conditions of vertical steady-state filtration of groundwater with a free surface from a reservoir into a coastal drain. In the case of averaged filtration velocity, we obtain an analytical solution of the considered problem. In the general case of variable filtration velocity, we present an algorithm for obtaining its numerical solution and the results of simulations.

Initial-Boundary Value Problem
We consider the problem of convective diffusion modeling with vertical steady-state filtration of groundwater from a reservoir to a buried semi-infinite drainage canal.
The filtration scheme corresponding to this problem has been studied quite well. For example, Lavryk et al. [5] provided a closed-form solution of the initial-boundary value problem obtained by conformal mapping of the filtration domain G z onto the domain G ω in the plane of complex flow potential ω = ϕ + iψ (where ϕ is the potential of filtration velocity and ψ is the flow function). As a result of filtration problem's solution, the characteristic flow function z = f (ω) was found and flow velocity field was determined. The corresponding relations are given in [5].
Within the framework of the original mathematical model based on Equation (4), the modeling of fractional-differential dynamics of convective diffusion in accordance with the considered filtration scheme can be performed finding in the domain G z × (0, ∞) the solutions of the equation with the following initial and boundary conditions: where C 1 is the given concentration of soluble substances on the filtration inflow y = 0, C 0 (x, y) is the given initial distribution function, n is the outer normal to the corresponding curve, γ 1 , γ 2 are the depression curves, H is the depth of the drainage canal, and υ x , υ y are the projection of the filtration velocity vector on the axes Ox and Oy.
Since the filtration domain G z is a non-canonical domain of complex configuration, an efficient way to solve initial-boundary value problems of this type is to make transition to new variables (ϕ, ψ)-the points of geometrically simpler domain of complex flow potential, which is [5] for the considered problem the rectangle G ω with sides ϕ 0 , Q (ϕ 0 = κH where κ is the soil's filtration coefficient and Q is the full filtration rate). Then, the problem in Equations (7)-(9) can be formulated in the domain of complex flow potential G ω in the form where Introducing the following variables and parameters: (υ 0 is the characteristic parameter of velocity), we rewrite the problem in Equations (10)- (12) in the form (the prime sign is further omitted to simplify the formulas) Therefore, the modeling of migration process dynamics using this fractional differential model reduces to solving the initial-boundary value problem in Equations (15)-(17) in the domain G ω × (0, +∞) and the following transition from the domain G ω = {(ϕ, ψ) : 0 < ϕ < ϕ 0 , 0 < ψ < 1} to the physical domain G z according to the solution of the corresponding filtration problem presented in [5]. Within the framework of the simplified mathematical model based on Equation (6), the modeling reduces to solving the initial-boundary value problem for the convective diffusion equation (16)- (17) are met.

Algorithm for Numerical Solution of the Initial-Boundary Value Problem in the Case of a Variable Filtration Velocity
In the case whenυ 2 (ϕ, ψ) = const, the approximated solution of the initial-boundary value problem can be found on the base of finite-difference approach as follows. Introducing the grid domain where h 1 , h 2 , τ are the grid steps with respect to the geometric variables ϕ, ψ and the time t, correspondingly, we associate with Equation (39) the following analog of the locally one-dimensional [26] finite-difference scheme of A.A. Samarskii: , Cφ ϕ , Cψ ψ are central and second finite-difference derivatives with respect to ϕ and ψ, correspondingly. Unwinding in Equation (43) the finite-difference operators taking introduced notations into consideration and collecting similar terms, we obtain on the half-integer time step t j+ 1 2 the following equations' system: where On the integer time step t j+1 from Equation (44) we get Let us note that the difference equations of the systems in Equations (45) and (48) are three-diagonal and, using the difference analogs of the boundary conditions in Equations (40) and (41), are effectively solved by the Thomas algorithm [26,27]. The stability of the Thomas algorithm for Equations (45) and (48) follows from the fact of diagonal prevalence in the coefficient matrices of these linear algebraic equations systems.
The brief description of the main stages of the simulation algorithm for solving the considered initial-boundary value problem based on the finite-difference technique using the locally one-dimensional scheme is as follows.
At the first stage, based on the analytical dependencies from Lavryk et al. [5], the elements of the array of velocity field values that correspond to the considered filtration scheme are calculated. Having a filtration velocity field, at the second stage, the systems of algebraic Equations (45)

The Results Of Simulations
Numerical modeling of the dynamics of the considered migration process using the original MIM model and the proposed simplified mathematical model, as well as a comparative analysis of simulation results according to both of these models, were performed with respect to dimensionless variables determined by the relations in Equation (14) under the condition C 0 (ϕ, ψ) ≡ 0.
Some obtained results for κ = 0.57 m/day, H = 0.45 m, Q = 6 m 2 /day, d = 0.1 m 2 /day, C 1 = ϕ 0 = 1, τ = 0.0005, and n = m = 30 are depicted in Figures 1-3.  Figure 1 shows the curves of concentration field distribution along the streamline ψ = 0.5 at a fixed time t = 4.2 h that, according to Equation (14), corresponds to the values t = 0.03 for α = 1, t = 0.212 for α = 0.84, and t = 0.0256 for α = 0.92. For the simulations, the results of which are shown in Figure 1, we used such values of fractional derivative exponents that yield accurate enough approximation of the initial model by the simplified one. Figure 2 shows the curves of the relative between the solutions according to the initial MIM model and the simplified model (here, C is the concentration that corresponds to the MIM model, C S is the concentration that corresponds to the simplified model, and · is the norm in L 2 ) depending on the value of the parameter γ for various fixed values of the parameter α. Similar curves for the E (α,γ) RE subject to the values of the parameter α for the fixed values of γ are presented in Figure 3.  The results of simulation using the simplified mathematical model are in good qualitative agreement with the corresponding results obtained using the original fractional-differential MIM model ( Figure 1).

2.
Simulation accuracy when the simplified model is used in comparison with the original MIM model is in many cases (e.g., when performing evaluative calculations in engineering practice) satisfactory. This can be illustrated by the following example: the inequality E The obtained data indicate that in some cases the simplified mathematical model of the geomigration process can be effectively used for carrying out evaluative calculations in engineering practice while developing constructive solutions in the field of designing environmentally hazardous engineering objects.

Conclusions
The use of a fractional-differential mathematical model of an anomalous geomigration process under the conditions of filtration taking into account the MIM approach in engineering practice presents certain severe difficulties associated with the complexity of simulation process. One of the effective approaches to simplify and accelerate the process of evaluative engineering calculations in this case is an approach based on the simplification of the original mathematical model and the transition to a simpler model based on the classical convective diffusion equation. We constructed such a simplified model and within its framework formulated the two-dimensional initial-boundary value problem of convective diffusion of soluble substances under the conditions of vertical steady-state groundwater filtration with a free surface from a reservoir into a coastal drain. In the case of averaged filtration velocity, an analytical solution of this problem is obtained. In the general case of a variable filtration velocity, an algorithm has been developed for obtaining its numerical solution. The results of simulations based on this algorithm indicate that the results of concentration fields modeling within the framework of the simplified mathematical model are in qualitative agreement with the corresponding results obtained using the original MIM model. The constructed model can be effectively used in carrying out evaluative computations in engineering practice for simplification and acceleration of simulation process. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.