MODELING OF THE FILTRATION PROCESSES IN A RECTANGULAR AREA SOILS USING THE DARCY

Modelling of filtration processes using a semi-empirical system of Darcy equations is characterized by the following [1]: – determining the type of boundary conditions which make it possible to consider water permeability of soil; – developing the algorithm and implementing a numerical method for solving the problem on filtration in soils; – designing a user-friendly software module interface for carrying out mathematical calculations based on the developed algorithm of implementation of the chosen mathematical model. Mathematical model of the filtration processes takes into consideration a special feature of outflow – permeability of the medium, rate of filtration, fluid pressure, peculiarities of configuration of the examined areas. The proposed simulation of filtration processes is a relevant task as it enables obtaining a detailed quantitative description of the processes, as well as exploring the patterns of fluid propagation in a porous medium. Results of modelling are important for forecasting the level of floodwater and soil water because they describe features of the filtration flow (the flow of liquid; dynamic viscosity and density of liquid).


Introduction
Modelling of filtration processes using a semi-empirical system of Darcy equations is characterized by the following [1]: -determining the type of boundary conditions which make it possible to consider water permeability of soil; -developing the algorithm and implementing a numerical method for solving the problem on filtration in soils; -designing a user-friendly software module interface for carrying out mathematical calculations based on the developed algorithm of implementation of the chosen mathematical model.
Mathematical model of the filtration processes takes into consideration a special feature of outflow -permeability of the medium, rate of filtration, fluid pressure, peculiarities of configuration of the examined areas. The proposed simulation of filtration processes is a relevant task as it enables obtaining a detailed quantitative description of the processes, as well as exploring the patterns of fluid propagation in a porous medium. Results of modelling are important for forecasting the level of floodwater and soil water because they describe features of the filtration flow (the flow of liquid; dynamic viscosity and density of liquid).

Literature review and problem statement
The tasks of modelling filtration processes were investigated by many authors, the basic patterns of groundwater flow were established and described [2], which is expressed by the parabolic Boussinesq equations: where h=h(x, y, t) is the function that assigns free surface of the groundwater; k=(S 0 ρg)/m is the coefficient of permeability, ρ is water density, g is the free fall acceleration, m is the volumetric mass, S 0 is the coefficient that is determined  (1) is the fact that equations (1) is nonlinear, it lacks an analytical solution. Determining boundary and initial conditions for the area under investigation is a complex task since the actual areas are large in geometrical dimensions and they have a significantly complicated geometric configuration [3].
Paper [4] addressed the issues of filtration flows and the influence of intensity and orientation of the filtration processes on the ecological-reclamation condition. Filtration processes that accompany a cyclic operation of the oil well whose downhole pressure is not constant but changes over time periodically, are described in [5]. The shortcomings of these technologies are the narrow scope of research because they focus on examining the filtration processes only for rice irrigation systems and in the oil industry.
There is a known algorithm for solving the problem on finding the unknown flow dividing lines and the points of flow "suspension" [6], which is employed to automatically build dynamic grids and to calculate filtration rate. However, the specified algorithms do not take into consideration water-physical properties of soils.
Article [7] describes basic processes and phases of the boundary stressed state of soils without taking into consideration the peculiarities of the filtration flow.
The problem on consolidation of the soil array with regard to mechanical suffusion at the border of contact between two soil layers with different fractions is examined in [8]. This problem does not take into consideration the compressibility of fluid and the porous skeleton of soil.
There is a known mathematical model of the two-dimensional problem on filtration consolidation of soil taking into consideration the influence of multicomponent salt solutions and non-isothermal conditions [9] that does not determine configuration of the free surface of groundwater at the initial point of time.
In [10], authors addressed the issue of soil preservation, improvement of their regime and properties, increase in fertility, with regard to only one water-physical property (water permeability).
Given the above, the following tasks remain unresolved: -not taking into consideration water-physical properties of soils; -configuration of the free surface of groundwater at the initial point of time; -narrow scope of research into filtration processes. There are also objective difficulties associated with the development of adequate mathematical models for filtration processes that would make it possible to estimate the level of floodwater and groundwater, to describe the features of the filtration flow as there is a lack of certainty in terms of boundary and initial conditions for the area of research into filtration processes.

