Numerical simulation of turbulent mixing and transport of biochemical substances in inland waters

Numerical modeling of inland water objects (lakes and reservoirs) and, in particular, mixing processes and transport of biochemical substances in these basins is considered. A three-dimensional hydrostatic model and a one-dimensional model LAKE based on averaging of three-dimensional equations over a horizontal section of a reservoir are used as tools for the calculations. In the LAKE model, seiche oscillations are taken into account by parameterizing the pressure gradient and horizontal eddy viscosity. A k-ε closure is used to describe the vertical exchange processes in both models. A parameterization of the Prandtl number is implemented in the 3D model which allows turbulence to persist at any values of the gradient Richardson number. The 3D model is also supplemented by equations for calculating biochemical substances by analogy with the one-dimensional biochemistry equations used in the LAKE model, which, in particular, allow us to study the transport of substances such as oxygen (O2) and methane (CH4). Numerical experiments aimed at studying turbulent mixing and transport of substances are carried out.


Introduction
This research is focused on the numerical simulation of inland water objects and, in particular, on the simulation of turbulent mixing and transport of biochemical substances in such objects. On the one hand, the study and modeling of lakes and water reservoirs is an important task, which is of interest for a wide range of problems of limnology and hydrology and, on the other hand, various approaches for the description of such objects are used in climate models [1][2][3][4][5]. The properties of inland waters play a significant role in the processes of general circulation of the atmosphere and hydrosphere, and also in climate change [6][7][8][9]. The influence of water bodies on climate is largely determined by biochemical substances, among which dissolved gases, in particular, greenhouse gases such as methane and carbon dioxide, play a special role [10]. These gases are produced during the decomposition of dead organic residues: carbon dioxide under aerobic conditions, methane under anaerobic conditions (with a very low oxygen level), and contribute to the creation of the greenhouse effect. With an increase in the concentration of greenhouse gases the chemical processes in the atmosphere change, which can lead to ENVIROMIS 2020 IOP Conf. Series: Earth and Environmental Science 611 (2020) 012013 IOP Publishing doi: 10.1088/1755-1315/611/1/012013 2 a deterioration of the environmental situation. In addition, the biochemical processes taking place in water bodies also affect water quality. Blue-green algae, or cyanobacteria, are a source of carbon dioxide, and their bloom significantly worsens the ecological conditions of water bodies.
To correctly take into account the two-way interaction of the hydrosphere and atmosphere, it is necessary to include the calculation of the thermohydrodynamic and biological characteristics of lakes and reservoirs in the climate models. Correct reproduction of the thermohydrodynamics of these objects is critical in large and mesoscale models, and the most important aspect in the modeling of inland waters is correct description of the turbulent mixing processes.

Model descriptions
There are mathematical models of various complexity that allow one to calculate the distribution of thermohydrodynamic and biochemical parameters of inland water objects. The most detailed description is given by complete three-dimensional models (see [11][12]) whose basis, in most cases, is the Reynolds-averaged Navier-Stokes system of equations under hydrostatics and Boussinesq approximations. In the present work, we used a three-dimensional hydrostatic model developed at the Research Computing Center of Moscow State University on the basis of a unified hydrodynamic code combining both RANS and DNS-(Direct Numerical Simulation), LES-(Large-Eddy Simulation) approaches for calculating geophysical turbulent flows at high spatial and time resolution (see [13][14][15]). The numerical model includes the equations of hydrodynamics in a stratified turbulent rotating fluid layer in the shallow water approximation, as well as an equation for heat transfer taking into account horizontal and vertical diffusion: Here is the velocity vector, is the deviation of the free surface from the equilibrium state, is the temperature, is the density; ( ) and ( ) are the coefficients of vertical (horizontal) turbulent viscosity and thermal diffusivity, respectively; ν, χ'are the coefficients of molecular viscosity and thermal diffusivity, f is the Coriolis parameter (assumed constant), g is the acceleration of gravity, z is the vertical coordinate (from the bottom of the reservoir z = -H (x, y) to the surface), and t is time.
is the advection operator: and and are the horizontal and vertical diffusion operators with coefficients λ and K, respectively. As for numerical methods, semi-implicit in time discretization of the equations is used, while the advection and horizontal diffusion terms are calculated explicitly.

