A method for simulating a compound flash flood disaster in a watershed based on a slope-gully topological structure

Abstract

The application discloses a kind of based on slope-gully topological structure's basin composite mountain torrent disaster distributed simulation method, comprising the following steps: step 1, slope-gully unit extraction and topological structure construction;Step 2, hydrological model calculation;Step 3, calculate slope-gully coupling transmission system quantitative transmission efficiency;Step 4, calculate soil saturation degree-slope risk index to assess landslide risk;Step 5, water and sediment evolution model calculation;Step 6, using multi-physical field adaptive time step strategy optimization calculation efficiency;Step 7, based on topological level parallel load balancing algorithm realizes distributed parallel computing.The present application realizes the close coupling of hydrological calculation, landslide starting judgment and water and sediment evolution calculation, can simulate the chain process of composite disaster, realizes the quantitative description of slope-gully transmission efficiency of slope runoff, improves the precision of slope-gully material exchange simulation, improves the accuracy and timeliness of mountain torrent disaster early warning and forecasting.

Description

Technical Field

[0001] This invention belongs to the field of hydrological and water resources and geological disaster simulation technology, and in particular relates to a distributed simulation method for watershed composite flash flood disasters based on slope-gully topology. Background Technology

[0002] Flash floods are one of the most serious natural disasters in mountainous areas of my country, characterized by their suddenness, destructive power, and short warning time. Under extreme rainstorm conditions, flash floods often exhibit a complex chain reaction: rainstorms trigger slope runoff and soil saturation, which in turn triggers landslides. After entering the gully, the landslide body mixes with the floodwater to form high-sediment-laden water flows or debris flows, causing even more severe disaster consequences.

[0003] Existing methods for simulating flash flood disasters mainly have the following problems:

[0004] 1) Traditional hydrological models mainly focus on rainfall-runoff processes and lack coupling mechanisms with geological disaster processes such as landslides and sediment transport;

[0005] 2) Landslide models and hydrological models usually run independently, making it impossible to achieve dynamic coupling where soil moisture content drives landslide initiation in real time;

[0006] 3) The lack of a unified topology to organize various simulation objects such as slope units, gully units, small watersheds, river sections, and nodes restricts the realization of large-scale parallel computing;

[0007] 4) Existing models lack a quantitative description of the material transport efficiency between slopes and gullies, making it difficult to accurately characterize the coupling strength under different terrain conditions;

[0008] 5) Landslide risk assessments often use a single indicator and lack a dynamic risk index that comprehensively considers multiple hydrological and mechanical factors;

[0009] 6) The load balancing strategy for parallel computing does not fully utilize the characteristics of the watershed topology, resulting in limited computational efficiency;

[0010] 7) The time step selection often uses a fixed value or a single CFL condition, failing to perform adaptive optimization for the characteristics of multi-physics coupling.

[0011] Therefore, there is an urgent need to develop a distributed simulation method for composite flash flood disasters that can couple hydrological processes, landslide processes, and water and sediment transport processes, in order to improve the accuracy and timeliness of flash flood disaster early warning and forecasting. Summary of the Invention

[0012] The purpose of this invention is to provide a distributed simulation method for watershed composite flash flood disasters based on slope-gully topology, so as to solve the above-mentioned technical problems.

[0013] To achieve the above objectives, the present invention provides the following technical solution:

[0014] This invention discloses a distributed simulation method for watershed composite flash flood disasters based on slope-gully topology, the method comprising the following steps:

[0015] Step 1: Slope-Valley Unit Extraction and Topology Construction: Based on the digital elevation model data of the study area, the D8 flow direction algorithm is used to calculate the flow direction grid and the runoff accumulation grid; a catchment area threshold is set. The process involves extracting the channel network and marking grids with cumulative runoff exceeding a threshold as channel units. Using the channel network as boundaries, the watershed is divided into several slope units, with each slope unit flowing into an adjacent channel unit, thus constructing a slope-channel binary tree topology. Nodes represent the confluence or branching points of channels. The relationship between slope units and channel units is established. Based on the small watershed-river segment-node topological relationship of the Chinese mountain flood hydrological model, an integrated topological structure is formed, encompassing five categories of objects: slope, channel, node, small watershed, and river segment.

[0016] Step 2, Hydrological Model Calculation: Perform hydrological model calculations, including rainfall interpolation, evapotranspiration calculation, soil moisture calculation, runoff generation calculation, slope runoff, and river network runoff evolution; output the results of the hydrological model calculations, including: Slope unit Soil saturation at time Slope unit outflow and small watershed outflow ;

[0017] Step 3: Calculate the slope-gully coupling transport coefficient to quantify transport efficiency: Based on the soil saturation obtained in Step 2, calculate the slope-gully coupling transport coefficient to quantify the efficiency of runoff transmission from the slope to the gully. This coefficient comprehensively considers topographic, soil, and vegetation factors, and the formula is as follows:

[0018] (5)

[0019] In the formula, For the first Slope unit The slope-channel coupling transmission coefficient at time t; For the first Topographic transfer factor for each slope unit; For the first Slope unit Soil transport factors at any given time; For the first Vegetation transport factor of each slope unit; The reference transmission coefficient;

[0020] The lateral inflow rate from the slope to the gully is then calculated as follows:

[0021] (9)

[0022] In the formula, The lateral inflow rate per unit channel length is expressed in m² / s. For the first The outflow rate of each slope unit is m³ / s; For corresponding channel units The length, in meters;

[0023] Step 4: Calculate the Soil Saturation-Landslide Risk Index to assess landslide risk: Define the Soil Saturation-Landslide Risk Index (LRI) to comprehensively assess the landslide risk of each slope unit.

[0024] (10)

[0025] In the formula, For the first Slope unit Soil saturation at any given time - landslide risk index; , , , These are the safety factor weight, saturation weight, terrain weight, and rainfall weight, respectively, satisfying... ; For the safety factor risk component; This represents the risk component of soil saturation. For terrain risk components; Accumulated risk component for rainfall;

[0026] Then, based on the soil saturation-landslide risk index (LRI), a two-level judgment of landslide risk was performed, divided into two levels: landslide risk warning and landslide initiation. Finally, the unit width flow rate of the landslide mass entering the gully was calculated. , ;

[0027] Step 5: Calculation of water and sediment evolution model: based on the outlet flow of the small watershed. Lateral inflow of coupled transmission and landslide inflow , A set of water-sediment coupling control equations was established, including the continuity equation, momentum equation, sediment transport equation, and riverbed deformation equation. The set of water-sediment coupling control equations was solved, and finally, the water depth of the channel unit was output. Water level Flow rate Unit width flow Riverbed changes and sediment concentration ;

[0028] Step 6: Optimize computational efficiency using a multiphysics adaptive time step strategy: Considering the coupling characteristics of hydrology, landslide, and sediment multiphysics, a comprehensive adaptive time step selection strategy is proposed to coordinate the time progression of each physical process; the multiphysics adaptive time step is:

[0029] (32)

[0030] In the formula, The adaptive time step for multiphysics is s; Let be the time step of the CFL condition constraint, in seconds; The time step is defined as s, which is a constraint on mass conservation. The time step for landslide risk constraints is s; Let be the time step for constraining riverbed deformation, in seconds; For safety factor;

