Simulation of Transfer and Spread of Pollutants into the Marine Environment Using a Probabilistic Model of Turbulence

The paper describes the problem and specifies the key sources of marine pollution. A mathematical mode of particle transfer and settling as well as the parameters of environment turbulence are proposed. A 3D diffusion-advective model discussed in the paper was developed on the basis of “wandering particles” method with the use of random numbers generator for nondeterministic processes simulation. But it has a number of specific features as compared with similar models in practical setting of boundary conditions at water-air boundaries and water-bottom boundaries, description of turbulence of different scales, in taking into account the effects of the initial immersion of a heavy jet, in the effect of flocculation as well as taking into account non-conservatism of the impurity. Different simulation results are presented.


Introduction
The methods for studying the problems of environmental safety are rapidly developing these days. The basis for identifying the environmental situation and developing measures to prevent and response pollution is comprehensive monitoring and detailed examination of polluted waters, unification of a number of parameters and methods that ensure their purification from pollution, including the construction of computer mathematical models to determine the degree of pollution.
The construction of computer mathematical models will allow predicting the formation of a pollution front and assessing the size of the polluted zone, which immeasurably expands the possibilities, both in its fundamental studies and in the field of its practical applications.
The engineering approach discussed below is used to calculate the distribution and settling of pollutants [1]. These days this problem is becoming topical for a promising area of the offshore oil and gas fields development. Later, mixing parameters will mean transfer characteristics (determined by the velocities of flows or particles) and design parameters of turbulence (vertical and horizontal turbulence caused by shear instability, surface impact, bottom friction, and other factors). In the first case, the transfer parameters determine the trajectory of the pollutant slick, which can be calculated with a sufficient degree of accuracy for a given hydrometeorological situation using the characteristics of wind, currents, and other factors. It is more difficult to determine the parameters of turbulence, that is, the characteristics of the evolution of a pollutant slick, for which there are no unified methods of description and much is determined by the approach proposed by Bogdanovsky A.A. and Kochergin I.E. [1].

