CC BY
6Original Russian Text © S.V. Kochergin, V.V. Fomin, 2019 published in MORSKOY GIDROFIZICHESKIY ZHURNAL, Vol. 35, Iss. 6 (2019)
Variational Identification of Input Parameters in the Model of Distribution of the Pollutants from the Underwater Source
S. V. Kochergin *, V. V. Fomin
Marine Hydrophysical Institute of RAS, Sevastopol, Russian Federation * vskocher@gmail.com
Purpose. The aim of the paper is to construct and to validate the variational procedure for determining the pollutant concentration and the water flow out velocity at the underwater source exit, as well as to analyze the algorithm sensitivity to the level of random noise in the measurement data. Methods and Results. The flow field was calculated using the three-dimensional baroclinic a-coordinate model of water circulation. Realization of the pollution transfer model included application of the TVD-type monotone schemes. The temperature and salinity initial profiles were preset based on the results of probing in the area of the underwater release, and the characteristic velocity of the background currents was defined using the data of the ADCP-measurements. The input parameters of the problem were identified by means of the iterative procedure for minimizing the quadratic functional. The numerical experiments on identifying parameters of the underwater pollution source showed that if noise was left out of account, the original parameters were reconstructed with a relative error < 1%. It is shown that the identification problem becomes of better conditionality in case the data from more informative points of the measurement scheme are assimilated.
Conclusions. Based on the analysis of the numerical experiments, the linearization algorithm is shown to be able to identify the parameters of the underwater source. The proposed algorithms can be used to solve a wide class of environmental problems, as well as to interpret and to plan the field experiments aimed at studying the wastewater distribution in the coastal waters.
Keywords: functional minimization, numerical modeling, parameter identification, problem in variations, assimilation of measurement data, linearization method.
Acknowledgements: the work was carried out within the framework of the state task on theme No. 0827-2018-0004 "Complex interdisciplinary studies of oceanological processes conditioning functioning and evolution of the ecosystems of the Black and Azov seas' coastal zones" and at partial support of the RFBR grant No. 18-45-920035 p_a.
For citation: Kochergin, S.V. and Fomin, V.V., 2019. Variational Identification of Input Parameters in the Model of Distribution of the Pollutants from the Underwater Source. Physical Oceanography, [e-journal] 26(6), pp. 547-556. doi:10.22449/1573-160X-2019-6-547-556
DOI: 10.22449/1573-160X-2019-6-547-556
© 2019, S. V. Kochergin, V. V. Fomin © 2019, Physical Oceanography
Introduction
Various environmental factors have a significant effect on the environmental situation in marine coastal waters, including the release of pollutants from underwater reservoirs. Availability of complete and reliable information about the situation in the water area of discharges is important for assessing the consequences and making managerial decisions. To obtain it, it is necessary to use contact and remote measurement methods [1-3], as well as the use of mathematical modeling methods [4-6].
Using mathematical modeling, the problem of identifying model parameters from measurement data naturally appears. Methods of determining the optimal parameters are most often based on minimizing the cost functionals, characterizing the difference between the model estimates of the impurity concentration and the measurement data. Assimilation algorithms [7-9] are constructed using adjoint ISSN 1573-160X PHYSICAL OCEANOGRAPHY VOL. 26 ISS. 6 (2019) 547
(® The content is available under Creative Commons Attribution-Noncommercial 4.0 TTT^ International (CC BY-NC 4.0) License
equations [10]. Such procedures are especially effective when the required parameters are variable in space and time. When some numerical modeling constants are subject to identification, one of the effective procedures for finding the optimal parameters is the estimation method based on the linearization [11, 12]. In this case, the problem in variations is solved once, which can sometimes be very effective for iterative identification of parameters.
The present paper is aimed to construct and validate the variational procedure for determining input parameters of a transfer model. A numerical hydrothermodynamic model is considered [13], in which the parameters to be identified are the pollution concentration in the underwater source (hereinafter Cp) and the outflow velocity at the release from it (hereinafter wp). It should be noted
that the identification of such parameters is an important task, since obtaining their real values is usually problematic.
Dynamical model
The nonlinear baroclinic model described in [13-15] is used to calculate the flow fields, impurity concentration, temperature and salinity of sea water in the area of underwater outlet. It is based on three-dimensional o-coordinate equations of ocean dynamics in the Boussinesq and hydrostatics approximation of the following form (summation over the indices a and P, which hereinafter change from 1 to 2, is assumed):
f ( Du. ) + A„. +E.p JDult +gD |l + DB. =A (^ ) + d0 (
d ( KM dUa
— +-(Du.) +-= 0,
dt cx do
dt
dt -(
dt
A9 = dt( Dup9 )+^( ),
d ( dT ^ at I dxP J d
dxp | do
d \ dS )
dxp dxn V P J + do
d dxp ( a dC' Ac — dxP \ p j d + — do
P = P (t , s ),
do I D do
do \ D do
f
D do
do
B =
a
Po
VdXa
D
•o
j pdo'
\
dD o-p
dXa J
(1) (2)
(3)
(4)
(5)
(6) (7)
t„„ = 2 A.
dUa dx„
_ AM
(dUP + dUa
dxa
(8)
where x\ = x, X2 = y; a - dimensionless vertical coordinate varying within -1 to 0; D = h0 + n - dynamic depth; xap - the turbulent stress tensor components; AM,
KM , AT, KT, AS, KS , AC, KC - coefficients of turbulent viscosity and diffusion; g - free fall acceleration; f - the Coriolis parameter; p0 - mean water density; eap=0 when a = P; e12 = -1; e21 = 1, C - contaminating impurity concentration; T, S - water temperature and salinity.
Sought variables of the system (1) - (8) are (u1,u2 ) = (u,v) - horizontal velocity components, w* - flow velocity component normal to the surfaces a = const, fields T, S, C and sea water density field p .
On the free surface (a = 0) the boundary conditions are of the following form
w* = 0, ^^ = 0, = 0, ^ = 0, = 0. (9)
D da D da D da D da
The boundary conditions at the bottom outside the source (a = -1, x ^ xp,y ^ yp) are written as follows:
KM 5un | | KT dT n KS dS n KC dC n .iri. w* = 0, —- = u uu, —— = 0, —— = 0, —— = 0, (10) D <3a l|a D a D a D a
where ^ - bottom friction coefficient; |u| =
-u 2
2 '
The boundary conditions at the bottom in the source action scope ( o = -1, x = xp, y = yp ) can be represented in the following form [14]:
= ' % IT = H u|Ua , (11)
D do
T-K^<K=wT wS-K^dS = wS wC-K^dc=wC pT D do plp ' ^ D do pSp ' pC D do pCp ■
At the initial moment of time, three-dimensional hydrothermodynamic fields and the concentration field of the pollutant are described by the following conditions:
u = U0, v = w = 0, n = 0 , T = T0 (a), S = S0 (a), C = 0. (13)
Here u, v, w - components of the flow velocity along the axes x, y, a correspondingly; U0 - depth constant background flow velocity; T<)(a) -background temperature distribution; S0(a) - background salinity distribution; wp = Qp/d2 - water outflow rate, Qp - water consumption; d - horizontal source
size; Tp , S p - temperature and salinity of outflowing water; Cp - impurity
concentration at the source exit.
The numerical procedure for solving the system of equations is based on the separation of the barotropic and baroclinic components in the velocity field and the use of explicit schemes for horizontal coordinates and implicit schemes for a-coordinate. The coefficients of horizontal turbulent diffusion are calculated by PHYSICAL OCEANOGRAPHY VOL. 26 ISS. 6 (2019) 549
the Smagorinsky formula [16]. The vertical turbulent viscosity and diffusion coefficients are determined by the Maylor - Yamada model [13, 17]. To approximate the advective terms in the model equations, the TVD schemes are used [18, 19]. They have the property of monotonicity, which is an important condition for adequate modeling of gradient fields. The algorithm for the numerical implementation of the model is described in detail in [14].
The problem is solved on a time interval [0,10 ] in a rectangular region Q = {0 < x < L; 0 < y < L; -1 < o <0} with liquid lateral borders and a free surface. When t > 0 at the basin bottom (a = -1) in the local area Q the source of mass begins to act.
Algorithm of the pollution concentration identification in the source
Consider the case when Cp is selected as the parameter to be identified. Such
dC
a problem can be solved based on the linearization method [11]. Let V =--is
dCp
C-function variation relative to the required parameter Cp. According to [20], for the equation (5) with boundary conditions from (9) - (13), the following problem is written in variations below:
d __ _ d
( dV ^
dt d x
— (DV) + AV = — Ac— +— AC —
d y \ d y
\ d x j
d
( . dV ^ d (Kr dV
do \ D do
(14)
Kc — = 0 when o = 0, (15)
D do
K dV
— = 0 when o = -1 and (x, y) gQs, (16)
D do
Kc- — = wp (V -1) when o = -1 and (x, y) eQp. (17)
D do
Let the problem (14) - (17) is solved within the time interval [0,10 ] in
the region Q. Then, to solve the problem of identifying the parameter Cp, it is necessary to find the minimum of the following functional:
J = \{P(RC - Cobs),P(RC - Cobs)), (18)
*0
where {a,b) = j jjjabdQdt - scalar product; Cohs - measured values of C at given
0 Q
points in the region Q at certain time points; R - projection operator in the observation point; P - operator of zero filling of the residual field in the absence of measurement data.
Suppose the variable C in the form
C = C + V(Cp - C*), (19)
where C - some assessment of pollution concentration, and Cp - its true meaning to be identified. After substituting (19) into (18) we obtain
J = 2(p(R(C+VC -cp))-),P(R(C+V(cp -C]))-Cofe)) . (20)
Further for simplicity we drop the dash above c . From the stationary
condition -dJ- = 0 we have dCp
= C
p p (PRV, PRV) '
Formula (21) is also valid for the case when there are measurements, for example, only on the sea surface. Then the problem of constructing the cost functional (18) is solved by the corresponding choice of the operators P and R. Identification algorithm of the parameter Cp is implemented iteratively. First,
the problem is solved once in variations (14) - (17) and the variable V is determined. Next, a certain initial approximation for the parameter Cp is specified
and the main problem (1) - (13) is solved. Subsequently, Cp is estimated by
the following formula
C- = C-1 + ^ (^ - ^ ) ,PRV) (22) p p (PRV, PRV) '
where n = 1, 2, .., N - iteration number.
Algorithm of source outflow velocity identification
Outflow velocity identification wp within the system (1) - (13) is significantly
more complicated than Cp identification, since it leads to a nonlinear problem in variations. This problem, even in the approximation of constant coefficients AM, KT, KS, KC is quite complex to solve. Therefore, the following simplified
dCp .
identification algorithm is proposed. At the first stage, the variation V = —- is
CWp
determined by averaging over the ensemble of solutions to a series of basic problems for various values of wp in the vicinity of the true value wp = w*p. At the second stage, the initial approximation is set for the outflow velocity and the value wp is refined by the following formula
(P(RCn -Cobs),PRV)
w"p = w;-1 -obs)\ ' . (23)
p p (PRV, PRV)
dJ
Expression (23) is obtained from the stationarity condition -= 0 for
dWp
the functional (18), where the value C is represented as C = C + V(wp - w"p).
Numerical experiments and discussion of the results
Model calculations were carried out for the area of the Golubaya Bay of Sevastopol, where the main outlet of urban wastewater is located. The computational domain had a horizontal size L = 2 km and a depth of h0 = = 30 m. A rectangular grid with a horizontal step of d = 20 m and a resolution of a 1 m was used. The time step of model integration was chosen to be 5 s. At the western boundary of the computational domain (x = 0), conditions of the form (1) were specified. On other liquid boundaries, for all variables, smooth continuation conditions were used: cfy/dx = 0 for x = L, 39/dy = 0 for y = 0 and y = L.
The pollution source center was located at a point with coordinates xp = 600 m, yp = 1000 m. The horizontal size of the source along the x and y axes was equal to the computational grid step. It was believed that the water flowing from the source has zero salinity (Sp = 0), and its temperature is equal to the ambient temperature, that is, Tp = T0 (-1). As T0 and S0, the results of sounding in the underwater release area on May 20, 2016 were used. According to the ADCP measurements in the area of the Heraclean Peninsula [21], the value of U0 was assumed to be 0.05 m/s.
First, a model calculation of the dynamic characteristics and values of the concentration of contaminants caused by the source action was carried out with a given "true" value C* = 0,5 kg/m3. The total integration time was 6 hours. During
the first hour of the model time, the water flow Qp in source linearly increased from 0 to 1.4 m3/s and then remained unchanged.
The spatial distribution of the polluted water field over a fixed horizon a = -0.3166 (z = -9.5 m) for t = 6 h is shown in Fig. 1. Here are the values of the relative concentration c = 100%-C/ C*p. The dimensionless horizontal coordinates have the form: x' = x/2d ; y' = y/2d . In Fig. 2, for t = 6 h, values of relative concentration c are shown in the section along the axis x' when y' = 25 . It can be seen from the figures that the field of polluted waters consists of a torch and a jet elongated in the direction of the background current and concentrated in the layer of the density jump (8-10 m). A detailed description of the effect of sea water density stratification and the background flow velocity on the impurity field is given in [14, 15].
Fig. 1. Horizontal distribution of the relative impurity concentration c (%) in the area of the underwater source at a = -0.3166 (z = -9.5 m)
Fig. 2. Vertical structure of the relative impurity concentration c (%) field in the area of the underwater source at y' = 25
In parallel with the "model" calculation, the problem was solved in variations (14) - (17). When solving it, the necessary vertical profiles V were preserved at the corresponding time moments of the measurements.
In the process of model integration, the values of concentration C were retained. It was assumed that some moving sensor sequentially measures the vertical profile C from the free surface to a depth of z = 25 m (aB = -0.8333). Measurements begin at time t\ = 2 hours in the point (x', y[) = (25, 20) and finish at time t15 = 5.5 hours in the point (xf5, y''5 ) = (35, 30). Time resolution of
measurements is 15 min. The indicated measurement scheme is selected from the considerations that in field studies, the identification of contaminated water is usually carried out by carrying out a series of vertical soundings perpendicular to the direction of jet propagation [3].
Next, a series of numerical experiments to test the above Cp parameter identification algorithm. Numerical experiments were carried out with a different number of stations. In addition, random noise £ e[- r/ 2, rj 2 ] was added to
the "measured" vertical profiles C(a), where, r is the given noise level, S e[ 0,1 ] is a normally distributed random variable. As an initial approximation, a value of Cp = 0.1 kg/m3 was set. It was established that for the convergence of the iterative process, 2-3 iterations are necessary.
The results of parameter recovery for various combinations of stations at 0 < r < C* /10 are given in Table 1. According to the comparison, without taking
into account random noise, the parameter C p is restored accurately, regardless of
the number and information content of the stations used. When random noise is taken into account, the best recovery occurs when the most informative stations are used (group A), indicated in Fig. 1 by white squares. Group B includes station contamination spots located on the periphery, the coordinates of which are indicated by black circles. Such a result agrees well with the work [22] devoted to the planning of experiments.
When testing the wp parameter identification algorithm, a series of numerical
experiments was carried out. For simplicity, it was believed that measurements of all 15 vertical profiles C, located according to the diagram in Fig. 1, were carried out simultaneously at t = 4 hours. Impurity concentration in the source C p = 0.5 kg/m3.
T a b l e 1
Results of reconstructing the pollutant concentration Cp at the underwater
source exit for different groups of the stations at different noise level r
in the input data
Stations r, kg/m3 Noise level, 100% • r/Cp ,% Reconstructed concentration Cp, kg/m3 Reconstruction error, 100 • |1 - CpjCp|, %
0.0000 0.0 0.5000 0.0
0.0125 2.5 0.4791 4.2
Group A
0.0250 5.0 0.5814 16.3
0.0500 10.0 0.8432 68.6
0.0000 0.0 0.5000 0.0
0.0125 2.5 0.5836 16.7
Group B
0.0250 5.0 0.6932 38.6
0.0500 10.0 1.0986 119.7
First, 8 main problems were solved (1 < j < 8 ), in which water consumption in the source Q varied from 0.8 to 2.0 m3/s with a step NQp =0.2 m3/s. The corresponding values of the outflow velocity were determined by the formula wpj = Qpj/d2. The "model" water consumption value Q* =1.4 m3/s, and
p
the corresponding "model" value of the outflow velocity w"p = 3.5-10 3 m/s.
554 PHYSICAL OCEANOGRAPHY VOL. 26 ISS. 6 (2019)
In the process of solving these problems, the fields C were saved at t = 4 hours. Then, for each wp j the difference derivative was calculated
V = Cv (wp,j) - Cp (wp,j-)/Awp , where Awp = AQp/d2; 2 < j < 8. After the mean
value V = 1/ l^Vj was determined. It was taken for variation in the formula (23).
j
As an initial approximation, wp = 10-3 m/s was used
The results of the reconstruction at different noise levels r in the input data are given in Tab. 2. As can be seen, at a noise level of 1.25% or lower, the value is reconstructed almost exactly. As with identification of Cp, the greatest decrease
required parameter is due to the assimilation of information from all measurement stations (including non-informative ones) and the approximate nature of determining the concentration variation. Naturally, when assimilating the most informative data, the accuracy of parameter reconstruction is increased due to an improvement in the conditionality of the solved identification problem.
T a b l e 2
Results of reconstructing velocity of the water release wp from the underwater source at different noise level r in the input data
r, kg/m3 Noise level, 100% • rjCp , % Reconstructed water release velocity wp, m/s Reconstruction error, 100 • 1 - wp!w'p\, %
0.0000 0.0 3.4991985-10-3 0.0
0.00675 1.25 3.5180813-10-3 0.5
0.0125 2.5 3.8568601-10-3 10.2
0.0250 5.0 4.6490142-10-3 32.8
Conclusions
Based on the analysis of numerical experiments, the applicability of the linearization algorithm for identifying the input parameters of the model of the distribution of pollutants from an underwater source is shown. The algorithms have good convergence of the iterative process, which permits quick evaluation of the pollution source concentration. The proposed algorithms can be used to solve a wide class of environmental problems, as well as to interpret and plan field experiments to study the distribution of wastewater in coastal waters.
REFERENCES
1. Bondur, V.G. and Grebenyuk, Yu.V., 2001. Remote Indication of Anthropogenic Influence on Marine Environment Caused by Depth Wastewater Plum: Modelling, Experiments. Issledovanie Zemli is Kosmosa, (6), pp. 49-67.
2. Bondur, V.G., 2005. Complex Satellite Monitoring of Coastal Water Areas. In: ISRSE, 2005. Proceedings of 31st International Symposium on Remote Sensing of Environment, ISRSE, 2005: Global Monitoring for Sustainability and Security. Saint Petersburg, 7 p. Available at: http://www.aerocosmos.info/pdf/1/2005_31_ISRSE_Bondur.pdf [Accessed: 20 January 2020].
3. Bondur, V.G., 2011. Satellite Monitoring and Mathematical Modelling of Deep Runoff Turbulent Jets in Coastal Water Areas. In: F. S. Garcia Einschlag, ed., 2011. Waste Water -Evaluation and Management. Croatia: InTech, pp. 155-180. http://dx.doi.org/10.5772/16134
4. Blumberg, A.F., Ji, Z.-G. and Ziegler, C.K., 1996. Modeling Outfall Plume Behavior Using Far Field Circulation Model. Journal of Hydraulic Engineering, 122(11), pp. 610-616. http://dx.doi.org/10.1061/(ASCE)0733-9429(1996)122:11(610)
5. Zhang, X.-Y. and Adams, E.E., 1999. Prediction of Near Field Plume Characteristics Using Far Field Circulation Model. Journal of Hydraulic Engineering, 125(3), pp. 233-241. https://doi.org/10.1061/(ASCE)0733-9429(1999)125:3(233)
6. Bondur, V.G., Zhurbas, V.M. and Grebenyuk, Yu.V., 2006. Mathematical Modeling of Turbulent Jets of Deep-Water Sewage Discharge into Coastal Basins. Oceanology, 46(6), pp. 757-771.
7. Marchuk, G.I. and Penenko, V.V., 1979. Application of Optimization Methods to the Problem of Mathematical Simulation of Atmospheric Processes and Environment. In: G. I. Marchuk, ed., 1979. Modelling and Optimization of Complex Systems. Lecture Notes in Control and Information Sciences, vol. 18. Berlin: Springer, pp. 240-252. https://doi.org/10.1007/BFb0004167
8. Marchuk, G.I., Agoskov, V.I. and Shutyaev, V.P., 1996. Adjoint Equations and Perturbation Algorithms in Nonlinear Problems. Boca Raton: CRC Press, 288 p.
9. Shutyaev, V.P., Le Dimet, F.-X. and Parmuzin, E., 2018. Sensitivity Analysis with Respect to Observations in Variational Data Assimilation for Parameter Estimation. Nonlinear Processes in Geophysics, [e-journal] 25(2), pp. 429-439. https://doi.org/10.5194/npg-25-429-2018
10. Marchuk, G.I., 1995. Adjoint Equations and Analysis of Complex Systems. Mathematics and Its Applications, Vol. 295. Dordrecht: Springer, 468 p. https://doi.org/10.1007/978-94-017-0621-6
11. Alifanov, O.M., Artyukhin, E.A. and Rumyantsev, S.V., 1995. Extreme Methods for Solving Ill-Posed Problems with Applications to Inverse Heat Transfer Problems. New York: Begell House, 306 p.
12. Gorsky, V.G., 1984. [Planning of Kinetic Experiments]. Moscow: Nauka, 241 p. (in Russian).
13. Ivanov, V.A. and Fomin, V.V., 2016. Numerical Simulation of Underwater Runoff Propagation in the Heraklean Peninsula Coastal Zone. Physical Oceanography, [e-journal] (6), pp. 82-95. https://doi.org/10.22449/1573-160X-2016-6-82-95
14. Ivanov, V.A. and Fomin, V.V., 2010. Mathematical Modeling of Dynamic Processes in the Sea-Land Area. Kiev: Akademperiodyka, 286 p.
15. Bondur, V.G., Ivanov, V.A. and Fomin, V.V., 2018. Peculiarities of Polluted Water Spreading from a Submarine Source in Stratified Coastal Environment. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 54(4), pp. 386-393. https://doi.org/10.1134/S0001433818040205
16. Smagorinsky, J., 1963. General Circulation Experiments with the Primitive Equations: I. The Basic Experiment. Monthly Weather Review, [e-journal] 91(3), pp. 99-164. https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2
17. Mellor, G.L. and Yamada, T., 1982. Development of a Turbulence Closure Model for Geophysical Fluid Problems. Review of Geophysics, [e-journal] 20(4), pp. 851-875. http://dx.doi.org/10.1029/RG020i004p00851
18. Harten, A., 1984. On a Class of High Resolution Total-Variation-Stable Finite-Difference Schemes. SIAM Journal on Numerical Analysis, [e-journal] 21(1), pp. 1-23. https://doi.org/10.1137/0721001
19. Fomin, V.V., 2006. Application of TVD Schemes to Numerical Modeling of Frontal Salinity Zones in a Shallow Sea. Russian Meteorology and Hydrology, (2), pp. 39-47.
About the authors:
Sergey V. Kochergin - Senior Research Associate, Marine Hydrophysical Institute of RAS (2 Kapitanskaya St., Sevastopol, 299011, Russian Federation), Ph. D. (Phys.-Math.), ORCID ID: 00000002-3583-8351, vskocher@gmail.com
Vladimir V. Fomin - Head of Department of Computer Engineering and Mathematical Modeling, Marine Hydrophysical Institute of RAS (2 Kapitanskaya St., Sevastopol, 299011, Russian Federation), Dr. Sci. (Phys.-Math.), ORCID ID: 0000-0002-9070-4460, v.fomin@ukr.net
Contribution of the co-authors:
Sergey V. Kochergin - problem statement, formulation of variational algorithm, analysis of the calculation result, writing and updating the paper text, discussion of the results, formulation of conclusions
Vladimir V. Fomin - problem statement, adjustment of the numerical model, carring out of the numerical calculations, editing of 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.
