Modularized distributed hydrothermal coupling hydrological model in cold region

Abstract

The application discloses a modularized distributed water and heat coupling hydrological model in cold regions, and comprises the following steps: S1, a basin is divided into a plurality of grids, and files required by a confluence process are prepared; S2, driving data and grid characteristic description data of each grid are prepared; S3, an energy balance module is used to calculate an energy balance process of soil and a vegetation surface; S4, an evapotranspiration module is used to calculate an evapotranspiration process; S5, an ice and snow melt water module is used, and a degree day factor method is used to calculate an ice and snow melting process; S6, a soil water and heat transmission and infiltration module is used to calculate soil water infiltration and internal water and heat transmission processes; S7, a runoff calculation and confluence module is used to calculate runoff; and S8, a model evaluation index is used to evaluate the simulation performance of the model. The application has a flexible modular structure, can automatically adjust a time step, improves the simulation capability of a runoff process in cold regions, and can provide scientific guidance for water resource management in a river source region.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of frozen soil hydrology and basin runoff generation, and particularly relates to a modularized distributed water and heat coupled hydrology model in cold regions. BACKGROUND

[0002] As a special regional "water-resisting layer" or "weak water-permeable layer", permafrost significantly reduces or hinders the hydraulic connection between precipitation, surface water body and groundwater at a certain spatiotemporal scale, and strongly affects the runoff generation process of surface runoff, the migration and distribution pattern of groundwater runoff. Under the background of climate change, the degradation of permafrost will have a profound impact on the regional hydrological conditions, including changes in soil moisture, seasonal characteristics of runoff, and distribution changes of surface and groundwater storage. Therefore, accurately simulating and predicting the changes of permafrost and its impact on hydrological processes under the background of climate change is of great importance to the management of water resources downstream.

[0003] Traditional hydrological models use some empirical formulas (such as Stefan equation, air temperature-water permeability index relationship, etc.) to reflect the influence of frozen soil on hydrological processes. Although this simple empirical method of linking air temperature and soil freezing and thawing cycle improves the accuracy of frozen soil basin runoff simulation to a certain extent, it cannot represent the mutual feedback of soil water and heat transfer processes, and is not suitable for studying the changes of frozen soil and hydrological processes under changing environment and the interaction between permafrost degradation and hydrological processes. Process-based hydrological models in cold regions, such as VIC, GBEHM, WEB-DHM-pf, etc., often have high requirements for input data, and have certain limitations in the application of data-poor Qinghai-Tibet Plateau region. In addition, most process-based models are designed with a specific problem-oriented or goal-oriented method, relying on a single mathematical theoretical basis, a unique parameterization scheme and fixed data requirements. This greatly limits the applicability of the model and its simulation / prediction ability under actual modeling conditions. SUMMARY

[0004] In order to solve the above problems, the present application provides a modularized distributed water and heat coupled hydrology model in cold regions to solve the above problems.

[0005] The present application discloses a modularized distributed water and heat coupled hydrology model in cold regions, which comprises an energy balance module, an evapotranspiration module, a snowmelt water module, a soil water and heat transport and infiltration module and a runoff calculation and confluence module, comprising the following steps:

[0006] S1, according to the elevation data and the basin outlet position, the GIS software is used to extract the basin range, and the basin is divided into several grids of different sizes, and the flow direction data, grid spherical distance and grid area proportion (the proportion of the grid inside the basin is 1, and the proportion of the grid on the boundary of the basin is less than 1) and other files required for the confluence process are prepared;

[0007] S2, the driving data (air temperature, precipitation, wind speed, relative humidity, incident shortwave radiation, incident longwave radiation, leaf area index, etc.) and the grid characteristic description data (elevation, soil layer thickness, soil texture, organic matter content, etc.) of each grid are prepared;

[0008] S3, the energy balance module is used to calculate the energy balance process of the soil and vegetation surface;

[0009] S4, the evapotranspiration module is used to calculate the evapotranspiration process;

[0010] S5, the ice and snow melt water module is used to calculate the ice and snow melting process by using the degree day factor method;

[0011] S6, the soil water and heat transport and infiltration module is used to calculate the soil water infiltration and internal water and heat transport process;

[0012] S7, the runoff calculation and confluence module is used to calculate the runoff;

[0013] S8, the simulation performance of the model is evaluated by the model evaluation index.

[0014] Preferably, the S3 comprises the following steps:

[0015] The net shortwave radiation is calculated as follows:

[0016] Ns veg =R s *FVC*[(1-α c )+τ c *(α soil / snow -1)];

[0017] Ns soil / snow =R s *(1-α soil / snow )*[(1-FVC)+τ c *FVC];

[0018] Wherein, Ns veg is the net shortwave radiation absorbed by the vegetation surface, Ns soil / snow is the net shortwave radiation absorbed by the bare soil / snow surface, R s is the incident shortwave radiation, τ c is the shortwave radiation transmission ratio of vegetation, α c is the albedo of the vegetation surface, and α soil / snowAlbedo of bare soil / snow surface, FVC is vegetation coverage;

[0019] Calculate net longwave radiation:

[0020] Nl veg = FVC * [L d + L soil / snow - 2 * L c ];

[0021] Nl soil / snow = (1 - FVC) * L d + FVC * L c - L soil / snow ];

[0022] Wherein, Nl veg is the net longwave radiation absorbed by the vegetation surface, Nl soil / snow is the net longwave radiation absorbed by the bare soil / snow surface, L c is the longwave radiation emitted by the canopy surface, L d is the incident longwave radiation, L soil / snow is the longwave radiation emitted by the bare soil / snow surface;

[0023] The net radiation of the surface is equal to the sum of the net shortwave radiation and the net longwave radiation:

[0024] R net = Ns + Nl;

[0025] Wherein, Ns is the net shortwave radiation, Nl is the net longwave radiation;

[0026] The formula for calculating the ground heat flux is as follows:

[0027] G0 = 0.35462 * R net - 47.79008;

[0028] Wherein, G0 is the ground heat flux.

[0029] Preferably, the S4 comprises the following steps:

[0030] The total evapotranspiration (E total ) of each grid is composed of three parts: vegetation transpiration (E trans ), interception evaporation (E int ) and soil evaporation (E soil ). Based on the calculated energy balance of the vegetation surface and the bare soil surface, transpiration, interception evaporation and soil evaporation are calculated respectively. The proportion of vegetation roots at different depths is calculated, and the calculated transpiration is distributed to each layer of soil according to the proportion of roots, and soil evaporation is completely distributed to the surface soil layer. The specific calculation process of each evaporation component is as follows:

[0031] Calculate potential evapotranspiration:

[0032]

[0033] where Δ is the slope of the curve of saturated water vapor pressure as a function of temperature, Q ne is the available energy, including latent heat consumed by ice-water phase change, p air is the density of air, c air is the specific heat of air, e s is the saturated water vapor pressure, e a is the actual water vapor pressure, r a is the aerodynamic impedance, r c is the canopy impedance, γ is the dry- wet constant, p liq is the density of liquid water, L v is the latent heat of vaporization;

[0034] The interception evaporation of the canopy surface is calculated as:

[0035]

[0036] where Int sto is the intercepted water output, Max it is the maximum interception water storage capacity, Δt is the time step;

[0037] The transpiration of the canopy surface is calculated as:

[0038]

[0039] E trans = E pot * f dry * root frac,i * Δt;

[0040] where f dry is the dryness ratio of the canopy surface, E trans is the transpiration of the canopy surface, root frac,i represents the root ratio in the i-th soil layer;

[0041] The soil evaporation is calculated as:

[0042]

[0043] where θ1is the liquid water content of the surface soil, η1is the porosity of the surface soil, K sat is the saturated hydraulic conductivity of the soil, m is the inverse of the Campbell pore size distribution exponent, Ψ e is the air entry potential.

[0044] Preferably, the S5 comprises the following steps:

[0045] The melting process of the accumulated snow is calculated by the following formula:

[0046]

[0047] wherein M snow is the snowmelt, DDF snow is the snowmelt degree-day factor, SRF snow is the shortwave radiation factor, Ns snow is the net shortwave radiation of the snowmelt surface, T air is the air temperature, T snow represents the highest temperature threshold when the precipitation is snowfall;

[0048] The glacier zone method is combined with the degree-day factor method and the area-volume-ratio method to simulate the glacier runoff process, including the following steps:

[0049] Each glacier is divided into several glacier zones, and the air temperature and precipitation data of each glacier zone are interpolated, and the area proportion of each glacier zone in each grid is calculated;

[0050] Rain and snow are divided based on the air temperature and precipitation data of each glacier zone obtained by interpolation, the rainfall and snowfall of each glacier zone are obtained, and the snow water equivalent of each glacier zone is calculated;

[0051] When the air temperature falls below the snowmelt threshold temperature, the potential snowmelt is calculated by the degree-day factor method, and compared with the current snow water equivalent of the glacier zone. If the potential snowmelt is less than the snow water equivalent, only snowmelt process occurs, otherwise the underlying glacier will start to melt:

[0052]

[0053] wherein M band represents the total snow and glacier meltwater, DDF ice is the glacier melting degree-day factor, SWE band is the current snow water equivalent of the glacier zone, M snow,band is the potential snowmelt, T band is the air temperature of the glacier zone, T melt is the melting temperature of the snow or glacier;

[0054] The runoff depth on each glacier zone is equal to the sum of rainfall and total snow and ice meltwater. On the grid scale, the runoff depth of the glacier area is equal to the product of the runoff yield and the area proportion of all glacier zones contained in the grid:

[0055] R band,i = P rain,band,i + M band,i ;

[0056]

[0057] where i is the number of the glacier zone, n is the number of glacier zones, R band,i is the runoff depth of the ith glacier zone, R glacier is the runoff depth of the glacier region, P rain,band,i is the rainfall of the ith glacier zone, M band,i is the snowmelt water of the ith glacier zone, A frac,i is the area ratio of the ith glacier zone in the grid.

[0058] Preferably, S6 uses a multi-layer soil structure to represent the water and heat transfer process inside the soil, and the division of soil thickness is determined by exponential layering method; within each time step, the calculation process of the water and heat transfer process inside the soil layer is as follows:

[0059] S61, based on the current soil liquid water content, ice content, and organic matter content, calculate the water and heat parameters of the soil, including thermal conductivity, water conductivity, volumetric heat capacity, and water flux;

[0060] S62, solve the heat transfer equation to obtain the temperature of each layer of soil;

[0061] S63, update the soil liquid water content and ice content, and the water conductivity and water flux based on the latest calculated soil temperature;

[0062] S64, solve the Richard equation considering ice-water phase change to obtain the water content and ice content of each layer of soil;

[0063] S65, repeat steps S61-S64 until the soil temperature and soil moisture error reaches the convergence criteria;

[0064] S66, use the multi-layer Green-Ampt algorithm to calculate the water infiltration process, and update the soil moisture and soil temperature.

[0065] Preferably, the water and heat transfer process inside the soil layer is calculated by the following formula:

[0066] For each layer of soil, calculate the change of soil temperature:

[0067]

[0068] where C s is the volumetric heat capacity of the soil, λ s is the thermal conductivity of the soil, ρ ice is the density of ice, θ ice is the volumetric ice content, q l is the liquid water flux, C liq is the volumetric heat capacity of liquid water, T is the temperature of the soil, z is the depth of the soil, L f is the latent heat of fusion;

[0069] The volumetric heat capacity of soil is expressed as the weighted sum of the volumetric heat capacities of all components in the soil:

[0070] C s =(1-organic)*ρ b *c m +organic*ρ b *c o +θ liq *ρ liq *c liq +θ ice *ρ ice *c ice +θ air *ρ air *c air ;

[0071] Where, ρ b For soil bulk density, c m For the specific heat of minerals, c o For the specific heat of organic matter, θ air θ represents the volume fraction of air. liq Liquid water content (m 3 ·m -3 ), where organic is the volume fraction of organic matter;

[0072] The thermal conductivity of the soil was obtained by weighted combination of dry-state thermal conductivity and saturated-state thermal conductivity using the Kosten number:

[0073] λ s =(λ sat -λ dry )*K e +λ dry ;

[0074] Among them, K e λ is a function of soil saturation and the volume fraction of organic matter. dry λ is the thermal conductivity of the soil in a dry state. sat Thermal conductivity of soil under saturated conditions;

[0075]

[0076] Among them, S r v represents soil saturation. sand v represents the volume fraction of sand grains. gravel The volume fraction of gravel is represented by α and β, which are adjustment factors; K is used to represent the volume fraction of gravel. e -S r The relationship between these parameters is used to predict the continuous variation of thermal conductivity across the entire soil saturation level and particle size distribution range, and is applicable to the calculation of thermal conductivity of permafrost in arid and semi-arid regions.

[0077] In cold regions, when the soil temperature is above freezing, the matric potential is expressed as a function of liquid water content, and when the temperature drops below freezing, the matric potential drops rapidly and is given by the freezing-point water potential equation:

[0078]

[0079] where Ψ is the soil matric potential, b is the Campbell pore-size distribution exponent, T f is the freezing point, θ sat is the saturated water content, and g is the gravitational acceleration;

[0080] The movement of water is described by the matric potential gradient:

[0081]

[0082] where K eff is the effective hydraulic conductivity;

[0083] When the soil begins to freeze, the retarding effect of soil ice on water transport is considered by introducing a retardation factor:

[0084]

[0085] where E is the retardation factor, m ice is the ratio of ice content to total water content;

[0086] Assuming that the soil medium is incompressible and ignoring water vapor migration, the basic equation of soil water movement is obtained according to Darcy's law and the principle of mass conservation:

[0087]

[0088] where S is the source-sink term in the water flux.

[0089] Preferably, the S7 comprises the following steps:

[0090] S71, calculate the surface runoff depth and the groundwater runoff depth of each grid.

[0091] Snowmelt and rain-through are the sources of infiltration water. The surface water storage is equal to the difference between the sum of rain-through and snowmelt and the total cumulative infiltration water, and the surface runoff is calculated by the following formula:

[0092] R sur = max(0, [(P rain + M snow )- Infil - Pond max ]*(1-glacier frac )+ glacier frac *Rglacier );

[0093] where P rain is the precipitation of the grid, Pond max is the maximum allowed ponding depth (mm), Infil is the total cumulative infiltration, glacier frac is the glacier area fraction of the grid.

[0094] During the soil freezing and thawing process, the freezing or thawing front temporarily acts as an impermeable layer, making it difficult to distinguish between the frozen layer water and interflow in the soil. Therefore, the subsurface runoff is used to describe the subsurface runoff process. The subsurface runoff of the ith layer of soil is:

[0095]

[0096] where θ f is the field capacity, Δz i is the thickness of the ith layer of soil, and α r is the subsurface runoff factor (-).

[0097] The total subsurface runoff is the sum of R sub,i of all soil layers above the bedrock.

[0098] S72, Confluence process calculation:

[0099] Based on the elevation data, the confluence file is prepared to simulate the output surface runoff depth and subsurface runoff depth as the driving force, and the Lohmann confluence model is used to calculate the slope confluence and river network confluence process to obtain the outlet runoff of the watershed. In the cold season, the river flow of permafrost watershed is mainly maintained by the water below the frozen layer, which is represented by Q90 in the flow duration curve. Based on the runoff variation characteristics of several rivers in the Siberian and Qinghai-Tibet permafrost regions in recent decades, an empirical relationship between annual average temperature, annual average precipitation, permafrost coverage and water below the frozen layer is established:

[0100]

[0101] where δ, β1, β2 and β3 are fitting parameters, is the multi-year moving average of air temperature (℃), is the multi-year moving average of precipitation (mm), The permafrost coverage (%) is a multi-year moving average. The Lohmann confluence model is used to calculate the confluence process, and the daily surface runoff depth and the groundwater runoff depth are used to drive the Lohmann confluence model to calculate the slope and river network confluence process. The Lohmann confluence model converges the runoff to the outlet of a given grid, and then adds it to the river network confluence model that couples all the grids together. The river network confluence process is estimated by solving the linear Saint-Venant equation to obtain the runoff at the outlet of the watershed. The final runoff at the outlet of the watershed is the sum of the output of the Lohmann confluence model and the output of the frozen layer water module, that is:

[0102]

[0103] wherein, is the final runoff at the outlet of the watershed (m 3 / s), is the watershed outlet runoff output by the Lohmann confluence model (m 3 / s), is the frozen layer water.

[0104] Preferably, the model evaluation index includes Nash efficiency coefficient, Kling-Gupta efficiency coefficient, percentage bias, root mean square error and determination coefficient.

[0105] The beneficial effects of the present application are:

[0106] (1) The model of the present application has a flexible modular structure, can be integrated with the most advanced methods or models, and can adjust the model structure and simulation strategy according to local characteristics and data availability.

[0107] (2) The model of the present application can automatically adjust the time step according to the convergence condition to ensure that the model always converges in the iterative calculation process of soil moisture and heat transport.

[0108] (3) The model of the present application comprehensively considers the influence of soil freeze-thaw cycle on infiltration, evapotranspiration and water-heat transport process in soil, greatly improving the simulation capability of runoff process in cold regions.

[0109] (4) The model of the present application can simulate and analyze the permafrost thermal conditions of the region, and is used to study the influence of permafrost degradation on hydrological processes under the background of climate change, and to provide scientific guidance for water resources management in the source region of rivers. BRIEF DESCRIPTION OF DRAWINGS

[0110] Figure 1 is a schematic diagram of the modular distributed cold region water-heat coupled hydrological model of the embodiment of the present application.

[0111] Figure 2A Yangtze River source area basin profile of an embodiment of the present application.

[0112] Figure 3 A Fenghuoshan small watershed multiple borehole different depth soil temperature simulation and measured comparison chart of an embodiment of the present application.

[0113] Figure 4 A Fenghuoshan small watershed multiple borehole different depth soil water content simulation and measured comparison chart of an embodiment of the present application.

[0114] Figure 5 A Yangtze River source area 2001-2017 evapotranspiration simulation and measured comparison chart of an embodiment of the present application.

[0115] Figure 6 A Yangtze River source area 2001-2017 snow depth simulation and measured comparison chart of an embodiment of the present application.

[0116] Figure 7 A Yangtze River source area 2001-2017 runoff simulation and measured comparison chart of an embodiment of the present application. DETAILED DESCRIPTION

[0117] In order to make the purpose, technical scheme and advantages of the present application more clear, the following refers to the drawings and embodiments, and further details of the present application are described.

[0118] The embodiment of the present application discloses a modular distributed water and heat coupled hydrological model in cold regions, including an energy balance module, an evapotranspiration module, a snow and ice melt water module, a soil water and heat transport and infiltration module, and a runoff calculation and confluence module, as shown in Figure 1 The following steps are included:

[0119] S1, according to the elevation data and the outlet position of the basin, the GIS software is used to extract the basin range, and the basin is divided into several grids of different sizes, and the files required for the confluence process such as flow direction data, grid spherical distance and grid area ratio (the grid ratio in the basin is 1, and the grid ratio on the boundary of the basin is less than 1) are prepared.

[0120] S2, prepare the driving data (air temperature, precipitation, wind speed, relative humidity, incident shortwave radiation, incident longwave radiation, leaf area index, etc.) and grid characteristic description data (elevation, soil layer thickness, soil texture, organic matter content, etc.) of each grid.

