Layered three-dimensional hydrological-hydrodynamic full-coupling simulation method for extra-small watershed

Abstract

The invention relates to the technical field of hydrological simulation, and particularly discloses an extra-small watershed layered three-dimensional hydrological-hydrodynamic full-coupling simulation method, which comprises the following steps: collecting multi-source data to construct a three-dimensional space discrete grid and dividing initial time units; the control equation set is converted into an algebraic system based on a Legendre polynomial and a Gaussian-Loabar node; adaptively adjusting a time unit length and a spectral basis function order according to the local time gradient change rate; a convergence spectral coefficient is obtained through spectral delay correction iteration solution; and finally, forcibly executing mass conservation constraint, controlling step recalculation according to a local time error, and gradually outputting a flow hydrograph and water depth and infiltration rate distribution according to time unit endpoints. According to the method, the dual self-adaption of the time step length and the polynomial order is realized, the numerical oscillation of the infiltration rate caused by the sudden change of the ponding depth can be effectively inhibited, the mass conservation is ensured, and the stability and the efficiency of the rainstorm flood simulation of the extra-small watershed are improved.

Description

Technical Field

[0001] This invention relates to the field of hydrological simulation technology, specifically to a layered, three-dimensional, fully coupled hydrological-hydraulic simulation method for very small watersheds. Background Technology

[0002] Small watersheds (closed units with an area of ​​a few hectares to a few square kilometers) are characterized by short runoff times, rapid rise and fall of flood peaks, and intense interaction between soil infiltration and surface runoff under short-duration heavy rainfall. Traditional hydrological models often use empirical runoff generation and runoff calculations (such as inference formulas and unit hydrographs) or sequential coupling methods (calculating runoff generation first and then calculating runoff), which are difficult to reflect the real-time two-way feedback between surface runoff, unsaturated soil infiltration, and shallow groundwater.

[0003] In the stratified three-dimensional hydrological-hydrodynamic fully coupled simulation of small watersheds, when torrential rain causes the surface water depth to jump dramatically from zero centimeters to several centimeters within minutes, traditional low-order time discretization methods based on fixed time steps or simple Coulomb number adjustments (such as backward Euler schemes) cannot smoothly handle such abrupt changes in water depth. This leads to non-physical sawtooth oscillations in infiltration rate calculations, characterized by steep drops followed by rebounds, which in turn causes divergence in the Newton iteration of the governing equations or spurious flood peak fluctuations in the simulation results. While existing high-order time schemes can improve accuracy, they lack adaptive and coordinated control of the time unit length and polynomial order, making it difficult to balance stability with computational efficiency. This invention introduces the time-domain spectral element method combined with a dual adaptive strategy of step size and order, and employs spectral delay correction iteration to achieve spectral accuracy tracking of the water depth abrupt change process, eliminating the infiltration rate oscillation phenomenon at the root of time discretization. Summary of the Invention

[0004] The purpose of this invention is to provide a layered, three-dimensional, fully coupled hydrological-hydraulic simulation method for very small watersheds to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions: A layered, three-dimensional, fully coupled hydrological-hydraulic simulation method for very small watersheds includes the following steps: S1 collects digital elevation data, soil vertical layer thickness data, land use type data, and time series rainfall data for the small watershed. Based on the collected data, a three-dimensional spatial discrete grid is constructed, and the total simulation duration is divided into multiple initial time units. The spatial discrete grid and initial time unit division results are output. S2 takes a three-dimensional spatial discrete grid and an initial time unit as input. Within each time unit, it uses a higher-order spectral basis function to approximate the temporal changes of surface water depth, soil pressure head, and groundwater level state variables. It also transforms the set of control equations describing surface water flow, soil water movement, and groundwater movement into a weak integral form based on higher-order spectral basis functions in the spatiotemporal domain. The output is an algebraic system containing unknown spectral coefficients in all time units. S3 takes the algebraic system as input, and dynamically adjusts the order of the higher-order spectral basis function and the length of the time unit based on the local time gradient change rate of the state variables of the previous time unit solution, generating an adaptive time discretization scheme, and outputting the adaptive time unit partitioning and corresponding order. S4 takes the adaptive time discretization scheme and algebraic system as input and uses spectral delay correction iterative solution: first, a low-order time scheme is used to obtain the predicted solution on the coarse time grid, then the residual is calculated using the spectral integral kernel and the predicted solution is corrected one by one until the residual meets the convergence threshold, and the convergence spectral coefficients of the state variables in each time unit are output. S5 takes the convergence spectrum coefficient as input, enforces the mass conservation integral constraint at the endpoint of each time unit, and accepts the current solution or recalculates by reducing the step size according to the local time error control, and advances to the next time unit. It outputs the flow process line of the outlet section of the small watershed and the spatiotemporal distribution of water depth and infiltration rate of the whole watershed successively according to the endpoint of the time unit.

[0006] As a further aspect of the present invention: the construction of the three-dimensional spatial discrete mesh specifically includes: A surface unstructured triangular grid is generated using digital elevation data. Then, based on the soil vertical layer thickness data, vertical layer nodes are inserted below each surface grid node according to the layer thickness ratio, so that the thin layer area is automatically densified. Land use type data is used to identify permeable and impermeable boundaries, and the surface grid is densified along the boundary lines. Based on the cumulative rainfall distribution of time series rainfall data, potential runoff paths are extracted, and the surface grid is linearly densified along the paths to obtain a three-dimensional spatial discrete grid.

