A slope surface runoff concentration method based on path-grid dual-domain coupling and local adaptive constraint

Abstract

The application discloses a kind of slope surface production concentration methods based on path-grid dual-domain coupling and local adaptive constraint, specifically related to hydrological simulation and basin flood forecasting technical field, including: by runoff obtaining surface runoff, interflow and groundwater recharge, and record its concentration path information;According to the start and end position of water flow, the path section and grid passed are determined, and the propagation time of water flow in each grid is calculated;Theoretical infiltration capacity is calculated by propagation time, and actual infiltration capacity is determined by two-stage correction of path water quantity constraint and grid capacity constraint;When competition or error is out of limit, affected section is reduced to recalculate time step;According to actual infiltration capacity, update path water quantity and soil moisture state, and the groundwater recharge is merged into linear reservoir to calculate base flow, and surface runoff, interflow are jointly used as the input of river concentration module.

Description

Technical Field

[0001] This invention relates to the field of hydrological simulation and watershed flood forecasting technology, and more specifically, to a slope runoff generation and confluence method based on path-grid dual-domain coupling and local adaptive constraints. Background Technology

[0002] Watershed runoff generated by rainfall typically consists of surface runoff, interflow, and groundwater baseflow. Surface runoff mainly flows rapidly along slopes and is a crucial component in the formation of flood peaks and the control of their timing. Interflow primarily moves laterally along the slope in shallow soil or weathered layers, playing a significant role in simulating flood rise and fall processes, as well as small to medium-sized floods. Deep seepage replenishes groundwater and ultimately enters the river channel as baseflow. Because the formation conditions, propagation media, and velocities of different runoff components vary significantly, runoff generation and confluence methods capable of describing the formation and confluence processes of these different runoff components are needed.

[0003] Existing watershed runoff simulation methods can be broadly categorized into two types. One type includes total runoff time methods, isochronous methods, and unit hydrograph methods. These methods are computationally efficient, but typically use total runoff time as the primary description of runoff propagation, simplifying the consideration of state changes along the propagation path. This makes it less suitable for further describing infiltration along the path, intermediate state updates, and the mutual influence of different paths within shared spatial units. The other type employs grid-by-grid distributed models or two-dimensional hydrodynamic models. These methods can comprehensively describe slope hydrodynamic processes, but involve numerous computational states, requiring significant computational time steps and resources. When the research objective primarily focuses on slope runoff propagation, infiltration along the path, and river inflow timing, using a grid-by-grid hydrodynamic solution across the entire watershed may impose unnecessary computational burden.

[0004] In addition, local time step technology has been studied and applied to hydrodynamic simulation, but there is no known technical solution for applying it to the accuracy control of infiltration calculation during slope runoff generation.

[0005] Therefore, existing runoff generation and confluence methods still face a certain balance between computational efficiency and process description along the slope. It is necessary to propose a slope runoff generation and confluence calculation method that takes into account both computational efficiency and process description capabilities. Summary of the Invention

[0006] To achieve the above objectives, the present invention provides the following technical solution: a slope runoff generation and confluence method based on path-grid dual-domain coupling and local adaptive constraints, comprising:

[0007] S1. Obtain basic watershed data and discretize it into calculation units. Obtain surface runoff, interflow recharge and groundwater recharge through runoff calculation.

[0008] S2. Establish orderly confluence paths from surface runoff and interflow recharge to the river channel, and record the path number, path progress and remaining water volume for each water volume unit in transit.

[0009] S3. Based on the start and end positions of the water flow along the confluence path during the current calculation period, determine the actual path segments traversed by the water flow, identify the grid cells corresponding to the path segments, and calculate the actual propagation time of the water flow in each grid cell.

[0010] S4. Calculate the theoretical infiltration amount corresponding to each path based on the actual propagation time; determine the path water volume constraint coefficient based on the current available infiltration water volume of the path, and correct the theoretical infiltration amount of each path so that the corrected infiltration amount of each path does not exceed its current available infiltration water volume; determine the grid capacity constraint coefficient based on the grid's maximum infiltration capacity in the current period, and uniformly correct the infiltration amount of each path in the same grid after the path water volume constraint correction so that the actual total infiltration of each path does not exceed the grid's maximum infiltration capacity in the current period, thereby determining the actual infiltration amount of each path.

[0011] S5. When the theoretical infiltration amount corresponding to multiple paths exceeds the infiltration capacity that the grid can accept, or the local calculation error exceeds the preset threshold, only the affected path segment is recalculated with a smaller time step for the water flow propagation process, grid infiltration loss and path remaining water volume in that segment. The remaining infiltration capacity of the grid is updated segment by segment. The refinement is terminated when the preset minimum subdivision time interval or maximum refinement iteration number is reached. If the error still exceeds the preset threshold when it is terminated, the current minimum step size calculation result is retained and the local error in this period is marked.

[0012] S6. Update the in-transit water volume, path progress, and grid soil moisture status according to the actual infiltration volume, collect the groundwater recharge into the groundwater linear reservoir to calculate the baseflow, and use surface runoff, interflow, and baseflow as inputs to the river confluence module.

[0013] Preferably, in step S1, the runoff calculation adopts one of the following: over-permeability runoff mechanism, full-storage runoff mechanism, or a combination of over-permeability and full-storage mechanism, or receives the runoff results output by an external hydrological model;

[0014] The amount of water entering the soil is divided into soil water storage increment, interflow recharge, and deep infiltration, with deep infiltration serving as groundwater recharge.

[0015] Preferably, in step S2, the surface runoff path is established based on the flow direction, slope, and river location of the digital elevation model;

