A method for inversion and prediction of reservoir sediment and temperature dual density flow

Abstract

The present invention belongs to the technical field of hydraulic engineering and discloses a method for inverting and predicting reservoir sediment and temperature dual density flows. First, a water-temperature-sediment-ice coupling model is established; the water-temperature-sediment-ice coupling model includes a hydrodynamic module, a sediment module, a water-temperature module, and an ice module. A cross-sectional model of the reservoir along the flow and water depth is then constructed, divided into several grid cells, and boundary and initial conditions are set. The water-temperature-sediment-ice coupling model is then solved to obtain the spatiotemporal distribution of sediment content, water temperature, and ice conditions in the reservoir area. Taking into account the hydrodynamic-water-temperature-ice-sediment coupling effect, the present invention accurately captures the sediment movement process in different types of sediment-laden river reservoirs in cold regions, as well as the thermal evolution process and patterns of reservoir water temperature and ice conditions.

Description

Technical Field

[0001] The invention belongs to the technical field of water conservancy projects, relates to the inversion and prediction of reservoir water temperature, sediment and ice conditions, and in particular to an inversion and prediction method for reservoir sediment and temperature dual density flow. Background Art

[0002] When sediment-laden water enters a reservoir, it experiences a significant density difference with the water within it, resulting in sediment hyperpycnal flow. This hyperpycnal flow not only alters the reservoir's sedimentation pattern but also affects the water temperature stratification within the reservoir, resulting in a temperature distribution that differs from that of reservoirs on rivers with lower sediment concentrations. These changes directly impact the thermal processes within the reservoir and its upstream and downstream channels. Rivers and reservoirs north of 30°N (referred to as cold-region reservoirs) experience varying degrees of ice conditions annually. During the ice-flood season, ice jams and ice dams are prone to form in the river channel, reducing the cross-sectional area of the river and causing water level rises, which can trigger ice-flood floods. The downward shift in the location of the maximum flow velocity during the ice-flood season causes intense riverbed scouring, leading to higher sediment concentrations during the ice-flood season. Reservoir areas typically experience vertical water temperature inversions in winter. Due to their higher density, the inflowing, low-temperature water, laden with sediment, can sink to the bottom of the higher-temperature water, altering the existing buoyancy pattern and indirectly affecting the spatial and temporal distribution of ice conditions within the reservoir.

[0003] It's worth noting that the current development of water conservancy and hydropower projects is concentrated in cold, high-altitude river source areas. For example, hydropower development in the upper reaches of the Yellow River has provided a large amount of clean energy to the northwest region, helping to alleviate energy shortages. However, building reservoirs on rivers with high sediment concentrations inevitably faces the coupling of sediment, ice conditions, and water temperature. For example, when the inflowing sediment concentration is high, the reservoir area may experience an inverted temperature distribution, sediment density flow may accelerate the discharge of water from the stagnant layer at the bottom of the reservoir, delay the formation and thawing of ice in the reservoir area or at the tail end, and reduce the range of ice formation.

[0004] The interaction between ice, water and sediment in a reservoir is mainly reflected in the impact of sediment density flow on the thermal state of the reservoir. Since the movement of sediment density flow itself is relatively complex, coupled with the influence of water temperature stratification, the movement of the dual density flow coupled with water temperature and sediment becomes even more complex. When high-sediment-content water enters the reservoir, it not only destroys the existing water temperature stratification structure, but also has a significant impact on the stratified flow of the reservoir and the material exchange within the water body. In this field, there are few existing methods for inverting the water temperature and ice conditions of reservoirs on high-sediment-content rivers. In particular, the sediment settling velocity is mostly calculated based on the median particle size of a single group. When the water flow has a high sediment concentration and the sediment particle size range is widely distributed, the error is large, affecting the water temperature inversion effect.

[0005] Water temperature is a key factor influencing fish reproduction and spawning. Different fish species have specific water temperature requirements during their reproduction. A suitable water temperature stimulates the maturation of fish gonads and encourages them to enter a reproductive state. Furthermore, sediment carries a rich supply of biogenic substances. The movement and distribution of sediment alters the physical environment and material fluxes within the water, in turn impacting phytoplankton, plants, and fish. Therefore, accurately predicting water temperature changes and sediment movement can provide a scientific basis for the conservation and management of aquatic ecosystems in rivers and reservoirs.

[0006] In summary, the presence of high sediment concentrations and ice conditions poses significant challenges to the inversion of reservoir water temperature and ice conditions, placing higher demands on the applicability of models. Furthermore, accurate prediction of water temperature and sediment in reservoirs along high-sediment-concentration rivers is crucial. Therefore, research on coupled water-ice-sediment models to invert water temperature and ice conditions in reservoirs along sediment-rich rivers in cold regions has significant scientific and engineering value. This research can provide a scientific basis for reservoir temperature and sediment control, ice disaster prevention and mitigation, and ecological protection along sediment-rich rivers in cold regions. Summary of the Invention

[0007] The purpose of the present invention is to overcome the deficiencies in the prior art and provide a method for inverting and predicting dual density flows of reservoir sediment and temperature, to accurately and efficiently simulate the water temperature, ice conditions and sediment content of the reservoir, and to invert, predict and evaluate the thermal spatiotemporal changes of the water temperature and ice conditions of high-sediment-content river reservoirs and the degree and scope of the influence of sediment on the water temperature and ice conditions.

[0008] The key link in this invention is to accurately simulate the sediment density flow by calculating the sediment grouping velocity. When the sediment density flow is moving in the water, it is subject to the resistance of the upper layer of clean water. Therefore, the density flow will also drive a part of the upper layer of clean water forward during the movement. At the same time, when the density flow is mixed with the clean water, it will also occupy the position of the original clean water, thus forming a clean water circulation above the density flow. Figure 2 As shown, the surface water near the point where the density current enters moves in the opposite direction to the muddy water, also known as reverse compensation flow. High sediment content in incoming flow causes it to submerge, forming a bottom-layer density current. The upper clear water near the reservoir tail forms vortices. However, the suction effect of the water intake in the area in front of the dam creates strong turbulence and a complex flow pattern. While the direction of the density current near the dam is the same as that of the water at the upper outlet, there is a significant velocity difference, resulting in an unbalanced force on the middle water layer, significantly disturbing the reservoir flow field and further altering the reservoir's thermal state.