[0031] By dynamically adjusting the time step based on local truncation error estimation, computational efficiency can be improved while maintaining accuracy.

[0032] (37)

[0033] In the formula, Tolerance; The current step is used to estimate the error; The order is in numerical format;

[0034] The time step nesting relationship between hydrological calculations and sediment evolution calculations is as follows:

[0035] (38)

[0036] In the formula, The nesting multiple of water and sand steps; The hydrological time step is the time step set during the operation of the hydrological model, expressed in seconds; this formula represents the time step at each hydrological time step. Internally, the calculation of water and sediment evolution is executed. Calculation of adaptive substep size;

[0037] Step 7: Implement distributed parallel computing based on topology-level parallel load balancing algorithm: Decompose the watershed space based on topology-level parallel load balancing algorithm, and implement distributed parallel computing using multi-process parallelism and GPU acceleration.

[0038] Furthermore, the nodes in step 1 use a binary tree encoding rule. Let the root node be encoded as 0. Then, for the node encoded as i:

[0039] (1)

[0040] (2)

[0041] (3)

[0042] The relationship between slope units and channel units is established as follows: each channel unit is associated with its left bank slope and right bank slope. The slope position is determined according to the left-right rule in the direction of water flow, that is, facing downstream, the left side is the left bank slope and the right side is the right bank slope.

[0043] Furthermore, in step 2, the first Slope unit Soil saturation at time Determined by the ratio of soil moisture content to watershed storage capacity:

[0044] (4)

[0045] In the formula, For the first Slope unit Total soil moisture content at any given time (mm); The watershed storage capacity (mm) of this slope unit.

[0046] Furthermore, the terrain transfer factor in step 3 Reflecting the impact of slope geometry on transmission efficiency:

[0047] (6)

[0048] In the formula, For the first The average slope of each slope unit; For reference slope; For the first Average flow length of each slope unit, in meters; The reference process length is in meters (m). The slope influence index;

[0049] Soil transport factors The dynamic impact of soil saturation state on transport efficiency:

[0050] (7)

[0051] In the formula, For reference saturation; For the first Soil saturated hydraulic conductivity of each slope unit, mm / h; The maximum saturated hydraulic conductivity is given in mm / h. , These are the saturation influence index and the hydraulic conductivity influence index, respectively.

[0052] Vegetation transport factors The effect of vegetation cover on transmission efficiency:

[0053] (8)

[0054] In the formula, For the first Vegetation coverage of each slope unit; For the first Leaf area index of each slope unit; Leaf area index for reference.

[0055] Furthermore, the formula for calculating the risk component of the safety factor in step 4 is as follows:

[0056] (11)

[0057] In the formula, For the first Slope unit The safety factor at any given time is greater than 1, indicating stability, and less than or equal to 1, indicating instability. This is the critical safety factor; The minimum safety factor threshold; It is a non-linear exponent;

[0058] Safety factor Calculations based on the Mohr-Coulomb criterion for infinite slopes:

[0059] (12)

[0060] In the formula, For the first Slope unit Effective cohesion at any given time, kPa; Let m be the depth of the sliding surface. Slope; The unit weight of saturated soil is kN / m³. The specific weight of water is kN / m³. For the first Slope unit The effective internal friction angle at any given time; where the soil strength parameter decreases with soil saturation:

[0061] (13)

[0062] (14)

[0063] In the formula, Effective cohesive strength in dry state, kPa; This is the effective internal friction angle in the dry state; , These are the cohesion attenuation coefficient and the internal friction angle attenuation coefficient, respectively.

[0064] The formula for calculating the soil saturation risk component is:

[0065] (15)

[0066] In the formula, Critical saturation; It is a non-linear exponent; h is the sensitivity coefficient for the rate of change of saturation. For the first The rate of change of saturation of each slope unit, 1 / h;

[0067] The formula for calculating the terrain risk component is:

[0068] (16)

[0069] In the formula, For the first The tangent of the slope of a slope unit; , These are the tangents of the minimum and maximum slopes at which the landslide occurred; For the first The relative elevation of each slope unit, in meters; The maximum relative elevation is in meters (m). This is a terrain nonlinearity index;

[0070] The formula for calculating the cumulative risk component of rainfall is:

[0071] (17)

[0072] In the formula, The critical cumulative rainfall is expressed in mm. For the first Slope unit The effective cumulative rainfall before the specified time, in mm, is calculated using exponential decay:

[0073] (18)

[0074] In the formula, For the first Each slope unit in Rainfall at any given time, in mm; The duration of the initial impact is in hours (h). h is the decay time constant.

[0075] Furthermore, the landslide risk warning in step 4 is as follows: a landslide risk warning is issued when the risk index of a slope unit reaches the critical risk index, and the warning conditions are as follows:

[0076] (19)

[0077] In the formula, This is a critical risk index;

[0078] Landslide initiation is defined as follows: Based on the met early warning conditions, when the safety factor of a slope element further decreases to an instability value, the landslide initiation is determined to have occurred in that slope element. The landslide initiation conditions are as follows:

[0079] (20)

[0080] Calculate the landslide volume, let Inversely calculate the instability thickness :

[0081] (twenty one)

[0082] The landslide volume is:

[0083] (twenty two)

[0084] In the formula, For the first The instability thickness of each slope element, in meters; For the first The area of ​​each slope unit, in m²; This is the shape correction factor; For the first The landslide volume of each slope unit, in m³;

[0085] Calculate the unit width flow rate of the landslide mass entering the gully. , The specific process is as follows: Assume the landslide body moves at a speed of... Angle along the direction of entry into the ditch During the import of history Inner width When flowing into a channel, the unit width flow rate is calculated using the following formula:

[0086] (twenty three)

[0087] (twenty four)

[0088] In the formula, , , respectively, are the unit width flow rates of the landslide inflow in the x and y directions, in m² / s; Let m be the volume of the landslide. The velocity of the sliding body as it enters the ditch is expressed in m / s. The angle of entry into the ditch; The width of the landslide entering the gully, in meters; The time it takes for the landslide mass to flow into the area is s.

[0089] Furthermore, the continuity equation in step 5 is:

[0090] (25)

[0091] The momentum equations include the momentum equations in the x-direction and the momentum equations in the y-direction;

[0092] The momentum equation in the x-direction is:

[0093] (26)

[0094] The momentum equation in the y-direction is:

[0095] (27)

[0096] The sediment transport equation is:

[0097] (28)

[0098] The equation for riverbed deformation is:

[0099] (29)

[0100] In the formula, The water depth is in meters (m). , The velocity in the x and y directions are respectively, in m / s; The acceleration due to gravity is expressed in m / s². For source and sink terms, m / s; Lateral inflow over the slope, m / s; , The bottom slope is in the x and y directions; , The friction gradients are in the x and y directions; This refers to the volume concentration of sediment. The concentration of sediment flowing in laterally; The riverbed scour rate is expressed in m / s. The sediment settling rate is expressed in m / s. The elevation of the riverbed is in meters (m). The porosity of the riverbed;