[0121] S3, the energy balance module is used to calculate the energy balance process of the soil and vegetation surface.

[0122] Net shortwave radiation:

[0123] Ns veg = R s * FVC * [(1-αc )+τ c *(α soil / snow -1)];

[0124] Ns soil / snow =R s *(1-α soil / snow )*[(1-FVC)+τ c *FVC];

[0125] where Ns veg is the net shortwave radiation absorbed by the vegetation surface (W·m -2 ), Ns soil / snow is the net shortwave radiation absorbed by the bare soil / snow surface (W·m -2 ), R s is the incident shortwave radiation (W·m -2 ), τ c is the proportion of shortwave radiation transmitted by the vegetation (-), α c is the albedo of the vegetation surface (-), α soil / snow is the albedo of the bare soil / snow surface (-), and FVC is the fractional vegetation cover (-).

[0126] Net longwave radiation:

[0127] Nl veg =FVC*[L d +L soil / snow -2*L c ];

[0128] Nl soil / snow =(1-FVC)*L d +FVC*L c -L soil / snow ];

[0129] where Nl veg is the net longwave radiation absorbed by the vegetation surface (W·m -2 ), Nl soil / snow is the net longwave radiation absorbed by the bare soil / snow surface (W·m -2 ), L c is the longwave radiation emitted by the canopy surface (W·m -2 ), L d is the incident longwave radiation (W·m -2 ), and L soil / snow is the longwave radiation emitted by the bare soil / snow surface (W·m -2 ).

[0130] The net radiation of a surface is the sum of the net shortwave and longwave radiations:

[0131] R net=Ns+Nl;

[0132] wherein Ns represents net shortwave radiation (W·m -2 ), and Nl represents net longwave radiation (W·m -2 ).

[0133] The ground surface heat flux is calculated according to the following formula:

[0134] G0=0.35462*R net -47.79008;

[0135] S4, calculating the evapotranspiration process by using the evapotranspiration module.

[0136] The total evapotranspiration (E total ) of each grid is composed of three parts, i.e., vegetation transpiration (E trans ), interception evaporation (E int ) and soil evaporation (E soil ). The transpiration, interception evaporation and soil evaporation are calculated based on the calculated energy balance of the vegetation surface and the bare soil surface. The calculated transpiration is distributed to each layer of soil according to the proportion of the vegetation root at different depths, and the soil evaporation is completely distributed to the surface soil layer. The specific calculation process of each evaporation component is as follows:

[0137] Calculating the potential evapotranspiration:

[0138]

[0139] wherein Δ is the slope of the curve of saturated water vapor pressure with temperature (Pa / K), Q ne is the available energy (W / m 2 ), including the latent heat consumed by ice-water phase change, ρ air is the density of air (kg / m 3 ), c air is the specific heat of air (J / (kg·K)), e s is the saturated water vapor pressure (Pa), e a is the actual water vapor pressure (Pa), r a is the aerodynamic resistance (s / m), r c is the crown layer resistance (s / m), γ is the dry and wet surface constant (Pa·K -1 ), ρ liq is the density of liquid water (kg / m 3 ), and L v is the latent heat of vaporization (MJ / kg).

[0140] Calculating the interception evaporation of the crown surface:

[0141]

[0142] where Int sto is the interception water yield (m), Max it is the maximum interception water storage capacity (m), and Δt is the time step (s). Here, E pot is calculated with r c = 0.

[0143] The transpiration of the canopy surface is calculated as:

[0144]

[0145] E trans = E pot * f dry * root frac,i * Δt;

[0146] where f dry is the dryness ratio of the canopy surface, E trans is the transpiration of the canopy surface (m), and root frac,i represents the root ratio in the i-th soil layer.

[0147] The soil evaporation is calculated as:

[0148]

[0149] where θ1 is the liquid water content of the surface soil, η1 is the porosity of the surface soil, K sat is the saturated hydraulic conductivity of the soil, m is the inverse of the Campbell pore-size distribution exponent, and Ψ e is the air entry potential (m). Here, E pot is calculated with r c = 0.

[0150] S5. The ice and snow melting module is used to calculate the process of ice and snow melting using the degree-day factor method.

[0151] The melting process of ice and snow is calculated using the degree-day factor method, in which the melting process of snow is calculated by the following formula:

[0152]

[0153] where M snow is the snowmelt, DDF snow is the snowmelt degree-day factor (mm·℃ -1 ·day -1 ), SRF snow is the shortwave radiation factor (mm·W -1 ·m 2 ·day -1 ), and Ns snowNet shortwave radiation of snow melting surface (W·m) -2 ), T air T represents the air temperature (°C). snow This indicates the highest temperature threshold (°C) at which precipitation occurs, assuming snowfall.

[0154] The "glacier zone" method, combined with the diurnal factor method and the area-volume ratio method, is used to simulate glacier runoff processes, including the following steps:

[0155] First, each glacier was divided into several "glacier zones" at intervals of 100m elevation difference. The MicroMet model was then used to calculate the air temperature (T) of each "glacier zone." band (℃) and precipitation (P) band Interpolate the data (mm) and calculate the area ratio (A) of each "glacier zone" within each grid. frac,i ).

[0156] The air temperature (T) of each "glacier zone" obtained by interpolation band (℃) and precipitation (P) band The data (mm) were divided into "rain and snow" categories to obtain the precipitation (P) for each "glacier zone". rain,band (mm) and snowfall (P) snow,band (mm), and calculate the snow water equivalent of each "glacier zone".

[0157] When the temperature drops below the snowmelt threshold temperature, the potential snowmelt amount (M) is calculated using the degree-day factor method. snow,band This is compared with the current snowmelt equivalent of the "glacier zone." If the potential snowmelt is less than the snowmelt equivalent, only snowmelt will occur; otherwise, the underlying glacier will begin to melt.

[0158]

[0159] Among them, M band DDF represents the total snow cover and glacial meltwater (mm). ice SWE is the glacier melt degree-day factor (mm / (℃·day)). band M represents the current glacial meltwater equivalent (mm). snow,band T represents potential snowmelt (mm). band T represents the air temperature (°C) in the glacier zone. melt The melting temperature (°C) of snow or glaciers is taken as 0°C in this embodiment.

[0160] The runoff depth in each "glacier zone" is equal to the sum of rainfall and total snowmelt. On a grid scale, the runoff depth in a glacier region is equal to the sum of the products of the runoff and the area percentage of all "glacier zones" included in that grid.

[0161] R band,i= P rain,band,i + M band,i ;

[0162]