[0009] The present invention provides an inversion and prediction method for reservoir sediment and temperature dual density flow, which comprises the following steps:

[0010] Step 1: establishing a water temperature-sediment-ice coupling model; the water temperature-sediment-ice coupling model includes a hydrodynamic module, a sand module, a water temperature module, and an ice module;

[0011] The hydrodynamic module includes the water flow continuity equation, momentum equation, turbulence equation, water state equation about water density, and free water surface equation obtained by vertical integration of the continuity equation;

[0012] The sand module includes suspended matter convection and diffusion equations after sediment settling;

[0013] The water temperature module includes a heat transfer equation;

[0014] The ice module includes an initial ice sheet formation equation and an ice sheet thermal generation and disappearance equation;

[0015] Step 2: Construct a cross-sectional model of the reservoir along the flow and depth directions, divide it into several grid cells, and set boundary conditions and initial conditions;

[0016] Step 3 solves the water temperature-sediment-ice coupling model to obtain the spatiotemporal distribution of sediment concentration, water temperature, and ice conditions in the reservoir area. This step includes the following sub-steps:

[0017] Step 3.1 Based on the spatial distribution of hydraulic pressure, sediment content, and ice conditions at the previous moment, the hydraulic module is used to preliminarily obtain the hydraulic spatial distribution at the current moment;

[0018] Step 3.2: Based on the spatial distribution of sediment content at the previous moment and the hydraulic spatial distribution at the current moment, the spatial distribution of sediment content at the current moment is obtained through the sediment content module;

[0019] Step 3.3: Based on the water temperature spatial distribution at the previous moment and the hydraulic spatial distribution at the current moment, obtain the water temperature spatial distribution at the current moment through the water temperature module;

[0020] Step 3.4 traverses the surface grid cells of the reservoir to determine whether the water temperature of the surface grid cells has dropped below freezing point at the current moment. If so, the initial ice sheet thickness is determined using the initial ice sheet formation equation, and the surface water temperature is set to freezing point and the properties of the surface grid cells are changed from ice-free to ice-covered. Then, proceed to step 3.5. Otherwise, proceed directly to step 3.5.

[0021] Step 3.5 traverses the surface grid cells of the reservoir and calculates the ice thickness change value of the surface grid cells with existing ice thickness to determine the growth and decline of the ice cover. If the ice thickness change value is positive, the ice thickness is growing; if the ice thickness change value is negative, the ice thickness is melting. When the ice thickness at the current moment is less than the melting ice thickness change value, the current surface grid cell state changes from ice-covered state to ice-free state.

[0022] Increase the time step and repeat the above steps S3.1-3.5 until the upper limit of the set time period is reached, and the reservoir hydraulic temporal and spatial distribution, sediment content temporal and spatial distribution, water temperature temporal and spatial distribution, and ice condition temporal and spatial distribution of the reservoir area are obtained.

[0023] The core idea of this invention is to further carry out the calculation of ice conditions based on the calculation of the impact of high-concentration sediment-laden water flow on reservoir hydrodynamics and water temperature. This invention proposes a water temperature-sediment-ice coupling model, such as Figure 1 As shown in Figure 2, it includes a hydrodynamic module, a sediment module, a water temperature module, and an ice module. The sediment module uses the group median particle size to determine the sediment settling velocity of the group, simulating the coupling effect of temperature-density flow and sediment-density flow, thereby inverting and predicting the spatiotemporal variation of water, ice, and sediment in the reservoir.

[0024] The key link in this invention is to accurately simulate the sediment density flow by calculating the sediment grouping velocity. When the sediment density flow is moving in the water, it is subject to the resistance of the upper layer of clean water. Therefore, the density flow will also drive a part of the upper layer of clean water forward during the movement. At the same time, when the density flow is mixed with the clean water, it will also occupy the position of the original clean water, thus forming a clean water circulation above the density flow. Figure 2 As shown, the surface water near the point where the density current enters moves in the opposite direction to the muddy water, also known as reverse compensation flow. High sediment content in incoming flow causes it to submerge, forming a bottom-layer density current. The upper clear water near the reservoir tail forms vortices. However, the suction effect of the water intake in the area in front of the dam creates strong turbulence and a complex flow pattern. While the direction of the density current near the dam is the same as that of the water at the upper outlet, there is a significant velocity difference, resulting in an unbalanced force on the middle water layer, significantly disturbing the reservoir flow field and further altering the reservoir's thermal state.

[0025] In step 1 above, for the hydrodynamic module, the water flow continuity equation is:

[0026]

[0027] Where: U is the longitudinal velocity, m / s; x is the flow direction; W is the vertical velocity, m / s; z is the vertical direction of the reservoir; q is the net inflow per unit width of the lateral grid unit, m 3 / s; B is the lateral width of the reservoir at different elevations along the flow direction, m.

[0028] The momentum equation includes the longitudinal momentum equation and the vertical momentum equation;

[0029] The longitudinal momentum equation is:

[0030]

[0031] Vertical momentum equation:

[0032]

[0033] Where: α is the riverbed slope, rad; τ xx is the turbulent normal stress, N / m 2 ; τ xz is the turbulent shear stress, N / m 2 ; η is the water surface elevation, m; U x is the x-direction component of the tributary velocity in the water flow direction, m / s; g is the acceleration due to gravity, m / s 2 ; ρ is the water density, kg / m 3 ; p is the hydrostatic pressure, N / m 2 .

[0034] Turbulence equation:

[0035]

[0036] Where: k is the turbulent kinetic energy, m 2 / s 2 ;ε is the turbulent kinetic energy dissipation rate; ν t is the turbulent eddy viscosity coefficient, ν t =C μ k 2 / ε,m 2 / s; ρ is the water density, kg / m 3 ; P is the turbulent kinetic energy generation term, m 2 / s 3 ; G is the buoyancy term, m 2 / s 3 ;P k 、P ε is the turbulent kinetic energy and dissipation rate generated by boundary friction; σ k , σ ε are the turbulent kinetic energy and the dissipation rate Prandtl number, which are 1.0 and 1.3 respectively; δ t is the turbulent Prandtl number, which is 1.0; C μ 、C ε1 and C ε2 are empirical constants, which are 0.09, 1.44 and 1.92 respectively.