[0101] The specific process for solving the water-sediment coupling control equations is as follows: spatial discretization is performed using the finite volume method, and the Riemann problem is solved using the HLLC or Roe scheme; the friction gradient is calculated separately for the x and y components using the Manning formula.

[0102] (30)

[0103] (31)

[0104] In the formula, The roughness coefficient is Manning's coefficient. .

[0105] Furthermore, the formula for calculating the time step of the CFL condition constraint in step 6 is as follows:

[0106] (33)

[0107] In the formula, CFL is the Courant number; Minimum grid size, m; , For grid cells The flow velocity in the x and y directions, in m / s; For grid cells The water depth, in meters; The acceleration due to gravity is expressed in m / s². Maximum wave speed, m / s;

[0108] The formula for calculating the time step of the mass conservation constraint is:

[0109] (34)

[0110] In the formula, This refers to the allowable percentage of quality variation. This is the minimum water depth threshold.

[0111] The landslide risk constraint time step is calculated more densely in high-risk areas, and the calculation formula is as follows:

[0112] (35)

[0113] In the formula, The baseline hydrological time step is s; Risk sensitivity coefficient; To adjust the index;

[0114] The formula for calculating the time step of riverbed deformation constraint is:

[0115] (36)

[0116] In the formula, This refers to the permissible proportion of riverbed changes; For reference, the riverbed thickness is measured in meters (m).

[0117] Furthermore, the specific process of performing watershed spatial decomposition using the topology-based parallel load balancing algorithm in step 7, and implementing distributed parallel computing through multi-process parallelism and GPU acceleration, is as follows:

[0118] A subdomain refers to the computational region allocated to a single parallel process after the watershed space is decomposed. This represents the subdomain allocated to the p-th process; first, define the computational load of the p-th subdomain. :

[0119] (39)

[0120] In the formula, Let p be the number of grid cells contained in the p-th subdomain; This represents the number of slope units; This represents the number of channel units; This represents the average number of substeps. , , These are the grid calculation weights, slope calculation weights, and gully calculation weights, respectively.

[0121] Define the cost of communication between subdomains :

[0122] (40)

[0123] In the formula, Let be the number of grid cells at the boundary between the p-th and q-th subdomains; Number of cross-boundary rivers; , This is the communication cost coefficient;

[0124] Then, the objective function for parallel load balancing optimization is constructed as follows:

[0125] (41)

[0126] In the formula, This represents the number of parallel processes. This is the subdomain allocated to the p-th process; For subdomain The boundary neighborhood; For communication - calculate the tradeoff coefficient;

[0127] To solve the above load balancing optimization objective function, a recursive decomposition algorithm based on binary tree hierarchy is used to divide the flow domain into parallel processes. The specific steps are as follows:

[0128] (a) Calculate the depth of the binary tree for:

[0129] (42)

[0130] (b) Starting from the root node, recursively divide according to the hierarchy:

[0131] (43)

[0132] (44)

[0133] (c) When partitioning at each level, adjust the boundaries to minimize the load difference between the left and right subdomains:

[0134] (45)

[0135] In the formula, The depth of the binary tree; For computing unit binary tree encoding; For subdomain Level; , These are the sets of left and right subdomains, respectively. The difference in load between the left and right subdomains; To allow for deviations in load balancing; Average subdomain load;

[0136] Computational Parallel Efficiency Evaluation Metrics for:

[0137] (46)

[0138] In the formula, The serial computation time is in seconds. The parallel computation time is in seconds (s). For subdomain load variance; Communication time, in seconds; For time calculation, s;

[0139] A dynamic load redistribution strategy is used to achieve runtime load balancing: when a load imbalance is detected ( Triggering dynamic adjustment of subdomain boundaries:

[0140] (47)

[0141] In the formula, , These are the subdomain boundaries before and after adjustment, respectively. The boundary adjustment amount is determined according to the proportion of the load difference between adjacent subdomains;

[0142] A pipelined parallel strategy based on river network confluence is adopted to achieve cross-level parallel computing:

[0143] (a) By topological level Grouping by river segment:

[0144] (48)

[0145] (b) River segments within the same level are calculated in parallel, while different levels are calculated sequentially from upstream to downstream;

[0146] (c) The pipeline speedup ratio is:

[0147] (49)

[0148] In the formula, This is the set of river segments in the Lth layer; The topological level at which river segment unit e resides; The speedup ratio of the production line; The maximum topology level; The number of the Lth layer river segment;

[0149] The water and sediment evolution is accelerated using GPUs, utilizing CUDA or OpenCL to achieve grid-level parallel computing, and employing a GPU thread block allocation strategy to map the grid to threads.

[0150] (50)

[0151] In the formula, Number of thread blocks; Number of grid cells; The number of threads per block.

[0152] The beneficial effects of the present invention are as follows:

[0153] 1) This invention constructs an integrated topological structure for five types of objects: slope, gully, node, small watershed, and river section. It achieves close coupling of hydrological calculation, landslide initiation judgment, and water and sediment evolution calculation, and can simulate the chain process of rainstorm-flash flood-landslide composite disaster.

[0154] 2) The slope-ditch coupling transport coefficient proposed in this invention takes into account topography, soil and vegetation factors, realizes the quantitative characterization of the transport efficiency of slope runoff to ditch, and improves the accuracy of slope-ditch material exchange simulation.

[0155] 3) The soil saturation-landslide risk index proposed in this invention integrates four dimensions: safety factor, soil saturation, topography, and rainfall, to achieve dynamic assessment of landslide risk. Compared with a single safety factor index, it has a better early warning lead time.

[0156] 4) The multiphysics adaptive time step strategy proposed in this invention comprehensively considers four types of constraints: CFL conditions, mass conservation, landslide risk, and riverbed deformation, achieving a balance between computational accuracy and efficiency.

[0157] 5) The parallel load balancing algorithm based on topology hierarchy proposed in this invention makes full use of the binary tree topology structure for spatial decomposition and load balancing, which improves the efficiency of distributed parallel computing and enhances the accuracy and timeliness of flash flood disaster early warning and forecast.

[0158] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments. Attached Figure Description

[0159] Figure 1 This is a schematic diagram of the method flow described in this invention. Detailed Implementation

[0160] This invention discloses a distributed simulation method for watershed composite flash flood disasters based on slope-gully topology, such as... Figure 1 As shown, the method includes the following steps:

[0161] Step 1: Slope-Valley Unit Extraction and Topology Construction: Based on the digital elevation model (DEM) data of the study area, the D8 flow direction algorithm is used to calculate the flow direction grid and the runoff accumulation grid. A catchment area threshold is set. (Value range: 0.05-0.5 km²) Extract the gully network and mark the grid cells with a cumulative runoff volume greater than a threshold as gully cells. Using the gully network as the boundary, divide the watershed into several slope cells. Each slope cell flows into the adjacent gully cell, thus constructing a slope-gully binary tree topology.

[0162] Nodes are the confluence or branching points of channels. Let the root node (watershed outlet) be coded as 0. Using the following binary tree coding rule, for nodes coded as... Nodes:

[0163] (1)

[0164] (2)

[0165] (3)

