A glacier meltwater runoff simulation method based on degree-day factor seasonal differentiation

Abstract

The application discloses a glacier meltwater runoff simulation method based on degree-day factor seasonal differentiation, which comprises the following steps: meteorological data processing, construction of a runoff calculation model with degree-day factor seasonal differentiation, self-adaptive calibration and verification of the runoff calculation model parameters, daily runoff process simulation of the glacier coverage grid scale, and quantitative analysis of the glacier meltwater contribution. In the application, the degree-day factor is improved from a fixed value in a year to a value taken according to the ablation stage, the degree-day factor in each stage is used as an independent parameter for calibration, and the glacier surface albedo is preferentially used as an ablation stage division index. A composite objective function containing a glacier meltwater contribution rate deviation penalty term is constructed, and a differentiated calibration strategy is designed for different hydrological years, so that the systematic deviation caused by the fixed degree-day factor can be eliminated, and the precision and rationality of the glacier meltwater runoff simulation are improved.

Description

Technical Field

[0001] This invention relates to a method for simulating glacial meltwater runoff based on seasonal variations of day-degree factors, belonging to the field of hydrological runoff simulation technology in cold regions. Background Technology

[0002] Generally speaking, the freezing and melting processes of glaciers and snow cover are an important source of river runoff and are crucial to the water resource security of watersheds. With global warming, the runoff contribution patterns of glacial meltwater are changing significantly.

[0003] Distributed hydrological (VIC) models are commonly used tools for simulating runoff in cold regions. The glacier diurnal factor is a key parameter in glacier meltwater simulation, its physical meaning being the glacier ablation depth corresponding to a unit of positive accumulated temperature. Hock R. (2003), after summarizing observations of multiple glaciers, pointed out that glacier ablation efficiency is not constant during the ablation period—in the early stages of ablation, the glacier surface is mostly covered by snow, albedo is high, and ablation efficiency is low; during the peak of ablation, bare ice is exposed, albedo is low, and ablation efficiency increases significantly; at the end of ablation, as solar radiation weakens, ablation efficiency decreases again.

[0004] Existing technologies have coupled the day-delay factor glacier ablation module into the VIC model, enabling grouped simulations of precipitation runoff, snowmelt runoff, and glacier runoff (SuT. et al., 2023). Furthermore, calculations of glacier area proportion at the grid scale based on glacier cataloging data, as well as differentiated runoff generation calculations between glacier-covered and non-glacier-covered grids, have also been applied in existing technologies, such as the technical solution disclosed in publication number CN116089765A.

[0005] However, existing technologies only describe the seasonal differences in day-degree factors from an observational perspective, without developing an engineering scheme that embeds phased day-degree factors into the VIC model, calibrates each phase's day-degree factor as an independent parameter, and optimizes parameters based on hydrological year-specific differences. Due to this lack of technical guidance, in existing VIC-day-degree factor coupled hydrological models, day-degree factors for the same glacier basin region are set to fixed values ​​throughout the year, failing to reflect the seasonal variation in ablation efficiency. In other words, simulating the entire ablation period with fixed day-degree factors will overestimate glacial meltwater volume in the early stages of ablation and underestimate it during the peak of ablation, leading to systematic biases in the simulation of the glacial meltwater process.

[0006] In addition, the existing model parameter calibration technology has two shortcomings: First, it adopts a uniform calibration strategy across the entire domain and does not set different search rules for different hydrological year types such as high water years, normal water years, and low water years; second, the objective function only relies on the total runoff fitting index and does not include the rationality of the contribution rate of glacial meltwater in the calibration constraints, which may lead to a good fit of the total runoff but a deviation in the runoff component classification.

[0007] Therefore, how to solve the problems of systematic bias caused by fixed degree-day factors, mismatch between calibration strategies and hydrological year patterns, and insufficient constraints of the objective function on runoff components without changing the basic structure of the VIC model is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0008] The purpose of this invention is to systematically solve the above three technical problems and provide a method for simulating glacial meltwater runoff based on the seasonal differentiation of the day-degree factor. By improving the day-degree factor from a fixed value throughout the year to a value differentiated according to the ablation stage, introducing an adaptive calibration strategy based on hydrological year type, and constructing a composite objective function that includes a penalty term for the deviation of glacial meltwater contribution rate, systematic biases in the simulation of glacial meltwater processes are avoided.

[0009] This invention achieves the above objective through the following technical solution: a method for simulating glacial meltwater runoff based on seasonal differentiation of day-degree factors, which includes the following steps: S1. Meteorological data processing: Collect multi-source meteorological data of the target area, and obtain a standardized meteorological dataset with grid scale through bias correction and interpolation processing; S2. Construct a runoff calculation model with seasonally differentiated day-degree factors: Based on the VIC model, a glacier meltwater calculation module is coupled to the glacier cover grid. The glacier ablation stage is divided into the initial ablation stage, the peak ablation stage, and the final ablation stage, and day-degree factors are set for different stages accordingly, thereby realizing the runoff calculation of the glacier cover grid. Glacier cover grids are obtained by spatially overlaying grids and glacier vector data of the target area. The diurnal factor of the glacier cover grid is determined according to the glacier basin zone to which the glacier cover grid belongs. The diurnal factor of the same glacier basin zone is determined according to the glacier ablation stage. The diurnal factor of each glacier basin zone is divided into multiple time windows according to the glacier ablation stage, and each time window takes a different value to form a seasonal differentiation scheme of diurnal factor. S3. Adaptive calibration and verification of runoff calculation model parameters: The hydrological year type is divided into high-water, normal-water, and low-water years using the guarantee rate division method. Calibration parameters are selected, and the runoff calculation model parameters for the corresponding hydrological year type are calibrated and verified through composite objective function and intelligent optimization algorithm to obtain the optimal runoff calculation model for the hydrological year type. S4. Simulation of daily runoff processes at the scale of glacier-covered grid: The standardized meteorological dataset is input into the optimal runoff generation calculation model of the hydrological year type to simulate the daily runoff processes of the glacier-covered grid and output the simulation results of the glacier-covered grid. S5. Quantitative analysis of the contribution of glacial meltwater: Based on the simulation results, the runoff contribution of glacial meltwater is quantified, and the glacial meltwater replenishment rate at different time periods is calculated.