[0007] As a further aspect of the present invention: dividing the total simulation duration into multiple initial time units specifically includes: The rate of change of rainfall intensity between adjacent moments is calculated based on time series rainfall data, and moments when the rate of change of rainfall intensity exceeds a preset threshold are set as candidate breakpoints. The length of the watershed confluence path is extracted using digital elevation data, and the water propagation time is estimated by combining it with surface roughness. The time interval between adjacent candidate breakpoints that is longer than the propagation time is taken as an initial time unit. During periods of stable rainfall intensity, the time units are divided into equal-length segments to ensure that each initial time unit contains at most one point of sudden change in rainfall intensity, and the division results are output.

[0008] As a further aspect of the present invention: the output process of the algebraic system specifically includes: Each time unit is mapped to a standard time interval, Legendre polynomials are selected as spectral basis functions within the standard time interval, and Gauss-Lobart nodes are defined as time nodes. The values ​​of surface water depth, soil pressure head, and groundwater level at each time point are used as unknown spectral coefficients. Integrating each term in the governing equations with the spectral basis functions over time units yields a system of algebraic equations with unknown spectral coefficients at all time points as variables, outputting an algebraic system.

[0009] As a further aspect of the present invention: the selection of Legendre polynomials as spectral basis functions within the standard time interval and the definition of Gauss-Lobart nodes as time nodes specifically include: Map the current time unit to a standard interval from negative one to one, and determine the order of the Legendre polynomial to be used based on the rate of change of the local time gradient of the previous time unit. Find the zeros of the derivative of the Legendre polynomial in the standard interval under the corresponding order, and combine these zeros with the two endpoints of the standard interval to form a Gauss-Lobart node sequence. Then, the corresponding node sequence is mapped back to the original time unit, so that the time nodes are distributed with denser ends and sparser middle within the time unit, and the output is a discrete time description with node positions.

[0010] As a further aspect of the present invention: S3 specifically includes: Extract the values ​​of each state variable at all time nodes from the solution of the previous time unit, calculate the absolute value of the difference between adjacent nodes, and take the maximum value as the local time gradient rate of change. Compare the local time gradient rate of change with a preset first threshold and a second threshold: If the local time gradient change rate is greater than the first threshold, the length of the current time unit is shortened to half the length of the previous time unit, and the order of the spectral basis function is increased by two. If the local time gradient change rate is less than the second threshold, the length of the current time unit is extended to twice the length of the previous time unit, and the order of the spectral basis function is reduced by one. If the local time gradient rate of change is between two thresholds, the time unit length remains unchanged, and the order of the spectral basis function is finely adjusted by one order only according to the rising or falling trend of the rate of change. Output the adaptive time unit division and its corresponding order.

[0011] As a further aspect of the present invention: S4 specifically includes: Even-numbered nodes in the adaptive time unit sequence are used as coarse time grid nodes, and the predicted spectral coefficients of each state variable are calculated using the backward Euler scheme on the coarse time grid nodes. Substitute the predicted spectral coefficients into the spectral integral kernel to calculate the weighted integral value of the algebraic system residual along the Legendre polynomial weight function within the coarse time grid interval; The weighted integral value is then used as a correction factor and superimposed onto the predicted spectral coefficients to generate the first corrected solution. Repeat the steps of calculating residuals and corrections until the maximum difference between two adjacent corrected solutions is less than the convergence threshold, and then output the final convergence spectral coefficients.

[0012] As a further aspect of the present invention: the calculation process of the predicted spectral coefficients is as follows: Extract the state variable values ​​of the current coarse time grid starting node from the convergence spectrum coefficients of the previous time unit, and use the state variable values ​​as the initial conditions for the backward Euler iteration; The time derivative term in the control equations is approximated by dividing the difference between the state variables of the starting node and the ending node by the coarse time step, and the spatial discrete term is expressed by the unknown spectral coefficient of the ending node, forming a closed algebraic equation system with the spectral coefficient of the ending node as the unknown. Then, Newton's iteration is used to solve the closed algebraic equation system. In each iteration, the nonlinear terms in the equation are updated using the endpoint node value obtained in the previous iteration until the absolute difference between two adjacent iteration solutions is less than a preset threshold. The predicted spectral coefficients of the endpoint node are then output.

[0013] As a further aspect of the present invention: S5 specifically includes: Substitute the state variables expressed by the convergence spectral coefficients in each time unit into the integral form of the mass conservation equation, calculate the difference between the net change in mass and the boundary flux integral in the entire time unit, and if the absolute value of the difference exceeds the preset tolerance, apply a correction flux at the endpoint of the corresponding time unit to force balance. Extract the lower-order part from the higher-order spectral coefficients within the same time unit, reconstruct the higher-order solution and the lower-order solution respectively, and calculate the maximum difference between the two within the time unit as the local time error; If the local time error is less than the allowable value, the current solution is accepted and the time unit is advanced to the next one; otherwise, the current solution is rejected and the length of the current time unit is reduced to half of the original length. The spectrum solution and error judgment are re-executed until the error meets the requirements and then the process is advanced. The convergence spectral coefficients at the endpoints of each time unit are converted into state variable values, and the flow process lines of the outlet section of the small watershed and the spatiotemporal distribution of water depth and infiltration rate of the entire watershed are output successively according to the endpoints of the time units.

[0014] The beneficial effects of this invention are: (1) By adaptively adjusting the length of the time unit and the order of the spectral basis function, combined with the spectral delay correction iteration, the numerical oscillation of infiltration rate caused by the drastic change in water depth during the rainstorm in the small watershed can be effectively suppressed, and the rapid water rise process can be stably simulated at a time step of minutes.

