Modularized distributed cold region hydrothermal coupling hydrological model

Abstract

The invention discloses a modularized distributed cold region hydrothermal coupling hydrological model which comprises the following steps: S1, dividing a drainage basin into a plurality of grids, and preparing a file required by a confluence process; s2, preparing driving data and grid feature description data of each grid; s3, calculating an energy balance process of the soil and vegetation surfaces by using an energy balance module; s4, calculating an evapotranspiration process by using an evapotranspiration module; s5, calculating an ice and snow melting process by using an ice and snow melting water module and adopting a degree-day factor method; s6, calculating soil moisture infiltration and an internal hydrothermal transmission process by utilizing a soil hydrothermal transmission and infiltration module; s7, calculating the runoff volume by using a runoff production calculation and confluence module; and S8, evaluating the simulation performance of the model through the model evaluation index. The device has a flexible modular structure, can automatically adjust the time step length, improves the simulation capability of the runoff process in the cold region, and can provide scientific guidance for water resource management in the river source region.

Description

Technical Field

[0001] The present invention relates to the technical field of frozen soil hydrology and basin runoff generation, and particularly to a modular distributed hydrothermal coupling hydrological model in cold regions. Background Art

[0002] As a special regional "impervious layer" or "weak pervious layer", permafrost with its significantly low hydraulic conductivity weakens or hinders the hydraulic connection between precipitation, surface water bodies and groundwater at a certain spatio-temporal scale, strongly affecting the runoff generation process of surface runoff, the migration and distribution patterns of subsurface 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, changes in the distribution of surface and groundwater storage, and so on. Therefore, accurately simulating and predicting the changes of permafrost and its impact on hydrological processes under the background of climate change is crucial for downstream water resources management.

[0003] Traditional hydrological models mostly use some empirical formulas (such as the Stefan equation, the relationship between air temperature and hydraulic conductivity index, etc.) to reflect the impact of permafrost on hydrological processes. Although this simple empirical method that links air temperature and soil freeze-thaw cycles improves the accuracy of runoff simulation in permafrost basins 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 permafrost and hydrological processes and the interaction between permafrost degradation and hydrological processes in a changing environment. Process-based cold region hydrological models, such as VIC, GBEHM, WEB-DHM-pf, etc., often have high requirements for input data and also have certain limitations in application in the Tibetan Plateau region with scarce data. In addition, most process-based models are designed by specific problem-oriented or goal-oriented methods, relying on a single mathematical theory basis, unique parameterization schemes and fixed data requirements. This greatly limits the applicability of the model and its simulation / prediction ability under actual modeling conditions. Summary of the Invention

[0004] In order to solve the above problems, the present invention proposes a modular distributed hydrothermal coupling hydrological model in cold regions to solve the above problems.

[0005] The present application discloses a modular distributed hydrothermal coupling 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 transfer and infiltration module, and a runoff generation calculation and confluence module, including the following steps:

[0006] S1. Based on the elevation data and the location of the watershed outlet, GIS software is used to extract the watershed range and divide the watershed into several grids of a specified size. The flow direction data, grid spherical distance, and grid area ratio (the ratio of grids inside the watershed is 1, and the ratio of grids at the watershed boundary is less than 1) and other files required for the confluence process are prepared;

[0007] S2. Prepare driving data (temperature, precipitation, wind speed, relative humidity, incident shortwave radiation, incident longwave radiation, leaf area index, etc.) and grid feature description data (altitude, soil layer thickness, soil texture, organic matter content, etc.) for each grid;

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

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

[0010] S5. Using the ice and snow meltwater module, the degree-day factor method is used to calculate the process of ice and snow melting;

[0011] S6. Calculate soil water infiltration and internal water and heat transfer process using the soil water and heat transfer and infiltration module;

[0012] S7, using the flow generation calculation and confluence module to calculate the runoff;

[0013] S8. Evaluate the simulation performance of the model through model evaluation indicators.

[0014] Preferably, S3 comprises the following steps:

[0015] Calculate the net shortwave radiation:

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

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

[0018] Among them, 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 proportion of shortwave radiation transmitted by vegetation, α c is the albedo of the vegetation surface, α soil / snowAlbedo is the albedo of the bare soil / snow surface, and FVC is the vegetation coverage;

[0019] Calculate the net long-wave 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] Where, Nl veg is the net long-wave radiation absorbed by the vegetation surface, Nl soil / snow is the net long-wave radiation absorbed by the bare soil / snow surface, L c is the long-wave radiation emitted by the canopy surface, L d is the incident long-wave radiation, L soil / snow is the long-wave radiation emitted by the bare soil / snow surface;

[0023] The net radiation of the surface is equal to the sum of the net short-wave radiation and the net long-wave radiation:

[0024] R net = Ns + Nl;

[0025] Where, Ns is the net short-wave radiation and Nl is the net long-wave radiation;

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

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

[0028] Where, G0 is the surface heat flux.

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

[0030] The total evapotranspiration (E total ) of each grid consists 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. Calculate the root proportion of the vegetation at different depths, and distribute the calculated transpiration amount to each soil layer according to the root proportion. The soil evaporation is completely distributed to the surface soil layer. The specific calculation processes of each evaporation component are as follows:

[0031] Calculate the potential evapotranspiration:

[0032]

