Anomalous Solute Transport in a Cylindrical Two-Zone Medium with Fractal Structure

Abstract

In this paper, a problem of anomalous solute transport in a coaxial cylindrical two-zone porous medium with fractal structure is posed and numerically solved. The porous medium is studied in the form of cylinder with two parts: macropore—with high permeability characteristics in the central part and micropore—with low permeability around it. Anomalous solute transport is modeled by differential equations with a fractional derivative. The solute concentration and pressure fields are determined. Based on numerical results, the influence of the fractional derivatives order on the solute transport process is analysed. It was shown that with a decrease in the order of the derivatives in the diffusion term of the transport equation in the macropore leads to a “fast diffusion” in both zones. Characteristics of the solute transport in both zones mainly depend on the concentration distribution and other hydrodynamic parameters in the macropore.

1. Introduction

Problems of solute transport and filtration of nonhomogeneous fluids are of great practical importance in many branches of engineering and technology. Oil reservoirs and aquifers may contain areas with low filtration characteristics where fluid filtration and solute transport occurs with arising of several specifics. Many works have been published where on the basis of experimental studies, several novel effects are revealed. In theoretical works, research is based mainly on phenomenological mathematical models.

Conceptually, models can be divided into two large groups. In the first group of models, the process is described from a microscopic point of view. The solute transport is considered in a single pore or channel with a certain geometry or in a hollow media between aggregates of a certain type. The transfer from macropores to the environment is described by the diffusion equation. Such types of models were analyzed, in particular, in [1,2,3,4,5,6,7,8]. In the second group of models, the geometry of macropores and their environment is not considered explicitly, but instead, channels of various sizes and surrounding rocks are considered to be whole and are studied from a macroscopic point of view. The medium is divided into two parts, in one of which the fluid is considered to be active, i.e., mobile, and in the other part—immobile or inactive. Mass transport between two parts (or zones) is usually described by a first order kinetic equation. Models of this type are commonly referred to as “mobile-immobile” type models. As one of the first models of this type, one can indicate the work [9], where the concepts of mobile and immobile parts (zones) of the environment were introduced. This approach was further developed in [10,11,12,13,14,15,16,17].

In [9], a porous medium was divided into two parts (zones): with a mobile and immobile fluid. The diffusion flow between the zones is considered proportional to the difference in concentration in the zones. Sorption processes in both zones were considered. Sorption was considered to be in equilibrium with a linear isotherm. The model, presented in the work, well describes the well-known phenomenon of “tailing”, which is characteristic for solute transport in aggregated media. In addition, an earlier manifestation of breakthrough curves has been established.

In a more general form, the issues of solute transport using the second group of models were analyzed in [18]. The porous medium was divided into two parts (zones): with a mobile and immobile fluid. The diffusion flow between the zones was considered proportional to the difference in concentration in the zones.

In [18], three cases were identified that can lead to “tailing”.

(1) Unsaturated conditions. At the same time, not all the void space of the medium is saturated with fluid, part of the space remains saturated with air. With a decrease in the water content of the media, significant tailing is observed. Part of the pores is excluded from the transport process and, accordingly, part of the fluid becomes immobile. This part of the media is called the “dead” zone or zone with a immobile fluid.

(2) Aggregated media. Such media consist of well and poorly permeable zones, i.e., parts. The poorly permeable zone contains many micropores, where diffusion is the main transport mechanism of solute, while convective transport is insignificant and can be neglected. This leads to the restriction of fluid mixing and, consequently, the appearance of tailing, even in the case of complete saturation of the media. In fairly large aggregates with an immobile fluid, the diffusion path increases, which leads to tailing.

(3) Filtration velocity. The results of some experimental studies show that tailing becomes more noticeable at low filtration velocity.

In [19], theoretical calculations obtained in [20] were compared with the results of experimental studies on the transport of tritium in an unsaturated sorbing porous medium. The experiments were carried out on the displacement of tritium through a 30 cm column with loam and the parameters of the model were evaluated. The experimental results show that tritium adsorption and ion exchange occur during movement through the medium. Ratio immobile water increases with a decrease in filtration rate and an increase in aggregate size, ranging from 6 to 45% of the total volume of water in the sample (medium). It was shown that the analytical solution [20] satisfactorily describes the experimental data, the effects of tailing can be well explained by the diffusion exchange of tritium between zones with a mobile and immobile fluid.

The “mobile-immobile” approach was used to describe various experimental data. In [21], the usual convective diffusion equation, as well as the equations of the “mobile-immobile” (MIM) approach, were used to describe the results of solute transport in a model of a homogeneous porous media filled with a sand in a fracture with smooth and non-smooth walls. The results show that the MIM model better describes the experimental breakthrough curves, especially their peak values and tailing. This model describes well the results for a coarse-wall fracture than a smooth-wall fracture.

The work [22] was also devoted to the experimental study of the solute transport in a laboratory sample of a porous media. A constant linearly dependent on distance, exponentially varying (distance) dispersion coefficient was used in the framework of the MIM model. It is shown that large values of the mass transport coefficient decrease the current concentration values in the mobile zone, which means a large mass transport from the mobile to the immobile zone of the medium.

In [23] experimentally investigated the solute transport in loam with simultaneous diffusion into the surrounding matrix. Organic hydrophobic substances such as phenanthrene, carbofuran, 1 and 2 dichlorobenzene, trichloroethane were used as substances. The results of experimental studies are compared with the calculated results based on the solution of known and new problems. The model took into account advective transport in macropore, diffusion transport into the surrounding micropore, as well as linear sorption in both zones.

The rocks of many oil fields, as a rule, are heterogeneous both on a microscopic and macroscopic scale. A typical example of nonhomogeneous porous media are aggregated and fractured-porous media.

The complex trajectory of fluid particles and substances in the inter-aggregate media, fractures and porous blocks causes anomalous transport, so that conventional convective transport equations cannot adequately describe solute transport, transport equations must take this anomalous into account. Such media can be considered to be fractals [24,25].

In fractals, one of the first, transport equations were proposed in [26]. In fractured-porous media, the transport equations were analyzed in [27,28,29]. It is shown that the order of the fractional derivative in the equations depends on the fractal dimension of the media.

It is known that for normal Fick diffusion, solute transport is described by the diffusion equation [30], which is invariant if the distance scales x and time t correlated as <$x^2$>∼t. Numerous experimental studies show that this relation is violated in fractal porous media. The mean square distance of a particle in a random walk obeys the relation <$x^2$>∼$t^{2/(2+\theta )}$, where $\theta$—index of anomalous diffusion. Therefore, the scales for the diffusion equation must be correlated as $x_0^2$∼$t_0^{2/(2+\theta )}$. It is clear that when $\theta =0$ we obtain the classical Fick’s law. The anomalous diffusion index is determined by the fractal dimension of the media $d_f$. For example, for the Koch curve we have $\theta =2d_f-2$, where $d_f=ln4/ln3$. Under these conditions, the equations of solute transport contain fractional derivatives both in time and in spatial coordinates.

One of the problems that arise when using fractional derivatives is that there is no unambiguous definition of them. Numerical methods for solving problems for equations with fractional derivatives depend on the selected derivative, so there is a need to analyze and compare the results obtained using different definitions and numerical methods [31].