[0010] As a further technical solution of the present invention, in S1, deviation correction refers to establishing a linear regression correction model based on meteorological data and reanalysis data and satellite data, and the correction formula is: Y_corr = a × Y_raw + b in, Y_corr For the corrected data, Y_raw The original data, a、b These are the regression coefficients; Interpolation processing refers to interpolating bias-corrected meteorological data onto a grid using Kriging interpolation to obtain a standardized meteorological dataset with grid scale.

[0011] As a further technical solution of the present invention, the glacial meltwater calculation module adopts the following formula: Glacier meltwater volume:

[0012] in, DDF j (t) For the first j The first glacier basin division in the 1st t The time factor corresponding to the time period T i (t) For the first i The glacier cover grid in the first t The average daily surface temperature of the glacier during that period; Snowmelt runoff:

[0013] in, S i ( t ) is the first i The grid in the first t Snowmelt runoff during the period Q melt Net energy available for snow melting p_ w The density of water, L_f The latent heat of melting ice; Total runoff:

[0014] in, R i (t) For precipitation, For the first i The percentage of glacier area covered by each glacier grid P i (t) For precipitation, M i,j (t) This refers to the amount of glacial meltwater. S i (t) For snowmelt runoff, R vic,i (t) The runoff depth in non-glacial areas simulated by the VIC model.

[0015] As a further technical solution of the present invention, in S2, the glacier ablation stage is divided into the initial ablation stage, the peak ablation stage, and the final ablation stage by using the glacier surface albedo: when the glacier surface albedo is higher than the first threshold, it is classified as the initial ablation stage; when the glacier surface albedo is lower than the second threshold, it is classified as the peak ablation stage; when the glacier surface albedo is between the first threshold and the second threshold, it is classified as the final ablation stage. The first and second thresholds are determined based on the actual albedo variation characteristics of the glaciers in the target area.

[0016] As a further technical solution of the present invention, in S2, for target areas lacking observation data on glacier surface albedo, the glacier ablation stages are divided using calendar months: the initial ablation period is from May to June, the peak ablation period is from July to August, and the final ablation period is in September.

[0017] As a further technical solution of the present invention, in S3, the composite objective function includes a runoff fitting error term and a glacier meltwater contribution rate deviation penalty term, specifically:

[0018] in, Object The objective function value, NSE The Nash efficiency coefficient. PBIAS This represents the percentage deviation in water volume. R glacier_sim To simulate the contribution rate of glacial meltwater, R glacier_ref Reference value for the contribution rate of glacial meltwater. This represents the penalty coefficient for each item.

[0019] As a further technical solution of the present invention, in S3, the intelligent optimization algorithm designs differentiated parameter adaptive calibration strategies for different hydrological year types as follows: Using the standard initial parameter settings of the intelligent optimization algorithm as the default settings, the parameter search step size is reduced by a smaller ratio compared to the default settings in dry years, while the search weight of the degree-day factor is increased by a larger ratio compared to the default settings; the infiltration parameter search weight is increased by a larger ratio compared to the default settings in wet years; and the default settings are used in normal years. The proportion of reducing the search step size for parameters in dry years and increasing the search weight for degree-day factors, as well as the proportion of increasing the search weight for infiltration parameters in wet years, are all determined based on the hydrological characteristics of the target area.

[0020] As a further technical solution of the present invention, in S5, quantifying the runoff contribution of glacial meltwater includes: calculating the proportion of glacial meltwater runoff at each measured station in wet years, normal years, and dry years from the perspective of hydrological year type.

[0021] As a further technical solution of the present invention, in S5, the formula for calculating the glacial meltwater replenishment rate is:

[0022] in, R recharge For glacial meltwater replenishment rate, R glacier This refers to the runoff of glacial meltwater. R total For total runoff, R precip This is to measure the precipitation flow rate.

[0023] The beneficial effects of this invention are: (1) This invention improves the degree-day factor from a fixed value throughout the year to a value differentiated according to the ablation stage, and uses the degree-day factor of each stage as an independent parameter in the calibration. It also prioritizes the use of glacier surface albedo as the indicator for dividing the ablation stage. This effectively eliminates the systematic bias of overestimating the fixed degree-day factor in the early stage of ablation and underestimating the glacier meltwater during the peak of ablation. The improvement effect of this bias is particularly significant in watersheds with a high proportion of glacier area. The examples show that the simulation accuracy of the peak of ablation is significantly improved.

[0024] (2) By introducing an adaptive calibration strategy based on hydrological year type, this invention strengthens the optimization of glacial meltwater parameters in dry years and the optimization of precipitation runoff parameters in wet years, thereby improving the simulation stability of the model under different hydrological year types.