[0033] Among them, Δ is the slope of the curve of saturated water vapor pressure varying with temperature, Q ne is the available energy, including the latent heat consumed by the ice - water phase change, ρ 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 psychrometric constant, ρ liq is the density of liquid water, L v is the latent heat of vaporization;

[0034] Calculate the interception evaporation amount on the canopy surface:

[0035]

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

[0037] Calculate the transpiration amount on the canopy surface:

[0038]

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

[0040] Among them, f dry is the dry proportion on the canopy surface, E trans is the transpiration amount on the canopy surface, root frac,i represents the root proportion in the i - th soil layer;

[0041] Calculate the soil evaporation amount:

[0042]

[0043] Among them, θ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 reciprocal of the Campbell pore - size distribution index, Ψ e is the air entry potential.

[0044] Preferably, the said S5 includes the following steps:

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

[0046]

[0047] where M snow is the snowmelt volume, DDF snow is the snowmelt degree-day factor, SRF snow is the shortwave radiation factor, Ns snow is the net shortwave radiation at the snowmelt surface, T air is the air temperature, and T snow represents the maximum temperature threshold when precipitation is snowfall;

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

[0049] Divide each glacier into several glacier zones, interpolate the air temperature and precipitation data for each glacier zone, and calculate the area ratio of each glacier zone within each grid;

[0050] Based on the interpolated air temperature and precipitation data of each glacier zone, conduct rain-snow division to obtain the rainfall and snowfall amounts of each glacier zone, and calculate the snow water equivalent of each glacier zone;

[0051] When the air temperature drops below the snowmelt threshold temperature, use the degree-day factor method to calculate the potential snowmelt volume, and compare it with the snow water equivalent of the current glacier zone. If the potential snowmelt volume is less than the snow water equivalent, only the snowmelt process occurs; otherwise, the underlying glacier will start to melt:

[0052]

[0053] where M band represents the total snow and glacier meltwater, DDF ice is the glacier melting degree-day factor, SWE band is the snow water equivalent of the current glacier zone, M snow,band is the potential snowmelt volume, T band is the air temperature of the glacier zone, and 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 the rainfall and the total ice and snow meltwater. On the grid scale, the runoff depth of the glacier area is equal to the sum of the products of the runoff yields of all glacier zones contained in the grid and their area ratios:

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

[0056]

[0057] Among them, i is the number of glacier zones, n is the number of glacier zones, and R band,i is the runoff depth of the i-th glacier zone, and R glacier is the runoff depth of the glacier area, and P rain,band,i is the rainfall of the i-th glacier zone, and M band,i is the ice and snow melt water volume of the i-th glacier zone, and A frac,i is the proportion of the area of the i-th glacier zone in the grid.

[0058] Preferably, the S6 uses a multi-layer soil structure to represent and reflect the water and heat transfer process inside the soil, and the division of the soil thickness is determined by the 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. The water and heat parameters include thermal conductivity, hydraulic conductivity, volumetric heat capacity, and water flux;

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

[0061] S63. Update the soil liquid water content, ice content, hydraulic 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 soil layer;

[0063] S65. Repeat steps S61 - S64 until the soil temperature and soil moisture errors reach the convergence criterion;

[0064] S66. Calculate the water infiltration process using the multi-layer Green-Ampt algorithm 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 soil layer, calculate the change in soil temperature:

[0067]

[0068] Among them, 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 soil temperature, z is the soil depth, and L f is the melting potential;

[0069] The volumetric heat capacity of the soil is expressed as the weighted sum of the volumetric heat capacities of the 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 is the bulk density of the soil, c m is the specific heat of the mineral matter, c o is the specific heat of the organic matter, θ air is the volume fraction of air, θ liq is the liquid water content (m 3 ·m -3 ), and organic is the volume fraction of the organic matter;

[0072] The thermal conductivity of the soil is obtained by weighted combination of the dry-state thermal conductivity and the saturated-state thermal conductivity with the Kostermann number:

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

[0074] Where K e is a function of the soil saturation and the volume fraction of the organic matter, λ dry is the dry-state thermal conductivity of the soil, and λ sat is the saturated-state thermal conductivity of the soil;

[0075]

[0076] Where S r is the soil saturation, v sand is the volume fraction of sand grains, v gravel is the volume fraction of gravel, and α and β are adjustment factors; the relationship between K e -S r is used to predict the continuous change of the thermal conductivity within the entire soil saturation level and particle size distribution range, and is applicable to the calculation of the thermal conductivity of frozen soil in arid and semi-arid regions;

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

[0078]

[0079] where Ψ is the soil matrix potential, b is the Campbell pore size distribution index, T f is the freezing point, θ sat is the saturated water content, and g is the acceleration due to gravity;

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

[0081]

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

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

[0084]

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

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

[0087]

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

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

[0090] S71. Calculate the surface runoff depth and the subsurface runoff depth for each grid.

[0091] Snowmelt and throughfall are the sources of infiltrating water. The surface water storage is equal to the difference between the sum of throughfall and snowmelt and the total cumulative infiltration water volume. 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] Among them, P rain is the rainfall of the grid, Pond max is the maximum allowable ponding depth (mm), Infil is the total cumulative infiltration volume, glacier frac is the proportion of the glacier area in the grid.

[0094] During the soil freezing and thawing process, the freezing front or thawing front will temporarily act as an aquitard, making it difficult to distinguish the perched water in the frozen layer and the subsurface flow in the soil. Therefore, the subsurface runoff is uniformly used to describe the subsurface runoff generation process. The subsurface runoff of the i-th layer of soil is:

[0095]

[0096] Among them, θ f is the field capacity, Δz i is the thickness of the i-th layer of soil, α 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. Calculation of the confluence process:

[0099] Based on the confluence file prepared from the elevation data, driven by the simulated surface runoff depth and subsurface runoff depth, the Lohmann confluence model is used to calculate the overland flow confluence and river network confluence processes to obtain the runoff at the basin outlet. In the cold season, the river flow in the permafrost basin is mainly maintained by the water under the frozen layer. The Q90 in the flow duration curve is used to represent the water under the frozen layer. Based on the runoff change characteristics of multiple rivers in the permafrost regions of Siberia and the Tibetan Plateau in recent decades, an empirical relationship between the annual average temperature, annual average precipitation, permafrost coverage rate and the water under the frozen layer is established:

[0100]

[0101] Among them, δ, β1, β2 and β3 are fitting parameters, is the multi-year moving average of the air temperature (°C), is the multi-year moving average of the precipitation (mm), is the multi-year moving average of the permafrost coverage rate (%). The Lohmann routing model is used to calculate the routing process. The daily surface runoff depth and subsurface runoff depth are used to drive the Lohmann routing model to calculate the hillslope and river network routing processes. The Lohmann routing model routes the runoff to the outlet of a given grid and then adds it to the river network routing model that couples all grids together. The river network routing process is estimated by solving the linear Saint-Venant equation to obtain the runoff at the basin outlet. The final runoff at the basin outlet is the sum of the output of the Lohmann routing model and the output of the sub-freezing groundwater module, i.e.,

[0102]

[0103] where, is the final runoff at the basin outlet (m 3 / s)), is the runoff at the basin outlet output by the Lohmann routing model (m 3 / s)), is the sub-freezing groundwater.

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

[0105] Advantages of the present invention:

[0106] (1) The model of the present invention has a flexible modular structure, can be integrated with state-of-the-art 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 invention can automatically adjust the time step according to the convergence situation to ensure that the model always converges during the iterative calculation process of soil moisture and heat transfer.

[0108] (3) The model of the present invention comprehensively considers the influence of soil freeze-thaw cycles on infiltration, evapotranspiration, and water and heat transfer processes inside the soil, greatly improving the simulation ability of the runoff process in cold regions.

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

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

[0111] Figure 2It is a general situation map of the Yangtze River source area basin in the embodiment of the present invention.

[0112] Figure 3 It is a comparison diagram of simulated and measured soil temperatures at different depths of multiple boreholes in the Fenghuoshan small watershed in the embodiment of the present invention.

[0113] Figure 4 It is a comparison diagram of simulated and measured soil water contents at different depths of multiple boreholes in the Fenghuoshan small watershed in the embodiment of the present invention.

[0114] Figure 5 It is a comparison diagram of simulated and measured evapotranspiration in the Yangtze River source area from 2001 to 2017 in the embodiment of the present invention.

[0115] Figure 6 It is a comparison diagram of simulated and measured snow depths in the Yangtze River source area from 2001 to 2017 in the embodiment of the present invention.

[0116] Figure 7 It is a comparison diagram of simulated and measured runoff in the Yangtze River source area from 2001 to 2017 in the embodiment of the present invention. Detailed implementation manners

[0117] To make the purpose, technical solutions and advantages of the present application clearer, the following takes examples with reference to the accompanying drawings and further elaborates on the present application in detail.

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

[0119] S1. According to the elevation data and the location of the basin outlet, use GIS software to extract the basin scope, divide the basin into several grids of a specified size, and prepare files required for the confluence process such as flow direction data, grid spherical distance, and grid area ratio (the ratio of internal basin grids is 1, and the ratio of basin boundary grids is less than 1).

[0120] S2. Prepare the driving data (air temperature, precipitation, wind speed, relative humidity, incident short-wave radiation, incident long-wave radiation, leaf area index, etc.) for each grid and the grid characteristic description data (altitude, soil layer thickness, soil texture, organic matter content, etc.).

[0121] S3. Use the energy balance module to calculate the energy balance process on the soil and vegetation surfaces.

[0122] Net short-wave 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 short - wave radiation absorbed by the vegetation surface (W·m -2 ), Ns soil / snow is the net short - wave radiation absorbed by the bare soil / snow surface (W·m -2 ), R s is the incident short - wave radiation (W·m -2 ), τ c is the proportion of short - wave 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 vegetation cover fraction (-).

[0126] Net long - wave 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 long - wave radiation absorbed by the vegetation surface (W·m -2 ), Nl soil / snow is the net long - wave radiation absorbed by the bare soil / snow surface (W·m -2 ), L c is the long - wave radiation emitted by the canopy surface (W·m -2 ), L d is the incident long - wave radiation (W·m -2 ), and L soil / snow is the long - wave radiation emitted by the bare soil / snow surface (W·m -2 ).

[0130] The net radiation at the surface is equal to the sum of the net short - wave radiation and the net long - wave radiation:

[0131] R net= Ns + Nl;

[0132] Among them, Ns represents the net short-wave radiation (W·m -2 ), and Nl represents the net long-wave radiation (W·m -2 ).

[0133] The calculation formula for the surface heat flux is as follows:

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