In [32], finite-difference approximation was proposed for the space fractional convective-diffusion model, which has the coefficients of space variable on the given bounded domain over time and space. It was shown that the Crank-Nicolson difference scheme based on the right shifted Grünwald-Letnikov difference formula is unconditionally stable and also has second-order consistency both in temporal and spatial terms with extrapolation to the limit approach.

An experimental study of the solute transport in a cylindrical zone with a central cylindrical macropore was the subject of work [33]. The influence of the fluid velocity on the transport characteristics is studied. In the experiments, tritium was used as a portable solute. The fluid velocity varied by two orders of magnitude, the feed had the character of a “short” and “long” pulse length. A review of laboratory and field experiments, where mass transfer rates are modelled with the first order mass transfer and the Fickian pore-diffusion models, is given in [34]. It was shown that one of the causes non-adequate description of the experimental data with the Fickian pore-diffusion model can be alternative mass transfer mechanisms, i.e., advection-diffusion processes then diffusion only. Therefore, in general, transport mechanism in both macro and micropores is advection-diffusion.

In this work, we consider an anomalous solute transport problem in a two-zone coaxial cylindrical porous medium. Intrinsic cylinder is treated as macropore with comparatively well permeability and other solute transport characteristics, and external cylindrical zone is treated as a micropore with comparatively small permeability and solute transport characteristics. Fist we present a mathematical model, i.e., fractional differential equations and corresponding initial and boundary conditions. Then a numerical procedure of the solution of posed problem is presented. After that some results are presented in graphical form and analysed. Influence of anomalous effects on solute transport characteristics is shown. Note, that in the work in distinguish with two-zones mobile-immobile solute transport problem we consider both zones as permeable media. So we have two zones with mobile liquid, but one of them is low permeable for the liquid. It considerably alters all filtration and solute transport characteristics.

2. Statement of Problem

A porous medium considered consists of two parts: (1) macroporous cylindrical medium (macropore) with a radius a (i.e., region ${\mathsf{\Omega }}_1\left\{0\le x<\infty ,0\le r\le a\right\}$), with large pores, characterized by a relatively high porosity and average fluid velocity in it, (2) surrounding cylindrical microporous medium (micropore) occupying the region ${\mathsf{\Omega }}_2\left\{0\le x<\infty ,a\le r\le b\right\}$, with low porosity (here and futher 1 stands for the macropore and 2 for the micropore) (Figure 1).

Macropore is considered to be a two-dimensional medium, i.e., pressure distributions and filtration velocity in this zone will be heterogeneous. To determine them, it is necessary to derive corresponding equations and solve together with the transport equations.

Assume that the outer cylindrical region ${\mathsf{\Omega }}_2$ has permeability $k_2$, and internal ${\mathsf{\Omega }}_1-k_1$, where $k_2<<k_1$. External surface of the cylindrical porous media ${\mathsf{\Omega }}_2$ is impermeable. To determine solute concentration fields, it is necessary to determine the distribution of pressure and filtration velocity in cylindrical regions.

Let the point $(0,0)$ is the inlet point for fluid. Inhomogeneous fluid inflows with constant pressure $p_c=$ const. Initially there was a constant pressure in the medium $p_0$, $p_c>p_0$.

Components of the filtration velocity in ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$ defined as

$$v_{1x}=-\frac{k_1}{\mu }\frac{\partial p_1}{\partial x},v_{1r}=-\frac{k_1}{\mu }\frac{\partial p_1}{\partial r},\left(r,x\right)\in {\mathsf{\Omega }}_1,$$

(1)

$$v_{2x}=-\frac{k_2}{\mu }\frac{\partial p_2}{\partial x},v_{2r}=-\frac{k_2}{\mu }\frac{\partial p_2}{\partial r},\left(r,x\right)\in {\mathsf{\Omega }}_2,$$

(2)

where $v_{1x},v_{2x}$—components of the filtration velocity in the x direction in ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$, $v_{1r},v_{2r}$—radial components of the filtration velocity in ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$, respectively, $p_1,p_2$—pressure in zones ${\mathsf{\Omega }}_1$, ${\mathsf{\Omega }}_2$, $\mu$—viscosity coefficient of the fluid.

To determine pressures in zones ${\mathsf{\Omega }}_1$, ${\mathsf{\Omega }}_2$ we use piezoconductivity equations

$$\frac{\partial p_1}{\partial t}={\chi }_1\left\{\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial p_1}{\partial r}\right)+\frac{{\partial }^2{p}_1}{\partial x^2}\right\},\left(r,x\right)\in {\mathsf{\Omega }}_1,$$

(3)

$$\frac{\partial p_2}{\partial t}={\chi }_2\left\{\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial p_2}{\partial r}\right)+\frac{{\partial }^2{p}_2}{\partial x^2}\right\},\left(r,x\right)\in {\mathsf{\Omega }}_2,$$

(4)

$${\chi }_1=\frac{k_1}{\mu {\beta }^{*}},{\chi }_2=\frac{k_2}{\mu {\beta }^{*}},$$

where ${\chi }_1,{\chi }_2$—coefficients of piezoconductivity, ${\beta }^{*}$—coefficient of elastic capacity of the media (${\beta }^{*}=m{\beta }_f+{\beta }_r$, ${\beta }_f$—elasticity coefficient of the fluid and ${\beta }_r$—of the rock of porous media). For the sake of simplicity here we assume that rocks of both media ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$ have same coefficient of elasticity—${\beta }_r$.

To describe the solute transport, we use the following equations:

$${\theta }_1\frac{\partial c_1}{\partial t}+\frac{\partial \left(v_{1x}{c}_1\right)}{\partial x}+\frac{\partial \left(v_{1r}{c}_1\right)}{\partial r}={\theta }_1\left\{D_{1r}\frac{1}{r^{{\gamma }_1}}\frac{{\partial }^{{\gamma }_1}}{\partial r^{{\gamma }_1}}\left(r^{{\gamma }_1}\frac{{\partial }^{{\gamma }_1}{c}_1}{\partial r^{{\gamma }_1}}\right)+D_{1x}\frac{{\partial }^{{\gamma }_1^{*}}{c}_1}{\partial x^{{\gamma }_1^{*}}}\right\},\left(r,x\right)\in {\mathsf{\Omega }}_1,$$

(5)

$${\theta }_2\frac{\partial c_2}{\partial t}+\frac{\partial \left(v_{2x}{c}_2\right)}{\partial x}+\frac{\partial \left(v_{2r}{c}_2\right)}{\partial r}={\theta }_2\left\{D_{2r}\frac{1}{r^{{\gamma }_2}}\frac{{\partial }^{{\gamma }_2}}{\partial r^{{\gamma }_2}}\left(r^{{\gamma }_2}\frac{{\partial }^{{\gamma }_2}{c}_2}{\partial r^{{\gamma }_2}}\right)+D_{2x}\frac{{\partial }^{{\gamma }_2^{*}}{c}_2}{\partial x^{{\gamma }_2^{*}}}\right\},\left(r,x\right)\in {\mathsf{\Omega }}_2,$$