[0166] A relationship between slope units and channel units is established. Each channel unit is associated with its left and right bank slopes. The slope position is determined according to the left-hand rule in the direction of water flow, i.e., facing downstream, the left side is the left bank slope and the right side is the right bank slope. Based on the watershed-segment-node topological relationship of the CNFF model (China Mountain Flood Hydrological Model), an integrated topological structure of five types of objects—slope, channel, node, watershed, and segment—is formed.

[0167] Step 2, Hydrological Model Calculation: Perform hydrological model calculations, including rainfall interpolation, evapotranspiration calculation, soil moisture calculation, runoff generation calculation, slope runoff, and river network runoff evolution. All calculations are performed using methods known in the field. Specifically: rainfall interpolation uses the Thiessen polygon method, inverse distance weighted method, or Kriging interpolation to interpolate rainfall station data to each calculation unit to obtain areal rainfall; evapotranspiration is calculated using a three-layer evaporation model; soil moisture is calculated using three-layer soil moisture content; runoff generation uses a full-storage runoff model, dividing it into three water sources: surface runoff, interflow, and groundwater runoff; in slope runoff, surface runoff uses the unit hydrograph method, while interflow and groundwater runoff use the linear reservoir method; and river network runoff evolution uses the Muskingen method or the dynamic Muskingen method.

[0168] The output results calculated by the hydrological model include: Slope unit Soil saturation at time Slope unit outflow and small watershed outflow wait.

[0169] Among them, the Slope unit Soil saturation at time (Dimensionless, 0-1), determined by the ratio of soil moisture content to watershed storage capacity:

[0170] (4)

[0171] In the formula, For the first Slope unit Total soil moisture content at any given time (mm); The watershed storage capacity (mm) of this slope unit.

[0172] Slope unit outflow (m³ / s), for the first Slope unit The slope outflow at a given time is used as input for the lateral inflow calculation in step 3; the small watershed outlet flow rate. (m³ / s), used as the upstream inflow condition for the river network confluence evolution and the water and sediment evolution in step 5.

[0173] Step 3: Calculate the slope-gully coupling transport coefficient to quantify transport efficiency: Based on the soil saturation obtained in Step 2, calculate the slope-gully coupling transport coefficient to quantify the efficiency of runoff transmission from the slope to the gully. This coefficient comprehensively considers topographic, soil, and vegetation factors, and the formula is as follows:

[0174] (5)

[0175] In the formula, For the first Slope unit Slope-channel coupling transmission coefficient at time (dimensionless). For the first Terrain transfer factor (dimensionless) for each slope unit. For the first Slope unit Soil transport factor at time (dimensionless). For the first Vegetation transport factor (dimensionless) for each slope unit. The reference transmission coefficient (dimensionless, with a value of 0.8-1.0).

[0176] Among them, terrain transmission factor Reflecting the impact of slope geometry on transmission efficiency:

[0177] (6)

[0178] In the formula, For the first Average slope (°) of each slope unit; For reference slope (°); For the first Average flow length (m) of each slope unit; Reference process length (m); The slope influence index (dimensionless, ranging from 0.5 to 1.0).

[0179] Soil transport factors The dynamic impact of soil saturation state on transport efficiency:

[0180] (7)

[0181] In the formula, For the first Slope unit Soil saturation at time (dimensionless, 0-1); Reference saturation (dimensionless); For the first Soil saturated hydraulic conductivity (mm / h) of each slope unit; The maximum saturated hydraulic conductivity (mm / h); , These are the saturation influence index and the hydraulic conductivity influence index (both dimensionless).

[0182] Vegetation transport factors The effect of vegetation cover on transmission efficiency:

[0183] (8)

[0184] In the formula, For the first Vegetation cover of each slope unit (dimensionless, 0-1). For the first Leaf area index (dimensionless) of each slope unit. Leaf area index (dimensionless) is used as a reference.

[0185] The lateral inflow rate from the slope to the gully is then calculated as follows:

[0186] (9)

[0187] In the formula, The lateral inflow rate per unit channel length (m² / s). For the first Outflow rate (m³ / s) of each slope unit. For corresponding channel units The length (m).

[0188] Step 4: Calculate the Soil Saturation-Landslide Risk Index to assess landslide risk: First, define the Soil Saturation-Landslide Risk Index (LRI) to comprehensively assess the landslide risk of each slope unit.

[0189] (10)

[0190] In the formula, For the first Slope unit Soil saturation at time - landslide risk index (dimensionless, 0-1). , , , These are the safety factor weight, saturation weight, terrain weight, and rainfall weight, respectively, satisfying... ; For the safety factor risk component; This represents the risk component of soil saturation. For terrain risk components; Accumulated risk components for rainfall (all risk components are dimensionless).

[0191] The formula for calculating the risk component of the safety factor is as follows:

[0192] (11)

[0193] In the formula, For the first Slope unit Safety factor at any given time (dimensionless, greater than 1 indicates stability, less than or equal to 1 indicates instability); The critical safety factor (dimensionless). The minimum safety factor threshold (dimensionless). It is a non-linear exponent (dimensionless).

[0194] Safety factor Calculations based on the Mohr-Coulomb criterion for infinite slopes:

[0195] (12)

[0196] In the formula, For the first Slope unit Effective cohesion at any given time (kPa); The depth of the sliding surface (m); Slope (°); The unit weight of saturated soil (kN / m³); The specific weight of water (kN / m³). For the first Slope unit The effective internal friction angle (°) at time t. The soil strength parameter decreases with soil saturation:

[0197] (13)

[0198] (14)

[0199] In the formula, Effective cohesive strength (kPa) in dry state; The effective internal friction angle (°) in the dry state; , These are the cohesion attenuation coefficient and the internal friction angle attenuation coefficient, respectively (both dimensionless).

[0200] The formula for calculating the soil saturation risk component is:

[0201] (15)

[0202] In the formula, Critical saturation (dimensionless); It is a non-linear exponent (dimensionless). The sensitivity coefficient for the rate of change of saturation (h); For the first The rate of change of saturation of each slope unit (1 / h).

[0203] The formula for calculating the terrain risk component is:

[0204] (16)

[0205] In the formula, For the first The tangent of the slope of a slope unit; , These are the tangents of the minimum and maximum slopes at which the landslide occurred; For the first Relative elevation (m) of each slope unit; Maximum relative elevation (m); This is a nonlinear topographic index (dimensionless).

[0206] The formula for calculating the cumulative risk component of rainfall is:

[0207] (17)

[0208] In the formula, Critical cumulative rainfall (mm, determined by region); For the first Slope unit The effective cumulative rainfall (mm) before the specified time is calculated using exponential decay:

[0209] (18)

[0210] In the formula, For the first Each slope unit in Rainfall at any given time (mm); The duration of the initial impact (h); The decay time constant is (h).

[0211] Then, based on the risk index, a two-level judgment of landslide risk is made, divided into two levels: landslide risk warning and landslide initiation.

[0212] (1) Landslide risk warning: When the risk index of a slope unit reaches the critical risk index, a landslide risk warning is issued. The warning conditions are as follows:

[0213] (19)