[0135] S4. Use the evapotranspiration module to calculate the evapotranspiration process.

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

[0137] Calculate the potential evapotranspiration:

[0138]

[0139] Among them, Δ is the slope of the curve of saturated water vapor pressure changing with temperature (Pa / K), Q ne is the available energy (W / m 2 ), including the latent heat consumed by the phase change of ice and water, ρ 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 impedance (s / m), r c is the canopy impedance (s / m), γ is the psychrometric constant (Pa·K -1 ), ρ liq is the density of liquid water (kg / m 3 ), L v is the latent heat of vaporization (MJ / kg).

[0140] Calculate the intercepted evaporation amount on the canopy surface:

[0141]

[0142] Among them, Int sto is the intercepted water discharge (m), Max it is the maximum intercepted water storage capacity (m), Δt is the time step (s), and here E pot needs to take r c as 0 during calculation.

[0143] Calculate the transpiration amount on the canopy surface:

[0144]

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

[0146] Among them, f dry is the dry ratio on the canopy surface, E trans is the transpiration amount on the canopy surface (m), root frac,i represents the root proportion in the i-th soil layer.

[0147] Calculate the soil evaporation amount:

[0148]

[0149] Among them, θ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 reciprocal of the Campbell pore size distribution index, Ψ e is the air entry potential (m), and here E pot needs to take r c as 0 during calculation.

[0150] S5. Use the snow and ice meltwater module and calculate the snow and ice melting process by the degree-day factor method.

[0151] The snow and ice melting process is calculated by the degree-day factor method. Among them, the snow melting process is calculated by the following formula:

[0152]

[0153] Among them, M snow is the snowmelt amount, 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 ),Ns snowFor the net shortwave radiation of the snowmelt surface (W·m -2 ), T air is the air temperature (°C), and T snow represents the maximum temperature threshold (°C) when precipitation is snowfall.

[0154] Combining the "glacier zone" method with the degree-day factor method and the area-volume ratio method to simulate the glacier runoff process, including the following steps:

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

[0156] Based on the air temperature (T band , °C) and precipitation (P band , mm) data of each "glacier zone" obtained by interpolation, conduct "rain-snow" division to obtain the rainfall (P rain,band , mm) and snowfall (P snow,band , mm) of each "glacier zone", and calculate the snow water equivalent of each "glacier zone".

[0157] When the air temperature drops below the snowmelt threshold temperature, the degree-day factor method is used to calculate the potential snowmelt (M snow,band ), and compare it with the snow water equivalent of the current "glacier zone". If the potential snowmelt is less than the snow water equivalent, only the snowmelt process occurs; otherwise, the underlying glacier will start to melt:

[0158]

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

[0160] The runoff depth on each "glacier zone" is equal to the sum of the rainfall and the total ice and snow meltwater. At the grid scale, the runoff depth of the glacier area is equal to the sum of the products of the runoff yields and area ratios of all "glacier zones" included in the grid:

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

[0162]

[0163] where i is the number of the "glacier zone", n is the number of "glacier zones", P rain,band,i is the rainfall (mm) of the i-th "glacier zone", R glacier is the runoff depth (mm) of the glacier area, M band,i is the ice and snow meltwater volume (mm) of the i-th "glacier zone", A frac,i is the area proportion of the i-th "glacier zone" in the grid, R band,i is the runoff depth (mm) of the i-th "glacier zone".

[0164] The key of the above ice and snow ablation simulation algorithm for "glacier zones" lies in the simulation of glacier retreat or advance. Currently, the methods for estimating the area evolution of glaciers can be roughly divided into two categories. One is the ice-flow model that predicts the movement and deformation of ice based on the physical properties of ice, mechanical processes, and surrounding environmental conditions. However, such methods require detailed geometric information of the glacier surface and base, so their application in areas with scarce data has limitations. 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 ), 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 annually according to the simulated glacier ablation situation, and the area of the glacier is calculated based on 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 altitude. When the glacier retreats, the area of the "glacier zone" with the lowest altitude will shrink. If the calculated shrinking area of the glacier is greater than the area of the lowest "glacier zone", the glacier area of this "glacier zone" will be set to 0, and the glacier area of the adjacent upper "glacier zone" will be reduced by the remaining shrinking area. On the contrary, when the glacier advances, the glacier area of the lowest "glacier zone" will increase and be limited within the initially input area, otherwise a new "glacier zone" will be introduced.

[0168] S6. Use the soil water heat transfer and infiltration module to calculate the soil water infiltration and the internal water heat transfer process.

[0169] A multi-layer soil structure is adopted to represent the hydrothermal transport process inside the soil. The division of the 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 scale 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, due to the high sensitivity to meteorological driving, the soil layer closer to the surface is thinner. The basic calculation process of the hydrothermal transport process inside the soil layer within each time step is as follows:

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

[0174] S62. Solve the heat transfer equation to obtain the temperature of each soil layer.

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

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

[0177] S65. Repeat steps S61 - S64 until the soil temperature and soil moisture errors reach the convergence criterion. In this embodiment, the convergence criterion is: ΔT < 0.001 °C, Δθ < 0.01 m 3 / m 3 , where ΔT is the soil temperature error between two consecutive iterations, and Δθ is the soil moisture error between two consecutive 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 hydrothermal transport process inside the soil layer is calculated by the following formula:

[0180] For each soil layer, use the heat conduction - convection equation and consider the ice-water phase change to calculate the change in soil temperature:

[0181]

[0182] Among them, C s is the volumetric heat capacity of the soil (J·kg -1 ·℃ -1 ), λ s is the thermal conductivity of the soil (W·m -1 ·℃ -1 ), ρ ice is the density of ice (kg·m -3 ), θ ice is the volumetric ice content (m 3 ·m -3 ), q l is the liquid water flux (m·s -1 ), C liq is the volumetric heat capacity of liquid water (J·kg -1 ·℃ -1 ), T is the soil temperature (°C), z is the soil depth (m), L f is the melting potential (3.35×10 5 J·kg -1 ).

[0183] The volumetric heat capacity of the soil is expressed as the weighted sum of the volumetric heat capacities of the 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] Among them, ρ b is the bulk density of the soil, c m is the specific heat of minerals (900 J / (kg·°C)), a o is the specific heat of organic matter (1920 J / (kg·°C)), θ air is the volume fraction of air, θ liq is the liquid water content (m 3 ·m -3 ), and organic is the volume fraction of organic matter.

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

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

[0188] where K e is a function of the soil saturation and the volume fraction of 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 volume fraction of sand grains, v gravel is the volume fraction of gravel, and α (taken as 0.24 in this embodiment) and β (taken as 18.1 in this embodiment) are adjustment factors. The relationship between K e -S r is used to predict the continuous variation of the thermal conductivity over the entire soil saturation level and particle size distribution range, and is applicable to the calculation of the thermal conductivity of frozen soil in arid and semi-arid regions.

[0191] 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. However, when the temperature drops below the freezing point, the matric potential (ignoring the osmotic potential) drops rapidly 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 index, T f is the freezing point (taken as 0 °C in this embodiment), θ sat is the saturated water content (m 3 / m 3 ), and g is the acceleration due to gravity (9.8 m / s 2 ).

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

[0195]

[0196] Among them, K eff is the effective hydraulic conductivity (m·s -1 ).

[0197] During the soil freezing process, since the flow velocity depends on the cross-sectional area of the flow path and the pore geometry, the hydraulic conductivity of the soil will decrease rapidly as soil ice accumulates in the pores. When the soil begins to freeze, the retardation effect of ice on water transport is considered by introducing a retardation factor E:

[0198]

[0199] Among them, E is the retardation factor, m ice is the ratio of the ice content to the total water content.

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

[0201]

[0202] Among them, S is the source-sink term in the water flux.

[0203] S7. Calculate the runoff using the runoff generation calculation and the confluence module.

[0204] S71. Calculate the surface runoff depth and the groundwater runoff depth for each grid.

[0205] Snowmelt and throughfall are the sources of infiltrating water. The surface water storage is equal to the difference between the sum of the throughfall and the snowmelt and the total cumulative infiltration water volume. 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] Among them, P rain is the rainfall (mm) of the grid, Pond max is the maximum allowable ponding depth (mm), Infil is the total cumulative infiltration water volume (mm), glacier frac is the proportion of the glacier area in this grid.

[0208] During the soil freeze-thaw process, the freezing front or thawing front will temporarily act as an aquiclude, making it difficult to distinguish between the perched water in the frozen layer and the subsurface flow in the soil. Therefore, subsurface runoff is uniformly used to describe the subsurface runoff process. The subsurface runoff of the i-th layer of soil is:

[0209]

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

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

[0212] S72. Calculation of the confluence process:

[0213] Based on the confluence file prepared from elevation data, driven by the simulated surface runoff depth and subsurface runoff depth, the Lohmann confluence model is used to calculate the overland flow confluence and river network confluence processes to obtain the runoff at the basin outlet. In the cold season, the river flow in the permafrost basin is mainly maintained by the groundwater below the frozen layer. Q90 in the flow duration curve is used to represent the groundwater below the frozen layer. Based on the runoff variation characteristics of multiple rivers in the permafrost regions of Siberia and the Qinghai-Tibet Plateau in recent decades, an empirical relationship between the annual average temperature, annual average precipitation, permafrost coverage rate, and the groundwater below the frozen layer is established:

[0214]

[0215] where δ, β1, β2, and β3 are fitting parameters, is the multi-year moving average of the air temperature (°C), is the multi-year moving average of the precipitation (mm), is the multi-year moving average of the permafrost coverage rate (%). 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 overland and river network confluence processes. 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 grids together. The river network confluence process is estimated by solving the linear Saint-Venant equation to obtain the runoff at the basin outlet. The final runoff at the basin outlet is the sum of the output of the Lohmann confluence model and the output of the groundwater below the frozen layer module, i.e.:

[0216]

[0217] where is the final runoff at the basin outlet (m3 / s)), is the runoff at the watershed outlet output by the Lohmann confluence model (m 3 / s)), is the subsurface flow below the frozen layer.

[0218] S8. Evaluate the simulation performance of the model through model evaluation indicators.

[0219] The model evaluation indicators include:

[0220] Nash efficiency coefficient:

[0221]

[0222] Kling-Gupta efficiency coefficient:

[0223]

[0224] Percentage bias:

[0225]

[0226] Root mean square error:

[0227]

[0228] Coefficient of determination:

[0229]

[0230] Among them, O i is the observed value, S i is the simulated value, is the average value of the observed values, is the average value of the simulated values, σ O is the standard deviation of the observed values, σ S is the standard deviation of the simulated values, n is the number of samples, and R is the correlation coefficient.