[0025] (3) This invention constructs a composite objective function that includes a penalty term for the deviation of glacial meltwater contribution rate. During the calibration process, it simultaneously constrains the fitting accuracy of total runoff and the rationality of runoff component classification. It also designs differentiated calibration strategies for different hydrological year types, avoiding the problem of "good fitting of total runoff but inaccurate simulation of glacial meltwater proportion". It can eliminate the systematic deviation caused by fixed degree-day factors, improve the simulation reliability of different water source component contribution rates, provide more accurate quantitative basis for the refined assessment of water resources in high-altitude and cold regions, and improve the accuracy and rationality of glacial meltwater runoff simulation. Attached Figure Description

[0026] Figure 1 This is a schematic diagram of the overall process of the present invention; Figure 2 This is a comparison chart of the simulation accuracy of the scheme in Embodiment 3 of the present invention and the scheme with fixed degree-day factor; Figure 3 This is a comparison chart of the simulated glacial meltwater volume deviation between the scheme of Embodiment 3 of the present invention and the scheme of fixed degree-day factor. Detailed Implementation

[0027] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0028] Example 1, as Figure 1 As shown, this embodiment provides a method for simulating glacial meltwater runoff based on seasonal differentiation of day-degree factors. This method specifically includes the following steps: First, meteorological data processing: collect available meteorological data from the target area (high-altitude cold region), and obtain a standardized meteorological dataset with grid scale after bias correction and interpolation.

[0029] Second, a runoff calculation model with seasonal differentiation of day-degree factors is constructed to realize runoff calculation for glacier-covered grids.

[0030] This runoff calculation model is based on the VIC model, retaining its energy balance snowmelt module, and coupling a glacier cover grid with a glacier meltwater calculation module. The glacier cover grid is obtained by spatially overlaying an N×N grid of data representing the target region with glacier vector data. The diurnal factor of the glacier cover grid is determined according to its glacier basin region. Within the same glacier basin region, the diurnal factor is divided into multiple time windows according to the glacier ablation stage, taking different values ​​(diurnal factors for different stages). Each diurnal factor for each ablation stage is further optimized as an independent parameter.

[0031] Based on the grid division of the target region, the glacier area percentage of each glacier-covered grid is calculated through spatial overlay analysis. ,Will The runoff calculation logic is embedded for the corresponding glacier-covered grid. The runoff calculation method for glacier-covered grids includes the contribution of glacier meltwater, while the runoff calculation for non-glacier-covered grids follows the original VIC model.

[0032] In this step, the glacial meltwater calculation module uses the following formula: Glacier meltwater volume:

[0033] in, DDF j (t) For the first j The first glacier basin division in the 1st t The time factor corresponding to the time period T i (t) For the first i The glacier cover grid in the first t The average daily surface temperature of the glacier during that period; Snowmelt runoff:

[0034] in, S i ( t ) is the first i The grid in the first t Snowmelt runoff (mm) during the time period. Q melt Net energy used for snow melting (MJ / m²) p_w The density of water (1000 kg / m³). L_f The latent heat of melting of ice (0.334 MJ / kg); Total runoff:

[0035] in, R i (t) For precipitation, For the first i The percentage of glacier area covered by each glacier grid P i (t) For precipitation, M i,j (t) This refers to the amount of glacial meltwater. S i (t)For snowmelt runoff, R vic,i (t) The runoff depth in non-glacial areas simulated by the VIC model.

[0036] In this step, there are two ways to divide the glacier ablation stages: The first method divides the ablation phase into three stages based on observed or remotely sensed albedo values ​​of the glacier surface: (1) Initial stage of melting: The albedo of the glacier surface is higher than the first threshold. At this time, the glacier surface is mostly covered by snow, the albedo is high, the melting efficiency is low, and the day factor takes a lower value. (2) Peak ablation period: The albedo of the glacier surface is below the second threshold. At this time, the bare ice is exposed, the albedo is low, the ablation efficiency is high, and the day factor takes a higher value. (3) End of ablation: The albedo of the glacier surface is between the first threshold and the second threshold, or the solar radiation is lower than the preset threshold, the ablation slows down, and the day factor takes the middle or lower value.

[0037] The second approach is to divide the glacier ablation phase into calendar months for target areas lacking glacier surface albedo observation data: the initial ablation period is from May to June, the peak ablation period is from July to August, and the final ablation period is in September.

[0038] It should be noted that the first and second thresholds are determined based on the actual albedo variation characteristics of glaciers in high-altitude and cold regions. For high-altitude and cold regions with available albedo observation data, the albedo index is preferentially used for stage division.

[0039] The core of this invention lies in the "seasonal differentiation" of the day-degree factor—that is, assigning different values ​​based on the natural stages of glacier melting, rather than being limited to specific calendar months. The division of melting stages should essentially be based on the physical indicator of glacier surface albedo. However, for high-altitude and frigid regions lacking continuous albedo observations or remote sensing inversion data, calendar months in a climatological statistical sense can be used as an approximate division method.

[0040] The core improvement in this step is that the day-temperature factor of each glacier basin is divided into multiple time windows according to the glacier ablation stage, and each time window takes a different value, thus forming a scheme of seasonal differentiation of day-temperature factor. At the same time, the day-temperature factor of each glacier ablation stage is optimized separately as an independent parameter in the next calibration process.