,
To study the hydrological and thermodynamic parameters of inland water bodies on seasonal, annual, and climatic scales, one-dimensional models, which are distinguished by computational simplicity, are of greatest interest today. In the present work, the one-dimensional model LAKE [16] was considered, based on a system of equations describing the vertical distribution of momentum and heat, which can be obtained by averaging the above three-dimensional equations (1-6) over the horizontal section of the reservoir [17]: Here A(z) is the area of the horizontal section of the reservoir, p is the hydrostatic pressure, and the horizontal line implies averaging over A(z). In accordance with the above, the heat flux at the lower boundary (at the bottom) is set equal to 0, and the momentum flux is considered constant at the boundary of each horizontal section (values indicated by the bot index). The LAKE model takes into account internal wave oscillations (seiches) due to the parameterization of the pressure gradient and horizontal viscosity. Seiches are not taken into account in most existing one-dimensional (vertical) models. However, when modeling lakes and reservoirs with horizontal dimensions much less than the Rossby internal radius of deformation L R , gravitational oscillations play a significant role compared to rotation (Coriolis force) in the dynamics of the mixed layer, and the one-dimensional models lacking seiche parameterization will overestimate the mixed layer thickness, for example, during the summer months.
When writing systems of equations (1-6) and (7)(8)(9), the validity of the gradient approximation for the description of turbulent flows is also assumed. For the problems considered in this paper, it seems important that vertical mixing in three-dimensional and one-dimensional models is represented in a similar way. Therefore, to calculate the coefficients and in both models, the k-ε closure (see, for example [18]) is used. It is based on prognostic equations for the kinetic energy of turbulence (k) and its dissipation rate.
Here the term P corresponds to the generation of turbulence energy due to the velocity shear, and B describes the generation or consumption of energy due to the action of buoyancy forces, are the turbulent Schmidt numbers for the turbulent kinetic energy and the dissipation rate, respectively, and , , are empirical constants, and are the stability functions for momentum and scalars. is the turbulent Prandtl number connecting the coefficients of turbulent viscosity and IOP Conf. Series: Earth and Environmental Science 611 (2020) 012013 IOP Publishing doi:10.1088/1755-1315/611/1/012013 4 eddy diffusivity, and it is assumed constant in the standard k-ε closure. However, there are different approaches to its parameterization, which will be discussed in the next subsection.

Parameterization of turbulent Prandtl number
There is a sufficient number of papers (see, for example, [19]) devoted to the need to clarify semiempirical closures for describing shear stratified turbulent flows. For the problems of describing the processes of transfer of heat, momentum, and concentrations of impurities in inland water bodies, it is important to consider the existence of turbulent mixing under strong stable stratification and large values of the gradient Richardson number ( where N is the Brunt-Väisälä frequency, and is the velocity shear). Numerous empirical dependences (for example, [18][19][20][21][22]) have been proposed, characterizing mainly the asymptotic behavior of the dependence , but not allowing one to establish with sufficient accuracy the dependence at intermediate values of Ri characterizing the transition from weak to strong turbulent stratified shear flows of a fluid. In [23][24], a theory of turbulent closure based on the balance equations for the kinetic, potential energies of turbulence, turbulent momentum fluxes, potential temperature and relaxation equation for a turbulent time scale is proposed. This model also allows one to remove the restrictions on the existence of turbulence at large Richardson numbers. In the framework of this study, we considered the parameterization obtained from the model proposed in [25]. The procedure for obtaining equations for the average velocity , density , as well as the kinetic energy of turbulence k and the variance of density fluctuations is similar to that used in the kinetic theory of gases; in this case the equation for the one-point distribution function is the so-called kinetic equation, and the Reynolds stresses are calculated from the known distribution function. The proposed approach made it possible to take into account some important but usually neglected effects, for example, the dependence of the vertical anisotropy of turbulence on stratification, a non-gradient correction to the classical expression for turbulent mass flux, the reverse transition from the potential energy of turbulent fluctuations to kinetic one, which leads to the elimination of restrictions on the existence of turbulence at large Richardson numbers.
The obtained dependence of the Prandtl number was implemented in a 3D hydrostatic model in order to take into account the peculiarities of turbulent mixing in stratified water bodies. The calculation results are presented in [26]. Taking into account the parameterization leads to smoothing of all sharp changes in the vertical distributions of turbulent kinetic energy, temperature, and the shock layer thickness. A significant decrease in the maximum of the temperature gradient in the thermocline was observed in [26]. The results obtained are determined by the possibility of describing turbulent mixing at values of Ri≫1 within the framework of the implemented approach.
The use of parameterization can influence the transport of biochemical substances in small inland water bodies, in particular, through the thermocline.