[0016] The interflow path is established based on the shallow aquifer interface, slope aspect, and interflow direction, or by using an equivalent interflow path with the same topology as the surface runoff path but with independently set propagation parameters.

[0017] Preferably, in step S2, the path progress is represented by one or more of the following: cumulative propagation distance, cumulative propagation time, current path segment number and position within the segment, remaining propagation distance, or remaining propagation time.

[0018] Preferably, in step S3, the actual path segment record of the water flow includes the grid cells that the water volume unit actually passes through during the current calculation period, as well as the actual propagation time, propagation length, and effective contact area in each grid cell.

[0019] Preferably, in step S4, the maximum infiltration capacity of the grid for surface runoff is determined by the remaining water storage space of the soil, the current infiltration capacity, the length of the calculation period, and the grid area.

[0020] The maximum infiltration capacity of the grid for interflow is determined by both the water storage space of the shallow soil and the permeability of the deep soil layer.

[0021] Preferably, in step S5, whether the theoretical infiltration amount corresponding to multiple paths exceeds the infiltration capacity that the grid can accept is determined based on the ratio of the total theoretical infiltration amount of multiple paths within the same grid to the maximum infiltration capacity of the grid.

[0022] Local calculation error is determined based on the difference between the calculation result of the current time step and the refined calculation result of the time step;

[0023] Local adaptation is only carried out in local areas where there is water competition or error exceeding the limit. Local areas include competing grids and related confluence paths or connected regions within their influence range.

[0024] Preferably, in step S6, the groundwater linear reservoir uses sub-basins as the basic calculation unit, sums up the deep seepage of each grid in the sub-basin as the groundwater recharge, and calculates the baseflow using a linear reservoir relationship that is proportional to the groundwater outflow and the groundwater storage.

[0025] The river confluence module uses the Muskingen method, linear reservoir capacity method, kinematic wave method, diffuse wave method, or Saint-Venant equation to calculate flood propagation.

[0026] The technical effects and advantages of this invention are as follows:

[0027] 1. This invention stores the water flow propagation state in the path domain and spatial attributes such as soil water storage capacity and infiltration capacity in the grid domain. The dynamic coupling relationship between the path and the grid is established only when the water flow actually passes through the corresponding grid, avoiding a large amount of invalid calculation caused by updating the entire watershed state grid by grid and time period by time. Thus, the computational efficiency is improved while ensuring the ability to describe the confluence process.

[0028] 2. This invention establishes a two-level coordination mechanism of path water volume constraint and shared grid capacity constraint, so that the actual infiltration volume of each path does not exceed the water volume it carries in transit, and the actual total infiltration volume received by each grid does not exceed its current maximum infiltration capacity. This effectively solves the problem of repeated infiltration and capacity overrun that may occur when multiple confluence paths share the same grid, and ensures water conservation and simulation stability during the calculation process.

[0029] 3. Based on the degree of infiltration competition among multiple paths within the grid and local calculation errors, this invention only reduces the time step for local areas where there is competition or errors exceed the limit, while maintaining the original calculation step size in other areas. This not only improves the accuracy of local infiltration calculations through time step refinement, but also avoids the significant increase in calculation costs caused by uniformly reducing the time step size across the entire watershed, thus achieving a reasonable balance between calculation accuracy and calculation efficiency. Attached Figure Description

[0030] Figure 1 This is a flowchart of the method of the present invention.

[0031] Figure 2 This is a general conceptual diagram of slope runoff generation and confluence coupling according to the present invention.

[0032] Figure 3 This is a schematic diagram illustrating the actual path interval and the spatiotemporal correspondence between the path and the grid in the current time period of this invention. Detailed Implementation

[0033] The present invention will be further described in detail below with reference to the accompanying drawings, formulas, and embodiments. In the specific embodiments described below, each formula is only used to illustrate one possible implementation; without changing the core technical solution of the present invention, formulas with equivalent physical meaning for flow generation, flow velocity, infiltration, soil-to-soil exchange, groundwater, and river channel calculations can be used instead. Those skilled in the art should understand that the following specific embodiments are only used to illustrate the present invention by way of example and should not be regarded as limiting the scope of protection of the present invention.

[0034] As attached Figure 1 As shown, this embodiment provides a slope runoff generation and confluence method based on path-grid dual-domain coupling and local adaptive constraints, including the following steps:

[0035] S1, Data Preprocessing and Flow Calculation

[0036] S1.1: Acquire the basic data for the target watershed, including digital elevation model, river network vector data, soil type distribution data, soil hydraulic parameters, land use type data, rainfall time series data, and evapotranspiration time series data. Soil hydraulic parameters include saturated hydraulic conductivity, porosity, field capacity, and wilting water content. Unify all the above data to the same coordinate system and spatial resolution. Discretize the watershed into regular raster computational units, irregular slope units, or computational units with unique spatial indexes.

[0037] Hydrological preprocessing is performed on the digital elevation model, including depression filling, river network correction, and flow direction calculation. Flow direction calculation methods can include the D8 single-flow-direction method, the D∞ multi-flow-direction method, the multi-flow-direction method, or manual correction of the flow direction. River grids are identified based on river network vector data, runoff accumulation thresholds, or a combination of both. For each non-river calculation unit, the flow direction is traced step-by-step to the first encountered river grid, which is then identified as the surface runoff target river unit for that calculation unit, thus forming the surface runoff confluence path topology for each calculation unit.

[0038] The confluence path of interflows is established based on soil structure information, distribution of shallow high-conductivity layers, and soil-bedrock interface morphology or slope aspect information. When subsurface data is insufficient, the same topological path as the surface runoff path can be used, but independent velocity parameters, water storage space parameters, and confluence parameters are set for the interflows.