[0214] In the formula, This is a critical risk index (dimensionless). Due to the risk index... This approach integrates four dimensions: safety factor, soil saturation, topography, and prior rainfall. It allows the system to exceed the warning threshold before the soil becomes saturated and the safety factor drops to an unstable level. Therefore, compared to relying solely on the safety factor, this method is more effective. A single criterion can issue early warnings and provide a longer lead time for warnings.

[0215] (2) Landslide initiation judgment: On the basis of meeting the early warning conditions, when the safety factor of the slope unit further decreases to the instability value, it is determined that the slope unit has started a landslide (instability). The landslide initiation conditions are as follows:

[0216] (20)

[0217] Calculate the landslide volume, let Inversely calculate the instability thickness :

[0218] (twenty one)

[0219] The landslide volume is:

[0220] (twenty two)

[0221] In the formula, For the first Instability thickness (m) of each slope element; For the first Area of ​​each slope unit (m²); , is the shape correction factor (dimensionless, considering the non-uniform thickness of the landslide body). For the first Landslide volume (m³) of each slope unit.

[0222] Then calculate the unit width flow rate of the landslide entering the gully. , This serves as the lateral inflow boundary condition for the water and sediment evolution model. Assume the landslide body moves at a velocity... Angle along the direction of entry into the ditch During the import of history Inner width When flowing into a channel, the unit width flow rate is calculated using the following formula:

[0223] (twenty three)

[0224] (twenty four)

[0225] In the formula, , These are the unit width flow rates (m² / s) of the landslide inflow in the x and y directions, respectively. The volume of the landslide is (m³). The velocity of the sliding body as it enters the ditch (m / s); The angle of entry into the ditch (°); The width of the landslide entering the gully (m); The duration (s) of the landslide mass inflow.

[0226] Step 5, Calculation of water and sediment evolution model: Perform water and sediment evolution calculation, receive the lateral inflow from the slope coupling transmission in Step 3 and the landslide inflow from Step 4, and solve the water and sediment coupling control equation set.

[0227] Specifically, based on the outflow of the small watershed Lateral inflow of coupled transmission and landslide inflow , Establish a set of water-sediment coupling control equations, including:

[0228] Continuity equation:

[0229] (25)

[0230] Momentum equation in the x-direction:

[0231] (26)

[0232] The momentum equation in the y-direction is:

[0233] (27)

[0234] Sediment transport equation:

[0235] (28)

[0236] Riverbed deformation equation:

[0237] (29)

[0238] In the formula, Water depth (m); , The values ​​are the flow velocities (m / s) in the x and y directions, respectively. The acceleration due to gravity (m / s²) Source and sink terms (m / s); Lateral inflow over the slope (m / s); , The bottom slope in the x and y directions (dimensionless); , The friction gradient in the x and y directions is dimensionless. The volume concentration of sediment (dimensionless). The lateral inflow sediment concentration (dimensionless). The riverbed scour rate (m / s); The sediment settling rate (m / s); Riverbed elevation (m); The porosity of the riverbed is dimensionless.

[0239] Spatial discretization is performed using the finite volume method, and the Riemann problem is solved using the HLLC or Roe scheme; the friction gradient is calculated separately for the x and y components using the Manning formula.

[0240] (30)

[0241] (31)

[0242] In the formula, Manning's roughness coefficient ( ).

[0243] Finally, output the water depth of the channel unit. Water level Flow rate Unit width flow Riverbed changes and sediment concentration .

[0244] Step 6: Optimize computational efficiency using a multiphysics adaptive time step strategy: Considering the coupling characteristics of hydrology, landslide, and sediment multiphysics, a comprehensive adaptive time step selection strategy is proposed to coordinate the time progression of each physical process. The multiphysics adaptive time step is:

[0245] (32)

[0246] In the formula, The adaptive time step (s) for multiphysics; The time step (s) for CFL condition constraints; The time step (s) is a mass conservation constraint. The time step (s) is the constraint for landslide risk. The time step (s) is the constraint time step for riverbed deformation. The safety factor (dimensionless).

[0247] The formula for calculating the time step of the CFL condition constraint is as follows:

[0248] (33)

[0249] In the formula, CFL is the Courant number (dimensionless, taken as 0.5-0.9); Minimum grid size (m); , For grid cells Flow velocity in the x and y directions (m / s); For grid cells Water depth (m); The acceleration due to gravity (m / s²) The maximum wave speed is (m / s).

[0250] The mass conservation constraint time step is calculated using the following formula and is used to limit the water depth variation within a single time step, ensuring the accuracy of mass conservation:

[0251] (34)

[0252] In the formula, The permissible proportion of mass variation (dimensionless). This is the minimum water depth threshold (m).

[0253] Landslide risk constraints time step, with more detailed calculations in high-risk areas:

[0254] (35)

[0255] In the formula, The baseline hydrological time step (s); The risk sensitivity coefficient (dimensionless). The adjustment index is dimensionless. When the maximum risk index approaches 1, the time step is significantly reduced to capture the moment of landslide initiation.

[0256] The time step for riverbed deformation constraint is calculated using the following formula, which is used to limit the magnitude of riverbed change within a single time step:

[0257] (36)

[0258] In the formula, The allowable proportion of riverbed variation (dimensionless); For reference, the riverbed thickness is (m).

[0259] By dynamically adjusting the time step based on local truncation error estimation, computational efficiency can be improved while maintaining accuracy.

[0260] (37)

[0261] In the formula, This is the allowable error (dimensionless). The current step estimation error (dimensionless); The numerical format order (dimensionless).

[0262] The time step nesting relationship between step 2 (hydrological calculation) and step 5 (water and sediment evolution calculation) is as follows:

[0263] (38)

[0264] In the formula, The nesting multiple of water and sand steps (dimensionless); This refers to the hydrological time step (s) (i.e., the time step set during the operation of the hydrological model). The formula represents the time step at each hydrological time step. Internally, the calculation of water and sediment evolution is executed. Calculate the adaptive substep size.

[0265] Step 7: Implement distributed parallel computing based on topology-level parallel load balancing algorithm: Decompose the watershed space based on topology-level parallel load balancing algorithm, and implement distributed parallel computing using multi-process parallelism and GPU acceleration.

[0266] Here, a subdomain refers to the computational region allocated to a single parallel process after the watershed space is decomposed. This represents the subdomain allocated to the p-th process. First, define the computational load of the p-th subdomain. :

[0267] (39)

[0268] In the formula, Let p be the number of grid cells contained in the p-th subdomain; This represents the number of slope units; This represents the number of channel units; The average number of substeps (dimensionless). , , These are the grid calculation weights, slope calculation weights, and gully calculation weights (dimensionless, calibrated based on actual tests).

[0269] Define the cost of communication between subdomains :

[0270] (40)

[0271] In the formula, Let be the number of grid cells at the boundary between the p-th and q-th subdomains; Number of cross-boundary rivers; , This is the communication cost coefficient (dimensionless).

[0272] Then, the objective function for parallel load balancing optimization is constructed as follows:

[0273] (41)

[0274] In the formula, This represents the number of parallel processes. This is the subdomain allocated to the p-th process; For subdomain The boundary neighborhood; For communication - calculate the tradeoff coefficient (dimensionless).