[0015] (2) Force the integral constraint of mass conservation at each time unit endpoint, and use the difference between high and low order solutions to evaluate the local time error. Automatically reduce the step size and recalculate for unqualified units. This reduces the water balance error of the whole basin, avoids invalid calculations, and improves the simulation efficiency under complex rainfall patterns. Attached Figure Description

[0016] The invention will now be further described with reference to the accompanying drawings.

[0017] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a flowchart of the process for constructing a three-dimensional spatial discrete mesh in this invention. Detailed Implementation

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

[0019] Please see Figure 1 As shown, this invention is a layered, three-dimensional, fully coupled hydrological-hydraulic simulation method for very small watersheds, comprising the following steps: S1 collects digital elevation data, soil vertical layer thickness data, land use type data, and time series rainfall data for the small watershed. Based on the collected data, a three-dimensional spatial discrete grid is constructed, and the total simulation duration is divided into multiple initial time units. The spatial discrete grid and initial time unit division results are output. S2 takes a three-dimensional spatial discrete grid and an initial time unit as input. Within each time unit, it uses a higher-order spectral basis function to approximate the temporal changes of surface water depth, soil pressure head, and groundwater level state variables. It also transforms the set of control equations describing surface water flow, soil water movement, and groundwater movement into a weak integral form based on higher-order spectral basis functions in the spatiotemporal domain. The output is an algebraic system containing unknown spectral coefficients in all time units. S3 takes the algebraic system as input, and dynamically adjusts the order of the higher-order spectral basis function and the length of the time unit based on the local time gradient change rate of the state variables of the previous time unit solution, generating an adaptive time discretization scheme, and outputting the adaptive time unit partitioning and corresponding order. S4 takes the adaptive time discretization scheme and algebraic system as input and uses spectral delay correction iterative solution: first, a low-order time scheme is used to obtain the predicted solution on the coarse time grid, then the residual is calculated using the spectral integral kernel and the predicted solution is corrected one by one until the residual meets the convergence threshold, and the convergence spectral coefficients of the state variables in each time unit are output. S5 takes the convergence spectrum coefficient as input, enforces the mass conservation integral constraint at the endpoint of each time unit, and accepts the current solution or recalculates by reducing the step size according to the local time error control, and advances to the next time unit. It outputs the flow process line of the outlet section of the small watershed and the spatiotemporal distribution of water depth and infiltration rate of the whole watershed successively according to the endpoint of the time unit.

[0020] Please see Figure 2 As shown, in S1, digital elevation data, soil vertical layer thickness data, land use type data, and time-series rainfall data of the small watershed are collected. A three-dimensional spatial discrete grid is constructed based on the collected data, and the total simulation duration is divided into multiple initial time units. The spatial discrete grid and initial time unit division results are output, specifically including: When collecting digital elevation data for the small watershed, UAV lidar measurements were used to acquire ground point clouds, which were then filtered and denoised to generate raster elevation values ​​with a resolution of 1m × 1m. Using these raster elevation values, a surface unstructured triangular mesh was generated using the Delaunay triangulation method. The side length of each triangle was controlled between 2m and 10m; for areas with a slope greater than 25 degrees, the side length was 2m, and for flat areas, it was 10m. Based on the soil vertical layer thickness data, obtained through borehole sampling, each borehole recorded the soil layer boundary every 10cm thickness between the surface and bedrock. The soil layer boundaries from all boreholes within the same watershed were averaged by depth to obtain 3 to 5 vertical layers, each with a thickness of 0-20cm, 20-50cm, 50-100cm, and below 100cm, respectively. Directly below each surface grid node, vertical layer nodes are inserted sequentially according to the proportion of each layer's thickness to the total vertical thickness. This automatically reduces the node spacing to 5 cm in shallow areas with a thickness less than 20 cm, and maintains a 20 cm spacing in deep areas with a thickness greater than 100 cm. Using land use type data derived from visual interpretation of Gaofen-2 satellite imagery, land cover is divided into impervious and permeable surfaces. Impervious surfaces include roofs, roads, and paved areas, while permeable surfaces include farmland, woodland, and grassland. Along the boundary between impervious and permeable surfaces, the triangles crossing this boundary in the unstructured triangular grid are densified by inserting new grid nodes along the boundary and re-subdividing adjacent triangles, ensuring that the grid side length on either side of the boundary does not exceed half the original side length. Based on time-series rainfall data recorded at 1-minute intervals by tipping bucket rain gauges located at the center of the watershed, the total rainfall distribution curve for a complete rainfall event is obtained. Starting from the initial moment recorded by the rain gauge, cumulative rainfall is calculated every minute. The surface runoff paths corresponding to the consecutive periods with the fastest increase in cumulative rainfall are identified as potential runoff paths. Specifically, the difference in cumulative rainfall between adjacent minutes is calculated as the rainfall per minute. Areas corresponding to periods with rainfall per minute greater than 5 mm are connected using flow direction extracted from digital elevation data to form a broken-line path from a high elevation to the watershed outlet. Along this broken-line path, all grid edges intersecting the path in the unstructured triangular mesh on the surface are linearly densified, i.e., new grid nodes are inserted within 0.5 meters on both sides of the path, ensuring that the grid edge length on the densified path does not exceed 1 meter. The aforementioned three-dimensional spatial discrete mesh is composed of the surface triangular mesh and the vertically layered nodes at each node below it. Each vertically layered node carries the three-dimensional coordinates of that node and the thickness information of the soil layer it belongs to.

