Flood forecasting method suitable for arid and semi-arid plain areas

Abstract

The application discloses a flood forecasting method suitable for arid and semi-arid plain areas, which comprises the following steps: constructing effective rainfall based on hydro-meteorological and underlying surface data, and carrying out super-permeation runoff generation discrimination to separate surface runoff and actual total infiltration; according to a dynamic bypass mechanism, the total infiltration is decomposed into matrix domain infiltration water and preferential flow infiltration water, the former is filled layer by layer to generate matrix storage fullness runoff, and the latter is divided into groundwater recharge and deep backfill based on penetration attenuation characteristics; river inflow is obtained by aggregating each runoff component, and after deducting the dynamic transmission loss determined by the initial wicking characteristics of the riverbed, confluence calculation is carried out to output the flood forecasting result. The application can realize full-path water closure, can depict the flood reduction effect of thickened aeration zone and dry river channel, and can improve the prediction accuracy of flood peak flow and peak occurrence time.

Description

Technical Field

[0001] This invention relates to the field of hydrological forecasting technology, and in particular to a flood forecasting method applicable to arid and semi-arid plains. Background Technology

[0002] Currently, the hydrological community has developed relatively mature forecasting methods for rainfall-runoff processes in humid and semi-humid regions. Modern hydrological forecasting models are typically constructed based on saturation-runoff mechanisms or infiltration-excess runoff mechanisms. These models assume that rainfall infiltration occurs uniformly layer by layer within the soil profile, and that a fixed loss coefficient is used to reduce water loss during the evolution of the flow in the river channel. This traditional model is generally applied to natural watersheds with shallow groundwater levels, uniform soil matrix, and perennial river flow, and its predicted flood volume and peak timing are relatively reliable.

[0003] However, some arid and semi-arid plains have long been affected by climate change and human activities, such as the over-extraction of deep groundwater, which has altered the underlying surface structure and runoff generation and confluence response characteristics of their watersheds. When faced with sudden rainstorms following such droughts, current forecasting models have limited capacity to characterize infiltration and diversion mechanisms in complex media, as well as the loss processes of water flow as it evolves on dry surfaces or in river networks. This leads to distortions in the predicted flow process lines in terms of peak flood magnitude and temporal phase. Summary of the Invention

[0004] Purpose of the invention: To propose a flood forecasting method applicable to arid and semi-arid plains, in order to solve the above-mentioned problems.

[0005] Technical solution: Flood forecasting methods applicable to arid and semi-arid plains, including:

[0006] Acquire hydrological and meteorological observation data and underlying surface data for the target watershed, and calculate the effective rainfall for the time period;

[0007] Based on the effective rainfall and underlying surface data for a given period, excess infiltration runoff is determined to obtain the excess surface runoff and the actual total infiltration.

[0008] Based on the dynamic bypass mechanism determined by the effective rainfall and hydrological and meteorological observation data for a given period, the actual total infiltration is decomposed into matrix domain infiltration and preferential inflow infiltration.

[0009] The infiltration water volume in the matrix domain was treated by layer-by-layer soil filling to obtain the matrix full runoff component;

[0010] Based on the underlying surface data, the penetration attenuation characteristics are determined, and the preferential inflow infiltration volume is divided into preferential groundwater recharge and deep recharge. The deep recharge volume is used for runoff calculation in subsequent periods.

[0011] The inflow rate of the river is obtained by summing the excess surface runoff, the matrix saturation runoff, and the preferential groundwater recharge. The dynamic transmission loss during the evolution process is calculated and deducted. Then, the river confluence calculation is performed, and the flood flow process forecast results of the forecast section are output.

[0012] Beneficial effects: This invention can achieve full-path water volume closure, characterize the flood reduction effect of thickened vadose zone and dried river channel, and improve the accuracy of flood peak flow and peak occurrence time forecasts. Attached Figure Description

[0013] Figure 1 This is a flowchart of the flood forecasting method applicable to arid and semi-arid plains regions according to the present invention.

[0014] Figure 2 This is a flowchart illustrating the calculation of effective rainfall over a given period, as presented in this invention.

[0015] Figure 3 This is a schematic diagram illustrating the dynamic discrimination of dual-super-flow generation according to the present invention.

[0016] Figure 4 This is a schematic diagram of the normalized distribution curve of the infiltration capacity of the present invention.

[0017] Figure 5 This is a flowchart illustrating the implementation method of the flood forecasting method of the present invention. Detailed Implementation

[0018] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. 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 should fall within the scope of protection of the present invention.

[0019] The applicant conducted an in-depth analysis and discovered the following problems with flood forecasting in arid and semi-arid plains: The continuous decline in groundwater levels in these plains leads to thickening of the vadose zone. After dehydration, the zone shrinks and forms rapid infiltration channels such as fissures, large pores, and root canals, causing rainfall infiltration to deviate from the assumption of uniform, layer-by-layer infiltration. In other words, some rainfall bypasses the layer-by-layer transport of the soil matrix, rapidly descending along infiltration channels and replenishing groundwater.

[0020] To solve these problems, combined with Figures 1 to 5 The present invention will be specifically described through the following embodiments. According to one aspect of this application, this embodiment provides a flood forecasting method, such as... Figure 1 As shown, it specifically includes:

[0021] Step 101: Obtain hydrological and meteorological observation data and underlying surface data of the target watershed, and calculate the effective rainfall for the time period.

[0022] In this embodiment, the effective rainfall for a target watershed is constructed based on hydrological and meteorological observation data. The target watershed can be an arid or semi-arid plain area and its surrounding piedmont transition zone. The hydrological and meteorological observation data includes multi-point rainfall sequences collected in real time by rain gauges and evaporation data obtained from meteorological stations. The underlying surface data includes topographic data, soil type distribution, land use types, and river network geometry within the target watershed. Based on this, through spatial interpolation of meteorological elements and evapotranspiration reduction processing, the net water volume that can participate in the runoff generation process at each calculation time step is obtained, i.e., the effective rainfall for that time period. This embodiment can eliminate water lost during rainfall due to vegetation interception, depression filling, etc.

[0023] Step 102: Based on the effective rainfall and underlying surface data for the time period, determine the infiltration runoff to obtain the infiltration surface runoff and the actual total infiltration.

[0024] In this embodiment, excess infiltration runoff is determined based on the effective rainfall over a given time period and the infiltration capacity represented by underlying surface data. Infiltration capacity refers to the maximum rate at which a unit area of ​​soil can absorb and infiltrate rainfall under the current soil moisture condition. Specifically, the surface runoff distribution is determined by comparing rainfall intensity with the current infiltration capacity. If the rainfall intensity exceeds the soil's infiltration capacity, the excess rainfall accumulates on the surface, generating excess surface runoff. The amount of water entering the soil is defined as the actual total infiltration, reducing flood peak forecasting errors caused by unclear runoff mechanisms. The formula for calculating the effective rainfall over a given time period in this embodiment is:

[0025] P _net =ΔR _s +ΔF _0 ;

[0026] Among them, P _net ΔR represents the effective rainfall over a given period. _s ΔF represents the excess surface runoff. _0 This represents the actual total infiltration volume.

[0027] Step 103: Based on the dynamic bypass mechanism determined by the effective rainfall amount over the time period and hydrological and meteorological observation data, or in other words, based on the dynamic bypass mechanism determined by the rainfall intensity obtained from the analysis of the effective rainfall amount over the time period and the drought state characterized by the hydrological and meteorological observation data, the actual total infiltration amount is decomposed into matrix domain infiltration amount and preferential inflow infiltration amount using the dynamic bypass mechanism.

[0028] In this embodiment, parameters representing the development degree of macropores or fissures in the soil are used to segment the path of water entering the soil, i.e., a dynamic bypass mechanism. The pre-drought state can be quantified by the number of rainless days before a rainfall event, which determines the development degree of drying shrinkage fissures formed by dehydration and shrinkage in cohesive soil.

[0029] Within the soil, the actual total infiltration will be divided. Part of the water moves slowly along the tiny pores of the soil, forming the matrix infiltration volume; the other part quickly flows downwards through cracks, root canals, etc., forming the preferential infiltration volume. The corresponding water volume relationship is as follows:

[0030] ΔF _0 =F _matrix +F _pref ;

[0031] Among them, F _matrix F represents the amount of water infiltrating into the matrix domain. _pref This is to prioritize the inflow of infiltrated water. This embodiment is used to address the problem that traditional hydrological models cannot reflect the phenomenon of rapid vertical transport of rainfall infiltrates due to fracture development under conditions of thickened vadose zones in the region.

[0032] Step 104: The infiltration water volume in the matrix domain is treated by layer-by-layer soil filling to obtain the matrix full runoff component.

[0033] Specifically, following the order from surface to depth, the matrix infiltration water sequentially fills the upper, lower, and deep water storage capacities of the soil. During the filling process, if the water storage capacity of each soil layer reaches saturation, the newly introduced water becomes free water, which signifies that the matrix is ​​fully saturated with runoff. This embodiment is used to assess the runoff response of a watershed, which together with infiltrated runoff.