[0275] To solve the above load balancing optimization objective function, this invention uses a recursive decomposition algorithm based on binary tree hierarchy to divide the flow domain into parallel processes. The specific steps are as follows:

[0276] (a) Calculate the depth of the binary tree for:

[0277] (42)

[0278] (b) Starting from the root node, recursively divide according to the hierarchy:

[0279] (43)

[0280] (44)

[0281] (c) When partitioning at each level, adjust the boundaries to minimize the load difference between the left and right subdomains:

[0282] (45)

[0283] In the formula, The depth of the binary tree (dimensionless); For computing unit binary tree encoding; For subdomain Level; , These are the sets of left and right subdomains, respectively. The difference in load between the left and right subdomains; The allowable deviation for load balancing (dimensionless). This represents the average subdomain load.

[0284] Computational Parallel Efficiency Evaluation Metrics for:

[0285] (46)

[0286] In the formula, The serial computation time is in seconds. The parallel computing time (s); For subdomain load variance; Communication time (s); To calculate the time (s).

[0287] Furthermore, a dynamic load redistribution strategy is adopted to achieve runtime load balancing: when a load imbalance is detected ( Triggering dynamic adjustment of subdomain boundaries:

[0288] (47)

[0289] In the formula, , These are the subdomain boundaries before and after adjustment, respectively. The boundary adjustment amount is determined according to the proportion of the load difference between adjacent subdomains, and boundary units with lower communication costs are adjusted first.

[0290] A pipelined parallel strategy based on river network confluence is adopted to achieve cross-level parallel computing:

[0291] (a) By topological level Grouping by river segment:

[0292] (48)

[0293] (b) River segments within the same level are calculated in parallel, while different levels are calculated sequentially from upstream to downstream;

[0294] (c) The pipeline speedup ratio is:

[0295] (49)

[0296] In the formula, This is the set of river segments in the Lth layer; The topological level at which river segment unit e resides; The flow rate is dimensionless. The maximum topology level; This represents the number of the Lth river segment.

[0297] The water and sediment evolution is accelerated using GPUs, utilizing CUDA or OpenCL to achieve grid-level parallel computing, and employing a GPU thread block allocation strategy to map the grid to threads.

[0298] (50)

[0299] In the formula, Number of thread blocks; Number of grid cells; The number of threads per block.

[0300] Example 1

[0301] This embodiment is an application example of the above method. This embodiment takes a mountainous watershed as the study area, with a watershed area of ​​approximately 200 km², characterized by large topographic relief, good vegetation cover, and historical occurrences of flash floods and debris flows.

[0302] This embodiment uses 30 m resolution DEM data and sets a catchment area threshold of 0.1 km² to extract the gully network. A total of 156 slope units, 89 gully units, and 45 confluence nodes were extracted. A binary tree coding method was used to establish the topological relationships, with the watershed outlet coded as 0 and recursively coded upstream. Each gully unit is associated with the left and right slopes, forming a slope-gully association table.

[0303] The parameters for the slope-channel coupling transmission coefficient in this embodiment are shown in Table 1.

[0304] Table 1. Slope-Ditch Coupling Transmission Coefficient Parameters

[0305]

[0306] The hydrological model parameter settings for this embodiment are shown in Table 2:

[0307] Table 2 Hydrological Model Parameters

[0308]

[0309] Among them, the watershed water storage capacity Wm(j) is determined according to the soil type and thickness of each slope unit, and its value ranges from 80 to 150 mm. In this embodiment, a representative value of 120 mm is uniformly used. The total soil moisture content W(j,t) is a state variable that is dynamically updated during the calculation process. In this embodiment, the initial condition is given as the initial soil saturation Sr(j,0)=0.45 for each slope unit, that is, W(j,0)=0.45·Wm(j)=54 mm, and the saturation calculation is started from this point.

[0310] The soil saturation-landslide risk index parameter settings in this embodiment are shown in Table 3:

[0311] Table 3 Soil saturation-landslide risk index parameters

[0312]

[0313] The soil mechanical parameters set in this embodiment are shown in Table 4:

[0314] Table 4 Soil mechanical parameters

[0315]

[0316] The parameters of the water and sediment evolution model in this embodiment are set as follows: the channel grid resolution is 2 m, the Manning roughness coefficient n is 0.035-0.05 (determined according to the channel type), and the riverbed porosity p = 0.4.

[0317] The adaptive time step parameter settings for this embodiment are shown in Table 5:

[0318] Table 5 Adaptive Time Step Parameters

[0319]

[0320] The parallel load balancing parameter settings for this embodiment are shown in Table 6:

[0321] Table 6 Parallel Load Balancing Parameters

[0322]

[0323] This embodiment sets the simulation duration to 24 hours and inputs a rainstorm event (maximum hourly rainfall intensity 50 mm / h, cumulative rainfall 180 mm) to simulate a complex flash flood disaster. The simulation results include:

[0324] 1) Hydrological response: The peak discharge at the basin outlet reached 1100 m³ / s, with a peak delay of approximately 3 hours;

[0325] 2) Coupled transport coefficient: The average coupled transport coefficient increased from 0.21 to 0.59 during the rainstorm, reflecting a significant increase in slope-ditch transport efficiency after soil saturation;

[0326] 3) Landslide risk early warning: Approximately 45 minutes before the actual landslide, the LRI of 8 slope units exceeded the warning threshold of 0.7, providing a greater lead time compared to a single landslide. The indicator is set approximately 30 minutes in advance;

[0327] 4) Landslide Initiation: A total of 12 slope units became unstable, mainly distributed in areas with slopes greater than 25°, with a total landslide volume of approximately [missing information]. m³;

[0328] 5) Water and sediment evolution: The maximum water depth in the gully reaches 2.8 m, the maximum flow velocity reaches 4.5 m / s, the peak sediment concentration reaches 15%, and the main channel riverbed is incised down by 0.3-0.8 m;

[0329] 6) Adaptive time step: The time step for water and sediment evolution is adaptively adjusted within the range of 0.05-0.8 s, and automatically densified to 0.05-0.1 s during the high-risk period of landslides;

[0330] 7) Parallel efficiency: The parallel efficiency of 8 processes reaches 0.82, and the load imbalance is controlled within 0.08;

[0331] 8) Computational efficiency: Using an 8-core CPU and GPU, it takes about 45 seconds to complete a 24-hour simulation, which meets the requirements for real-time early warning.

[0332] Finally, it should be noted that the above description is only used to illustrate the technical solution of the present invention and not to limit it. Although the present invention has been described in detail with reference to the preferred arrangement, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solution of the present invention without departing from the spirit and scope of the technical solution of the present invention.

Claims