Materials and methods
The classical problems of simulation of the suspended substances propagation in a turbulent flow are of interest from the two points of view. First of all, from a theoretical point of view, they are the simplest example of the phenomena associated with such an important and yet not fully understood phenomenon as turbulence. Secondly, from a practical point of view, the interest in such problems has significantly increased recently, in particular, because of the need to assess the influence of various anthropogenic impacts on aquatic biological resources.
The need for such assessments arises, for example, when planning the construction of drilling platforms in the ocean shelf, when laying subsea pipelines, when dredging the bottom during repair works to clean up intra-port waters and shipping channels of ports from sediments.
Theoretical studies of the process of turbulent dissipation of suspended substances are carried out using the apparatus of stochastic differential equations [2] and even the methods of quantum field theory [3]. However, when addressing practical problems in which it is necessary to simulate the process of propagation of various suspended substances in an aquatic environment, more traditional equations of mathematical physics are taken as the basis. When describing the distribution of suspended substances, two qualitatively different areas can be distinguished, namely the "near zone", the spatial scale of which correlates with the size of the object polluting the water area.
All the algorithms discussed in the paper are implemented in the form of computer programs in the Wolfram Mathematica programming languages.
The main impact of the considered sources of pollution of the aquatic environment is characterized by the formation of a plume of turbidity (suspensions) from the place of works in the direction of the currents as well as by zones of sediments falling on the seabed from settling particles of the solid phase of grounds in the area of the works.
Excavation of concrete from under the crib structure is planned to be carried out using a Primorets dredger. Excavation of concrete shall last for 2.92 hours. The volume of the concrete excavated will comprise 105 m 3 . The speed of concrete excavation will comprise 35.9 m 3 per hour.
When removing ground with a dredger, the following main types of operations lead to the formation of suspended substances in seawater: − loosening of the seabed with a bucket accompanied by bottom detachment of the excavated ground; − continuous lifting of the bucket with ground with leakage of water and ground particles from the bucket, the so-called ground spillage.
The calculation of the total amount of suspended substances entering the aquatic environment during ground excavation with a dredger is carried out as follows [5]: where -the total volume of excavated ground; -dredger spillage factor; , -the contents of fine fractions in the ground. The average rate of suspended matters entry into the water layer is calculated according to the formula: where Qdredging is dredger performance. According to the technical design, it is planned to remove 105 m 3 of concrete and free the seabed within the area of 555 m 2 .
According to studies carried out by DNIIMF, spillage for multi-bucket dredgers comprises about 10% for homogenous muddy bottom. According to [4], the spillage during dredging with multi-bucket dredgers can range from 0 to 8.5% and averages 0.8%. In this paper we adopt one of the most conservative variants of spillage as kspillage = 0.05. According to the analysis of file materials and similar projects, the factor kfine = 0.05.
Dredger performance in case of concrete excavation Qdredging = 35.9 m 3 /h. The particle-size composition of the grounds corresponds to the composition given in Table 1.
On the basis of the above formulas and factors, the technical parameters were calculated and a model of the source of impact on the marine environment during the excavation of concrete from under the crib structure was prepared ( The shape of the source is distributed in the overlying water in the form of a rectangular parallelepiped 2 × 2 × 6 m in size and occupies the initial position in the considered coordinate system.

Mathematical model
In the model description, the pollutants entering marine environment from the effluents are simulated in the form of a set of markers that retain the properties of mass, buoyancy, fractional distribution, etc. In this case, the markers interact with marine environment and the water-air and water-bottom boundaries. At the initial stage of mixing the markers can interact with each other. In some cases, their destruction can occur due to the non-conservatism of the impurity. The basic equations describing the trajectories of the markers in the water layer are represented as: where -the present coordinates of i-th marker; u, v, w -the components of advective transfer determined by the fluid velocity (tidal and nontidal currents); , , -turbulent pulsation velocity components, while for , in shelf waters the following spectrum ranges of turbulent pulsations are distinguished (depending on the area under study) -0.5 -15 minutes, 15 -60 minutes and 60 -360 minutes; -the velocity of vertical movement of markers under the action of gravitational force, describing substances with neutral, positive or negative buoyancy, at this, in relation to the discharges of sticking particles, the correction to the settling rate at the expense of the effect of "sticking" or "flocculation" is taken into account. So, it should be taken into account that even for the markers with neutral buoyancy the vertical movement does not equal to zero in some cases.
-additional velocity of vertical movement of markers at the initial stage of jet destruction, arising from the difference in the density of the jet and seawater. With regard to the unit discharge of heavy effluents, the jet descends more quickly in accordance with the "ball" effect The initial conditions for the system (Eq. 04) determine the spatial coordinates of markers outgoing from their source , , , and the corresponding turn-on time . The boundary conditions determine the effects that are possible at the boundaries of the media and at the boundary of the computational region and can be represented as a logical operator . When going beyond the boundary of the calculated horizontal region [0,X], [0,Y], the operator fulfills the condition of excluding markers from further calculation; Boundary conditions at the boundary with air Z≤0 for solid phase of total reflection. In case of non-vaporization of the marker, operator requires the fulfillment of the conditions for vertical velocities: (5) At solid boundaries z => h( , ) (where h ( , ) is the depth of the basin at point , ,), the conditions for the complete reflection of dissolved substances are set. In case of reflection of the marker, operator requires that the horizontal velocities correspond to those calculated for the bottom and that the following conditions are met:

Characteristics of transfer caused by particles settling
For the solid phase, the settling rate of the markers during simulation is determined by two mechanisms -natural sedimentation due to gravitational forces and resistance, and in the case of a certain type of grounds -by the addition due to the effects of jet lowering and flocculation.
The rate of natural settling for particles of the solid phase in (2.1) is found from the solution of the balance equation of three forces, for spherical particles of the reduced equivalent radius (diameter), which can be written in the form of the formula by Todes O.M., Rosenbaum R.B., taking into account the influence of the irregularity of the shape of the falling particle [6]: (7) where g -downward acceleration; -effective particle diameter; -particle specific gravity; ρ(z) -fluid density at the depth of the marker location; k -geometrical factor of the shape; μ -viscosity factor determined by the empirical relation (for fresh water): (8) where t -Celsius temperature.
Mechanism for determining in (2.4) the equivalent diameter of a ball-shaped particle for irregular-shaped particles with a certain degree of smoothness is described in [6] (Kurganov, Fedorov, 1973). For engineering calculations the shape irregularity factor ]. The effect of "stickiness" (flocculation) takes into account the adhesion of solid phase particles to each other and the formation of bonded combinations of particles. In foreign literature, to describe this phenomenon, the term flocculation is used (floc is translated as snowflake) [7]. Sticking is characteristic of clay fractions of ground or colloidal solutions, while suspension particles are enveloped in a suspension, which facilitates their bonding for some time after mixing with seawater [8]. Such coarse particle (floc) is settled quicker.
The effect of jet, or "ball", sinking that is reflected in equation (Eq. 04) with the parameter , manifests itself during intensive discharges and consists in taking the total mass of pollutants to the bottom at a much higher speed than during low-intensity discharge [6]. It happens at the expense of a combined action of several processes. First of all, the jet is destroyed due to the shear instability, turbulent pulsations, and other factors. Secondly, the jet falls down because its density is higher than the density of surrounding seawater, which is considered to be a "ball' effect. Thirdly, when the structure of the flow is destroyed with a decrease in the concentration of the solid phase in the slick, the effect of jet falling stops and only the settling of solid phase particles occurs according to the balance of gravitational and other forces.
On average, depending on the speed of currents, the rate of discharge and other factors, the duration of the "ball" effect is characterized by a range of 0.5-1.5 minutes. This effect manifests itself most significantly during a unit discharge.
Thus, to calculate the settling velocity of markers and , expressions (Eq. 07, 08) are used taking into account the effect of increasing the mass and equivalent radius due to the effect of flocculation and faster jet settling due to the "ball" effect.