[0021] Based on time-series rainfall data, which records rainfall values ​​at 1-minute intervals, the difference in rainfall between two adjacent 1-minute records is calculated to obtain the rate of change of rainfall intensity per minute, in millimeters per minute (mm / min). A preset threshold of 2 mm / min is used. When the rate of change of rainfall intensity within a minute exceeds 2 mm / min, the corresponding moment of that minute is recorded as a candidate breakpoint. If multiple consecutive minutes exceed the threshold, only the moment of the first minute is taken as the candidate breakpoint. The watershed runoff path length is extracted using digital elevation data. Specifically, starting from the furthest grid node in the watershed, the path is traced grid by grid along the direction of surface slope descent to the outlet section. The side lengths of all grids during the tracing process are accumulated to obtain the longest runoff path length, denoted as L, in meters. The surface roughness is determined using the Manning coefficient, with the average value obtained from a table based on land use type data; for example, 0.15 is used for paddy fields and 0.03 for dry land. The water propagation time T is calculated as follows: divide the longest confluence path length L by the average flow velocity. The average flow velocity is calculated by multiplying the reciprocal of the Manning coefficient by the square of the water depth (3 / 4) and then by the square of the slope (1 / 2). The water depth is taken as an empirical value of 0.1 meters, and the slope is taken as the average slope of the watershed. Output the value of T in seconds. For any two adjacent candidate breakpoints, calculate the time interval between them. If the time interval is greater than the water propagation time T, then the time interval is used as an initial time unit; if the time interval is less than or equal to T, then the two candidate breakpoints are merged, and no separate time unit is defined. For periods of stable rainfall intensity, i.e., periods where the rate of change in rainfall intensity is less than 0.5 mm / min for more than 10 consecutive minutes, supplementary divisions are made using an equal-length method. The length of the equal-length division is half of the water propagation time T, but does not exceed 60 seconds. After the above division, each initial time unit contains at most one point of sudden change in rainfall intensity (i.e., the moment when the rate of change in rainfall intensity exceeds 2 mm / min). Output the start and end times of each initial time unit, as well as the corresponding surface grid and vertical layer node numbers.

[0022] In S2, a three-dimensional spatial discrete grid and initial time units are used as input. Within each time unit, higher-order spectral basis functions are used to approximate the temporal changes of surface water depth, soil pressure head, and groundwater level state variables. The governing equations describing surface water flow, soil water movement, and groundwater movement are uniformly transformed into weak integral forms based on higher-order spectral basis functions in the spatiotemporal domain. The output is an algebraic system containing unknown spectral coefficients in all time units, specifically including: The three-dimensional spatial discrete grid and the start and end times of each initial time unit are used as inputs. For each initial time unit, its time interval is first linearly transformed to a standard interval from -1 to 1: Let the start time of the initial time unit be T0 and the end time be T1. Then the relationship between the standard time variable τ and the original time t is: τ equals twice the amount of t minus T0 divided by the amount of T1 minus T0 minus 1, thus τ = -1 when t = T0 and τ = 1 when t = T1. Within the standard interval [-1, 1], Legendre polynomials are selected as the spectral basis functions, and the order of the Legendre polynomial is denoted as . Define Gaussian-Lobart nodes as time nodes; this sequence of nodes includes the two endpoints of the interval and the zeros of the Legendre polynomial derivative in the open interval. For a given order... Legendre polynomial derivative yes A polynomial of degree 1, which has the following expression in (-1, 1) A number of distinct zeros are used to form a combination of these zeros and two endpoints. Nodes Surface water depth, soil pressure head, and groundwater level were measured at each time point. The values ​​on are taken as unknown spectral coefficients, and are denoted as follows: , and For each state variable, its change within a time unit is obtained by summing the spectral coefficients at all nodes by multiplying them by the corresponding Lagrange interpolation basis function, which is defined at each node. The value is 1 at one node and 0 at other nodes. Each term in the governing equations describing surface water flow, soil water movement, and groundwater movement is multiplied by each spectral basis function within a time unit, and then Gaussian-Lobart numerical integration is performed over a standard interval, utilizing the same set of nodes. and the corresponding weight coefficients The integral is approximated as the sum of the products of the function values ​​at each node and the weight coefficients. (Weight coefficients) The calculation method is as follows: for endpoint nodes and The weighting coefficient equals Multiply Divide by the reciprocal of the square of the value of the Legendre polynomial at the endpoints; for internal nodes The weighting coefficient is equal to 1 minus the weighting coefficient inside the parentheses. The derivative of the square multiplied by the Legendre polynomial is in The reciprocal of the product of the squares of the values. After the above integration operation, the governing equations within each initial time unit are transformed into an algebraic equation system with the unknown spectral coefficients at all time points as variables. The algebraic equation systems of all initial time units are combined in chronological order to form a large algebraic system, and the coefficient matrix and right-hand side vector of this algebraic system are output.