[0037] The present invention takes into account the change in density in the gravity term in the control equation, and the state equation of water temperature, total dissolved solids and suspended solids, i.e., the water state equation, which is in the following form:

[0038]

[0039] Where: T w is the water temperature; ρ T is the density of water under the influence of temperature; △ρ ss,iis the water density increment caused by the i-th level of suspended sediment; C i is the sediment content per unit volume of water; γ s,i is the dry bulk density of suspended sediment at the i-th level; n represents the number of sediment level groups.

[0040] In order to simulate the change of water level, the free water surface equation is also needed. By integrating the continuity equation from the bottom to the surface, the free water surface equation is obtained:

[0041]

[0042] Where: B η is the lateral width of the reservoir surface along the flow direction, m; η is the water surface elevation, m; h is the riverbed elevation, m.

[0043] For the sand module, the present invention uses the suspended sediment convection diffusion equation to represent the sediment transport process in the reservoir, while taking into account the influence of sediment sedimentation. The suspended sediment is divided into i levels, and the suspended sediment convection diffusion equation for each level is:

[0044]

[0045] Where: C i is the sediment content per unit volume of water; ω s,i is the sediment settling velocity; D x 、D z are the longitudinal turbulent diffusion coefficient and the vertical turbulent diffusion coefficient respectively;

[0046] The formula for the sediment settling velocity of a single particle uses the commonly used Stokes formula as follows:

[0047]

[0048] Where: d i is the median particle size of the suspended sediment particles at level i, m; υ is the kinematic viscosity of water, m 2 / s; γ T is the bulk density of water without sand, γ T =ρ T g, N / m 3 ; γ s,i is the dry bulk density of suspended sediment at level i, γ s,i =ρ s,i g, N / m 3 ,ρ s,i is the density of suspended sediment at level i, kg / m 3 .

[0049] Generally speaking, density currents can only carry fine sediment particles with a diameter less than 25 μm. After entering the reservoir, medium-coarse sediment with larger particle sizes settles quickly under the influence of gravity, mainly settling at the end of the reservoir. Fine sediment has better tracking ability. The finer suspended sediment can even reach the front of the dam under the action of density currents and be transported downstream. Taking sediment particles with a particle size of 25 μm as an example, in still water at 20°C, according to the Stokes formula, the settling velocity will reach 38.2 m / day (4.4×10 -4 m / s). Therefore, this paper ignores the impact of medium and coarse sand on the flow temperature field near the dam, and mainly considers the fine sand (d < 25 μm). Fine sand is further divided into three levels: 0-5 μm, 5-10 μm, and 10-25 μm. The settling velocity corresponding to the median particle size is calculated using formula (2-2).

[0050] For the water temperature module, the heat transfer equation is:

[0051]

[0052] Where: T w is the water temperature, ℃; D x is the longitudinal diffusion coefficient, m 2 / s;D z is the vertical diffusion coefficient, m 2 / s; Q is the rate of lateral heat flux in and out of the grid unit, Q = qT q , T q is the water temperature of the lateral flow, °C; S is the source-sink term, including surface heat exchange, riverbed heat exchange, ice-water heat exchange, and / or solar radiation absorption at depth z, J / m 3 / s.

[0053] Water surface heat exchange mainly consists of three parts: radiation, evaporation and heat conduction, which can be expressed by the following formula:

[0054] φ n =H s +H a -H b -H e -H c (3-2)

[0055] Where, φ n Net heat flux through the water surface, W / m 2 ;H s is the solar shortwave radiation absorbed within 1m of the surface water, W / m 2 ;H a is the longwave radiation returned by the atmosphere, W / m 2 ;H b is the long-wave radiation from the water surface, W / m 2 ;H eis the evaporation heat flux, W / m 2 ;H c is the heat conduction flux, W / m 2 .

[0056] (1) Solar shortwave radiation H s

[0057] H s =H s0 (1-γ)β(1-0.65C 2 ) (3-3)

[0058] Where H s0 is the shortwave solar radiation reaching the water surface on a sunny day, W / m 2 ; γ is the reflectivity of the water surface; β is the degree of solar radiation absorption within 1m of the water surface; C is the cloud cover rate.

[0059] Solar radiation entering a body of water is attenuated according to Beer's law:

[0060] H s (z) = (1-β)(1-γ)H s0 e -ζz (3-4)

[0061] Where H s (z) is the solar radiation at water depth z, W / m 2 ;ζ is the solar radiation attenuation coefficient, 1 / m.

[0062] (2) Atmospheric longwave radiation H a

[0063] The solar energy absorbed by the atmosphere is emitted to the ground in the form of long waves. The intensity of the long-wave radiation depends on the temperature and cloud cover and can be calculated using the Stefan-Boltzman law:

[0064] H a =(1-γ a )σε a (273.15+T a ) 4 (3-5)

[0065] Where, T a (℃) is the air temperature at 2m above the water surface, ℃; γ a is the long-wave reflectivity; σ is the Stefan-Boltzman constant, which is 5.67×10 -8 W / (m 2 ·K 4 );for ε a Atmospheric emissivity, atmospheric emissivity on a clear day ε ac It can be calculated using Idso and Jackson formulas.

[0066] ε ac =1-0.261·exp(-7.77×10 -4 T a 2 ) (3-6)

[0067] In cloudy weather, the Bolz formula can be used to correct it:

[0068] ε a =ε ac (1+K c ·C 2 ) (3-7)

[0069] The parameter K c Related to cloud height, the Tennessee Engineering Administration recommends an average value of 0.17.

[0070] (3) Water body long wave return radiation H b

[0071] The long-wave return radiation of water is an important part of the heat loss of water, which can be calculated using the Stefan-Boltzman law:

[0072] H b =σε w (273.15+T ws ) 4 (3-8)

[0073] Where, ε w is the long-wave emissivity of water; T ws is the surface water temperature of the reservoir, ℃.

[0074] (4) Evaporation heat loss H e

[0075] The evaporation process of water requires the absorption of heat, and the evaporation heat loss expression is:

[0076] H e =f(W)(e s -e a ) (3-9)

[0077] Where f(W) is the wind function:

[0078] f(W)=9.2+0.46W 2 (3-10)

[0079] Where W is the wind speed at 2m above the ground, m / s; e s is the saturated evaporation pressure of the water surface, mmHg; e a is the atmospheric vapor pressure, mmHg.

[0080] (5) Heat conduction flux H c