[0034] Step 105: Based on the underlying surface data, determine the penetration attenuation characteristics, and accordingly divide the preferential inflow infiltration volume into preferential groundwater recharge and deep recharge; the deep recharge is used for runoff calculation in subsequent periods. That is, based on the groundwater depth in the underlying surface data, determine the penetration attenuation characteristics, and according to these characteristics, divide the preferential inflow infiltration volume into preferential groundwater recharge and deep recharge; then, add the deep recharge to the deep soil water storage to participate in the runoff calculation for subsequent periods.

[0035] In this embodiment, penetration attenuation characteristics are used to characterize the attenuation of preferential inflow seepage volume as it passes through the thick vadose zone and reaches the groundwater surface, due to absorption by the surrounding matrix. In areas with significant groundwater depth, the preferential flow channel is unlikely to penetrate the entire vadose zone. Therefore, this embodiment calculates the penetration efficiency based on the current groundwater depth in the underlying surface data. The portion of the preferential inflow seepage volume that reaches the groundwater surface is taken as the preferential flow groundwater recharge, while the portion that does not reach the surface is taken as the deep recharge and recharged into the deep soil. This embodiment corrects the assumption that the preferential flow recharges the entire groundwater by using a depth response factor, making the groundwater recharge prediction closer to reality.

[0036] Step 106: Add up the excess surface runoff, the matrix saturation runoff, and the preferential groundwater recharge to obtain the river inflow.

[0037] River inflow refers to the total amount of water flowing into the river network within a unit of calculation, after deducting losses such as evaporation. It is calculated by summing the surface, subsurface, and groundwater runoff components to form the input boundary for river channel calculation, thus connecting slope runoff to river confluence.

[0038] Step 107: Calculate and deduct the dynamic transmission loss during the evolution process. That is, calculate the dynamic transmission loss of the river inflow during the evolution process based on the initial infiltration characteristics of the riverbed in the underlying surface data. After deducting the dynamic transmission loss from the river inflow, perform river confluence calculation and output the flood flow process forecast results of the forecast section.

[0039] To address the long-term dryness of rivers in arid and semi-arid plains during the non-flood season, the initial infiltration characteristics of the riverbed are used to characterize the strong absorption capacity of the dry riverbed due to matrix suction in the initial stage of flood contact. The dynamic transport loss of the riverbed changes dynamically with the evolution of the riverbed from dry to wet, i.e., calculating the energy loss and infiltration loss during the inflow process, reflecting the reduction characteristics of flood propagation in dry riverbeds, and used to improve the accuracy of downstream section flood volume forecasting. Using the effective water volume after deducting losses as the driving force, a discrete confluence model is used to simulate the spatiotemporal evolution of floods in the river network. Flood flow process forecasting results include the flow sequence values ​​of the forecast section over time, from which hydrological characteristic indicators such as peak flow, total flood volume, and peak occurrence time are extracted. By calculating the energy loss and infiltration loss during the inflow process, the reduction phenomenon of flood propagation in dry riverbeds is reflected, improving the accuracy of downstream section flood volume forecasting.

[0040] Based on the above embodiments, the construction of effective rainfall boundary conditions is further explained. In one possible implementation, the effective rainfall for a given period is calculated based on hydrological and meteorological observation data. This process specifically includes:

[0041] Step 201: Extract multi-point rainfall sequences and meteorological evaporation elements from hydrological and meteorological observation data.

[0042] Multiple hydrological and meteorological stations were deployed within the target watershed to simultaneously collect raw datasets, namely hydrological and meteorological observation data. From this dataset, multi-site rainfall sequences characterizing precipitation processes and meteorological evaporation elements for evaporation calculations were extracted. The multi-site rainfall sequences included precipitation depth records from each rain gauge within a set time step, the duration of which was predetermined. Meteorological evaporation elements included at least the following meteorological parameters: air temperature, net radiation flux, relative humidity, wind speed, and soil heat flux. Through this data extraction, the mixed raw observation data could be decoupled into independent meteorological inputs driving the hydrological model.

[0043] Step 202: Spatial weights are used to spatially weight the multi-point rainfall sequence to obtain the average rainfall over the watershed.

[0044] Since rainfall records from a single rain gauge only reflect precipitation at a local location, they need to be expanded to represent the average rainfall across the entire watershed. In this embodiment, based on the geographic coordinates of each rain gauge, the Thiessen polygon method is used to assign a corresponding control area proportion, i.e., a pre-configured spatial weight, to each rain gauge. The time-period rainfall at each station is multiplied by its corresponding spatial weight, and then summed to obtain the watershed average rainfall. This spatial weighting process can suppress the interference of local rainfall extremes and reflect the macroscopic distribution characteristics of rainfall at the watershed scale.

[0045] In other alternative implementations, if the watershed topography is highly undulating, the pre-configured spatial weights can also be pre-calculated using the inverse square distance method that takes into account the influence of elevation, or the Kriging interpolation method.

[0046] Step 203: Calculate the potential evapotranspiration based on meteorological evaporation factors, and reduce the potential evapotranspiration by combining the soil moisture stress state characterized by the underlying surface data to obtain the actual evapotranspiration.

[0047] Specifically, the Priestley-Taylor formula is used, with net radiation, air temperature, and soil heat flux from the meteorological evaporation elements input, to calculate the potential evapotranspiration under conditions of unrestricted water supply. Based on the actual water content of the current soil layer, the potential evapotranspiration is converted into actual evapotranspiration. Due to the vertical differences in soil moisture distribution, this embodiment employs a stratified calculation logic: evaporation demand is subtracted from the water storage of the upper soil layer; when the upper layer water is insufficient, the remaining evaporation demand is transferred to the lower layers and reduced using the lower soil moisture stress coefficient; if the lower layer water is still insufficient, the demand continues to be transferred to deeper layers and reduced using the deep soil moisture stress coefficient. The sum of the reduced evaporations from each layer is the actual evapotranspiration. This embodiment reflects the resistance of deep soil to the evaporation process during drought periods.

[0048] Step 204: Determine the effective rainfall for the time period based on the difference between the average rainfall and the actual evapotranspiration in the watershed.

[0049] After obtaining the average rainfall and actual evapotranspiration over the watershed, the difference between the two is calculated. When the rainfall is greater than the actual evapotranspiration, the difference is positive, indicating that there is net water input into the watershed. This difference is taken as the effective rainfall for the current time period and included in runoff calculation. When the rainfall is less than or equal to the actual evapotranspiration, it indicates that the precipitation has been consumed by evaporation, and the effective rainfall for the time period is set to zero. Before calculating the difference, the input rainfall and evapotranspiration must be consistent in terms of time scale and numerical units.

[0050] Based on the overall forecasting process provided in the above embodiments, this embodiment further explains the dual-overflow differentiation mechanism and its spatial probability allocation. Here, dual-overflow refers to the over-permeability flow mechanism and the saturation flow mechanism.

[0051] In one possible implementation, excess infiltration runoff is determined based on effective rainfall over a given time period and underlying surface data to obtain excess surface runoff and actual total infiltration, including:

[0052] Step 301: Analyze the rainfall intensity corresponding to the effective rainfall amount during the time period, and calculate the ratio of rainfall intensity to the current infiltration capacity to obtain the water supply degree. When the water supply degree indicates that infiltration overflow has occurred, construct a normalized distribution feature of infiltration capacity reflecting spatial heterogeneity based on the underlying surface data. Combine the water supply degree with the normalized distribution feature of infiltration capacity for spatial probability integration. Alternatively, construct the effective water supply degree based on the effective rainfall amount during the time period and the maximum cumulative infiltration capacity of the hypothetical unit, and combine the effective water supply degree with the normalized distribution feature of infiltration capacity for spatial probability integration.

[0053] The average infiltration rate of the watershed is calculated; the actual total infiltration is determined based on the average infiltration rate of the watershed, and the excess water volume in the effective rainfall of the period that exceeds the actual total infiltration is taken as the surface runoff of the infiltrate; when it is determined that no infiltrate runoff has occurred in the water supply indicator, the entire effective rainfall of the period is taken as the actual total infiltration, and the surface runoff of the infiltrate is taken as zero.

[0054] In this embodiment, the analytical operation refers to converting the effective rainfall within a time interval into an intensity value per unit time. Water supply is a parameter used to determine whether surface runoff exceeds infiltration. Current infiltration capacity refers to the upper limit of the soil matrix's absorption rate at the current moment, taking into account the effects of previous soil wetting.

[0055] Specifically, the infiltration rate is calculated using a modified single-point infiltration formula, which is related to the soil water content before infiltration, soil pore size distribution parameters, and saturated hydraulic conductivity. After determining that excessive infiltration has occurred, this embodiment introduces the concept of a fictitious unit to address the spatial heterogeneity of infiltration capacity within the watershed. This divides the watershed into a set of countless micro-elements with different infiltration characteristics. The normalized distribution characteristic of infiltration capacity is characterized by establishing a functional relationship between relative infiltration capacity and cumulative probability, thus representing the statistical distribution law of infiltration capacity at each point within the watershed.