[0231] In a specific embodiment, the Yangtze River Source Basin with the largest permafrost coverage rate on the Qinghai-Tibet Plateau is selected as the study area, and the soil temperature and soil water content data of the measured stations in the basin are used to evaluate the simulation ability of the soil hydrothermal process of the modular distributed cold-region hydrothermal coupling hydrological model proposed in this application; the GLEAM evapotranspiration data is selected as the reference data, and the simulation results of the model proposed in this application and the data results of the mainstream GLDAS / Noah, GLDAS / VIC, and PML_V2 models are compared with GLEAM at the same time to evaluate the evapotranspiration simulation ability of the model proposed in this application; the snow depth satellite inversion products such as SSM / I, SSMIS, AMSR-E, and AMSR2 are selected to evaluate the snowmelt process simulation ability of the model proposed in this application; the runoff data (2001 - 2017) of the Zhimenda Hydrological Station at the basin outlet is selected to evaluate the runoff simulation ability of the model proposed in this application.

[0232] The Yangtze River Source area is located in the permafrost region in the central part of the Qinghai-Tibet Plateau. The basin area is 158,096 km², and the altitude ranges from 3,500 to 6,500 m. The permafrost coverage rate in the basin is over 80%, the permafrost thickness ranges from 10 to 120 m, and the active layer thickness ranges from 1 to 4 m, belonging to a typical permafrost basin. The annual average temperature in the Yangtze River Source area ranges from -5.5 to -1.7 °C, showing an overall trend of higher in the southeast and lower in the northwest in spatial distribution, with typical inland plateau climate characteristics. The annual precipitation varies between 260 - 496 mm, and mainly (80%) is concentrated in summer and autumn, with very little snowfall in winter. The runoff in the basin is the largest from May to September every year, lagging behind the precipitation peak by 1 month, accounting for about 80% of the total annual runoff. During 1986 - 2009, the annual average runoff of the Zhimenda Hydrological Station at the basin outlet was 13.8×10⁹ m³, and the average runoff coefficient was 0.20.

[0233] The implementation steps are as follows:

[0234] Step 1: Collect and organize as Figure 2The data required to run the model proposed in this application in the Yangtze River source basin shown include elevation data, precipitation data, air temperature data, wind speed data, relative humidity data, downward short-wave radiation and downward long-wave radiation data, and leaf area index data, and the software required for data processing and model operation, such as ArcGIS, Python, etc., are installed. Among them, the air temperature, wind speed, incident short-wave radiation and incident long-wave radiation data are from the China Meteorological Forcing 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 meteorological stations of the China Meteorological Administration to obtain the precipitation data of the entire Yangtze River source area, and the triple collocation method is used to fuse the interpolated precipitation data with the precipitation data in the CMFD and ERA5-Land datasets for finally driving the model. The leaf area index data in the GLASS dataset and the GIMMS dataset are averaged to be used as the leaf area index data required for the final operation of the model. The elevation data, vegetation type data and glacier distribution data are from the National Tibetan Plateau Data Center of China. The soil texture data and organic matter content data, etc., are from the Basic Attribute Dataset of the Chinese High-Resolution National Soil Information Grid

[0235] Step 2: Using the Nash-Sutcliffe Efficiency (NSE), Kling-Gupta Efficiency (KGE), Percent Bias (Pbias), Root Mean Square Error (RMSE), and Coefficient of Determination (R 2 ) as the objective functions, calibrate the model parameters and evaluate the simulation performance of the model from multiple aspects such as soil temperature and humidity, evapotranspiration, snow depth, and runoff

[0236] The simulation results of the soil temperature and soil liquid water content after the implementation of the technical solution are respectively as shown in Figure 3 and Figure 4 shown. The model proposed in this application accurately captures the dynamic changes of soil temperature (average R 2 =0.98, average RMSE = 0.63°C) and liquid water content (average R 2 =0.87, average RMSE = 0.03 m 3 / m -3 ) under different underlying surface conditions, including the latent heat effect during the freeze-thaw process and the ice-water phase change process. The simulation result of the total evapotranspiration is as shown in Figure 5 shown. The average RMSE between the monthly evapotranspiration simulated by the model proposed in this application and the GLEAM data is only 15.7 mm, and the average R 2 is as high as 0.9. The overall performance is superior to the GLDAS / Noah, GLDAS / VIC, and PML_V2 models, indicating the superiority of the model proposed in this application in the simulation of evapotranspiration in cold regions. The simulation result of the snow depth is as shown in Figure 6As shown, the average snow depths simulated by the model proposed in this application and retrieved by satellite are 1.13 cm and 1.31 cm respectively, and the average RMSE between the simulated snow depth and the satellite retrieved value is only 0.86 cm, which proves that the model proposed in this application can reasonably simulate and reproduce the spatio-temporal variations of snow depth in cold region basins. The results of runoff simulation are as Figure 7 shown. On the monthly scale, the model of the present invention obtained relatively accurate simulation results during the calibration period, with NSE being 0.88, KGE being 0.87, R 2 being 0.89, and Pbias being -6.86%. After parameter calibration, the simulation ability of the model is still quite excellent, and the values of NSE, KGE, R 2 and Pbias are 0.85, 0.90, 0.86, and -5.86% respectively. Generally speaking, considering the complex cold region environment, high spatial heterogeneity, and scarce data, the runoff simulation ability of the modular distributed cold region water-heat coupled hydrological model proposed in this application is quite excellent and can be used for simulating the runoff process in complex cold region environments. The evaluation index results of a modular distributed cold region water-heat coupled model proposed in this application are shown in Table 1.

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

