CC BY
7Original Russian Text © The Authors, 2020. published in MORSKOY GIDROFIZICHESKIY ZHURNAL, Vol. 36, Iss. 2 (2020)
Determination of Location of the Concentration Initial Field of a Possible Contamination Source in the Black Sea Water Area near the Gerakleisky Peninsula Based on the Adjoint Equations Method
S. G. Demyshev, O. A. Dymova, V. S. Kochergin, S. V. Kochergin*
Marine Hydrophysical Institute of RAS, Sevastopol, Russian Federation * vskocher@gmail.com
Purpose. The aim of the present article is to study propagation of contamination in the Sevastopol coastal zone using the dynamic model of high spatial discretization and which takes into account actual atmospheric forcing. Location of possible contamination sources is to be calculated and, hence, determined by applying the adjoint equations method using the measurement data assimilation. Methods and Results. The currents' field was calculated using the three-dimensional baroclinic model of the Black Sea water circulation developed in the Marine Hydrophysical Institute, RAS. The model is characterized by high degree of spatial discretization. The simulation included application of the SCIRON data on actual atmospheric forcing. The adjoint equations method was used for defying location of possible contamination sources. The model currents in the coastal zone were analyzed for the chosen time period. The admixture distribution from the initial contamination location in the Black Sea water area nearby the Gerakleisky Peninsula was calculated. Having been integrated, a series of the adjoint problems permitted to determine location of the contamination source. Conclusions. Based on the results of the numerical experiments on integrating a series of the adjoint problems, defined is the area where a possible contamination source can be located. The location of this area agrees well with the initial distribution of concentration that, in case the information on the admixture character is available, permits to define the place where it enters the marine environment. The proposed approaches can be used for solving various ecological problems as well as for interpreting and planning the field experiments on studying sewage propagation in the coastal zone.
Keywords: simulation, identification of contamination source, adjoint problem, admixture transport.
Acknowledgments: the study was carried out within the framework of the state task on theme No. 0827-2018-0004 "Complex interdisciplinary investigations of the oceanologic processes conditioning functioning and evolution of the Black and Azov seas' coastal zones".
For citation: Demyshev, G.S., Dymova, O.A., Kochergin, V.S. and Kochergin, S.V., 2020. Determination of Location of the Concentration Initial Field of a Possible Contamination Source in the Black Sea Water Area near the Gerakleisky Peninsula Based on the Adjoint Equations Method. Physical Oceanography, [e-journal] 27(2), pp. 210-221. doi: 10.22449/1573-160X-2020-2-210-221
DOI: 10.22449/1573-160X-2020-2-210-221
© S. G. Demyshev, O. A. Dymova, V. S. Kochergin, S. V. Kochergin, 2020 © Physical Oceanography, 2020
Introduction
The intensive development of the Azov-Black Sea region infrastructure in recent years has led to a significant increase in the load on its ecological system. Therefore, the development of computer systems for assessing the ecological state of the studied basins is an important and urgent task. The main elements of such systems are numerical models of water dynamics and contamination distribution models. The task of determining the location of contamination sources and their output is also an important problem.
The adjoint equations method applied in [1, 2] is used for solving this problem. Such an approach was first proposed in [1] when analyzing the function of meteorological parameters' impact on the formation of trends in geopotential, temperature and the field of vertical movements. Further, the solution of adjoint problem was successfully used to determine the "neutron cost" in the calculations of nuclear reactors [2]. In [3], the effect functions for the model of the atmosphere thermal interaction with the ocean and continents were first constructed. The results of this work were actively used in the planning of field experiments and expeditionary research in the oceans. A wide class of environmental problems based on the adjoint equations method is considered in the monograph [4]. An approach based on solving related problems provides the construction of the effect functions of transport model input parameters on the admixture concentration in the area under study and determination of the location of possible pollution sources. The effect functions indicate those areas that significantly affect the admixture concentration in a given area. Currently, the method has been successfully used for determining the sources of various types of contamination, including oil ones [5]. A similar approach has also been successfully used to solve various problems for the Black Sea and the Sea of Azov [6-9]. In this work, the adjoint equations method is used in the model of passive admixture transport in the Black Sea near the Crimean Peninsula in order to determine possible contamination sources in the coastal zone. The model of Marine Hydrophysical Institute (MHI) of the Russian Academy of Sciences [10, 11] with 1.6 km spatial resolution was applied as a dynamic model.
Thermohydrodynamic model
Numerical experiment is carried out by using three-dimensional non-linear MHI model. Model equations are written in the Cartesian coordinate system in the Boussinesq approximation, hydrostatics and incompressibility of the sea water:
ut - (4 + f )v + wuz = -gqx --(P1 + E)x + (vvuz)z-vHV4u , (1)
Po
Vt + (4 + f )u + wvz = -gqy - — (P1 + E)y + (VvVz) z - vH V4v, (2)
y Po y
z
P = gpo? + g j"p^=gpoq+p\ (3)
0
ux + vy + wz = 0, (4)
H
qt + J (ux + vy )dz = (Pr - Ev)/ (5)
0
Tt+(uT)x + (vT)y + (wT)z = -khV4T + (ktTz)z, (6)
St+(uS)x + (vS)y + (wS)z = -KH V4S + (KSSz)z, (7)
P = p(T, S), (8)
PHYSICAL OCEANOGRAPHY VOL. 27 ISS. 2 (2020) 211
where vh and kh are the coefficients of horizontal turbulent viscosity and diffusion, respectively; vv and kt' S are the coefficients of vertical turbulent viscosity and diffusion, respectively; Ev is evaporation; Pr is precipitation; the remaining notations are generally accepted. In the equations (1), (2)
dv du „ u2 + v2 TTT , ... . , .
q =---, E = p0-. We point out that the equation (5) is obtained under
dx dy 2
the assumption that the linearized kinematic condition is satisfied in Pr — Ev
the form w = —qt +--. Let M be an integration domain of the model on the
Pi
time interval [0, tL ], r be a solid boundary of the domain M ; M-t = M x[0, tL ], n
is a normal to r boundary.
For the system of equations (1) - (8), we write down the boundary conditions on the surface z = 0 in the form
PiVvUz =—Tx, PiVvVz =—Ty, PKTTZ = QT, KSSz = Pr—Ev So + P(Scl — So), (9)
Pi
where tx, ty a tangential wind stress; QT is a heat flux; So is a surface salinity; Scl is a climatic salinity; p1 is a sea water density in the upper layer; p is a relaxation parameter.
At the bottom at z = H(x, y)
u = 0, v = 0, w = 0, Tz = 0, Sz = 0. (10)
On the solid sidewalls: - for meridional sections of the boundary
u = 0, V2ux = 0, vx = 0, V2vx = 0,
Tx = 0, (V2 T) x = 0, Sx = 0, (V2 S) x = 0, - for zonal sections of the boundary
(11)
v = 0, V vy = 0, uy = 0, V uy = 0, Tv = 0, (V2T) v = 0, Sv = 0, (V2S) v = 0.
(12)
On the boundary sections where the water flows in, the Dirichlet conditions are used:
- for meridional sections
u = up, V2u = 0, V = 0, V2v = 0,
xx
T = Tp, S = Sp, (V2T)x = 0, (V2S)x = 0,
(13)
- for zonal sections
v = vp, V2vy = 0, uy = 0, V uy = 0,
2 2 (14)
T = Tp, S = Sp, (V2T)y = 0, (V2S)y = 0. 212 PHYSICAL OCEANOGRAPHY VOL. 27 ISS. 2 (2020)
For the Upper Bosphorus current and for the Kerch Strait, when the current is directed from the Black Sea to the Sea of Azov,
v = vs, V2vy = 0, uv = 0, V2uy = 0,
y y2 y 2 (15)
Ty = 0, Sy = 0, (V2T)y = 0, (V2S)y = 0.
In (13) - (15) the following notations were introduced: up, vp, vs are
the velocity in the estuaries and straits; Tp, Sp are temperature and salinity in the rivers. In the initial moment of time at t = 0 the following conditions are satisfied:
u = u 0( X, y, z), v = v0(x, y, z), g = g°( x, y),
0 0 (16)
T = T 0( x, y, z), S = S 0( x, y, z).
The coefficients of vertical turbulent exchange and vertical turbulent diffusion of heat and salt are calculated in accordance with the Mellor - Yamada 2.5 parameterization [12].
The finite-difference system of equations approximating equations (1) - (8) with the corresponding boundary (10) - (15) and initial (16) conditions is constructed on grid C according to the Arakawa classification [13]. The spatial resolution of the model is 1.6 km horizontally, 27 z-horizons are set vertically and the time step is 1.5 minutes. A detailed description of the model is given in [11].
The model of passive admixture transport
The equation of a passive admixture transport in Cartesian coordinates has the following form:
dC d(uC) d(vC) d(wC) , d „ dC
— + ^—- + ——- + —-- = AHV2C +—AV—, (17)
dt dx dy dz dz dz
where C is an admixture concentration; Ah is a coefficient of horizontal turbulent diffusion; Av is a coefficient of vertical turbulent diffusion.
Boundary condition on a free surface (z = 0) is set by the absence of flux:
Av ^ = 0. (18)
dz
On solid boundaries, at the bottom and in the area of river estuaries the conditions for the absence of a substance flux are set. In the Bosphorus Strait and the Kerch Strait regions the Dirichlet condition of the first kind is realized (in this calculation it is zero). At the initial moment of time the concentration field
C0 (x, y, z) is set.
The finite-difference form of differential equation (17) is written out on the grid C [13]. To approximate the advective terms, the TVD-scheme is used [14], the horizontal diffusion transport is parameterized by a harmonic operator with the coefficient Ah = 105 cm2/s. The values of vertical turbulent diffusion coefficient are predetermined vertically within the range of 2.5-0.03 cm2/s [15].
Adjoint problem
We associate the following adjoint [4] problem to the transport equation of the passive admixture (17) with the boundary conditions (18) and the initial data:
dC d(uC*) d(vC*) d(wC*) = A y2c, + C (i9) dt dx dy dz dz dz
z = 0: AV —— wC * = 0, V dz
z = H : AV — = 0
dz (20)
*
r: AH C = 0, dn
C* = h.
tL
In (20) h is an arbitrary function so far, it will be defined below; r is solid boundary of the region. In the region of estuaries and straits in these calculations C* = 0. Multiplying (17) by C* and integrating by parts, taking into account (4), (19) and boundary conditions, we obtain
J hCtdM = J C0C*dM . (21)
M M
We are to consider a certain domain Q e M, in which we are interested in certain functionals subject to control of admixture concentration. We choose h in the form
— within the domain Q, h = ^ m(Q) (22)
0 — outside the domain Q,
where m is a measure of the selected domain. Then, on the left side of expression (21), we obtain the average concentration Ct in Q at a finite point in time tL.
By choosing Q as a cell of a computational net, we can assess the concentration in it using the initial data and the corresponding solution to the adjoint problem. Thus, from (21) and (22) it can be seen that the solution of the adjoint problem C* is an effect function of the initial data on the value on the left side of expression (21). Therefore, by analyzing the spatial structure C0for various points of the computational domain, at which the concentration values significantly exceed the background values, we can determine the region of the contamination source possible location.
The results of numerical experiments
For the period of October 12-18, 2016, current velocity fields, obtained on the basis of the MHI thermohydrodynamic model with 1.6 km step with regard to real atmospheric forcing for 2016 according to SKIRON data with 0.1° spatial
resolution of [16], were selected. As the initial fields, the data of the Black Sea Marine Forecast Center were used [17]. In order to assess the sea dynamics effect on the problem results, we consider the features of currents in the analyzed period. The surface field of currents is shown in Fig. 1. At the initial moment of time (October 12), the Rim Current (RC) velocity off the Crimean coast reached 27 cm/s. Two anticyclonic gyres with ~ 50 km dimensions and up to 20 cm/s orbital velocities were observed westwards of Sevastopol. The velocity of currents in the vicinity of Sevastopol was low and did not exceed 6 cm/s. As a result of wind direction change from north-eastern to western, a current in the surface velocity field, which spread northwards along the coast from Khersones Cape to Tarkhankut Cape and was observed on October 13-15, occurred. The current velocity near the coast reached 25 cm/s. Westwards of Sevastopol the eddies weakened, and two new anticyclones arose on the left of the coastal current. The Rim Current in this period was characterized by a velocity of more than 40 cm/s. By October 14, an anticyclone with 25 km dimension was formed between the Rim Current and Khersones Cape. Over time, the orbital velocity in the anticyclone increased to 20 cm/s due to interaction with the Rim Current. With no resizing, the eddy moved towards the Rim Current motion and was observed until the end of the integration period.
5 10 15 20 25 30 35
F i g. 1. Currents' field (cm/s) on the sea surface on October 14, 2016
After October 16, a moderate northeastern wind prevailed in the area under study. An increase in the Rim Current was observed and the coastal current in the area of the Crimea western part weakened. Its velocity decreased to 5-10 cm/s. The structure of the velocity field northwestwards of Sevastopol was irregular. By October 18, the Rim Current velocity increased to 45 cm/s, which, apparently, again led to the appearance of a coastal current directed to the north near
the Khersones Cape. According to the modeling data, current velocity here was 25 cm/s.
Further, using the calculated current velocities as input ones, an experiment for simulating the impurity propagation for October 12-18, 2016 period was carried out. The initial impurity concentration field was set in the coastal part of the Gerakleisky Peninsula by the formed contamination spot, indicated in Fig. 2 by square. In this region, a unit conditional concentration was set from the surface to the bottom. The result of modeling the impurity concentration field on October 18 is shown in Fig. 3. It can be seen that pollution from the source for a given time period at the circulation that took place can propagate to a wide area. Considering the points on the periphery of the pollution area (white squares in Fig. 3) as coordinates for setting the initial data when solving the adjoint problems (19) -(22), we get a series of solutions given in Fig. 4-6. Calculations using the passive admixture transport model were performed on a grid with a spatial resolution of 1.6 km. The same spatial discretization was applied to solve the adjoint problem (19), (20).
F i g. 2. Concentration initial field and area ¥
In Fig. 4 the normalized values of initial data effect function on the impurity concentration at point 1 (see Fig. 3) are represented. It can be seen that the concentration of impurities at a given point (for a selected period of integration of the model equations and input parameters) is mainly affected by the concentration of impurities in the coastal region eastwards of Balaklava and in the more seaward zone located in the Rim Current passage area. If the first zone is characterized by the presence of an alongcoast current on October 13-15, then the second, in all probability, depends on the intensity of the anticyclonic formation that arose at that time southwards of the Gerakleisky Peninsula. A similar double structure of the effect functions is obtained for points 2-4.
I I I I
32(7 33. (7 34.0 35.0
F i g. 3. Model field of concentration
Black Sea
3Z0= 33. cf 34.0= 35.0= £
F i g. 4. Normalized effect function for point 1
The difference between the effect functions for points 4 and 5 is interesting. If the concentration at the point 4 located in the coastal zone is affected by the RC, then the concentration at point 5 is affected by the vast coastal zone from Foros to Yevpatoriya and the region of anticyclonic formation. When setting the initial perturbation at point 6 (Fig. 5), the dynamics of the anticyclonic eddy has a greater effect than the coastal zone. The study of the effect functions of the initial data on the concentration at the peripheral point 7 (Fig. 6) showed that the coastal zone and the region of anticyclonic formation have the main effect on the concentration at this point. A similar situation is obtained for point 8.
PHYSICAL OCEANOGRAPHY VOL. 27 ISS. 2 (2020) 217
—r--1—:-r--i—
32.0 33.0 34.0 35.0
F i g. 5. Normalized influence function for point 6
Black Sea
-r--1---I—-I—-
32.0" 33.0 34.0 35.0 £
F i g. 6. Normalized influence function for point 7
If the location of the initial pollution spot is to be determined, then by solving N of adjoint problems (19), (20), we can determine the region ¥ = fl^„, n = 1,...,
N, where ^n are the areas of essential values C* . In Fig. 2, such area ¥ is marked by dark gray. It can be seen that the coordinates of the initially specified contamination spot correlate well with the location of this region.
Conclusion
As a result of numerical experiments on integrating a series of adjoint problems, the region in which a possible source of pollution can be located is determined. The location of this region is in good agreement with the initial concentration distribution, which, if information on the impurity nature is available, provides the localization of the place where it enters the aquatic medium. The calculations showed that on the used time interval the concentration of impurities in the studied area of Sevastopol coastal zone is mainly affected by the values of the initial concentration field in the coastal zone itself and in the areas adjacent to the southern coast of the Gerakleisky Peninsula and the Balaklava Bay. The results of numerical experiments have shown that under certain meteorological and hydrological conditions prevailing in the Sevastopol region, the situations when pollution from treatment facilities and associated sources can reach the main urban coastal infrastructure in a fairly short time are possible. Therefore, the solution of such a problem is important from an environmental point of view. In the presence of necessary information to fill in the dynamic model, transport model and the associated problem, it becomes possible to determine the location of pollution sources and analyze their effect on the concentration of impurities in the coastal zone of the selected area. Thus, based on the integration of adjoint equations, it is possible to realize the task of assessing the controlled values of the impurity concentration in the sea region of interest and identifying the location of the initial contamination spot. Such information may be useful in making decisions for reducing the anthropogenic load on the Sea of Azov-Black Sea basin ecosystem.
REFERENCES
1. Buleev, N.I. and Marchuk, G.I., 1958. [Dynamics of Large-Scale Atmospheric Processes]. Trudy Instituta Fiziki Atmosfery AN SSSR, (2), pp. 66-104 (in Russian).
2. Marchuk, G.I., 1959. Numerical Methods for Nuclear Reactor Calculations. New York: Consultants Bureau, Inc., 295 p.
3. Marchuk, G.I. and Skiba Yu.N., 1976. Numerical Calculation of the Conjugate Problem for the Model of Thermal Interaction of the Atmosphere, Oceans and Continents. Izvestiya of Academy of Sciences, USSR, Atmospheric an Oceanic Physics, 12(5), pp. 459-469 (in Russian).
4. Marchuk, G.I., 1986. Mathematical Models in Environmental Problems. Amsterdam: Elsevier Science Publishers B.V., 216 p.
5. Skiba, Y.N. and Parra-Guevara, D., 2017. Application of Adjoint Approach to Oil Spill Problems. Enviromental Modeling & Assessment, 22(4), pp. 379-395. doi:10.1007/s10666-016-9540-4
6. Ruabtzev, Yu.N. and Shapiro, N.B., 2009. Determination of the Initial Position Discovered in Open Part of Sea Lenses of Small Salinity. In: MHI, 2009. Ekologicheskaya Bezopasnost' Pribrezhnykh i Shel'fovykh Zon i Kompleksnoe Ispol'zovanie Resursov Shel'fa [Ecological Safety of Coastal and Shelf Zones and Comprehensive Use of Shelf Resources]. Sevastopol: ECOSI-Gidrofizika. Iss. 18, pp. 141-157 (in Russian).
7. Kochergin, V.S. and Kochergin, S.V., 2017. Identification of the Parameters of the Instantaneous Point Pollution Source in the Azov Sea Based on the Adjoint Method. Physical Oceanography, (1), pp. 62-67. doi:10.22449/1573-160X-2017-1-62-67
8. Demyshev, S.G., Kochergin, V.S. and Kochergin, S.V., 2018. Using the Variational Approach and Adjoint Equations Method Under the Identification of the Input Parameter of the Passive Admixture Transport Model. In: V. Karev, D. Klimov and K. Pokazeev, eds., 2018. Physical and Mathematical Modeling of Earth and Environment Processes. PMMEEP 2017. Cham: Springer, pp. 51-61. https://doi.org/10.1007/978-3-319-77788-7_7
9. Kochergin, V.S., 2011. Determining the Concentration Field with Use the Data Measurements Based on Solving of Adjoint Problems. In: MHI, 2014. Ekologicheskaya Bezopasnost' Pribrezhnykh i Shelfovykh Zon i Kompleksnoe Ispol'zovanie Resursov Shel'fa [Ecological Safety of Coastal and Shelf Zones and Comprehensive Use of Shelf Resources]. Sevastopol: ECOSI-Gidrofizika. Iss. 25(2), pp. 370-376 (in Russian).
10. Demyshev, S.G., 2012. A Numerical Model of Online Forecasting Black Sea Currents. Izvestiya, Atmospheric and Oceanic Physics, 48(1), pp. 120-132. https://doi.org/10.1134/S0001433812010021
11. Demyshev, S.G. and Dymova, O.A., 2013. Numerical Analysis of the Mesoscale Features of Circulation in the Black Sea Coastal Zone. Izvestiya, Atmospheric and Oceanic Physics, 49(6), pp. 603-610. https://doi.org/10.1134/S0001433813060030
12. Mellor, G.L. and Yamada, T., 1982. Development of a Turbulence Closure Model for Geophysical Fluid Problems. Reviews of Geophysics, 20(4), pp. 851-875. doi:10.1029/RG020i004p00851
13. Arakawa, A. and Lamb, V.R., 1981. A Potential Enstrophy and Energy Conserving Scheme for the Shallow Water Equations. Monthly Weather Review, 109(1), pp. 18-36. https://doi.org/10.1175/1520-0493(1981)109<0018:APEAEC>2.0.CO;2
14. Harten, A., 1997. High Resolution Schemes for Hyperbolic Conservation laws.
Journal of Computational Physics, 135(2), pp. 260-278. https://doi.org/10.1006/jcph.1997.5713
15. Demyshev, S.G., Zapevalov, A.S., Kubryakov, A.I. and Chudinovskikh, T.V., 2001. Evolution of the 137Cs Concentration Field in the Black Sea after the Passage of the Chernobyl Cloud. Meteorologiya i Gidrologiya, (10), pp. 49-61 (in Russian).
16. Kallos, G., Nickovic, S., Papadopoulos, A., Jovic, D., Kakaliagou, O., Misirlis, N., Boukas, L., Mimikou, N., Sakellaridis, G., Papageorgiou, J., Anadranistakis, E. and Manousakis, M., 1997. The Regional Weather Forecasting System SKIRON: An Overview. In: University of Athens, 1997. Proceedings of the Symposium on Regional Weather Prediction on Parallel Computer Environments, Vol. 15. Athens, Greece: University of Athens Greece, pp. 109-122.
17. Korotaev, G.K., Ratner, Y.B., Ivanchik, M.V., Kholod, A.L. and Ivanchik, A.M., 2016. Operational System for Diagnosis and Forecast of Hydrophysical Characteristics of the Black Sea. Izvestiya, Atmospheric and Oceanic Physics, 52(5), pp. 542-549. https://doi.org/10.1134/S0001433816050078
About the authors:
Sergey G. Demyshev - Chief Research Associate, Head of Wave Theory Department, Marine Hydrophysical Institute of RAS (2 Kapitanskaya Str., Sevastopol, 299011, Russian Federation), Dr. Sci. (Phys.-Math.), ORCID ID: 0000-0002-5405-2282, demyshev@gmail.com
Ol'ga A. Dymova - Senior Research Associate, Marine Hydrophysical Institute of RAS (2 Kapitanskaya Str., Sevastopol, 299011, Russian Federation), Ph. D. (Phys.-Math), SPIN-code: 796-208, ORCID ID: 0000-0003-4036-2447, ResearcherID: P-9669-2015, olgadym@yahoo.com
Vladimir S. Kochergin - Junior Research Associate, Marine Hydrophysical Institute of RAS (2 Kapitanskaya St., Sevastopol, 299011, Russian Federation), ORCID ID: 0000-0002-6767-1218, ResearcherID: AAG-4209-2020, vskocher@gmail.com
Sergey V. Kochergin - Senior Research Associate, Marine Hydrophysical Institute of RAS (2 Kapitanskaya St., Sevastopol, 299011, Russian Federation), Ph. D. (Phys.-Math.), ORCID ID: 0000-0002-3583-8351, ResearcherID: AAG-4206-2020, ko4ep@mail.ru
Contribution of the co-authors:
Sergey G. Demyshev - statement of the problem, analysis of the results, preparation of the paper text
Ol'ga A. Dymova - carrying out numerical experiments, analysis and interpretation of research results
Vladimir S. Kochergin - carrying out numerical experiments, analysis of the study results, preparation of the paper text
Sergey V. Kochergin - problem statement, realization of adjoint model, analysis of the calculation result, writing and updating the paper text, discussion of the results, formulation of conclusions
All the authors have read and approved the final manuscript.
The authors declare that they have no conflict of interest.