(6)

where ${\theta }_1$ and ${\theta }_2$—porosity; $c_1,c_2$—volume concentrations; $D_1,D_2$—diffusion coefficients in ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$, respectively, ${\gamma }_1,{\gamma }_1^{*},{\gamma }_2,{\gamma }_2^{*}$—orders of derivatives, t— time.

Initial and boundary conditions for Equations (3) and (4) have the form:

$$p_1\left(0,r,x\right)=p_2\left(0,r,x\right)=p_0,p_0=\mathrm{const},$$

(7)

$$p_1\left(t,0,0\right)=p_c,p_c=\mathrm{const},p_c>p_0,$$

(8)

$$\frac{\partial p_1}{\partial x}\left(t,r,0\right)=0,0<r\le a,$$

(9)

$$\frac{\partial p_1}{\partial x}\left(t,r,\infty \right)=0,0<r\le a,$$

(10)

$$\frac{\partial p_1}{\partial r}\left(t,0,x\right)=0,0<x<\infty ,$$

(11)

$$\frac{\partial p_2}{\partial x}\left(t,r,0\right)=0,a<r\le b,$$

(12)

$$\frac{\partial p_2}{\partial x}\left(t,r,\infty \right)=0,a<r\le b,$$

(13)

$$\frac{\partial p_2}{\partial r}\left(t,b,x\right)=0,0<x<\infty ,$$

(14)

$$k_1\frac{\partial p_1}{\partial r}\left(t,a,x\right)=k_2\frac{\partial p_2}{\partial r}\left(t,a,x\right),0\le x<\infty ,$$

(15)

$$p_1\left(t,a,x\right)=p_2\left(t,a,x\right),0\le x<\infty .$$

(16)

Initial and boundary conditions for (5)–(6) have the form:

$$c_1\left(0,r,x\right)=c_2\left(0,r,x\right)=0,$$

(17)

$$c_1\left(t,0,0\right)=c_0,c_0=\mathrm{const},$$

(18)

$$\frac{\partial c_1}{\partial x}\left(t,r,0\right)=0,0<r\le a,$$

(19)

$$\frac{\partial c_1}{\partial x}\left(t,r,\infty \right)=0,0<r\le a,$$

(20)

$$\frac{\partial c_1}{\partial r}\left(t,0,x\right)=0,0<x<\infty ,$$

(21)

$$\frac{\partial c_2}{\partial x}\left(t,r,0\right)=0,a<r\le b,$$

(22)

$$\frac{\partial c_2}{\partial x}\left(t,r,\infty \right)=0,a<r\le b,$$

(23)

$$\frac{\partial c_2}{\partial r}\left(t,b,x\right)=0,0<x<\infty ,$$

(24)

$$D_{1r}\frac{{\partial }^{{\gamma }_1}{c}_1}{\partial r^{{\gamma }_1}}\left(t,a,x\right)=D_{2r}\frac{{\partial }^{{\gamma }_2}{c}_2}{\partial r^{{\gamma }_1}}\left(t,a,x\right),0\le x<\infty ,$$

(25)

$$c_1\left(t,a,x\right)=c_2\left(t,a,x\right),0\le x<\infty .$$

(26)

The system of Equations (1)–(4) with initial and boundary conditions (7)–(16) allows determining the pressure and filtration velocity distributions, and system (5)–(6) with conditions (17)–(26)—solute concentration in both zones.

3. Numerical Solution of the Problem

To solve the problem (1)–(26), we apply the finite difference method [35,36,37]. For fractional derivatives we use Coputo’s defination [37]. In the region ${\mathsf{\Omega }}_1\bigcup {\mathsf{\Omega }}_2$ we introduce the grid

$${\mathsf{\Omega }}_{\tau h_1{h}_2}^{}=\left\{\left(t_k,x_i,r_j\right),t_k=\tau k,\tau =\frac{T}{K},k=\overline{0,K},x_i=i{h}_1,h_1=\frac{L}{I},i=\overline{1,I},\right.$$

$$\left.r_j=j{h}_2,h_2=\frac{a}{N},j=\overline{1,N},r_j=j{h}_2,h_2=\frac{b-a}{J-N},j=\overline{N+1,J}\right\},$$

(27)

where $\tau$—grid step in time, $h_1$—grid step in x direction, $h_2$—grid step in r direction, T—maximum time, during which the process is investigated, N—number of intervals along the radius in the macropore, L—length of cylinder, K—number of grid intervals in time, I—number of intervals along the length, J—number of intervals along the radius.

Approximating the Equations (1)–(4) in the grid zone ${\mathsf{\Omega }}_{\tau h_1{h}_2}$ with finite differences, we get

$${(v_{1x})}_{i,j}=-\frac{k_1}{\mu }\frac{{\left(p_1\right)}_{i+1,j}-{\left(p_1\right)}_{i,j}}{h_1},i=\overline{0,I-1},j=\overline{0,N},$$

(28)

$${(v_{1r})}_{i,j}=-\frac{k_1}{\mu }\frac{{\left(p_1\right)}_{i,j+1}-{\left(p_1\right)}_{i,j}}{h_2},i=\overline{0,I},j=\overline{0,N-1},$$

(29)

$${(v_{2x})}_{i,j}=-\frac{k_2}{\mu }\frac{{\left(p_2\right)}_{i+1,j}-{\left(p_2\right)}_{i,j}}{h_1},i=\overline{0,I-1},j=\overline{N+1,J},$$

(30)

$${(v_{2r})}_{i,j}=-\frac{k_2}{\mu }\frac{{\left(p_2\right)}_{i,j+1}-{\left(p_2\right)}_{i,j}}{h_2},i=\overline{0,I},j=\overline{N+1,J-1},$$

(31)

