Method for simulating rainfall runoff based on nonlinear confluence mechanism of underground base flow

Abstract

The application discloses a rainfall runoff simulation method based on a nonlinear confluence mechanism of underground base flow, and comprises the following steps: obtaining precipitation and evapotranspiration data of a basin, generating a basin average input sequence after spatial weighting and quality control; calculating total yield flow based on the sequence, and combining a free water storage structure to divide the total yield flow into surface runoff, interflow and underground runoff; converting the surface runoff and the interflow into flow of the same time step and performing linear regulation and confluence calculation; constructing an underground water storage balance model, updating the storage based on the underground runoff, and adopting a state self-adaptive mechanism to smoothly adjust a nonlinear index according to a storage saturation state, and outputting nonlinear base flow; and superimposing the three water sources to obtain a total runoff process through river network calculation. The application overcomes the limitation of a traditional single fixed attenuation mechanism, and can improve the comprehensive simulation precision of the model in the flood peak period and the runoff recession stage in the dry season.

Description

Technical Field

[0001] This invention relates to the field of hydrological and water resources simulation and forecasting technology, and in particular to a rainfall-runoff simulation method based on the nonlinear confluence mechanism of underground baseflow. Background Technology

[0002] Rainfall-runoff simulation is a core foundation of hydrological science. Its technical significance lies in quantitatively characterizing the transport and distribution patterns of atmospheric precipitation transformed into surface and groundwater flows through mathematical and physical models. High-precision runoff simulation can reveal the multidimensional dynamic mechanisms of the watershed water cycle, providing crucial data support for flood peak projection, flood resource utilization, and baseflow maintenance during the dry season. It is a fundamental engineering tool for watershed water resource allocation and water environment capacity assessment.

[0003] Current hydrological models primarily rely on empirical conceptual structures to simulate runoff processes. In handling groundwater runoff, existing conventional models generally employ a linear reservoir theory with constant coefficients; that is, a fixed recession coefficient is set, causing the calculated baseflow to decay exponentially over time. After parameter calibration, these methods perform calculations using an invariant physical hysteresis mechanism, regardless of whether the aquifer is in a high- or low-storage state, lacking dynamic response to the deep groundwater migration state and tracking of storage balance.

[0004] However, due to the highly complex pore structure and permeability characteristics of natural aquifers, groundwater runoff exhibits a strong state dependence throughout the entire hydrological cycle. Existing methods, employing fixed-parameter linear assumptions, struggle to reconcile the accelerated discharge mechanisms during aquifer saturation with the slow-release mechanisms during extreme drought. This leads to simulation systems frequently underestimating flood peaks or experiencing severe imbalances in water volume during the receding period when handling different hydrological years. Therefore, there is an urgent need to develop a runoff calculation method that can adapt to different hydrological physical states and achieve high-precision simulations of both high and low flow rates, thereby improving the reliability and robustness of hydrological forecasts under extreme climatic conditions. Summary of the Invention

[0005] The purpose of this invention is to provide a rainfall-runoff simulation method based on the nonlinear confluence mechanism of underground baseflow, in order to solve the above-mentioned problems in the prior art.

[0006] To achieve the above objectives, the present invention provides the following technical solution, including: Acquire watershed precipitation and evapotranspiration data, perform spatial weighting and quality control processing on them, and obtain the watershed average input sequence; The total runoff is calculated based on the average input sequence of the watershed, and the total runoff is divided into surface runoff, interflow and groundwater runoff based on the free water storage structure. Based on the dimensional transformation of surface runoff and interflow, and using linear regulation and storage mechanisms for runoff calculation, the surface runoff flow and interflow flow were obtained. A groundwater storage balance model is constructed, and the groundwater storage is dynamically updated based on the groundwater runoff. A state adaptive mechanism is adopted to smoothly adjust the nonlinear exponent of the baseflow calculation according to the saturation state of the groundwater storage, and the nonlinear groundwater baseflow is calculated. The groundwater storage at the previous moment is recursively updated based on the difference between the groundwater runoff and the nonlinear groundwater baseflow to obtain the groundwater storage at the current moment. The total runoff process of the watershed is obtained by superimposing the surface runoff, the interflow runoff, and the nonlinear groundwater baseflow, and then processing them through river network lag propagation calculation.

[0007] Among them, a dimensional conversion was performed based on surface runoff and interflow, and a linear regulation and storage mechanism was used for runoff calculation to obtain surface runoff and interflow runoff, including: Surface runoff and interflow, measured in units of runoff depth, are uniformly converted into flow rates within the corresponding time step. By introducing a surface flow recession coefficient, a linear reservoir model is used to iteratively calculate the surface runoff converted into flow rate, and the surface runoff is output. By introducing the interflow recession coefficient, a linear reservoir model is used to iteratively calculate the interflow converted into flow rate, and the interflow confluence flow rate is output.

[0008] The computational parameters required for method execution are automatically calibrated through the following pre-construction steps: Acquire historical measured flow data of the watershed and corresponding historical input meteorological data; Based on historical meteorological data, the preceding rainfall-runoff simulation process is executed to calculate historical simulated flow data; Based on historical measured flow data of the basin and combined with historical simulated flow data, a multi-objective evaluation function is constructed. The multi-objective evaluation function includes the Nash efficiency coefficient index for evaluating the process fitting accuracy and the flood volume relative error index for evaluating the water balance state. A pre-defined global composite evolutionary optimization algorithm is adopted to automatically iterate and optimize the parameter space involving the flow generation mechanism and the nonlinear confluence mechanism, with the goal of satisfying the comprehensive optimal convergence condition of the multi-objective evaluation function, and outputs a pre-calibrated model parameter set.

[0009] Before iteratively summing the historical groundwater storage from the previous moment based on the water volume difference, the process also includes the step of determining the initial groundwater storage: Acquire historical groundwater depth data and measured soil moisture content data for the target watershed; Based on historical groundwater depth data and measured soil moisture content data, the initial groundwater storage is determined by water volume conversion, which serves as the starting point for the historical groundwater storage at the first time step of the groundwater storage balance model.

[0010] After obtaining the total runoff process of the watershed, it also includes: The total runoff process of the basin is compared and analyzed with the preset flood control warning threshold of the basin; When the peak flow in the total runoff process of the basin is determined to exceed the basin flood control warning threshold, a flood warning signal containing the peak arrival time and the estimated flow is automatically generated, or the flood control scheduling strategy of the corresponding water conservancy hub in the basin is triggered.

[0011] Beneficial effects: By constructing an independent differential groundwater storage balance model, the macroscopic watershed water volume status is quantified into a standard continuous relative saturation control variable, enabling physical tracking of the evolution of the aquifer's abundance and scarcity, and overcoming the technical deficiency of traditional runoff simulation methods in lacking tracking of deep water transport status.

[0012] Building upon this, to resolve the contradiction between the fixed-parameter linear assumption and the accelerated discharge mechanism during flood season and the slow-release mechanism during dry season, a state-adaptive nonlinear baseflow calculation mechanism with bidirectional asymptotic characteristics is introduced. This allows the exponent controlling the discharge rate to be smoothly and dynamically adjusted according to the current relative saturation of the aquifer, correcting the underestimation of peak flow during flood season and reducing the calculation error of receding water volume during dry season, thus achieving high-precision dynamic fitting of high and low flow rates throughout the entire hydrological cycle.

[0013] Furthermore, to address the challenges of numerical divergence and model parameter calibration inherent in nonlinear hydrological systems, a pre-defined reference constant is introduced to perform dimensionless base calculations to eliminate dimensional drift. Additionally, a three-dimensional physical extreme boundary constraint encompassing storage, exponent, and baseflow is applied at the end of the nonlinear calculation. Simultaneously, a multi-objective evaluation function is constructed by combining process fit and water balance indices for offline global optimization. In summary, this enhances the underlying numerical stability of the nonlinear calculations, enabling the model to maintain good forecast robustness and spatial generalization ability even under extreme climate conditions. Attached Figure Description

[0014] Figure 1 This is a schematic diagram of the overall process of the rainfall-runoff simulation method based on the nonlinear confluence mechanism of underground baseflow of the present invention.

[0015] Figure 2 This is a schematic diagram of the method for calculating and classifying total output flow and its types according to the present invention.

[0016] Figure 3 This is a schematic diagram of the method for constructing a groundwater storage balance model and updating groundwater storage according to the present invention.

[0017] Figure 4 This is a schematic diagram of the calculation method for nonlinear underground base flow of the present invention. Detailed Implementation

[0018] Example 1: like Figure 1 As shown in the figure, this embodiment provides a rainfall-runoff simulation method based on the nonlinear confluence mechanism of underground baseflow, which mainly includes the following steps: Step 101: Obtain watershed precipitation and evapotranspiration data, perform spatial weighting and quality control processing on them, and obtain the watershed average input sequence.

[0019] Precipitation and evapotranspiration, as the main driving factors of the hydrological cycle, often exhibit significant spatial non-uniformity in their spatial distribution.