[0163] where i is the index of the "glacier zone", n is the number of "glacier zones", P rain,band,i is the precipitation of the i-th "glacier zone" (mm), R glacier is the runoff depth of the glacier region (mm), M band,i is the snowmelt of the i-th "glacier zone" (mm), A frac,i is the area ratio of the i-th "glacier zone" in the grid, R band,i is the runoff depth of the i-th "glacier zone" (mm).

[0164] The key of the above-mentioned "glacier zone" snow and ice ablation simulation algorithm lies in the simulation of glacier retreat or advance. At present, the glacier area evolution estimation method can be roughly divided into two categories. One is the ice-flow model based on the physical properties, mechanical processes and surrounding environmental conditions of ice to predict the movement and deformation of ice. However, this method requires detailed information of the glacier surface and base geometry, so it has limitations in application in areas with scarce data. The volume-area ratio method with lower parameter requirements is widely used to simulate the change of glacier area:

[0165]

[0166] where A glacier is the area of the glacier (km 2 ), V glacier is the volume of the glacier (km 3 ), and the constant c is taken as 0.0365 and the scale factor γ is taken as 1.375.

[0167] The volume of the glacier is updated every year according to the simulated glacier ablation, and the area of the glacier is calculated according to the updated glacier volume. The area of each "glacier zone" is updated according to the adjusted glacier area, and it is assumed that the retreat or advance of the glacier first occurs in the "glacier zone" with the lowest elevation. When the glacier retreats, the area of the "glacier zone" with the lowest elevation will decrease. If the calculated glacier shrinkage area is greater than the area of the lowest "glacier zone", the glacier area of the "glacier zone" will be set to 0, and the glacier area of the adjacent upper "glacier zone" will be reduced by the remaining shrinkage area. Conversely, when the glacier advances, the glacier area of the lowest "glacier zone" will increase and be limited within the initial input area, otherwise a new "glacier zone" will be introduced.

[0168] S6, calculate the soil water infiltration and internal water and heat transfer process by using the soil water and heat transfer and infiltration module.

[0169] A multi-layer soil structure is used to represent the water and heat transport process inside the soil. The division of soil thickness can be specified by the user or determined by the exponential layering method:

[0170] z i =f s {e (0.5(i-0.5)) -1};

[0171]

[0172] where z i is the depth of the i-th soil layer, f s is the scaling factor, which is taken as 0.025 in this embodiment, and Δz i is the thickness of the i-th soil layer. In the exponential layering method, the soil layer thickness is thinner closer to the surface due to the high sensitivity to meteorological driving. The basic calculation process of the water and heat transport process inside the soil layer in each time step is as follows:

[0173] S61, calculate the water and heat parameters of the soil based on the current soil liquid water content, ice content, and organic matter content, including thermal conductivity, hydraulic conductivity, volumetric heat capacity, and water flux.

[0174] S62, solve the heat transport equation to obtain the temperature of each layer of soil.

[0175] S63, update the soil liquid water content and ice content, as well as the hydraulic conductivity and water flux, based on the newly calculated soil temperature.

[0176] S64, solve the Richard equation considering ice-water phase change to obtain the water content and ice content of each layer of soil.

[0177] S65, repeat steps S61-S64 until the soil temperature and soil moisture error reaches the convergence criteria, in this embodiment, the convergence criteria is: ΔT < 0.001℃, Δθ < 0.01m 3 / m 3 , ΔT is the soil temperature error of the previous two iterations, and Δθ is the soil moisture error of the previous two iterations.

[0178] S66, calculate the water infiltration process using the multi-layer Green-Ampt algorithm, and update the soil moisture and soil temperature.

[0179] The water and heat transport process inside the soil layer is calculated by the following formula:

[0180] For each layer of soil, the heat conduction-convection equation is used to calculate the change of soil temperature, considering the ice-water phase change:

[0181]

[0182] Among them, C s The volumetric heat capacity of soil (J·kg) -1 ·℃ -1 ), λ s Thermal conductivity of soil (W·m) -1 ·℃ -1 ), ρ ice The density of ice (kg·m) -3 ), θ ice Ice content by volume (m³) 3 ·m -3 ), q l Liquid water flux (m·s) -1 ), C liq The volumetric heat capacity of liquid water (J·kg⁻¹) -1 ·℃ -1 T is the soil temperature (°C), z is the soil depth (m), and L is the soil depth. f To melt potential (3.35×10) 5 J·kg -1 ).

[0183] The volumetric heat capacity of soil is expressed as the weighted sum of the volumetric heat capacities of all components in the soil:

[0184] C s =(1-organic)*ρ b *c m +organic*ρ b *c o +θ liq *ρ liq *c liq +θ ice *ρ ice *c ice +θ air *ρ air *c air ;

[0185] Where, ρ b For soil bulk density, c m The specific heat of the mineral (900 J / (kg·℃)), a o The specific heat of organic matter (1920 J / (kg·℃)), θ air θ represents the volume fraction of air. liq Liquid water content (m 3 ·m -3 ), where organic is the volume fraction of organic matter.

[0186] The thermal conductivity of the soil is obtained by a Kersten number (Ke) weighted combination of the dry and saturated state thermal conductivities:

[0187] λ s = (λ sat - λ dry ) * K e + λ dry ;

[0188] where K e is a function of the volumetric fraction of soil saturation and organic matter, λ dry is the dry state thermal conductivity of the soil, and λ sat is the saturated state thermal conductivity of the soil.

[0189]

[0190] where S r is the soil saturation, v sand is the volumetric fraction of sand, v gravel is the volumetric fraction of gravel, and a (taken as 0.24 in this example) and β (taken as 18.1 in this example) are adjustment factors. The relationship between K e and S r is used to predict the continuous change in thermal conductivity across the range of soil saturation and particle size distributions, which is applicable to the calculation of the thermal conductivity of frozen soils in arid and semi-arid regions.

[0191] In cold regions, the matric potential is expressed as a function of the liquid water content when the soil temperature is above freezing. However, when the temperature drops below freezing, the matric potential (neglecting osmotic potential) rapidly decreases and is given by the freezing point water potential equation:

[0192]

[0193] where Ψ is the soil matric potential (m), b is the Campbell pore-size distribution exponent, T f is the freezing point (taken as 0°C in this example), θ sat is the saturated water content (m 3 / m 3 ), and g is the gravitational acceleration (9.8 m / s 2 ).

[0194] During the freezing process of the soil, the matric potential at the freezing front rapidly decreases, and the unfrozen soil water migrates to the freezing front under the action of the matric potential gradient. The discontinuity of the soil water content at the soil interface between the soil layers makes the traditional water flow description based on the water content gradient unsuitable. Therefore, the matric potential gradient is used to describe the movement of water:

[0195]

[0196] where K eff is the effective hydraulic conductivity (m s -1 ).

[0197] During the freezing process of soil, the hydraulic conductivity of soil will decrease rapidly with the accumulation of ice in the pore, because the flow rate depends on the cross-sectional area of the flow path and the pore geometry. When the soil starts to freeze, the resistance of ice to water transport is considered by introducing a resistance factor E:

[0198]

[0199] where E is the resistance factor, m ice is the ratio of ice content and total water content.

[0200] Assuming that the soil medium is incompressible and ignoring water vapor migration, the basic equation of soil water movement is obtained according to Darcy's law and the principle of mass conservation:

[0201]

[0202] where S is the source and sink term in the water flux.

[0203] S7, calculating runoff using the runoff calculation and confluence module.

[0204] S71, calculating the surface runoff depth and groundwater runoff depth of each grid.

[0205] Snowmelt and rainforest are the sources of infiltration water. The surface water storage is equal to the difference between the sum of rainforest and snowmelt and the total cumulative infiltration water, and the surface runoff is calculated by the following formula:

[0206] R sur = max(0, [(P rain +M snow )-Infil-Pond max ]*(1-glacier frac )+glacier frac *R glacier );

[0207] where P rain is the rainfall of the grid (mm), Pond max is the maximum allowed water depth (mm), Infil is the total cumulative infiltration water (mm), and glacier frac is the proportion of ice area in the grid.

[0208] During the freeze-thaw process of soil, the freezing front or thawing front will temporarily act as an impermeable layer, making it difficult to distinguish the frozen layer water and interflow in the soil. Therefore, the subsurface runoff is used to describe the process of groundwater runoff. The subsurface runoff of the ith layer of soil is:

[0209]

[0210] where θ f is the field capacity, Δz i is the thickness of the ith layer of soil, and α r is the subsurface runoff factor (-).

[0211] The total subsurface runoff is the sum of R sub,i of all soil layers above the bedrock.

[0212] S72, Confluence process calculation:

[0213] The confluence file prepared based on the elevation data is driven by the simulated surface runoff depth and subsurface runoff depth to calculate the slope confluence and river network confluence process using the Lohmann confluence model to obtain the runoff at the outlet of the watershed. In the cold season, the river flow of permafrost watershed is mainly maintained by the water under the frozen layer, which is represented by Q90 in the flow duration curve. Based on the runoff variation characteristics of several rivers in the Siberian and Qinghai-Tibet permafrost regions in recent decades, an empirical relationship between annual average temperature, annual average precipitation, permafrost coverage and water under the frozen layer is established:

[0214]

[0215] where δ, β1, β2 and β3 are fitting parameters, is the multi-year moving average of air temperature (℃), is the multi-year moving average of precipitation (mm), is the multi-year moving average of permafrost coverage (%). The Lohmann confluence model is used to calculate the confluence process. The daily surface runoff depth and subsurface runoff depth are used to drive the Lohmann confluence model to calculate the slope and river network confluence process. The Lohmann confluence model converges the runoff to the outlet of a given grid, and then adds it to the river network confluence model that couples all the grids together. The river network confluence process is estimated by solving the linear Saint-Venant equation to obtain the runoff at the outlet of the watershed. The final runoff at the outlet of the watershed is the sum of the output of the Lohmann confluence model and the output of the water under the frozen layer module, i.e.

[0216]

[0217] where, is the final runoff at the outlet of the watershed (m3 / s)), The watershed outlet runoff (m³) output by the Lohmann confluence model 3 / s)), Water is supplied to the frozen layer.

[0218] S8. Evaluate the simulation performance of the model using model evaluation metrics.

[0219] Model evaluation metrics include:

[0220] Nash efficiency coefficient:

[0221]

[0222] Klin-Gupta efficiency coefficient:

[0223]

[0224] Percentage deviation:

[0225]

[0226] Root mean square error:

[0227]

[0228] Coefficient of determination:

[0229]

[0230] Among them, O i For the observed value, S i These are simulated values. The average of the observed values. σ is the average of the simulated values. O σ is the standard deviation of the observed values. S denoted as the standard deviation of the simulated values, n is the sample size, and R is the correlation coefficient.

[0231] In a specific embodiment, the Yangtze River source basin with the largest permafrost coverage on the Qinghai-Tibet Plateau is selected as the research area, and the soil temperature and soil moisture data of the measured stations in the basin are selected to evaluate the soil water and heat process simulation capability of the modular distributed water and heat coupled hydrological model in cold regions proposed in the application; the GLEAM evapotranspiration data are selected as the reference data, the simulation results of the model proposed in the application and the data results of mainstream GLDAS / Noah, GLDAS / VIC and PML_V2 model are compared and analyzed at the same time with GLEAM, and the evapotranspiration simulation capability of the model proposed in the application is evaluated; the snow depth satellite inversion products such as SSM / I, SSMIS, AMSR-E and AMSR2 are selected to evaluate the snowmelt process simulation capability of the model proposed in the application; the runoff data (2001-2017) of the outlet Zhimenda hydrological station are selected to evaluate the runoff simulation capability of the model proposed in the application.

[0232] The Yangtze River source area is located in the permafrost region in the middle of the Qinghai-Tibet Plateau, with a basin area of 158096km2, an elevation range of 3500-6500m, a permafrost coverage rate of more than 80% in the basin, a permafrost thickness of 10-120m, an active layer thickness of 1-4m, and belongs to a typical permafrost basin. The annual average temperature of the Yangtze River source area is between-5.5 and-1.7℃, and the overall trend of spatial distribution is generally high in the southeast and low in the northwest, with typical inland plateau climate characteristics. The annual precipitation varies between 260-496mm, and mainly (80%) concentrates in summer and autumn, with little snow in winter. The runoff in the basin is the largest from May to September each year, with a lag of one month from the precipitation peak, accounting for about 80% of the total annual runoff. From 1986 to 2009, the annual average runoff of the Zhimenda hydrological station at the outlet of the basin was 13.8x109m3, and the runoff coefficient was 0.20 on average.

[0233] The implementation steps are as follows:

[0234] Step one: collect and organize as Figure 2The data required for running the model in the Yangtze River source basin shown in the application include elevation data, precipitation data, air temperature data, wind speed data, relative humidity data, downward shortwave radiation and downward longwave radiation data, and leaf area index data, and the software required for data processing and model running is installed, such as ArcGIS, Python, etc. Among them, the air temperature, wind speed, incident shortwave radiation and incident longwave radiation data are from the China Regional Surface Meteorological Element Driving Dataset (CMFD), and the relative humidity data is from the ERA5-Land dataset. The thin plate spline curve method is used to interpolate the data of the China Meteorological Bureau station to obtain the precipitation data of the entire Yangtze River source area, and the triple configuration method is used to fuse the interpolated precipitation data with the precipitation data in the CMFD and ERA5-Land datasets for the final driving model. The leaf area index data in the GLASS dataset and the GIMMS dataset are ensemble averaged as the leaf area index data required for the final running of the model. The elevation data, vegetation type data and glacier distribution data are from the National Qinghai-Tibet Plateau Scientific Data Center. The soil texture data and organic matter content data are from the China High-Resolution National Soil Information Grid Basic Attribute Dataset

[0235] Step two: taking the Nash efficiency coefficient (NSE), Kline-Gupta efficiency coefficient (KGE), percentage bias (Pbias), root mean square error (RMSE), determination coefficient (R 2 ) as the objective function, the model parameters are calibrated, and the simulation performance of the model is evaluated from the aspects of soil temperature and humidity, evapotranspiration, snow depth and runoff, etc.

[0236] The simulation results of soil temperature and soil liquid water content after implementation of the technical scheme are shown in Figure 3 and Figure 4 , the model proposed in the application accurately captures the dynamic changes of soil temperature (average R 2 = 0.98, average RMSE = 0.63℃) and liquid water content (average R 2 = 0.87, average RMSE = 0.03m 3 / m -3 ) under different underlying surface conditions, including the latent heat effect in the freezing and thawing process and the ice-water phase change process. The simulation results of total evapotranspiration are shown in Figure 5 , the average RMSE between the monthly evapotranspiration simulated by the model proposed in the application and the GLEAM data is only 15.7mm, and the average R 2 is as high as 0.9, the overall performance is better than GLDAS / Noah, GLDAS / VIC and PML_V2 models, indicating the superiority of the model proposed in the application in the simulation of evapotranspiration in cold regions. The simulation results of snow depth are shown in Figure 6As shown, the average snow depth simulated by the model and satellite inversion in the application is 1.13 cm and 1.31 cm respectively, and the average RMSE between the simulated snow depth and the satellite inversion value is only 0.86 cm, which proves that the model proposed in the application can reasonably simulate and reproduce the temporal and spatial variation of snow depth in the cold region basin. The results of runoff simulation are as shown in the following table 1. Figure 7 As shown, on a monthly scale, the model of the application obtained relatively accurate simulation results during calibration, with NSE of 0.88, KGE of 0.87, R 2 of 0.89, and Pbias of-6.86%. After parameter calibration, the simulation capability of the model is still good, with NSE, KGE, R 2 and Pbias values of 0.85, 0.90, 0.86 and-5.86% respectively. Overall, considering the complex environment, high spatial heterogeneity and data scarcity in the cold region, the runoff simulation capability of the modular distributed cold region water-heat coupled hydrological model proposed in the application is relatively good, and can be used for runoff process simulation in the complex cold region environment. The evaluation index results of the modular distributed cold region water-heat coupled model proposed in the application are shown in table 1.

[0237] Table 1 Evaluation index of the modular distributed cold region water-heat coupled model in runoff simulation in the source region of the Yangtze River

[0238]

[0239] The above shows and describes the basic principles, main features and advantages of the application. It should be understood by those skilled in the art that the application is not limited by the above examples, and the above examples and descriptions in the specification are only to illustrate the principles of the application, and various changes and improvements can be made to the application without departing from the spirit and scope of the application, and these changes and improvements all fall within the scope of the claimed application. The scope of protection of the application is defined by the appended claims and their equivalents.

Claims

1. A modularized distributed hydrothermal coupling hydrological model construction method, characterized in that, The model comprises an energy balance module, an evapotranspiration module, a snow and ice melt water module, a soil water and heat transport and infiltration module, and a runoff calculation and confluence module, and the model construction comprises the following steps: S1. According to the elevation data and the outlet position of the basin, the range of the basin is extracted, and the basin is divided into a plurality of grids, and files required for the confluence process are prepared; S2. Preparing driving data and grid characteristic description data for each grid; S3. Calculating the energy balance process of the soil and vegetation surface by using the energy balance module; S4. Calculating the evapotranspiration process by using the evapotranspiration module; S5. Calculating the ice and snow melting process by using the ice and snow melting module and adopting the degree-day factor method; The melting process of snow is calculated by the following formula: wherein, is the snowmelt amount, is the snowmelt degree-day factor, is the shortwave radiation factor, is the net shortwave radiation on the snowmelt surface, is the air temperature, denotes the maximum temperature threshold when precipitation is snow. The glacier zone method is combined with the degree-day factor method and the area-volume-ratio method to simulate the glacier runoff process, comprising the following steps: Each glacier is divided into a plurality of glacier zones, the temperature and precipitation data of each glacier zone are interpolated, and the area proportion of each glacier zone in each grid is calculated; Based on the temperature and precipitation data of each glacier zone obtained by interpolation, rain and snow are divided, the rainfall and snowfall of each glacier zone are obtained, and the snow water equivalent of each glacier zone is calculated; When the temperature drops below the snowmelt threshold temperature, the potential snowmelt amount is calculated by using the degree-day factor method, and compared with the current snow water equivalent of the glacier zone, if the potential snowmelt amount is less than the snow water equivalent, only snowmelt process occurs, otherwise the underlying glacier will start to melt: wherein, represents the total snow and glacier melt water, is the degree-day factor for glacier melting, is the current snow water equivalent of the glacier, is the potential snow melt, is the air temperature of the glacier, is the melting temperature of the snow or glacier; The runoff depth on each glacier zone is equal to the sum of rainfall and total ice and snow melt, and on the grid scale, the runoff depth of the glacier area is equal to the product of the runoff of all the glacier zones contained in the grid and the area proportion: wherein, is the number of glacier zones, is the number of glacier zones, is the runoff depth of the glacier zone, is the runoff depth of the glacier region, is the precipitation of the glacier zone, is the snowmelt water of the glacier zone, is the area proportion of the glacier zone in the grid; S6. Calculating the soil water infiltration and internal water and heat transport process by using the soil water and heat transport and infiltration module; S7. Calculating the runoff by using the runoff calculation and confluence module; S8. Evaluating the simulation performance of the model by using the model evaluation index.