[0039] S1.2: Calculate the runoff of each calculation unit based on rainfall, evapotranspiration, and soil water storage status. This invention does not limit the specific runoff model; any runoff method capable of outputting three types of water sources—surface runoff, interflow, and groundwater recharge—is applicable.

[0040] The following explanation uses the Xin'anjiang Three-Source Model as an example to illustrate the runoff calculation process. The Xin'anjiang Model, based on the theory of runoff generation at full capacity, is suitable for hydrological simulation in humid and semi-humid regions. The model's runoff calculation consists of two steps: first, calculating the total runoff; then, using a free-water reservoir, dividing the total runoff into three sources: surface runoff, interflow, and groundwater runoff.

[0041] A three-layer evapotranspiration calculation model was adopted, dividing soil tension moisture into three layers: upper, lower, and deep. The watershed evaporation capacity was calculated from measured water surface evaporation, as shown in formula (1):

[0042] Formula (1): ,

[0043] In formula (1), For grid In the Evaporation capacity of the watershed over a given period, expressed in millimeters; This is the evapotranspiration conversion factor, which is dimensionless; For grid In the The measured water surface evaporation during the specified time period is expressed in millimeters.

[0044] Evapotranspiration in three layers is calculated using formulas (2) to (5):

[0045] Formula (2): ,

[0046] Formula (3): ,

[0047] Formula (4): ,

[0048] Formula (5): ,

[0049] In the above formulas, , , The first Evapotranspiration of upper, lower, and deep layers during a given time period, in millimeters; , The first Initial soil moisture content of upper and lower layers at the beginning of the time period, in millimeters; The lower layer tension water capacity is expressed in millimeters. is the deep evapotranspiration coefficient, which is dimensionless; For the first Total evapotranspiration over a period of time, expressed in millimeters, when the upper layer moisture content is... hour, ,at this time Evaporation in the lower and deeper layers automatically becomes zero.

[0050] Net rainfall during a given period is calculated using formula (6):

[0051] Formula (6): ,

[0052] In formula (6), For the first Time Grid Net rainfall, in millimeters; For the first Time Grid Rainfall amount, in millimeters; For the first Time Grid Total evapotranspiration, expressed in millimeters. Net rainfall and precipitation are expressed in... and This indicates that the same symbol should not be used to represent two physical quantities.

[0053] The assumption of full runoff generation is adopted, and the spatial distribution of watershed water storage capacity is described by tension water storage capacity curve, as shown in formula (7):

[0054] Formula (7): ,

[0055] In formula (7), For water storage capacity less than or equal to The area, in square kilometers; The total area of ​​the basin is expressed in square kilometers. The water storage capacity at a point within the basin, expressed in millimeters; The maximum point water storage capacity within the basin, expressed in millimeters; It is a dimensionless index representing the uneven distribution of water storage capacity.

[0056] The relationship between the average water storage capacity of the basin and the maximum water storage capacity is shown in formula (8):

[0057] Formula (8): ,

[0058] In formula (8), The average tension water storage capacity of the grid watershed is expressed in millimeters.

[0059] The ordinate of the curves showing the initial average water storage capacity of the basin and the preceding impact rainfall is calculated using formula (9):

[0060] Formula (9): ,

[0061] In formula (9), The preceding rainfall is measured in millimeters and is equal to the average soil water storage in the watershed at the beginning of the period. ; The vertical axis of the water storage capacity curve corresponds to the amount of rainfall that affected the previous period, and the unit is millimeters.

[0062] Total output flow is calculated in two scenarios:

[0063] Partial area runoff, i.e. When, calculate according to formula (10):

[0064] Formula (10): ,

[0065] Whole-basin runoff, i.e. When, calculate according to formula (11):

[0066] Formula (11): ,

[0067] when hour, .

[0068] After the calculation is completed, update the soil moisture content at the end of the time period according to formula (12):

[0069] Formula (12): ,

[0070] Formula (12) indicates that the soil moisture content at the end of the period is equal to the soil moisture content at the beginning of the period plus the net rainfall and then minus the runoff, which conforms to the principle of water balance.

[0071] The total runoff is divided into three sources: surface runoff, interflow, and groundwater runoff. The Xin'an River three-source model incorporates a free-water reservoir for source classification.

[0072] The proportion of runoff area is calculated according to formula (13):

[0073] Formula (13):

[0074] In terms of runoff generation area, the water balance of the free water reservoir is as shown in formula (14):

[0075] Formula (14): ,

[0076] In formula (14), , The first Free water depth at the beginning and end of the time period (over runoff area), in millimeters; To convert the average runoff of the entire basin into the inflow depth over the runoff area; , , These represent the depths of surface runoff, interflow, and groundwater runoff over the runoff-producing area, respectively, in millimeters.

[0077] The three outflow components are calculated according to formulas (15), (16), and (17) (all representing water depths over the runoff area):

[0078] Formula (15): ,

[0079] Formula (16): ,

[0080] Formula (17): ,

[0081] Among them, surface runoff is calculated according to the overflow mechanism of free water reservoirs: when the free water storage on the runoff-producing area... Exceeding the free water storage capacity At that time, the excess overflows and forms surface runoff; The free water storage after surface runoff overflow and before subsurface runoff outflow is calculated according to formula (18):

[0082] Formula (18): ,

[0083] In formulas (16) and (17), , , respectively, are the outflow coefficients of interflow and groundwater runoff, dimensionless; Free water storage capacity, in millimeters. Free water storage capacity at the end of the period is... .

