Two-dimensional-one-dimensional splitting scheme for numerical solution of problems of suspended matter transport in coastal systems

. This article discusses the problems of numerical solution of non-stationary convection-diffusion-reaction problems using the model problem of suspended matter transport as an example. In the difference scheme proposed by the authors, at each time layer, the original spatial-three-dimensional problem is split along horizontal directions into a chain of two-dimensional and one-dimensional problems. In order to ensure the unconditional skew-symmetry of the convective transfer operator and its energy neutrality, the convective terms are written in symmetric form (half the sum of the non-divergent and divergent forms). The approximation of the initial boundary value problem, to which the suspended matter transport model is reduced, is considered in the Hilbert space of grid functions, which in subsequent discussions will allow us to focus on the use of general results of the theory of stability (correctness) of operator-difference schemes.


Introduction
Coastal systems occupy an exceptional role in the life of the sea and represent a special area of application for human economic activity.As urbanization increases, anthropogenic pressure on the ecosystem of coastal areas increases.Without taking into account its influence, it is impossible to objectively assess the state of waters and effectively implement a policy of rational environmental management [1][2][3].
Knowledge of the spatial variability of suspended particles is an important prerequisite for describing and predicting environmental conditions in coastal zones.This knowledge is based on the simultaneous study of the dynamics and movement of water masses, the direction and intensity of currents, the characteristics of suspended matter, etc.To conduct research in this area, reliable quantitative forecasting methods are needed, based on modern non-stationary models in a complete formulation [4][5][6].An important role in the development of ideas about these processes at different times was played by the works of Soviet and Russian scientists, including Longinova V.V., Zenkovicha V.P., Antsyferova S.M., Kosyana R.D., Pykhova N.V., Pavlidis Yu.A., Aibulatova N.A., Aksenova A.A., Leontyeva I.O., Debolsky V.K., Safyanova G.A. and many others.They are joined by the works of foreign specialists Nielsen P., Horikawa K., Bruun P., Deigaart R., Hanes D.M., Vincent C.E., Huntley D.H., Clark T.L., Hayes M.O., Hom-ma M., Horikawa K., Kana T.W., Kantardgi I.G. and many others.The accumulation of new knowledge and experimental data encourages us to obtain new results on the problem that interests us.This paper presents a three-dimensional spatial model of suspended matter transport, which allows one to reliably describe the physics of processes and calculate suspended matter concentrations.For the numerical study, calculation schemes are constructed, obtained by splitting the original problem into a chain of two-dimensional problems in horizontal directions and one-dimensional problems in the vertical direction.Particular attention is paid to approximations in space of convective terms.Writing the convective terms in a symmetric form allows us to subsequently construct monotonic approximations (fulfillment of the maximum principle at a discrete level) and ensure the stability of schemes in the Hilbert space of grid functions [7,8].

Initial-boundary value problem of suspension transport
Let  ⊂ ℝ 3 is the area in which the process occurs, and which is a parallelepiped Closing an area G is carried out in an obvious way and consists of joining to internal points G all faces of the parallelepiped.We believe that in the region G there are suspended particles that are at the point   ,, x y z and at the moment t have concentration

 
, , , q q x y z t  ; t is a temporary variable.The equation describing the behavior of suspended particles will look like this [9,10]: The following notations are used here: ,, u v w is a components of the velocity vector U of the aquatic environment; g w is a hydraulic particle size of particles; , h   is a coefficients of horizontal and vertical turbulent diffusion of particles; F is a function of external sources.
It is assumed that the rate of descent of particles g w in an aquatic environment, with equal gravity and resistance, is constant, i.e. g w const  .
The solution to equation ( 1) is found in a certain given area arguments, which is a four-dimensional cylinder with generatrices parallel to the time axis Ot heights T . 1) we add the initial conditions at time 0 t  :

 
, , ,0 , , , , , ; and boundary conditions: where n is the outer normal to the boundary of the surface l  , Г U is the velocity vector of fluid motion at the surface l  boundary, q there are known concentration values; : 3 Construction of a difference scheme To equation (7) we add initial conditions of the form: q x y z q x y z q x y z t q x y z t n N x y z G  (9) : : .
Further, the set of internal grid nodes ,, x y z    we will denote accordingly as , , , as well as their Cartesian products as On a space-time grid    let us approximate problems ( 8) -( 13) on grids adopted in computational fluid dynamics with setting velocities (convective transfer coefficients) at grid nodes shifted by half a step: for u , grid Oz .This result can be arrived at using the integro-interpolation method (more generally known as the finite volume method).
Next, taking into account expression (13), we will construct an additive unconditionally stable splitting scheme along geometric directions, focusing on a chain of continuous problems of the form -two-dimensional and one-dimensional [7]: , .
Next, the symbol "-" is above the functions (1), ( 2), q q F F will mean that they belong to the class of grid functions.
For the convective transport operator, we have The approximation of terms describing diffusion processes has the form:  q x h y z q x y z Dq x h y z x h y z hh q x y h z q x y z q x y z q x h y z x y h z h h h q x y z q x y h z x y h z x y z h nn Taking into account the constructed approximations of the form ( 17) -(20), we come to the approximation of equations ( 14), (15) at internal grid nodes  -to implicit schemes of the form:

qq w x y z h q x y z h w x y z h q x y z h h q
x y z h q x y z x y z h hh To the system of difference equations it is necessary to add initial conditions for a chain of initial-boundary value problems of the form (8) for   ,, x y z   , as well as approximation of boundary conditions ( 9) - (12).If, with respect to the external normal, the components   u,v,w of the velocity vector U is less than zero, then the boundary conditions (10) are approximated for the surface l  in the following way: 0, , , 0,5 , , 0, , , , 0,5 , , 0, , , ; ,0, , ,0,5 , 0, , , , , 0,5 , 0, , , .
To set boundary conditions (10) on the surface l  for positive values of the normal component of the velocity vector, as well as for setting boundary conditions on surfaces For a given grid, grid  is internal.
We will assume that where \   there are boundary grid nodes   .
In addition, we will assume that the values of the velocity vector components of the aquatic environment and the hydraulic size of suspended particles at the grid nodes are known with fractional index values: , For those grid nodes \   , which are outside the reservoir, the values of the velocity vector components are assumed to be zero.Let us present separately the approximations of the convective terms for the nodes on each of the surfaces ,, i In the case of flows on the side faces of the area G , coinciding in direction with external normals to the surface l  , i.e., if the conditions are met 0,5 , , 0,5 , , 0, 0,5 , , 0,5 , , 0, , , ; ,0,5 , , 0,5 , 0, , 0,5 , , 0,5 , 0, , ,

y z u h y z u L h y z u L h y z x y z v x h z v x
Neumann boundary conditions apply.
Considering the composition of the surface l  , which is formed by dots which can be considered a difference approximation of the convective term at 0 x with an error   Oh : from which we have , , , , .
From ( 27) and ( 28) we obtain Arguing in a similar way, we get

1
, , 0,5 , , 0,5 , , . 2 For surface f  we will use boundary conditions (11) For surface f  we will use boundary conditions (12) of the third kind when approximating a one-dimensional (vertical) difference problem.We have its difference approximation with an error   On the other hand, we have the formal equality Expressing ,, from equality (34), we have and, substituting it into equality (35), we obtain Перейдем к аппроксимации диффузионных членов на границе области G .Let 0 x .Then formally on the extended grid   can be written down nn Taking into account the Neumann condition at the boundary 0 x if the first of conditions (25) is satisfied, we can write 2 0 1 0,5 , , 0,5 , , , , 0, , .
D q h y z h y z q h y z q y z h Likewise, for x xL  taking into account equality , , , , q L h y z =q L h y z  when the second of conditions ( 25) is satisfied, we obtain x Arguing in a similar way, we have for 0, ,0,5 , , 0,5 , , , ,0, ; Finally, when 0 z on a surface f  we have , ,0,5 , , 0,5 , , , ,0 .

D q x y h x y h q x y h q x y h
Note that in this relation   , , 0,5 0 vy x y h   (there is no turbulent diffusion outside the reservoir), so one should accept On a surface f  condition (12) must be taken into account.Let us formally write the expression for the diffusion term in difference form: and, using the difference boundary condition (36), we arrive at the relation In relation (44) the value of the coefficient   , , 0,5 is determined not only by the turbulent mixing of the aquatic environment at the bottom, but also by the physical properties of bottom sediments, i.e. particle size distribution, porosity, etc. Consideration of this issue is beyond the scope of this study and is not given here.
Regarding the values of the velocity vector components     0,5 , , , 0,5 , , , it is assumed that in expanding the area G in horizontal directions (in   ) there is an aquatic environment and these quantities can be determined in the hydrodynamic block of the combined model "hydrodynamics -sediments and suspensions".
The constructed splitting scheme is monotonic, and each of the difference equations is non-degenerate if the grid Peclet number is less than one.

Conclusions
This article considers the approximation of the three-dimensional problem of transport of suspended matter in coastal systems based on two-dimensional and one-dimensional splitting schemes.During approximation, splitting along horizontal directions into a chain of twodimensional and one-dimensional problems was performed at each time layer.The constructed splitting scheme is monotonic, and each of the difference equations is nondegenerate if the grid Peclet number is less than one.This was achieved by writing the convective terms in a symmetric form and choosing a sufficiently small time step.The numerical implementation of the resulting system of difference equations, with a non-self-adjoint operator (matrix) in the Hilbert space of grid functions, and which, as a rule, has a high dimension (10 6 -10 10 ), is a separate non-trivial problem.In a numerical implementation for a two-dimensional problem, iterative methods can be used.In particular, you can use the previously developed adaptive alternating triangular method [8], as well as the algorithmically simpler Seidel method.In this case, with the choice of space-time grid parameters described above and taking into account the guaranteed diagonal dominance of the system matrix operator, it is possible to ensure the convergence of the Seidel method at a geometric progression rate with a denominator of 0.7 -0.9.For a one-dimensional threepoint problem, the sweep method should be used.Additional advantages of this difference scheme appear when it is necessary to solve the problem in parallel for operational forecasts of the spread of pollution in coastal systems.Algorithms based on these schemes have lower time costs for performing exchanges between computing nodes of multiprocessor systems compared to traditional splitting schemes.Compared to full approximation schemes, these schemes, when implemented by iterative methods, require a significantly smaller number of iterations.However, consideration of these issues is beyond the scope of this article.

E3S
Web of Conferences 458, 03017 (2023) EMMFT-2023 https://doi.org/10.1051/e3sconf/202345803017 , bf  convenient to enter into the extended grid: of the area G .Let us first return to equalities (23), which can be refined taking into account the introduction of the grid   :