Method for simulating heat flux of flooding lake wetland in watershed hydrological model

Abstract

The invention discloses a flood lake wetland heat flux simulation method in a watershed hydrological model, and relates to the technical field of heat flux simulation, and the method comprises the following steps: dispersing a watershed space into a slope river channel and a lake wetland, simulating the water movement of a slope river channel unit, and simulating the water movement of the lake wetland; the method comprises the following steps: simulating surface energy exchange for a lake wetland unit on the basis of simulating a hydrological process, simulating land hydrology and lake inflow runoff of a sub-basin by adopting a grid-based distributed hydrological model, respectively dividing simulation grid units for a slope river channel and a lake wetland, establishing evapotranspiration, infiltration and surface runoff models, and establishing an evapotranspiration, infiltration and surface runoff model; a two-dimensional shallow water equation is used for processing slope confluence, a one-dimensional hydrodynamic model is used for calculating river channel confluence, a lake water quantity balance model is established for calculating the lake water level, the lake submerging process is simulated, and the lake wetland distributed heat flux process is simulated. And the technical blank of flooding wetland hydrothermal simulation is filled.

Description

Technical Field

[0001] The present invention relates to the technical field of heat flux simulation, and particularly to a method for simulating heat flux of floodplain lake wetlands in a watershed hydrological model. Background Technique

[0002] Wetlands are special types of underlying surfaces in the transitional zone between land and water, known as the "kidneys of the earth". Their water and heat flux processes are of great significance for regulating regional climate, influencing the internal water and heat circulation in wetlands, and maintaining the health of the ecosystem. Floodplain lake wetlands formed due to the periodic wet-dry alternation of the hydrological rhythm account for about 15% of the total wetland area in the world. They are extremely important wetland types in the wetland ecosystem and also key interface systems for the exchange of water and energy between the earth's surface and the atmosphere, playing a key role in the material and energy balance of the land-air coupling system. The large amplitude of the high water level causes drastic changes in the wetland inundation dynamics, resulting in more complex interface properties and land surface parameters of floodplain lake wetlands compared to ordinary wetlands. Accurately simulating the dynamic change process of the heat flux of floodplain lake wetlands is of great significance for deeply understanding the interaction between the earth and the atmosphere, exploring the regulation mechanism of floodplain wetlands on regional climate, and formulating wetland protection strategies.

[0003] As an important tool for studying land surface hydrological processes, hydrological models can simulate the water and energy processes at the watershed scale through modules such as surface water, groundwater, soil water, and lake hydrology. Most hydrological models focus on the runoff generation and confluence processes, with the simulation target being the flow process at the outlet river cross-section. For complex lake basins composed of sub-watersheds - rivers - lakes, especially large floodplain lakes, there is still a lack of a general hydrological model with strong applicability, high stability, and certain accuracy. Although some domestic and foreign hydrological models have developed modules for water and heat fluxes of lake wetlands, these parameterization schemes are all for static wetlands with long-term inundation, and the underlying surface parameterization schemes still use fixed water surfaces or land as boundaries, unable to depict the impact of large-scale inundation and exposure dynamics caused by the floodplain process of large lakes on the regional energy balance. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for simulating heat flux of floodplain lake wetlands in a watershed hydrological model to solve the problems proposed in the above background technique.

[0005] To solve the above technical problems, the present invention provides the following technical solutions: A method for simulating heat flux of floodplain lake wetlands in a watershed hydrological model, comprising the following steps:

[0006] S1: Discretize the watershed space into two parts: slope channels and lake wetlands;

[0007] S2: Divide grid cells for slope channels and lake wetlands respectively;

[0008] S3: Establish evapotranspiration, infiltration, and surface runoff models, and analyze the main hydrological processes of evaporation, infiltration, runoff generation, and confluence based on the models;

[0009] S4: Establish a lake water balance model to calculate the lake water level and calculate the flow rate at the lake outlet discharge;

[0010] S5: Simulate the lake inundation process and conduct a distributed heat flux process simulation for the lake wetland.

[0011] Furthermore, in step S1, the watershed space is discretized into two parts: slope channels and lake wetlands. For slope channel units, simulate the water movement process; for lake wetland units, on the basis of simulating the hydrological process, simultaneously simulate the surface energy exchange process. The slope channel units and lake wetland units exchange water volume through in-lake or out-lake channels, and there will be no overland flow across the lake boundary.

[0012] Furthermore, in step S2, grid cells are divided for the three parts of slopes, channels, and lake wetlands respectively. For slopes and lake wetlands, square grid cells are divided according to the established spatial resolution (1 km - 25 km); for channels, they are divided into rectangular grid cells according to the channel shape. Generally, the average channel width is set as the grid width, and the grid length can be set according to the simulation accuracy of the confluence process (0.1 - 5 km); the slopes and channels are divided into coarse grid cells, and the lake wetland part is divided into fine grid cells. The resolutions of the coarse grid cells and fine grid cells are in an integer multiple relationship, so as to achieve the nesting between simulation regions.

[0013] Furthermore, in step S3, the runoff generation and confluence process of the sub-watershed includes the runoff generation process at the slope scale and the confluence process at the slope channel scale. The runoff generation process includes precipitation, canopy interception, vegetation transpiration, soil evaporation, surface water infiltration, vertical movement of soil moisture, and soil water - groundwater exchange. The confluence process includes overland flow and subsurface flow sub - modules at the slope scale, and a channel confluence sub - module at the channel scale, and finally the inflow into the lake is obtained. The slope is directly connected to the channel, and the slope provides lateral inflow for channel confluence directly through overland flow and subsurface flow; based on net radiation, saturation vapor pressure, actual vapor pressure, actual pressure of water vapor, and wind speed, establish an evapotranspiration process model for the evapotranspiration amount ET:

[0014]

[0015] where Δ is the slope of the saturation vapor pressure curve (kPa / ℃), which is related to temperature, R n is the net radiation, with the unit of MJ / m 2·d represents the net radiant energy per unit area; γ is the dry air constant, with the unit of kPa / ℃, and usually takes the value of 0.066 kPa / ℃; es is the saturated water vapor pressure, with the unit of kPa, representing the maximum pressure of water vapor in the air at a specific temperature; ea is the actual water vapor pressure, with the unit of kPa, representing the actual pressure of water vapor in the air under specific conditions; R is the wind speed;

[0016] Based on the saturated hydraulic conductivity of the soil, the soil water head, the pressure head of water in the soil, the soil water potential, and the infiltration surface area, an infiltration process model for the total infiltration amount I(t) is established:

[0017]

[0018] where f(t) is the infiltration rate varying with time; K s is the saturated hydraulic conductivity of the soil (m / s); h is the soil water head (m), representing the pressure head of water in the soil; φ is the soil water potential (m), usually representing the water absorption capacity of the soil; A1 is the infiltration surface area; I(t) is the total infiltration amount during the time period, obtained by integration;

[0019] Based on the precipitation, evapotranspiration, change in soil moisture, and exchange between soil water and groundwater, the surface runoff Q1 is calculated:

[0020] Q1 = P1 - ET - ΔS1 + L;

[0021] where Q1 is the surface runoff, P1 is the precipitation, ET is the evapotranspiration, ΔS1 is the change in soil moisture, and L is the exchange between soil water and groundwater;

[0022] The two-dimensional shallow water equation is used to handle the calculation of overland flow concentration, and its vector form is:

[0023]

[0024] where U1 is the unknown quantity to be solved in the two-dimensional shallow water equation, E1 and G1 are the two-dimensional directional fluxes; S1 is the source term in the two-dimensional shallow water equation;

[0025] The one-dimensional hydrodynamic model is used to calculate the river flow concentration, and the governing equations adopt the Saint-Venant equations, and its vector form is:

[0026]

[0027] In the formula, U2 is the unknown quantity to be solved, F is the flux, and S2 is the source term. The finite volume method is used for second-order accurate spatial discretization of the governing equations, and the second-order Runge-Kutta explicit scheme is used for time discretization. The numerical flux between grids is obtained by calculating with an approximate Riemann solver. The bottom slope term in the source term is calculated based on the hydrostatic reconstruction method and central difference discretization. The friction term is treated fully implicitly and calculated using the Newton-Raphson iteration to obtain better numerical stability. The model calculation step size adopts a conditional adaptive time step size to ensure the stability of the calculation.

[0028] Further, in step S4, considering the lake water volume change, lake surface precipitation, lake surface evaporation, lake water surface area, total inflow river flow, and flow at the lake outlet, a lake water volume balance model regarding the lake water level in the previous simulation period and the lake water level in the current period is established:

[0029] ΔV = (P - E) * A + R in +R out

[0030]

[0031] where ΔV is the lake water volume change, P1 is the lake surface precipitation, E is the lake surface evaporation, A2 is the lake water surface area, R in is the total inflow river flow, and R out is the flow at the lake outlet. h1 and H2 are the lake water levels in the previous simulation period and the current period respectively. The relationship between the lake surface area and the water level can be obtained through the water level - area curve; the total inflow river flow includes the sum of the runoff of all inflow rivers;

[0032] R out can be calculated by the following formula. Based on the assumption of broad-crested weir flow, assuming that the velocity head can be ignored, the discharge R out is calculated:

[0033]

[0034] where b is the flow width (m), g is the acceleration due to gravity, z is the current lake depth (m), and z min is the elevation (m) above the weir or the lake outlet. The discharge coefficient c d is used to consider the inflow velocity, non-parallel streamlines at the top, and energy losses. The value of c d varies approximately between 0.8 and 1.2.

[0035] Further, in step S5, the seed flooding algorithm is as follows: The seed flooding algorithm starts from a seed point and marks adjacent pixel regions as the same type or color by gradually diffusing and filling pixels. The seed flooding algorithm is used in combination with DEM to simulate the active flooding process of a shallow lake. The steps of simulating the lake flooding process based on the seed flooding algorithm include: determining the location and flow data of the river flowing into the lake, initialization, starting the seed flooding, and updating the flooding process;

[0036] Determining the location and flow data of the river flowing into the lake: Determine the injection point location of the river flowing into the lake and the corresponding flow data. The injection point is regarded as the seed of the lake basin, and the flooding process is simulated starting from the injection point;

[0037] Initialization: Take the injection point as the seed point, and take the flow rate at the seed point as the initial flooding flow rate. At the same time, load the lake basin DEM data and the initial water level data into the computing environment;

[0038] Starting the seed flooding: Starting from the seed point, search and diffuse to the surrounding 8 neighborhood grids according to the eight-neighborhood connectivity domain, and layer by layer judge whether each grid meets the flooding conditions, and add the qualified grids to the queue to be processed; Continuously iterate this process until no new grids are flooded;

[0039] Updating the flooding process: According to the flooding elevation values of each time period, a dynamic image of the flooding process can be generated to show the process of the flooding range gradually expanding, so as to update the flooding value H new :