[0023] After mapping the current initial time unit to the standard interval [-1, 1], the order of the Legendre polynomial to be used is determined based on the spectral energy distribution of the state variables of the solution from the previous time unit. Within the previous time unit, the change in surface water depth over time was expanded using Legendre polynomials to obtain the sequence of spectral coefficients of each order. Let the th... The spectral coefficients at the order are in the first order. The contribution value at each time point ,in From 1 to , Let be the total number of Gaussian-Lobart nodes in the previous time unit. Then, the energy of the k-th order spectral coefficients... Calculate using the following formula: ; Calculate the sum of the energies for all orders to obtain the total energy. Take the sum of the energies of the two highest orders. ,like Divide by If the ratio is greater than 0.1, then the order of the current time unit will be adjusted. Set it to the order of the previous unit plus 2 (but not exceeding 20); if the ratio is less than 0.01, then Set the order to the previous unit order minus 1 (but not lower than 2); otherwise, keep the order unchanged. Determine the order. Then, the Legendre polynomials need to be solved. derivative exist The zero point within the interval. Solving using Newton's iteration method: uniformly select within the interval... Initial guess value For each guessed value, iteratively update using the following formula: ; in For the second derivative of the Legendre polynomial. Repeat the iteration until... Less than To obtain zero point to Connect the two endpoints and These zeros together form the Gauss-Lobart node sequence. This node sequence is then transformed back to the original time unit using an inverse mapping, the formula for which is: Because the Gaussian-Lobart nodes are dense at the endpoints and sparse in the middle, the transformed time nodes also exhibit a dense-at-the-ends and sparse-in-the-middle distribution within the original time unit. Output the node position sequence for each initial time unit and its corresponding Legendre polynomial order. This serves as a time-discrete description of the time unit, for use in subsequent steps.

[0024] In S3, the algebraic system is taken as input. Based on the local time gradient rate of change of the state variables of the previous time unit solution, the order of the higher-order spectral basis functions and the time unit length of the current time unit are dynamically adjusted to generate an adaptive time discretization scheme. The adaptive time unit partitioning and corresponding order are output, specifically including: Using an algebraic system as input, the values ​​of three state variables—surface water depth, soil pressure head, and groundwater level—at all Gauss-Lobart time nodes are extracted from the solution of the previous time unit. For each state variable, the absolute value of the difference between the values ​​at adjacent nodes is calculated sequentially, and then divided by the time interval between adjacent nodes to obtain the average rate of change for each sub-interval. The maximum value among all the average rates of change in sub-intervals is taken as the local temporal gradient rate of change for that state variable. After calculating the local temporal gradient rates of change for each of the three state variables, the maximum value is selected as the local temporal gradient rate of change for the current time unit, denoted as the gradient value.

[0025] The first threshold is preset to 0.8 mm / min, and the second threshold is 0.08 mm / min. The calculated gradient value is compared with the first and second thresholds. If the gradient value is greater than the first threshold, it indicates that the state variable is changing drastically. In this case, the length of the current time unit is set to half the length of the previous time unit, and the order of the spectral basis function is increased by two orders compared to the previous time unit, but the order shall not exceed twenty. If the gradient value is less than the second threshold, it indicates that the state variable is changing gradually. In this case, the length of the current time unit is set to twice the length of the previous time unit, and the order of the spectral basis function is decreased by one order, but the order shall not be lower than two orders. If the gradient value is between the second threshold and the first threshold, the length of the current time unit is kept the same as the previous time unit. At this time, it is necessary to further determine the trend of gradient value: compare the currently calculated gradient value with the gradient value calculated in the previous time unit. If the current gradient value is greater than the gradient value in the previous time unit, it indicates that the trend is upward, so the order of the spectral basis function is increased by one; if the current gradient value is less than the gradient value in the previous time unit, it indicates that the trend is downward, so the order of the spectral basis function is decreased by one; if the two are equal, the order remains unchanged.

[0026] When adjusting the time unit length, it is also necessary to check whether the length exceeds a reasonable range: the shortest time unit length is set to 0.01 seconds, and the longest time unit length is set to 60 seconds. If the length after shortening is less than 0.01 seconds, then take 0.01 seconds; if the length after extending is greater than 60 seconds, then take 60 seconds. After adjustment, output the start time, end time, and length of the current time unit, as well as the order of the spectral basis function corresponding to that time unit, as part of the adaptive time discretization scheme. Repeat the above process, processing time units one by one, until the entire simulation duration is covered, and finally output the division results of all time units and the order corresponding to each unit.

[0027] In S4, the adaptive time discretization scheme and the algebraic system are taken as input, and the spectral delay correction iterative solution is used: first, a low-order time scheme is used to obtain the predicted solution on a coarse time grid; then, the residual is calculated using a spectral integral kernel, and the predicted solution is corrected successively until the residual meets the convergence threshold. The convergence spectral coefficients of the state variables in each time unit are output, specifically including: The adaptive time discretization scheme and the algebraic system are used as input. First, from the adapted time unit sequence, a node is selected for every other time unit; that is, the even-numbered time in the order of the start and end times of all time units is taken as the coarse time grid node. The interval between two adjacent coarse time grid nodes is called the coarse time grid interval, and its length is denoted as the coarse time step. For each coarse time grid interval, the predicted spectral coefficients of each state variable within the interval are calculated using the backward Euler scheme.

[0028] The specific calculation process for the predicted spectral coefficients is as follows: From the spectral coefficients that have converged in the previous time unit, extract the surface water depth, soil pressure head, and groundwater level values ​​at the starting node of the current coarse time grid. Use these values ​​as the initial conditions for the backward Euler iteration. In the backward Euler scheme, the time derivative term in the governing equations is approximated as: the state variable value at the endpoint minus the state variable value at the starting node, divided by the coarse time step. The spatial discrete terms in the governing equations are all expressed using the unknown spectral coefficients at the endpoint, including the water depth gradient in the surface water flow term, the second derivative of the pressure head in the soil water movement term, and the Laplace operator of the groundwater movement term. This constitutes a closed algebraic equation system with the unknown spectral coefficients of the three state variables at the endpoint as the solution object. The closed algebraic equation system is solved using Newton's iterative method: First, an initial guess value is assigned to the unknown spectral coefficients of the endpoint node, which is taken as the state variable value at the starting node. Then, the residual vector of the algebraic equation system under the current guess value and the partial derivative matrix of the residual with respect to each unknown spectral coefficient are calculated. A correction is obtained by solving the linear equation system, and the guess value is added to the correction to obtain the updated guess value. This process is repeated until the maximum absolute value of the difference between the spectral coefficients of the endpoint node obtained from two adjacent iterations is less than a preset threshold, which is set to 10 to the power of -8. After the iteration converges, the predicted spectral coefficients of the endpoint node of the coarse time grid interval are output. The above process is performed sequentially for all coarse time grid intervals to obtain the predicted spectral coefficients at all coarse time grid nodes.