[0041] It should be noted that the physical basis for the seasonal differentiation of the day-day factor lies in the seasonal variation of glacier surface albedo, rather than the seasonal variation of air temperature. Specifically, in the day-day factor model, the formula for calculating glacier meltwater is:

[0042] wherein, PDD is the positive degree-day, which reflects the energy available for ablation; DDF is the degree-day factor, which reflects the amount of ice melted per unit energy.

[0043] The seasonal variation of positive degree-days is already embodied in PDD. The independent physical significance of the seasonal differentiation of the degree-day factor lies in: in the initial ablation period, the snow cover on the glacier surface leads to high albedo, and with the same input of positive degree-days, the proportion of energy actually absorbed by the glacier for ablation is lower, so DDF should take a lower value; in the peak ablation period, the exposure of bare ice leads to low albedo, and with the same input of positive degree-days, the proportion of energy actually absorbed is higher, so DDF should take a higher value. Therefore, the seasonal differentiation of the degree-day factor is not a double counting of the seasonal effect of temperature, but an independent parameterization of the difference in energy utilization efficiency caused by the seasonal variation of albedo on the glacier surface.

[0044] Third, adaptive calibration and verification of parameters of the runoff generation calculation model.

[0045] Based on the measured runoff data, the guarantee rate division method is adopted to complete the division of hydrological year types into wet years, normal years and dry years. Soil parameters and the degree-day factors of each ablation stage in each glacial watershed zone are selected as calibration parameters, and the composite objective function and intelligent optimization algorithm are used to complete the calibration and verification of the parameters of the runoff generation calculation model, so as to obtain the optimal runoff generation calculation model for the hydrological year types.

[0046] In this step, a manner of completing the division of hydrological year types into wet years, normal years and dry years by adopting the guarantee rate division method based on measured runoff data is provided: sort the measured annual total runoff of n years in descending order, and calculate the empirical frequency P of each year as P=m / (n+1)×100%, where m is the sorting rank. Wet years correspond to the years with P≤37.5%, normal years correspond to the years with 37.5%<P≤62.5%, and dry years correspond to the years with P>62.5%.

[0047] In this step, there are two core improvements: First, a composite objective function including a runoff fitting error term and a deviation penalty term for the glacial meltwater contribution rate is constructed, which is specifically as follows:

[0048] wherein, Object is the objective function value, NSE is the Nash-Sutcliffe efficiency coefficient, PBIAS is the percentage of water volume deviation (substituted in decimal form, for example, when the water volume error is 15%, PBIAS it is taken as 0.15), R glacier_sim is the simulated glacial meltwater contribution rate, R glacier_ref is the reference value of the glacial meltwater contribution rate obtained by a reference method independent of simulation, and the reference method is an isotope tracer method or a water balance method, These are the penalty coefficients for the total runoff trend error term, the water balance error term, and the glacier meltwater contribution rate deviation term, respectively. They are determined based on the hydrological characteristics of the target area (high-altitude cold region) and the magnitude differences of each error term, and their values ​​range from 0.1 to 1.0.

[0049] Further explanation of the above formula: (1- NSE ) and | PBIAS | This represents the runoff fitting error term, which measures the model's accuracy in simulating total runoff and water balance;| R glacier_sim -R glacier_ref The term "penalty term" is used to further constrain the simulated glacial meltwater contribution rate from deviating significantly from the reference value, in addition to the runoff fitting objective.

[0050] Regarding NSE, PBIAS , R glacier_sim The calculation formula is:

[0051] in, Q obs (t) For the first t Measured runoff during the time period Q sim (t) For the first t Simulated runoff over a given period of time Q obs_m (t) The average measured runoff; the contribution rate of simulated glacial meltwater over a single time period. R glacier_sim (t) Calculated from the ratio of glacial meltwater flow to total watershed outflow during the same period. Q glacier (t) For the first t The time-period model outputs glacial meltwater runoff. Q total ( t ) is the first t The total outflow of all water sources in the basin during the specified time period is calculated; the average annual contribution rate of glacial meltwater under this condition is obtained by summing the integrals of all time periods for a single hydrological year. R glacier_sim Substitute the parameters into the upper-level objective function to complete the parameter optimization with component penalty constraints.

[0052] Consider (1- NSE), | PBIAS |、| R glacier_sim - R glacier_ref The fluctuation ranges of the three factors differ, and their contribution percentages in the parameter optimization process are determined by balancing the runoff fitting error, water balance error, and meltwater component deviation using three sets of independent weights; among them... R glacier_sim - R glacier_ref The term "|" represents a penalty constraint for meltwater components, used to add additional loss constraints to parameter solutions that deviate from the measured component reference values. This invention divides long-sequence runoff into three hydrological scenarios: high-water years, normal-water years, and low-water years. Weight combinations are configured to address the differences in runoff recharge sources under different scenarios; the value range for all three weights is set to 0.1. Version 1.0 allows for differentiated value assignment based on the runoff recharge characteristics of different hydrological scenarios, such as high-water years, normal-water years, and low-water years.

[0053] The weighting coefficients for each component are set differently based on the hydrological year type: Taking a normal water year as the baseline scenario, λ1=λ2=λ3=1 / 3, ensuring balanced weights for all components. In dry years, glacial meltwater accounts for a high proportion of water supply, so λ3 is increased to 0.50, while λ1 and λ2 are each decreased to 0.25, strengthening the constraint on glacial meltwater components. In wet years, rainfall is the main source of water supply, so λ2 is increased to 0.50, while λ1 and λ3 are each decreased to 0.25, strengthening the constraint on water balance. The adjustment range for each weight does not exceed ±50% of the baseline value, ensuring that each error always retains a perceptible contribution to the objective function.