[0081] H c =C c f(W)(T ws -T a ) (3-11)

[0082] Where C c is the Bowen constant, which is taken as 0.47 mmHg / ℃.

[0083] Compared with surface heat exchange, riverbed heat exchange is smaller in magnitude and many models do not consider riverbed heat exchange. Research shows that riverbed heat exchange must be considered to accurately simulate the water temperature of the stagnant layer at the bottom of the reservoir. The riverbed heat exchange formula is:

[0084] H bw =K bw (T b -T wb ) (3-12)

[0085] Where: H bw is the heat exchange flux between the riverbed and the bottom water, W / m 2 ;K bw is the heat exchange coefficient, W / (m 2 ·℃), generally take 0.3; T wb is the water temperature at the bottom of the reservoir, ℃; T b is the riverbed temperature, ℃, which is generally estimated using the multi-year average temperature.

[0086] When covered by ice, the heat flux form of the reservoir surface changes from water-air heat exchange to ice-water heat exchange. The ice-water heat exchange calculation can be expressed by the following formula:

[0087] H wi =h wi (T ws -T m ) (3-13)

[0088] Where: H wi is the heat exchange flux between the reservoir surface and the bottom of the ice sheet, W / m 2 ;h wi is the ice-water heat exchange coefficient, Wm -2 ℃ -1 ;T m is the ice-water interface temperature, 0℃.

[0089] For the ice module, the formation and melting of ice cover is a complex heat exchange process, influenced by factors such as reservoir surface weather, subglacial water temperature, and ice thickness. The ice cover is also affected by the reservoir's hydrodynamic factors, which influence the ice-water heat exchange coefficient. Heat exchange during freezing and melting primarily involves ice-air heat exchange, ice-water heat exchange, and internal heat conduction within the ice.

[0090] Before the initial ice sheet is formed, the water surface temperature drops to freezing point. As the water surface continues to lose heat, the initial ice sheet begins to form. The thickness of the initial ice sheet is calculated by the following formula:

[0091]

[0092] Where: θ0 is the initial ice thickness, m; C p is the constant-pressure specific heat capacity of water (constant), J / (kg·℃); θ s is the vertical grid thickness of the surface layer, m; ρ is the water density, kg / m 3 ρ i is the density of ice, kg / m 3 ;L f is the latent heat of freezing (constant), J / kg.

[0093] The ice sheet thermal generation and dissipation equation for ice thickness θ is given by the following equation:

[0094]

[0095] Where: θ represents ice thickness, m; T a is the air temperature, ℃; k i is the thermal conductivity of ice, W / (m·℃); h ai is the ice-air heat exchange coefficient, W / (m 2 ·℃).

[0096] In step 2 above, the cold-region reservoir waters were studied and divided into a number of triangular or rectangular units (i.e., grid cells) along a specified rectangular coordinate system. The hydraulic pressure, sediment concentration, water temperature, and ice conditions at each grid cell's node were used to characterize the hydraulic pressure, sediment concentration, water temperature, and ice conditions at that location. The spatiotemporal distribution of hydraulic pressure includes the temporal variations of flow velocity and water surface elevation; the spatiotemporal distribution of sediment concentration is the temporal variation of sediment concentration in the reservoir area; the spatiotemporal distribution of water temperature is the temporal variation of water temperature in the reservoir area; and the spatiotemporal distribution of ice conditions is the temporal variation of surface ice thickness in the reservoir area.

[0097] In the above step 3, the present invention discretizes all the models constructed in the above step 1 using the finite volume method to obtain a discrete set of equations for each model on the grid unit, thereby solving the discrete set of equations.

[0098] Based on the above analysis, it can be seen that the present invention takes into account the impact of sediment particle size grouping on water density. Using the suspended matter convection and diffusion equation, different levels of sediment are calculated and the corresponding sediment content distribution is obtained, further updating the water density.

[0099] The hydraulic, sediment content, water temperature, and ice condition parameters obtained at the previous moment are assigned to each variable as the initial value of the hydraulic, sediment content, water temperature, and ice condition parameter variables at the current moment. The hydraulic, sediment content, water temperature, and ice condition of the reservoir area at each moment are solved through time iteration, thereby obtaining the spatiotemporal distribution of hydraulic, sediment content, water temperature, and ice condition of the reservoir area.

[0100] The convergence conditions of the hydrodynamic field, temperature field, and sediment field are:

[0101] (1) Under the condition of accurately given initial fields, the various field data obtained by calculation are the time series results to be obtained;

[0102] (2) When the initial field cannot be accurately given, the initial conditions can be determined by random initialization, and then multiple time cycles are used. When the result of the mth cycle and the result of the m-1th cycle are within the allowable difference range, the result of the mth cycle is considered to be the intended time series result.

[0103] The reservoir boundary conditions include inflow flow, inflow water temperature, inflow sediment volume and gradation, meteorological conditions, and outflow flow. The initial conditions are a given flow field, temperature field, and sediment field. The initial ice thickness is zero, meaning there is no ice cover on the reservoir surface. The water temperature at the inlet boundary is the actual or given water temperature at the reservoir tail, and the velocity is assumed to be uniform. Fully developed turbulence is assumed at the dam body and outlet section, and a no-slip boundary condition is used on the dam body surface, which is also an adiabatic boundary. A no-slip boundary condition is used at the reservoir bottom.

[0104] Stability is an important aspect to consider when solving numerical simulations, and can be verified by the time step. In numerical calculations, the time step of the grid unit must meet the following formula conditions:

[0105]

[0106] Where △t is the time step, s; △x and △z are the lengths of the grid cell along the x and z directions respectively; A x is the longitudinal eddy viscosity, m 2 / s;A z is the vertical eddy viscosity, m 2 / s;Q e is the total flow into or out of the grid cell, m 3 / s; V is the volume of the grid unit, m 3; H is the depth of the water column where the unit grid is located, m; ρ is the water density; △ρ is the maximum density difference of the water column where the unit grid is located, kg / m 3 .

[0107] The reservoir sediment and temperature dual density flow inversion and prediction method provided by the present invention has the following beneficial effects:

[0108] (1) The present invention establishes a water temperature-sediment-ice coupling model from the perspective of thermodynamics and kinetics, which can simulate the coupling effect of temperature and sediment density flow, and can scientifically and reasonably reflect the spatiotemporal evolution of the water temperature structure in the reservoir area of high-sediment-content river reservoirs in cold regions and quantify the impact period and degree of low-temperature and high-temperature water discharged after reservoir operation;