[0056] When performing spatial probability integral processing, the average watershed water supply is substituted into the integral function of the normalized distribution curve of infiltration capacity to calculate the average watershed infiltration rate. That is, the coverage area of ​​the distribution curve within the water supply range is numerically integrated.

[0057] In one alternative implementation, the expression for the normalized distribution curve of infiltration capacity is:

[0058] G(f / f _m )=1-(1-(f / f _m -f _min / f _m ) / (1-f _min / f _m ) 1 / (1+b) ;

[0059] Among them, f _min / f _m Let f be the minimum relative infiltration capacity within the watershed, b be the shape parameter, and f be the infiltration capacity. _m represents the maximum infiltration capacity of a hypothetical unit. f represents the actual infiltration capacity of a single point or unit, which is the instantaneous infiltration capacity variable at a certain location within the watershed.

[0060] Define the effective water supply quantity X as the effective rainfall P during the current period. _net Maximum cumulative infiltration capacity Δ of the virtual unit _FmThe ratio of Δ to infiltration rate, which differs from the ratio of rainfall intensity to infiltration rate used to determine infiltration overshoot, are independent parameters used in different calculation stages. This yields the actual infiltration ratio, reflecting the average level across the entire watershed. Wherein, Δ... _Fm =f _m ×Δt; Δt is the calculation time step length.

[0061] The actual total infiltration is obtained by multiplying the average infiltration rate of the watershed by the maximum infiltration capacity of the hypothetical unit. The portion of the effective rainfall during a given period that fails to infiltrate is converted into excess surface runoff. The corresponding calculation logic is as follows:

[0062] X=P _net / Δ _Fm ;

[0063] ΔF _0 =η(X)×Δ _Fm ;

[0064] ΔR _s =P _net -ΔF _0 ;

[0065] Where X represents the effective water supply degree, and P _net Δ represents the effective rainfall amount for that period of time. _Fm Let η(X) be the maximum cumulative infiltration capacity of the hypothetical unit, and ΔF be the average infiltration rate function of the watershed. _0 ΔR represents the actual total infiltration. _s This refers to the permeable surface runoff.

[0066] Step 302: During the runoff determination process, the infiltration excess surface runoff and the matrix saturation runoff component dynamically switch or coexist within the same calculation period in the following ways: The infiltration excess runoff mechanism is triggered in real time based on the magnitude relationship between rainfall intensity obtained from the analysis of effective rainfall during the period and the infiltration capacity represented by the underlying surface data; the saturation runoff mechanism is triggered in real time based on the soil tension water saturation state determined by the current soil water storage, generating the matrix saturation runoff component; when the rainfall intensity is greater than the infiltration capacity, the infiltration excess runoff mechanism dominates the runoff generation; when the rainfall intensity is less than or equal to the infiltration capacity and the soil water storage reaches its maximum storage capacity, the saturation runoff mechanism dominates the runoff generation. That is, when heavy rainfall occurs and the soil moisture content is low, the infiltration excess runoff mechanism dominates the runoff generation; when moderate to light rainfall occurs and the soil moisture content is high, the saturation runoff mechanism dominates the runoff generation.

[0067] In this embodiment, dynamic switching or coexistence means that within the same calculation step, there can be surface runoff due to excessive rainfall intensity and runoff due to soil saturation within the watershed.

[0068] Specifically, the rainfall intensity and current soil infiltration capacity are monitored in real time. When the rainfall intensity exceeds the infiltration capacity, the excess surface runoff is calculated. If the tension water storage capacity of the soil layer has reached its upper limit, the newly infiltrated water is no longer retained by the soil pores, thus triggering the saturation runoff mechanism and generating a matrix saturation runoff component. Through this two-way determination and dynamic switching, the runoff component is deconstructed to provide a realistic flow composition for runoff calculation.

[0069] Based on the above embodiments, this embodiment further explains the determination of the dynamic bypass coefficient and the allocation of infiltration volume. In one possible implementation, based on the effective rainfall over a period of time and a dynamic bypass mechanism determined by hydrological and meteorological observation data, the actual total infiltration volume is decomposed into matrix domain infiltration volume and preferential inflow infiltration volume, including the following steps:

[0070] Step 401: Based on the rainfall intensity of the current period obtained from the analysis of effective rainfall in the current period, the underlying surface data and the soil moisture state characterized by the current soil water storage, and the macropore development state determined based on hydrological and meteorological observation data, determine the dynamic bypass coefficient.

[0071] Because the soil layer is not a homogeneous structure, water entering the soil follows different transport paths. In this embodiment, the system acquires rainfall dynamic characteristics in real time and assesses the current water-holding capacity of the soil. Accordingly, a dynamic development index of macropores is introduced, and a dynamic bypass coefficient is calculated through multi-factor coupling. This coefficient is a dynamic distribution ratio between 0 and an upper limit constraint value, used to replace the fixed distribution constant in traditional models, providing a high-resolution numerical basis for the diversion of infiltrated water.

[0072] Step 401a: Extract the duration of continuous no effective rainfall before the occurrence of a rainfall event from hydrological and meteorological observation data.

[0073] In this step, duration is a timescale variable used to quantify the previous drought state. By tracing historical meteorological observation data of the target watershed, starting from the beginning of the current rainfall event and searching backwards along the historical timeline, the number of days with daily rainfall consistently below a preset invalid rainfall threshold is counted. This threshold indicates that when daily rainfall falls below this threshold, moisture is rapidly lost, making it difficult to increase soil moisture content, leading to a continuous accumulation of drought days. In this embodiment, duration is extracted to obtain long-range memory features of soil drought history.

[0074] Step 401b: Construct an early drought-induced crack development factor based on the duration and a preset crack development characteristic time scale to characterize the internal crack development state of the soil caused by long-term dehydration and shrinkage.

[0075] Early-stage drought fissure development factors are used to characterize changes in soil physical structure caused by alternating wet and dry conditions. For example, when clay or loam soils in arid or semi-arid regions lack water replenishment for a long period, they will undergo volume shrinkage, forming a network of drying shrinkage fissures within the soil. These fissures constitute physical channels for the rapid migration of infiltrated water.

[0076] Accordingly, this embodiment establishes a mapping relationship between the number of consecutive rainless days and this factor, transforming meteorological variables into parameters characterizing the connectivity of micropores. The corresponding calculation relationship is as follows:

[0077] Ψ _Tdry =Ψ _min +(1-Ψ _min )×(1-exp(-T _dry / T _c ));

[0078] Among them, Ψ _Tdry Ψ is a factor contributing to the development of early-stage drought-induced fissures. _min Let T be the minimum fracture development factor, exp be an exponential function with the natural constant as the base, and T be the minimum fracture development factor. _dry T represents the duration of the period without effective rainfall in the preceding period. _c This represents the preset timescale for fracture development characteristics. The fracture development factor during the early stages of drought exhibits a non-linear positive correlation with the duration of drought. That is, as the drought duration increases, the negative exponential decay term approaches zero, and the factor shows a non-linear monotonically increasing trend, converging to a saturated upper limit.

[0079] In this embodiment, the minimum fracture development factor is used to characterize the bypass capacity of a permanent macroporous substrate unaffected by wet-dry cycles. The fracture development characteristic timescale reflects the typical time required for soil to evolve from a saturated wet state to fully fractured development.

[0080] Step 401c: The dynamic bypass coefficient is calculated by combining the early drought fissure development factor as the modulation weight with the rainfall intensity and soil moisture status obtained from the analysis of effective rainfall during the time period.

[0081] The proportion of water entering the soil that can bypass the matrix pores is called the dynamic bypass coefficient. In this embodiment, a combination operator is used to map the instantaneous driving force of rainfall, the soil water holding capacity, and the degree of pore cracking. The corresponding calculation logic is as follows:

[0082] β=β _max ×(1-exp(-k _p ×(P-P _0 ) / I _s ))×(W / W _M ) σ ×Ψ _Tdry P>P _0 ;

[0083] β=0, P≤P _0 ;

[0084] Where β is the dynamic bypass coefficient, β _max The maximum bypass coefficient is the upper limit, exp is an exponential function with the natural constant as the base, and k is the number of decimal places. _p P is the rainfall intensity sensitivity coefficient, where P is the rainfall intensity. _0 The rainfall intensity parameter, i.e. the effective rainfall threshold, is used to determine the priority flow. _s The baseline permeability intensity for the watershed reflects the response threshold of the soil matrix to rainfall input. Its value is on the same order of magnitude as saturated hydraulic conductivity and is determined through calibration as a pre-configured parameter for the model. W represents the current soil water storage. _M The maximum tensile water storage capacity is σ, which is the soil wetting state response parameter. Its value can be determined by combining it with other model characteristic parameters based on the soil texture characteristics and historical flood data of the target watershed.

[0085] This calculation limits the maximum bypass coefficient to the maximum flow rate of the preferential flow channel. Rainfall intensity driving term (1 - exp(-k) _p ×(P-P _0 ) / I _s This indicates that the greater the rainfall intensity, the higher the bypass allocation ratio. Soil moisture state term (W / W) _M ) σ This indicates that under sandy loam conditions, as the matrix gradually becomes saturated, newly infiltrated water is squeezed into the macropores, and the two are positively correlated. The pre-drought fissure development factor, as a multiplicative modulating weight, controls the baseline level of the bypass coefficient.

[0086] In some alternative implementations, for clay strata with expansive properties, the rapid volume expansion after water absorption leads to fissure closure, and the promoting effect of soil wetting state on preferential flow channels exhibits a reversed characteristic. The wetting state term exponent σ can be modified to a negative value, making the dynamic bypass coefficient exhibit a non-linear negative correlation with the current soil water storage, thereby adapting to the physical-hydrological response characteristics of this soil composition. When σ < 0, the bypass coefficient still needs to satisfy β ≤ β _max The upper limit constraint is used to avoid numerical overflow under extremely dry soil conditions.

[0087] Furthermore, the rainfall intensity sensitivity coefficient k _p Response parameters σ of soil moisture state in the early stage and time scale T of crack development characteristics _c It can be determined based on historical flood data through subsequent multi-objective calibration methods. During the calibration process, k should be ensured... _p It satisfies the physical constraint that the bypass coefficient monotonically increases with increasing rainfall intensity, and σ satisfies the sign consistency constraint that the runoff response changes with soil moisture.

[0088] Furthermore, when calculating the dynamic bypass coefficient, the initiation of the priority flow is constrained by combining a preset effective rainfall threshold.

[0089] In step 401d, when the rainfall intensity is less than or equal to the effective rainfall threshold, the dynamic bypass coefficient is set to zero, and the actual total infiltration is entirely classified as matrix domain infiltration.

[0090] Specifically, the minimum kinetic energy threshold required to overcome the suction of the soil surface granular matrix is ​​defined as the effective rainfall threshold. When the rainfall intensity does not exceed this threshold, infiltrated water, dominated by the suction gradient, is intercepted by surface pores and struggles to accumulate into gravity flow entering the vertical macropores. Therefore, the dynamic bypass coefficient is set to 0, cutting off the computational branch leading to the preferred flow channel. This logic avoids computational bias caused by erroneously activating bypass channels under low-intensity rainfall conditions.

[0091] Step 401e: When the rainfall intensity is greater than the effective rainfall threshold, the dynamic bypass coefficient calculated by the modulation weight satisfies the preset maximum bypass upper limit constraint.

[0092] When rainfall intensity exceeds this threshold, the preferential flow channel is activated. Due to the physical extremities of the volume and instantaneous hydraulic conductivity of macropores, the bypass allocation ratio cannot be unlimited. After completing the multi-factor coupling calculation, the system needs to perform a truncation check. If the calculated coefficient value exceeds the maximum bypass upper limit, it is replaced with this upper limit value for calculations under stable extremely heavy rainfall conditions. The maximum bypass upper limit is used to characterize the ultimate hydraulic conductivity of regional tubular pores.

[0093] Step 402: The actual total infiltration volume is proportionally allocated using the dynamic bypass coefficient. The water volume allocated to the macroporous channels is taken as the priority infiltration volume, and the remaining infiltration volume is taken as the matrix domain infiltration volume.

[0094] After obtaining the dynamic bypass coefficient after boundary constraints, multiply it by the actual total infiltration. The product is the preferential inflow infiltration volume into the deep fracture network. Subtract this preferential inflow infiltration volume from the actual total infiltration volume, and the difference is the matrix domain infiltration volume that slowly seeps along the gaps between conventional soil particles.

[0095] Based on the above embodiments, this embodiment further explains the vertical distribution of moisture in the matrix domain and the correction of water balance. In one possible implementation, the infiltration water in the matrix domain is treated by layer-by-layer soil filling to obtain the matrix-saturated runoff component, including the following steps:

[0096] The system obtains the current water storage capacity of each soil layer maintained in real time. The initial water storage capacity is determined by the soil's early water content in the underlying surface data. That is, the current water storage capacity of each soil layer is obtained, and the difference between the current water storage capacity and the maximum water storage capacity of each soil layer is calculated in the order of upper soil layer, lower soil layer, and deep soil layer. This difference is used as the recharge capacity of each layer. The infiltration water in the matrix domain is used to fill the recharge capacity layer by layer in sequence, and the current water storage capacity of each soil layer is updated accordingly and saved to the system status for use in the next calculation period. The remaining water after the upper soil layer, lower soil layer, and deep soil layer are filled is input into the preset free water storage module, and the matrix full-storage flow component is calculated according to the preset outflow distribution rules.

[0097] In this embodiment, the layered soil filling treatment is used to simulate the process of rainfall infiltration water gradually wetting the soil from the surface to deeper layers. The infiltration water in the matrix domain first enters the upper soil water storage layer. Correspondingly, the current water storage capacity of the upper soil is extracted and subtracted from the preset maximum water storage capacity of the upper layer to calculate the amount of water the upper soil can further absorb, i.e., the recharge capacity. The corresponding aboveground recharge amount is:

[0098] Δ _WU =min(F _matrix ,WUM-WU _t );

[0099] Where, Δ _WU F represents the actual replenishment amount of the upper soil layer. _matrix WU represents the infiltration volume of the matrix domain, WUM represents the maximum water storage capacity of the upper layer, and WU represents the maximum water storage capacity of the upper layer. _t is the current water storage capacity of the upper soil layer at the start of the current step, and min is the function for finding the minimum value.

[0100] After completing the upper layer recharge, update the upper soil moisture content and calculate the amount of water that has migrated to subsequent layers:

[0101] WU _next =WU _t +Δ _WU ;

[0102] P _net1 =F _matrix -Δ _WU ;

[0103] Among them, WU _next For the updated water storage capacity of the upper soil, P _net1 This represents the remaining water volume after replenishment from the upper layer.

[0104] Based on this, the remaining water volume P _net1 Following the same logic, the water enters the lower soil water storage layer and the deeper soil water storage layer in sequence for replenishment.

[0105] In this embodiment, a closed-loop correction for water balance also needs to be considered. Based on physical mechanisms, some preferential flow may be absorbed by the surrounding dry matrix during its penetration through the thick vadose zone. Therefore, a correction term from the preferential flow channel is introduced into the water storage update of deep soil layers, and the corresponding formula is:

[0106] WD _next =min(WD _t +ΔWD _t +R _retain WDM);

[0107] R _retain =F _pref ×(1-κ);

[0108] Among them, WD _next For the updated deep soil water storage, WD _t ∆WD represents the current water storage in the deep soil at the start of the current step. _t R represents the deep replenishment amount from the layer-by-layer filling of the matrix domain. _retain WDM represents the deep recharge volume that does not penetrate to the groundwater surface, and is the maximum deep water storage capacity.

[0109] When the infiltration water in the matrix zone is replenished by the upper, lower, and deep soil layers, and the deep recharge is taken into account, if there is still excess water exceeding the soil's water storage capacity, this portion of the water will be converted into gravity flow, i.e.:

[0110] P _net3 =max(0,WD _t +ΔWD _t +R _retain -WDM) + Existing excess water volume in each layer before replenishment;

[0111] Among them, P _net3 This represents the total amount of water entering the free water storage module.

[0112] Accordingly, the free water storage module, based on a preset outflow distribution rule—that is, using pre-configured lateral outflow coefficients and deep permeability coefficients—divides the total water volume into an interflow component and a component of the model's original groundwater. These components aggregate to form the matrix-filled runoff component. The free water storage module adopts a linear reservoir outflow model, whereby after the surplus water is input into the free water storage reservoir, it is proportionally distributed according to pre-configured surface runoff outflow coefficients and groundwater runoff outflow coefficients to obtain the interflow and groundwater runoff components. This can be implemented with reference to the standard outflow distribution rule of the free water storage reservoir in the Xin'anjiang model. The sum of the interflow and groundwater runoff components constitutes the matrix-filled runoff component.

[0113] Based on the above embodiments, this embodiment provides a detailed description of the preferential flow penetration depth attenuation and water volume classification process based on groundwater depth. In one possible implementation, when determining the penetration attenuation characteristics based on the groundwater depth in the underlying surface data, and classifying the preferential inflow seepage volume into preferential groundwater recharge and deep recharge volume according to the penetration attenuation characteristics, the process includes:

[0114] Step 501: Obtain the representative groundwater burial depth and preferential flow characteristic penetration depth of the target area.

[0115] In this embodiment, groundwater depth refers to the vertical distance from the surface to the groundwater level. In arid and semi-arid regions, due to long-term groundwater extraction, this vertical distance often exhibits spatial heterogeneity. This groundwater depth can be obtained through regional hydrogeological survey data or collected in real-time by groundwater monitoring wells deployed within the target area.

[0116] The preferential flow characteristic penetration depth is a preset physical parameter used to characterize the average vertical depth of the effective extension of preferential flow channels within a watershed. This parameter reflects geological and ecological characteristics such as root canal length, the depth of desiccation fracture development, and macropore connectivity. In practical applications, the value of the preferential flow characteristic penetration depth depends on the local soil profile structure.

[0117] Step 502: Based on the ratio of groundwater burial depth to the penetration depth of preferential flow characteristics, a penetration attenuation feature with an exponentially decreasing law is constructed.

[0118] Penetration attenuation characteristics were used to quantify the proportion of water loss caused by the lateral suction of the unsaturated soil matrix during the vertical infiltration of the preferential flow. Physical experiments and tracer detection showed that the preferential flow channels are not absolutely closed, and water is continuously absorbed by the pore walls during transport, with the loss proportion increasing the longer the transport distance.

[0119] Therefore, this embodiment utilizes the ratio of groundwater burial depth to the penetration depth of preferential flow characteristics to establish the expression corresponding to the mathematical attenuation model:

[0120] κ=exp(-D _w / D _ref );

[0121] Where κ is the penetration attenuation characteristic, representing the percentage of water that can penetrate the vadose zone and reach the groundwater surface, and D _w D is the depth of groundwater burial. _ref The penetration depth is the characteristic of the priority flow.

[0122] It should be understood that when the groundwater depth approaches zero, the penetration attenuation characteristic approaches 1, indicating that water reaches the water surface with almost no loss; when the groundwater depth is greater than the penetration depth of the preferential flow characteristic, the characteristic value decays exponentially and approaches zero, indicating that water is absorbed by the vadose zone matrix during long-distance transport.

[0123] Step 503: The preferential inflow infiltration volume is reduced and extracted using penetration attenuation characteristics. The effective penetration volume reaching the groundwater surface after reduction is taken as the preferential groundwater recharge volume, and the remaining water absorbed during vertical transport is taken as the deep recharge volume. In other words, the remaining water absorbed by the vadose zone matrix during vertical transport that does not reach the groundwater surface is taken as the deep recharge volume. This is then overlaid and updated with the current deep soil water storage volume maintained in real time in the system, and the updated deep soil water storage volume is saved to the system status for subsequent runoff calculations.

[0124] In this embodiment, the separated macropore water flow is secondary-distributed using penetration attenuation characteristics. Specifically, the preferential inflow seepage volume is multiplied by the penetration attenuation characteristics to obtain the actual amount of water that crosses the vadose zone and flows into the groundwater reservoir, i.e., the preferential groundwater recharge volume. According to the law of conservation of mass, the preferential inflow seepage volume is subtracted from the preferential groundwater recharge volume to obtain the amount of water intercepted by the surrounding matrix, i.e., the deep recharge volume.

[0125] Accordingly, the deep recharge amount is input as a compensation term into the deep water storage state calculation module of the matrix domain, and is superimposed with the existing water content of the deep soil to update the current water storage of the deep soil. If the total water volume after superposition exceeds the maximum water holding capacity of the deep soil, the excess water will enter the pre-set free reservoir to participate in the lateral outflow distribution.

[0126] Furthermore, for some karst basins, such as those with extremely shallow groundwater levels or bedrock fissures directly connecting to groundwater, the absorption effect of the vadose zone matrix is ​​weak. Therefore, a fixed preferential flow groundwater recharge coefficient needs to be pre-configured to replace the exponential decay function. This coefficient is then multiplied by the preferential inflow seepage volume to obtain the preferential flow groundwater recharge amount.

[0127] For arid and semi-arid riverbeds characterized by prolonged dryness or water deficit, this embodiment improves the accuracy of flood evolution forecasting by introducing an initial infiltration enhancement mechanism. One possible implementation specifically includes:

[0128] Based on the initial infiltration characteristics of the riverbed in the underlying surface data, or based on historical river flow records in hydrological and meteorological observation data, the dynamic transmission loss of river inflow during its evolution is calculated. This includes: extracting the maximum potential infiltration capacity of the riverbed from the underlying surface data and calculating the ratio of river inflow to the maximum potential infiltration capacity of the riverbed to obtain the generalized water supply of the river; assessing the early wet state of the riverbed based on historical river flow records and water evaporation and drying characteristics, and updating the early wet state of the riverbed synchronously according to the continuous inflow of river inflow during the current calculation period; and performing spatial infiltration reduction calculation by combining the generalized water supply of the river and the dynamically updated early wet state of the riverbed to obtain the dynamic transmission loss of the river.

[0129] Specifically, the maximum potential infiltration capacity of a riverbed refers to the maximum volume of water that the riverbed material can accept per unit time. The corresponding calculation formula is:

[0130] f _max_ch =f _0_ch ×B _w ×L _r ;

[0131] Among them, f _max_ch The maximum potential infiltration capacity of the riverbed, expressed in meters (m). 3 / s;f _0_ch B represents the maximum infiltration rate of the riverbed material, expressed in m / s. _w The representative wet width of the river section, in meters (m); L _r The length of the river section is expressed in meters (m).

[0132] The generalized water supply ratio of a river channel is used to characterize the supply-demand ratio between the current amount of water entering the river channel and the riverbed's water absorption capacity. The corresponding calculation formula is:

[0133] X _riv =Q _in / f _max_ch ;

[0134] Among them, X _riv For the generalized water supply degree of the river channel, Q _in This represents the inflow rate into the river during the current period.

[0135] The initial wetting state of the riverbed is used to characterize the water saturation level of the surface material. Specifically, a preset riverbed wetting response coefficient and riverbed drying attenuation coefficient are obtained. Based on the continuous inflow of river flow during the current calculation period, this state is updated step-by-step. The dynamic update logic is as follows:

[0136] θ _ch (t+Δ _t )=θ _ch (t)+C _w ×(1-θ _ch(t))×(Q _in / (Q _in +Q _ref ))-C _d ×θ _ch (t)×(Q _ref / (Q _in +Q _ref ));

[0137] Where, θ _ch (t+Δ _t ) represents the updated early wet state of the riverbed, θ _ch (t) represents the current state of the riverbed's initial wetting condition, C _w C is the preset riverbed wetting response coefficient. _d Q is the preset riverbed drying attenuation coefficient. _in Q represents the inflow rate into the river channel. _ref The pre-configured reference flow can be the multi-year average flow or the base flow level of the river section. It can be determined based on the historical flow statistics of the target river section. Those skilled in the art can directly obtain it based on the statistical mean of the measured flow data of the river section hydrological station over many years or the base flow segmentation results.

[0138] This formula reflects the physical competition between the wetting and drying terms, that is, when the river flow Q... _in When the flow rate is much higher than the reference flow rate, the wetting term dominates and the wetting state tends to be 1; when the riverbed is dry or the flow rate is low, the drying term dominates and the wetting state tends to be 0. Through this mechanism, the continuous changes in the dryness and wetness of the riverbed can be tracked in real time.

[0139] To ensure the effectiveness of the initial moist state of the riverbed, C _w With C _d All should satisfy C _w ∈(0,1], and C _d ∈(0,1]; after each update calculation, if the result exceeds the interval [0,1], it is truncated to 0 or 1 respectively. C _w With C _d The constraint conditions are determined by combining multi-objective calibration methods with historical field data, and are used to make θ _ch It belongs to [0,1] during the simulation.

[0140] In this embodiment, when Q _in Much larger than Q _ref When the riverbed is dry or the flow is low, the drying term dominates and the wet state tends to be 1; when the riverbed is dry or the flow is low, the drying term dominates and the wet state tends to be 0.

[0141] In one optional implementation, spatial infiltration reduction calculations are performed by combining the generalized water supply of the river channel with the dynamically updated pre-wetting state of the riverbed to obtain the dynamic transport loss of the river channel, including:

[0142] Based on the degree of dryness characterized by the early wet state of the riverbed, an initial infiltration enhancement factor exhibiting a nonlinear decay law was constructed to quantify the rapid infiltration capacity driven by matrix suction when the dry riverbed is initially irrigated.

[0143] By incorporating the initial infiltration enhancement factor into the transmission loss calculation process, a dynamic reduction weight higher than the steady-state permeability is assigned to the riverbed in a dry state, thus obtaining the dynamic reduction weight for the current time period.

[0144] As the initial wet state of the riverbed gradually approaches saturation under the replenishment of water flow, the dynamic reduction weight is smoothly reduced with the increase of wetness by the initial infiltration enhancement factor until the steady state corresponding to the stable infiltration state of the riverbed is reached. Combined with the generalized water supply of the river channel, the final dynamic transmission loss of the river channel is calculated.

[0145] To address the significant initial water loss in dried-up river channels in arid and semi-arid plains, this embodiment introduces a nonlinear enhancement mechanism, with the initial infiltration enhancement factor formulated as follows:

[0146] ξ _θ =1+λ _s ×(1-θ _ch ) n _s ;

[0147] Where, ξ _θ λ is the initial permeation enhancement factor. _s λ is a preset initial absorption enhancement coefficient used to characterize the matrix suction strength of dry riverbed materials. _s The constraint condition is λ _s ≥0, i.e., ξ _θ ≥1 ensures that the infiltration capacity of the dry riverbed is not lower than the steady-state level. Its value can be determined using the multi-objective calibration method of this invention, utilizing measured transmission loss data from the first flood passage period in historical data. θ _ch The riverbed was in its early, moist state; n _s The preset absorption attenuation index is used to characterize the rate attenuation of the control enhancement effect as the degree of wetting changes.

[0148] Accordingly, after obtaining the enhancement factor, it is incorporated into the reduction operator of river transport loss. When the riverbed is in an extremely dry state, i.e., θ... _ch When ξ approaches 0 _θ The maximum value is taken, and the initial reduction weight is amplified. As the flood progresses and the riverbed gradually becomes wet, the enhancement factor quickly drops back to 1, and the transmission loss returns to a steady-state level.

[0149] In one alternative implementation, the river dynamic transport loss Q _loss The expression is:

[0150] Q _loss =η _ch (X _riv ,θ _ch )×ξ _θ ×f _max_ch ;

[0151] Where, η _ch It is the infiltration reduction function of the riverbed space, and its form is similar to the average infiltration rate function η(X) of the watershed.

[0152] ξ _θ As a dynamic amplification term of the maximum potential infiltration capacity of the riverbed, the effective upper limit of riverbed infiltration f for the current period is obtained. _eff_ch The calculation formula is as follows:

[0153] f _eff_ch =ξ _θ ×f _max_ch ;

[0154] The corresponding dynamic transport loss in the river channel is:

[0155] Q _loss =η _riv (X _riv * )×f _eff_ch ;

[0156] X _riv * =Q _in / f _eff_ch ;

[0157] Among them, X _riv * For the corrected river water supply; when Q _loss ≥Q _in When, take Q _loss =Q _in The effective inflow into the river channel is zero.

[0158] Based on the river loss calculation results obtained from the above embodiments, this embodiment further explains the implementation process of river confluence calculation and the output of flood forecast results. In one possible implementation, the river confluence calculation is performed after deducting the river dynamic transmission loss from the river inflow, and the flood flow process forecast results for the forecast section are output, including:

[0159] Step 601: The dynamic transport loss of the river channel is used as the source-sink dissipation term in the river continuity equation and deducted from the river inflow to obtain the effective river inflow.

[0160] Specifically, river inflow refers to the total volume of water that converges into the river cross-section after surface runoff, interflow, and corrections. Effective river inflow refers to the net flow component that actually participates in downstream propagation calculations, taking into account riverbed infiltration. Correspondingly, the dynamic transport loss of the river is treated as a negative source term in the river water balance equation, and water dissipation is handled through subtraction. The corresponding formula is:

[0161] Q _e =Q-Q _loss ;

[0162] Among them, Q _e The effective inflow rate is Q, where Q is the inflow rate of the river during the current time period. _loss This is the calculated dynamic transmission loss in the river channel. During the continuous evolution of river floods, if the loss in the current period is greater than or equal to the inflow, the effective river inflow is zero, indicating that the flood is intercepted or absorbed in this section of the river and is difficult to propagate downstream.

[0163] Step 602: Based on the hydrodynamic storage and discharge characteristics parameters of the target river section, the discrete confluence method is used to perform time-by-time flow propagation calculations on the effective river inflow to deduce the flood flow process forecast results of the forecast section.

[0164] The target river section refers to the river channel interval from the runoff inflow point to the predicted control section. Hydrodynamic storage and discharge characteristics parameters include storage and discharge parameters reflecting the river channel's flood channel storage capacity, and weighting coefficients reflecting the distribution of discharge weights. The Muskingan method can be used for calculations in the discrete runoff confluence method.

[0165] Based on the calculation time step, storage and release parameters, and weighting coefficients, the Muskingan calculation coefficients are pre-calculated, namely:

[0166] C _0 =(0.5×Δ _t -K×ε) / (K-K×ε+0.5×Δ _t );

[0167] C _1 =(0.5×Δ _t +K×ε) / (K-K×ε+0.5×Δ _t );

[0168] C _2 =(K-K×ε-0.5×Δ _t ) / (K-K×ε+0.5×Δ _t );

[0169] Among them, C _0 C _1 C _2 Calculate the coefficients for Muskingen, Δ_t For the calculation of the time step, K is the channel storage and discharge parameter, and ε is the weighting coefficient. Using this coefficient, a time-series calculation of the effective channel inflow is performed to obtain the predicted outlet flow at the cross-section. The corresponding formula is:

[0170] O _t +1=C _0 ×I _e,t +1+C _1 ×I _e,t +C _2 ×O _t ;

[0171] Among them, O _t +1 represents the outflow rate of the forecast section at the current moment, O _t For the outflow rate of the cross-section predicted at the previous moment, I _e,t +1 represents the current effective river inflow, I _e,t This represents the effective inflow rate into the river channel at the previous moment.

[0172] Repeating this calculation process for each calculation period yields a complete flood flow forecast. Based on this flow process curve, key flood indicators for the forecast section are extracted, and the corresponding calculation logic is as follows:

[0173] W _f =∑(Q _t ×Δ _t );

[0174] Q _max =max(Q _t );

[0175] T _p =argmax _t (Q _t );

[0176] Among them, W _f To predict the total flood volume at the cross-section, Q _max For peak flow, T _p For peak time, Q _t This is the time series value of the flow at the forecast section. This embodiment introduces the confluence calculation after dynamic transmission loss, which can avoid the problems of early peak time and excessive flood volume that occur when traditional confluence methods are applied to rivers, thus improving the rationality of flood forecasting.

[0177] Based on the above embodiments, this embodiment further explains the offline calibration process of pre-configured model feature parameters involved in flood forecasting methods. The pre-configured model feature parameters required for the first-step dynamic bypass mechanism, penetration attenuation characteristics, and initial infiltration characteristics of the riverbed are obtained through the following offline construction steps:

[0178] A sample of historical typical flood events in the target watershed is obtained, which includes historical meteorological input data and actual observed flood flow events. Forecast calculations are performed using the model feature parameters to be calibrated and the historical meteorological input data to obtain the historical simulated flood flow events. A multi-objective optimization function is constructed based on the differences between the actual observed flood flow events and the historical simulated flood flow events. The multi-objective optimization function is used to iteratively optimize the model feature parameters to be calibrated until the preset convergence conditions are met. The final parameter combination is extracted as the pre-configured model feature parameters for use in real-time flood forecasting.

[0179] In this embodiment, the acquisition of historical typical flood process samples is the data basis for parameter calibration.

[0180] For example, taking a watershed in the western piedmont transition zone as an example, the catchment area of ​​this watershed is over 2000 square kilometers. Hourly rainfall observation records from 12 rain gauge stations within the watershed are collected and used as historical meteorological input data. Hourly measured flow data from the control section of one of the hydrological stations are obtained as the actual observed flood flow process. In this embodiment, 10 typical flood events from a certain period are selected as samples, and they are divided into calibration phase samples and verification phase samples.

[0181] The model characteristic parameters to be calibrated include the upper limit of the maximum bypass coefficient, the rainfall intensity sensitivity coefficient, the benchmark reference infiltration intensity, the rainfall intensity parameter for preferential flow initiation, the time scale of fracture development characteristics, the penetration depth of preferential flow characteristics, the riverbed wetting response coefficient, the riverbed drying attenuation coefficient, the initial infiltration enhancement coefficient, and the infiltration attenuation index.

[0182] During forecast calculations, historical meteorological data is input into the constructed hydrological model, and runoff generation and concentration are extrapolated using current initial parameter values ​​to obtain a continuously distributed historical simulated flood flow process along the time axis. By comparing the deviations between the simulated and measured flow processes in terms of peak flow, flood volume, and phase, the discrepancies are extracted.

[0183] In this embodiment, the multi-objective optimization function is used to evaluate the Nash efficiency coefficient, relative error of flood peak, relative error of total flood volume, and peak occurrence time error of the model forecast. The peak occurrence time error is expressed in a normalized form, i.e.:

[0184] dt _norm =∣T _p , _sim -T _p , _obs | / T _ref ;

[0185] Among them, T _ref For reference time scales, the duration of a flood event or the number of calculation steps can be used, and the formula is as follows:

[0186] F _obj =w _1 ×(1-NSE)+w _2 ×RE _q +w _3 ×RE _w +w _4 ×dt _norm ;

[0187] Among them, F _obj For the multi-objective optimization function value, w _1 w _2 w _3 w _4 The weighting coefficients are preset, NSE is the Nash efficiency coefficient, used to reflect the degree of fit between the simulated and measured flow processes, and RE is the weighting coefficient. _q RE represents the relative error of the peak flow rate. _w The relative error of the total flood volume is dt _norm This represents the peak occurrence time error.

[0188] The iterative optimization process can employ a hybrid evolutionary algorithm or a genetic algorithm. By continuously adjusting the model feature parameters to be calibrated, the calculation is repeated to determine the optimal function value. Iteration stops when the optimal function value reaches the global minimum or the rate of change is less than a preset convergence threshold. The final parameter combination obtained is the pre-configured model feature parameter. To ensure the robustness of the forecast, in application, these parameters are limited to a reasonable physical range.

[0189] According to another aspect of this application, this application also discloses another specific implementation of flood forecasting applicable to drought and semi-arid regions. It is understood that, based on the above embodiments, some implementation processes in this embodiment directly adopt the implementation methods of the above embodiments, and will not be specifically described in subsequent embodiments.

[0190] Specifically, hourly rainfall data from multiple rain gauge stations were selected, and the Thiessen polygon method was used to calculate the average rainfall P over the watershed. _average The expression is as follows:

[0191] P _average =∑ _(j=1) n (ω _j ×P _j );

[0192] ∑ _(j=1) n ω _j =1;

[0193] Among them, P _j Let ω be the rainfall at the j-th rain gauge during the specified time period. _jThe corresponding weights are given, and n represents the total number of rain gauge stations. The potential evapotranspiration E is calculated using the Priestley-Taylor formula. _M :

[0194] E _M =α _PT ×Δ / (Δ+γ)(R _n -G);

[0195] Where, α _PT The coefficients for the Priestley-Taylor formula can be selected based on the experience of those skilled in the art; Δ is the slope of the saturated vapor pressure as a function of temperature, γ is the hygrometer constant, and R... _n G represents net radiation, which is the net solar radiation received by the Earth's surface, and G represents soil heat flux.

[0196] Accordingly, the actual evapotranspiration E is calculated using the water stress function. _P :

[0197] E _P =K _C ×E _M ;

[0198] Among them, K _C This is due to water stress. According to P... _average With E _P Calculate the effective rainfall P for each time period. _net .

[0199] Furthermore, the water supply specificity X is defined as the ratio of rainfall intensity P to the current infiltration capacity f, specifically:

[0200] X = P / f;

[0201] When X > 1, it is determined that infiltration excess runoff has occurred in the current time period; when X ≤ 1, it is determined that rainfall in the current time period preferentially infiltrates and enters the soil water distribution process. The pre-infiltration soil water content B is defined. _0 for:

[0202] B _0 =θ _0 / Φ;

[0203] Where, θ _0 Where t is the initial water content and Φ is the porosity. Single-point infiltration rate f(t, B) _0 ) and cumulative infiltration F(t,B _0 They are respectively:

[0204] f(t,B _0 ) = 1 / 2 × S _r (1-B _0 (c+1) )×t (－1 / 2) +K_s (1-B _0 (2c+1) );

[0205] F(t,B _0 )=S _r ×(1-B _0 (c+1) )×t (1 / 2) +K _s (1-B _0 (2c+1) )×t;

[0206] Among them, S _r To ensure the macroscopic absorption rate of the soil is fully dried, c represents the soil pore size distribution parameter, and K... _s Let t be the saturated hydraulic conductivity, and t be the corresponding time.

[0207] Considering the complex slopes of the watershed and the spatial differences in infiltration capacity, this embodiment introduces a hypothetical maximum infiltration capacity f of a unit. _m And its cumulative value, the corresponding formula is:

[0208] f _m (t,(B _0 *)=1 / 2×S _r (1-B _0 *) (c+1) )×t (－1 / 2) +K _s ×(1-(B _0 *) (2c+1) );

[0209] F _m (t,(B _0 *)=S _r ×(1-(B _0 *) (c+1) )×t (1 / 2) +K _s ×(1-(B _0 *) (2c+1) )t;

[0210] Among them, B _0 * Represents the average pre-infiltration soil moisture content of the watershed. Relative infiltration capacity α is defined. _i for:

[0211] α _i =f _i / f _m ;

[0212] Among them, f _i This represents the actual infiltration capacity of the i-th unit.

[0213] Establish a normalized distribution curve for infiltration capacity, specifically as follows:

[0214] β(α)=0, α≤α _0 ;

[0215] β(α)=1-((1-α) / (1-α _0 )) b α _0 <α≤1;

[0216] Where α represents relative infiltration capacity, and its meaning is equivalent to α _i ;α _0 This refers to the minimum relative infiltration capacity within the watershed, which is equivalent to the aforementioned f. _min / f _m , where b is the shape parameter. Based on the effective water supply X at the watershed scale, the average infiltration rate function η(X) of the watershed is obtained:

[0217] η(X)=X, X≤α _0 ;

[0218] η(X)=α _0 +(1-α _0 ) / (b+1)[1-((1-X) / (1-α _0 )) (b+1) ], α _0 <X<1;

[0219] η(X)=α _0 +(1-α _0 ) / (b+1), X≥1;

[0220] The actual total infiltration amount ΔF during the calculation period _0 and excess permeable surface runoff ΔR _s :

[0221] ΔF _0 =η(X)×ΔF _m ;

[0222] ΔR _s =P _net -ΔF _0 ;

[0223] In this embodiment, the soil water structure still adopts the original three-layer soil water storage structure of the Xin'anjiang model, including an upper soil water storage layer, a lower soil water storage layer, and a deep soil water storage layer. The infiltration water volume in the matrix domain is replenished layer by layer according to the three soil layers of the Xin'anjiang model; the preferential inflow infiltration water volume does not participate in the layer-by-layer filling process of the three soil layers, bypasses the three soil water storage structure, and is input into the groundwater module as groundwater recharge.

[0224] Alternatively, the preferential inflow of seepage water is reduced based on the groundwater depth and the penetration attenuation characteristic κ, resulting in the effective seepage water Q reaching the groundwater surface. _pref_gw As the preferred groundwater recharge, i.e. F _pref The product of κ and the remaining water volume R. _retain As a deep recharge amount, it is added to the deep soil water storage (WD).

[0225] Corrected groundwater content:

[0226] Q _GS_corr =Q _GS +Q _pref_gw ;

[0227] Unlike the above embodiments, this embodiment redefines the dynamic bypass coefficient β as:

[0228] β=0, P≤f _mat ;

[0229] β=β _max / (1+exp(-k _p (P-P _0 )))×(W / W _M ) σ , P>f _mat ;

[0230] Among them, f _mat This represents the threshold for soil matrix infiltration capacity. The actual total infiltration amount ΔF is then used. _0 Decomposed into matrix domain infiltration water volume F _matrix and preferential inflow infiltration volume F _pref :

[0231] F _matrix =(1-β)×ΔF _0 ;

[0232] F _pref =β×ΔF _0 ;

[0233] For F _matrix The soil is replenished layer by layer, in the order of upper soil, lower soil and deep soil.

[0234] Let W be the current water storage capacity of the upper soil layer, lower soil layer, and deep soil layer. _U W _L and W _D The corresponding maximum water storage capacities are W _UM W _LM and W _DM Then we have:

[0235] ΔW _U =min(F_matrix W _UM -W _U );

[0236] ΔP _1 =F _matrix -ΔW _U ;

[0237] ΔW _L =min(ΔP _1 W _LM -W _L );

[0238] ΔP _2 =ΔP _1 -ΔW _L ;

[0239] ΔW _D =min(ΔP _2 W _DM -W _D );

[0240] ΔP _3 =ΔP _2 -ΔW _D ;

[0241] The water storage updates for the upper, lower, and deep soil layers are as follows:

[0242] W _U (t+Δt)=W _U +ΔW _U ;

[0243] W _L (t+Δt)=W _L +ΔW _L ;

[0244] W _D (t+Δt)=W _D +ΔW _D ;

[0245] Preferred inflow of infiltration volume F _pref Based on the penetration attenuation characteristic κ, the preferential flow groundwater recharge Q is then classified. _pref_gw and deep replenishment amount R _retain .

[0246] In a preferred embodiment, a penetration attenuation characteristic is constructed based on the groundwater burial depth and the penetration depth of the preferential flow characteristics. The preferential flow groundwater recharge and the remaining water volume are superimposed on the deep soil water storage as the deep recharge amount.

[0247] In another alternative implementation, for karst watersheds with shallow groundwater depth or weak absorption effect of the vadose zone matrix, the following alternative calculation method can be used:

[0248] Q _pref_gw =F _pref ;

[0249] Q _pref_gw =g _p ×F _pref ;

[0250] For karst watersheds with shallow groundwater depth or weak absorption effect of the vadose zone matrix, the entire preferential inflow infiltration volume can be used as groundwater recharge, or a pre-configured fixed preferential flow groundwater recharge coefficient g can be adopted. _p , where g _p It can be determined through multi-objective calibration.

[0251] Furthermore, the groundwater recharge from the preferential flow is superimposed on the original groundwater component of the model to obtain the corrected groundwater component.

[0252] In this embodiment, the concept of water supply degree in the slope runoff process is extended to the river confluence process, defining the generalized water supply degree X of the river. _riv And calculate the dynamic transport loss Q in the river channel. _loss :

[0253] Q _loss =Q _in X _riv ≤α _0 * ;

[0254] Q _loss =f _max_ch ×[X _riv -1 / (τ+2)×(X) _riv -α _0 * ) (τ+2) / ((1-α _0 * ) (τ+1) )], α _0 * <X _riv <1;

[0255] Q _loss =f _max_ch ×(α _0 * +(1-α _0 * ) / (τ+2)), X _riv ≥1;

[0256] Where, α _0 * Let τ be the minimum relative infiltration capacity of the riverbed, and τ be the shape parameter of the infiltration capacity distribution of the riverbed.

[0257] Incorporate the dynamic transport loss of the river channel as a source-sink term into the river continuity equation:

[0258] (dW _riv ) / dt=Q _in -O-Q _loss ;

[0259] Among them, W _riv O represents the water storage capacity of the river channel, and O represents the outflow of the river channel.

[0260] The effective inbound traffic after deducting dynamic transmission losses is:

[0261] Q _e =Q _in -Q _loss ;

[0262] Next, the Muskingen method is used to calculate river confluence, outputting flood forecast results. The matrix saturation runoff component (including interflow and matrix groundwater components) and the groundwater runoff component updated by merging with the preferential groundwater recharge are combined to form the aggregated result of the matrix saturation runoff component and the preferential groundwater recharge. The obtained excess surface runoff component ΔR _s The components of interflow and subsurface runoff are input into the river confluence calculation module to obtain the flood discharge process Q(t) at the control section of the hydrological station, and the total flood volume W is calculated. _f Peak flow Q _max And peak time T _p In this embodiment, parameter α can be determined by manual calibration, automatic optimization, or a combination of both. _PT S _r c, K _s α _0 b, β _max g _p P _0 α _0 * τ, K, and ε, etc. The calibration target can comprehensively consider the relative error of the flood peak, the Nash efficiency coefficient, the total flood volume error, and the peak occurrence time error, so as to improve the comprehensive forecasting performance of the method of this invention in typical flood processes in arid and semi-arid plain areas.

[0263] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A flood forecasting method applicable to arid and semi-arid plains, characterized in that, include: Acquire hydrological and meteorological observation data and underlying surface data for the target watershed, and calculate the effective rainfall for the time period; Based on the effective rainfall and underlying surface data for a given period, excess infiltration runoff is determined to obtain the excess surface runoff and the actual total infiltration. Based on the dynamic bypass mechanism determined by the effective rainfall and hydrological and meteorological observation data for a given period, the actual total infiltration is decomposed into matrix domain infiltration and preferential inflow infiltration. The infiltration water volume in the matrix domain was treated by layer-by-layer soil filling to obtain the matrix full runoff component; Based on the underlying surface data, the penetration attenuation characteristics are determined, and the preferential inflow infiltration volume is divided into preferential groundwater recharge and deep recharge. The deep recharge volume is used for runoff calculation in subsequent periods. The river inflow is obtained by summing the excess surface runoff, the matrix saturation runoff component, and the preferential groundwater recharge. The dynamic transmission loss during the evolution process is calculated and deducted. Then, the river confluence calculation is performed, and the flood flow process forecast results of the forecast section are output. The actual total infiltration is decomposed into matrix domain infiltration and preferential inflow infiltration, including: The dynamic bypass coefficient is determined based on the current rainfall intensity obtained from the analysis of effective rainfall over a period of time, the underlying surface data and the soil moisture state characterized by the current soil water storage, and the macropore development state determined based on hydrological and meteorological observation data. The actual total infiltration volume is proportionally allocated using a dynamic bypass coefficient. The water volume allocated to the macroporous channels is taken as the priority infiltration volume, and the remaining infiltration volume is taken as the matrix domain infiltration volume. Determining the dynamic bypass coefficient includes: Extract the duration of continuous, ineffective rainfall preceding a rainfall event from hydrological and meteorological observation data; An early drought-induced fracture development factor was constructed based on the duration and a pre-defined timescale of fracture development characteristics. The dynamic bypass coefficient was calculated by using the drought fissure development factor as the modulation weight and combining it with the rainfall intensity and soil moisture status obtained from the analysis of effective rainfall over time. The preferential inflow of seepage water is divided into preferential groundwater recharge and deep recharge, including: Obtain representative groundwater depth and preferential flow penetration depth in the target area; Based on the ratio relationship between the two, a penetration attenuation characteristic with an exponential decreasing law is constructed; The preferential inflow infiltration volume is extracted by reducing the penetration attenuation characteristics. The effective penetration volume reaching the groundwater surface after reduction is taken as the preferential groundwater recharge volume. The remaining water volume absorbed during vertical transport is taken as the deep recharge volume.

2. The method according to claim 1, characterized in that, The effective rainfall for the calculation period includes: Extracting multi-point rainfall sequences and meteorological evaporation elements from hydrological and meteorological observation data; Spatial weights are used to spatially weight the multi-location rainfall sequences to obtain the average rainfall over the watershed area; The potential evapotranspiration is calculated based on meteorological evaporation factors, and the actual evapotranspiration is obtained by reducing the potential evapotranspiration based on the soil moisture stress state characterized by the underlying surface data. The effective rainfall for a given period is determined based on the difference between the average rainfall over the watershed and the actual evapotranspiration.

3. The method according to claim 1, characterized in that, When calculating the dynamic bypass coefficient, a preset effective rainfall threshold is also used to constrain the initiation of the preferential flow: When the rainfall intensity is less than or equal to the effective rainfall threshold, the dynamic bypass coefficient is set to zero, the actual total infiltration is divided into matrix domain infiltration, and the preferential infiltration is set to zero. When the rainfall intensity is greater than the effective rainfall threshold, the dynamic bypass coefficient calculated by using the previous drought fissure development factor as the modulation weight satisfies the preset maximum bypass upper limit constraint.

4. The method according to claim 1, characterized in that, The infiltration water volume in the matrix domain was treated with layer-by-layer soil filling to obtain the matrix-saturated runoff component, including: Obtain the current water storage of each soil layer, and calculate the difference between the current water storage and the maximum water storage capacity of each soil layer in the order of upper soil layer, lower soil layer and deep soil layer, as the replenishment capacity of each layer. The infiltration water volume of the matrix domain is replenished layer by layer in sequence, the current water storage of each soil layer is updated accordingly, and the data is saved to the system status for use in the next calculation period. After the top layer, bottom layer and deep layer of soil are filled in sequence, the remaining water is input into the preset free water storage module, and the matrix full-storage flow component is calculated according to the preset outflow distribution rules.

5. The method according to claim 1, characterized in that, Based on historical river flow records from hydrological and meteorological observation data, the dynamic transmission loss of river inflow during its evolution is calculated, including: The maximum potential infiltration capacity of the riverbed is extracted from the underlying surface data, and the ratio of the river inflow to the maximum potential infiltration capacity of the riverbed is calculated to obtain the generalized water supply of the river. The early moist state of the riverbed is assessed based on historical river flow records and water flow evaporation and drying characteristics, and the early moist state of the riverbed is updated synchronously according to the continuous inflow of river flow during the current calculation period. By combining the generalized water supply of the river channel with the previously moist state of the riverbed after dynamic updates, spatial infiltration reduction calculations are performed to obtain the dynamic transmission loss.

6. The method according to claim 1, characterized in that, Output the flood discharge process forecast results for the forecast section, including: The dynamic transmission loss is treated as a source-sink dissipation term in the river continuity equation and deducted from the river inflow to obtain the effective river inflow. Based on the hydrodynamic storage and discharge characteristics of the target river section, the discrete confluence method is used to calculate the flow propagation of the effective river inflow in time intervals, and the flood flow process forecast results of the forecast section are obtained.

7. The method according to claim 1, characterized in that, During the runoff determination process, the excess permeability surface runoff and the matrix saturation runoff component are dynamically switched or coexist in the same calculation period through the following methods: The order of magnitude of rainfall intensity obtained from the analysis of effective rainfall over a period of time and the infiltration capacity characterized by underlying surface data are used to trigger the infiltration runoff generation mechanism in real time, generating infiltration surface runoff. Based on the soil tension water full state determined by the current soil water storage, the full water production mechanism is triggered in real time to generate the matrix full water production component. When rainfall intensity exceeds infiltration capacity, runoff is generated primarily by the excess infiltration runoff generation mechanism; When rainfall intensity is less than or equal to infiltration capacity and soil water storage reaches its maximum storage capacity, runoff is generated primarily by the full-storage-and-run generation mechanism.

Images