[0040]

[0041] The H new is the updated flooding elevation value, H old is the flooding elevation value of the current grid, Q3 is the injection flow data, and A3 is the influence range of the injection flow;

[0042] The unit grids are divided into two categories: lake grids and non-lake grids, where the lake grids have seasonal changes in the flooding range during the simulation period;

[0043] Calculate the surface energy components of each grid using the energy balance equation and analyze the net radiation R n :

[0044] R n =H + ρ w *λ v *E0 + G;

[0045] Among them, R n is the net radiation (W*m -2 ), H is the sensible heat flux (W*m -2 ), ρ w *λv *E is the latent heat flux (W·m -2 )(ρ w is the liquid water density, with the unit of kg·m -3 ; λ v is the latent heat of vaporization of water, with the unit of J·kg -1 ), E0 is the evapotranspiration, and G is the geothermal flux (W·m -2 );

[0046] Based on the surface albedo, downward shortwave radiation, surface emissivity, and downward longwave radiation, the net radiation R input to the grid is obtained n :

[0047]

[0048] where α is the surface albedo of this land cover type, R s is the downward shortwave radiation (W·m -2 ), ε is the surface emissivity of this land cover type, R L is the downward longwave radiation (W·m -2 ), T s is the surface temperature (K), and σ is the Stefan-Boltzmann constant (5.67×10-8W·m -2 ·K -4 );

[0049] Based on the aerodynamic resistance of the heat flux, air density, specific heat capacity of air, surface temperature, and surface air temperature, the sensible heat flux H is calculated as follows:

[0050]

[0051] where α is the grid inundation area ratio, T W is the water surface temperature, T g is the soil surface temperature of the exposed part of the grid, T a is the air temperature, ρ a is the air density, c p is the specific heat capacity of air, r h,w is the aerodynamic resistance of the water surface, r h,g is the aerodynamic resistance of the land surface, r s,w is the surface resistance of the water body, r s,g is the soil evaporation impedance.

[0052] Achieve continuous transition of water-land energy fluxes through the submergence ratio, and accurately characterize the heat exchange behavior of "partially submerged" grids; the water temperature calculation includes a heat storage term to reflect the heat buffering effect during flood retention; the introduced term includes the influence of the air pressure gradient force to avoid overestimation of latent heat fluxes in high-humidity environments; the land surface resistance term is related to soil moisture, and the model's feedback on evaporation suppression during dry periods is more significant.

[0053] Derivation process: Assume that there are two types of land surfaces within the grid, with the water area proportion being α and the surface temperature being T W ; the soil area proportion is 1 - α, and the surface temperature is T g . The total sensible heat flux needs to sum the two according to the area weights. According to the law of conservation of energy, the greater the resistance, the smaller the sensible heat flux. The sensible heat flux of the water surface is driven by the temperature difference T W -T a , and the resistance includes the aerodynamic resistance r h,w (reflecting the influence of air turbulence on heat transfer) and the water surface resistance r s,w (possibly introduced due to evaporation suppression or water surface characteristics). The formula for calculating the sensible heat flux of water is:

[0054]

[0055] Similarly, for the sensible heat flux of the soil, there are two resistance terms: the aerodynamic resistance r h,g (affected by soil roughness on turbulence) and the soil evaporation impedance r s,g (related to soil moisture, with higher impedance for dry soil).

[0056]

[0057] By combining the contributions of water and soil according to the area weights, the formula for calculating the sensible heat flux of the water-soil mixed grid can be obtained:

[0058]

[0059] Based on the soil thermal conductivity, the soil temperature between the first and second soil layers, and the thickness of the first soil layer, the geothermal flux G of the topsoil is calculated as:

[0060]

[0061] where k is the soil thermal conductivity (W*m -1 *K -1 ), T1 is the soil temperature between the first and second soil layers (K), and D1 is the thickness of the first soil layer (m);

[0062] Based on the sensible heat flux, latent heat flux, changes in the heat storage in the overlying water body and soil of the floodplain wetland, the advective heat flux carried by water flow, and the heat flux entering the underlying soil, the surface net radiation R is calculated. n :

[0063] R n = H + LE + ΔS2 + A4 + Q B ;

[0064] Among them, R n is the surface net radiation, H is the sensible heat flux, LE is the latent heat flux, ΔS2 is the change in the heat storage in the overlying water body and soil of the floodplain wetland, A4 is the advective heat flux carried by water flow, and Q B is the heat flux entering the underlying soil. Among them, H and LE are obtained from flux observation data, and ΔS2 and A4 need to be estimated based on key parameters such as specific heat capacity and flow velocity in combination with the energy equation;

[0065] Based on the soil depth, soil temperature, soil specific heat capacity, and the heat flux at the reference depth z ref at, by integrating the one-dimensional heat diffusion equation, the change in the heat storage of the exposed soil ΔS(z) can be obtained:

[0066]

[0067] Among them, z is the soil depth (m), T is the soil temperature (K), ρ s c s is the soil specific heat capacity (J*kg -1 *K -1 ), and S(z ref ) is the heat flux at the reference depth z ref . If z ref is greater than 1 m and S(z ref ) is less than 1% of the surface soil heat flux, it can be assumed that S(z ref ) = 0;

[0068] Given the temperature profile T(z i ), based on the time interval of the discrete time step, the reference depth, the interval in the vertical direction, and the depth coordinate in the vertical direction, the soil heat flux G is calculated:

[0069]

[0070] Among them, Δt represents the time interval for the discrete time step; t represents the current moment; z ref represents the reference depth, Δz represents the interval in the vertical direction, and z i represents the depth coordinate in the vertical direction, and T (zi,t) is the temperature at a depth of Z i, the soil temperature corresponding to time t, the soil specific heat capacity can be calculated based on soil water content and soil porosity, and the temperature profile is estimated using the soil heat diffusion equation;

[0071] Based on the total water depth, average water temperature, and the change in average water temperature, the change in water body heat storage ΔS during the flooding period is calculated w :

[0072]

[0073] where ρ w c p is the specific heat capacity of water (J*kg -1 *K -1 ), z w is the total water depth, T w is the average water temperature, ΔT w is the change in average water temperature within the time period Δt. According to water temperature observations, when the water depth is less than 1 m, it is considered that the lake water temperature is evenly mixed; when the water depth is greater than 1 m, ΔT w profile is estimated through water thermal conductivity and surface temperature from remote sensing. The change in heat storage is estimated by weighted averaging the water heat storage change (ΔS w ) during the flood season and the soil heat flux (G) during the dry season; due to the drastic change in water conditions and the relatively fast surface water flow velocity of the overlying water body in the floodplain wetland, the horizontal advection flux of the overlying water body is calculated by the following formula:

[0074]

[0075] where u is the water flow velocity along the direction of temperature gradient measurement, x is the horizontal distance, is the horizontal gradient of water temperature, ρ w c p is the specific heat capacity of water.

[0076] Compared with the prior art, the beneficial effects achieved by the present invention are as follows: It overcomes the defect that traditional hydrological models are insufficient in depicting the hydrothermal processes of large lake wetlands, makes up for the problem that the parameterization scheme for simulating the hydrothermal processes of lake wetlands in existing hydrological models relies too much on complex hydrodynamic equations, realizes multi-level nested simulation of the basin-river-lake and fine grid encryption of lake wetlands, uses the seed spreading algorithm to simulate the flood inundation process, and drives the energy balance equation with the inundation dynamic variables, improving the model's simulation ability for the heat flux of flooded wetlands, quantifying the hydrothermal balance process of large lake wetlands, and enriching and perfecting the theoretical system of lake-atmosphere energy exchange and regional atmospheric hydrological cycle response. From the perspective of practical application, the present invention has the ability to improve the accuracy of regional hydrometeorological simulation and prediction from the physical mechanism. Compared with the existing hydrodynamic simulation scheme, it has the advantages of low calculation cost, less data requirement, simple model structure, stable calculation performance, and improved efficiency, filling the technical gap in the hydrothermal simulation of flooded wetlands. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] The accompanying drawings are used to provide a further understanding of the present invention, and constitute a part of the specification. They are used together with the embodiments of the present invention to explain the present invention, but do not constitute a limitation to the present invention. In the accompanying drawings:

[0078] Figure 1 is a flowchart of a method for simulating the heat flux of flooded lake wetlands in a basin hydrological model of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0079] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0080] Please refer to Figure 1 , the present invention provides a technical solution: a method for simulating the heat flux of flooded lake wetlands in a basin hydrological model, including the following steps:

[0081] S1: Discretize the basin space into two parts: the slope channel and the lake wetland;

[0082] S2: Divide grid cells for the slope channel and the lake wetland respectively;

[0083] S3: Establish evapotranspiration, infiltration, and surface runoff models, and analyze the main hydrological processes of evaporation, infiltration, runoff generation, and confluence according to the models;

[0084] S4: Establish a lake water balance model to calculate the lake water level and calculate the flow rate at the lake outlet discharge;

[0085] S5: Simulate the lake inundation process and conduct a simulation of the distributed heat flux process in lake wetlands.

[0086] In step S1, the watershed space is discretized into two parts: the hillslope channel and the lake wetland. For the hillslope channel unit, simulate the water movement process; for the lake wetland unit, on the basis of simulating the hydrological process, simultaneously simulate the surface energy exchange process. The hillslope channel unit and the lake wetland unit exchange water volume through the inflow or outflow channels, and there will be no overland flow across the lake boundary.

[0087] In step S2, grid cells are divided for the hillslope, channel, and lake wetland parts respectively. For the hillslope and lake wetland parts, square grid cells are divided according to the established spatial resolution (1 km - 25 km); for the channel part, it is divided into rectangular grid cells according to the channel shape. Generally, the average channel width is set as the grid width, and the grid length can be set according to the simulation accuracy of the confluence process (0.1 - 5 km); the hillslope and the channel are divided into coarse grid cells, and the lake wetland part is divided into fine grid cells. The resolutions of the coarse grid cells and the fine grid cells are in an integer multiple relationship, so as to achieve the nesting between the simulation areas.

[0088] In step S3, the runoff generation and concentration process of the sub - watershed includes the runoff generation process at the hillslope scale and the concentration process at the hillslope channel scale. The runoff generation process includes precipitation, canopy interception, vegetation transpiration, soil evaporation, surface water infiltration, vertical movement of soil moisture, and soil - water - groundwater exchange. The concentration process includes the overland flow and subsurface flow sub - modules at the hillslope scale, and the channel concentration sub - module at the channel scale. Finally, the runoff into the lake is obtained. The hillslope is directly connected to the channel, and the hillslope provides lateral inflow for the channel concentration directly through overland flow and subsurface flow; based on net radiation, saturation vapor pressure, actual vapor pressure, actual pressure of water vapor, and wind speed, establish an evapotranspiration process model for evapotranspiration amount ET:

