Numerical models for calculating hydrologic processes in river and lake-river systems

We use one-dimensional (1D) and two-dimensional (2D) longitudinal-vertical mathematical models and their 2D+1D combination as well as numerical methods to study unsteady processes in the complex open channel systems under the influence of water management measures. The analysis shows the economic feasibility and efficiency of using the developed mathematical models to study hydrological process in water bodies. The study of the physical processes in complex water body, consisting of significantly different components, based on the use of only one chosen mathematical model, is uneconomical and inefficient from the viewpoint of computational expense.


Introduction
In recent decades, in our country, as well as all over the world, there has been a general trend in changing climatic conditions, leading to an increase in the number of natural disasters associated with flooding and elevated groundwater levels as a result of catastrophic floods, as well as drying up of watercourses due to lack of precipitation. In addition to natural factors, significant role in the formation of hydrothermal and environmental regimes is plaid by anthropogenic factors, such as the construction of hydropower plants on the rivers, various hydrotechnical engineering structures on the watercourses of lake-river systems, the construction of the water utilization systems in the territories adjacent to the waterways, etc. No less complexity is observed in relation to the hydrological regime in estuaries of the Siberian rivers, where due to overregulation or inadequate rainfall along with the shortage of river runoff and the effect of wind-induced surge, we can observe not only salinization of water in the estuary ducts, but also the increase in the penetration range of salty sea water into fresh river water.
To study the hydrological processes in water bodies under the impact to different natural and anthropogenic factors, one can use the set of mathematical models of different dimensionality developed at the Institute of Hydrodynamics SB RAS (IH SB RAS). These models are based on oneand two-dimensional lateral-vertical shallow water equations (Saint-Venant equations) and their combinations. Models are designed for investigation of heat and mass transfer processes, which are taking place at water bodies under the effects of different anthropogenic and natural factors.
It is assumed that the water body may represent a water utilization system, including in the general case the objects strongly differing from each other in terms of morphometric and hydraulic characteristics (separate watercourses, ponds, lakes, and their associated systems, e.g., lake-river system with tributaries or intake structures, etc.). The possible diversity of such components in the studied water utilization system or its fragment is the main problem not only in terms of choosing relevant mathematical models appropriately describing investigated physical processes, but also mathematical simulation of the above mentioned processes according to a certain selected mathematical model. The set of computer programs used in the current research is a combination of mathematical models different in terms of setting up the problems (one-dimensional, two-dimensional longitudinalvertical, etc.) with different capabilities of taking into account the effects of the various physical or any other external factors on water body [1][2][3][4]. Flexibility in the use of a set of programs is specified by provided set of corresponding correlations (linking and conjugation conditions) between individual parts of the studied water body.