[0109] (2) The present invention adopts the method of using the median particle size of the group to determine the sediment settling velocity of the group, which better reflects the effect of fine-grained sediment on the water temperature structure of the reservoir area and the downstream water temperature process; the finer-grained sediment has a small settling velocity and good flow-following property, and produces a greater disturbance to the flow field and temperature field in the reservoir and in front of the dam by diving, forming a large circulation inside the reservoir and producing a "crowding-out" effect on the water body retained at the bottom of the reservoir, and the water temperature in the reservoir area shows an inverted temperature distribution structure of "the bottom is higher than the middle and upper layers in spring and summer"; in the winter when there is sediment inflow in cold-region reservoirs, the cold water entering the reservoir will dive due to the presence of sediment, rather than the surface buoyancy flow process of the cold water, further affecting the spatiotemporal process of the ice condition in the reservoir; in rivers with high sediment content and wide particle size distribution, using the median particle size of the group to calculate the sediment settling velocity can more truly reflect the spatiotemporal process of reservoir heat transfer under the coupling of water and sediment;

[0110] (3) The present invention can accurately obtain the sediment movement process in different types of sand-containing river reservoirs in cold regions and the thermal evolution process and laws of reservoir water temperature and ice conditions by considering the coupling effect of hydrodynamics, water temperature, ice and sand. BRIEF DESCRIPTION OF THE DRAWINGS

[0111] Figure 1 This is a block diagram of the water temperature-sediment-ice coupling model of the present invention;

[0112] Figure 2 Schematic diagram of the sediment density flow simulation process involved in the present invention;

[0113] Figure 3 This is a schematic diagram of the reservoir model grid division;

[0114] Figure 4 In the calculation process of a test to simulate the impact of sediment inflow into the reservoir, the sediment treatment used a single group and three groups of water temperature changes over time and compared them with the situation without sediment;

[0115] Figure 5In the calculation process of a test on the impact of sediment inflow into the reservoir, the sediment treatment adopts the change of sediment content over time in a single group and three groups;

[0116] Figure 6 The simulation results show that the reservoir water temperature changes in the whole year in three groups are compared with the situation when there is no sediment (representative months in each season);

[0117] Figure 7 The simulation results of the reservoir sediment grouping throughout the year were obtained by using three groups of water temperature changes in front of the dam throughout the year and the comparison with the sediment-free period (15th of each month);

[0118] Figure 8 The simulated reservoir sediment inflow was divided into three groups throughout the year, and the discharge water temperature process without sediment was compared with the natural water temperature process at the dam site (without a reservoir).

[0119] Figure 9 To simulate the reservoir sediment entering the reservoir throughout the year, three groups of water temperature changes during the winter ice period were used to compare with the situation when there was no sediment. DETAILED DESCRIPTION

[0120] The present invention is described in detail below through embodiments and application examples. It is necessary to point out that the present embodiments are only used to further illustrate the present invention and cannot be understood as limiting the scope of protection of the present invention. Those skilled in the art in this field can make some non-essential improvements and adjustments based on the above-mentioned contents of the present invention.

[0121] Example 1

[0122] This embodiment uses the inversion and prediction method of reservoir sediment and temperature dual density flow provided by the present invention to predict the water temperature distribution and ice cover distribution after the construction of a reservoir on the Yellow River.

[0123] The inversion and prediction method of reservoir sediment and temperature dual density flow provided in this embodiment includes the following steps:

[0124] Step 1: Establish a water temperature-sediment-ice coupling model, such as Figure 1 The water temperature-sediment-ice coupled model includes a hydrodynamic module, a sediment module, a water temperature module, and an ice module.

[0125] The hydrodynamic module includes the water flow continuity equation (Formula (1-1)), the momentum equation (Formulas (1-2)-(1-3)), the turbulence equation (Formulas (1-4)-(1-5)), the water state equation for water density (Formulas (1-6)-(1-8)), and the free water surface equation obtained by vertical integration of the continuity equation (Formula (1-9)).

[0126] The sand module includes suspended matter convection and diffusion equations after considering sediment settling (Formulas (2-1)-(2-2));

[0127] The water temperature module includes the heat transfer equation (Formula (3-1));

[0128] The ice module includes the initial ice sheet formation equation (Formula (4-1)) and the ice sheet thermal generation and disappearance equation (Formula (4-2)).

[0129] Step 2: Construct a cross-sectional model of the reservoir along the flow and depth directions, divide it into several grid cells, and set boundary conditions and initial conditions.

[0130] In this embodiment, unequally spaced rectangular grids are used to perform grid division based on the reservoir's large cross-section, cross-section spacing, horizontal grid spacing distribution, and vertical grid spacing distribution. The rectangular grid size is set according to the reservoir's characteristics. In this embodiment, the reservoir is discretized into 382×62 rectangular grids, with a longitudinal (along the mainstream direction) grid size of 20 to 800 m and a vertical (along the water depth direction) grid size of 2 m. The grids after division are as follows: Figure 3 The water temperature and surface ice cover conditions of each grid cell are used to represent the water temperature and ice cover conditions at the location.

[0131] Reservoir boundary conditions include inflow flow, inflow water temperature, inflow sediment load and gradation, meteorological conditions, and outflow flow. Initial conditions are given flow, temperature, and sediment fields. The initial ice thickness is zero, meaning there is no ice cover on the reservoir surface. The water temperature at the inlet boundary is the measured or given temperature at the reservoir tail, and the velocity is assumed to be uniform. Fully developed turbulent flow is assumed at the dam body and outlet section, and no-slip boundary conditions are applied to the dam surface, which is adiabatic. A no-slip boundary condition is applied to the reservoir bottom.

[0132] Step 3 solves the water temperature-sediment-ice coupling model to obtain the spatiotemporal distribution of sediment content, water temperature, and ice conditions in the reservoir area.

[0133] The spatiotemporal distribution of hydraulic power includes the temporal changes of flow velocity and water surface elevation; the spatiotemporal distribution of sediment content is the temporal changes of sediment content in the reservoir area; the spatiotemporal distribution of water temperature is the temporal changes of water temperature in the reservoir area; the spatiotemporal distribution of ice conditions is the temporal changes of surface ice thickness in the reservoir area.