[0084] Formulas (15) to (17) above calculate the water depth over the runoff-generating area. In actual runoff calculations, the depth needs to be converted to the entire catchment area and calculated using formulas (19), (20), and (21):

[0085] Formula (19): ,

[0086] Formula (20): ,

[0087] Formula (21): ,

[0088] In the above formulas, The average surface runoff depth for the entire basin; The average depth of water flowing through the soil across the entire basin; The average groundwater runoff depth for the entire basin is expressed in millimeters. Formulas (15) to (17) are not multiplied by the runoff area ratio. The conversion is only performed uniformly in formulas (19) to (21) to avoid It is multiplied repeatedly.

[0089] The amount of water entering the soil is divided into soil water storage increment, interflow recharge, and deep infiltration, with deep infiltration serving as groundwater recharge.

[0090] In practical applications, other water source classification methods can be used to replace the above-mentioned free water reservoir model, such as the Horton super-infiltration runoff model, the Green-Ampt infiltration model, the SCS curve number method, or other distributed runoff models. It is only necessary to ensure that the output consists of three water sources: surface runoff, interflow, and groundwater runoff.

[0091] like Figure 2 As shown, the surface runoff obtained from runoff generation calculations evolves towards the river channel along the surface runoff confluence path, the interflow recharge evolves towards the river channel along the interflow confluence path, and the groundwater recharge enters the groundwater linear reservoir. The three types of water sources enter the river channel with different propagation speeds and path forms, and ultimately serve as the input to the river confluence module.

[0092] S2. Establish the confluence path and calculate the propagation speed.

[0093] S2.1: Construct ordered confluence paths for each in-transit water volume unit of surface runoff and interflow. An in-transit water volume unit refers to a water volume component that has not yet reached the river channel and is moving along the confluence path. Surface runoff paths are established based on the flow direction, slope, and river channel location using a digital elevation model; interflow paths are established based on the shallow aquifer interface, slope aspect, and interflow direction, or using an equivalent interflow path with the same topology as the surface runoff path but with independently set propagation parameters.

[0094] S2.2: The velocities of surface runoff and interflow are calculated according to formula (22):

[0095] Formula (22): ,

[0096] In formula (22), Runoff type Water volume unit Path segment The equivalent transport velocity of water flow, expressed in meters per second; Represents surface runoff, Represents the flow of the soil; For grid Corresponding runoff type The velocity coefficient; For grid Corresponding runoff type The hydraulic gradient is dimensionless. The velocity index; Runoff type The minimum constraint velocity is expressed in meters per second. When the velocity calculated from the slope and roughness is less than the minimum constraint velocity, the minimum constraint velocity is taken to ensure that the velocity in formula (23) is not zero and to avoid the denominator being zero.

[0097] Flow rate coefficient and flow velocity index The parameters for surface runoff can be determined by comprehensively considering runoff type, land use type, soil texture, and surface roughness. For example, the Manning-type slope flow parameters can be used. Generally taken Inter-soil flow is controlled by seepage through porous media. Generally taken The above parameters can also be obtained by calibration from measured rainfall-runoff data.

[0098] For surface runoff or interflow, the propagation time of the path segment is uniformly expressed by formula (23):

[0099] Formula (23): ,

[0100] In formula (23), Water volume unit Path segment Runoff type The single-segment propagation time, in seconds, is calculated once during the model initialization phase and remains constant throughout the entire process. For path segment The length is expressed in meters.

[0101] The assumption that the propagation time remains constant throughout the entire single segment applies to slope confluence scenarios where the flow rate variation is small and the flow velocity is mainly controlled by the slope and surface roughness. When the flow rate variation is large, the velocity coefficient can be dynamically adjusted based on the remaining water volume during the local adaptive time period described in S5, or the kinematic wave velocity formula can be used instead. (q is the unit width flow rate) to take into account the impact of flow rate on flow velocity.

[0102] The cumulative propagation time along the path is calculated using formula (24):

[0103] Formula (24): ,

[0104] In formula (24), Water volume unit Runoff type ,forward Total cumulative propagation time of the segment path, in seconds; This refers to the sequence number of the path segment participating in the runoff calculation along the runoff path. With the unique grid unit Correspondingly, path segment length That is, the water flow in the corresponding grid. The propagation length within; the path segment index in subsequent formulas. With grid subscript They represent the same spatial location and are only used to distinguish between path attributes and raster attributes.

[0105] S3. Determine the path segments and establish the spatiotemporal correspondence between the path and the raster.

[0106] This step describes the complete process of updating the propagation progress of water units in transit along the confluence path in the path domain and recording the correspondence between the actual paths they traverse and the grid.

[0107] S3.1: For each water volume unit in transit within the path domain, record the following status information: flow type, path number, current path progress, remaining water volume, target inflow unit, and index of the current path segment. Path progress is represented by one or more of the following: cumulative propagation distance, cumulative propagation time, current path segment number and position within the segment, remaining propagation distance, or remaining propagation time.

[0108] S3.2: Calculate the time interval step based on the current model, and calculate the prediction end progress according to formula (25):

[0109] Formula (25): ,

[0110] In formula (25), For the first After time period projection, runoff type Path unit The cumulative propagation progress of the predicted path, in seconds; For the first Initial time period, runoff type Path unit The cumulative propagation progress along the path, in seconds; Runoff type Path unit The total cumulative propagation time to the end of the complete path, i.e., the upper limit of the total duration of the path convergence, is expressed in seconds.

[0111] like If the water volume unit reaches the end of the path within this time period, it will flow into the corresponding river channel unit after exchange and settlement along the route; if If the updated progress and remaining water volume are saved, the process will continue from the current position in the next time period. This distinguishes between water volume that has already entered the river and water volume that is still en route.