The aim and objectives of the study
The aim of present work is to develop a mathematical model and a corresponding numerical method for modelling the process of fluid filtration flow in the region with a discharge of fluid through the surface. This will make it possible to obtain new results in solving the tasks on the evaluation of environmental risks in the case of emergencies at the sites of oil and gas industry. Such a technique could also be applied when predicting the development of flood phenomena in mountain rivers.
To achieve the set aim, the following tasks must be solved: -to improve the models of filtration flow based on the Darcy and Forchheimer laws taking into consideration the specificity of area configurations; -to develop and explore convergence of the numerical method for solving the system of equations with partial derivatives that describes a filtration flow in the examined rectangular plot of soil using the Darcy equations; -to conduct testing and simulation calculations employing a software module.

Mathematical fundamentals of modelling filtration processes in soils. Numerical method for solving the Dirichlet problem
One of the most effective methods for predicting the level of flood waters is the method of mathematical modelling of processes that occur when a significant amount of liquid penetrates the soil at large volume of precipitation. In this case, the most appropriate model is a model of filtration processes that employs the Darcy equations, which form a system of differential equations with partial derivatives: where u, v are the components of rate of the fluid filtered through the medium, ρ is the fluid density, p is pressure, µ is the coefficient of dynamic viscosity of the fluid, g is the free fall acceleration, x, y are coordinates, K is the coefficient of filtration (permeability). It is important when solving the task of modeling filtration process to assign appropriate boundary conditions. Boundary conditions are assigned under assumption about the fact that the motion of fluid is predetermined by the difference in pressure between the area and the border.
When constructing a mathematical model, we accept the following assumptions: -the flow of liquid in a cylindrical form with a bottom radius R and height H is symmetric relative to the vertical axis of the area, that is, we actually consider a flow in a rectangular region; -dynamic viscosity µ and fluid density ρ fluid are considered constant; -temperature of the medium and its permeability are considered constant; -pressure at points is determined only by the atmospheric pressure and the pressure of the column of fluid.
A problem on the estimation of filtration rate in the medium with resistance in a rectangular region: under condition V>>V 2 , where V is the module of velocity vector V(u; v), is reduced to the problem of solving a system of the Darcy equations: where u, v are the components of the filtration velocity vector; p is the pressure of the liquid. In the case when condition V>>V 2 is not satisfied, the last two equations of system (4) are recorded using the Forchheimer filtration law: where ß is the dimensionless coefficient, which depends on the structure of the porous medium, ß≈1; k is the coefficient, which coincides with the coefficient of permeability by dimensionality; V is the characteristic filtration rate. Systems (4) and (5) are supplemented by boundary conditions of the form: where g(x, y) is the piecewise-continuous function. Conditions (6) describe the case when the walls or borders of region G are freely permeable to liquid. In the case of applying conditions (7), we solve the problems on exploring filtration flow under condition that part of region G has borders that are impenetrable for liquid. In addition, solution to problems (4) and (5) under conditions (6) allows us to examine the size of the zone of influence of the presence of outflows from the area on the flow configuration as a whole.
In the process of modelling, we made the following assumptions: -modelling problem is solved for the area that has the shape of a rectangle with a given height H and width L ( Fig. 1) whose bottom part has openings, through which the fluid is filtered; -fluid is characterized by constant viscosity µ (Pa•s) and density ρ; -temperature of the medium is constant; -filtration rate is such that the magnitudes of order V 2 can be neglected; -pressure at all points of the area is determined by the atmospheric pressure and the pressure of the column of liquid, which is due to the smallness of magnitude V 2 ; -resistance of the walls of the area is considered to be zero.
In this case, boundary conditions for pressure are recorded in the form: .  (7) take into consideration that the liquid flows through the area at section [y i , y i+1 ]. The introduction of indexes means that there can be several such zones of outflow. By assigning boundary conditions in the form (8), we model a configuration of soil classes at which the flow of fluid occurs only at certain depths: A is the zone of impenetrable rock, B is the zone of penetrable rock. Under condition µ=const, K=const, system (2) is substantially simplified and is reduced to a Dirichlet problem on the pressure of fluid, recording (2) in the form: hence, we obtain with regard to the first equation of the system: which is solved with conditions (7). A Dirichlet problem (or the first boundary value problem) is the task on finding a field-regular solution to the second-order elliptic equation, which accepts the preset values at the area's borders.
Because the boundary conditions are discontinuous, in order to find a solution to problem (9) with conditions (8), the numerical methods are applied [8]. To solve the Dirichlet problem with appropriate boundary conditions numerically, the upper relaxation method by rows is used [8]. The mathematical features of solving a problem employing the upper relaxation method are considered in paper [9] where the relationship between the values of relaxation parameters and the convergence of an iterative process was established.
We select for equation (9) a computational grid of the area (Fig. 2) with a similar step for the x and y coordinates, which is equal to h. I is the number of partition points along x, J -the number of partition points along y.