Numerical experiments. The role of horizontal dimensions of the reservoir on mixing
A numerical implementation of the classical Kato-Phillips laboratory experiment [27] which serves as a basic material for the calibration of turbulent closure has been carried out. The obtained results were compared with calculations of the changes in the mixed layer depth in idealized reservoirs with different horizontal dimensions. In the Kato-Phillips experiment, a homogeneous stratified fluid of sufficiently large depth is considered, and vertical boundaries are absent. The initial temperature profile is linear and the only source of turbulence is a wind with a constant velocity. In a classical statement described in [27], a ring tank was considered whose surface was influenced by friction stress in the direction of the circle. The inner and outer diameters were 106.7 and 152.4 cm, respectively, so that the channel width was 22.8 cm. The depth of the tank was 28 m.
The results of this experiment are well described by the theoretical formula for the change in the mixed layer thickness with time [28]:  (14) where is the mixed layer thickness, and is the friction velocity on the surface ( ). Also, idealized rectangular reservoirs with various horizontal sizes (10 and 1000 meters) were considered. It was demonstrated that the models similarly describe the dynamics of the mixed layer thickness, and a good agreement with the theoretical formula by Price [28] for the Kato-Phillips experiment was obtained (see Figure 1). It should be noted that with increasing size of the reservoir, the mixing speed also increases, and the larger the reservoir, the closer the result will be to the result of the classical Kato-Phillips experiment. The next step was to compare the various characteristics linked with the mixing processes in lakes, and we analyzed the vertical temperature distribution and the distribution of the U-velocity on the 7th day of the calculations (see Figure 2a   Agreement in temperature between the 1d and 3d models is extremely good for a scale of several days, and this is precisely the case where it is important, for example, for the tasks of forecasting. As for the velocity profiles, there are certain differences in the results, which can be related to the features of the models: the one-dimensional model is based on the averaging of the three-dimensional equation, and the parameterization of seiches is still a simplification of the influence of the pressure gradient on the flow dynamics. Note that the decrease in the velocity fluctuations below the thermocline is influenced by molecular viscosity.

Biochemical substances in inland waters and simulation of their transport
In this study, a code section for calculating the biochemical characteristics of the lake was developed and implemented within the framework of a three-dimensional hydrostatic model of the reservoir. The model was supplemented with equations for calculating biochemical substances by analogy with the one-dimensional equations of biochemistry used in the one-dimensional LAKE model [16]. The equations describe the transport, diffusion, and reactions for such substances as methane (CH4), oxygen (O2), carbon dioxide (CO2), nitrogen (N), argon (Ar), living and dead particles of phyto-and zooplankton, and so on, and are as follows: where is the concentration of substances, and are the coefficients of turbulent and molecular diffusion, respectively, and the term describes the reactions. The main reactions that are considered primarily at this stage of study are photosynthesis and methane oxidation: Using the developed biochemistry code, a numerical experiment was carried out, similar to the above-described numerical formulation of the laboratory experiment of Kato-Phillips. Boundary conditions have been added for substances such as oxygen and methane. The methane flow at the bottom is set constant: the oxygen concentration at the bottom is assumed to be zero: . The fluxes into the atmosphere are specified through the gas exchange coefficient calculated according to the so-called "surface renewal model" [29][30], which takes into account the dissipation rate of the kinetic energy of turbulence, which is explicitly calculated in the k-ε-closure used in both the 3D and 1D models of the reservoir.  It should be noted that the most important characteristic in the study of methane is the flux of this gas into the atmosphere, because an increase in the methane concentration in the atmosphere enhances the greenhouse effect [9], since methane intensively absorbs the thermal radiation of the Earth. A time series of values of methane flux into the atmosphere obtained using the three-dimensional model is shown in Figure 4.

Conclusions
We studied the processes of turbulent mixing and transport of biochemical substances using numerical models. The three-dimensional hydrostatic model and the one-dimensional LAKE model, which are based on averaging of three-dimensional equations over a horizontal section of a reservoir, were chosen for carrying out the calculations. A k-ε closure was used to describe the vertical exchange processes in both models. In the LAKE model, seiche oscillations are taken into account by parameterizing the pressure gradient and taking into account horizontal eddy viscosity. A parameterization of the turbulent Prandtl number obtained on the basis of a closure that takes into account a two-sided transformation of the kinetic and potential energies of turbulent pulsations was implemented in the 3D model. Also, the 3D model was supplemented by equations for calculating biochemical substances by analogy with the one-dimensional biochemistry equations used in the LAKE model. The equations describe transport, diffusion, and reactions for substances such as: methane (CH4), oxygen (O2), carbon dioxide (CO2), nitrogen (N), argon (Ar), living and dead particles of phyto-and zooplankton. The authors carried out a number of numerical experiments based on the laboratory Kato-Phillips experiment. The first series of numerical experiments demonstrated that for lakes with horizontal sizes much less than the internal Rossby radius seiches should be taken into account. The agreement between 1D and 3D models was the best for a scale of several days, and this is precisely the case where it is important, for example, for forecasting. The second numerical experiment was carried out to study the dynamics of two dissolved gases: methane (CH4) and oxygen (O2).