[0020] Specifically, watershed precipitation and evapotranspiration data can be acquired in real time through remote sensing stations deployed within the target watershed, or directly retrieved from the meteorological department's historical data database.

[0021] In some application scenarios, since single-point meteorological data cannot directly represent the macroscopic physical state of the entire basin, spatial interpolation technology is needed to transform it into a surface average, i.e., spatial weighting.

[0022] Furthermore, due to potential sensor malfunctions or communication interruptions, the original meteorological sequence often contains outliers and data gaps. By performing quality control processing, outliers exceeding statistical boundary extremes can be identified and removed. Gaps can be filled by utilizing the spatiotemporal correlations between adjacent time periods or stations, thus constructing a watershed-averaged input sequence with spatiotemporal continuity and numerical stability. This provides noise-free boundary conditions for subsequent hydrophysical calculations.

[0023] As an alternative implementation method, evapotranspiration data can be indirectly calculated by collecting measured meteorological elements, such as temperature, wind speed, and radiation, and combining them with preset empirical physical formulas (such as the Penman-Monteith formula or the Hargreaves formula), or by directly loading measured observation data from the evaporation pan.

[0024] Step 102: Calculate the total runoff based on the average input sequence of the watershed to obtain the total runoff, and divide the total runoff into surface runoff, interflow and groundwater runoff in combination with the free water storage structure.

[0025] Runoff calculation is a primary response process that converts meteorological driving forces into displacements of surface and groundwater bodies.

[0026] After obtaining the continuous average input sequence of the watershed, the system uses the preset soil moisture content curve and evaporation conversion factor to perform time-series iteration. Based on whether the current soil moisture has reached saturation, the total runoff generated within the current time step is calculated.

[0027] The hydrological model no longer treats the flow as a single entity, but rather breaks it down based on the stratified water conductivity characteristics of the soil profile.

[0028] Specifically, based on the free water storage structure and its dynamic distribution ratio mechanism, the total yield is forcibly divided into surface runoff with extremely fast response, interflow with moderate lateral migration speed, and groundwater runoff that slowly seeps into the aquifer.

[0029] The three-source water distribution model can accurately depict the physical distribution path of rainwater on surfaces with different porosities, avoiding the confusion between the rapid flow at the beginning of a flood and the slow flow during the receding period.

[0030] Step 103: Based on the unified dimensional conversion of surface runoff and interflow, and using the linear regulation and storage mechanism to perform runoff calculations, the surface runoff flow rate and interflow flow rate are obtained.

[0031] In the hydrophysical field, the surface runoff and interflow obtained in the previous step are usually represented by runoff depth, while the observation benchmark for the watershed outlet control section is the flow rate.

[0032] To prevent dimensional conflicts when multiple water sources overlap in subsequent calculations, a unified dimensional conversion needs to be performed before confluence calculations. Specifically, by combining the total watershed area parameter and the current time step, surface runoff and interflow are uniformly converted into volumetric flow rate.

[0033] For these two types of water sources that migrate in the shallow layer, different receding constants are introduced, and independent iterative calculations are performed using a linear regulation and storage mechanism.

[0034] This mechanism utilizes the smoothing and peak-shaving effect of reservoirs on inflow to calculate the surface runoff and interflow runoff of surface water and shallow soil water, which evolve after physical delay and water regulation driven by gravity and topographic slope.

[0035] Step 104: Construct a groundwater storage balance model. Based on the groundwater storage at the previous moment, use a state adaptive mechanism to smoothly adjust the nonlinear exponent of the baseflow calculation according to the saturation state of the groundwater storage at the previous moment, and calculate the nonlinear groundwater baseflow (groundwater discharge). Based on the difference between the groundwater runoff and the nonlinear groundwater baseflow, recursively update the groundwater storage at the previous moment to obtain the groundwater storage at the current moment.

[0036] To address the technical problem of groundwater runoff exhibiting a single, fixed exponential decay in traditional hydrological models.

[0037] Specifically, an independent groundwater storage balance model is extracted and constructed within the system. This model receives the groundwater runoff generated in the aforementioned steps as the sole positive recharge source.

[0038] By accumulating the difference between replenishment and discharge, the groundwater storage is updated differentially over time, transforming it into a dynamic indicator reflecting the abundance or scarcity of deep groundwater in the basin.

[0039] Furthermore, by adopting a state-adaptive mechanism, the nonlinear exponent used to control the baseflow discharge velocity is no longer a constant, but becomes a dependent variable of the groundwater storage saturation state.

[0040] When the aquifer is close to saturation, this mechanism can adaptively increase the nonlinear exponent to stimulate a rapid discharge response during the flood peak; when the aquifer is extremely dry, the exponent is smoothly decreased, causing the baseflow to enter a very slow decay mode.

[0041] Based on this dynamically adjusted nonlinear exponential output, the nonlinear groundwater base flow rate is mathematically reproduced to represent the complex transport characteristics of groundwater throughout the entire hydrological cycle.

[0042] Step 105: The surface runoff, interflow runoff, and nonlinear groundwater baseflow are superimposed and processed by river network lag propagation calculation to obtain the total runoff process of the watershed.

[0043] Through the above steps, shallow linear calculations and deep nonlinear calculations are completed. Subsequently, under the same volumetric flow rate dimensions, a simultaneous summation operation can be performed on surface runoff, interflow runoff, and nonlinear subsurface baseflow.

[0044] In some scenarios, considering the resistance factors such as riverbed roughness and channel tortuosity that affect water flow transmission in natural river networks, the superimposed composite flow data is fed into a pre-defined river network lag propagation calculation module for waveform evolution processing. This module can employ the Muskingan method or its equivalent linear channel confluence method.

[0045] Through the above calculation process, the propagation time delay and wave crest flattening phenomenon of flood waves advancing downstream in the physical river channel can be reproduced, and the total runoff process of the basin can be obtained that matches the physical state of the actual control section at the basin outlet.

[0046] Optionally, the river network lag propagation calculation can be implemented using the Muskingan method. This method is based on the wedge-shaped and prismatic water storage characteristics of the river channel for flood waves. It uses the river section storage volume as a linear weighted function of the inflow and outflow within the time step, and iteratively calculates the physical delay and peak flattening of the flood wave as it propagates downstream along the river network.

[0047] Specifically, two parameters, flow proportion coefficient and storage propagation time, are introduced. The inflow at the previous time step and the current time step are weighted and combined using a differential iteration method to output the outflow at the downstream end of the river section.

[0048] Example 2: The specific steps involved in acquiring watershed precipitation and evapotranspiration data, performing spatial weighting and quality control processing, and obtaining the watershed average input sequence are further explained.

[0049] That is, in one possible implementation, watershed precipitation data and evapotranspiration data are acquired, and spatial weighting and quality control processing is performed on them to obtain the watershed average input sequence, including the following steps: Step 201: Spatial weighting calculation is performed on the watershed precipitation data and evapotranspiration data using the Thiessen polygon method or the weighted average method to obtain spatial weighted data.

[0050] Precipitation and evapotranspiration data provided by meteorological stations are discrete, point-like spatial characteristics, while watershed hydrological simulation requires driving data at the areal scale. To eliminate model input errors caused by uneven spatial distribution, it is necessary to convert station observations into watershed mean values.

[0051] Specifically, after obtaining a base map of the watershed containing multiple observation stations, a spatial geometric segmentation network is constructed using the Thiessen polygon method. This method uses each meteorological station as the centroid of a polygon, divides the effective representative area of ​​each station according to the geometric bisector of the spatial distance, and calculates the ratio of this effective representative area to the total watershed area to obtain the corresponding area weight coefficient. Subsequently, the measured data of each station are multiplied by its corresponding area weight coefficient and summed to output spatially weighted data.

[0052] In some alternative implementations, if the target watershed area is small and the terrain is flat with uniform station distribution, a weighted average method can be used. This involves assigning equal weights to the observations at all stations and directly calculating the arithmetic mean as the spatially weighted data. The above processing is a mature spatial interpolation method in the industry, and the specific algebraic derivation formulas are not listed here.

[0053] Step 202 involves imputing missing data and removing outliers from the spatially weighted data, and identifying potential mutations through a consistency check to generate a stable and reliable watershed average input sequence.

[0054] Due to limitations such as power outages for field sensors, communication link failures, or environmental noise interference, the initially obtained spatially weighted data sequences often contain missing values ​​or physically meaningless extrema. Directly inputting data with these defects into a hydrological model can lead to divergence in iterative calculations or overflow of boundary conditions.

[0055] Therefore, quality control of the sequence execution system is necessary. By designing a data repair channel that includes interpolation and elimination logic, and combining it with a consistency verification method based on the statistical characteristics of time series, a data chain that meets the requirements of continuous computation can be reconstructed and solidified into the watershed average input sequence that ultimately participates in the model's runoff calculation.

[0056] Step 203: For missing data, if it is determined to be a short-term missing data, then linear interpolation is used to interpolate and obtain short-term interpolated data. If it is determined to be a long-term missing data, the measured data of the adjacent observation station is obtained as the neighbor station data, and the neighbor station data substitution method is used to interpolate to obtain the long-term interpolated data.