[0238]

[0239] The above has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited by the above embodiments. What is described in the above embodiments and the specification only illustrates the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and these changes and improvements all fall within the scope of the present invention claimed. The scope of protection claimed by the present invention is defined by the appended claims and their equivalents.

Claims

1. A modular distributed cold region water-heat coupling hydrological model, characterized in that: It includes energy balance module, evapotranspiration module, snow and ice melt water module, soil water and heat transfer and infiltration module, and runoff calculation and confluence module, including the following steps: S1. Extract the watershed range according to the elevation data and the location of the watershed outlet, divide the watershed into several grids, and prepare the files required for the confluence process; S2, preparing driving data and grid feature description data for each grid; S3, using the energy balance module to calculate the energy balance process of the soil and vegetation surface; S4, using the evapotranspiration module to calculate the evapotranspiration process; S5. Using the ice and snow meltwater module, the degree-day factor method is used to calculate the process of ice and snow melting; S6. Calculate soil water infiltration and internal water and heat transfer process using the soil water and heat transfer and infiltration module; S7, using the flow generation calculation and confluence module to calculate the runoff; S8. Evaluate the simulation performance of the model through model evaluation indicators.

2. The modular distributed cold region water-heat coupling hydrological model according to claim 1 is characterized in that: The S3 comprises the following steps: Calculate the net shortwave radiation: Ns veg =R s *FVC*[(1-a c )+τ c *(a soil / snow -1)]; Ns soil / snow =R s *(1-a soil / snow )*[(1-FVC)+τ c *FVC]; Among them, 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 proportion of shortwave radiation transmitted by vegetation, α c is the albedo of the vegetation surface, α soil / snow is the albedo of bare soil / snow surface, and FVC is the vegetation cover; Calculate the net longwave radiation: Nl veg =FVC*[L d +L soil / snow -2*L c ]; Nl soil / snow =(1-FVC)*L d +FVC*L c -IT soil / snow ]; Among them, 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 long-wave radiation emitted by the canopy surface, L d is the incident long-wave radiation, L soil / snow is the longwave radiation emitted by the bare soil / snow surface; The net radiation at a surface is equal to the sum of the net shortwave radiation and the net longwave radiation: R net =Ns+Nl; Among them, Ns is the net shortwave radiation, and Nl is the net longwave radiation; The surface heat flux calculation formula is as follows: G0=0.35462*R net -47.79008; Where G0 is the surface heat flux.

3. The modular distributed cold region water-heat coupling hydrological model according to claim 2 is characterized in that: The S4 comprises the following steps: Calculate potential evapotranspiration: Where Δ is the slope of the curve of saturated water vapor pressure changing with temperature, Q ne is the available energy, including the latent heat consumed by the ice-water phase transition, ρ 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 psychrometric constant, ρ liq is the density of liquid water, L v is the latent heat of vaporization; Calculate the intercepted evaporation at the canopy surface: Among them, Int sto To intercept the water volume, Max it is the maximum interception and storage capacity, Δt is the time step; Calculate transpiration from the canopy surface: E trans =E pot *f dry *root frac,i *Δt; Among them, f dry is the canopy surface dryness ratio, E trans is the transpiration of the canopy surface, root frac,i represents the proportion of roots in the i-th soil layer; Calculate soil evaporation: Among them, θ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 index, Ψ e For the air to enter.

4. The modular distributed cold region water-heat coupling hydrological model according to claim 3 is characterized in that: The S5 comprises the following steps: The melting process of snow is calculated by the following formula: Among them, M snow is the amount of snowmelt, DDF snow is the snowmelt degree day factor, SRF snow is the shortwave radiation factor, Ns snow is the net shortwave radiation from the snowmelt surface, T air is the temperature, T snow Indicates 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, including the following steps: Each glacier is divided into several 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 interpolated temperature and precipitation data of each glacier belt, rain and snow are divided, the rainfall and snowfall in each glacier belt are obtained, and the snow water equivalent of each glacier belt is calculated; When the temperature drops below the snowmelt threshold, the degree-day factor method is used to calculate the potential snowmelt and compare it with the snow water equivalent of the current glacier zone. If the potential snowmelt is less than the snow water equivalent, only snowmelt will occur, otherwise the underlying glacier will begin to melt: Among them, M band represents the total snow cover and glacier meltwater, DDF ice is the glacier melting degree-day factor, SWE band is the current snow water equivalent in the glacier belt, M snow,band is the potential snowmelt, T band is the temperature in the glacier zone, T melt is the melting temperature of snow or glacier; The runoff depth of each glacier belt is equal to the sum of rainfall and total snowmelt water. At the grid scale, the runoff depth of the glacier area is equal to the sum of the product of the runoff of all glacier belts contained in the grid and the area proportion: R band,i =P rain,band,i +M band,i ; Where i is the number of the glacier belt, n is the number of glacier belts, R band,i is the runoff depth of the ith glacier zone, R glacier is the runoff depth in the glacier area, P rain,band,i is the rainfall in the ith glacier zone, M band,i is the amount of meltwater from ice and snow in the ith glacier zone, A frac,i is the area proportion of the i-th glacier belt in the grid.