[0029] After obtaining the predicted spectral coefficients, they are substituted into the spectral integral kernel. The spectral integral kernel is defined as follows: within a coarse time grid interval, the residual function of the algebraic system is multiplied by the Legendre polynomial weight function, and the result is integrated over the entire interval. Specifically, the coarse time grid interval is mapped to a standard interval from -1 to 1. Using the Gauss-Lobart nodes and corresponding weight coefficients within this interval, the residual function is calculated at each node, multiplied by the corresponding node's weight coefficient, and summed to obtain the weighted integral value. This weighted integral value is used as a correction and back-stacked onto the predicted spectral coefficients; that is, the predicted spectral coefficients are subtracted from the correction value to generate the first corrected solution. After one correction, the corrected solution is used as the new predicted value. The weighted integral value obtained from the spectral integral kernel is calculated again, and back-stacked again to obtain the second corrected solution. The above steps of calculating the residual weighted integral and back-stacking correction are repeated. After each iteration, the current corrected solution is compared with the previous corrected solution, and the absolute value of the maximum difference between the two at all nodes is taken. When the maximum difference is less than the convergence threshold (set to the power of 10 - 10), the iteration stops, and the current corrected solution is output as the final convergence spectral coefficients. If convergence is not achieved after more than 20 iterations, the length of the coarse time grid interval is halved, and prediction and correction are performed again. Finally, the convergence spectral coefficients of each state variable in each time unit at all Gauss-Lobart nodes are output.

[0030] In S5, the convergence spectral coefficients are used as input. Mass conservation integral constraints are enforced at the endpoints of each time unit. Based on local time error control, the current solution is accepted or the step size is reduced for recalculation, proceeding to the next time unit. The flow process curves at the outlet section of the ultra-small watershed and the spatiotemporal distribution of water depth and infiltration rate across the entire watershed are output sequentially according to the endpoints of the time units. Specifically, this includes: The convergence spectral coefficients of the state variables within each time unit are used as input. For each time unit, the functions of surface water depth, soil pressure head, and groundwater level as expressed by the convergence spectral coefficients over time are first substituted into the integral form of the mass conservation equation. The integral form of the mass conservation equation is: within a time unit, the net change in water mass within the watershed equals the amount of water flowing into the boundary of the unit minus the amount of water flowing out of the boundary, plus the rainfall input. In specific calculations, the cumulative amount of mass change over time within the entire time unit is obtained using the spectral coefficient values ​​at all Gauss-Lobart nodes within the time unit and a numerical integration method. Simultaneously, based on the three-dimensional spatial discrete grid output from step one, the total water storage of the entire watershed at the start and end times of the time unit is calculated, and the difference between the two is the net change in mass. The boundary flux integral is calculated based on surface boundary conditions (such as no-flow boundary or given head boundary) and subsurface boundary conditions. The calculated net change in mass is subtracted from the boundary flux integral to obtain the difference, and the absolute value of this difference is taken. The preset tolerance value is one-thousandth of the total water storage of the entire basin. If the absolute value of this difference exceeds the preset tolerance, it indicates that the mass conservation within that time unit has not been satisfied. In this case, a correction flux is applied at the endpoint of that time unit. The correction flux is calculated by dividing the difference between the net change in mass and the integral of the boundary flux by the length of the time unit to obtain the correction flow rate per unit time. This correction flow rate is then uniformly superimposed on all surface grid nodes at the endpoint of the time unit to force the balance of mass error.

[0031] After completing the forced correction based on mass conservation, local time error is assessed. Low-order components are extracted from the higher-order convergent spectral coefficients within the same time unit. Specifically, if the order of the spectral basis function for that time unit is P, the first Q order spectral coefficients are used to construct the low-order solution, where Q equals P minus 2, but Q must not be lower than 2. The state variable variation curve within that time unit is reconstructed using the low-order spectral coefficients to obtain the low-order solution; simultaneously, the higher-order solution is reconstructed using all P order spectral coefficients. Twenty sampling points are uniformly selected within the time unit, and the absolute value of the difference between the higher-order and lower-order solutions at each sampling point is calculated. The maximum difference among all sampling points is taken as the local time error. The preset allowable values ​​for local time error are 0.001 mm (for water depth) and 0.0001 mm water column (for pressure head). The calculated local time error is compared with the allowable value: if the local time error is less than the allowable value, the solution for the current time unit is accepted, the time unit is marked as converged, and the process continues to the next time unit, continuing steps three through five; if the local time error is greater than or equal to the allowable value, the solution for the current time unit is rejected, the calculation result for that unit is not saved, and the length of the current time unit is reduced to half of its original length. After reduction, starting from step three again, the dynamic adjustment, spectral delay correction solution, and mass conservation and error judgment of this step are performed again for this shortened time unit until the local time error meets the requirements. Only when the solution of a time unit passes both the mass conservation constraint and the local time error control is its spectral coefficients retained as a valid result and advanced to the next time unit.