[0134] In this step, for all the models constructed in the above step 1, the finite volume method is used to discretize each model to obtain a discrete equation group of each model on the grid unit, and then the discrete equation group is solved.

[0135] In this step, the hydraulic, sediment content, water temperature, and ice condition parameters obtained at the previous moment are assigned to various variables as the initial values of the hydraulic, sediment content, water temperature, and ice condition parameters at the current moment (for the initial moment, the hydraulic, sediment content, water temperature, and ice condition at the current moment correspond to the initial conditions). The hydraulic, sediment content, water temperature, and ice condition of the reservoir area at each moment are solved through time iteration, thereby obtaining the spatiotemporal distribution of hydraulic, sediment content, water temperature, and ice condition of the reservoir area.

[0136] This step includes the following sub-steps:

[0137] Step 3.1: Based on the spatial distribution of hydraulic pressure, sediment content, and ice conditions at the previous moment, the hydraulic module is used to preliminarily obtain the hydraulic spatial distribution at the current moment.

[0138] In this step, the initial values of the hydraulic, sediment content, water temperature, and ice condition parameter variables at the current moment are set based on the spatial distribution of hydraulic, sediment content, and ice conditions at the previous moment. Combined with the boundary conditions, the water flow continuity equation, momentum equation, turbulence equation, water body state equation, and free water surface equation of the hydrodynamic module are discretely solved to obtain the flow velocity (U and W) and water surface elevation (η) at the current moment, that is, the hydraulic spatial distribution at the current moment.

[0139] Step 3.2 obtains the spatial distribution of sediment content at the current moment through the sediment content module based on the spatial distribution of sediment content at the previous moment and the hydraulic spatial distribution at the current moment.

[0140] In this step, the initial value of the sediment content at the current moment is set according to the spatial distribution of the sediment content at the previous moment. Combined with the boundary conditions and the hydraulic spatial distribution at the current moment, the suspended sediment convection and diffusion equations at each level of the sand module are discretely solved to obtain the spatial distribution of sediment at different levels at the current moment (C i ), and then the spatial distribution of sediment content at the current moment is obtained.

[0141] Step 3.3 obtains the water temperature spatial distribution at the current moment through the water temperature module based on the water temperature spatial distribution at the previous moment and the hydraulic spatial distribution at the current moment.

[0142] In this step, the initial value of the water temperature at the current moment is set according to the spatial distribution of the water temperature at the previous moment. Combined with the boundary conditions and the hydraulic spatial distribution at the current moment, the heat transfer equation of the water temperature module is discretely solved to obtain the spatial distribution of the water temperature at the current moment (T w ).

[0143] Step 3.4 traverses the surface grid cells of the reservoir to determine whether the water temperature of the surface grid cells has dropped below freezing at the current moment. If so, the initial ice sheet thickness is determined using the initial ice sheet formation equation, and the surface water temperature is set to freezing point and the properties of the surface grid cells are changed from ice-free to ice-covered. Then, proceed to step 3.5. Otherwise, proceed directly to step 3.5.

[0144] In this step, each grid cell on the surface of the reservoir is traversed to determine whether the water temperature of the surface grid cell has dropped below the freezing point at the current moment. If so, the initial ice sheet thickness is calculated according to the initial ice sheet formation equation given above, and the surface water temperature is set to the freezing point and the properties of the surface grid cell are changed from the ice-free state to the ice-covered state. Then, the process proceeds to step 3.5. Otherwise, the process proceeds directly to step 3.5.

[0145] Step 3.5 traverses the surface grid cells of the reservoir, calculates the ice thickness change value of the surface grid cells with existing ice thickness, and determines the growth and decline process of the ice cover. If the ice thickness change value is positive, the ice thickness is growing; if the ice thickness change value is negative, the ice thickness is melting. When the ice thickness at the current moment is less than the melting ice thickness change value, the current surface grid cell state changes from ice-covered state to ice-free state.

[0146] In this step, the initial ice thickness for the current moment is set based on the spatial distribution of ice conditions at the previous moment. Combined with boundary conditions and the spatial distribution of water temperature at the current time, the ice module's ice cover thermal generation and dissipation equation is discretely solved to obtain the surface ice thickness (θ) for each surface grid cell with existing ice thickness at the current moment. The ice thickness change value is then subtracted from the current surface ice thickness of each surface grid cell from the corresponding surface ice thickness at the previous moment. The ice thickness change value is used to determine the ice cover's growth and dissipation process: a positive ice thickness change indicates growth; a negative ice thickness change indicates ablation. When the current ice thickness is less than the ablation ice thickness change value, the current surface grid cell transitions from ice-covered to ice-free.

[0147] Increase the time step and repeat the above steps S3.1-3.5 until the upper limit of the set time period is reached, and the reservoir hydraulic temporal and spatial distribution, sediment content temporal and spatial distribution, water temperature temporal and spatial distribution, and ice condition temporal and spatial distribution of the reservoir area are obtained.

[0148] The convergence conditions of the hydrodynamic field, temperature field, and sediment field are:

[0149] (1) Under the condition of accurately given initial fields, the various field data obtained by calculation are the time series results to be obtained;

[0150] (2) When the initial field cannot be accurately given, the initial conditions can be determined by random initialization, and then multiple time cycles are used. When the result of the mth cycle and the result of the m-1th cycle are within the allowable difference range, the result of the mth cycle is considered to be the intended time series result.

[0151] The time step size is set to meet the requirements of formula (5) given above.

[0152] The following is a detailed explanation of the water temperature distribution of the Yellow River reservoir and the ice cover distribution inversion and situation mentioned above.