[0112] S3.3: Determine the correspondence between the actual path traversed by the water volume unit and the grid during the current time period. The actual path segment record of the water flow includes the grid cells that the water volume unit actually traversed during the current calculation time period, as well as the actual propagation time, propagation length, and effective contact area in each grid cell.

[0113] like Figure 3 As shown, based on the propagation time range of the water flow along the confluence path within the current calculation period, that is... And the propagation time intervals corresponding to each path segment. The propagation time is the time range within which the water flow completes its journey along the path segment. It is a pre-calculated static attribute inherent to the path. The actual propagation time of the water flow in each grid cell is calculated, and its value is the length of the intersection of the two time intervals. When the two intervals do not overlap, the propagation time is zero.

[0114] Path propagation information only records the actual paths and grid combinations that the water flow passes through during the current calculation period, and does not record combinations that have not propagated, thereby reducing the amount of computation and improving computational efficiency.

[0115] S4. Calculation of Theoretical Infiltration and Dual Constraints

[0116] S4.1: Calculate the theoretical infiltration rate for each path based on the actual propagation time of water flow within the grid. For surface runoff, This indicates the current infiltration capacity of the grid; for interflow, This indicates the current groundwater recharge capacity of the grid, i.e., the deep seepage capacity. The theoretical infiltration rate is calculated using formula (26):

[0117] Formula (26): ,

[0118] In formula (26), For the first Time period, runoff type ,path Passing through the grid Unconstrained theoretical infiltration rate, in cubic meters; For the first Time Grid Corresponding runoff type The infiltration capacity, measured in millimeters per hour; For the first Time-of-day runoff type ,path In the grid The actual transmission time within the time frame, in hours; For grid Area, in square meters; Runoff type ,path In the grid The effective contact correction factor is dimensionless; it is divided by 1000 as a unit conversion factor. When the infiltration capacity is expressed in millimeters per hour, the propagation time in hours, and the grid area in square meters, the calculation result is converted to cubic meters.

[0119] Effective contact correction factor The actual contact degree of water flow within the grid is reflected and calculated according to formula (27):

[0120] Formula (27): ,

[0121] In formula (27), The effective contact area of ​​the water flow within the grid is determined by the width of the water flow and its propagation length within the grid. Sure. The range of values ​​is When measured data on water flow width are lacking, empirical values ​​can be obtained based on land use type and degree of gully formation, or calibrated using measured infiltration data.

[0122] S4.2: Calculate the path water volume constraint coefficient based on the theoretical infiltration amount calculated for each path in the grid cells along the path. When the sum of the theoretical infiltration amounts of a path in multiple grid cells exceeds the current in-transit water volume it carries, the theoretical infiltration amounts of each grid cell are uniformly adjusted according to the path water volume constraint coefficient to ensure that the total infiltration amount of the path does not exceed its current available infiltration water volume. The path water volume constraint coefficient is calculated according to formula (28):

[0123] Formula (28): ,

[0124] In formula (28), For the first Time-of-day runoff type ,path The water quantity constraint coefficient is dimensionless. For the first Time period path Total water volume carried in transit, in cubic meters; For path The total theoretical infiltration volume of the entire grid, in cubic meters; This is a correction term for the denominator of the minimum value to prevent the denominator from being zero and dimensionless.

[0125] After path flow constraint correction, runoff type ,path In the grid The infiltration rate is calculated according to formula (29):

[0126] Formula (29): ,

[0127] S4.3: Summarize the theoretical infiltration amounts of each path of the same runoff type in the same grid, and calculate the theoretical infiltration capacity of the grid for the current period according to formula (30):

[0128] Formula (30): ,

[0129] In formula (30), For the first Time Grid After superimposing all internal paths, the runoff type Total theoretical infiltration volume, in cubic meters; For grid Internal flow type The set of paths.

[0130] For surface runoff, the maximum capacity that the grid can accept in the current time period is limited by the remaining water storage space of the soil and the maximum infiltration capacity of the time period, and is calculated according to formula (31):

[0131] Formula (31):

[0132] In formula (31), For the first Time Grid The maximum amount of infiltrated water that surface runoff can accept, in cubic meters; The maximum water storage capacity of the grid soil is expressed in millimeters, to distinguish it from the watershed average tension water storage capacity in formula (8). ; The output of formula (12) is the grid soil water storage after S1 runoff update and before friction exchange, in millimeters; Ensure that the acceptable capacity of the soil is not negative when it is full. For the first Time Grid Maximum soil infiltration capacity, measured in millimeters per hour; The time step for model calculation is in hours.

[0133] For interflows, the maximum acceptable capacity of shallow interflow units is expressed by formula (32):

[0134] Formula (32): ,

[0135] In formula (32), For the first Time Grid The maximum amount of water that can be received by the soil flow, in cubic meters; This represents the maximum water storage capacity in the shallow soil layer of the grid, in cubic meters. For the first Initial water storage in shallow soil of the grid during a given time period, in cubic meters; For the first Vertical saturated hydraulic conductivity (i.e., deep permeability) of shallow soil in a time-segmented grid, in millimeters per hour. This refers to the grid area, in square meters. The time period is expressed in hours; divide by 1000 to convert to cubic meters. Interflow seepage into deeper layers is a vertical movement controlled by vertical hydraulic conductivity; therefore, deep seepage capacity, rather than lateral hydraulic capacity, is used as the limiting factor.

[0136] S4.4: Calculate the grid capacity constraint coefficient based on the grid's maximum infiltration capacity, using formula (33):

[0137] Formula (33): ,