[0089]

[0090] The Δ is the slope of the saturation vapor pressure curve (kPa / °C), which is related to temperature, R n is the net radiation, with the unit of MJ / m 2 ·d, representing the net radiation energy per unit area; γ is the dry air constant, with the unit of kPa / °C, usually taking a value of 0.066 kPa / °C; es is the saturation vapor pressure, with the unit of kPa, representing the maximum pressure of water vapor in the air at a specific temperature; ea is the actual vapor pressure, with the unit of kPa, representing the actual pressure of water vapor in the air under specific conditions; R is the wind speed;

[0091] Based on the saturated hydraulic conductivity of the soil, the water head of the soil, the pressure head of the water in the soil, the water potential of the soil, and the infiltration surface area, an infiltration process model for the total infiltration amount I(t) is established:

[0092]

[0093] where f(t) is the infiltration rate varying with time; K s is the saturated hydraulic conductivity of the soil (m / s); h is the water head of the soil (m), representing the pressure head of the water in the soil; φ is the water potential of the soil (m), usually representing the water absorption capacity of the soil; A1 is the infiltration surface area; I(t) is the total infiltration amount during the time period, obtained by integration;

[0094] Based on the precipitation, evapotranspiration, change in soil moisture, and exchange between soil water and groundwater, the surface runoff Q1 is calculated:

[0095] Q1 = P1 - ET - ΔS1 + L;

[0096] where Q1 is the surface runoff, P1 is the precipitation, ET is the evapotranspiration, ΔS1 is the change in soil moisture, and L is the exchange between soil water and groundwater;

[0097] The two-dimensional shallow water equation is used to handle the overland flow calculation, and its vector form is:

[0098]

[0099] where U1 is the unknown quantity to be solved in the two-dimensional shallow water equation, E1 and G1 are the two-dimensional directional fluxes; S1 is the source term in the two-dimensional shallow water equation;

[0100] The one-dimensional hydrodynamic model is used to calculate the river channel flow, and the governing equations adopt the Saint-Venant equations, and its vector form is:

[0101]

[0102] In the formula, U2 is the unknown quantity to be solved, F is the flux, and S2 is the source term. The finite volume method is used to discretize the governing equations with second-order accuracy in space, and the time discretization adopts the second-order Runge-Kutta explicit format. The numerical flux between grids is calculated by an approximate Riemann solver. The bottom slope term in the source term is discretized based on the hydrostatic reconstruction method and central difference. The friction term is treated fully implicitly and calculated using Newton-Raphson iteration to obtain better numerical stability. The model calculation step size adopts a conditional adaptive time step size to ensure the stability of the calculation.

[0103] In step S4, considering the lake water volume change, lake surface precipitation, lake surface evaporation, lake water surface area, total inflow of inflowing rivers, and the flow rate at the lake outlet, a lake water balance model regarding the lake water level in the previous simulation period and the lake water level in the current period is established:

[0104] ΔV = (P - E) * A + R in +R out

[0105]

[0106] where ΔV is the lake water volume change, P1 is the lake surface precipitation, E is the lake surface evaporation, A2 is the lake water surface area, R in is the total inflow of inflowing rivers, and R out is the flow rate at the lake outlet. h1 and H2 are the lake water levels in the previous simulation period and the current period respectively. The relationship between the lake surface area and the water level can be obtained through the water level - area curve; the total inflow of inflowing rivers includes the sum of the runoff of all inflowing rivers;

[0107] R out can be calculated by the following formula. Based on the assumption of broad - crested weir flow and assuming that the velocity head can be ignored, the discharge R out is calculated as:

[0108]

[0109] where b is the flow width (m), g is the acceleration due to gravity, z is the current lake depth (m), and z min is the elevation (m) above the weir or the lake outlet. The discharge coefficient c d is used to consider the inflow velocity, non - parallel streamlines at the top, and energy losses. The value of c d varies approximately between 0.8 and 1.2.

[0110] In step S5, the seed flooding algorithm is as follows: The seed flooding algorithm starts from the seed point and marks adjacent pixel regions as the same type or color by gradually spreading and filling pixels. The seed flooding algorithm is used in combination with the DEM to simulate the active flooding process of a shallow lake. The steps of the lake flooding process simulation based on the seed flooding algorithm include: determining the location and flow rate data of the inflowing river channels, initialization, starting the seed flooding, and updating the flooding process;

[0111] Determining the location and flow rate data of the inflowing river channels: Determine the injection point location and the corresponding flow rate data of the inflowing river. The injection point is regarded as the seed of the lake basin, and the flooding process is simulated starting from the injection point;

[0112] Initialization: Take the injection point as the seed point, and take the flow rate at the seed point as the initial inundation flow rate. At the same time, load the lake basin DEM data and the initial water level data into the calculation environment;

[0113] Seed inundation start: Starting from the seed point, search for the eight-neighborhood connected domain and diffuse to the surrounding 8 neighborhood grids. Layer by layer, judge whether each grid meets the inundation conditions, and add the grids that meet the conditions to the queue to be processed; Continuously iterate this process until no new grids are inundated;

[0114] Update the inundation process: According to the inundation elevation values of each time period, a dynamic image of the inundation process can be generated to show the process of the gradually expanding inundation range, thereby updating the inundation value H new :

[0115]

[0116] The H new is the updated inundation elevation value, H old is the inundation elevation value of the current grid, Q3 is the injection flow rate data, and A3 is the injection flow rate influence range;