1. A distributed simulation method for watershed composite flash flood disasters based on slope-gully topology, characterized in that, The method includes the following steps: Step 1: Slope-Valley Unit Extraction and Topology Construction: Based on the digital elevation model data of the study area, the D8 flow direction algorithm is used to calculate the flow direction grid and the runoff accumulation grid; a catchment area threshold is set. The process involves extracting the channel network and marking grids with cumulative runoff exceeding a threshold as channel units. Using the channel network as boundaries, the watershed is divided into several slope units, with each slope unit flowing into an adjacent channel unit, thus constructing a slope-channel binary tree topology. Nodes represent the confluence or branching points of channels. The relationship between slope units and channel units is established. Based on the small watershed-river segment-node topological relationship of the Chinese mountain flood hydrological model, an integrated topological structure is formed, encompassing five categories of objects: slope, channel, node, small watershed, and river segment. Step 2, Hydrological Model Calculation: Perform hydrological model calculations, including rainfall interpolation, evapotranspiration calculation, soil moisture calculation, runoff generation calculation, slope runoff, and river network runoff evolution; output the results of the hydrological model calculations, including: Slope unit Soil saturation at time Slope unit outflow and small watershed outflow ; Step 3: Calculate the slope-gully coupling transport coefficient to quantify transport efficiency: Based on the soil saturation obtained in Step 2, calculate the slope-gully coupling transport coefficient to quantify the efficiency of runoff transmission from the slope to the gully. This coefficient comprehensively considers topographic, soil, and vegetation factors, and the formula is as follows: （5） In the formula, For the first Slope unit The slope-channel coupling transmission coefficient at time t; For the first Topographic transfer factor for each slope unit; For the first Slope unit Soil transport factors at any given time; For the first Vegetation transport factor of each slope unit; The reference transmission coefficient; The lateral inflow rate from the slope to the gully is then calculated as follows: （9） In the formula, The lateral inflow rate per unit channel length is expressed in m² / s. For the first The outflow rate of each slope unit is m³ / s; For corresponding channel units The length, in meters; Step 4: Calculate the Soil Saturation-Landslide Risk Index to assess landslide risk: Define the Soil Saturation-Landslide Risk Index (LRI) to comprehensively assess the landslide risk of each slope unit. （10） In the formula, For the first Slope unit Soil saturation at any given time - landslide risk index; , , , These are the safety factor weight, saturation weight, terrain weight, and rainfall weight, respectively, satisfying... ; For the safety factor risk component; This represents the risk component of soil saturation. For terrain risk components; Accumulated risk component for rainfall; Then, based on the soil saturation-landslide risk index (LRI), a two-level judgment of landslide risk was performed, divided into two levels: landslide risk warning and landslide initiation. Finally, the unit width flow rate of the landslide mass entering the gully was calculated. , ; Step 5: Calculation of water and sediment evolution model: based on the outlet flow of the small watershed. Lateral inflow of coupled transmission and landslide inflow , A set of water-sediment coupling control equations was established, including the continuity equation, momentum equation, sediment transport equation, and riverbed deformation equation. The set of water-sediment coupling control equations was solved, and finally, the water depth of the channel unit was output. Water level Flow rate Unit width flow Riverbed changes and sediment concentration ; Step 6: Optimize computational efficiency using a multiphysics adaptive time step strategy: Considering the coupling characteristics of hydrology, landslide, and sediment multiphysics, a comprehensive adaptive time step selection strategy is proposed to coordinate the time progression of each physical process; the multiphysics adaptive time step is: （32） In the formula, The adaptive time step for multiphysics is s; Let be the time step of the CFL condition constraint, in seconds; The time step is defined as s, which is a constraint on mass conservation. The time step for landslide risk constraints is s; Let be the time step for constraining riverbed deformation, in seconds; For safety factor; By dynamically adjusting the time step based on local truncation error estimation, computational efficiency can be improved while maintaining accuracy. （37） In the formula, Tolerance; The current step is used to estimate the error; The order is in numerical format; The time step nesting relationship between hydrological calculations and sediment evolution calculations is as follows: （38） In the formula, The nesting multiple of water and sand steps; The hydrological time step is the time step set during the operation of the hydrological model, expressed in seconds; this formula represents the time step at each hydrological time step. Internally, the calculation of water and sediment evolution is executed. Calculation of adaptive substep size; Step 7: Implement distributed parallel computing based on topology-level parallel load balancing algorithm: Decompose the watershed space based on topology-level parallel load balancing algorithm, and implement distributed parallel computing using multi-process parallelism and GPU acceleration.

2. The distributed simulation method for watershed composite flash flood disasters based on slope-gully topology as described in claim 1, characterized in that, The nodes in step 1 use a binary tree encoding rule. Let the root node be encoded as 0. Then, for the node encoded as i: （1） （2） （3） The relationship between slope units and channel units is established as follows: each channel unit is associated with its left bank slope and right bank slope. The slope position is determined according to the left-right rule in the direction of water flow, that is, facing downstream, the left side is the left bank slope and the right side is the right bank slope.

3. The distributed simulation method for watershed composite flash flood disasters based on slope-gully topology as described in claim 1, characterized in that, Step 2 Slope unit Soil saturation at time Determined by the ratio of soil moisture content to watershed storage capacity: （4） In the formula, For the first Slope unit Total soil moisture content at any given time (mm); The watershed storage capacity (mm) of this slope unit.

4. The distributed simulation method for watershed composite flash flood disasters based on slope-gully topology as described in claim 1, characterized in that, Terrain transfer factor in step 3 Reflecting the impact of slope geometry on transmission efficiency: （6） In the formula, For the first The average slope of each slope unit; For reference slope; For the first Average flow length of each slope unit, in meters; The reference process length is in meters (m). The slope influence index; Soil transport factors The dynamic impact of soil saturation state on transport efficiency: （7） In the formula, For reference saturation; For the first Soil saturated hydraulic conductivity of each slope unit, mm / h; The maximum saturated hydraulic conductivity is given in mm / h. , These are the saturation influence index and the hydraulic conductivity influence index, respectively. Vegetation transport factors The effect of vegetation cover on transmission efficiency: （8） In the formula, For the first Vegetation coverage of each slope unit; For the first Leaf area index of each slope unit; Leaf area index for reference.

5. The distributed simulation method for watershed composite flash flood disasters based on slope-gully topology according to claim 1, characterized in that, The formula for calculating the risk component of the safety factor in step 4 is: （11） In the formula, For the first Slope unit The safety factor at any given time is greater than 1, indicating stability, and less than or equal to 1, indicating instability. This is the critical safety factor; The minimum safety factor threshold; It is a non-linear exponent; Safety factor Calculations based on the Mohr-Coulomb criterion for infinite slopes: （12） In the formula, For the first Slope unit Effective cohesion at any given time, kPa; Let m be the depth of the sliding surface. Slope; The unit weight of saturated soil is kN / m³. The specific weight of water is kN / m³. For the first Slope unit The effective internal friction angle at any given time; where the soil strength parameter decreases with soil saturation: （13） （14） In the formula, Effective cohesive strength in dry state, kPa; This is the effective internal friction angle in the dry state; , These are the cohesion attenuation coefficient and the internal friction angle attenuation coefficient, respectively. The formula for calculating the soil saturation risk component is: （15） In the formula, Critical saturation; It is a non-linear exponent; h is the sensitivity coefficient for the rate of change of saturation. For the first The rate of change of saturation of each slope unit, 1 / h; The formula for calculating the terrain risk component is: （16） In the formula, For the first The tangent of the slope of a slope unit; , These are the tangents of the minimum and maximum slopes at which the landslide occurred; For the first The relative elevation of each slope unit, in meters; The maximum relative elevation is in meters (m). This is a terrain nonlinearity index; The formula for calculating the cumulative risk component of rainfall is: （17） In the formula, The critical cumulative rainfall is expressed in mm. For the first Slope unit The effective cumulative rainfall before the specified time, in mm, is calculated using exponential decay: （18） In the formula, For the first Each slope unit in Rainfall at any given time, in mm; The duration of the initial impact is in hours (h). h is the decay time constant.