Turbulence parameters
The turbulent pulsation velocity components in (Eq. 04) are determined using mean-square deviations , , , calculated from the experimental data of observations over currents or constructed using semi-empirical relations, and the mean values of the velocity modulus. Based on the hypothesis of normal distribution of the ocean turbulence pulsation spectrum [9], we can write down the expressions to determine the pulsation velocity components in (2.1) for each statistical test: (9) where -a random variable determined by the formula: A random variable takes equiprobable values of -1 or +1, has function (x)=0.5 and is determined on a set of two points: x={-1;+1}.
A random variable is distributed according to the normal law with parameters μ (expectation) and σ (mean-square deviation) and has a distribution function: (11) Let us apply the formula to obtain a random variable distributed according to the normal law (2.8) as per [10]: (12) where -a random variable distributed according to the normal distribution law that takes the values within the range . The velocity components of turbulent pulsations of the minimum scale in (2.1), that is comparable with the calculation discreteness in case of absence or lack of experimental data, are determined on the basis of the theory of shear instability that connects velocity components dispersion with Richardson criterion [11]. The relation for the calculation of the relations of mean-square deviation: Mathematic expectation of velocity modulus for relations (2.9) is determined through relations of mean-square deviation taking into account empiric-formula dependences obtained for large scales The exponential form for stable stratification [12] has proven itself the best, which is typical for most of the time period. So, we adopted the convention that .

The simulation results
As a result of the program operation the following results were obtained. Taking into account that initially the source was at the beginning of the system of coordinates under study, we can see that in 5 minutes suspended material consisting of 105 markers was shifted by 140 meters from its initial place and located within 140-160 meters. Moreover, figure 2.1 shows that suspended material is not homogenous and consists of the markers of different diameter, density and mass, which is reflected by their location in terms of depth: the more the density is, the higher the marker is located. In Figure 2.2, the dispersion of the suspended material of markers consisting of 105 is considered in 15 minutes. The graph shows that the suspended material was shifted by 420 meters and located within 420-480 meters.

Results and discussion
The result of the simulation of points with different diameter where the first point has diameter of 0.1-0.05 mm and, being on the bottom, is slowly moving to the surface ( Figure  1), the second point with diameter <0.005 mm is moving near the bottom (Figure 2), the third point with diameter >10 mm is located on the water surface (Figure 3), and the fourth point with diameter 2-1 mm is located on the water surface ( Figure 4). The graphs show that the trajectory of a particle depends on its diameter.

Conclusions
A three-dimensional diffusion-advective model was calculated for 105 markers within time intervals of 5 and 15 minutes was calculated in the computer algebra system Wolfram Mathematica. The equations describing the trajectories of the markers in the water layer are solved. Boundary conditions at the solid boundaries and the boundaries with air are determined. The settling velocities of markers under the action of gravitational force at the initial stage of jet destruction are calculated. Turbulence parameters are determined.
It is also worth noting that on the basis of the obtained data it will be possible to calculate: − average volumes of water layer pollution with suspended substances over the entire period of the source of impact operation, subject to pollution with different gradations of pollutant concentrations − average areas of water surface polluted with suspended substances over the entire period of the source of impact operation, with different gradations of pollutant concentrations − total volumes of polluted seawater for the entire period of operation − zones of settling of the solid phase of suspended substances on the seabed with the areas of sedimentation zones for different thicknesses of sediments − characteristics of quasisteady plumes of suspended substances in the water layer, including medium and maximum lengths as well as the time of plumes existence with different gradations of pollutant concentrations.