[0153] Since the sediment settling velocity characteristics will significantly affect the sediment distribution in the reservoir area and thus affect the water temperature structure, it is particularly critical to reasonably deal with the sediment settling velocity of the group. Figure 4 The paper presents the comparison of the water temperature changes with time in the process of the single sediment inflow test calculation, using a single sediment treatment group (0-25μm sediment median particle size) and three sediment treatment groups (0-5μm, 5-10μm, 10-25μm three-level sediment median particle size, the mass ratio of the three levels of sediment being 3:1:3) and the case without sediment. Figure 5 A two-dimensional distribution diagram of reservoir sediment content is given using the sediment grouping (3 groups) median particle size calculation method and the single group median particle size calculation method. The calculation starts from July, the flood season with higher sediment content (there is no sediment influence before July), under the given same initial and boundary conditions, and the simulation results on the 6th, 12th, 18th and 24th days are given. The median settling velocity of a single group of heterogeneous sand was calculated to be 3 m / day. The numerical results show that the vertical diffusion of sediment is inhibited. The large density difference between clear and turbid water causes the clear water layer to form a continuous vortex during the infiltration of the density flow. The vortex drives the turbid water with higher local temperature and higher sediment content to flow in the non-mainstream direction. The sediment is divided into three groups according to particle size from small to large, and the corresponding settling velocity of each group is calculated. Among them, the sediment group with smaller particle size has a low settling velocity (0.38 m / day), good flow-following ability, and is easy to diffuse vertically, which inhibits the formation of vortices. The temperature field and sediment field are more uniform. At the same time, fine-grained sediment is more likely to reach the front of the dam, and the vertical water temperature and sediment content in front of the dam increase faster. Due to the varying particle sizes of sediment entering a reservoir, it is sorted and deposited under the influence of water flow. The median particle size of the sediment decreases from the reservoir tail to the dam front. Simultaneously, the density flow is blocked by the dam, resulting in a turbid reservoir. Generally, the interface between clear and turbid water is clearly defined and nearly parallel to the water surface. The results calculated after particle size grouping generally conform to these characteristics. Fine sediment is the primary factor affecting the water temperature and flow field structure in front of the dam. Therefore, in this paper, the median particle size of the sediment particle size grouping is used to calculate the settling velocity.

[0154] Figure 6The comparison between the two-dimensional distribution of water temperature in the reservoir area predicted by the present invention and the one without sediment is given. Figure 7 To correspond to the vertical water temperature distribution in front of the dam (the 15th of each month). Overall, the present invention better simulates the coupling effect of reservoir sediment density flow and temperature density flow. The water temperature structure of the reservoir is mainly controlled by the sediment content of the incoming water body. Except for May and June, the inflowing water body in each month is mainly diving. During the flood season, the incoming high-sediment hot water basically replaces the cold water in the reservoir area. The suspended sediment increases the density of the water body at the bottom of the reservoir, causing the vertical mixing in front of the dam to weaken in winter, and the hot water stagnates in front of the dam. During the warming period, the high-temperature inflow water dives, supporting the low-temperature water at the bottom of the reservoir to exit the reservoir, and forming an inverted temperature distribution. The overall vertical stratification is weak, and the vertical temperature difference is 1.2 to 13.7 ° C. When there is no sediment, the water temperature in the reservoir area undergoes a major reversal in March and November, forming an isothermal distribution. When there is sediment, in addition to forming an isothermal distribution in March, the high-temperature water entering in August during the flood season gradually sinks and replaces the low-temperature water in the reservoir area, forming an isothermal distribution. In November, due to the influence of the sinking of high-temperature water during the flood season, a certain vertical temperature difference is always maintained in front of the dam.

[0155] Figure 8 The predicted reservoir discharge temperature process is compared with the natural water temperature process at the dam site. The present invention effectively simulates the flattening and delayed effects of the discharge water temperature under the influence of dual hyperpycnal flow. Compared with the natural water temperature process at the dam site, low-temperature water discharge occurs from March to July, and high-temperature water discharge occurs from September to February. Under the influence of sediment, the maximum low-temperature water amplitude is 10.2°C. The low-temperature water discharge is more significant from May to August, with the amplitude increasing in July. This is mainly due to the increased sediment content of the incoming water during this period, the downward flow of hyperpycnal flow, and the accelerated discharge of low-temperature water stored at the reservoir bottom, resulting in a significant drop in the discharge water temperature. Because the high-sediment inflow changes the water temperature structure of the reservoir, especially in the summer, the water in the stagnant temperature layer at the reservoir bottom is essentially replaced by high-temperature water. After entering the cooling period, the inflow and sediment content decrease, the high-temperature water discharge is slow, and the reservoir bottom water temperature is relatively high in winter, which leads to a large amplitude of high-temperature water in winter, with a maximum amplitude of 11.1°C.

[0156] Figure 9 The present study compares the reservoir's ice-age temperature process predicted using this method with that observed in the absence of sediment. Winter ice conditions in the reservoir are weakened due to the inflow of sediment-laden water. Sediment alters the water's buoyancy pattern, causing low-temperature water to sink with the sediment's density flow. Vertical flow contributes less to surface water cooling, resulting in slower surface water temperature drop compared to the absence of sediment, a slower ice cover advance, and a smaller ice-covered area.

Claims

1. A method for inversion and prediction of reservoir sediment and temperature dual density flow, characterized by: The following steps are involved: Step 1: establishing a water temperature-sediment-ice coupling model; the water temperature-sediment-ice coupling model includes a hydrodynamic module, a sand module, a water temperature module, and an ice module; The hydrodynamic module includes the water flow continuity equation, momentum equation, turbulence equation, water state equation about water density, and free water surface equation obtained by vertical integration of the continuity equation; The sand module includes suspended matter convection and diffusion equations after sediment settling; The water temperature module includes a heat transfer equation; The ice module includes an initial ice sheet formation equation and an ice sheet thermal generation and disappearance equation; Step 2: Construct a cross-sectional model of the reservoir along the flow and depth directions, divide it into several grid cells, and set boundary conditions and initial conditions; Step 3 solves the water temperature-sediment-ice coupling model to obtain the spatiotemporal distribution of sediment concentration, water temperature, and ice conditions in the reservoir area. This step includes the following sub-steps: Step 3.1 Based on the spatial distribution of hydraulic pressure, sediment content, and ice conditions at the previous moment, the hydraulic module is used to preliminarily obtain the hydraulic spatial distribution at the current moment; Step 3.2: Based on the spatial distribution of sediment content at the previous moment and the hydraulic spatial distribution at the current moment, the spatial distribution of sediment content at the current moment is obtained through the sediment content module; Step 3.3: Based on the water temperature spatial distribution at the previous moment and the hydraulic spatial distribution at the current moment, obtain the water temperature spatial distribution at the current moment through the water temperature module; Step 3.4 traverses the surface grid cells of the reservoir to determine whether the water temperature of the surface grid cells has dropped below freezing point at the current moment. If so, the initial ice sheet thickness is determined using the initial ice sheet formation equation, and the surface water temperature is set to freezing point and the properties of the surface grid cells are changed from ice-free to ice-covered. Then, proceed to step 3.

