Comparison of two methods of numerical tracking of the soil contamination dynamics during a leak from a pipeline

The situation of leakage of a polluting liquid from a longitudinal crack of the pipeline lying on the ground surface is considered. The two-dimensional nonstationary mathematical model is based on the mass balance equation in terms of pressure, which is satisfied in a domain with an unknown moving boundary. This area corresponds to the area of contaminated zone. A function characterizing the region of action of the equation is introduced, which makes it possible to obtain the formulation of the problem in a fixed domain. Two types of finite-difference approximation of the problem statement are proposed. They differ by approximation of the convective term. Counter-current approximation and approximation along characteristics are used. The results of computational experiments, which are in favor of using the method of characteristics, are presented. The methods application is illustrated by an example of spread of oil pollution.


Introduction
Man-caused pollution of soil and water resources, including groundwater, is constantly in the spotlight all over the world. For instance, various aspects of occurrence and propagation of pollution of soil and water resources in various regions of the world (e.g., Europe, China, India, Brazil, Africa) are considered in the works [1][2][3][4][5][6][7]. Methods of modeling and numerical simulation of the contaminant transport are studied, for example, in [8][9][10]. Special attention is paid to modeling of the consequences of catastrophic failures of oil pipelines [3].
This work is devoted to modeling the spread of pollution into the ground from its surface caused by an oil pipeline crack.
The problem under consideration is characterized by the fact that the equation is written in a domain which boundary is mobile and must also be found simultaneously with the solution of the equation inside the domain with this boundary. Such tasks can not be solved by any known software package. The construction of methods for solving free-boundary problems and their investigation is constantly in development [11][12][13][14][15][16][17][18][19][20][21][22][23].

The problem statement
The situation is considered when a longitudinal crack is formed in the pipe lying on the ground surface and leakage of the polluting liquid occurs through it. As a result, a stain of contamination is formed in the soil (figure 1). The segment [ ] corresponds to the crack in the pipe through which liquid flows. It is assumed that the soil is saturated inside the stain. Outside the contaminated zone pressure and saturation are equal to zero. The boundary of the contaminated zone ( ) changes its form with time.
Since the contamination zone is symmetrical, one-half of this zone is considered and referred to as .
The mathematical model has the form: where is the porosity of the soil, is the fluid pressure, is the filtration velocity vector according to the Darcy's law [8]: Here, and are hydraulic conductivity for and directions. To describe the compaction of soil with depth, the model dependence ( ) is introduced in the direction of the axis, and ( ) is in the direction. Constants , , mean viscosity, fluid density and acceleration of gravity respectively.
In order to formulate the problem by dimensionless quantities, new variables ̃ ̅ , ̅ are the scale of the problem for , , , , respectively.
Neglecting the tilde sign above dimensionless variables, a dimensionless formulation of the problem is obtained: where , ̅ ̅ ( ̅ ) , ̅ ̅ ( ̅ ) and is the saturation of the soil.

Material and methods
То formulate the problem in a fixed domain, the is enclosed into the rectangle [ ] [ ], the size of which is sufficient to ensure that the contamination never reaches the boundaries of the rectangle (figure 2). The function ( ) is introduced for the fixed domain . It is the Heaviside function [24]. Thus, the statement of the problem in a fixed domain is obtained: Further, the soil is assumed to be homogeneous, i.e. ( ) ( ) and . An implicit finite-difference approximation is constructed for the differential problem in a fixed domain. The solution is found by relaxation method [25].

Results and discussion
A set of computational experiments is carried out. The pressure decreases with depth for each fixed and decreases for each in the direction from the symmetry axis to the contaminated zone boundary.
Comparison of the computational results for the counter-current approximation and approximation along the characteristics is made.
For the calculations, the results of which are presented below as typical, the following values of the parameters are taken: the hydraulic conductivity , the fluid pressure in the pipe , , , . It should be noted that the grid steps in space for all presented results are . At the same time, the observed pointthe maximum depth of contaminant penetration at each time point moves during one time step to a distance of approximately (greater than ) for the time step , and a distance less than for . In the latter case, the boundary turns out to be smoother ( figure 3). At such an average flow velocity, the contamination will reach the groundwater table in ̅ ̅ ̅ . It should be noted that the average filtration rate differs from the actual fluid velocity. The time step ̅ . The estimation time in 20 steps will be 1.12 days = 27 hours. The model pressure in the pipe can be estimated by selecting the pressure measurement scale as ̅ ̅ (technical atmosphere). For the computational experiment the pressure in the pipe is set 10 times less than ̅ , which can correspond to an oil storage. In the real pipeline, the liquid moves on a pressure of several atmospheres. Then the calculation area in the program must be set to a much larger size. It can be seen on the figure 3 that during the time in 20 time steps (27 hours) the oil reaches a depth of (line (c)), which is (about ). Thus, the choice of approximation along the characteristics gives a better result than the countercurrent approximation. In this case, there is no scattering of the values of the function , which is responsible for determining the boundary of the contaminated zone, if the grid steps in space and time are chosen so that the characteristics pass through grid nodes.