[0054] This composite objective function ensures that the calibration process not only pursues the fitting accuracy of the total runoff volume, but also constrains the simulated runoff component classification to approach the measured reference value.

[0055] Reference value R of glacial meltwater contribution rate glacier_ref The methods for obtaining the reference values ​​are as follows: The isotope tracing method involves collecting river water and samples from various water sources, determining the stable isotopic composition, and using an isotope mixing model to segment the runoff contribution ratio of each water source; the water balance method, based on measured glacier mass balance data, precipitation data, and evaporation data, estimates the glacier meltwater volume and its contribution ratio to the total runoff using a water balance equation. The acquisition of these reference values ​​is independent of the simulation process in this embodiment and can be obtained through literature review or field sampling analysis before parameter calibration.

[0056] Secondly, the intelligent optimization algorithm employs a parameter adaptive calibration strategy tailored to the different hydrological year types: The standard initial parameter settings of the intelligent optimization algorithm are used as default settings. In dry years, the parameter search step size is reduced relative to the default settings, while the degree-day factor search weight is increased relative to the default settings to enhance the optimization of glacial meltwater parameters during dry years. In wet years, the infiltration parameter search weight is increased relative to the default settings to enhance the optimization of precipitation runoff parameters during wet years. In normal years, the default settings (default search step size and default parameter weights) are used to optimize all parameters in a balanced manner. It is important to note that the reduction in the parameter search step size and the increase in the degree-day factor search weight during dry years, as well as the increase in the infiltration parameter search weight during wet years, are all determined based on the hydrological characteristics of the target area.

[0057] It should be noted that the calibration parameters include: variable infiltration curve parameters, nonlinear baseflow rate ratio, maximum baseflow rate, nonlinear baseflow soil moisture content ratio, second soil aquifer thickness, third soil aquifer thickness, and degree-day factor for each ablation stage in each glacier basin region.

[0058] Furthermore, the runoff calculation model is validated periodically and during the verification period to obtain the optimal parameter set of the runoff calculation model for the corresponding hydrological year type, and thus obtain the optimal runoff calculation model for the hydrological year type.

[0059] Fourth, simulation of daily runoff processes at the glacier cover grid scale.

[0060] A standardized meteorological dataset is input into the optimal runoff generation calculation model based on the hydrological year type to simulate the daily runoff process of each glacier cover grid and output the simulation results for each glacier cover grid. The simulation results include snowmelt runoff, glacier meltwater runoff, and total runoff.

[0061] Fifth, a quantitative analysis of the contribution of glacial meltwater.

[0062] Based on the simulation results of each glacier cover grid, the runoff contribution of glacier meltwater is quantified, and the glacier meltwater replenishment rate at different times under different hydrological year types is calculated.

[0063] Quantifying the runoff contribution of glacial meltwater involves calculating the proportion of glacial meltwater runoff at each measured station in wet, normal, and dry years from the perspective of hydrological year type.

[0064] The formula for calculating the glacial meltwater replenishment rate is:

[0065] in, R recharge For glacial meltwater replenishment rate, R glacier This refers to the runoff of glacial meltwater. Rtotal For total runoff, R precip This represents the precipitation runoff. The glacier meltwater replenishment rate. R recharge Unlike the full-caliber glacial meltwater contribution rate used in the parameter calibration stage. R glacier_sim To eliminate the interference of rainfall runoff, this study is used to refine the analysis of the contribution of glaciers to non-rainfall runoff at the hydrological annals level.

[0066] Example 2: This example provides a method for simulating glacial meltwater runoff based on the seasonal variation of day-degree factors, applied to a watershed upstream of a hydrological station in the Yangtze River source region (geographical location 90°43′E~96°45′E, 32°30′N~35°35′N, watershed area approximately 137,000 km²). The specific steps are as follows: (I) Multi-source meteorological data processing The following data were collected: daily precipitation and temperature observation data from 1961 to 2018 from three meteorological stations: Tuotuohe, Wudaoliang, and Qumalai; CMFD reanalysis dataset from 1979 to 2018 (spatial resolution 0.1°×0.1°); AVHRR land cover data from the University of Maryland at 1km resolution; FAO global soil data; daily runoff data from a hydrological station from 1968 to 2017; daily runoff data during the flood season from the Tuotuohe hydrological station from 1959 to 2017; and glacier cataloging and fusion data from the Third Pole region.

[0067] Bias correction is performed on multi-source meteorological data. The bias correction method used in this embodiment is as follows: A linear regression correction model is established using station observation data as a benchmark, combining reanalysis data and satellite data. The correction formula is: Y_corr = a × Y_raw + b in, Y_corr For the corrected data, Y_raw The original data, a、b is the regression coefficient.

[0068] The corrected data was interpolated to a 9km×9km grid using Kriging interpolation to obtain a standardized meteorological dataset.

[0069] (ii) Constructing a runoff calculation model with seasonal differentiation of day-degree factors

[0070] The study area was divided into 2217 9km × 9km grids. Glacier vector data was spatially overlaid with the grids to obtain 163 glacier-covered grids, and the glacier area percentage of each grid was calculated. f i .