. Sequential upper relaxation in rows
Using the method of finite differences, we shall obtain the following estimation formula for implementing an iterative process: In relation (11), ß=∆x/∆y=1-y in the case of the same step of the grid, ω≤ß 2 +1 is the parameter of upper relaxation, k is the number of an iteration, + is the pressure at partition point х і and у і at the step of iterative process number k+1. Equation (11)  is taken from the preceding step of the iterative process, magnitude + -1 , 1 k i j p is taken from the preceding layer at the step of iterative process number k+1. To implement the specified method, it is required to assign the initial approximation of pressure distributionfor example, at the first step of the iterative process pressure distribution is assigned by dependence (7). System (9) is solved at each step by the Thomas algorithm.
The Thomas algorithm is applied to matrices, in which most of the elements are equal to zero. In this case, the structure of the matrix is not chaotic, but fully determined, specifically, the matrix is tridiagonal, nonzero elements are located on the main diagonal and the two adjacent to it. The Thomas algorithm is a particular case of the Gaussian method and it also consists of direct and reverse sweep. To solve the system, the matrix must first be reduced to the bidiagonal form: Dividing the first row of the matrix above by -b 1 , it is obvious that P 1 =c 1 /b 1 ; Q 1 =-d 1 /b 1 and we can derive a formula for the direct sweep: Next, it is necessary to perform the reverse sweep -to find vector X -it follows from the last row of the transformed matrix that x n =Q n . In this case, other elements of the vector are calculated from formula: The Thomas algorithm is stable if │P i │<1, i=1,…, N (which follows from the diagonal dominance in the matrix) and correct if (b i -a i )P i-1 )≠1 (otherwise the formulae of the direct sweep do not make any sense).
By assigning the initial approximation of function p 0 =(x, y) taking into consideration the boundary conditions, we obtain a system of linear algebraic equations with a tridiagonal matrix: where A=ω/(2(1+ß 2 )); C=1-ω; B=Aß 2 ; ß=∆x/∆y; N is the number of partition points along x; s is the number of layer along y.
With this approach, sequential upper relaxation is included in the algorithm for solution in every row rather than as a separate step of the procedure of solution. Since, when using the program, it is desirable to ensure the diagonal dominance, then condition ω≤1+ß 2 must be maintained in the algorithm. Iterative process is interrupted when the difference between pressure values at the adjacent steps of the iterative process does not exceed a given level of accuracy ε. Upon finding pressure distribution, components of the velocity vector are calculated from the second and third equation of system (2).