[0138] In formula (33), For the first Time Grid Runoff type The capacity constraint coefficient is dimensionless. For grid Runoff type Maximum acceptable infiltration rate, in cubic meters; For grid Runoff type Total theoretical infiltration volume, in cubic meters; This is a correction term for the denominator of the minimum value to avoid division by zero.

[0139] This invention establishes a two-tiered coordination mechanism: path water volume constraint and shared grid capacity constraint. First, a path water volume constraint coefficient is determined based on the current available infiltration volume of each path, performing a first-level correction to the theoretical infiltration volume of each path. Then, a grid capacity constraint coefficient is determined based on the maximum infiltration capacity of the grid during the current time period, performing a second-level unified correction to the infiltration volume of each path, ultimately determining the actual infiltration volume of each path. This coefficient is used to uniformly adjust the infiltration volume of each path within the grid, resulting in dual constraints: path water volume constraint and grid capacity constraint. Subsequently, the runoff type... ,path In the grid The actual infiltration amount is calculated according to formula (34):

[0140] Formula (34): ,

[0141] The process ensures two conditions simultaneously: first, the total actual infiltration volume of each path does not exceed the current in-transit water volume carried by that path; second, the total actual infiltration volume received by each grid does not exceed the maximum infiltration capacity of that grid in the current time period.

[0142] This invention employs a sequence of path constraint followed by grid constraint: path constraint is performed independently within each path to ensure water conservation within a single path; grid constraint applies the path-constrained infiltration volume to each path within the same grid at the same ratio. Scaling does not change the relative distribution ratio between paths. Under this order, the water allocation result is unique and satisfies global water conservation.

[0143] S5, Local Time Step Adaptive Calculation

[0144] S5.1: Calculate the grid capacity competition index according to formula (35):

[0145] Formula (35): ,

[0146] In formula (35), For the first Time period, runoff type Grid The capacity competition metric is dimensionless. For grid The total theoretical infiltration volume after superimposing all internal pathways, in cubic meters; For grid Corresponding runoff type The maximum acceptable infiltration volume, in cubic meters; This is a minimum value correction term to avoid the denominator being zero.

[0147] Whether the theoretical infiltration volume corresponding to multiple paths exceeds the grid's acceptable infiltration capacity is determined by the ratio of the total theoretical infiltration volume of multiple paths within the same grid to the grid's maximum infiltration capacity. When competition indicators... When this occurs, it indicates that the maximum infiltration capacity of the grid in the current time period can meet the infiltration demand of all paths; when This indicates that the theoretical infiltration rate of multiple paths exceeds the grid's capacity limit, based on a preset competition threshold. Determine whether to initiate local time-step adaptive computation.

[0148] S5.2: To improve the accuracy of local infiltration calculations, the original time step is used. With halved time step The actual infiltration rate of the grid is calculated separately, and the relative error index is calculated based on the two sets of results according to formula (36). The local calculation error is determined based on the difference between the calculation result of the current time step and the refined calculation result of the time step. When the error exceeds the preset threshold... Local time step adaptive calculation is initiated at that time.

[0149] Formula (36): ,

[0150] In formula (36), For the first Time period, runoff type Grid The local infiltration calculation error index is dimensionless; This represents the total actual infiltration volume of the grid calculated at the original time step, in cubic meters. The sum of the actual infiltration amounts in the two sub-steps under a halved time step is expressed in cubic meters.

[0151] S5.3: The execution rules for local time step adaptive computation are as follows:

[0152] First, the calculation is only performed on local areas where there is water competition or error exceeding the limit. The local areas include the competing grid and the relevant confluence paths or connected areas within its influence range. The remaining areas maintain the original calculation step size and do not require full recalculation.

[0153] Second, a fixed bisection ratio is used for refinement. Let the number of refinement layers be l, and the substep size be determined according to formula (37a):

[0154] Formula (37a): ,

[0155] For each substep The following recursive calculations are performed sequentially:

[0156] The substep propagation progress update is determined according to formula (37b):

[0157] Formula (37b): ,

[0158] Determine the substep using the time interval intersection method in S3.3. The grid through which the internal water flows and the actual propagation time And calculate the theoretical infiltration amount of the sub-step according to formula (26). ;

[0159] The sub-step path water volume constraint coefficient is determined according to formula (37c):

[0160] Formula (37c): ,

[0161] The remaining capacity of the sub-step grid is deducted step by step according to formula (37d):

[0162] Formula (37d): ,

[0163] in That is, the initial capacity of the first sub-step is the maximum capacity that the grid can accept in the original time period;

[0164] The substep grid capacity constraint coefficient is determined according to formula (37e):

[0165] Formula (37e): ,

[0166] The actual infiltration rate of the sub-step is determined according to formula (37f):

[0167] Formula (37f): ,

[0168] The remaining water volume in the sub-step is updated according to the formula (37g):

[0169] Formula (37g): ,

[0170] When multiple paths compete within the same grid, each path is determined according to the grid capacity constraint coefficient. The water is allocated proportionally, rather than on a first-come, first-served basis, to ensure the conservation of water volume within each step.

[0171] Third, the refined results are summarized back to the original time period according to formula (38):

[0172] Formula (38): ,

[0173] Fourth, set the minimum subdivision time interval. Maximum number of refinement iterations As an efficiency constraint, when the time step is refined... Or refine the number of iterations Refinement is terminated at a certain time. If there is an error at the termination time... Still exceeds the preset threshold If the current minimum step size calculation result is retained, the local calculation error for this time period is marked.

[0174] S6. Calculation of in-transit water volume and soil water storage and renewal, groundwater baseflow and river runoff