[0071] Based on the spatial distribution of glaciers, they are divided into three glacier basin regions. Each glacier basin region adopts a scheme based on seasonal variation of the degree-day factor. According to the actual climatic characteristics of the Yangtze River source region (a high-altitude, cold region with a long ablation period), the ablation period (May to September) is divided into three stages: (1) Early stage of melting (May to June): The glacier surface is mostly covered with snow, with high albedo and low melting efficiency, and the day factor takes a low value. (2) Peak ablation period (July to August): Bare ice is exposed, albedo is low, ablation efficiency is high, and the day factor takes a higher value; (3) End of ablation (September): Solar radiation weakens, ablation slows down, and the day factor takes the middle value.

[0072] The range of DDF (Day Degree Factor) values ​​for each stage in each glacier basin region is shown in Table 1. Table 1 shows the day-degree factor values ​​for each stage in each glacier watershed region.

[0073] Note: Units are mm / (℃·d)

[0074] It should be noted that for high-altitude and cold regions with continuous observation or remote sensing inversion data of glacier surface albedo, more precise stage division can be performed based on albedo thresholds. For example, a first threshold of 0.45 can be set, with an albedo above 0.45 classified as the initial stage of ablation; a second threshold of 0.30, with an albedo below 0.30 classified as the peak stage of ablation; and an albedo between 0.30 and 0.45 or a significant decrease in solar radiation classified as the final stage of ablation. The calendar month division method used in this embodiment is merely an exemplary approach under conditions where albedo data is lacking.

[0075] The day-of-day factor for each ablation stage is optimized separately as an independent parameter in the next calibration step. The original energy balance snowmelt module of the VIC model is retained and coupled with the glacial meltwater calculation module.

[0076] The calculation process for the glacial meltwater calculation module is as follows: Glacier meltwater volume:

[0077] in, DDF j (t) For the first j The first glacier basin division in the 1st t The degree-day factor corresponding to the time period is determined according to the seasonal differentiation scheme mentioned above, which varies with the ablation stage. T i (t) For the firsti The glacier cover grid in the first t The average daily surface temperature of the glacier during that period; Snowmelt runoff:

[0078] in, S i ( t ) is the first i The grid in the first t Snowmelt runoff (mm) during the time period. Q melt Net energy available for snow melting (MJ / m²). p_w The density of water (1000 kg / m³). L_f The latent heat of melting of ice (0.334 MJ / kg); The snowmelt module built into the VIC model can effectively simulate the accumulation and melting of snow. This embodiment retains this module for snowmelt runoff calculation. The core logic of the snowmelt module is to calculate the energy budget of the snow surface at each time step within each glacier cover grid. When snow is present and the surface temperature reaches 0°C, the net energy input exceeding the latent heat of melting drives snowmelt. The liquid water produced by snowmelt mixes with concurrent liquid precipitation and participates in soil water movement and runoff generation according to the stratified infiltration mechanism of the VIC model, in conjunction with precipitation. P i (t) and glacial meltwater M i,j ( t Together, they constitute the total runoff input of the glacier cover grid.

[0079] Total runoff:

[0080] in, R i (t) For precipitation, For the first i The percentage of glacier area covered by each glacier grid P i (t) For precipitation, M i,j (t) This refers to the amount of glacial meltwater. S i (t) For snowmelt runoff, R vic,i (t) The runoff depth in non-glacial areas simulated by the VIC model.

[0081] The glacier-covered grid calculates the total runoff using the above-mentioned area-weighted method, and the non-glacier-covered grid follows the original runoff yield calculation of the VIC model. The cross-section where this hydrological station is located is the confluence outlet of the entire watershed.

[0082] (3) Adaptive calibration and verification of runoff yield calculation model parameters

[0083] Based on the measured runoff data of a certain hydrological station from 1980 to 2017, the guarantee rate division method is adopted to classify hydrological year types: 12 wet years (guarantee rate P≤37.5%), 13 normal years (37.5%<P≤62.5%), and 13 dry years (P>62.5%).

[0084] Selected calibration parameters include: the variable infiltration curve parameter INFILT, the nonlinear baseflow rate ratio Ds, the maximum baseflow rate Ds_MAX, the nonlinear baseflow soil moisture content ratio Ws, the thickness of the second soil aquifer DEPTH_2, the thickness of the third soil aquifer DEPTH_3, and the degree-day factor DDF for each ablation stage in each glacial watershed subregion.

[0085] A composite objective function is constructed, and an intelligent optimization algorithm is used for parameter calibration. The composite objective function is as follows:

[0086] where NSE is the Nash-Sutcliffe efficiency coefficient, PBIAS is the percentage of water volume deviation, R glacier_sim is the simulated glacier meltwater contribution rate, R glacier_ref is the reference value of the glacier meltwater contribution rate obtained by a reference method independent of simulation.

[0087] In this embodiment, the reference value of the glacier meltwater contribution rate R glacier_ref adopts the isotope tracing result of a certain hydrological station in the source region of the Yangtze River recommended in the literature, which is about 22%.

[0088] In this embodiment, l the value rule is as follows: (1) Normal year: l 1= l 2= l 3=1 / 3, the weights of the three are relatively balanced, and the total runoff, water balance, and glacier meltwater component are constrained synchronously; (2) Wet year: Rainfall is the main recharge source, focusing on ensuring the accuracy of the annual water balance, and increasing l 2, set l 1=0.25, l 2=0.5, l3 = 0.25; (3) High proportion of glacial meltwater supply in dry years: Glacial meltwater is the core supply source, and it is necessary to strengthen the fitting effect of meltwater components and adjust upward. l 3. Take l 1 = 0.25 l 2 = 0.25 l 3=0.5, which amplifies the penalty for loss due to component deviation.

[0089] Different parameter adaptive calibration strategies are designed for different hydrological year types: the standard initial parameter settings of the CMA-ES algorithm are used as the default settings, the parameter search step size is reduced by 50% and the search weight of the degree-day factor is increased by 30% in dry years; the parameter search step size is kept at the default in wet years and the search weight of INFILT is increased by 30%; and the default search step size and default parameter weights are used in normal years.

[0090] Two independent calibration periods and validation periods were set up: the first calibration period was from 2001 to 2005, and the validation period was from 2006 to 2009; the second calibration period was from 2013 to 2015, and the validation period was from 2016 to 2017. After calibration, the optimal parameter set for the corresponding hydrological year type was obtained, and then the optimal runoff calculation model for the hydrological year type was obtained.

[0091] (iv) Simulation of daily runoff processes at the glacier cover grid scale

[0092] Standardized meteorological datasets were input into the optimal runoff calculation model corresponding to the hydrological year type to simulate the daily runoff process of 2217 grids from 1980 to 2017. Each glacier-covered grid outputs three types of data: snowmelt runoff, glacier melt runoff, and total runoff. The daily runoff process at the watershed scale of a certain hydrological station and the Tuotuohe station was obtained through the confluence model.

[0093] (v) Quantitative analysis of the contribution of glacial meltwater.

[0094] Based on the simulation results, the proportion of glacier runoff and snowmelt runoff were calculated.

[0095] Glacier runoff percentage = (Glacier meltwater runoff / Total runoff) × 100%

[0096] Snowmelt runoff percentage = Snowmelt runoff / Total runoff × 100%

[0097] The results show that at this hydrological station, glacial runoff accounts for approximately 26.9% and snowmelt approximately 9.8% in dry years, glacial runoff approximately 22.6% and snowmelt approximately 9.9% in normal years, and glacial runoff approximately 21.0% and snowmelt approximately 10.0% in wet years. At the Tuotuohe station, glacial runoff accounts for approximately 62.0% and snowmelt approximately 5.2% in dry years.

[0098] Simultaneously, the glacial meltwater recharge rate during the non-flood season in a dry year is calculated using the following formula: Glacier meltwater replenishment rate = glacier meltwater runoff / (total runoff - precipitation yield) × 100%.

[0099] The results showed that the glacial meltwater replenishment rate during the non-flood season in dry years was about 58.7%, indicating that glacial meltwater during the non-flood season is the main source of river runoff.

[0100] Example 3: To verify the technical effect of Example 2, this example is compared with a comparative example using a fixed degree-day factor.

[0101] like Figure 2 and Figure 3 As shown, the comparative design is as follows: the degree-day factor for each glacier watershed zone is set to a fixed value throughout the year, without distinguishing between ablation stages. To ensure fairness in the comparison, the fixed degree-day factor is selected according to conventional practices in existing technologies. During the calibration process, the DDF is optimized as a single constant parameter, and the fixed DDF value obtained from the calibration is the median of the calibration DDF at the initial and peak ablation stages of each zone. The calibration of the comparison group adopts a conventional unified search strategy and an objective function that only includes runoff fitting error. Other conditions (meteorological driving data, VIC model basic parameters, glacier area proportion, and division of calibration and validation periods) are exactly the same for both schemes, and the results are shown in Tables 2 and 3.

[0102] Table 2 shows the comparison of simulation accuracy.

[0103] Table 3 shows the analysis of the deviation direction in the simulation of glacier meltwater volume.

[0104] The above results indicate that the fixed degree-day factor scheme suffers from significant systematic biases—overestimation during the initial stage of ablation and underestimation during the peak of ablation. This proposed scheme effectively eliminates these systematic biases and significantly improves the accuracy of glacial meltwater simulation by combining three techniques: seasonal differentiation of degree-day factors, adaptive calibration of hydrological year patterns, and a composite objective function.

[0105] The above are merely preferred embodiments of the present invention and do not limit the scope of the patent. The specific DDF value ranges, month division methods, and accuracy comparison data in the above embodiments are merely exemplary implementations of the present invention in the Yangtze River source area and do not constitute a limitation on the scope of protection of the claims. The application scenarios of the present invention are not limited to the Yangtze River source area. For other high-altitude cold regions (such as the Tianshan Mountains, Kunlun Mountains, Karakoram Mountains, Himalayas, Altai Mountains, Qilian Mountains, and other Asian high-altitude cold basins, as well as the Alps, Andes Mountains, Arctic Islands, and other global high-altitude cold basins), the albedo index division, alternative index division, or calendar month division method can be flexibly selected according to the actual ablation period length of local glaciers, albedo variation characteristics, and data availability, and the degree-day factor values ​​for each stage can be determined accordingly. All equivalent structural transformations made using the content of the present invention under the inventive concept of the present invention, or direct / indirect applications in other related technical fields, are included within the scope of patent protection of the present invention.

[0106] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A method for simulating glacial meltwater runoff based on the seasonal variation of day-degree factors, characterized in that, The method for simulating glacial meltwater runoff includes the following steps: S1. Meteorological data processing: Collect meteorological data of the target area and obtain a standardized meteorological dataset with grid scale through bias correction and interpolation; S2. Construct a runoff calculation model with seasonally differentiated day-degree factors: Based on the VIC model, a glacier meltwater calculation module is coupled to the glacier cover grid. The glacier ablation stage is divided into the initial ablation stage, the peak ablation stage, and the final ablation stage, and day-degree factors are set for different stages accordingly, thereby realizing the runoff calculation of the glacier cover grid. S3. Adaptive calibration and verification of runoff calculation model parameters: The hydrological year type is divided into high-water, normal-water, and low-water years using the guarantee rate division method. Calibration parameters are selected, and the runoff calculation model parameters for the corresponding hydrological year type are calibrated and verified through composite objective function and intelligent optimization algorithm to obtain the optimal runoff calculation model for the hydrological year type. S4. Simulation of daily runoff process at the scale of glacier cover grid: Input the standardized meteorological dataset into the optimal runoff generation calculation model of the hydrological year type to simulate the daily runoff process of the glacier cover grid and output the simulation results of the glacier cover grid. S5. Quantitative analysis of glacial meltwater contribution: Based on the simulation results, the runoff contribution of glacial meltwater is quantified, and the glacial meltwater replenishment rate at different time periods is calculated.

2. The method for simulating glacial meltwater runoff according to claim 1, characterized in that, In S1, the deviation correction refers to establishing a linear regression correction model based on meteorological data and reanalysis data and satellite data. The correction formula is as follows: Y_corr = a × Y_raw + b in, Y_corr For the corrected data, Y_raw This is the original data. a、b These are the regression coefficients; The interpolation process refers to interpolating the bias-corrected meteorological data to a grid using Kriging interpolation to obtain a standardized meteorological dataset with a grid scale.

3. The method for simulating glacial meltwater runoff according to claim 1, characterized in that, In S2, the glacier cover grid is obtained by spatially superimposing the grid and glacier vector data of the target area. The day-degree factor of the glacier cover grid is determined according to the glacier watershed partition to which the glacier cover grid belongs. The day-degree factor of the same glacier watershed partition is determined according to the glacier ablation stage.

4. The method for simulating glacial meltwater runoff according to claim 3, characterized in that, The glacial meltwater calculation module uses the following formula: Glacier meltwater volume: in, DDF j (t) For the first j The first glacier basin division in the 1st t The time factor corresponding to the time period T i (t) For the first i The glacier cover grid in the first t The average daily surface temperature of the glacier during that period; Snowmelt runoff: in, S i ( t ) is the first i The grid in the first t Snowmelt runoff during the period Q melt Net energy used for snow melting ρ_w The density of water, L_f The latent heat of melting ice; Total runoff: in, R i (t) For precipitation, For the first i The percentage of glacier area covered by each glacier grid P i (t) For precipitation, M i,j (t) This refers to the amount of glacial meltwater. S i (t) For snowmelt runoff, R vic,i (t) The value represents the runoff depth in non-glacial areas simulated by the VIC model.

5. The method for simulating glacial meltwater runoff according to claim 1, characterized in that, In S2, the glacier ablation stage is divided into the initial ablation stage, the peak ablation stage, and the final ablation stage by using the glacier surface albedo: when the glacier surface albedo is higher than the first threshold, it is classified as the initial ablation stage; when the glacier surface albedo is lower than the second threshold, it is classified as the peak ablation stage; when the glacier surface albedo is between the first and second thresholds, it is classified as the final ablation stage. The first threshold and the second threshold are determined based on the actual albedo variation characteristics of the glacier in the target area.

6. The method for simulating glacial meltwater runoff according to claim 1, characterized in that, In S2, for target areas lacking glacier surface albedo observation data, the calendar months are used to divide the glacier ablation stages: the initial ablation period is from May to June, the peak ablation period is from July to August, and the final ablation period is in September.

7. The method for simulating glacial meltwater runoff according to claim 1, characterized in that, In S3, the composite objective function includes a runoff fitting error term and a glacier meltwater contribution rate deviation penalty term, specifically: in, Obj The objective function value, NSE The Nash efficiency coefficient. PBIAS This represents the percentage deviation in water volume. R glacier_sim To simulate the contribution rate of glacial meltwater, R glacier_ref Reference value for the contribution rate of glacial meltwater. This represents the penalty coefficient for each item.

8. The method for simulating glacial meltwater runoff according to claim 1, characterized in that, In S3, the intelligent optimization algorithm designs differentiated parameter adaptive calibration strategies for different hydrological year types as follows: Using the standard initial parameter settings of the intelligent optimization algorithm as the default settings, the parameter search step size is reduced by a smaller ratio compared to the default settings in dry years, while the search weight of the degree-day factor is increased by a larger ratio compared to the default settings; the infiltration parameter search weight is increased by a larger ratio compared to the default settings in wet years; and the default settings are used in normal years. The proportion of reducing the search step size for parameters in dry years and increasing the search weight for degree-day factors, as well as the proportion of increasing the search weight for infiltration parameters in wet years, are all determined based on the hydrological characteristics of the target area.

9. The method for simulating glacial meltwater runoff according to claim 1, characterized in that, In S5, the quantitative contribution of glacial meltwater to runoff includes: calculating the proportion of glacial meltwater runoff at each measured station in wet, normal, and dry years from the perspective of hydrological year type.

10. The method for simulating glacial meltwater runoff according to claim 1, characterized in that, In S5, the formula for calculating the glacial meltwater replenishment rate is: in, R recharge For glacial meltwater replenishment rate, R glacier This refers to the runoff of glacial meltwater. R total For total runoff, R precip This is to measure the flow rate of precipitation.

Images