Results of studying the algorithm and the implementation of a numerical method for solving the problems on filtration in soils
Studying the algorithm and the implementation of a numerical method for solving the problems on filtration in soils were carried out based on the algorithms for solving hydro-and gas-dynamic tasks of hydrodynamics and heat exchange [9].
The algorithm of solving the problem implies the following: -initial approximation of solution p 0 =(x, y) is assigned; -at each layer, along the y coordinate, system (13) is solved, as a consequence of which, using known p K =(x, y), approximation p K+1 =(x, y) is found; -if the following condition is satisfied: where ε is the assigned level of accuracy of the solution to the problem, then the iterative process is terminated. The important point is that systems (3) and (4) can be reduced to the form: and the form of constants C 1 and C 2 depends only on the type of model ((3) or (4)). Thus, it will suffice to determine the distribution of pressure at the examined plot and its partial derivatives along x and y. Absolute values of components of the velocity vector can be determined as a linear combination of the found derivatives with appropriate coefficients. The coefficients characterize the environment of filtration (coefficient of permeability) and the substance that is filtered (dynamic viscosity and density).
Employing the specified algorithms that implement a numerical procedure of the simulation of filtration process, we developed a software module in the Delphi programming environment (Fig. 3).
The main software interface is shown in Fig. 3. The developed software module enables carrying out test calculations of pressure distribution over the examined area with assigning various initial and boundary conditions.
The results of modelling the location of penetrable zone are shown in Fig. 4: Fig. 4, a -on the left border of width 10 points, and on the right border of width 3 points (the total number of points of the computational grid is 50×50), Fig. 4, b -in the absence of outflow on the left border.
The source data for the results shown in Fig. 4, b are given in Fig. 3, the field "Borders of the outflow".
The results shown in Fig. 4 are obtained by implementing the above-described computational algorithm: initial conditions were assigned based on dependence (8), an iterative procedure to calculate the pressure field is performed by formula (11).
The result of the simulation of the horizontal velocity component for the same cases of the presence of outflow is shown in Fig. 5.
Exploring a horizontal component of velocity for the height of the region revealed that the impact of the presence of these zones is felt only in the small vicinity of these zones, that is, the presence of outflows for region's height almost does not affect the parameters of flow at the bottom of this region. a b Fig. 4. Pressure distribution: a -outflow on the right and the left border; b -outflow on the right border A criterion for terminating an iterative process is the fulfillment of condition (12). The estimate of the magnitude of error is given graphically using the dependence shown in Fig. 6.
In addition, the software module makes it possible to study stability of the iterative procedure depending on the parameter of relaxation. Results of calculations are given in Table 1. The data given in Table 1 were obtained while applying a computational grid 50×50 and at the magnitude of calculation error of 0.0001. An analysis of the obtained data allows us to choose a relaxation factor in order to optimize conducting the calculations. Thus, when implementing a numerical method, we obtained the following results: -we created a numerical algorithm and a program for the implementation of the algorithm, which makes it possible to calculate the pressure field in the examined region. We conducted test calculations and established that the convergence of the iterative process is achieved in 4-5 steps at optimal selection of the relaxation parameter; -we obtained a distribution of velocity fields for the model region taking into consideration the dependence of medium on the y coordinate; -it was established that under boundary conditions, recorded in form (7), the filtration process is carried out in such a way that one does not observe any accumulation of fluid on the height of the region.
We performed calculations for the model regions which are penetrated by water from outside and where the phenomenon of filtration is observed.
The developed mathematical procedure and software module could be a component of the system for simulating a filtration process with a known value of the dynamic viscosity of the filtered liquids and the coefficient of permeability [11,12].
As a result of modelling the processes of filtration in soils, the following was established and implemented: -we improved a numerical algorithm for the calculation of pressure field in the examined region by assigning the kind of boundary conditions and building a converging iterative process at optimal selection of the relaxation parameter; -we conducted a qualitative analysis of the size of the fluid propagation zone when implementing technological processes of shale gas extraction.
The advantages of the present research over analogues include ensuring the stability of solving a numerical scheme through selection of the optimal relaxation parameter and taking into consideration properties of the medium that borders the region of simulation. An alternative mathematical model is the model of filtration consolidation of soils taking into consideration multifractional suffusion [9] that does not establish configuration of the free surface of groundwater at the initial point of time.
It should be noted that a shortcoming of the study is the underestimation of morphological characteristics of the river and the river system. That is why continuation of the present research implies engagement of technical means in order to obtain operational data on a change in the level of groundwater and flood waters in real time to ensure ecological safety of the environment.

Conclusions
1. We performed modelling of filtration flow based on a system of the Darcy equations by taking into consideration a wide class of boundary conditions (a velocity vector module) when simulating the flow with an outflow of fluid through the surface of the examined region, which would make it possible to prevent the rising of water levels in rivers during floods.
2. We implemented a numerical method to solve the Dirichlet problem using a difference method of sequential upper relaxation and taking into consideration different values of the relaxation parameters and properties of the examined medium (temperature, permeability) in order to study behavior of ground water in terms of ecological (floods, territory flooding) and technological (shale gas extraction) aspects.
3. An analysis of the obtained numerical results was performed with regard to the properties of soils in the flooded areas and zones with possible emergency outflows, the volume of fluid that penetrates the area, as well as hydrological, climatic and natural conditions of the regions for which using the proposed software module is recommended. When implementing the numerical method, we obtained the following results: distribution of velocity fields for the model region taking into consideration the dependence of environment on coordinate; it was also found that the process of filtration occurs when fluid is not accumulating on the height of the area.