$$\frac{{\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(p_1\right)}_{i,j}^k}{0,5\tau }={\chi }_1\left\{{\mathsf{\Lambda }}_2{\left(p_1\right)}_{i,j}^k+\frac{1}{j{h}_2}\frac{{\left(p_1\right)}_{i,j}^k-{\left(p_1\right)}_{i,j-1}^k}{h_2}+{\mathsf{\Lambda }}_1{\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right\},$$

(32)

$$i=\overline{1,I-1},j=\overline{1,N-1},$$

$$\frac{{\left(p_1\right)}_{i,j}^{k+1}-{\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{0,5\tau }={\chi }_1\left\{{\mathsf{\Lambda }}_2{\left(p_1\right)}_{i,j}^{k+1}+\frac{1}{j{h}_2}\frac{{\left(p_1\right)}_{i,j}^{k+1}-{\left(p_1\right)}_{i,j-1}^{k+1}}{h_2}+{\mathsf{\Lambda }}_1{\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right\},$$

(33)

$$i=\overline{1,I-1},j=\overline{1,N-1},$$

$$\frac{{\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(p_2\right)}_{i,j}^k}{0,5\tau }={\chi }_2\left\{{\mathsf{\Lambda }}_2{\left(p_2\right)}_{i,j}^k+\frac{1}{j{h}_2}\frac{{\left(p_2\right)}_{i,j}^k-{\left(p_2\right)}_{i,j-1}^k}{h_2}+{\mathsf{\Lambda }}_1{\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right\},$$

(34)

$$i=\overline{1,I-1},j=\overline{N+1,J-1},$$

$$\frac{{\left(p_2\right)}_{i,j}^{k+1}-{\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{0,5\tau }={\chi }_2\left\{{\mathsf{\Lambda }}_2{\left(p_2\right)}_{i,j}^{k+1}+\frac{1}{j{h}_2}\frac{{\left(p_2\right)}_{i,j}^{k+1}-{\left(p_2\right)}_{i,j-1}^{k+1}}{h_2}+{\mathsf{\Lambda }}_1{\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right\},$$

(35)

$$i=\overline{1,I-1},j=\overline{N+1,J-1}.$$

where ${\mathsf{\Lambda }}_1{p}_{i,j}^{k+1}=\frac{p_{i-1,j}^{k+1}-2p_{i,j}^{k+1}+p_{i+1,j}^{k+1}}{h_1^2}$, ${\mathsf{\Lambda }}_2{p}_{i,j}^{k+1}=\frac{p_{i,j-1}^{k+1}-2p_{i,j}^{k+1}+p_{i,j+1}^{k+1}}{h_2^2}$.

Initial and boundary conditions (7)–(16) approximated as

$${\left(p_1\right)}_{i,j}^k=p_0,{\left(p_2\right)}_{i,j}^k=p_0,i=\overline{0,I},j=\overline{0,J},k=0,p_0>0,$$

(36)

$${\left(p_1\right)}_{i,j}^k=p_c,i=0,j=0,k=\overline{0,K},$$

(37)

$$\frac{{\left(p_1\right)}_{i,j}^k-{\left(p_1\right)}_{i-1,j}^k}{h_1}=0,i=1,j=\overline{0,N},k=\overline{0,K},$$

(38)

$$\frac{{\left(p_1\right)}_{i,j}^k-{\left(p_1\right)}_{i-1,j}^k}{h_1}=0,i=I,j=\overline{0,N},k=\overline{0,K},$$

(39)

$$\frac{{\left(p_1\right)}_{i,j}^k-{\left(p_1\right)}_{i,j-1}^k}{h_2}=0,i=\overline{0,I},j=1,k=\overline{0,K},$$

(40)

$$\frac{{\left(p_2\right)}_{i,j}^k-{\left(p_2\right)}_{i-1,j}^k}{h_1}=0,i=1,j=\overline{N,J},k=\overline{0,K},$$

(41)

$$\frac{{\left(p_2\right)}_{i,j}^k-{\left(p_2\right)}_{i-1,j}^k}{h_1}=0,i=I,j=\overline{N,J},k=\overline{0,K},$$

(42)

$$\frac{{\left(p_2\right)}_{i,j}^k-{\left(p_2\right)}_{i,j-1}^k}{h_2}=0,i=\overline{0,I},j=J,k=\overline{0,K},$$

(43)

$$k_1\frac{{\left(p_1\right)}_{i,j}^k-{\left(p_1\right)}_{i,j-1}^k}{h_2}=k_2\frac{{\left(p_2\right)}_{i,j+1}^k-{\left(p_2\right)}_{i,j}^k}{h_2},i=\overline{0,I},j=N,k=\overline{0,K},$$

(44)

$${\left(p_1\right)}_{i,j}^k={\left(p_2\right)}_{i,j}^k,i=\overline{0,I},j=N,k=\overline{0,K}.$$

(45)

Conditions (17)–(26) we write in finite difference form

$${\left(c_1\right)}_{i,j}^k=0,{\left(c_2\right)}_{i,j}^k=0,i=\overline{0,I},j=\overline{0,J},k=0,$$

(46)

$${\left(c_1\right)}_{i,j}^k=c_0,i=0,j=0,k=\overline{0,K},$$

(47)

$$\frac{{\left(c_1\right)}_{i,j}^k-{\left(c_1\right)}_{i-1,j}^k}{h_1}=0,i=1,j=\overline{1,N},k=\overline{0,K},$$

(48)

$$\frac{{\left(c_1\right)}_{i,j}^k-{\left(c_1\right)}_{i-1,j}^k}{h_1}=0,i=I,j=\overline{0,N},k=\overline{0,K},$$

(49)

$$\frac{{\left(c_1\right)}_{i,j}^k-{\left(c_1\right)}_{i,j-1}^k}{h_2}=0,i=\overline{1,I},j=1,k=\overline{0,K},$$

(50)

$$\frac{{\left(c_2\right)}_{i,j}^k-{\left(c_2\right)}_{i-1,j}^k}{h_1}=0,i=1,j=\overline{N,J},k=\overline{0,K},$$

(51)

$$\frac{{\left(c_2\right)}_{i,j}^k-{\left(c_2\right)}_{i-1,j}^k}{h_1}=0,i=I,j=\overline{N,J},k=\overline{0,K},$$

(52)

$$\frac{{\left(c_2\right)}_{i,j}^k-{\left(c_2\right)}_{i,j-1}^k}{h_2}=0,i=\overline{0,I},j=J,k=\overline{0,K},$$

(53)

$$D_{1r}{\theta }_1\frac{{\left(c_1\right)}_{i,j}^k-{\gamma }_1{\left(c_1\right)}_{i,j-1}^k}{\mathsf{\Gamma }\left(2-{\gamma }_1\right)h_2}=D_{2r}{\theta }_2\frac{{\left(c_2\right)}_{i,j+1}^k-{\gamma }_2{\left(c_2\right)}_{i,j}^k}{\mathsf{\Gamma }\left(2-{\gamma }_2\right)h_2},$$

(54)

$$i=\overline{0,I},j=N,k=\overline{0,K},$$

$${\left(c_1\right)}_{i,j}^k={\left(c_2\right)}_{i,j}^k,i=\overline{0,I},j=N,k=\overline{0,K}.$$

(55)

where $\mathsf{\Gamma }(·)$—gamma function.

The computation sequence is as follows: first we define $p_1$,$p_2$ from the grid Equations (32)–(35), then the filtration velocity components—from (28)–(31).

Equations (32) and (33) are solved using the one-dimensional Thomas’ algorithm and are given the following form

$$A_1{\left(p_1\right)}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-B_1{\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+E_1{\left(p_1\right)}_{i+1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}=-{\left(F_1\right)}_{i,j}^k,$$

$$i=\overline{1,I-1},j=\overline{0,N-1},k=\overline{0,K-1},$$

(56)

$$A_1^{*}{\left(p_1\right)}_{i,j-1}^{k+1}-B_1^{*}{\left(p_1\right)}_{i,j}^{k+1}+E_1^{*}{\left(p_1\right)}_{i,j+1}^{k+1}=-{\left(F_1^{*}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.},$$

$$i=\overline{0,I},j=\overline{1,N-1},k=\overline{0,K-1},$$

(57)

where

$$A_1=\frac{0,5\tau {\chi }_1}{h_1^2},B_1=\frac{\tau {\chi }_1}{h_1^2}+1,E_1=\frac{0,5\tau {\chi }_1}{h_1^2},$$

$${\left(F_1\right)}_{i,j}^k={\left(p_1\right)}_{i,j}^k+\frac{0,5\tau {\chi }_1}{h_2^2}\left({\left(p_1\right)}_{i,j+1}^k-{\left(p_1\right)}_{i,j}^k\right),j=0,$$

$${\left(F_1\right)}_{i,j}^k={\left(p_1\right)}_{i,j}^k+\frac{0,5\tau {\chi }_1}{h_2^2}\left(\frac{{\left(p_1\right)}_{i,j}^k-{\left(p_1\right)}_{i,j-1}^k}{j}+{\left(p_1\right)}_{i,j-1}^k-2{\left(p_1\right)}_{i,j}^k+{\left(p_1\right)}_{i,j+1}^k\right),$$

$$j=\overline{1,N-1},$$

$$A_1^{*}=\frac{0,5\tau {\chi }_1}{h_2^2}\left(1-\frac{1}{j}\right),B_1^{*}=\frac{\tau {\chi }_1}{h_2^2}\left(1-\frac{0,5}{j}\right)+1,E_1^{*}=\frac{0,5\tau {\chi }_1}{h_2^2},$$

$${\left(F_1^{*}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}={\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+\frac{0,5\tau {\chi }_1}{h_1^2}\left({\left(p_1\right)}_{i+1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right),i=0,$$

$${\left(F_1^{*}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}={\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+\frac{0,5\tau {\chi }_1}{h_1^2}\left({\left(p_1\right)}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-2{\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+{\left(p_1\right)}_{i+1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right),i=\overline{1,I-1},$$

$${\left(F_1^{*}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}={\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+\frac{0,5\tau {\chi }_1}{h_1^2}\left({\left(p_1\right)}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right),i=I.$$

Equations (34) and (35) are reduced to the similar form

$$A_2{\left(p_2\right)}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-B_2{\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+E_2{\left(p_2\right)}_{i+1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}=-{\left(F_2\right)}_{i,j}^k,$$

$$i=\overline{1,I-1},j=\overline{N+1,J},k=\overline{0,K-1},$$

(58)

$$A_2^{*}{\left(p_2\right)}_{i,j-1}^{k+1}-B_2^{*}{\left(p_2\right)}_{i,j}^{k+1}+E_2^{*}{\left(p_2\right)}_{i,j+1}^{k+1}=-{\left(F_2^{*}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.},$$

$$i=\overline{0,I},j=\overline{N+1,J-1},k=\overline{0,K-1},$$

(59)

where

$$A_2=\frac{0,5\tau {\chi }_2}{h_1^2},B_2=\frac{\tau {\chi }_2}{h_1^2}+1,E_2=\frac{0,5\tau {\chi }_2}{h_1^2},$$

$${\left(F_2\right)}_{i,j}^k={\left(p_2\right)}_{i,j}^k+\frac{0,5\tau {\chi }_2}{h_2^2}\left(\frac{{\left(p_2\right)}_{i,j}^k-{\left(p_2\right)}_{i,j-1}^k}{j}+{\left(p_2\right)}_{i,j-1}^k-2{\left(p_2\right)}_{i,j}^k+{\left(p_2\right)}_{i,j+1}^k\right),$$

$$j=\overline{N+1,J-1},$$

$${\left(F_2\right)}_{i,j}^k={\left(p_2\right)}_{i,j}^k+\frac{0,5\tau {\chi }_2}{h_2^2}\left({\left(p_2\right)}_{i,j-1}^k-{\left(p_1\right)}_{i,j}^k\right),j=J;$$

$$A_2^{*}=\frac{0,5\tau {\chi }_2}{h_2^2}\left(1-\frac{1}{j}\right),B_2^{*}=\frac{\tau {\chi }_2}{h_2^2}\left(1-\frac{0,5}{j}\right)+1,E_2^{*}=\frac{0,5\tau {\chi }_2}{h_2^2},$$

$${\left(F_2^{*}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}={\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+\frac{0,5\tau {\chi }_2}{h_1^2}\left({\left(p_2\right)}_{i+1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right),i=0,$$

$${\left(F_2^{*}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}={\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+\frac{0,5\tau {\chi }_2}{h_1^2}\left({\left(p_2\right)}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-2{\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+{\left(p_2\right)}_{i+1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right),i=\overline{1,I-1},$$

$${\left(F_2^{*}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}={\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+\frac{0,5\tau {\chi }_2}{h_1^2}\left({\left(p_2\right)}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right),i=I.$$

Equations (56) and (58) are solved by the Thomas’ algorithm in the directions i. For this we use the following relations

$${\left(p_m\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}={\alpha }_{i+1,j}^m{\left(p_m\right)}_{i+1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+{\xi }_{i+1,j}^m,i=\overline{0,I-1},j=\overline{0,J},m=1,2,$$

(60)

where ${\alpha }_{i,j}^m,{\xi }_{i,j}^m$—coefficients of Thomas algorithm.

To solve Equations (57) and (59), Thomas algorithm method is used

$${\left(p_1\right)}_{i,j}^{k+1}={\alpha }_{i,j+1}^{}{\left(p_1\right)}_{i,j+1}^{k+1}+{\xi }_{i,j+1}^{},i=\overline{0,I},j=\overline{0,N-1},$$

(61)

$${\left(p_2\right)}_{i,j}^{k+1}={\alpha }_{i,j-1}^{*}{\left(p_2\right)}_{i,j-1}^{k+1}+{\xi }_{i,j-1}^{*},i=\overline{0,I},j=\overline{N+1,J},$$

(62)

where ${\alpha }_{i,j}^{},{\xi }_{i,j}^{},{\alpha }_{i,j}^{*},{\xi }_{i,j}^{*}$—coefficients of Thomas’ algorithm.

Substituting (60) in (56), (58) we get the following recurrent formulas for determining coefficients ${\alpha }_{i,j}^m,{\xi }_{i,j}^m$, $m=1,2$ ($m=1$ by $j=\overline{0,S}$, $m=2$ by $j=\overline{N+1,J}$ )

$${\alpha }_{i+1,j}^m=\frac{E_m}{B_m-A_m{\alpha }_{i,j}^m},{\xi }_{i+1,j}^m=\frac{{\left(F_m\right)}_{i,j}^k+A_m{\xi }_{i,j}^m}{B_m-A_m{\alpha }_{i,j}^m},$$

$$i=\overline{1,I-1},j=\overline{0,J},m=1,2.$$

(63)

To determine the Thomas’ algorithm coefficients in (61), (62), we obtain the following recurrent formulas

$${\alpha }_{i,j+1}^{}=\frac{E_3}{B_3-A_3{\alpha }_{i,j}^{}},{\xi }_{i,j+1}^{}=\frac{{\left(F_3\right)}_{i,j}^k+A_3{\xi }_{i,j}^{}}{B_3-A_3{\alpha }_{i,j}^{}},i=\overline{0,I},j=\overline{1,N-1},$$

(64)

$${\alpha }_{i,j-1}^{*}=\frac{A_4}{B_4-E_4{\alpha }_{i,j}^{*}},{\xi }_{i,j+1}^{*}=\frac{{\left(F_4\right)}_{i,j}^k+E_4{\xi }_{i,j}^{*}}{B_4-E_4{\alpha }_{i,j}^{*}},i=\overline{0,I},j=\overline{J,N+1}.$$

(65)

Starting values of Thomas’ algorithm coefficients ${\alpha }_{i,j}^m,{\xi }_{i,j}^m$, $m=1,2$ determined from boundary conditions (36), (38) and (41), i.e.

$$\left\{\begin{array}{c}{\alpha }_{i,j}^m=0,{\xi }_{i,j}^m=p_c\mathrm{at}i=1,j=0,m=1,\hfill \\ {\alpha }_{i,j}^m=\frac{\tau {\chi }_m}{\tau {\chi }_m+h_1^2},{\xi }_{i,j}^m=\frac{{\left(F_m\right)}_{i,j}^k}{B_m}\mathrm{at}i=1,j=\overline{1,J},m=1,2.\hfill \end{array}\right.$$

(66)

Starting values ${\alpha }_{i,j}^{},{\xi }_{i,j}^{},{\alpha }_{i,j}^{*},{\xi }_{i,j}^{*}$ in (64), (65) determined from the conditions (37), (40), (43)

$$\left\{\begin{array}{c}{\alpha }_{i,j}^{}=0,{\xi }_{i,j}^{}=p_c\mathrm{at}j=1,i=0,\hfill \\ {\alpha }_{i,j}^{}=\frac{\tau {\chi }_1}{\tau {\chi }_1+h_2^2},{\xi }_{i,j}^{}=\frac{{\left(F_3\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{B_3}\mathrm{at}j=1,i=\overline{1,I},\hfill \end{array}\right.$$

(67)

$${\alpha }_{i,j}^{*}=\frac{\tau {\chi }_2}{\tau {\chi }_2+h_2^2},{\xi }_{i,j}^{*}=\frac{{\left(F_4\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{B_4}atj=J-1,i=\overline{0,I}.$$

(68)

To determine the values ${\left(p_1\right)}_{i,j}^{k+1}$ and ${\left(p_2\right)}_{i,j}^{k+1}$ ($i=\overline{0,I},j=S$) from the boundary conditions (44), (45) we obtain the expression

$${\left(p_1\right)}_{i,j}^{k+1}={\left(p_2\right)}_{i,j}^{k+1}=\frac{{\xi }_{i,j}^{}+k{\xi }_{i,j}^{*}}{1-{\alpha }_{i,j}^{}-k\left({\alpha }_{i,j}^{*}-1\right)},i=\overline{0,I},j=N,$$

(69)

where $k=\raisebox{1ex}{$k_2$}\!\left/ \!\raisebox{-1ex}{$k_1$}\right.$.

The following computing sequence is set. From (63) taking into account (66), we determine the values of the Thomas’ algorithm coefficients ${\alpha }_{i,j}^m,{\xi }_{i,j}^m$ ($m=1,2$, $i=\overline{1,I},j=\overline{0,J}$). Then, taking into account (39), (42)

$${\left(p_m\right)}_{I,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}={\xi }_{I+1,j}^m=\frac{{\left(F_m\right)}_{I,j}^k+A_m{\xi }_{I,j}^m}{B_m-A_m{\alpha }_{I,j}^m},j=\overline{0,J},m=1,2$$

(70)

are calculated.

From relations (60) we determine ${\left(p_1\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ and ${\left(p_2\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ in directions i. In view of (69), by (61), (62) ${\left(p_1\right)}_{i,j}^{}$ and ${\left(p_2\right)}_{i,j}^{}$ are calculated on the $k+1$ layer.

After determining the pressure distribution using Equations (28)–(31) the filtration velocity field is determined on ${\mathsf{\Omega }}_{\tau h_1{h}_2}$.

Equation (5) is approximated as follows

$${\theta }_1\frac{c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-c_{ij}^k}{0,5\tau }+\frac{{\left(v_{1x}\right)}_{i,j}^k{c}_{ij}^k-{\left(v_{1x}\right)}_{i-1,j}^k{c}_{i-1,j}^k}{h_1}+\frac{{\left(v_{1r}\right)}_{i,j}^k{c}_{ij}^k-{\left(v_{1r}\right)}_{i,j-1}^k{c}_{i,j-1}^k}{h_2}=$$

$$={\theta }_1\left\{D_{1r}\frac{1}{{\left(j{h}_2\right)}^{{\gamma }_1}}\frac{h_2^{{\gamma }_1}{\left(j+\frac{1}{2}\right)}^{{\gamma }_1}-{\gamma }_1{h}_2^{{\gamma }_1}{\left(j-\frac{1}{2}\right)}^{{\gamma }_1}}{\mathsf{\Gamma }\left(2-{\gamma }_1\right)h_2^{{\gamma }_1}}\frac{c_{ij}^k-{\gamma }_1{c}_{i,j-1}^k}{\mathsf{\Gamma }\left(2-{\gamma }_1\right)h_2^{{\gamma }_1}}+\right.$$

$$+\frac{D_{1r}}{\mathsf{\Gamma }\left(3-2{\gamma }_1\right)h_2^{2{\gamma }_1}}\sum _{l=0}^{j-1}\left(c_{i,j-(l-1)}^k-2c_{i,j-l}^k+c_{i,j-\left(l+1\right)}^k\right)\left({\left(l+1\right)}^{2-2{\gamma }_1}-l^{2-2{\gamma }_1}\right)$$

$$\left.+\frac{D_{1x}}{\mathsf{\Gamma }\left(3-{\gamma }_1^{*}\right)h_1^{{\gamma }_1^{*}}}\sum _{l=0}^{i-1}\left(c_{i-(l-1),j}^k-2c_{i-l,j}^k+c_{i-\left(l+1\right),j}^k\right)\left({\left(l+1\right)}^{2-{\gamma }_1^{*}}-l^{2-{\gamma }_1^{*}}\right)\right\},$$

(71)

$$i=\overline{1,I-1},j=\overline{1,N-1},k=\overline{0,K-1},$$

$${\theta }_1\frac{c_{ij}^{k+1}-c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{0,5\tau }+\frac{{\left(v_{1x}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(v_{1x}\right)}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{h_1}+\frac{{\left(v_{1r}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(v_{1r}\right)}_{i,j-1}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{i,j-1}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{h_2}=$$

$$={\theta }_1\left\{D_{1r}\frac{1}{{\left(j{h}_2\right)}^{{\gamma }_1}}\frac{h_2^{{\gamma }_1}{\left(j+\frac{1}{2}\right)}^{{\gamma }_1}-{\gamma }_1{h}_2^{{\gamma }_1}{\left(j-\frac{1}{2}\right)}^{{\gamma }_1}}{\mathsf{\Gamma }\left(2-{\gamma }_1\right)h_2^{{\gamma }_1}}\frac{c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\gamma }_1{c}_{i,j-1}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{\mathsf{\Gamma }\left(2-{\gamma }_1\right)h_2^{{\gamma }_1}}+\right.$$

$$+\frac{D_{1r}}{\mathsf{\Gamma }\left(3-2{\gamma }_1\right)h_2^{2{\gamma }_1}}\sum _{l=0}^{j-1}\left(c_{i,j-(l-1)}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-2c_{i,j-l}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+c_{i,j-\left(l+1\right)}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)\left({\left(l+1\right)}^{2-2{\gamma }_1}-l^{2-2{\gamma }_1}\right)$$

$$\left.+\frac{D_{1x}}{\mathsf{\Gamma }\left(3-{\gamma }_1^{*}\right)h_1^{{\gamma }_1^{*}}}\sum _{l=0}^{i-1}\left(c_{i-(l-1),j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-2c_{i-l,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+c_{i-\left(l+1\right),j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)\left({\left(l+1\right)}^{2-{\gamma }_1^{*}}-l^{2-{\gamma }_1^{*}}\right)\right\},$$

(72)

$$i=\overline{1,I-1},j=\overline{1,N-1},k=\overline{0,K-1}.$$

Equation (6) is approximated in the form

$${\theta }_2\frac{c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-c_{ij}^k}{0,5\tau }+\frac{{\left(v_{2x}\right)}_{i,j}^k{c}_{ij}^k-{\left(v_{2x}\right)}_{i-1,j}^k{c}_{i-1,j}^k}{h_1}+\frac{{\left(v_{2r}\right)}_{i,j}^k{c}_{ij}^k-{\left(v_{2r}\right)}_{i,j-1}^k{c}_{i,j-1}^k}{h_2}=$$

$$={\theta }_2\left\{D_{2r}\frac{1}{j{h}_2}\frac{h_2^{{\gamma }_2}{\left(j+\frac{1}{2}\right)}^{{\gamma }_2}-{\gamma }_1{h}_2^{{\gamma }_2}{\left(j-\frac{1}{2}\right)}^{{\gamma }_2}}{\mathsf{\Gamma }\left(2-{\gamma }_2\right)h_2^{{\gamma }_2}}\frac{c_{ij}^k-{\gamma }_2{c}_{i,j-1}^k}{\mathsf{\Gamma }\left(2-{\gamma }_2\right)h_2^{{\gamma }_2}}+\right.$$

$$+\frac{D_{2r}}{\mathsf{\Gamma }\left(3-2{\gamma }_2\right)h_2^{2{\gamma }_2}}\sum _{l=0}^{j-1}\left(c_{i,j-(l-1)}^k-2c_{i,j-l}^k+c_{i,j-\left(l+1\right)}^k\right)\left({\left(l+1\right)}^{2-2{\gamma }_2}-l^{2-2{\gamma }_2}\right)$$

$$\left.+\frac{D_{2x}}{\mathsf{\Gamma }\left(3-{\gamma }_2^{*}\right)h_1^{{\gamma }_2^{*}}}\sum _{l=0}^{i-1}\left(c_{i-(l-1),j}^k-2c_{i-l,j}^k+c_{i-\left(l+1\right),j}^k\right)\left({\left(l+1\right)}^{2-{\gamma }_2^{*}}-l^{2-{\gamma }_2^{*}}\right)\right\},$$

(73)

$$i=\overline{1,I-1},j=\overline{N+1,J-1},k=\overline{0,K-1},$$

$${\theta }_2\frac{c_{ij}^{k+1}-c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{0,5\tau }+\frac{{\left(v_{2x}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(v_{2x}\right)}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{i-1,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{h_1}+\frac{{\left(v_{2r}\right)}_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\left(v_{2r}\right)}_{i,j-1}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{i,j-1}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{h_2}=$$

$$={\theta }_2\left\{D_{2r}\frac{1}{j{h}_2}\frac{h_2^{{\gamma }_2}{\left(j+\frac{1}{2}\right)}^{{\gamma }_2}-{\gamma }_1{h}_2^{{\gamma }_2}{\left(j-\frac{1}{2}\right)}^{{\gamma }_2}}{\mathsf{\Gamma }\left(2-{\gamma }_2\right)h_2^{{\gamma }_2}}\frac{c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-{\gamma }_2{c}_{i,j-1}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{\mathsf{\Gamma }\left(2-{\gamma }_2\right)h_2^{{\gamma }_2}}+\right.$$

$$+\frac{D_{2r}}{\mathsf{\Gamma }\left(3-2{\gamma }_2\right)h_2^{2{\gamma }_2}}\sum _{l=0}^{j-1}\left(c_{i,j-(l-1)}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-2c_{i,j-l}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+c_{i,j-\left(l+1\right)}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)\left({\left(l+1\right)}^{2-2{\gamma }_2}-l^{2-2{\gamma }_2}\right)$$

$$\left.+\frac{D_{2x}}{\mathsf{\Gamma }\left(3-{\gamma }_2^{*}\right)h_1^{{\gamma }_2^{*}}}\sum _{l=0}^{i-1}\left(c_{i-(l-1),j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-2c_{i-l,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+c_{i-\left(l+1\right),j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)\left({\left(l+1\right)}^{2-{\gamma }_2^{*}}-l^{2-{\gamma }_2^{*}}\right)\right\},$$

(74)

$$i=\overline{1,I-1},j=\overline{N+1,J-1},k=\overline{0,K-1}.$$

We introduce the following notation

$$V_{mx}^k=\frac{{\left(v_{mx}\right)}_{i,j}^k{c}_{ij}^k-{\left(v_{mx}\right)}_{i-1,j}^k{c}_{i-1,j}^k}{h_1},$$

$$V_{mr}^k=\frac{{\left(v_{mr}\right)}_{i,j}^k{c}_{ij}^k-{\left(v_{mr}\right)}_{i,j-1}^k{c}_{i,j-1}^k}{h_2},$$

$$G_{mr}^k=D_{mr}\frac{1}{j{h}_2}\frac{h_2^{{\gamma }_m}{\left(j+\frac{1}{2}\right)}^{{\gamma }_m}-{\gamma }_m{h}_2^{{\gamma }_m}{\left(j-\frac{1}{2}\right)}^{{\gamma }_m}}{\mathsf{\Gamma }\left(2-{\gamma }_m\right)h_2^{{\gamma }_m}}\frac{c_{ij}^k-{\gamma }_m{c}_{i,j-1}^k}{\mathsf{\Gamma }\left(2-{\gamma }_m\right)h_2^{{\gamma }_m}},$$

$$S_{mr}^k=\frac{D_{mr}}{\mathsf{\Gamma }\left(3-2{\gamma }_m\right)h_2^{{\gamma }_m}}\sum _{l=0}^{i-1}\left(c_{i-(l-1),j}^k-2c_{i-l,j}^k+c_{i-\left(l+1\right),j}^k\right)\left({\left(l+1\right)}^{2-2{\gamma }_m}-l^{2-2{\gamma }_m}\right),$$

$$S_{mx}^k=\frac{D_{mx}}{\mathsf{\Gamma }\left(3-{\gamma }_m^{*}\right)h_1^{{\gamma }_m^{*}}}\sum _{l=0}^{i-1}\left(c_{i-(l-1),j}^k-2c_{i-l,j}^k+c_{i-\left(l+1\right),j}^k\right)\left({\left(l+1\right)}^{2-{\gamma }_m^{*}}-l^{2-{\gamma }_m^{*}}\right),$$

$$m=1,2.$$

Given these notation, Equations (71)–(74), have the form respectively

$$c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}=c_{ij}^k-\frac{0,5\tau }{{\theta }_1}\left(V_{1x}^k+V_{1r}^k\right)+0,5\tau \left(G_{1r}^k+S_{1r}^k+S_{1x}^k\right),$$

(75)

$$i=\overline{1,I-1},j=\overline{1,N-1},k=\overline{0,K-1},$$

$$c_{ij}^{k+1}=c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-\frac{0,5\tau }{{\theta }_1}\left(V_{1x}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+V_{1r}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)+0,5\tau \left(G_{1r}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+S_{1r}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+S_{1x}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right),$$

(76)

$$i=\overline{1,I-1},j=\overline{1,N-1},k=\overline{0,K-1},$$

$$c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}=c_{ij}^k-\frac{0,5\tau }{{\theta }_2}\left(V_{2x}^k+V_{2r}^k\right)+0,5\tau \left(G_{2r}^k+S_{2r}^k+S_{2x}^k\right),$$

(77)

$$i=\overline{1,I-1},j=\overline{N+1,J-1},k=\overline{0,K-1},$$

$$c_{ij}^{k+1}=c_{ij}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-\frac{0,5\tau }{{\theta }_2}\left(V_{2x}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+V_{2r}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)+0,5\tau \left(G_{2r}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+S_{2r}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+S_{2x}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right),$$

(78)

$$i=\overline{1,I-1},j=\overline{N+1,J-1},k=\overline{0,K-1}.$$

From (54) we define $c_{i,j}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ and $c_{i,j}^{k+1}$ ($i=\overline{0,I-1},j=N$)

$$c_{i,S}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}=\frac{D{c}_{i,S+1}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+{\gamma }_1{c}_{i,S-1}^{k+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{1+D{\gamma }_2},$$

(79)

$$c_{i,S}^{k+1}=\frac{D{c}_{i,S+1}^{k+1}+{\gamma }_1{c}_{i,S-1}^{k+1}}{1+D{\gamma }_2},$$

(80)

where $D=\frac{D_{2r}\mathsf{\Gamma }\left(2-{\gamma }_1\right)h_2^{{\gamma }_1}}{D_{1r}\mathsf{\Gamma }\left(2-{\gamma }_2\right)h_2^{{\gamma }_2}}$.

The concentration field is determined from (75), (80) with conditions (46)–(53).

4. Results and Discussion

The proposed methodology serves to carry out computational experiments, where we used the following values of the initial parameters: $\mu ={10}^{-1}$ Pa·s, ${\beta }^{*}={10}^{-9}$ Pa ${}^{-1}$, $k_1={10}^{-12}÷2.5·{10}^{-12}$ m${}^2$, $k_2=5·{10}^{-13}÷·{10}^{-12}$ m${}^2$, $p_c=1000$ Pa, $p_0=0$ Pa, $T=3600$ s, $h_1=0.01$ m, $h_2=0.005$ m, $\tau =0.1$ s, $a=0.5$ m, $b=1.5$ m.

In Figure 2 and Figure 3 we present surfaces $c_1,c_2$ and p for the classic case, i.e., ignoring anomalous effects. It is easy to notice the spread of surfaces along the directions x and r. On the common border of the two zones, ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$, i.e., at $r=a$ one can observe a kink of surfaces, which means a discontinuity in concentration gradients on $r=a$. It means that ${\left.\frac{\partial c_1}{\partial r}\right|}_{r=a-0}\ne {\left.\frac{\partial c_2}{\partial r}\right|}_{r=a+0}$. This phenomenon is characteristic for macroscopically inhomogeneous media, where on the boundary of two media gradients of many characteristics break. One can notice the progress of the surfaces in the directions x and r with increasing the time t. The same phenomena occur with the pressure distribution (Figure 3).

Using pressure distribution (Figure 3) according to (1), (2) filtration velocity fields are determined. On Figure 4 we present filtration velocity contour lines. Filtration velocity is expressed through $v_1=\sqrt{v_{1r}^2+v_{1x}^2}$, $v_2=\sqrt{v_{2r}^2+v_{2x}^2}$, where $v_1$ and $v_2$ are filtration velocity modules in ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$, respectively.

One can observe that filtration velocities in ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$ are considerably differ. On the common boundary of ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$ we have a drastically change of the filtration velocity due to different values of permeabilities in ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$.

Figure 5 and Figure 6 show the surfaces of $c_1$ and $c_2$ for values ${\gamma }_1^{*}$ and ${\gamma }_2^{*}$, less than 2. From the presented material, one can notice an increase in diffusion effects in the direction x in both zones, when ${\gamma }_1^{*}$ takes values less than 2. When ${\gamma }_2^{*}$ decreases from 2, “fast diffusion” is observed mainly in the zone ${\mathsf{\Omega }}_2$. In addition, surface $c_1$ in ${\mathsf{\Omega }}_1$ do not undergo much change. Thus, the results obtained show that the transport characteristics of the solute mainly depend on the distribution of concentration and other hydrodynamic parameters of the zone ${\mathsf{\Omega }}_1$. The transport characteristics in ${\mathsf{\Omega }}_2$ are determined depending on the characteristics of ${\mathsf{\Omega }}_1$.

Numerical experiments show that the decreasing of ${\gamma }_1$ and ${\gamma }_2$ from 1 leads to the “fast diffusion” in the radial direction similar to that as in x direction. The smaller values of ${\gamma }_2$ than 1 that ${\gamma }_1=1$ leads to the “fast diffusion” almost only in ${\mathsf{\Omega }}_2$, whearase in decreasing of ${\gamma }_1$ from 1 we can observe “fast diffusion” in both zones ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$. In other words, solute transport characteristics in the micropore are influenced by the corresponding characteristics of the macropore.

5. Conclusions

In this paper, an anomalous solute transport problem in a two-zone coaxial cylindrical porous medium is considered. The medium is treated as a non-homogenous porous medium of fractal structure consists of macropore (inner cylinder) and micropore (external cylinder). To solve the problem a numerical algorithm is developed on the basis of finite difference schemes. Solute transport is modelled by partial differential equations of fractional order in both zones. To approximate the fractional derivatives Caputo’s definition is used. On the common boundary of the zones solute transfer can occur. On the basis of computer experiments it is shown that the decrease in the order of the derivative in the diffusion term of the transport equation in the macropore leads to “fast diffusion” in both zones, while a decrease in the order of the derivative in the diffusion term of the transport equation in the micropore leads to “fast diffusion”, mainly in the micropore. It was found that the characteristics of the solute transport in both zones mainly depend on the concentration distribution and other hydrodynamic parameters in the macropore. A decrease in the order of the derivative in the diffusion term of the transport equation in the macropore from 2 leads to an increase in the relative solute mass transfer through the common boundary of the zones. Received results show us that solute transport model in two-zone media with different flow and transport characteristics describes the process much better than conventional mobile-immobile solute transport models. As a consequence, one can hope that the presented model can more adequately describe solute transport processes in two-zone media. Because, as mentioned above, two-zone mobile-immobile model with the Fickian pore-diffusion solute transfer schematization can be unable to describe experimental data.