The description of mathematical models
The numerical modeling of hydrological processes in water bodies was done using two types of mathematical models different in terms of their dimensionality.
2.1. One -dimensional (1D) model of hydrologic processes in watercourses and their systems. This model is based on Saint-Venant equations averaged over channel cross-section. The system allows describing dynamics of hydrological processes in water body channels with adequate accuracy [1,2] where tis time; xis the axial channel coordinate; is the lateral inflow distributed along the channel length; are the width of free surface, the cross-section area of channel, and the rate of discharge at the depth of h, respectively; w Rare the external factors (for example, wind-induced surge (negative surge), or the barometric pressure gradient); gis the acceleration gravity.
This mathematical model, being a kind of classic in the field of hydraulics, describes a wide class of problems on unsteady water flows in open channels. In this case, any of the following factors may become the source of the unsteady water flow. This can be abundant rainfall or snowmelt, the most relevant for the North-Eastern regions of Siberia and the Far East, or wind-induced surges (storms) at lakes, reservoirs, and estuarine areas of the rivers in the Far North. Moreover, this can be uneven operation mode of hydrotechnical structures, industrial enterprises, and agricultural plants; accidental discharges of industrial and domestic wastewater, as well as partial or complete destruction of hydraulic structures (hydro-electric plants, pumping stations), etc.
The solutions to the above mentioned different problems have a common goal consisting in the definition of discharge rates ) , ( t x Q and water levels ) , ( t x Z along the length of the watercourse (reservoir) at different points in time.
The given model can help to obtain the quantitative estimate of water resources availability in a water body. This is either a direct quantification, as for example in the problem of water resources efficient distribution and if necessary their internal territorial redistribution, or additional assessment, as for example in the problem of determining damage from flooding and elevated groundwater levels at possible destruction of hydraulic works. Such universality of generalized Saint-Venant model provides the opportunity to use it in the development of mathematical models to describe a very wide range of applied (engineering) hydro-economic problems, such as calculation of water flow in channels with floodplains [5], or for example, the below-mentioned study of salt wedge task.
Task. The sea-water penetration into river-mouth sections [2]. The feature of unsteady flow calculations on estuary sections of rivers is associated not only with necessity of taking into account the river water flow dynamics, but also to reflect the direction and strength of wind, atmospheric pressure and density inhomogeneities. In this regard, to study the salt water distribution process in the estuarine sections of rivers, we can use a mathematical model (1-2), represented as: The equation of the state [6]: are determined from balance dispersion equation (the conservative additive and heat balance equation): Here tis time; xis the axial channel coordinate; are the heat and salt inflows coming together with lateral inflow is the longitudinal dispersion coefficient, which is calculated by Harleman formula [7]: The mathematical model (1D) used in the research, is a part of "Complex hydraulic models" designed at the IH SB RAS to address a wide range of open channel hydraulics problems. After thorough testing and successful adjustment to the real conditions, the latter was included to the State Fund of Algorithms and Programs (USSR State Fund of Algorithms and Programs) [8].

The two-dimensional (longitudinal-vertical) model (2D) of temperature-stratified flows in deep elongated reservoirs.
The given model is based on two-dimensional equations obtained from threedimensional equations of hydrodynamics, and averaged over channel width or watercourse, as well as assumptions about the hydrostatic law of pressure [3,4]: where Here tis time; z x,are horizontal and vertical coordinates; are horizontal and vertical velocity components; is the characteristic density value; t is the coefficient equal to the sum of molecular and turbulent viscosities; is the friction resistance (shear stress) on lateral surfaces; is the friction coefficient.
Heat and salinity transfer equations were obtained by averaging the relevant three-dimensional equations: According to this model, the following functions are unknown: the time-distance distribution of water discharge rates and levels (continuity and dynamic equations are described in variables of level and discharge rate), vertical and horizontal velocity fields, as well as temperature and salinity distribution fields (in terms of space and time) and dependent density fields. The mathematical model was developed for the study of hydro-thermal and mass transfer processes in low flow densitystratified narrow deep reservoirs.

Comprehensive 2D-1D mathematical model for calculation of hydrological processes in open channel systems.
This mathematical model is developed on the basis of models being produced above, and represents the combination of one-dimensional and two-dimensional (longitudinal-vertical) models [2,3,4,10].
This mathematical model can be used to study the dynamics of wave processes occurring in the water bodies complicated in structure, having components greatly differing from each other in terms of their geometric, hydraulic, and morphometric characteristics. When building one-dimensional and two-dimensional models, the real morphometric and hydraulic characteristics of the channels and adjacent floodplains were taken into account, as well as their interaction, and, if necessary, the impact of meteorological factors (wind, atmospheric pressure) on the wave processes.
Numerical modeling was carried out step by step according to the technique developed for the above mathematical models. It is described in detail in [10]. So, when modeling, the original physical object is represented schematically by relevant (single, multiple) physical model (topological diagram), i.e. the analogue of a real object. Then the selected physical model is associated with a corresponding mathematical model (task). Compliance of both models can be achieved using sets of boundary and matching conditions of different type provided by the computer programs [1][2][3][4].
Consider the example of some applied hydrodynamic problem, which was addressed using set of computer programs, created on the basis of the aforementioned mathematical models, numerical methods and algorithms, designed at the IH SB RAS [8].

Example. The calculation of the unsteady flow in flow-through water system of Lama Lake -Lama River [4].
The main objective of the calculations was to determine the possibility of using the developed comprehensive (2D+1D) mathematical model to study wave processes in the system of "river-lakeriver" type. To this end, we considered a fragment of a lake-river system in the Putorana Lake Province located in the North-West of the Central Siberian plateau in Krasnoyarsk Territory of the Russian Federation [11]. The considered system fragment includes Lake Lama, 80 km long, extending to East-West, and rivers flowing in and out of the lake (Figure 1). Lake Lama represents a narrow deep reservoir 3.7 km wide in its middle part with the depth of 254 m, and 10-13 km wide at the west edge (with the depth of 2.5-9 m ).
The river flowing into the lake played a role of the cumulative inflow with elevation of water level ) (t Z Z pritok  ; the river flowing out of the lake was considered as a channel with arbitrary crosssection (20 km long with maximum width up to 500 m and depth of 4-7 m, Figure 1). For numerical modeling, we have built two different topological diagrams of lake-river system fragment. According to the first diagram, the lake-river system (both lake and river) were considered as two interrelated two-dimensional domains extended to ( 20 2  L km). The calculations according to the second scheme were produced on the basis of combined (2D+1D) mathematical model. In both cases the lake and the river were described by corresponding hydraulic and morphometric characteristics. In this case, the same roughness coefficient equal to 0.030 was used in both calculating domains. Numerical calculations were based on the average annual flood, which reflects real flood waves and may well agree the possible flood waves in the territory under consideration.
The boundary conditions, when carried out calculations according to both schemes were as follows: • At the left boundary (entering , corresponding to a rise of the water of the river flowing into the lake, based on the descriptions of repeatedly observed of the water raisings and recession of flood in the spring and summer periods (from late May to early August) during the 70 days. According to the observations, water rise above the low water level in the lake was at the maximum of 4.5 m (in mid-June), which then was gradually decreased (to early or mid-August) down to the original level. Besides, at this boundary (in both cases), vertical distribution of longitudinal velocity component (elliptic law) was set [12].
• At the right boundary (the closing Lama River cross section • At the joint of calculation domains (cross section of the river flowing out the lake ) the balance of discharges and equality of water levels were assigned.
The initial conditions were taken from numerical calculations aimed at determination of steady condition, where the rates of water discharges and levels corresponded to hydraulic regime ( As noted above, the numerical simulations were aimed at comparing the relevance of both mathematical models in terms of describing the flood wave propagation process in the real lake-river system, since surveillance data is not available. In passing, we investigated the problem of choosing the method of determining the vertical distribution of the horizontal velocity ) , , ( ) on the left boundary of the two-dimensional domain. The study showed that setting it according to elliptic law does not violate the general flow pattern in the lake in both stationary and non-stationary conditions. The velocity set at the border, even at the first computational step ( x x   ) starting from the border, adjusts to the velocity field of the lake in both flow regimes ( Figure 2). Possibility to assign a constant value for the vertical longitudinal velocity component, when solving the two-dimensional problem, was also noted in [13]. The calculations of the unsteady flow propagation in the system according to both models were performed with a constant time step equal to 30 s, while the total duration of the calculations of the physical process of passing flood in the lakeriver system was 70 days.
In the absence of observational data, the analysis was made to compare the results of calculations performed according to both mathematical models. The analysis has shown that: 1) the water levels in the lake in both steady and unsteady motion practically coincide (max divergence was equal to 0.14%); for the river, the divergence reaches a maximum of 2.03-2.5%; 2) when calculating water discharge values, divergence reaches the maximum of 9.0-9.5% at the closing cross section of river-lake system; 3) the computing cost for the calculation of 70-days flood according to two-dimensional (2D) model is almost four times higher than that of the same flood calculated based on integrated (2D+1D) model, other conditions being equal.

Conclusion
Experience of solving similar tasks using the set of computer programs and techniques developed at the IH SB RAS, has shown that [4, 10, 14-16]: a) the choice of the appropriate mathematical model to investigate hydrological processes occurring in a water body, as a rule, is dictated not only by the goals of the conducted study and structure of the research object, but also the opportunity to fulfill the most economical (in terms of computational costs) and exact numerical calculations without loss of relevancy of models to real object. b) the developed one-dimensional and two-dimensional mathematical models, as well as numerical techniques are efficient and therefore useful for solving a broad range of problems of applied hydromechanics. Particularly, their advantage over three-dimensional models is evident in cases where the water bodies, which represent systems of quite extended real channels and reservoirs, need consideration of the specifics of the morphometry, operation hydraulic structures, local resistances, lateral inflow, etc.