[0117] Divide the unit grids into two categories: lake grids and non-lake grids, where the lake grids have seasonal changes in the inundation range during the simulation period;

[0118] Use the energy balance equation to calculate the surface energy components of each grid and analyze the net radiation R n :

[0119] R n = H + ρ w * λ v * E0 + G;

[0120] Where R n is the net radiation (W * m -2 ), H is the sensible heat flux (W * m -2 ), ρ w * λ v * E is the latent heat flux (W * m -2 )(ρ w is the density of liquid water, with the unit of kg * m -3 ; λ v is the latent heat of vaporization of water, with the unit of J * kg -1 ), and G is the geothermal flux (W * m -2 );

[0121] Based on the surface albedo, downward shortwave radiation, surface emissivity, and downward longwave radiation, obtain the net radiation R input to the grid n :

[0122]

[0123] where α is the surface albedo of the land cover type, R s is the downward shortwave radiation (W*m -2 ), ε is the surface emissivity of the land cover type, R L is the downward longwave radiation (W*m -2 ), σ is the Stefan-Boltzmann constant (5.67×10-8W*m -2 *K -4 );

[0124] Based on the aerodynamic resistance of heat flux, air density, specific heat capacity of air, surface temperature and surface air temperature, the sensible heat flux H is calculated as follows:

[0125]

[0126] where α is the grid inundation area ratio, T W is the water surface temperature, T g is the soil surface temperature of the exposed part of the grid, T a is the air temperature, ρ a is the air density, c p is the specific heat capacity of air, r h,w is the aerodynamic resistance of the water surface, r h,g is the aerodynamic resistance of the land surface, r s,w is the surface resistance of the water body, r s,g is the soil evaporation impedance.

[0127] The continuous transition of water-land energy flux is realized through the inundation ratio, accurately depicting the heat exchange behavior of "partially inundated" grids; the water temperature calculation includes a heat storage term, reflecting the heat buffering effect during flood retention; the introduced term includes the influence of the pressure gradient force, avoiding the overestimation of the latent heat flux in a high-humidity environment; the land surface resistance term is related to soil moisture, and the model's feedback on evaporation inhibition during the dry period is more significant.

[0128] where r h is the aerodynamic resistance of heat flux (s*m -1 ), ρ a is the air density, c p is the specific heat capacity of air, T s is the surface temperature, T a is the surface air temperature;

[0129] Based on the soil thermal conductivity, the soil temperature between the first and second soil layers, and the thickness of the first soil layer, the geothermal flux G of the topsoil is calculated as follows:

[0130]

[0131] where k is the soil thermal conductivity (W*m -1 *K -1 ), T1 is the soil temperature (K) between the first and second soil layers, and D1 is the thickness (m) of the first soil layer;

[0132] Based on the sensible heat flux, latent heat flux, changes in the heat storage in the overlying water body and soil of the floodplain, the advective heat flux carried by water flow, and the heat flux entering the underlying soil, the net surface radiation R n is calculated as:

[0133] R n = H + LE + ΔS2 + A4 + Q B ;

[0134] where R n is the net surface radiation, H is the sensible heat flux, LE is the latent heat flux, ΔS2 is the change in the heat storage in the overlying water body and soil of the floodplain, A4 is the advective heat flux carried by water flow, and Q B is the heat flux entering the underlying soil. Among them, H and LE are obtained from flux observation data, and ΔS2 and A4 need to be estimated based on key parameters such as specific heat capacity and flow velocity in combination with the energy equation;

[0135] Based on the soil depth, soil temperature, soil specific heat capacity, and the heat flux at the reference depth z ref , by integrating the one-dimensional heat diffusion equation, the change in the heat storage of the exposed soil ΔS(z) can be obtained:

[0136]

[0137] where z is the soil depth (m), T is the soil temperature (K), ρ s c s is the soil specific heat capacity (J*kg -1 *K -1 ), and S(z ref ) is the heat flux at the reference depth z ref . If z ref is greater than 1 m and S(z ref ) is less than 1% of the surface soil heat flux, it can be assumed that S(z ref ) = 0;

[0138] Given the temperature profile T(z i ), based on the time interval of the discrete time step, the reference depth, the interval in the vertical direction, and the depth coordinates in the vertical direction, the soil heat flux G is calculated as:

[0139]

[0140] Among them, Δt represents the time interval for discrete time steps; t represents the current moment; z ref represents the reference depth, Δz represents the interval in the vertical direction, and z i represents the depth coordinate in the vertical direction. The specific heat capacity of the soil can be calculated based on the soil water content and soil porosity, and the temperature profile is estimated using the soil heat diffusion equation;

[0141] Based on the total water depth, average water temperature, and the change in the average water temperature, the change in water body heat storage ΔS during the flooding period is calculated w :

[0142]

[0143] Among them, ρ w c p is the specific heat capacity of water (J*kg -1 *K -1 ), z w is the total water depth, T w is the average water temperature, and ΔT w is the change in the average water temperature within the time period Δt. According to water temperature observations, when the water depth is less than 1 m, it is considered that the lake water temperature is evenly mixed; when the water depth is greater than 1 m, ΔT w profile is estimated through the water thermal conductivity and the surface temperature from remote sensing. The change in heat storage is estimated by the weighted average of the water heat storage change (ΔS w ) during the flood season and the soil heat flux (G) during the dry season; due to the drastic change in water conditions and the relatively fast surface flow velocity of the water body in the floodplain wetland, the horizontal advection flux of the overlying water body is calculated by the following formula:

[0144]

[0145] Among them, u is the water flow velocity along the measurement direction of the temperature gradient, x is the horizontal distance, is the horizontal gradient of the water temperature, and ρ w c p is the specific heat capacity of water.

[0146] Example 1: In the basin, there are slopes and river channels formed by undulating mountains, as well as large areas of lake wetlands, which are of great significance to regional ecology and water resource regulation.

[0147] First, the basin space is discretized into two parts: the hillslope-channel and the lake-wetland. For the hillslope-channel unit, the focus is on simulating the water movement process, such as how the water flows on the hillslope after precipitation and the convergence and transmission of water in the channel. For the lake-wetland unit, while simulating the hydrological process, such as the inflow and outflow of lake water and the water level change, the surface energy exchange process also needs to be simulated because the energy exchange in the lake-wetland has a key impact on its ecosystem. It is assumed that the hillslope and channel unit and the lake-wetland unit exchange water only through the inlet or outlet channels of the lake, and overland flow across the lake boundary is prohibited.

[0148] Next, simulation grid units are divided for the hillslope, channel, and lake-wetland parts respectively. For the hillslope and lake-wetland parts, they are divided into square grid units according to the established spatial resolution, which can more accurately simulate their hydrological and energy processes. The channel part is divided into rectangular grids according to the channel shape to adapt to the morphological characteristics of the channel. Moreover, the hillslope and channel are divided into coarse grid units, and the lake-wetland part is divided into fine grid units. The resolutions of the coarse grid units and the fine grid units are in an integer multiple relationship, so as to effectively simulate at different scales.

[0149] In the process of runoff generation and concentration in the sub-basin, the runoff generation process at the hillslope scale includes precipitation, canopy interception, vegetation transpiration, soil evaporation, surface water infiltration, vertical movement of soil moisture, and soil water-groundwater exchange, etc. The runoff concentration process at the hillslope-channel scale covers the overland flow and subsurface flow sub-modules at the hillslope scale, and the channel runoff concentration sub-module at the channel scale, and finally the lake inflow runoff is accurately calculated.

[0150] Considering factors such as lake water volume change, lake surface precipitation, lake surface evaporation, lake water surface area, total inflow of inflowing rivers, and flow at the lake outlet, a lake water balance model is established to accurately grasp the dynamic change of the lake water level. Using the seed flooding algorithm, combined with DEM to simulate the active flooding process of shallow lakes, vividly demonstrating the dynamic process of the gradually expanding flooding range. At the same time, the distributed heat flux process is simulated, and through methods such as the energy balance equation, key energy components such as net radiation, sensible heat flux, and geothermal flux are comprehensively analyzed, providing strong support for in-depth understanding of the basin hydrology and energy conditions and facilitating the scientific management and protection of the basin.

[0151] For those skilled in the art, it is obvious that the present invention is not limited to the details of the above exemplary embodiments, and can be implemented in other specific forms without departing from the spirit or basic characteristics of the present invention. Therefore, in any regard, the embodiments should be regarded as exemplary and non-limiting. The scope of the present invention is defined by the appended claims rather than the above description. Therefore, all changes falling within the meaning and scope of the equivalent elements of the claims are intended to be included in the present invention. Any reference signs in the claims should not be regarded as limiting the claimed rights.

Claims

1. A method for simulating heat flux of flooded lake wetlands in a watershed hydrological model, characterized by: The method comprises the following steps: S1: Discretize the watershed space into two parts: slope river channel and lake wetland; S2: Divide the grid cells for the slope river channel and lake wetland respectively; S3: Establish evapotranspiration, infiltration and surface runoff models, and analyze the main hydrological processes of evaporation, infiltration, runoff generation and confluence based on the models; S4: Establish a lake water balance model to calculate the lake water level and the flow at the lake outlet; S5: Simulate the lake flooding process and simulate the distributed heat flux process in lake wetlands.

2. The method for simulating heat flux of flooded lake wetlands in a watershed hydrological model according to claim 1, characterized in that: In step S1, the watershed space is discretized into two parts: slope river channel and lake wetland. For the slope river channel unit, the water movement process is simulated; for the lake wetland unit, the surface energy exchange process is simulated on the basis of simulating the hydrological process. The slope river channel unit and the lake wetland unit exchange water through the river entering or leaving the lake, and no overflow across the lake boundary will occur.

3. The method for simulating heat flux of flooded lake wetlands in a watershed hydrological model according to claim 2, characterized in that: In step S2, the three parts of the slope, river channel and lake wetland are divided into grid units respectively. For the slope and lake wetland parts, square grid units are divided according to the established spatial resolution; for the river channel part, it is divided into rectangular grid units according to the shape of the river channel; the slope and river channel are divided into coarse grid units, and the lake wetland part is divided into fine grid units, and the resolution of the coarse grid units and the fine grid units are in integer multiples.