6. The distributed simulation method for watershed composite flash flood disasters based on slope-gully topology according to claim 5, characterized in that, The landslide risk warning in step 4 is as follows: a landslide risk warning is issued when the risk index of a slope unit reaches the critical risk index. The warning conditions are as follows: （19） In the formula, This is a critical risk index; Landslide initiation is defined as follows: Based on the met early warning conditions, when the safety factor of a slope element further decreases to an instability value, the landslide initiation is determined to have occurred in that slope element. The landslide initiation conditions are as follows: （20） Calculate the landslide volume, let Inversely calculate the instability thickness : （21） The landslide volume is: （22） In the formula, For the first The instability thickness of each slope element, in meters; For the first The area of ​​each slope unit, in m²; This is the shape correction factor; For the first The landslide volume of each slope unit, in m³; Calculate the unit width flow rate of the landslide mass entering the gully. , The specific process is as follows: Assume the landslide body moves at a speed of... Angle along the direction of entry into the ditch During the import of history Inner width When flowing into a channel, the unit width flow rate is calculated using the following formula: （23） （24） In the formula, , , respectively, are the unit width flow rates of the landslide inflow in the x and y directions, in m² / s; Let m be the volume of the landslide. The velocity of the sliding body as it enters the ditch is expressed in m / s. The angle of entry into the ditch; The width of the landslide entering the gully, in meters; The time it takes for the landslide mass to flow into the area is s.

7. The distributed simulation method for watershed composite flash flood disasters based on slope-gully topology according to claim 1, characterized in that, The continuity equation in step 5 is: （25） The momentum equations include the momentum equations in the x-direction and the momentum equations in the y-direction; The momentum equation in the x-direction is: （26） The momentum equation in the y-direction is: （27） The sediment transport equation is: （28） The equation for riverbed deformation is: （29） In the formula, The water depth is in meters (m). , The velocity in the x and y directions are respectively, in m / s; The acceleration due to gravity is expressed in m / s². For source and sink terms, m / s; Lateral inflow over the slope, m / s; , The bottom slope is in the x and y directions; , The friction gradients are in the x and y directions; This refers to the volume concentration of sediment. The concentration of sediment flowing in laterally; The riverbed scour rate is expressed in m / s. The sediment settling rate is expressed in m / s. The elevation of the riverbed is in meters (m). The porosity of the riverbed; The specific process for solving the water-sediment coupling control equations is as follows: spatial discretization is performed using the finite volume method, and the Riemann problem is solved using the HLLC or Roe scheme; the friction gradient is calculated separately for the x and y components using the Manning formula. （30） （31） In the formula, The roughness coefficient is Manning's coefficient. .

8. The distributed simulation method for watershed composite flash flood disasters based on slope-gully topology according to claim 1, characterized in that, The formula for calculating the time step of the CFL condition constraint in step 6 is: （33） In the formula, CFL is the Courant number; Minimum grid size, m; , For grid cells The flow velocity in the x and y directions, in m / s; For grid cells The water depth, in meters; The acceleration due to gravity is expressed in m / s². Maximum wave speed, m / s; The formula for calculating the time step of the mass conservation constraint is: （34） In the formula, This refers to the allowable percentage of quality variation. This is the minimum water depth threshold; The landslide risk constraint time step is calculated more densely in high-risk areas, and the calculation formula is as follows: （35） In the formula, The baseline hydrological time step is s; Risk sensitivity coefficient; To adjust the index; The formula for calculating the time step of riverbed deformation constraint is: （36） In the formula, This refers to the permissible proportion of riverbed changes; For reference, the riverbed thickness is measured in meters (m).

9. A distributed simulation method for watershed composite flash flood disasters based on slope-gully topology as described in claim 1, characterized in that, The specific process of performing watershed spatial decomposition using the topology-based parallel load balancing algorithm in step 7, and implementing distributed parallel computing through multi-process parallelism and GPU acceleration, is as follows: A subdomain refers to the computational region allocated to a single parallel process after the watershed space is decomposed. This represents the subdomain allocated to the p-th process; first, define the computational load of the p-th subdomain. : （39） In the formula, Let p be the number of grid cells contained in the p-th subdomain; This represents the number of slope units; This represents the number of channel units; This represents the average number of substeps. , , These are the grid calculation weights, slope calculation weights, and gully calculation weights, respectively. Define the cost of communication between subdomains : （40） In the formula, Let be the number of grid cells at the boundary between the p-th and q-th subdomains; Number of cross-boundary rivers; , This is the communication cost coefficient; Then, the objective function for parallel load balancing optimization is constructed as follows: （41） In the formula, This represents the number of parallel processes. This is the subdomain allocated to the p-th process; For subdomain The boundary neighborhood; For communication - calculate the tradeoff coefficient; To solve the above load balancing optimization objective function, a recursive decomposition algorithm based on binary tree hierarchy is used to divide the flow domain into parallel processes. The specific steps are as follows: (a) Calculate the depth of the binary tree for: （42） (b) Starting from the root node, recursively divide according to the hierarchy: （43） （44） (c) When partitioning at each level, adjust the boundaries to minimize the load difference between the left and right subdomains: （45） In the formula, The depth of the binary tree; For computing unit binary tree encoding; For subdomain Level; , These are the sets of left and right subdomains, respectively. The difference in load between the left and right subdomains; To allow for deviations in load balancing; Average subdomain load; Computational Parallel Efficiency Evaluation Metrics for: （46） In the formula, The serial computation time is in seconds. The parallel computation time is in seconds (s). For subdomain load variance; Communication time, in seconds; For time calculation, s; A dynamic load redistribution strategy is used to achieve runtime load balancing: when a load imbalance is detected ( Triggering dynamic adjustment of subdomain boundaries: （47） In the formula, , These are the subdomain boundaries before and after adjustment, respectively. The boundary adjustment amount is determined according to the proportion of the load difference between adjacent subdomains; A pipelined parallel strategy based on river network confluence is adopted to achieve cross-level parallel computing: (a) By topological level Grouping by river segment: （48） (b) River segments within the same level are calculated in parallel, while different levels are calculated sequentially from upstream to downstream; (c) The pipeline speedup ratio is: （49） In the formula, This is the set of river segments in the Lth layer; The topological level at which river segment unit e resides; The speedup ratio of the production line; The maximum topology level; The number of the Lth layer river segment; The water and sediment evolution is accelerated using GPUs, utilizing CUDA or OpenCL to achieve grid-level parallel computing, and employing a GPU thread block allocation strategy to map the grid to threads. （50） In the formula, Number of thread blocks; Number of grid cells; The number of threads per block.

Images