Simulation of the process of plane-radial filtration of fluid in a porous medium considering the migration

. Modeling plane-radial fluid filtration in a porous medium is considered in the article, taking into account the condition of fluid flow to a perfect well with a constant flow rate drilled in the center of the formation; the problem of identifying the parameters, entering the mathematical formulation to the problem, is solved. Solving the problem, the values of the following parameters were determined: T - the coefficient of water conductivity, which ranges from 15% to 27%, and 𝜇 ∗ − the coefficient of water loss (water filtration), which ranges from 22% to 31% relative to the exact values used in solving the forward problem. In the course of the study, it was revealed that one of the main parameters in the process of ground and underground water filtering in porous media is the identification of evaporation parameters. To identify their values, the following methods were applied in the study: the method of barycentric coordinates; the least squares method; and the gradient method. Numerical experiments were conducted at a time when the critical depth of groundwater was a known value, while 𝑞 0 , 𝑛 were the unknown parameters. The analysis of the studies showed that the "discrepancy" takes the minimum value when the found parameters (taken as zero approximations) are used in the least squares method.


Introduction
An analysis of long-term statistical data showed that the deterioration of the hydrodynamic regime of groundwater and the process of their mineralization led to a decrease in the fertility of coastal pastures and lands of the Central Asian region.In particular, the decrease in the level of the Aral Sea has significantly changed the aquifers, the exploitation of groundwater in the coastal zone, and had a negative impact on the environment.
The total natural groundwater reserves in Uzbekistan are estimated at 24.35 km 3 .Of this amount, 20.79 km 3 is in the Quaternary layer, 2.92 km 3 is in the Upper Pliocene-Quaternary layer, and 0.46 km 3 is in the Upper Cretaceous layer.Fresh groundwater is concentrated mainly in the Fergana Valley (34.5%) and Tashkent (25.7%),Samarkand (18%), Surkhandarya (9%), and Kashkadarya (5.5%) regions, and the rest of groundwater reserves are saline or brackish, their use is low.Fresh water in other areas amounts only to 7%.
In the western part of the country (the lower reaches of the Zarafshan and the basins of the Kashkadarya, Syrdarya, Amudarya Rivers and Central Kyzylkum), groundwater is highly mineralized and hard.Subsurface lens waters used to supply drinking water to the Khorezm region and the Republic of Karakalpakstan, formed by the flow of large rivers (Amu Darya and irrigation canals), over the past 10-15 years do not meet the requirements of national standards due to increased mineralization and hardness (due to irrigation effect).
From this point of view, it is possible to predict changes in the level and salinity of groundwater, as well as to make decisions and prepare recommendations by developing mathematical models that describe the process under study.
The studies on the state of groundwater have shown that in most cases the main aquifer from which pumping is done is blocked from above by a low-permeable cover layer, and from below, it is underlain by a low-permeable layer through which the seeping into the underlying aquifer occurs.Based on the above, the low-permeability interlayers containing groundwater were considered non-compressible, or at best, they were considered to contain elastic reserves, though insignificant ones.Consequently, when studying the process of groundwater filtration, Moore, Shilsueiz, Hurst, Masket, and others took into account only the elasticity of the fluid, but did not consider the elasticity of layers, and the discrepancy between theoretical calculations and actual data led to the need to create a more accurate notion of the elastic regime of filtration in a porous medium.As practice has shown, under natural conditions, the layers are in a compressed state, and with the start of pumping groundwater by a well, a decrease in pressure near the wells occurs, and in the process of pumping out, pressure increases in the zone of drawdown.The released water due to elastic reserves reduces the efficiency of the well, as shown by the study; accounting for released water makes it possible to more accurately reflect the physical pattern of the groundwater flow in layered soils.
Research conducted by V.V. Shchelkachev, V.M. Shestakov, N.N.Verigin, I.A. Charny et al. was aimed at solving this problem; they obtained significant results of a fundamental and applied nature, which took into account the regime of groundwater filtration in multilayer media.The studies of Hantush and his students were devoted to the process of groundwater filtration, considering the elastic regime in low-permeabe layers.
The main empirical law of groundwater filtration is Darcy's law, the experimental method obtained in 1856 by the French engineer Henri Darcy; this law, later called Darcy's law, connects the flow rate of the underground flow with energy losses during its flow.
F. Forchheimer, J.Dupuy.N.E.Zhurkovsky, F.B. Abutaliev, P.Ya.Polubarinova-Kochina, V.I.Aravin, S.N.Numerov, N.N.Verigin, G.N. Kamensky and others played a significant role in the development of mathematical models, effective numerical methods and algorithms for a comprehensive study of the process of groundwater filtration.
It was suggested in [1,2], that the interaction between surface and ground waters is important in land reclamation, engineering hydrology and hydrogeology; mathematical modeling of the process of non-stationary fluid filtration in a non-homogeneous porous medium was considered.The mathematical model was developed on the basis of partial differential equations of parabolic type with shifted initial and boundary conditions, and an analytical solution to the problem was found using the Laplace transform.Computational experiments were performed to determine the change in pressure along the length of the filtration layers without considering the elastic regime.Numerical calculations have shown that the pressure in both layers grows exponentially, and the fluid flow through the interface between the filter layers depends significantly on the piezoconductivity coefficient of the well-permeable layer.The authors of the article obtained a new generalized formula for controlling well adits and proposed a mathematical apparatus that allows drawing up layouts and capacity of vertical drainage wells to protect irrigated and non-irrigated areas from flooding, to protect groundwater from pollution sources, and isolate already polluted areas.
In [3], the authors stated that arid and semi-arid regions face serious challenges in managing scarce freshwater resources due to economic development and climate change.As noted by the authors, accurate groundwater level prediction is an essential component of efficient water management, and physically based models are often used to model and predict groundwater flows.The focus of this study is on the application and comparison of three data-driven models to predict short-term groundwater levels.The paper evaluates and compares a set of popular data-driven models, including artificial neural networks (ANNs), support vector machines (SVMs), and the M5 model tree.The feasibility and capabilities of these models are demonstrated by predicting groundwater levels five days in advance in arid and semi-arid basins located in northwest China.Encouraging modeling results show that the techniques can simplify and improve the procedure for predicting groundwater levels.
In the dissertation work [4], fluid-dynamic processes in the underground hydrosphere were studied by the methods of computational experiments; according to field studies, models were developed for calculating physical quantities characterizing these processes, namely, fields of pressure, temperature, and deformations in porous and fractured-porous layers of the upper lithosphere.The results obtained link together thermal, filtration and deformation processes; this is an important step toward creating complex dynamic models of the evolution of the upper lithosphere under the influence of a combination of natural and anthropogenic factors of varying intensity and duration.In practical terms, all this is the basis for building permanent models for the development of hydrosphere resources, as well as for creating systems for environmental monitoring of near-surface waters.
An exact two-parameter solution to the Cauchy problem on the flow of two-layer shallow water under a solid cover of two infinite contacting liquid layers with a small difference in density moving at different velocities in a horizontal channel with solid walls, was constructed in [5], where the distortion of the interface between layers occurs due to the Kelvin-Helmholtz instability.The process is described by a system of two quasi-linear partial differential equations of the hyperbolic type of the first order, and to construct the solution, an option of the hodograph method based on the conservation law was used; it makes it possible to transform the system of partial quasi-linear equations of the first order to a partial linear equation with variable coefficients of the second order.
Modeling the process of groundwater filtration in multilayer porous media is reflected in [6][7][8], in these publications, the level of groundwater in multilayer reservoirs, as well as salt concentrations in them are predicted.The authors of these publications have developed a mathematical tool for designing and developing (in hydro-technical structures) regulated groundwater flows to prevent flooding, salinization and waterlogging of lands, which cause enormous damage to the national economy.The research works related to this problem were analyzed in detail in these articles.To derive a mathematical model of the process, Darcy's law was used considering infiltration (rain, watering), and evaporation.To solve problems that describe the process of nonlinear partial differential equations, a numerical algorithm was developed using finite-difference schemes, besides, an iterative scheme was used for nonlinear equations in order to check their convergence.
To obtain a one-dimensional analytical solution for the conservative transport of dissolved substances along an unsteady groundwater flow in a semi-infinite aquifer, a comparative study of the Laplace transform and Fourier transform methods was conducted [9].The time-dependent source of pollutant concentration was calculated at the start; at the other end of the aquifer, it was assumed zero.Initially, the aquifer is not free of solutes, which means that the concentration of solutes enters the groundwater system, and is uniform.The aquifer is assumed homogeneous and semi-infinite, and time-dependent velocity expressions are considered.
The study in [10,11] solved the problem related to the process of changing the level of groundwater and the transfer of mineral salts in soils; the problem is described by a system of partial differential equations and the corresponding initial, internal, and boundary conditions.To conduct a comprehensive study of the filtration process and changes in the salt regime of groundwater, mathematical models and an efficient numerical algorithm were proposed, taking into account external sources and evaporation.
In [12], the process of predicting changes in the level of groundwater and pressure water was considered.For a comprehensive study of the problem under consideration, the authors developed a mathematical model that took into account the external source, evaporation, filtration coefficients, active porosity, filtration rate, and bilateral boundary conditions.An efficient numerical algorithm was developed for predicting changes in the groundwater level using a combination of finite difference schemes and run methods.Changes in the level of ground and pressure waters and filtration permeability were studied.
Reference [13,14] proposes a mathematical model for conducting a comprehensive study of groundwater filtration, taking into account the clogging of soil pores by fine particles over time; changes in the coefficient of water permeability of soil, water loss and filtration coefficient; changes in the initial porosity and porosity of the settled mass and an efficient numerical algorithm based on the Samarsky-Fryazinov vector scheme.In that article, to construct a mathematical model of salt transfer, the authors assumed that the pressure gradient in the channel is constant and equal to atmospheric pressure.It was established that with poor irrigation, the maximum absorption of water and the transfer of accumulated salts occur in the upper layers of soil.Numerical calculations showed that the change in the rate of water transfer in soil depends on porosity, water permeability of soil, filtration coefficient, composition and structure of soil, and porosity of the settled mass.
In [15], the concentration profile is studied for a solution that saturates poorly permeable soil.The simulation results showed the presence of three flow regimes.The accumulation of salt near the phase transition boundary increases the density of the solution and can lead to the development of natural concentration convection interacting with the upward flow (forced convection).The stability threshold of the forced flow and the influence of the emerging natural convection on it are determined.It is shown that with an intense flow to the surface of evaporation, the impurity concentration at this boundary rapidly increases and reaches the saturation concentration.In this case, the impurity is precipitated out.In the mode of slow evaporation, the impurity diffuses from the region of high concentration, which prevents the development of the convective flow.
The study in [16] presents a mathematical model for calculating evaporation from bare saline soils.The equations of water and heat transfer in a zone of partially saturated soil surface are related to the model of the boundary layer of the lower atmosphere.The transport model considers the multicomponent migration of salts.The flow equation takes into account the osmotic flow due to the presence of ions in the soil solution.Nonlinear one-dimensional equations were solved using fully implicit difference schemes.The soil of the model was tested using experimental data on the redistribution of water and dissolved substances under temperature gradients in sealed soil columns.To demonstrate the effect of soil salinity on evaporation, the problem of water infiltration into soil with initially low moisture content and its subsequent evaporation was considered.
In [17], an alternative numerical upward weighting method was proposed for modeling the transport process in groundwater under various conditions.The discretized system of solute transport equations was derived from the integral form expressing the mass balance in each element.The upward weight dependent on the local Peclet number was directly added to the basis functions.When convective transport dominates over dispersive transport, this method can efficiently eliminate solution oscillations with only a nominal increase in computational power compared to Galerkin's general linear triangular finite element method.The accuracy of the proposed numerical model was tested on two classical problems for which there were analytical solutions.