4. The method for simulating heat flux of flooded lake wetlands in a watershed hydrological model according to claim 3, characterized in that: In step S3, the runoff generation and confluence process of the sub-basin includes the runoff generation process at the slope scale and the confluence process at the slope river scale. The runoff generation process includes precipitation, canopy interception, vegetation transpiration, soil evaporation, surface water infiltration, vertical movement of soil moisture and soil water-groundwater exchange. The confluence process includes the overland flow and soil midflow submodules at the slope scale, and the river confluence submodule at the river scale, and finally the runoff into the lake is obtained. The slope is directly connected to the river, and the slope directly provides lateral inflow for the river confluence through the overland flow and soil midflow; based on net radiation, saturated water vapor pressure, actual water vapor pressure, actual water vapor pressure and wind speed, an evapotranspiration process model about evapotranspiration ET is established; Based on the saturated hydraulic conductivity of soil, the hydraulic head of soil, the pressure head of water in soil, the water potential of soil, and the infiltration surface area, an infiltration process model for the total infiltration volume I(t) is established. Based on precipitation, evapotranspiration, infiltration, and soil water and groundwater exchange, a surface runoff model for surface runoff Q1 is established; The two-dimensional shallow water equation is used to handle the slope runoff calculation, and the one-dimensional hydrodynamic model is used to calculate the river runoff.

5. The method for simulating heat flux of flooded lake wetlands in a watershed hydrological model according to claim 4, characterized in that: In step S4, a lake water balance model is established based on the lake water level in the previous simulation period and the lake water level in the current period, taking into account changes in lake water volume, lake surface precipitation, lake surface evaporation, lake surface area, total flow of rivers entering the lake and flow at the lake outlet.

6. The method for simulating heat flux of flooded lake wetlands in a watershed hydrological model according to claim 5, characterized in that: In step S5, the lake flooding process is simulated based on the seed flooding algorithm to simulate the distributed heat flux process of the lake wetland. The seed flooding algorithm is as follows: the seed flooding algorithm takes the seed point as the starting point, and marks the adjacent pixel areas as the same type or color by gradually diffusing and filling the pixels. The seed flooding algorithm cooperates with the DEM to simulate the active flooding process of the shallow lake. The steps of simulating the lake flooding process based on the seed flooding algorithm include: determining the location and flow data of the river entering the lake, initialization, starting the seed flooding and updating the flooding process; Determine the location and flow data of the river entering the lake: Determine the location of the injection point of the river entering the lake and the corresponding flow data. The injection point is regarded as the seed of the lake basin, and the flooding simulation process starts from the injection point; Initialization: The injection point is used as the seed point, and the flow at the seed point is used as the initial flooding flow. At the same time, the lake basin DEM data and initial water level data are loaded into the computing environment; Seed flooding starts: starting from the seed point, the eight-neighborhood connected domain search spreads to the surrounding eight neighboring grids, judging layer by layer whether each grid meets the flooding condition, and adding the grid that meets the condition to the queue to be processed; this process is iterated continuously until no new grid is flooded; Update the flooding process: According to the flooding elevation values ​​at each time period, a dynamic image of the flooding process can be generated to show the gradual expansion of the flooding range, thereby updating the flooding value H new .

7. The method for simulating heat flux of flooded lake wetlands in a watershed hydrological model according to claim 6, characterized in that: The cell grids are divided into two categories: lake grids and non-lake grids, where the lake grids have seasonal inundation extent changes during the simulation period.

8. The method for simulating heat flux of flooded lake wetlands in a watershed hydrological model according to claim 7, characterized in that: The energy balance equation is used to calculate the surface energy components of each grid and analyze the net radiation R n ; Based on the surface albedo, downward shortwave radiation, surface emissivity and downward longwave radiation, the net radiation R input to the grid is obtained. n ; The sensible heat flux H is calculated based on the aerodynamic resistance to heat flow, air density, air specific heat capacity, ground surface temperature, and ground-surface air temperature: Where α is the grid flooded area ratio, T W is the water surface temperature, T g is the soil surface temperature of the exposed part of the grid, T a is the air temperature, ρ a is the air density, c p is the specific heat capacity of air, r h,w is the aerodynamic resistance on the water surface, r h,g is the aerodynamic resistance on the land surface, r s,w is the water surface resistance, r s,g is the soil evaporation resistance; Based on the soil thermal conductivity, the soil temperature between the first and second soil layers and the thickness of the first soil layer, the geothermal flux G of the top soil layer is calculated: where k is the thermal conductivity of the soil, T1 is the soil temperature between the first and second soil layers, and D1 is the thickness of the first soil layer.

9. The method for simulating heat flux of flooded lake wetlands in a watershed hydrological model according to claim 7, characterized in that: The net surface radiation R is calculated based on the sensible heat flux, latent heat flux, changes in heat storage in the overlying water and soil of the floodplain, the advective heat flux carried by the water flow, and the heat flux entering the underlying soil. n ; Based on soil depth, soil temperature, soil specific heat capacity and reference depth z ref The heat flux at the point can be integrated into the one-dimensional heat diffusion equation to obtain the change in heat storage of the soil during the exposure period, ΔS(z): Given a temperature profile T(z i ), based on the time interval of discrete time steps, reference depth, vertical interval, and vertical depth coordinate, the soil heat flux G is calculated: Where Δt represents the time interval, which is used for discrete time steps; t represents the current moment; z ref represents the reference depth, Δz represents the interval in the vertical direction, z i It represents the depth coordinate in the vertical direction. The soil specific heat capacity can be calculated based on the soil water content and soil porosity. The temperature profile is estimated using the soil heat diffusion equation. Based on the total water depth, average water temperature, and the change in average water temperature, the change in water heat storage during the flooding period, ΔSw, was calculated; the horizontal advection flux A of the overlying water body was calculated by the following formula: Where u is the water flow velocity along the temperature gradient measurement direction, x is the horizontal distance, is the horizontal gradient of water temperature, ρ w c p is the specific heat capacity of water.

Images