[0057] Because the duration of missing data varies greatly, using a single imputation algorithm would cause data distortion. Therefore, this step introduces a conditional branch decision mechanism based on time period length comparison.

[0058] A threshold for the duration of missing data is pre-configured in the system, such as 24 hours, or other values ​​can be designed as needed. When the system detects a continuous missing value window in the spatially weighted data, it calculates the time span of the window and compares it with the preset time period threshold.

[0059] In the first execution path, if the missing data is determined to be short-term, linear interpolation is used for interpolation.

[0060] Specifically, when the gap span is less than 24 hours, it is determined to be a short-term missing data period. At this time, the physical processes of rainfall or evapotranspiration have not changed drastically. The system extracts the previous and next valid observation values ​​before the missing data gap, calculates the slope of the numerical change between the two, and distributes the slope proportionally to each time point within the gap according to the time step to generate fill-in values.

[0061] In the second execution path, if the missing data is determined to be missing for a long period of time, the data from neighboring stations will be used for interpolation.

[0062] Specifically, when the gap span is ≥24 hours, it is determined to be a long-term missing data period. Since the autocorrelation of meteorological elements decays over long periods, the system extracts several candidate stations closest to the target station and calculates the correlation coefficient between the target station and the candidate stations during their historical synchronous observation periods. The neighboring station with the highest correlation coefficient is selected as the reference benchmark. A scaling factor is calculated based on the ratio of their historical averages. The measured data of the reference neighboring station during the same missing data period are multiplied by this scaling factor to complete the long-term missing data filling for the target station segment.

[0063] Step 204: Use a statistical boundary detection method based on the standard deviation principle to identify outliers in the spatially weighted data and correct them to obtain the corrected spatially weighted data.

[0064] In addition to missing data, outliers with anomalous jumps in the sequence also need to be intercepted.

[0065] In practice, the system calculates the overall mean and standard deviation of the data within a fixed-time sliding window. Based on the Gaussian distribution characteristics, the tolerance boundary is set to the mean plus or minus three times the standard deviation. The measured values ​​are compared hourly. When the absolute value of the data at a certain time exceeds the statistical tolerance boundary, the system determines it as an outlier caused by instrument malfunction and removes it from the original sequence, generating an artificial gap.

[0066] The system invokes the interpolation path described in step 203 above to recalculate the values ​​for the artificially created window. Further, after interpolation and correction, the Mann-Kendall nonparametric statistical test is introduced to assess the mean stability of the long sequence, identify trend abrupt changes caused by climate change or station migration, perform piecewise smoothing, and finally output the watershed average input sequence to drive the hydrological model.

[0067] Example 3: Based on the above embodiment one, as follows Figure 2 As shown, this embodiment further explains how to calculate the total runoff based on the average input sequence of the watershed, and how to divide the total runoff into specific stages of surface runoff, interflow, and groundwater runoff in combination with the free water storage structure.

[0068] In one possible implementation, the following steps are included: Step 301: The net rainfall is obtained by calculating the evapotranspiration capacity in the average input sequence of the watershed, and the soil moisture status is updated according to the preset three-layer soil structure.

[0069] Evapotranspiration is the main pathway for water consumption in a watershed. Potential evapotranspiration data is obtained from the average input sequence of the watershed, and this data is multiplied by a pre-configured conversion factor to calculate the actual evapotranspiration capacity.

[0070] Subtracting the actual evapotranspiration capacity from the average precipitation in the watershed yields the net rainfall, which characterizes the effective water input to the surface. To depict the characteristics of water transport in the vertical soil profile, the system pre-defines a three-layer soil structure consisting of an upper layer, a lower layer, and a deep layer.

[0071] After obtaining the net rainfall, following the physical principle of prioritizing consumption in the upper layer and progressive replenishment in the lower and deeper layers, the evapotranspiration consumption and water replenishment of each layer are calculated, and the soil moisture status of each layer is dynamically updated accordingly, providing real-time soil bottom water conditions for subsequent runoff calculations.

[0072] Step 302: If the net rainfall is less than or equal to zero, the total runoff is set to zero. Under the condition of having net rainfall, the total runoff is calculated based on the full runoff generation mechanism and combined with the soil moisture status, and the current runoff generation area ratio is determined according to the full runoff generation mechanism.

[0073] The physical significance of the saturation runoff mechanism lies in the fact that runoff will only occur after the water storage capacity of the soil vadose zone reaches saturation.

[0074] Conditions are determined based on the net rainfall calculated above. When the calculated net rainfall is less than or equal to 0, it is determined that the current period does not meet the conditions for runoff generation, and the total runoff is directly set to 0.

[0075] When the net rainfall condition is met (i.e., net rainfall is greater than 0), the system extracts the current watershed water storage from the updated soil moisture state and adds it to the net rainfall to obtain the theoretical total water storage. The theoretical total water storage is then compared with the preset average watershed water storage capacity for judgment.

[0076] For example, if the theoretical total water storage is less than or equal to the average water storage capacity of the basin, the basin is determined to be under unsaturated conditions. In this case, only some areas within the basin generate runoff due to full storage. The system, taking into account the non-uniform distribution characteristics of the basin's water storage capacity curve, calculates the runoff values ​​corresponding to some runoff-producing areas as the total runoff.

[0077] If the theoretical total water storage capacity exceeds the average water storage capacity of the basin, the basin is determined to be under full saturation conditions. At this time, all areas within the basin are at full capacity, and the system converts the incoming net rainfall and the remaining water volume after deducting soil water shortage into total runoff.

[0078] Step 303: Based on the current free water storage state obtained from the previous time step iteration output or preset initialization, the total runoff is dynamically allocated, and combined with the runoff area ratio, the total runoff is divided into surface runoff, interflow and groundwater runoff.

[0079] After obtaining the total production flow, it needs to be separated into runoff components with different response rates according to the pore structure of the underlying surface.

[0080] A free water storage structure was constructed to characterize the total amount of gravity water in shallow soil pores. The total runoff generated at the current time step was added to the initial free water storage at the previous time step to calculate the updated free water storage state. The proportion of surface runoff generation is positively correlated with the degree of soil free water saturation.

[0081] Based on the ratio of the current free water storage state to the preset free water storage capacity, combined with the curve index describing the spatial non-uniformity of the watershed's underlying surface water storage capacity, and combined with the proportion of runoff-producing area, the surface runoff is calculated and allocated.

[0082] For the amount of water infiltrating into the soil, the system calculates and separates the interflow moving laterally along the shallow soil layer based on the preset interflow outflow coefficient and the current free water storage state. Similarly, the system calculates and separates the groundwater runoff that infiltrates vertically into the deep aquifer based on the preset groundwater runoff outflow coefficient and the current free water storage state.

[0083] Through the dynamic allocation process described above, the total runoff is physically isolated into three groups of water sources with independent flow velocities, providing input data for subsequent classification and runoff calculations. The runoff and water source division calculations involved in this step are based on well-known algorithms in existing hydrological models, and their core logic has been described in the textual action above; therefore, the specific algebraic mathematical equations will not be listed here.

[0084] Example 4: Further explanation of the specific steps involved in using a linear regulation and storage mechanism to calculate the runoff and obtain the surface runoff and interflow runoff.

[0085] In one possible implementation, the following steps are included: Step 401: Surface runoff and interflow, measured in units of runoff depth, are uniformly converted into flow rates within the corresponding time step.

[0086] In hydrophysical calculations, surface runoff and interflow obtained through runoff generation mechanisms are typically measured using water depth as the unit of measurement. However, subsequent nonlinear baseflow calculations and multi-source superposition processes require flow rate per unit time as the fundamental physical quantity. If different water sources are summed while maintaining their own independent units, it will lead to computational divergence and physical inconsistencies in the system.

[0087] Therefore, the current model calculation time step parameter is extracted, and the surface runoff and interflow in the water depth dimension are converted into flow parameters constrained by the calculation time step.

[0088] Specifically, this conversion process maintains the spatial area attribute of the water volume parameter, distributing it evenly across the current time step, and outputting a combination of the watershed area parameter and the calculation time step. For example, the runoff depth, measured in millimeters, is uniformly converted into a flow rate value measured in cubic meters per second within the corresponding time step, providing a standardized data interface for subsequent dimensional superposition of the three water sources.

[0089] Step 402: Introduce the surface flow recession coefficient, use a linear reservoir model to iteratively calculate the surface runoff converted into flow rate, and output the surface runoff flow rate.

[0090] As surface runoff converges on slopes and moves toward rivers, it is affected by a combination of factors, including micro-topographical undulations, vegetation blocking water, and slope friction, resulting in physical lag and storage characteristics.

[0091] Obtain a pre-configured surface flow retreat coefficient, which characterizes the retention capacity of the underlying slope of the watershed for surface water flow.

[0092] The surface runoff, converted to flow rate, is used as the system inflow at the current time step, and a linear reservoir model is invoked for calculation.

