УДК 519.6
One-dimensional Model for Studying Seasonal Changes of Vertical Structure of Salt Lake Uchum
Victor M. Belolipetskii* Svetlana N. Genova^
Institute of Computational Modelling SB RAS Akademgorodok, 50/44, Krasnoyarsk, 660036
Russia
Andrey G. Degermendzhy* Vladimir V. Zykov§
Institute of Biophysics SB RAS Akademgorodok, 50/50, Krasnoyarsk, 660036
Russia
Denis Yu. Rogozin^
Institute of Biophysics SB RAS Akademgorodok, 50/50, Krasnoyarsk, 660036 Siberian Federal University Svobodny, 79, Krasnoyarsk, 660041
Russia
Received 15.06.2018, received in revised form 08.09.2018, accepted 10.11.2018 Many salt lakes are meromictic, in which the water column is not mixed to the bottom for at least one year. In the stratified lake, the upper (epilimnion) and deep (hypolimnion) layers are distinguished, in which the density gradients are small. Between them is a layer of water (metalimnion), within which the density gradient is great. In the near-bottom layer, hydrogen sulphide accumulates and there is no oxygen. One-dimensional mathematical models are used to determine the dynamics of the vertical thermohaline structure of the salt lake Uchum. Two periods in the annual regime are considered: the period of the absence of the ice cover (IC) and the period with the existence of the IC. The vertical distributions of temperature and water salinity in Lake Uchum have been calculated in different seasons. The results of the calculations are consistent with the data of field measurements.
Keywords: salt lake, one-dimensioanal vertical model, temperature and salinity profiles of water. DOI: 10.17516/1997-1397-2019-12-1-100-108.
Introduction
Many salt lakes are meromictic, which means that no complete mixing of the water column occurs down to the bottom for at least one year. In a stratified lake, an upper (epilimnion) and a deep (hypolimnion) layer are distinguished where density gradients are small. In the intermediate layer of water (metalimnion) there is a sharp density gradient. The degree of mixing of the water masses depends on the vertical density gradient. The stability of the density stratification depends on the vertical gradients of temperature and salinity. In the absence of
* belolip@icm.krasn.ru
tsv@icm.krasn.ru
iibp@ibp.ru
§ zykovvv@yandex.ru ^ rogozin@ibp.ru © Siberian Federal University. All rights reserved
ice cover, vertical mixing of stratified lakes is obtained by vertical turbulent exchange and wind currents. To determine the dynamics of the thickness of the ice cover, a simplified model based on a quasi-stationary temperature regime in the solidified region can be applied. Due to ice formation as a result of water crystallization, salt is released in salty lakes during the winter period. Thus, unstable density stratification leads to intensive vertical circulation. A layer of convective mixing (LCM) is formed.
Mathematical models of different levels for study hydrophysical processes in stratified reservoirs are considered. It is necessary to take into account the main features of the study object when choosing mathematical models for the study of hydrophysical processes. For the seas and deep lakes there were developed models based on the equations of hydrodynamics of turbulent motion and propagation of salinity and temperature, taking into account the influence of Coriolis forces, wind stress on the surface and friction at the bottom [1,2]. To calculate the convection of deep lakes of the Baikal type, a quasi-two-dimensional model is used, which takes into account the compressibility of water and two Coriolis parameters [3]. Two-dimensional vertical plane (x, z) model with two-layer temperature stratification is used for elongated pools in the longitudinal direction [4]. Two-dimensional x, z-models are often used to model stratified flows in water bodies [5-8]. One-dimensional models are used to estimate the vertical temperature distribution in low-flow reservoirs [9-11].
The most studied by the meromictic lake of the south of Siberia is the salt lake Shira. The maximum depth of the lake is above 20 m. The mineralization of Lake Shira in the epilimnion is about 15 g/l, in the hypolimnion — 19 g/l. A special issue of the journal Aguatic Ecology (44(2010), no. 3) is devoted to the results of research on this lake.
The closed-basin salt lake named Lake Uchum is located in the South of Siberia (Krasnoyarsk Krai, Russia) at 90 km North of lake Shira. This small lake has an oval form (1.5x4 km); the maximum depth does not exceed 7.9 m (2015). The total salinity in the epilimnion is about 24 g/L, in the hypolimnion — about 34 g/L. There are data of the vertical distributions of water temperature and salinity in Lake Uchum. The thermal vertical structure is characterized by pronounced seasonal adjustment, which is due to the annual cycle of heat exchange between the atmosphere and the lake surface. Unlike Lake Shira in which water temperature in the bottom layer changes only little in the winter (less than 0.5 degrees), in Lake Uchum the temperature of the bottom layer changes by several degrees in winter. In spring, due to the melting of ice, a thin surface layer of low salinity is formed. In the process of wind mixing, the salinity in the upper layer increases. In spring and summer due to warming of the top water layer the density stratification increases. In autumn cooling of the reservoir takes place. The water temperature remains homogeneous down to a depth of 5 m before the freeze-up. Below his depth temperature inversion occurs and stratification is produced by a salinity gradient.
1. A one-dimensional model (in vertical direction) of a salt lake
Changes of temperature and salinity of water in stratified lakes mainly occur in two dimensions: vertical and temporal; in horizontal directions the changes are small. Therefore, to study the dynamics of the vertical hydrophysical structure of a stratified lake a horizontally homogeneous model can be used [11,12].
For the period of ice cover absence a one-dimensional model is based on the solution of one-dimensional (in the vertical direction) diffusion equations for temperature and salinity:
where t is time, 2 is the vertical coordinate directed downwards, H is the lake depth, Kc is the vertical turbulent exchange, C(t, z) is the temperature T (salinity S) of water. The boundary conditions are
dC dC
Kcdz lz=0 = -Fc , Kcdz \z=H = -Fch (or C\z=h = Ch),
the initial condition is
C(0, z) = Co(z).
Flows are set on the water surface and on the bottom. The heat flow on the water surface is calculated from known empirical formulas [13].
The state equation of salt water is adopted in the Boussinesq approximation for seawater [14]:
P = Po(£l + £2 TO + e3 SO), (1)
where T is the temperature, S is the water salinity p0 = 1.0254 g/cm3, e1 =0.9753, e2 = -0.00317, e3 = 0.02737, T0 = 17.5°C, S0 = 35 g/L.
A more accurate state equation of sea water [15]
P(S,T) • 103 = Pw+
+ S • (0.824493 - 4. • 10~3T + 7.6438 • 10~5T2 - 8.2467 • 10~7T3 + 5.3875 • 10~9T4)+
+ S3/2 • (-5.72466 • 10~3 + 1.0227 • 10~4T - 1.6546 • 10~6T2) + 4.8314 • 10~4S2, (2)
where pw (the density of fresh water) is determined by the formula
pw = 999.842594 + 6.793952 • 10~2T - 9.09529 • 10~3T2 +
+ 1.001685 • 10~4T3 - 1.120083 • 10~6T4 + 6.536332 • 10~9T5.
The wind shear stress on water surface is defined by the Saimons formula [16] t = 1.5-10~6w2, here w is wind speed (cm/s), t is wind stress (g/(cms2).
Parametrization of the coefficient of vertical turbulent exchange.
The coefficient of vertical turbulent exchange Kz is defined by the Prandtl-Obukhov formula [11]:
K f (0.05h)2VB + Kmin if b > 0, (3)
Z \ Kmin if B< 0, ()
f dU \ 2 g d
B = ( —— )--- —, U(z) is the horizontal velocity of water flow, h is the depth of quasi-
\dz J po dz
homogeneous layer, Kmin is the background value of the coefficient of vertical turbulent exchange.
The horizontal velocity of water flow for shallow lakes is determined by a stationary hydrostatic model of slow flows of a homogeneous fluid:
^ d2U dz
KO-TY + g^ = 0,
dz2 dx
Ko dU = - — for z = 0,
I Udz = 0.
o
Po
U = 0 for z = H
dz
A slope of the water surface is — = const. From (3), (4) formulas for determination of the
ax
coefficient of the vertical turbulent exchange are obtained:
Ko = 0.5Kmin + ^/0.25Klin + (0.05h)2r/po,
(I-H - 1). '
L0
dU t ( z \ 2 „
— =- 1.5--1, h =-H.
dz p0K0\ H Г 3
2. Vertical mixing in salt lakes in winter
In intermediate latitudes in winter, an ice cover (IC) is formed on the surface of the reservoir. The algorithm for determining the vertical distributions of water temperature and salinity consists of three stages.
Stage 1: the thickness of the IC is determined. To determine the dynamics of the thickness of the ice cover, a simplified model based on the quasi-stationary temperature regime in the solidified region is applied [11,12]. In this case, from the balance of heat flow at the water-ice interface, taking into account the latent heat of the phase transition, we obtain the equation:
d^ice = ^ Tph - Tice F _ Fc f (7 )
LPice j, лгсв ^ Fw-i Fi • f \£ice):
dt £ice
1-(1 + exp(-KV£~) _
f (7ice) = -27--exp(-«V£ice).
к £ice
The heat flow at the water-ice interface Fw-i is determined from the formula [17]:
Fw-i = ^T(7k) - Tph),
where t is time, z is the vertical coordinate (axis Oz is directed downwards), £ice is the ice thickness, £k is the depth of convection distribution, H is the reservoir depth, Xice is the coefficient of ice thermal conductivity, pice is the ice density, L is the specific heat of the ice melting, T is the water temperature, Tice is the ice temperature, Tph is the water crystallization temperature which depends on the water salinity S, Ff = Fi exp(-0.47s„), £sn is the snow layer thickness, к is the coefficient of attenuation of solar radiation in the ice cover. To determine the attenuation of solar radiation in the ice layer, the empirical dependence of the square root is used [18], a = 0.00832(——) is the coefficient of heat transfer from water to the bottomice surface [17].
V см2 c-град/ 1 J
The crystallization temperature of salt water is determined by the formula [15]: Tph = -0.0575 • S + 1.710523 • 10-3 • S3/2 - 2.154996 • 10-4 • S2.
Stage 2: the thickness of a layer of convective mixing is estimated. In salty lakes, salt is released as a result of the ice formation by the crystallization of water. Unstable density stratification is formed, leading to intense vertical circulation and the formation of a layer of a convective mixing. There is an alignment of temperature and salinity in this layer. It is assumed that convective mixing extends to a horizon at which the density of water becomes equal to the density of the underlying layer of water. Since in winter the water temperature varies little in depth, the water density mainly depends on the salinity: p = p(S,T*), T* is the characteristic water temperature under the ice cover.
A simplified formulation for the assessment of the vertical structure of a stratified reservoir during the winter period is considered. Vertical distributions of the temperature and salinity of water before the formation of IC are assumed to be three-layer, T = T*. In the upper layer of water adjacent to the free surface, the water density is remains at a constant value because of
mixing (epilimnion). In the pycnocline the density varies linearly in the depth (metalimnion). In the near-bottom layer the density varies little, p « const (hypolimnion), d1, d2, d3 are the thicknesses of the corresponding layers, (di + d2 + d3 = H is the reservoir depth). The water density is assumed to be continuous at the interfaces of the layers.
Before the ice formation, we consider the vertical distribution of salinity in the form:
{So, 0 < z < d1,
So + (SM - So) --—dl■, di < z < d,
Sbt, d < z < H.
Here S0 = const is water salinity in upper, Sbt is salinity near-bottom layer, H is the lake depth, d = d1 + d2. Let's assume that the thickness of the ice cover is known. When salt water freezes, an unstable density stratification is obtained and a layer of convective mixing (LCM) is formed. The salinity in this layer S is found by the formula:
Sk = So + (Sbt — So) —1.
d2
The thickness of LCM (Ck — Cice) is determined from the condition of equality of the salt concentration in the layer of convective mixing Sk and the salinity of water at the depth Ck [19]:
Ck = Cice + \jCtce + 2d2^ice ^ —^ + d^l — 2Cice).
So, in winter period, the thickness of the convective mixing layer depends on the autumn values of the So, Sbt, d1, d2, and on Cice, Sice.
If Ck > d1 + d2 then the convective mixing layer extends to the bottom. The ice thickness Cm at excess of which the LCM will spread to the bottom is determined by the formula:
d1 + 0.5d2
1 +
So—SiCe
Sbt —So
Thus, during the winter period, the whole water column can be mixed to the bottom if
Cice > Cm.
Stage 3. Determination of the water temperature between the ice cover and the bottom. The temperature of the water before the formation of ice is presented in the form:
'T°, 0 < z < d1,
T0(-) = ^T0° + (Tbt — T0°)-—-d1, d1 < z < d,
Tt, 2 d < z < H.
Here T° is the water temperature of the upper layer, Tbt is the temperature of the near-bottom layer. The water temperature in the mixing layer on (n +1) time step T^+1 is equal to the average temperature in the layer £"+1 < z < Cfc + 1:
Tn+1 = (Cn — £+1) • Tkn + (cn+1 — cn) [ to — (d1 — 0.5(gn+1 + cn)) (Tbt — T0o)/d2]
k cn+1 _ Cn+1 *
ck cice
Determination of water temperature in the bottom layer.
It is assumed that the thickness of the bottom layer does not change. The change in the water temperature in the bottom layer is determined from the solution of the equation:
dTbt d2Tbt . dTbt
a „ „ , Aw
dt dz2 ' w dz
, Tbt — Tk dTbt
Aw—;-z—, A,
ww z=d d — Ck dz
= 0, (5)
z=H
where Xw is the coefficient of water thermal conductivity, a =
Aw
is the water thermal diffu-
Cp Po
sivity, cp is the specific heat capacity of water, p0 is density. It is assumed that the temperature of the water in the bottom layer does not vary in depth, but changes with time. From equation (5) we obtain the equation for computation Tbt:
dTbt dt
_ l(Tk - Tbt),
satisfying the initial condition Tbt \ =0 = T°t. Here p =
Aw
Cp po (H - d)(d - ) '
The numerical solution of the obtained problem is based on an implicit scheme:
+ 0.5 At-
rpn Tbt
Tn+1 _ _
bt _ 1 + 0.5 At in+1
,.n(Tn Tn\ I .,n+lT>n+1 ^ (Tk Tbt) + 1 Tk
Aw
ln
cppo(H - d)(d - en)'
1 + 0.5At|n+1
rp \ _ rp0
±bt\t=0 _ Tbt.
Refinement of the water temperature in the layer of convective mixing. Taking into account the heat transfer at the boundaries of the convective mixing layer, we obtain the equation:
dTk dt
1
Cp P0(ek - ^ice)
Tbt - Tk
Aw ;-r--a (lk - Tph)
d - ek
The formula for recalculation of Tk follows from this equation:
At
T^n+l _
'Poiin-aj
Tn + 1 I 1 '
A Lbt__L aTn+x
Aw J r-n + a T ph
■d-in
1+
At
Aw
d-H
+a
The described algorithm for determining the vertical distributions of salinity and water temperature in the winter, there is no dependence on the type of the state equation of salt water.
Tn +
c
1
Cp Poiin-iZe)
3. Results and discussion
The parametrization of vertical turbulent exchange for a shallow lake is proposed. The modification of the one-dimensional model, which allows to estimate the temperature change of the bottom layer in winter is realized.
Measurements of vertical distributions of water temperature and salinity in the deep-water area of lake Uchum (in 2015-2016) were made. The measurements were carried out with immersed multi-channel probes Data-Sonde 4a (Hydrolab, Austin, Texas, USA) and YSI 6600 (Yellow Springs, Ohio, USA).
Calculations of the annual dynamics of the vertical thermohaline structure using the equation of state were carried out (2). The results were almost identical to the calculations using the Boussinesq approximation (1).
Numerical calculations of water temperature and salinity profiles using meteorological data for the corresponding periods were made. The results of calculations (solid lines) and data of measurements (points) are presented in Figs. 1, 2. The theoretical results in different seasons are consistent with measurement data. In spring, when the inflow of fresh water prevails over evaporation, the water level in the lake increases and there is a desalination of the lake. In summer-autumn evaporation exceeds precipitation, the water level in the lake decreases, and salinity increases. It is possible take into account these processes by the balance equations for
known changes in the depth of the spring (Delta N is greater than 0) and autumn The calculated position of the lower boundary of the convective mixing layer (5.45 m) is close to the real (5.64 m). The measured thickness of the ice cover on March 11, 2016 is 72 cm, calculated one is 72.35 cm.
Fig. 1. Vertical distribution of the water salinity in Lake Uchum: solid lines — calculations, points — measurements
Fig. 2. Vertical distribution of the water temperature in Lake Uchum: solid lines — calculations, points — measurements
The variants of climate change are considered: the average annual air temperature for the "cold " year decreased by 2 degrees compared to the meteorological data of 2015-2016 hydro-logical year and "warm" (the average annual air temperature increased by 2 degrees) with the same distributions of initial temperature and salinity. In the "warm" year, water temperatures in summer epilimnion increases by about 2 degrees, in the "cold" one is reduced by 4 degrees. The beginning-end of freeze-up shifts accordingly, unlike the real one. The period with the presence of ice cover for the "cold" year increases by 16 days, and for the "warm" year decreases by 8 days. The maximum thickness of the ice: in "cold" is 84 cm, in "warm" is 61 cm, in 2016 it is 72 cm.
The main result of this work is a one-dimensional model of the thermohaline regime of a shallow salt lake.It considers the processes of diffusion of heat and salt, heat exchange with the atmosphere, the processes of evolution of the ice cover, convective mixing under the ice cover. Thus, in the first approximation, the model takes into account all the main processes that form the seasonal and interannual variability of the lake state.
The proposed simple model makes it possible to predict the long-term dynamics of the vertical thermohaline structure of a shallow salt lake.
The work was supported by the RFBR grants (16-05-00091) and Russian Foundation for Basic Research and the government of the region of the Russian Federation, grant 18-45-243002.
References
[1] B.V.Arkhipov, P.P.Koryavov, Three-Dimensional unsteady model of water dynamics and changes in salinity in the reservoir. In: Modeling and experimental studies of hydrological processes in lakes, L., Nauka, Leningrad branch, 10-13, 1986 (in Russian).
[2] G.P.Astrakhantsev, L.A.Oganesyan, L.A.Rukhovets, Three-Dimensional mathematical model of hydrothermics a closed body of water. In: Modeling and experimental studies of hydrological processes in lakes, L., Nauka, Leningrad branch, 13-16, 1986 (in Russian).
[3] E.A.Tsvetova, The Effect of the Coriolis force on convection in a deep lake: a computational experiment, J. Appl. Mech. and Tech. Phys., 39(1981), no. 4, 127-134 (in Russian).
[4] Z.N.Dobrovol'skaya, G.P.Epihov, P.P.Koryavov, N.N.Moiseev, A Mathematical model for calculating the dynamics and quality of complex water systems, Water resources, (1981), no. 3, 32-51 (in Russian).
[5] B.V.Arkhipov, V.V.Solbakov, D.A.Shapochkin, A two-dimensional vertical model of the temperature regime of the reservoir-cooler, Water resources, 22(1995), no. 6, 653-666 (in Russian).
[6] O.F.Vasiliev, A.F.Voevodin, V.S.Nikiforovskaya, Numerical simulation of temperature-stratified flows in deep lake systems, Computational Technologies, 10(2005), no. 5, 29-38 (in Russian).
[7] D.S.Dunbar, R.W.Burling, A numerical model of stratified circulation in Indian Arm. British Columbia, J. Geophys Res., 92(1987), no. C12, 13075-13105.
[8] J.M.Lavelle, E.D.Cokelet, G.A.Cannon, A model study of density intrusion into and circulation within a deep, silled estuary: Puget Sound, J.Geographycal Research, 96(1991), no. C9, 16779-16800.
[9] G.Sh.Ignatov, I.V.Kwon, One-Dimensional model of seasonal thermocline in lakes, Vodnye resursy, (1979), no. 6, 119-126 (in Russian).
[10] A.T.Zinoviev, Yu.N.Kopylov, A.A.Kuzmin, One-dimensional vertical model of the sedimentation process in a deep reservoir, Water resources, 22(1995), no. 6, 676-683 (in Russian).
[11] S.N.Genova, V.M.Belolipetskii, D.Y.Rogozin, A.G.Degermendzhi, W.M.Mooij, A one-dimensional model of vertical stratification of Lake Shira focussed on winter conditions and ice cover, Aquat. Ecol., 44(2010), 571-584.
[12] V.M.Belolipetskii, S.N.Genova, Numerical simul, Computational technologies, 13(2008), no. 4, 34-43 (in Russian).
[13] V.M.Belolipetskii, S.N.Genova, V.B.Tugovikov, Y.I.Shokin, Numerical modelling of problems of hydro-icethermics of currents, SB RAS, Novosibirsk, 1994 (in Russian).
[14] V.P.Kochergin, Theory and methods of calculation of ocean currents, Moscow, Science, 1978 (in Russian).
[15] A.E.Gill, Atmosphere - Ocean Dynamics, Department of Applied Mathematics and Theoretical Physics University of Cambridge, Cambridge, England, Academic Press, 1982.
[16] N.N.Filatov, Dynamics of lakes, Leningrad, Gidrometeoizdat, 1983 (in Russian).
[17] S.D.Vinnikov, B.V.Proskurjakov, Hydrophysics, L., Gidrometeoizdat, 1988 (in Russian).
[18] M.S.Krass, V.G.Merzlikin, Radiation physics of snow and ice, L., Gidrometeoizdat, 1990 (in Russian).
[19] V.M.Belolipetskii, S.N.Genova, A.G.Degermendzhi, D.Yu.Rogozin, The Change of the circulation regime in the salt lake Shira (Siberia, Republic of Khakassia), Reports of the Academy of Sciences, 474(2017), no. 4, 486-489 (in Russian).
Одномерная модель для исследования сезонных изменений вертикальной структуры соленого озера Учум
Виктор М. Белолипецкий Светлана Н. Генова
Институт вычислительного моделирования СО РАН Академгородок, 50/44, Красноярск, 660036
Россия
Андрей Г. Дегерменджи Владимир В. Зыков
Институт биофизики СО РАН Академгородок, 50/50, Красноярск, 660036
Россия
Денис Ю. Рогозин
Институт биофизики СО РАН Академгородок, 50/50, Красноярск, 660036 Сибирский федеральный университет Свободный, 79, Красноярск, 660041
Россия
Многие соленые озера относятся к меромиктическим, в которых в течение как минимум одного года водная толща не перемешивается до дна. В стратифицированном озере выделяют верхний (эпилимнион) и глубинный (гиполимнион) слои, в которых градиенты плотности малы. Между ними располагается слой воды (металимнион), в пределах которого градиент плотности велик. В придонном слое накапливается сероводород и отсутствует кислород. Для определения динамики вертикальной термохалинной структуры соленого озера Учум применяются одномерные математические модели. Рассматриваются два периода в годовом режиме: период отсутствия ледяного покрова (ЛП) и период с существованием ЛП. Выполнены расчеты вертикальных распределений температуры и солености воды в озере Учум в различные сезоны. Результаты расчетов согласуются с данными натурных измерений.
Ключевые слова: соленое озеро, одномерная вертикальная модель, профили температуры и солености воды.