Formulation of the problem
Problem a).The problem of fluid inflow to a perfect well with a constant flow rate, drilled in the center of the formation, taking into account migration, is considered; it satisfies the following differential equation: (1) with initial and boundary conditions: Here:  is the radial distance from the well;  1 − the level;  2 -the pressure in the aquifer;  2 0 -the initial pressure; Q -the well flow rate;   -the well radius;  * -the coefficient of elastic water loss; T -the coefficient of water conductivity; t is time;  2 =     ;   -the interlayer thickness,   -the interlayer filtration coefficient.
Equation (1) also includes coefficient α, which is equal to one when migration is taken into account, otherwise, it equals zero.
From the statement of problem ( 1)-( 2), it follows that it is difficult to obtain an explicit analytical solution.For a numerical solution to problem (1) -(2), we introduce the following substitutions: Then we obtain: (3) For convenience, dashes are omitted and (3), (4) take the following form: To solve problem ( 5)-( 6), we introduce the change of variables  =   2 and obtain: Applying a finite-difference approximation to ( 7) and ( 8), the solution is sought in the following form: where Coefficients entering   ; and   are calculated as follows: Then, replacing the boundary conditions of the problem posed ( 7)-( 8) by finite-difference ones, we obtain the relation for determining coefficients  0 and  0 : Having determined the values of running coefficients   , and   , we calculate the values of heads  2, according to the recursive relation (9).
To check the calculation of the head values  2, , we integrate equation ( 5) over t and S and obtain: or replacing with a finite-difference approximation and integrating numerically, we obtain: where ∆ characterizes the accuracy of calculations.For the task of identifying the parameters included in equation (1) with condition (2), the solution to problem (1) with condition (2) is used as information for conducting numerical calculations on a computer.The analysis of computational experiments showed that the found values of parameter T -the water conductivity coefficient, range from 15% to 27%, and parameter  * -the water loss coefficient varies from 22% to 31% relative to the exact values determined in the solution to the forward problem.
Problem b).The conducted studies have shown that one of the main parameters in the process of filtration of ground and underground waters in porous media is the identification of evaporation parameters.To identify their values, we use the method of barycentric coordinates, the least squares method, and the gradient method.
For this purpose, we take the solution to the following problem as a mathematical tool: with initial and boundary conditions: Here:   is the free water loss or lack of saturation;   is the filtration coefficient of the cover stratum;  0 is the intensity of evaporation on the day surface;  1  is the critical depth of standing groundwater; n is the exponent.
As follows from the statement of the problem, the evaporation parameters are  0 , ,  1  , to determine their values, the following functional is minimized: To solve the above problem, a numerical algorithm and a software tool were developed to conduct numerical calculations using a computer system.The results of numerical calculations are given in Tables 1, 2.  The technique for calculating the parameters is performed in the following sequence: using the method of barycentric coordinates, calculating and obtaining the corresponding values of the parameters, the target discrepancies are checked and, depending on this, the calculation process either ends or transfers the values of the found parameters as zero approximations to the next algorithm -the least squares method or the gradient method.
Numerical experiments were performed when the critical depth of groundwater was a known value, and  0 , , were the unknown parameters.The results of numerical calculations are given in Table 2.In all numerical calculations performed using the method of barycentric coordinates, strongly perturbed zero approximations were taken.
The numerical calculations performed for various values of parameters  0 , n showed that the "discrepancy" took the minimum values when the found parameters as zero approximations were used in the least squares method (Table 2).

Table 1 .
The results of numerical calculations 1.

Table 2 .
The results of numerical calculations 2.