[0093] In practice, the historical confluence state value calculated from the previous time step is extracted, and the historical confluence state value is weighted and attenuated using the surface flow decay coefficient. The attenuated value is then summed with the portion of the system inflow in the current time step that was not reduced by the decay coefficient.

[0094] By using the above differential iteration method, the evolution waveform of surface water flow on the slope is simulated, and the surface runoff flow affected by physical regulation is output.

[0095] Step 403: Introduce the interflow recession coefficient, use a linear reservoir model to iteratively calculate the interflow converted into flow rate, and output the interflow confluence flow rate.

[0096] Compared to surface runoff, interflows move laterally through shallow soil pores, encountering greater medium resistance, having a slower flow velocity, and exhibiting a stronger regulation and storage effect in their confluence paths.

[0097] Obtain the pre-configured interflow recession coefficient. Due to the strong damping characteristics of porous media, the value of this interflow recession coefficient is set to be greater than that of the surface flow recession coefficient.

[0098] A linear reservoir model architecture, identical to that used for surface runoff, is adopted, with interflow converted to flow rate as the input boundary condition. The historical confluence state values ​​of interflow at the previous time step are weighted using the interflow recession coefficient and linearly combined with the current inflow variable.

[0099] This iterative calculation process transforms the pulsed water infiltration into the shallow soil into a flow sequence with a gentle receding tail, outputting the interflow runoff. Through the independent linear regulation calculation in this embodiment, the two hydrophysical processes with different response rates at the surface and in the shallow soil are successfully decoupled, and data preparation under a unified dimension is completed.

[0100] Example 5: Figure 3 As shown in the figure, this embodiment further details the specific steps of constructing a groundwater storage balance model and dynamically updating groundwater storage based on groundwater runoff.

[0101] In one possible implementation, the following steps are included: Step 501: Obtain historical groundwater depth data and measured soil moisture content data for the target watershed.

[0102] When starting the groundwater storage balance model, a prerequisite condition reflecting the initial geological water storage state of the watershed is required, namely, cold start parameters.

[0103] Historical groundwater depth data was acquired through a network of groundwater level monitoring wells deployed within the target watershed, while measured soil moisture content data from the deep vadose zone was read using a soil moisture sensor array. This well and sensor data objectively recorded the physical state of the aquifer at the start of the model calculation, providing boundary constraints for subsequent continuous water volume calculations.

[0104] Step 502: Based on historical groundwater depth data and measured soil moisture content data, the initial groundwater storage is determined by water volume conversion, which serves as the starting point for the historical groundwater storage at the first time step of the groundwater storage balance model.

[0105] Directly obtained groundwater depth and soil moisture content are usually expressed in spatial distance or volume percentage dimensions, and cannot be directly used in water balance calculations. Therefore, by combining the specific yield and porosity parameters of the aquifer in the target watershed, discrete historical groundwater depth data and measured soil moisture content data are integrated by volume and converted to equivalent water depth.

[0106] Through the above conversion process, the initial groundwater storage is obtained in terms of physical depth. The specific value range is usually set to 50 to 200 mm, and it is loaded into the computing memory as the starting point of the historical groundwater storage for the first time step of the groundwater storage balance model.

[0107] Step 503: In the groundwater storage balance model, a differential form of the water balance equation is used, and groundwater runoff is taken as the groundwater recharge.

[0108] Traditional hydrological models often rely on surface runoff generation modules for baseflow calculations, lacking an independent description of groundwater transport. This embodiment constructs an independent groundwater storage control volume within the system and performs calculations using a differential form of the water balance equation for this control volume.

[0109] At the input boundary of the equation, the groundwater runoff obtained from the aforementioned steps is redirected and set as the positive groundwater recharge of the groundwater reservoir, thereby realizing the physical channel connection for the migration of surface seepage water to deep aquifers in the model architecture.

[0110] Step 504: Calculate the nonlinear groundwater base flow rate at the current time step based on the groundwater storage at the previous time step; calculate the water volume difference between the groundwater recharge and the nonlinear groundwater base flow rate at the current time step.

[0111] After establishing the input and output terms, the system performs an algebraic comparison of the water volume entering and leaving the control volume within one computational time step. It extracts the groundwater recharge allocated for the current time step and simultaneously acquires the nonlinear groundwater baseflow calculated by the nonlinear equation.

[0112] The water volume difference is calculated by subtracting the discharge from the replenishment amount. This water volume difference characterizes the net water storage change of the aquifer within the current time window. A positive difference indicates that the aquifer is in the water storage period, while a negative difference indicates that the aquifer is in the water loss period.

[0113] Step 505: Based on the water volume difference, iteratively accumulate the historical groundwater storage at the previous moment to realize the recursive update of groundwater storage.

[0114] The state variables are updated by calling the differential form of the water balance formula.

[0115] Specifically: W a (t)=W a (t-1)+Q r (t)-Q b (t); In the formula, W a (t) represents the updated groundwater storage at the current moment, W a (t-1) represents the historical groundwater storage at the previous moment, Q r (t) represents the current groundwater recharge, Q b (t) represents the nonlinear underground base flow rate at the current moment.

[0116] This formula is used to algebraically add the calculated water volume difference to the historical groundwater storage, and output the updated groundwater storage parameters to ensure the water volume conservation of the groundwater reservoir in long-term series simulation.

[0117] Step 506: Obtain the preset boundary values ​​of minimum and maximum groundwater storage.

[0118] The water storage capacity of an aquifer is limited by its geological substrate structure and maximum pore space, exhibiting a definite physical limit. The system pre-determines the minimum and maximum groundwater storage boundary values ​​in the parameter configuration table.

[0119] The minimum groundwater storage corresponds to the dead water level of the aquifer during the dry season; for example, a reasonable value range is between 30 and 80 millimeters. The maximum groundwater storage corresponds to the water volume when the entire aquifer profile reaches saturation.

[0120] The specific values ​​of the above parameters can be determined by those skilled in the art based on the hydrogeological survey data of the target watershed. The above boundary values ​​serve as a reference for subsequent data standardization processing.

[0121] Step 507: Use the minimum and maximum groundwater reserves to perform extreme value normalization on the groundwater reserves at the current moment, and extract the relative saturation of groundwater in the range of zero to one.

[0122] Due to the significant differences in absolute water content across different watersheds, directly using the absolute value of groundwater storage in nonlinear calculations would cause the model parameters to lose their spatial generalization ability. Therefore, an extreme value normalization algorithm is employed to perform dimensionality reduction mapping of the storage state.

[0123] Specifically: S(t) = (W a (t)-W a_min ) / (W a_max -W a_min ); Where S(t) is the relative saturation of groundwater, W a (t) represents the current groundwater storage, W a_min W is the preset minimum groundwater storage. a_max The maximum groundwater storage is preset. The system calculates the difference between the current storage and the minimum storage as the numerator, and the difference between the maximum storage and the minimum storage as the denominator. The dimensionless relative saturation of groundwater is obtained by calculating the ratio.

[0124] The relative saturation is strictly limited to a continuous closed interval between zero and one. When the value approaches zero, it indicates that the watershed is in an extremely dry state; when the value approaches one, it indicates that the groundwater aquifer is in a state of near saturation. This state parameter serves as the standardized control variable for the subsequent state-driven adaptive nonlinear mechanism.

[0125] Example 6: Figure 4 As shown in the figure, this embodiment further explains the specific steps of using a state adaptive mechanism to smoothly adjust the nonlinear exponent of the baseflow calculation according to the saturation state of the groundwater storage, and to calculate the nonlinear groundwater baseflow.

[0126] In one possible implementation, the following steps are included: Step 601: Set the nonlinear index as a dynamic state variable that depends on the relative saturation of groundwater. Using a preset S-shaped continuous function, with the relative saturation of groundwater as the input variable, perform a smooth transition mapping between the preset minimum nonlinear index and the maximum nonlinear index to calculate the nonlinear index that dynamically adjusts with the relative saturation of groundwater.

[0127] Traditional linear reservoir models employ constant decay parameters, which cannot accommodate the accelerated discharge state of groundwater aquifers at high water levels and the slow decay state at low water levels. To resolve this physical conflict, the system sets the exponential parameter controlling the discharge rate as a dynamic state variable. Specifically, the system uses an sigmoid continuous function with bidirectional asymptotic properties to construct a state mapping operator. This mapping operator receives the relative saturation of groundwater calculated in the preceding steps as an independent variable, and controls the degree of nonlinear distortion of the curve through internal sensitivity coefficients and inflection point parameters, enabling the output variable to smoothly transition between preset upper and lower boundaries.

[0128] The specific calculation formula is as follows: α(S) = α min +(α max -α min ) / (1+exp[-β*(S-S0)]); In the formula, α(S) is a dynamically adjusted nonlinear exponent, α min α is the preset minimum nonlinear exponent. max exp represents the natural exponential function operation, β is the state sensitivity coefficient used to control the steepness of the nonlinear exponent change with saturation, S is the relative saturation of groundwater at the current moment, and S0 is the preset saturation state inflection point.