5. Otherwise, proceed directly to step 3.

5. Step 3.5 traverses the surface grid cells of the reservoir and calculates the ice thickness change value of the surface grid cells with existing ice thickness to determine the growth and decline process of the ice cover; if the ice thickness change value is positive, the ice thickness is growing; If the ice thickness change value is negative, the ice thickness is melting. When the ice thickness at the current moment is less than the melting ice thickness change value, the current surface grid cell state changes from ice-covered state to ice-free state. Increase the time step and repeat steps 3.1-3.5 above until the upper limit of the set time period is reached to obtain the spatiotemporal distribution of hydraulic pressure, sediment content, water temperature, and ice conditions in the reservoir area.

2. The method for inversion and prediction of reservoir sediment and temperature dual density flow according to claim 1 is characterized in that: In step 1, the water flow continuity equation is: （1-1） Where: U is the longitudinal velocity; x is the flow direction; W is the vertical velocity; z is the vertical direction of the reservoir; q is the net inflow flow per unit width of the lateral grid unit; B is the horizontal width of the reservoir at different elevations along the flow direction; The momentum equation includes the longitudinal momentum equation and the vertical momentum equation; The longitudinal momentum equation is: （1-2） Vertical momentum equation: （1-3） Where: α is the riverbed slope; is the turbulent normal stress, is the turbulent shear stress; η is the water surface elevation; U x is the x-direction component of the tributary velocity in the direction of water flow; g is the acceleration of gravity; ρ is the water density; p is the hydrostatic pressure; Turbulence equation: （1-4） （1-5） Where: k is the turbulent kinetic energy; ε is the turbulent kinetic energy dissipation rate; ν t is the turbulent eddy viscosity coefficient, ν t =C μ k 2 / ε; ρ is the water density; P is the turbulent kinetic energy generation term, ; G is the buoyancy term, ;P k 、P ε is the turbulent kinetic energy and dissipation rate generated by boundary friction; σ k , σ ε are the turbulent kinetic energy and the dissipation rate Prandtl number, respectively; is the turbulent Prandtl number; C μ 、C ε1 and C ε2 is an empirical constant; The specific form of the water state equation is as follows: (1-6) (1-7) (1-8) Where: T w is the water temperature; is the density of water under the influence of temperature; is the water density increment caused by the i-th level of suspended sediment; C i is the sediment content per unit volume of water; γ s,i is the bulk density of the ith level of sediment particles; n represents the number of sediment level groups; The free water surface equation is: (1-9) Where: is the lateral width of the reservoir surface along the flow direction; η is the water surface elevation; h is the riverbed elevation.

3. The inversion and prediction method of reservoir sediment and temperature dual density flow according to claim 2 is characterized in that: The suspended sediment is divided into i levels, and the suspended sediment convection and diffusion equations at each level are: （2-1） Where: C i is the sediment content per unit volume of water; ω s,i is the settling velocity of a single sediment particle; 、 are the longitudinal turbulent diffusion coefficient and the vertical turbulent diffusion coefficient respectively; The formula for the sediment settling velocity of a single particle uses the commonly used Stokes formula as follows: （2-2） Where: ω s,i is the settling velocity of a single sediment particle; d i is the median particle size of suspended sediment particles at level i; υ is the kinematic viscosity of water; is the bulk density of water without sand, ; is the dry bulk density of suspended sediment at the i-th level, , is the density of suspended sediment at the i-th level.

4. The method for inversion and prediction of reservoir sediment and temperature dual density flow according to claim 3 is characterized in that: The heat transfer equation is: （3-1） Where: T w is the water temperature; Q is the rate of lateral heat flux in and out of the grid unit, , is the water temperature of the lateral flow; S is the source-sink term, including surface heat exchange, riverbed heat exchange, ice-water heat exchange and / or solar radiation absorption at water depth z.

5. The inversion and prediction method of reservoir sediment and temperature dual density flow according to claim 4 is characterized in that: The initial ice sheet thickness is calculated as follows: （4-1） Where: θ0 is the initial ice thickness; is the surface water temperature of the reservoir; C p is the specific heat capacity of water at constant pressure; θ s is the vertical grid thickness of the surface layer; ρ is the water density; ρ i is the density of ice; L f is the latent heat of freezing; The ice sheet thermal generation and dissipation equation for ice thickness θ is given by the following equation: （4-2） Where: θ is the ice thickness; T m is the ice-water interface temperature; T a is the air temperature; k i is the thermal conductivity of ice; h ai is the ice-air heat exchange coefficient; is the heat exchange flux between the reservoir surface and the bottom of the ice sheet.

6. The method for inversion and prediction of reservoir sediment and temperature dual density flow according to claim 1, characterized in that: In step 2, the boundary conditions include inflow flow, inflow water temperature, inflow sediment volume and gradation, meteorological conditions, and outflow flow. The initial conditions are the given flow field, temperature field, and sediment field. The initial ice thickness is 0, that is, there is no ice cover on the reservoir surface.

7. The method for inversion and prediction of reservoir sediment and temperature dual density flow according to claim 1, characterized in that: The spatiotemporal distribution of hydraulic power includes the temporal changes of flow velocity and water surface elevation; the spatiotemporal distribution of sediment content is the temporal changes of sediment content in the reservoir area; the spatiotemporal distribution of water temperature is the temporal changes of water temperature in the reservoir area; the spatiotemporal distribution of ice conditions is the temporal changes of surface ice thickness in the reservoir area.

8. The method for inversion and prediction of reservoir sediment and temperature dual density flow according to claim 1, characterized in that: The time step meets the following formula conditions: ； Where, is the time step; is the longitudinal eddy viscosity; is the vertical eddy viscosity; is the total flow into or out of the grid cell; V is the volume of the grid cell; H is the depth of the water column where the grid cell is located; ρ is the water density; is the maximum density difference in the water column where the grid cell is located.

9. The method for inversion and prediction of reservoir sediment and temperature dual density flow according to any one of claims 1 to 8, characterized in that: In step 3.2, the suspended matter convection-diffusion equation is used to calculate the sediment content distribution of different levels and further update the water density.

Images