5. The modular distributed cold region water-heat coupling hydrological model according to claim 4 is characterized in that: The S6 uses a multi-layer soil structure to reflect the water and heat transfer process inside the soil. The division of soil thickness is determined by the exponential layering method. The calculation process of the water and heat transfer process inside the soil layer in each time step is as follows: S61. Calculate the hydrothermal parameters of the soil based on the current soil liquid water content, ice content, and organic matter content. The hydrothermal parameters include thermal conductivity, hydraulic conductivity, volumetric heat capacity, and water flux. S62, solving the heat transfer equation to obtain the temperature of each layer of soil; S63, updating soil liquid water content and ice content as well as hydraulic conductivity and water flux based on the latest calculated soil temperature; S64. Solve the Richard equation considering the ice-water phase transition to obtain the water content and ice content of each layer of soil; S65, repeating steps S61-S64 until the soil temperature and soil moisture errors reach the convergence standard; S66. Use the multi-layer Green-Ampt algorithm to calculate the water infiltration process and update the soil moisture and soil temperature.

6. The modular distributed cold region water-heat coupling hydrological model according to claim 5 is characterized in that: The water and heat transfer process inside the soil layer is calculated by the following formula: For each soil layer, calculate the change in soil temperature: Among them, C s is the volume heat capacity of the soil, λ s is the thermal conductivity of soil, ρ ice is the density of ice, θ ice is the volume ice content, q l is the liquid water flux, C liq is the volume heat capacity of liquid water, T is the temperature of the soil, z is the depth of the soil, L f To melt potential; The volumetric heat capacity of soil is expressed as the weighted sum of the volumetric heat capacities of the components in the soil: C s =(1-organic)*ρ b *c m +organic*p b *c o +θ liq *r liq *c liq +θ ice *r ice *c ice +θ air *r air *c air ; Among them, ρ b is the soil bulk density, c m is the specific heat of the mineral, c o is the specific heat of organic matter, θ air is the volume fraction of air, θ liq is the liquid water content (m 3 ·m -3 ), organic is the volume fraction of organic matter; The thermal conductivity of the soil is obtained by weighting the dry state thermal conductivity and the saturated state thermal conductivity through the Kosten number: l s =(λ sat -l dry )*K e +λ dry ; Among them, K e is a function of soil saturation and volume fraction of organic matter, λ dry is the thermal conductivity of soil in dry state, λ sat is the thermal conductivity of soil in saturated state; Among them, S r is the soil saturation, v sand is the volume fraction of sand particles, v gravel is the volume fraction of gravel, α and β are adjustment factors; K e -S r The relationship between the thermal conductivity of the soil and the particle size distribution is used to predict the continuous change of thermal conductivity in the whole soil saturation level and particle size distribution range, which is suitable for the calculation of thermal conductivity of frozen soil in arid and semi-arid areas. In cold regions, when soil temperature is above freezing, the matrix potential is expressed as a function of liquid water content. When the temperature drops below freezing, the matrix potential drops rapidly and is given by the freezing point water potential equation: Where Ψ is the soil matrix potential, b is the Campbell pore size distribution index, T f is the freezing point, θ sat is the saturated water content, g is the gravitational acceleration; The movement of water is described using the matrix potential gradient: Among them, K eff is the effective hydraulic conductivity; When the soil begins to freeze, the retardation effect of soil ice on water transport is taken into account by introducing a retardation factor: Among them, E is the blocking factor, m ice is the ratio of ice content to total water content; Assuming that the soil medium is incompressible and ignoring water vapor migration, the basic equation of soil moisture movement is obtained according to Darcy's law and the principle of mass conservation: Where S is the source and sink term in the moisture flux.

7. The modular distributed cold region water-heat coupling hydrological model according to claim 6 is characterized in that: The S7 comprises the following steps: S71, calculating the surface runoff depth and underground runoff depth of each grid; Surface runoff is calculated using the following formula: R sur =max(0,[(P rain +M snow )-Infil-Pond max ]*(1-glacier frac )+glacier frac *R glacier ); Among them, P rain is the rainfall of the grid, Pond max is the maximum allowable water depth, Infil is the total accumulated water seepage, glacier frac is the percentage of glacier area in the grid; The underground runoff of the i-th soil layer is: Among them, θ f is the field water capacity, Δz i is the thickness of the i-th soil layer, α r is the underground runoff factor; Total groundwater runoff is the R of all soil layers above the bedrock. sub,i The sum of S72, confluence process calculation: Based on the confluence file prepared by elevation data, the Lohmann confluence model is used to calculate the slope confluence and river network confluence process driven by the simulated surface runoff depth and underground runoff depth, and the basin outlet runoff is obtained; Establish the relationship between annual average temperature, annual average precipitation, permafrost coverage and water below the frozen layer; The final runoff at the basin outlet is the sum of the output of the Lohmann confluence model and the output of the frozen layer water module, that is: in, is the final runoff at the basin outlet, is the outlet runoff of the basin output by the Lohmann confluence model, Water for the frozen layer.

8. The modular distributed cold region water-heat coupling hydrological model according to claim 7 is characterized in that: The model evaluation indicators include Nash efficiency coefficient, Kling-Gupta efficiency coefficient, percentage deviation, root mean square error and determination coefficient.

Images