[0129] For example, in the first scenario, after continuous rainfall in the target watershed, the aquifer is close to saturation. At this point, the relative saturation of the groundwater is in a numerical range far above the inflection point of saturation. Due to the characteristic of the S-shaped continuous function to rapidly approach the upper asymptote above the inflection point, the calculated nonlinear exponent will be close to the preset maximum nonlinear exponent, causing the baseflow to exhibit a high-order accelerated release response, consistent with the physical phenomenon of accelerated discharge caused by the increased water head driving force when the aquifer is saturated.

[0130] In the second scenario, when the target watershed is at the end of the dry season, the aquifer level drops significantly to an extremely arid state. At this time, the relative saturation of groundwater is far below the inflection point of saturation. The S-shaped continuous function slowly approaches the lower asymptote below the inflection point, and the calculated nonlinear exponent is close to the preset minimum nonlinear exponent. The baseflow decay rate slows down significantly, and it enters a slow and continuous release state, avoiding premature cessation of baseflow.

[0131] For example, assuming the system configuration parameters are a minimum nonlinearity index of 0.8, a maximum nonlinearity index of 2.0, a state sensitivity coefficient of 10, and a saturation inflection point of 0.5, when the aquifer is close to saturation (relative groundwater saturation of 0.8), substituting these parameters into the above formula yields a dynamic adjustment nonlinearity index of approximately 1.94. This value is much greater than 1, indicating that the baseflow will exhibit accelerated release characteristics. When the aquifer is extremely dry (relative groundwater saturation of 0.2), the calculated dynamic adjustment nonlinearity index is approximately 0.86. This value is less than 1, indicating that the baseflow decay rate slows down, transitioning to a slow release state.

[0132] In other alternative implementations, the sigmoid continuous function can also be replaced with other continuous monotonically increasing mapping functions with bidirectional asymptotic properties, such as linear transformations of the arctangent function or hyperbolic tangent function. Those skilled in the art can select a suitable mapping function form based on the response characteristics of the actual watershed.

[0133] Step 602: Obtain the pre-configured baseflow coefficient and the preset reference constant for balancing dimensions. Use the preset reference constant to perform dimensionless calculation on the groundwater storage to obtain the base term. Perform a power function operation on the base term and the nonlinear exponent, and multiply the result of the power function operation with the baseflow coefficient to obtain the nonlinear groundwater baseflow.

[0134] After introducing a dynamic nonlinear exponent, directly performing a power function calculation on groundwater storage with spatial depth dimensions will cause the comprehensive dimensions on the right side of the baseflow equation to continuously drift as the nonlinear exponent changes hourly, resulting in an inability to balance the dimensions on both sides of the equation. To eliminate this algebraic logic conflict, the system extracts a preset reference constant from a pre-injected parameter table. The current groundwater storage is divided by this preset reference constant, and this ratio, which has physical dimensions, is defined as the dimensionless base term. Subsequently, using this base term as the base and the dynamically adjusted nonlinear exponent as the exponent, a power function evaluation is performed. Finally, the baseflow coefficient, characterizing the permeability of the aquifer, is obtained, multiplied by the power function evaluation result, and the nonlinear groundwater baseflow rate within the current time step is output.

[0135] Specific: Q b (t)=K g *(W a (t) / W ref )α(S) ; Among them, Q b (t) represents the nonlinear underground baseflow at the current moment, K g For the pre-configured base current coefficient, W a (t) represents the current groundwater storage, W ref α(S) is a preset reference constant used for balancing dimensions, and α(S) is a dynamically adjusted nonlinear exponent.

[0136] In actual time-by-time iterative calculations, the groundwater storage W at the current moment... a (t) Take the storage value W updated after water balance at the previous moment. a (t-1), or the intermediate state of the reserves before discharge after replenishment at the current time can be used. Those skilled in the art can choose the specific time series processing method according to the calculation accuracy requirements.

[0137] In one optional implementation, the preset reference constant used for balancing dimensions is set to a fixed constant unit with a value of 1 and the same depth dimension, so that the value of the base term is equal to the absolute value of the groundwater storage, thereby ensuring that the base is greater than 1, and ensuring that the result of the power function operation increases accordingly when the nonlinear exponent increases, which is consistent with the physical mechanism of accelerated discharge under the high storage state of the aquifer.

[0138] Step 603: Apply extreme value constraints to the groundwater storage, so that its calculated value is limited to the minimum and maximum groundwater storage, and obtain the constrained groundwater storage; apply extreme value constraints to the nonlinear groundwater baseflow, so that its output value is not less than the preset minimum positive baseflow and not greater than the preset maximum baseflow, and obtain the constrained nonlinear groundwater baseflow.

[0139] Alternatively, one can apply an exponential extremum constraint to the nonlinear exponent, limiting its calculated value to between the minimum and maximum nonlinear exponents, thus obtaining a constrained nonlinear exponent.

[0140] To ensure that all parameters remain within a reasonable hydrophysical range during the nonlinear calculation process, the system embeds three-dimensional hydrophysical boundary constraints at the end of the nonlinear base current calculation cycle.

[0141] If the temporary calculated value of groundwater storage is determined to be lower than the numerical limit used to limit the minimum water storage during the dry season, i.e. the preset minimum groundwater storage, then the value is forcibly assigned to the minimum storage. If the value is determined to be higher than the numerical limit used to limit the saturated water storage of the aquifer, i.e. the preset maximum groundwater storage, then the overflow water volume is forcibly stripped and assigned to the maximum storage.

[0142] In one possible implementation, the hydrophysical boundary constraints applied to each parameter include a minimum nonlinear exponent ranging from 0.6 to 1.0 and a maximum nonlinear exponent ranging from 1.2 to 2.5. Furthermore, the minimum groundwater storage capacity is set as the numerical limit of the minimum water storage capacity during the corresponding dry season, and the maximum groundwater storage capacity is set as the numerical limit of the saturated water storage capacity of the corresponding groundwater aquifer.

[0143] Specifically, for the nonlinear exponent, the system extracts pre-configured intervals for the minimum and maximum nonlinear exponents. In actual parameter space calibration, the minimum nonlinear exponent is limited to a range of 0.6 to 1.0, and the maximum nonlinear exponent is limited to a range of 1.2 to 2.5. The system detects the output value of the sigmoid continuous function, and performs boundary truncation when it exceeds these set intervals to ensure the physical rationality of the exponent values.

[0144] For the final output of nonlinear subsurface baseflow, the system sets an absolute flow rate threshold. The minimum positive baseflow and maximum baseflow, calibrated based on historical extreme hydrological events in the watershed, are obtained, and bidirectional truncation filtering is applied to the calculated baseflow values. When the calculated value is biased towards the negative domain, it is reset to the minimum positive baseflow; when the value exceeds the historical extreme value limit, it is suppressed to the maximum baseflow. Through this three-dimensional boundary interception operation, the system maintains its underlying numerical stability while leveraging the fitting advantages provided by the nonlinear mechanism.

[0145] Example 7: Providing an offline build branch independent of the online process.

[0146] Before the rainfall-runoff simulation model is officially put into real-time prediction application, this embodiment further pre-constructs and calibrates the various calculation parameters required for the method execution.

[0147] In one possible implementation, the following steps are included: Step 701: Obtain historical measured flow data of the watershed and the corresponding historical input meteorological data.

[0148] Hydrological models contain a large number of empirical and conceptual parameters that characterize the physical properties of the underlying surface and the hydrogeological conditions of aquifers. These parameters cannot be obtained directly through physical measurements and need to be derived by inversion using historical observation sequences.

[0149] By connecting to a historical hydrological and meteorological database, a representative long-term historical time window is defined, which typically includes a complete hydrological cycle of multiple high-water years, normal-water years, and low-water years.

[0150] Historical measured flow data of the watershed outlet control section within the time window are extracted. At the same time, historical input meteorological data such as precipitation and evapotranspiration covering the entire watershed are also extracted. These historical records are loaded into the memory space of the offline calibration environment as the data basis for model parameter inversion.

[0151] Step 702: Based on historical input meteorological data, execute the preceding rainfall-runoff simulation process to calculate historical simulated flow data.

[0152] In the offline calibration environment, the system initializes the hydrological model using an initial set of random parameters. Subsequently, the system uses the acquired historical meteorological data as the driving condition to initiate the complete physical calculation process in Examples 1 to 6, including runoff generation classification, shallow linear storage and regulation, and deep dynamic nonlinear baseflow calculation.

[0153] Through this forward physics calculation process, the system outputs hourly runoff forecasts for a given time period, which are defined as historical simulated flow data. This data reflects the theoretical response of the current parameter set to specific historical meteorological inputs.

[0154] Step 703: Based on historical measured flow data of the watershed and combined with historical simulated flow data, construct a multi-objective evaluation function. The multi-objective evaluation function includes the Nash efficiency coefficient index for evaluating the process fitting accuracy and the flood volume relative error index for evaluating the water balance state.