[0175] S6.1: Based on the actual total infiltration within the corresponding grid of each path segment, update the remaining in-transit water volume of the runoff unit, and calculate according to formula (39):

[0176] Formula (39): ,

[0177] In formula (39), For the first At the end of the period, runoff type Path unit The remaining water volume in transit, in cubic meters; For the first The total water volume carried by the initial path unit during the time period, in cubic meters; This represents the total actual infiltration volume generated by this path unit across all grid path segments, expressed in cubic meters.

[0178] when When the water volume unit in transit is depleted and its propagation ceases, it will no longer participate in the calculation of subsequent time periods.

[0179] Surface runoff infiltrates into the soil aquifer along its course and is then redistributed within the soil. A grid system is installed. The actual total infiltration along the surface path in various parts of the country is As shown in formula (40a):

[0180] Formula (40a): ,

[0181] The water volume first replenishes the soil tension water storage, and the remaining part enters the shallow free water layer and is divided into interflow and deep infiltration according to the proportion, as shown in formulas (40b) to (40d):

[0182] Formula (40b): ,

[0183] Formula (40c): ,

[0184] Formula (40d): ,

[0185] In formulas (40b) to (40d), The increment of tensile water storage is expressed in millimeters. The grid soil water storage capacity after runoff renewal and before friction exchange, calculated using formula (12), is expressed in millimeters. The volume of water entering the shallow free aquifer is expressed in cubic meters. , These are the recharge from interflow formed by the conversion of infiltrated water along the course and the deep infiltration, respectively, in cubic meters; , These are the interflow distribution coefficient and deep seepage distribution coefficient for water exchange along the flow path, respectively. They are dimensionless, determined by soil hydraulic parameters, and satisfy the following conditions: Unconverted portion It remains in the shallow soil.

[0186] The soil water storage status of the grid is updated according to formula (41):

[0187] Formula (41): ,

[0188] In formula (41), The soil moisture content after updating the S1 runoff output by formula (12) (after deducting evapotranspiration during this period) and local production flow ), as the initial value for the along-the-path exchange update; For the first Initial grid soil water storage for the time period, in millimeters. Updating data along the route only includes upstream infiltration recharge and outflow from interbedded flow and deep seepage; evapotranspiration is not deducted again because evapotranspiration... It has been passed in formula (12) Deduct. Formula (12) updates the soil moisture content change caused by rainfall runoff in this grid, and formula (41) updates the soil moisture content change caused by infiltration recharge along the upstream path; both act on the same soil moisture state variable, and in the same time period, S1 runoff update is performed first, and then formula (41) is performed for exchange update along the path, so as to jointly ensure water balance. This is used as the initial soil water storage for the next period S1. The amount of shallow free water that is not drained in time is conservatively included in the soil water storage, which will gradually decrease in subsequent periods through evapotranspiration and exchange along the way.

[0189] This formula enables continuous water volume conversion among four types of water bodies: surface runoff, soil water storage, interflow, and groundwater.

[0190] Interflow in the soil can directly enter the river channel at the foot of the slope, in the valley, or at the exposure point of the shallow aquifer, or it can be transformed into a new surface water volume unit. The exposure determination is made according to formula (42):

[0191] Formula (42): ,

[0192] In formula (42), This refers to the thickness of the grid soil layer, in meters. The exposed thickness threshold is expressed in meters. The saturated water content is dimensionless. The unit is cubic meters. The unit is square meters. The dimension of is cubic meters, and Consistent; Pre-defined grid sets for outcropping areas such as slope toes and gullies. If any of the above conditions are met, it is determined that interflow is exposed in that grid.

[0193] The values ​​for the exposed water volume and the new surface water volume in transit are assigned according to formulas (43) and (44):

[0194] Formula (43): ,

[0195] The amount of water exposed is deducted from the original flow path, and the remaining water volume in the original path is updated according to formula (43a):

[0196] Formula (43a): ,

[0197] when At that time, the water volume unit in the soil flow is marked as exhausted and the propagation stops.

[0198] Formula (44): ,

[0199] In formulas (43), (43a), (44), The exposure ratio (determined by soil layer thickness and saturation level; the slope toe grid can be set to 1.0). The initial water volume of the transformed surface runoff; the water flow type of the new water volume unit is set to surface runoff ( ), the path number is set to the exposed grid Corresponding surface runoff path The initial path progress is set as the cumulative propagation time of the exposed raster in the path represented at that location. The target river entry unit has the same path as the local area.

[0200] S6.2: Summarize the deep seepage of each grid within each sub-basin to obtain the groundwater recharge. The groundwater storage of the sub-basin is updated according to formula (45):

[0201] Formula (45): ,

[0202] In formula (45), , The first , Groundwater storage in the sub-basin during this period, in cubic meters; This represents the total groundwater recharge of the sub-basin during this period, expressed in cubic meters. This refers to the outflow rate of groundwater, expressed in cubic meters per second. The duration is expressed in seconds.

[0203] The outflow of groundwater is calculated using a linear reservoir relationship, according to formula (46):

[0204] Formula (46): ,

[0205] In formula (46), This is the groundwater outflow coefficient, in units of... ; For the first The total groundwater storage in the sub-basin during this period is expressed in cubic meters.

[0206] This method ultimately identifies three types of river inflow components: residual water that has completed the surface runoff path propagation forms surface lateral inflow; water that has completed the interflow runoff path propagation or has emerged and transformed at the slope toe forms interflow lateral inflow; and outflow from the groundwater linear reservoir forms the river baseflow. These three types of inflow are aggregated according to the target river grid or river nodes and used as input conditions for river runoff calculation.