[0032] After completing the calculations for all time units, the convergence spectral coefficients at the endpoints of each time unit are concatenated in chronological order. At each endpoint, the surface water depth and infiltration rate at all three-dimensional discrete grid nodes across the entire watershed are calculated based on the spectral coefficients. The infiltration rate is calculated using Darcy's law to determine the vertical flux at the surface using the spectral coefficients of soil pressure head. Simultaneously, the flow rate at the watershed outlet is calculated using Manning's formula based on the surface water depth and topographic slope. The calculation results are output sequentially at the endpoints of each time unit, yielding the flow process curve at the outlet of the small watershed, as well as the water depth and infiltration rate distributions at different times for each grid node across the entire watershed. The output time interval is consistent with the endpoint times of each time unit, forming a complete spatiotemporal distribution result.

[0033] The working principle of this invention is as follows: First, digital elevation data, soil vertical layer thickness, land use type, and time-series rainfall data are collected for a small watershed. Based on these data, a three-dimensional spatial discrete grid is constructed, including a surface unstructured triangular grid and vertical layer nodes. The total simulation duration is divided into initial time units, each containing at most one rainfall intensity abrupt change point, according to the rainfall intensity variation rate and confluence propagation time. Then, using the spatial discrete grid and initial time units as input, each time unit is mapped to a standard interval. Legendre polynomials are selected as the spectral basis functions, and Gauss-Lobart nodes are defined. The values ​​of surface water depth, soil pressure head, and groundwater level at each node are used as unknown spectral coefficients. The governing equations are transformed into an algebraic system with all node spectral coefficients as variables using the Galerkin weighted residual method. Then, based on… The local time gradient change rate of the solution in the previous time unit is used to dynamically adjust the order of the spectral basis function and the time unit length of the current time unit to generate an adaptive time discretization scheme. Then, spectral delay correction is used for iterative solution. First, the predicted spectral coefficients are obtained on the coarse time grid using the backward Euler scheme. Then, the weighted integral value of the residual is calculated using the spectral integral kernel and repeatedly back-stacked for correction until the residual converges. The convergent spectral coefficients of the state variables in each time unit are output. Finally, the mass conservation integral constraint is enforced at the endpoint of each time unit, and the local time error is evaluated by the difference between the higher-order solution and the lower-order solution. If the error exceeds the limit, the step size is reduced and recalculated. Otherwise, the current solution is accepted and the process is advanced to the next time unit. Finally, the flow process line of the outlet section of the small watershed and the spatiotemporal distribution of water depth and infiltration rate of the whole watershed are output successively according to the endpoint of the time unit.

[0034] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A layered, three-dimensional, fully coupled hydrological-hydraulic simulation method for very small watersheds, characterized in that, Includes the following steps: S1 collects digital elevation data, soil vertical layer thickness data, land use type data, and time series rainfall data for the small watershed. Based on the collected data, a three-dimensional spatial discrete grid is constructed, and the total simulation duration is divided into multiple initial time units. The spatial discrete grid and initial time unit division results are output. S2 takes a three-dimensional spatial discrete grid and an initial time unit as input. Within each time unit, it uses a higher-order spectral basis function to approximate the temporal changes of surface water depth, soil pressure head, and groundwater level state variables. It also transforms the set of control equations describing surface water flow, soil water movement, and groundwater movement into a weak integral form based on higher-order spectral basis functions in the spatiotemporal domain. The output is an algebraic system containing unknown spectral coefficients in all time units. S3 takes the algebraic system as input, and dynamically adjusts the order of the higher-order spectral basis function and the length of the time unit based on the local time gradient change rate of the state variables of the previous time unit solution, generating an adaptive time discretization scheme, and outputting the adaptive time unit partitioning and corresponding order. S4 takes the adaptive time discretization scheme and algebraic system as input and uses spectral delay correction iterative solution: first, a low-order time scheme is used to obtain the predicted solution on the coarse time grid, then the residual is calculated using the spectral integral kernel and the predicted solution is corrected one by one until the residual meets the convergence threshold, and the convergence spectral coefficients of the state variables in each time unit are output. S5 takes the convergence spectrum coefficient as input, enforces the mass conservation integral constraint at the endpoint of each time unit, and accepts the current solution or recalculates by reducing the step size according to the local time error control, and advances to the next time unit. It outputs the flow process line of the outlet section of the small watershed and the spatiotemporal distribution of water depth and infiltration rate of the whole watershed successively according to the endpoint of the time unit.

2. The method for stratified three-dimensional hydrological-hydraulic fully coupled simulation of a very small watershed according to claim 1, characterized in that, The construction of the three-dimensional spatial discrete mesh specifically includes: A surface unstructured triangular grid is generated using digital elevation data. Then, based on the soil vertical layer thickness data, vertical layer nodes are inserted below each surface grid node according to the layer thickness ratio, so that the thin layer area is automatically densified. Land use type data is used to identify permeable and impermeable boundaries, and the surface grid is densified along the boundary lines. Based on the cumulative rainfall distribution of time series rainfall data, potential runoff paths are extracted, and the surface grid is linearly densified along the paths to obtain a three-dimensional spatial discrete grid.

3. The method for stratified three-dimensional hydrological-hydraulic fully coupled simulation of a very small watershed according to claim 1, characterized in that, The process of dividing the total simulation duration into multiple initial time units specifically includes: The rate of change of rainfall intensity between adjacent moments is calculated based on time series rainfall data, and moments when the rate of change of rainfall intensity exceeds a preset threshold are set as candidate breakpoints. The length of the watershed confluence path is extracted using digital elevation data, and the water propagation time is estimated by combining it with surface roughness. The time interval between adjacent candidate breakpoints that is longer than the propagation time is taken as an initial time unit. During periods of stable rainfall intensity, the time units are divided into equal-length segments to ensure that each initial time unit contains at most one point of sudden change in rainfall intensity, and the division results are output.