[0155] Since a single evaluation index cannot adequately account for both the high peak flow during flood season and the low flow tail during dry season, the system employs a multi-objective evaluation framework to quantitatively assess the performance of the model's current parameter set. The system introduces the Nash efficiency coefficient to evaluate the overall fit between the simulated and measured flow sequences; its calculation formula is as follows: NSE=1-∑((Q obs (t)-Q sim (t)) 2 ) / ∑((Q obs (t)-Q obs_mean ) 2 ); Where NSE is the Nash efficiency coefficient, Q obs (t) represents the historical measured flow data of the basin at the current moment, Q sim (t) represents the historical simulated flow data at the current moment, Q obs_mean ∑() represents the arithmetic mean of the measured flow data of the historical watershed within a given time period, and ∑() represents the summation operation of the continuous time series within the time period.

[0156] Meanwhile, the system introduces the relative error index of flood volume to evaluate the overall water volume conservation of the model within a time period, and its calculation formula is as follows: RE=(∑(Q sim (t))-∑(Q obs (t))) / ∑(Q obs (t))*100%; Here, RE stands for Relative Error of Flood Volume. The system linearly weights and combines the two indicators, NSE and RE, according to pre-configured weighting coefficients to construct a comprehensive multi-objective evaluation function. This function, when guiding parameter evolution, both constrains the model to capture the temporal pattern of the flood peak and forcibly controls the total water volume deviation during the long-term recession process.

[0157] Step 704: A preset global composite evolutionary optimization algorithm is adopted to automatically iterate and optimize the parameter space involving the flow generation mechanism and nonlinear confluence mechanism, with the goal of satisfying the comprehensive optimal convergence condition of the multi-objective evaluation function, and outputs a pre-calibrated model parameter set.

[0158] Due to the introduction of a state-adaptive nonlinear baseflow mechanism, the parameter space dimension of the model increases and the degree of nonlinearity improves, making traditional local gradient descent algorithms prone to getting trapped in local optima. The system employs a pre-defined global composite evolutionary optimization algorithm to drive parameter search. Specifically, the system randomly generates multiple initial parameter populations within the parameter space to be optimized. This parameter space encompasses the evapotranspiration conversion coefficient, interflow runoff coefficient, surface flow recession coefficient, interflow recession coefficient, and baseflow coefficient K. g The core variables include the state sensitivity coefficient β, the saturation state inflection point S0, and the maximum and minimum nonlinear exponents that control the nonlinear morphology.

[0159] The system calls the aforementioned multi-objective evaluation function to assess the fitness of individuals in the population, and performs compound parallel search, subpopulation information exchange, and mutation operations based on the evaluation results. As the number of iterations increases, the parameter population gradually converges towards a high-fitness region. When the system determines that the change in the multi-objective evaluation function value over several consecutive generations is less than a preset convergence threshold, or the number of iterations reaches a preset computational upper limit, the system determines that the comprehensive optimal convergence condition is met. The system stops the iterative search, extracts the individual with the highest current fitness as the global optimal solution, decodes and solidifies it, and outputs a pre-calibrated model parameter set. This parameter set is then loaded into the online hydrological forecasting system to support real-time rainfall-runoff simulation calculations.

[0160] Example 8: Based on Examples 1 to 7 above, this example provides a branch of flood control scheduling application based on runoff simulation.

[0161] In one possible implementation, after obtaining the total runoff process of the watershed, the following steps are also included: Step 801: Compare and analyze the total runoff process of the basin with the preset flood control warning threshold of the basin to obtain the result of the flood peak exceeding the limit.

[0162] The total runoff process of a watershed is a discrete sequence of flow data continuously distributed along a time axis, reflecting the flow change trend at the control sections of the watershed within the prediction time window. To transform the predicted data into flood control indicators with engineering guidance significance, preset watershed flood warning thresholds are retrieved in memory.

[0163] The preset flood control warning threshold for the basin is determined by the water resources management department based on the historical flow capacity of the downstream river and the design standards of the dikes along the river. During the comparative analysis, the system employs a time-series traversal algorithm to extract the simulated flow value corresponding to each time point in the total runoff process of the basin, and then compares the difference between this value and the basin flood control warning threshold.

[0164] In some alternative implementations, the flood control warning threshold for the watershed can be configured as a single constant flow extreme value parameter, or as a dynamic water level curve that varies with the flood season. The system dynamically loads the corresponding upper limit of the threshold based on the attribute of the current calculation time step.

[0165] Step 802: When the flood peak exceedance determination result indicates that the flood peak flow in the total runoff process of the basin exceeds the flood control warning threshold of the basin, a flood control warning signal containing the flood peak arrival time and the estimated flow is automatically generated, or the flood control scheduling strategy of the corresponding water conservancy hub in the basin is triggered.

[0166] Based on the comparison results of the aforementioned steps, the execution status is determined and subsequent service routing is distributed.

[0167] When it is detected that the simulated flow exceeds the basin flood control warning threshold during a specific period within the prediction time window, the warning condition is determined to be met.

[0168] The peak flow rate is extracted from the data range exceeding the limit using a peak-finding algorithm, and the timestamp corresponding to this maximum flow rate is recorded as the peak arrival time. Based on this condition, the system is configured with two response execution paths that can be selected concurrently.

[0169] For example, in the first execution path, an automatic early warning mechanism is executed.

[0170] The extracted peak arrival time and estimated flow rate are encapsulated to generate a structured flood warning signal. This flood warning signal can take the form of a warning communication message pushed to mobile terminals, or a time-series line graph rendered in red on the hydrological dispatch and monitoring console interface. This flood warning signal provides a data foundation for flood control duty personnel.

[0171] In the second execution path, the underlying control network of the water conservancy hub is connected to trigger the flood control scheduling strategy of the corresponding water conservancy hub in the basin.

[0172] Specifically, the system feeds data packets containing estimated flow and peak arrival time into the reservoir dispatch control center, which then converts the inflow into corresponding floodgate opening control parameters based on the flood control dispatch strategy.

[0173] The control parameters are transmitted to the field controller nodes of the corresponding reservoir via an industrial Ethernet interface. By opening the spillway or flood discharge channel in advance, the flood control capacity of the physical reservoir is emptied, thereby using pre-emptive storage to physically intercept the upcoming flood peak. Through this business extension, the numerical simulation results are transformed into physical control actions for the actual water conservancy project.

[0174] This invention provides a rainfall-runoff simulation method based on a nonlinear confluence mechanism of underground baseflow, comprising the following steps: Step S1: Rainfall and runoff data collection and preprocessing.

[0175] For the target watershed, basic data required for hydrological simulation are collected, mainly including precipitation, potential evapotranspiration, and measured runoff data. Precipitation data uses daily or hourly observations from rain gauges within the watershed, combined with satellite inversion data when necessary to enhance spatial representativeness; evapotranspiration data is calculated based on meteorological factors or directly uses observational data; runoff data selects measured flow rates at the watershed outlet control section for model calibration and validation.

[0176] After data compilation, spatial integration and quality control processing were performed. First, spatial weighting was applied to precipitation and evapotranspiration data. The basin-wide average precipitation was calculated using the Thiessen polygon or weighted average method. ; in t represents the spatially weighted average precipitation of the entire watershed at time t; n represents the total number of rain gauges (or precipitation data units) participating in the calculation within the watershed; i represents the index of the rain gauge (or precipitation data unit). Let be the precipitation observation value at time t for the i-th rain gauge station; Let be the watershed area represented by the i-th rain gauge station; The total area of ​​the target watershed is studied.

[0177] Evapotranspiration data were spatially weighted using the same method to ensure the spatial representativeness of the input data.

[0178] Secondly, data quality control and consistency processing were carried out: for missing data, linear interpolation was used for short periods and neighbor station substitution was used for long periods; outliers were identified and corrected based on the 3σ principle to reduce the impact of extreme errors; the data time resolution was unified; and the Mann-Kendall method was used to check the consistency of the sequences, identify potential mutations, and ensure that the data is stable and reliable.

[0179] Step S2: Runoff calculation and division of the three water sources.

[0180] Step S2.1: Calculate the potential evapotranspiration capacity based on the reference evapotranspiration rate and then calculate the net rainfall: ； ； in: This refers to the actual evaporation capacity; Potential evaporation; This is the conversion factor; Net rainfall; Let be the precipitation at time t. The evapotranspiration process follows the distribution principle of "priority to the upper layer, with progressive replenishment from the lower and deeper layers." The evapotranspiration of each layer is calculated, and the soil moisture state is updated.

[0181] Step S2.2: Under net rain conditions ( Based on the full-saturation runoff mechanism, the total runoff is calculated according to two scenarios based on the soil water storage status. :when When the condition is unsaturated, it corresponds to partial flow generation; when When the time is saturation condition, it corresponds to full-scale flow generation; when No runoff is generated at this time. ).in This represents the maximum soil water storage capacity. This process reflects the physical mechanism by which soil goes from water storage to saturation and then to runoff.