[0207] S6.3: Perform river-course evolution calculations for all inflow volumes. This implementation method uses the Muskingan method discrete calculation scheme to perform river confluence calculations, as shown in formula (47):

[0208] Formula (47): ,

[0209] The Muskingan model coefficients are calculated according to formulas (48), (49), and (50), respectively:

[0210] Formula (48): ,

[0211] Formula (49): ,

[0212] Formula (50): ,

[0213] The model coefficients satisfy the normalization constraint: .

[0214] In the above formulas, , The river section is the first , Outflow rate during a given time period, expressed in cubic meters per second; , The river section is the first , Inflow rate during a given time period, expressed in cubic meters per second; This is the river channel storage constant, expressed in seconds. The Muskingan flow weighting coefficient is dimensionless and is determined by a combination of river propagation time, calculation time step, and river storage characteristics.

[0215] It should be noted that the river confluence module of this invention is not limited to the Muskingen method. It can also use the linear reservoir method, the moving wave method, the diffuse wave method, or the complete Saint-Venant equation to carry out river flood propagation calculations. The appropriate calculation method can be flexibly selected according to the actual watershed scale, simulation accuracy requirements, and calculation efficiency requirements.

[0216] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A slope runoff generation and confluence method based on path-grid dual-domain coupling and local adaptive constraints, characterized in that: include: S1. Obtain basic watershed data and discretize it into calculation units. Obtain surface runoff, interflow recharge and groundwater recharge through runoff calculation. S2. Establish orderly confluence paths from surface runoff and interflow recharge to the river channel, and record the path number, path progress and remaining water volume for each water volume unit in transit. S3. Based on the start and end positions of the water flow along the confluence path during the current calculation period, determine the actual path segments traversed by the water flow, identify the grid cells corresponding to the path segments, and calculate the actual propagation time of the water flow in each grid cell. S4. Calculate the theoretical infiltration amount corresponding to each path based on the actual propagation time; determine the path water volume constraint coefficient based on the current available infiltration amount of the path, and correct the theoretical infiltration amount of each path so that the corrected infiltration amount of each path does not exceed its current available infiltration amount. The grid capacity constraint coefficient is determined based on the grid's maximum infiltration capacity during the current period. The infiltration amount of each path within the same grid, after path water volume constraint correction, is uniformly corrected so that the actual total infiltration amount of each path does not exceed the grid's maximum infiltration capacity during the current period, thereby determining the actual infiltration amount of each path. S5. When the theoretical infiltration amount corresponding to multiple paths exceeds the infiltration capacity that the grid can accept, or the local calculation error exceeds the preset threshold, only the affected path segment is recalculated with a smaller time step for the water flow propagation process, grid infiltration loss and path remaining water volume in that segment. The remaining infiltration capacity of the grid is updated segment by segment. The refinement is terminated when the preset minimum subdivision time interval or maximum refinement iteration number is reached. If the error still exceeds the preset threshold when it is terminated, the current minimum step size calculation result is retained and the local error in this period is marked. S6. Update the in-transit water volume, path progress, and grid soil moisture status according to the actual infiltration volume, collect the groundwater recharge into the groundwater linear reservoir to calculate the baseflow, and use surface runoff, interflow, and baseflow as inputs to the river confluence module.

2. The method according to claim 1, characterized in that: In step S1, the runoff calculation adopts one of the following mechanisms: infiltration runoff mechanism, storage runoff mechanism, and combination of infiltration and storage runoff mechanism, or receives the runoff results output by an external hydrological model; The amount of water entering the soil is divided into soil water storage increment, interflow recharge, and deep infiltration, with deep infiltration serving as groundwater recharge.

3. The method according to claim 1, characterized in that: In step S2, the surface runoff path is established based on the flow direction, slope, and river location of the digital elevation model; The interflow path is established based on the shallow aquifer interface, slope aspect, and interflow direction, or by using an equivalent interflow path with the same topology as the surface runoff path but with independently set propagation parameters.

4. The method according to claim 1, characterized in that: In step S2, the path progress is represented by one or more of the following: cumulative propagation distance, cumulative propagation time, current path segment number and position within the segment, remaining propagation distance, or remaining propagation time.

5. The method according to claim 1, characterized in that: In step S3, the actual path segment record of the water flow includes the grid cells that the water volume unit actually passes through during the current calculation period, as well as the actual propagation time, propagation length and effective contact area in each grid cell.

6. The method according to claim 1, characterized in that: In step S4, the maximum infiltration capacity of the grid for surface runoff is determined by the remaining water storage space of the soil, the current infiltration capacity, the length of the calculation period, and the grid area. The maximum infiltration capacity of the grid for interflow is determined by both the water storage space of the shallow soil and the permeability of the deep soil layer.

7. The method according to claim 1, characterized in that: In step S5, whether the theoretical infiltration amount corresponding to multiple paths exceeds the infiltration capacity that the grid can accept is determined based on the ratio of the total theoretical infiltration amount of multiple paths within the same grid to the maximum infiltration capacity of the grid. Local calculation error is determined based on the difference between the calculation result of the current time step and the refined calculation result of the time step; Local adaptation is only carried out in local areas where there is water competition or error exceeding the limit. Local areas include competing grids and related confluence paths or connected regions within their influence range.

8. The method according to claim 1, characterized in that: In step S6, the groundwater linear reservoir uses sub-basins as the basic calculation unit, sums up the deep seepage of each grid in the sub-basin as the groundwater recharge, and calculates the base flow using a linear reservoir relationship that is proportional to the groundwater outflow and the groundwater storage. The river confluence module uses the Muskingen method, linear reservoir capacity method, kinematic wave method, diffuse wave method, or Saint-Venant equation to calculate flood propagation.

Images