2. The modular, distributed, coupled-hydraulic-hydrothermal model construction method according to claim 1, wherein, The S3 comprises the following steps: Calculating the net shortwave radiation: wherein, net shortwave radiation absorbed by the vegetated surface, net shortwave radiation absorbed by the bare soil / snow surface, incident shortwave radiation, fraction of shortwave radiation transmitted by the vegetation, albedo of the vegetated surface, albedo of the bare soil / snow surface, fraction of vegetation cover; Calculating the net longwave radiation: wherein, net longwave radiation absorbed by the vegetated surface, net longwave radiation absorbed by the bare soil / snow surface, longwave radiation emitted by the canopy surface, incident longwave radiation, longwave radiation emitted by the bare soil / snow surface; The net radiation of the surface is equal to the sum of the net shortwave radiation and the net longwave radiation: wherein, is the net shortwave radiation, is the net longwave radiation; The formula for calculating the surface heat flux is as follows: wherein, is the surface heat flux.

3. The modular, distributed, cold region hydrothermal coupled hydrological model construction method of claim 2, wherein, The S4 comprises the following steps: Calculating the potential evapotranspiration: wherein is the slope of the curve of the saturation water vapor pressure as a function of temperature, is the available energy, including the latent heat consumed by the ice-water phase change, is the density of air, is the specific heat of air, is the saturation water vapor pressure, is the actual water vapor pressure, is the aerodynamic impedance, is the canopy impedance, is the psychrometric constant, is the density of liquid water, is the latent heat of vaporization; Calculating the interception evaporation of the canopy surface: wherein, is the amount of water retained, is the maximum water retention capacity, is the time step; Calculating the transpiration of the canopy surface: wherein, is the canopy surface drying ratio, is the canopy surface transpiration, denotes the first root fraction in the layer soil layer; Calculating the soil evaporation: wherein, is the liquid water content of the topsoil, is the porosity of the topsoil, is the saturated hydraulic conductivity of the soil, is the inverse of the Campbell pore-size distribution exponent, is the air entry potential.

4. The modular, distributed, cold region hydrothermal coupled hydrological model construction method of claim 3, wherein, The S6 adopts a multi-layer soil structure to reflect the water and heat transport process in the soil, and the division of soil thickness is determined by the exponential layering method; in each time step, the calculation process of the water and heat transport process in the soil layer is as follows: S61. Based on the current soil liquid water content, ice content and organic matter content, the water and heat parameters of the soil are calculated, including thermal conductivity, water conductivity, volume heat capacity and water flux; S62. Solving the heat transport equation to obtain the temperature of each layer of soil; S63. Updating the soil liquid water content and ice content, and the water conductivity and water flux based on the latest calculated soil temperature; S64. Solving the Richard equation considering the ice and water phase change to obtain the water content and ice content of each layer of soil; S65. Repeating steps S61-S64 until the soil temperature and soil water error reaches the convergence standard; S66, the water infiltration process is calculated by using the multilayer Green-Ampt algorithm, and the soil moisture and soil temperature are updated.

5. The modular, distributed, cold region hydrothermal coupled hydrological model construction method of claim 4, wherein, The water and heat transfer process inside the soil layer is calculated by the following formula: For each layer of soil, the change of soil temperature is calculated: wherein, is the volumetric heat capacity of the soil, is the thermal conductivity of the soil, is the density of ice, is the volumetric ice content, is the liquid water flux, is the volumetric heat capacity of liquid water, is the temperature of the soil, is the depth of the soil, is the latent melting potential; The volumetric heat capacity of soil is expressed as the weighted sum of the volumetric heat capacities of the components in the soil: wherein, is the bulk density of the soil, is the specific heat of the mineral matter, is the specific heat of the organic matter, is the volume fraction of air, is the liquid water content (m 3 ·m -3 ), is the volume fraction of organic matter; The thermal conductivity of soil is obtained by the weighted combination of the dry state thermal conductivity and the saturated state thermal conductivity through the Kostensen number: wherein, is a function of the volumetric fraction of soil saturation and organic matter, is the soil dry state thermal conductivity, is the soil saturated state thermal conductivity; where, is the soil saturation, is the volume fraction of sand, is the volume fraction of gravel, and is the adjustment factor; the relationship between is used to predict the continuous change of the entire soil saturation level and the thermal conductivity within the particle size distribution range, which is suitable for the calculation of the thermal conductivity of frozen soil in arid and semi-arid regions; In cold regions, when the soil temperature is above the freezing point, the matric potential is expressed as a function of the liquid water content, when the temperature drops below the freezing point, the matric potential drops rapidly, and is given by the freezing point water potential equation: wherein, is the soil matric potential, is the Campbell pore-size distribution exponent, is the freezing point, is the saturated water content, is the gravitational acceleration; The movement of water is described by the matric potential gradient: wherein, Keffis the effective hydraulic conductivity; When the soil starts to freeze, the retardation effect of soil ice on water transport is considered by introducing a retardation factor: wherein is a retardation factor, is the ratio of ice content and total water content; Assuming that the soil medium is incompressible, ignoring water vapor migration, according to Darcy's law and the principle of mass conservation, the basic equation of soil water movement is obtained: wherein, is the source-sink term in the water flux.

6. The modular, distributed, cold region hydrothermal coupled hydrological model construction method of claim 5, wherein, The S7 includes the following steps: S71, calculate the surface runoff depth and the groundwater runoff depth of each grid; The surface runoff is calculated by the following formula: wherein, is the amount of rainfall for the grid, is the maximum allowed depth of water accumulation, is the total accumulated amount of water seepage, is the proportion of glacier area in the grid; The Subsurface runoff of the layer of soil is: wherein, is the field capacity, is the first is the thickness of the layer of soil, is the subsurface runoff factor; Total runoff from the ground is the sum of all soil layers above the bedrock ; S72, the confluence process calculation: Based on the elevation data prepared, the Lohmann confluence model is used to calculate the slope confluence and river network confluence process, and the outlet runoff of the basin is obtained, driven by the simulated surface runoff depth and groundwater runoff depth; The relationship between the annual average temperature, the annual average precipitation, the multi-year permafrost coverage rate and the water under the frozen layer is established; The final runoff at the outlet of the basin is the sum of the output of the Lohmann confluence model and the output of the water under the frozen layer module, that is: wherein, is the final runoff at the outlet of the catchment, is the runoff at the outlet of the catchment output by the Lohmann flow model, is the water under the frozen layer.

7. The modular, distributed, cold region hydrothermal coupled hydrological model construction method of claim 6, wherein, The model evaluation indexes include Nash efficiency coefficient, Kling-Gupta efficiency coefficient, percentage bias, root mean square error and determination coefficient.

Images