[0182] Step S2.3: After the runoff calculation is completed, based on the free water storage structure and the runoff area ratio... Total output flow Further subdivided into surface runoff ( ), soil flow ( ) and underground runoff ( The runoff is divided into three parts to reflect the formation mechanisms and confluence characteristics of different runoff components. The division process takes into account both net rainfall periods and periods without net rainfall periods, and dynamically allocates each runoff component based on the current free water storage status, thereby achieving a reasonable characterization of the watershed runoff composition and drainage process.

[0183] The formula for updating free water storage is: ; in, This represents the initial free water storage capacity. This is the updated free water storage capacity; Total output flow; , and These represent surface runoff, interflow, and groundwater runoff, respectively.

[0184] Based on this, the simplified formula for dividing the three water sources is as follows.

[0185] (1) Surface runoff; ; in, This refers to the free water storage capacity. The free water capacity curve index is used to describe the spatial non-uniformity of the water storage capacity of the underlying surface of a watershed. This formula shows that as the degree of soil free water storage increases, the proportion of surface runoff gradually increases, meaning that the closer the watershed is to saturation, the easier it is for surface runoff to form.

[0186] (2) Flowing through the soil; ; in, This is the interflow coefficient. Interflow refers to the discharge of a portion of the water in a free reservoir as shallow lateral flow. Its response rate is between that of surface runoff and groundwater runoff, and it has a significant impact on the flood recession process.

[0187] (3) Groundwater runoff; ; in, This is the outflow coefficient of groundwater runoff. Groundwater runoff refers to the slow flow process formed after free water seeps into deep soil and underground aquifers. Its main function is to maintain the base flow of the watershed and the river runoff during the dry season.

[0188] Step S3: Linear confluence of surface runoff and subsurface runoff.

[0189] Step S3.1: Convert the water volume in units of runoff depth to flow rate for subsequent runoff calculations. The conversion formula is: ; in, This corresponds to the flow rate of the water source; For runoff depth; The drainage area; The time step is 3.6; the coefficient 3.6 is the unit conversion factor.

[0190] Step S3.2: Surface runoff responds rapidly, mainly propagating along slopes and shallow river networks. Its confluence process can be described using a linear reservoir model. The calculation formula is: ; in, Let t be the surface runoff input flow rate; Let t be the output flow rate after merging. Output flow for the previous moment; It is the surface runoff concentration coefficient (dimensionless), reflecting the slope's water storage capacity.

[0191] Step S3.3: Intersoil flow moves laterally in the shallow soil layer. Its propagation speed is slower than surface runoff, but it can still be approximated as a linear storage process. A linear reservoir model is also used for calculation: ; in, Let t be the soil inflow rate. The output flow rate after the confluence of the soil streams; The interflow confluence coefficient (dimensionless) is usually greater than the surface runoff coefficient, reflecting its stronger lag.

[0192] Step S4: Nonlinear convergence of underground base flow.

[0193] Step S4.1: Groundwater storage reservoir Used to characterize the water storage capacity of groundwater aquifers in a watershed, its dynamic changes follow the principle of water balance, and the core water balance equation is in differential form: ; In the formula, The groundwater storage at time t reflects the water storage state of the groundwater aquifer, initial value. The value is determined by conversion based on historical groundwater depth data and soil moisture content data in the watershed, and the range is 50~200mm; The groundwater recharge at time t is the groundwater runoff obtained from step S3. Provided directly, i.e. To ensure consistency between the supply volume and the production flow process; : Groundwater discharge at time t (i.e., baseflow), which is the result of the improved nonlinear confluence calculation and is the main component of runoff during the dry season; t: time, consistent with the time step Δt of the model calculation; The time derivative represents the tiny increment of time t; This represents the rate of change of groundwater reserves over time.

[0194] To facilitate numerical calculations in the model, the differential form of the water balance equation is converted into a difference form, resulting in a recursive update formula for groundwater storage: ; In the formula: Δt represents the groundwater storage at time t-1; Δt is the time step of the model calculation to ensure the stability of the calculation results.

[0195] To characterize the saturation state of groundwater aquifers, a groundwater relative saturation index is introduced. : ; in: This represents the relative saturation of groundwater. and These represent the minimum and maximum groundwater reserves, respectively. After normalization... .when When the value approaches 0, it indicates that the aquifer is in a state of low water content and low storage capacity; when... When the value approaches 1, it indicates that the aquifer is close to saturation.

[0196] Step S4.2: To overcome the problem that traditional fixed exponential parameters cannot simultaneously adapt to the confluence characteristics of flood and dry seasons, nonlinear exponential parameters are used... Defined as a continuous function of groundwater storage state This allows the nonlinearity of the flow convergence to automatically adjust according to the aquifer saturation state. The improved baseflow expression is: ; In the formula, The baseflow coefficient characterizes the discharge capacity of a groundwater aquifer and is related to the hydrogeological conditions of the watershed (aquifer permeability, porosity, etc.). The higher the permeability coefficient, the greater the baseflow coefficient. The larger the value; It is a state-adaptive nonlinear exponential function. To ensure that the parameter changes are continuous, stable, and have a clear physical meaning, a sigmoid-type continuous function is used to construct it. : ; in: It is the minimum nonlinear exponent; The maximum nonlinear exponent; This is the state sensitivity coefficient, used to control the steepness of the nonlinear exponent as it changes with saturation. This is the inflection point for saturation, typically taken as 0.4 to 0.6. This function exhibits the following hydrophysical characteristics: (1) When the groundwater storage is low, that is hour, At this time, the baseflow decay rate slows down, which can more realistically depict the slow release process of groundwater during the dry season and improve the simulation capability of low flow rate. (2) When the groundwater storage is relatively high, that is hour, At this time, the baseflow shows an accelerating growth trend with the increase of storage, which can effectively describe the rapid discharge process under high groundwater storage during the flood season and improve the accuracy of flood peak simulation; (3) When the aquifer is in a state of medium reserves, exist and The transitions between them are continuous and smooth, thus avoiding the parameter mutation problem caused by traditional piecewise functions and improving the numerical stability and physical consistency of the model.

[0197] Compared to the traditional fixed parameter format: This method enables dynamic adaptive adjustment of the nonlinearity of underground baseflow, allowing the model to automatically switch between "slow-release" and "rapid-discharge" confluence mechanisms based on the groundwater storage status, thereby simultaneously improving the simulation capabilities during both flood and dry seasons.

[0198] Step S4.3: To ensure the reasonableness of the calculation results, the groundwater storage... and base current Set constraints: (1) Groundwater storage constraints: ,in This represents the minimum groundwater storage capacity, ranging from 30 to 80 mm, corresponding to the minimum groundwater storage capacity during the dry season. This represents the maximum groundwater storage capacity, ranging from 150 to 300 mm, corresponding to the saturated water storage capacity of the groundwater aquifer. (2) Nonlinear exponential constraint: ,in: The value range is 0.6 to 1.0; The value range is 1.2 to 2.5.

[0199] (3) Base current constraint: ,in The minimum base flow rate is set to a value ranging from 0.1 to 0.5 mm / Δt to avoid calculating a negative base flow rate. The maximum base flow rate is determined based on the historical maximum base flow data of the basin to ensure consistency with the actual hydrological process.

[0200] The nonlinear fundamental current at time t is obtained through the above calculations. This study simulates the groundwater runoff process, overcoming the limitations of traditional linear runoff models.

[0201] Step S5: Calculation of multiple water source confluence and river network confluence.

[0202] Step S5.1: Calculate the obtained surface runoff runoff concentration. , soil confluence flow , and nonlinear base current By superimposing the results, the total runoff of the basin can be obtained using the following formula: ; In the formula: Let t be the total runoff of the basin at time t, representing the total runoff after all water sources converge.

[0203] Step S5.2: River network runoff calculation. Based on the multi-source runoff, and considering the river channel regulation and propagation processes, the runoff is calculated using the lag calculation method: ; in, Let be the flow rate at the outlet section of the basin at time t; The river network water storage and receding coefficient reflects the river's regulation and storage capacity; River network confluence delay indicates the propagation time of runoff in the river channel.

[0204] Step S6: Model parameter calibration and result output.

[0205] Step S6.1: The parameter system of the model includes four parts: evapotranspiration, runoff generation, water source delineation, linear runoff parameters in the surface soil, and nonlinear runoff parameters in the groundwater baseflow.

[0206] Step S6.2: The SCE-UA global optimization algorithm is used to automatically optimize the parameters. During the calibration process, the measured flow rate at the watershed outlet section is used. For reference, a multi-objective evaluation function is constructed to comprehensively evaluate the model performance from three aspects: overall accuracy, water balance, and performance during the dry season. Specifically, it includes: (1) Nash efficiency coefficient (NSE). Used to evaluate the degree of fit between simulated and measured flow rates; (2) Relative error of flood volume (RE). Used to evaluate the deviation between simulated total runoff and measured total runoff.

[0207] The closer RE is to 0, the better the model's water conservation. Furthermore, the relative error of runoff during the dry season is specifically used to evaluate the model's simulation capability during the dry season.

[0208] The simulation results of the method of this invention are compared and analyzed with the simulation results of the traditional linear baseflow Xin'anjiang model (α=1). The focus is on evaluating the improvement effect of the model in terms of peak flow simulation, recession process simulation and dry season runoff reproduction. The effectiveness and superiority of the method are verified by comparison.

[0209] Based on the same inventive concept as the above embodiments, this embodiment provides a verification process for evaluating the performance of rainfall-runoff simulations. This verification process, independent of the steps defined in the claims, is used to quantitatively analyze the adjustment state of the simulation output data after introducing a state-adaptive nonlinear baseflow calculation mechanism, providing physical and statistical support data that sufficiently demonstrates the present invention. Specifically, it includes the following data extraction and comparison calculation steps.

[0210] Meteorological and runoff data for the target watershed within a preset time period are acquired and input into a first model configured with a linear baseflow mechanism and a second model configured with a state-adaptive nonlinear baseflow mechanism, respectively, for calculation to obtain first and second simulated data. In practice, historical rainfall-runoff observation sequences are extracted from a specific humid southern watershed, and these sequences are divided into typical flood periods and typical dry season periods based on hydrological characteristics. The nonlinearity exponent of the baseflow calculation module in the first model is fixed at a constant of 1.0. The second model is configured to execute the methods described in Embodiments 1 to 6 of this invention. The first and second models are run, and the first and second simulated data corresponding to the typical flood periods and typical dry season periods are output, respectively.

[0211] The Nash efficiency coefficient and relative error of flood volume were used as quantitative evaluation indicators to calculate and compare the output data of the first and second simulations. Simulation results for typical flood periods show that the multi-period Nash efficiency coefficient of the first model ranges from 0.47 to 0.73; while, under the same input boundary conditions, the multi-period Nash efficiency coefficient of the second model ranges from 0.77 to 0.87. Furthermore, the absolute value of the relative error of flood volume output by the second model tends to decrease. For example, in the simulation of a typical flood period, the relative error of flood volume for the first model was -8.63%, while the output value for the relative error of flood volume for the second model was -2.70%.

[0212] The same evaluation index calculation rules were applied to simulated data from the dry season in typical years. In the simulation calculations during the dry season, the Nash efficiency coefficient of the first model fluctuated across multiple years, ranging from -2.13 to 0.63. The second model, receiving the same boundary conditions, shifted its output Nash efficiency coefficient distribution range to 0.10 to 0.76, reducing the frequency of negative coefficients. Analyzing from the perspective of water balance assessment, in a simulation of an extremely dry year, the relative error of the flood volume output by the first model was -52.72%, indicating an underestimation of the actual water volume; for the same period, the relative error of the flood volume calculated by the second model was 1.02%.

[0213] The comparative calculation results output by the above verification process reflect, from a mathematical statistical perspective, the calibration effect of the dynamically adjusted nonlinear exponent in the state-adaptive mechanism on the low-flow characteristic data sequence. When the aquifer is in a low-storage state, the minimum nonlinear exponent generated by the sigmoid continuous function mapping slows down the decay rate of the baseflow, thereby adjusting the model's flow output value during the dry season and reducing the difference in water volume calculation between the model and the measured runoff sequence.

Claims

1. A rainfall-runoff simulation method based on a nonlinear confluence mechanism of underground baseflow, characterized in that, include: Acquire watershed precipitation and evapotranspiration data, perform spatial weighting and quality control processing on them, and obtain the watershed average input sequence; The total runoff is calculated based on the average input sequence of the watershed, and the total runoff is divided into surface runoff, interflow and groundwater runoff based on the free water storage structure. Based on the dimensional transformation of surface runoff and interflow, and using linear regulation and storage mechanisms for runoff calculation, the surface runoff flow and interflow flow were obtained. A groundwater storage balance model is constructed. Based on the groundwater storage at the previous moment, a state adaptive mechanism is adopted to smoothly adjust the nonlinear exponent of the baseflow calculation according to the saturation state of the groundwater storage at the previous moment, and the nonlinear groundwater baseflow is calculated. The groundwater storage at the previous moment is recursively updated based on the difference between the groundwater runoff and the nonlinear groundwater baseflow to obtain the groundwater storage at the current moment. The total runoff process of the watershed is obtained by superimposing the surface runoff, the interflow runoff, and the nonlinear groundwater baseflow, and then processing them through river network lag propagation calculation.

2. The method according to claim 1, characterized in that, Watershed precipitation and evapotranspiration data were acquired, and spatial weighting and quality control processing was performed on them to obtain the watershed average input sequence, including: Spatial weighted data are obtained by using the Thiessen polygon method or the weighted average method to spatially weight the watershed precipitation data and evapotranspiration data. Spatial weighted data are imputed for missing data and outliers are removed. Potential mutations are identified through consistency checks to generate a stable and reliable watershed average input sequence.

3. The method according to claim 2, characterized in that, The specific steps in imputing missing data and removing outliers from spatially weighted data include: For missing data, if it is determined to be a short-term missing data, linear interpolation is used to interpolate and obtain short-term interpolated data; If it is determined to be a long-term missing data, the measured data of the adjacent observation station is obtained as the neighbor station data, and the neighbor station data substitution method is used to interpolate to obtain the long-term interpolated data; A statistical boundary detection method based on the standard deviation principle is used to identify outliers in spatially weighted data and correct them, resulting in corrected spatially weighted data.

4. The method according to claim 1, characterized in that, The total runoff is calculated based on the average input sequence of the watershed, and then divided into surface runoff, interflow, and groundwater runoff, including: Net rainfall is obtained by calculating the evapotranspiration capacity in the average input sequence of the watershed, and the soil moisture content is updated according to the preset three-layer soil structure. If the net rainfall is less than or equal to zero, the total runoff is set to zero. Under the condition of net rainfall, the total runoff is calculated based on the full runoff generation mechanism and combined with the soil moisture status, and the current runoff generation area ratio is determined according to the full runoff generation mechanism. The total runoff is dynamically allocated based on the current free water storage state obtained from the previous time step iteration output or preset initialization, and the total runoff is divided into surface runoff, interflow and groundwater runoff based on the runoff area ratio.

5. The method according to claim 1, characterized in that, Construct a groundwater storage balance model and dynamically update groundwater storage based on groundwater runoff, including: In the groundwater storage balance model, a differential form of the water balance equation is used, and groundwater runoff is taken as the groundwater recharge. Calculate the water volume difference between the groundwater recharge and the nonlinear groundwater baseflow within the current time step; Based on the difference in water volume, the historical groundwater storage at the previous moment is iteratively accumulated to achieve the recursive update of groundwater storage.

6. The method according to claim 5, characterized in that, After obtaining the groundwater reserves, the following is also included: Obtain the preset boundary values ​​between the minimum and maximum groundwater storage. The groundwater storage at the current moment is normalized by using the minimum and maximum groundwater storage to extract the relative saturation of groundwater in the range of zero to one.

7. The method according to claim 6, characterized in that, A state-adaptive mechanism is employed to smoothly adjust the nonlinear exponent of the baseflow calculation according to the saturation state of the groundwater storage, including: The nonlinear exponent is set as a dynamic state variable that depends on the relative saturation of groundwater; Using a preset S-shaped continuous function and groundwater relative saturation as the input variable, a smooth transition mapping is performed between the preset minimum nonlinear index and the maximum nonlinear index to calculate the dynamically adjusted nonlinear index.

8. The method according to claim 7, characterized in that, The nonlinear underground baseflow was calculated, including: Obtain the pre-configured base current coefficient and the preset reference constants used for balance dimensions; The base term is obtained by performing dimensionless calculation of groundwater storage using a preset reference constant; The nonlinear underground baseflow is obtained by performing a power function operation on the base term and the nonlinear exponent, and then multiplying the result of the power function operation with the baseflow coefficient.

9. The method according to claim 8, characterized in that, The nonlinear calculation process also includes applying hydrophysical boundary constraints to each parameter, specifically including: By imposing extreme value constraints on groundwater reserves, the calculated value is limited to the range between the minimum and maximum groundwater reserves, resulting in constrained groundwater reserves. By imposing exponential extremum constraints on the nonlinear exponent, its calculated value is restricted to between the minimum and maximum nonlinear exponents, resulting in a constrained nonlinear exponent. By applying a baseflow extreme value constraint to the nonlinear underground baseflow, the output value is ensured to be no less than the preset minimum positive baseflow and no greater than the preset maximum baseflow, thus obtaining the constrained nonlinear underground baseflow.

10. The method according to claim 9, characterized in that, When applying hydrophysical boundary constraints to each parameter, the minimum nonlinear index used to limit the nonlinear index ranges from 0.6 to 1.0, and the maximum nonlinear index ranges from 1.2 to 2.

5. Furthermore, the minimum groundwater storage capacity is set as the numerical limit of the minimum water storage capacity during the corresponding dry season, and the maximum groundwater storage capacity is set as the numerical limit of the saturated water storage capacity of the corresponding groundwater aquifer.

Images