4. The method for stratified three-dimensional hydrological-hydraulic fully coupled simulation of a very small watershed according to claim 1, characterized in that, The output process of the algebraic system specifically includes: Each time unit is mapped to a standard time interval, Legendre polynomials are selected as spectral basis functions within the standard time interval, and Gauss-Lobart nodes are defined as time nodes. The values ​​of surface water depth, soil pressure head, and groundwater level at each time point are used as unknown spectral coefficients. Integrating each term in the governing equations with the spectral basis functions over time units yields a system of algebraic equations with unknown spectral coefficients at all time points as variables, outputting an algebraic system.

5. The layered three-dimensional hydrological-hydraulic fully coupled simulation method for very small watersheds according to claim 4, characterized in that, The selection of Legendre polynomials as spectral basis functions within the standard time interval and the definition of Gauss-Lobart nodes as time nodes specifically include: Map the current time unit to a standard interval from negative one to one, and determine the order of the Legendre polynomial to be used based on the rate of change of the local time gradient of the previous time unit. Find the zeros of the derivative of the Legendre polynomial in the standard interval under the corresponding order, and combine these zeros with the two endpoints of the standard interval to form a Gauss-Lobart node sequence. Then, the corresponding node sequence is mapped back to the original time unit, so that the time nodes are distributed with denser ends and sparser middle within the time unit, and the output is a discrete time description with node positions.

6. The method for stratified three-dimensional hydrological-hydraulic fully coupled simulation of a very small watershed according to claim 1, characterized in that, S3 specifically includes: Extract the values ​​of each state variable at all time nodes from the solution of the previous time unit, calculate the absolute value of the difference between adjacent nodes, and take the maximum value as the local time gradient rate of change. Compare the local time gradient rate of change with a preset first threshold and a second threshold: If the local time gradient change rate is greater than the first threshold, the length of the current time unit is shortened to half the length of the previous time unit, and the order of the spectral basis function is increased by two. If the local time gradient change rate is less than the second threshold, the length of the current time unit is extended to twice the length of the previous time unit, and the order of the spectral basis function is reduced by one. If the local time gradient rate of change is between two thresholds, the time unit length remains unchanged, and the order of the spectral basis function is finely adjusted by one order only according to the rising or falling trend of the rate of change. Output the adaptive time unit division and its corresponding order.

7. The method for stratified three-dimensional hydrological-hydraulic fully coupled simulation of a very small watershed according to claim 1, characterized in that, S4 specifically includes: Even-numbered nodes in the adaptive time unit sequence are used as coarse time grid nodes, and the predicted spectral coefficients of each state variable are calculated using the backward Euler scheme on the coarse time grid nodes. Substitute the predicted spectral coefficients into the spectral integral kernel to calculate the weighted integral value of the algebraic system residual along the Legendre polynomial weight function within the coarse time grid interval; The weighted integral value is then used as a correction factor and superimposed onto the predicted spectral coefficients to generate the first corrected solution. Repeat the steps of calculating residuals and corrections until the maximum difference between two adjacent corrected solutions is less than the convergence threshold, and then output the final convergence spectral coefficients.

8. The method for stratified three-dimensional hydrological-hydraulic fully coupled simulation of a very small watershed according to claim 7, characterized in that, The calculation process for the predicted spectral coefficients is as follows: Extract the state variable values ​​of the current coarse time grid starting node from the convergence spectrum coefficients of the previous time unit, and use the state variable values ​​as the initial conditions for the backward Euler iteration; The time derivative term in the control equations is approximated by dividing the difference between the state variables of the starting node and the ending node by the coarse time step, and the spatial discrete term is expressed by the unknown spectral coefficient of the ending node, forming a closed algebraic equation system with the spectral coefficient of the ending node as the unknown. Then, Newton's iteration is used to solve the closed algebraic equation system. In each iteration, the nonlinear terms in the equation are updated using the endpoint node value obtained in the previous iteration until the absolute difference between two adjacent iteration solutions is less than a preset threshold. The predicted spectral coefficients of the endpoint node are then output.

9. The method for stratified three-dimensional hydrological-hydraulic fully coupled simulation of a very small watershed according to claim 1, characterized in that, S5 specifically includes: Substitute the state variables expressed by the convergence spectral coefficients in each time unit into the integral form of the mass conservation equation, calculate the difference between the net change in mass and the boundary flux integral in the entire time unit, and if the absolute value of the difference exceeds the preset tolerance, apply a correction flux at the endpoint of the corresponding time unit to force balance. Extract the lower-order part from the higher-order spectral coefficients within the same time unit, reconstruct the higher-order solution and the lower-order solution respectively, and calculate the maximum difference between the two within the time unit as the local time error; If the local time error is less than the allowable value, the current solution is accepted and the time unit is advanced to the next one; otherwise, the current solution is rejected and the length of the current time unit is reduced to half of the original length. The spectrum solution and error judgment are re-executed until the error meets the requirements and then the process is advanced. The convergence spectral coefficients at the endpoints of each time unit are converted into state variable values, and the flow process lines of the outlet section of the small watershed and the spatiotemporal distribution of water depth and infiltration rate of the entire watershed are output successively according to the endpoints of the time units.

Images